cloudy
trunk
|
00001 /* This file is part of Cloudy and is copyright (C)1978-2008 by Gary J. Ferland and 00002 * others. For conditions of distribution and use see copyright notice in license.txt */ 00003 /*ContSetIntensity derive intensity of incident continuum */ 00004 /*extin do extinction of incident continuum as set by extinguish command */ 00005 /*sumcon sums L and Q for net incident continuum */ 00006 /*ptrcer show continuum pointers in real time following drive pointers command */ 00007 /*conorm normalize continuum to proper intensity */ 00008 /*qintr integrates Q for any continuum between two limits, used for normalization */ 00009 /*pintr integrates L for any continuum between two limits, used for normalization */ 00010 #include "cddefines.h" 00011 #include "physconst.h" 00012 #include "iso.h" 00013 #include "extinc.h" 00014 #include "noexec.h" 00015 #include "ionbal.h" 00016 #include "hextra.h" 00017 #include "trace.h" 00018 #include "dense.h" 00019 #include "oxy.h" 00020 #include "prt.h" 00021 #include "heavy.h" 00022 #include "rfield.h" 00023 #include "phycon.h" 00024 #include "called.h" 00025 #include "hydrogenic.h" 00026 #include "timesc.h" 00027 #include "secondaries.h" 00028 #include "opacity.h" 00029 #include "thermal.h" 00030 #include "ipoint.h" 00031 #include "atmdat.h" 00032 #include "rt.h" 00033 #include "radius.h" 00034 #include "geometry.h" 00035 #include "grainvar.h" 00036 #include "continuum.h" 00037 00038 /* these are weights used for continuum integration */ 00039 static double aweigh[4]={-0.4305682,-0.1699905, 0.1699905, 0.4305682}; 00040 static double fweigh[4]={ 0.1739274, 0.3260726, 0.3260726, 0.1739274}; 00041 00042 /*conorm normalize continuum to proper intensity */ 00043 STATIC void conorm(void); 00044 00045 /*pintr integrates L for any continuum between two limits, used for normalization */ 00046 STATIC double pintr(double penlo, 00047 double penhi); 00048 00049 /*qintr integrates Q for any continuum between two limits, used for normalization */ 00050 STATIC double qintr(double *qenlo, 00051 double *qenhi); 00052 00053 00054 /*sumcon sums L and Q for net incident continuum */ 00055 STATIC void sumcon(long int il, 00056 long int ih, 00057 realnum *q, 00058 realnum *p, 00059 realnum *panu); 00060 00061 /*extin do extinction of incident continuum as set by extinguish command */ 00062 STATIC void extin(realnum *ex1ryd); 00063 00064 /*ptrcer show continuum pointers in real time following drive pointers command */ 00065 STATIC void ptrcer(void); 00066 00067 /* called by Cloudy to set up continuum 00068 * return value is zero if ok, 1 if no radiation - main would then stop */ 00069 int ContSetIntensity(void) 00070 { 00071 bool lgCheckOK; 00072 00073 long int i, 00074 ip, 00075 j, 00076 n; 00077 00078 realnum EdenHeav, 00079 ex1ryd, 00080 factor, 00081 occ1, 00082 p, 00083 p1, 00084 p2, 00085 p3, 00086 p4, 00087 p5, 00088 p6, 00089 p7, 00090 pgn, 00091 phe, 00092 pheii, 00093 qgn; 00094 00095 realnum xIoniz; 00096 00097 double HCaseBRecCoeff, 00098 wanu[4], 00099 alf, 00100 bet, 00101 fntest, 00102 fsum, 00103 ecrit, 00104 tcompr, 00105 tcomp, 00106 RatioIonizToRecomb, 00107 r3ov2; 00108 00109 double amean, 00110 amean2, 00111 amean3, 00112 peak, 00113 wfun[4]; 00114 00115 /* fraction of beamed continuum that is varies with time */ 00116 double frac_beam_time , frac_beam_time1; 00117 /* fraction of beamed continuum that is constant */ 00118 double frac_beam_const , frac_beam_const1; 00119 /* fraction of continuum that is isotropic */ 00120 double frac_isotropic , frac_isotropic1; 00121 00122 long int nelem , ion; 00123 00124 DEBUG_ENTRY( "ContSetIntensity()" ); 00125 00126 /* set continuum */ 00127 if( trace.lgTrace ) 00128 { 00129 fprintf( ioQQQ, " ContSetIntensity called.\n" ); 00130 } 00131 00132 /* find normalization factors for the continua - this decides whether continuum is 00133 * per unit area of luminosity, and which is desired final product */ 00134 conorm(); 00135 00136 /* define factors to convert rfeld.flux array into photon occupation array OCCNUM 00137 * by multiplication */ 00138 factor = (realnum)(EN1RYD/PI4/FR1RYD/HNU3C2); 00139 00140 /*------------------------------------------------------------- */ 00141 lgCheckOK = true; 00142 00143 for( i=0; i < rfield.nupper; i++ ) 00144 { 00145 /* this was original anu array with no continuum averaging */ 00146 rfield.anu[i] = rfield.AnuOrg[i]; 00147 rfield.ContBoltz[i] = 0.; 00148 fsum = 0.; 00149 amean = 0.; 00150 amean2 = 0.; 00151 amean3 = 0.; 00152 frac_beam_time = 0.; 00153 frac_beam_const = 0.; 00154 frac_isotropic = 0.; 00155 00156 for( j=0; j < 4; j++ ) 00157 { 00158 /* aweigh is symmetric about 0.5 widflx */ 00159 wanu[j] = rfield.anu[i] + rfield.widflx[i]*aweigh[j]; 00160 /* >>chng 02 jul 16, add test on continuum limits - 00161 * this was exceeded when resolution set very large */ 00162 wanu[j] = MAX2( wanu[j] , rfield.emm ); 00163 wanu[j] = MIN2( wanu[j] , rfield.egamry ); 00164 /* >>chng 06 feb 03, the continuum binning can change dramatically 00165 * at some energies - make sure that this cell does not overextend the 00166 * boundaries of its neighbors */ 00167 if( i > 0 && i < rfield.nupper-1 ) 00168 { 00169 wanu[j] = MAX2( wanu[j] , rfield.anu[i-1] + 0.5*rfield.widflx[i-1] ); 00170 wanu[j] = MIN2( wanu[j] , rfield.anu[i+1] - 0.5*rfield.widflx[i+1] ); 00171 } 00172 00173 wfun[j] = fweigh[j]*ffun(wanu[j] , 00174 &frac_beam_time1 , 00175 &frac_beam_const1 , 00176 &frac_isotropic1 ); 00177 /*if( i==76 ) 00178 fprintf(ioQQQ,"DEBUG ffun %li %.4e at %.4e R\n", j , 00179 ffun(wanu[j]) , wanu[j] );*/ 00180 fsum += wfun[j]; 00181 amean += wanu[j]*wfun[j]; 00182 amean2 += wanu[j]*wanu[j]*wfun[j]; 00183 amean3 += wanu[j]*wanu[j]*wanu[j]*wfun[j]; 00184 frac_beam_time += fweigh[j]*frac_beam_time1; 00185 frac_beam_const += fweigh[j]*frac_beam_const1; 00186 frac_isotropic += fweigh[j]*frac_isotropic1; 00187 } 00188 00189 ASSERT( fabs( 1.-frac_beam_time-frac_beam_const-frac_isotropic)< 00190 10.*FLT_EPSILON); 00191 /* This is a fix for ticket #1 */ 00192 if( fsum*rfield.widflx[i] > BIGFLOAT ) 00193 { 00194 fprintf( ioQQQ, "\n Cannot continue. The continuum is far too intense.\n" ); 00195 for( j=0; j < rfield.nspec; j++ ) 00196 { 00197 if( (wfun[j]*rfield.widflx[i] > BIGFLOAT) && ( rfield.nspec > 1 ) ) 00198 { 00199 fprintf( ioQQQ, " Problem is with source number %li\n", j ); 00200 break; 00201 } 00202 } 00203 cdEXIT(EXIT_FAILURE); 00204 } 00205 00206 rfield.flux[i] = (realnum)(fsum*rfield.widflx[i]); 00207 00208 /* save separation into isotropic constant and beamed, and possibly 00209 * time-variable beamed continua */ 00210 rfield.flux_beam_const[i] = rfield.flux[i] * (realnum)frac_beam_const; 00211 rfield.flux_beam_time[i] = rfield.flux[i] * (realnum)frac_beam_time; 00212 rfield.flux_isotropic[i] = rfield.flux[i] * (realnum)frac_isotropic; 00213 00214 if( rfield.flux[i] > 0. ) 00215 { 00216 /*fprintf(ioQQQ,"DEBUG %4li %e %e %.3e %.3e\n", 00217 i , rfield.anu[i] , (amean2/amean) , amean2 , amean );*/ 00218 rfield.anu[i] = (realnum)(amean2/amean); 00219 rfield.anu2[i] = (realnum)(amean3/amean); 00220 /* mesh must be strictly monotonically increasing - make it so */ 00221 if( i > 0 && rfield.anu[i] <= rfield.anu[i-1] ) 00222 { 00223 /* prevent roundoff from allowing i cell to lie below i-1 00224 * cell when continuum mesh is very fine. */ 00225 /* use 2*epsilon to protect against unusual rounding modes */ 00226 rfield.anu[i] = rfield.anu[i-1]*(1.f+2.f*FLT_EPSILON); 00227 rfield.anu2[i] = pow2(rfield.anu[i]); 00228 } 00229 ASSERT( i==0 || rfield.anu[i] > rfield.anu[i-1] ); 00230 /* define array of LOG10( nu(ryd) ) */ 00231 rfield.anulog[i] = (realnum)log10(rfield.anu[i]); 00232 } 00233 00234 else if( rfield.flux[i] == 0. ) 00235 { 00236 rfield.anu2[i] = rfield.anu[i]*rfield.anu[i]; 00237 rfield.anulog[i] = (realnum)log10(rfield.anu[i]); 00238 } 00239 00240 else 00241 { 00242 rfield.anu2[i] = rfield.anu[i]*rfield.anu[i]; 00243 fprintf( ioQQQ, " negative continuum returned at%6ld%10.2e%10.2e\n", 00244 i, rfield.anu[i], rfield.flux[i] ); 00245 lgCheckOK = false; 00246 } 00247 rfield.anu3[i] = rfield.anu2[i]*rfield.anu[i]; 00248 00249 rfield.ConEmitReflec[i] = 0.; 00250 rfield.ConEmitOut[i] = 0.; 00251 rfield.convoc[i] = factor/rfield.widflx[i]/rfield.anu2[i]; 00252 00253 /* following are Compton exchange factors from Tarter */ 00254 alf = 1./(1. + rfield.anu[i]*(1.1792e-4 + 7.084e-10*rfield.anu[i])); 00255 bet = 1. - alf*rfield.anu[i]*(1.1792e-4 + 2.*7.084e-10*rfield.anu[i])/ 00256 4.; 00257 rfield.csigh[i] = (realnum)(alf*rfield.anu2[i]*3.858e-25); 00258 rfield.csigc[i] = (realnum)(alf*bet*rfield.anu[i]*3.858e-25); 00259 } 00260 00261 /*i = 76; 00262 fprintf(ioQQQ,"\nDEBUG %4li %e \n", 00263 i , rfield.anu[i] );*/ 00264 00265 /* >>chng 05 jul 01 add this 00266 * finished with stored continua - return these vectors */ 00267 #if 0 00268 /* commented out since we must conserve energy, and continuum was set with old widflx */ 00269 /* now fix widflx array so that it is correct */ 00270 for( i=1; i<rfield.nupper-1; ++i ) 00271 { 00272 /*rfield.widflx[i] = rfield.anu[i+1] - rfield.anu[i];*/ 00273 rfield.widflx[i] = ((rfield.anu[i+1] - rfield.anu[i]) + (rfield.anu[i] - rfield.anu[i-1]))/2.f; 00274 } 00275 #endif 00276 00277 if( !lgCheckOK ) 00278 { 00279 ShowMe(); 00280 cdEXIT(EXIT_FAILURE); 00281 } 00282 00283 if( trace.lgTrace && trace.lgComBug ) 00284 { 00285 fprintf( ioQQQ, "\n\n Compton heating, cooling coefficients \n" ); 00286 for( i=0; i < rfield.nupper; i += 2 ) 00287 { 00288 fprintf( ioQQQ, "%6ld%10.2e%10.2e%10.2e", i, rfield.anu[i], 00289 rfield.csigh[i], rfield.csigc[i] ); 00290 } 00291 fprintf( ioQQQ, "\n" ); 00292 } 00293 00294 /* option to check frequencies in real time, drive pointers command, 00295 * routine is below, is file static */ 00296 if( trace.lgPtrace ) 00297 ptrcer(); 00298 00299 /* extinguish continuum if set on */ 00300 extin(&ex1ryd); 00301 00302 /* now find peak of hydrogen ionizing continuum - for PDR calculations 00303 * this will remain equal to 1 since the loop will not execute */ 00304 prt.ipeak = 1; 00305 peak = 0.; 00306 00307 for( i=iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH1s]-1; i < rfield.nupper; i++ ) 00308 { 00309 if( rfield.flux[i]*rfield.anu[i]/rfield.widflx[i] > (realnum)peak ) 00310 { 00311 /* prt.ipeak points to largest f_nu at H-ionizing energies 00312 * and is passed to other parts of code */ 00313 /* i+1 to keep ipeak on fortran version of energy array */ 00314 prt.ipeak = i+1; 00315 peak = rfield.flux[i]*rfield.anu[i]/rfield.widflx[i]; 00316 } 00317 } 00318 00319 /* find highest energy to consider in continuum flux array 00320 * peak is the peak product nu*flux */ 00321 peak = rfield.flux[prt.ipeak-1]/rfield.widflx[prt.ipeak-1]* 00322 rfield.anu2[prt.ipeak-1]; 00323 00324 /* say what type of cpu this is, if desired */ 00325 if( trace.lgTrace ) 00326 { 00327 fprintf( ioQQQ, " ContSetIntensity: The peak of the H-ion continuum is at %.3e Ryd - its value is %.2e\n", 00328 rfield.anu[prt.ipeak-1] , peak); 00329 } 00330 00331 if( peak > 1e38 ) 00332 { 00333 fprintf( ioQQQ, " PROBLEM DISASTER The continuum is too intense to compute. Use a fainter continuum. (This is the nu*f_nu test)\n" ); 00334 fprintf( ioQQQ, " Sorry.\n" ); 00335 cdEXIT(EXIT_FAILURE); 00336 } 00337 00338 /* FluxFaint set in zero.c; normally 1e-10 */ 00339 /* this will be faintest level of continuum we want to consider. 00340 * peak was set above, is peak of hydrogen ionizing radiation field, 00341 * and is zero if no H-ionizing radiation */ 00342 fntest = peak*rfield.FluxFaint; 00343 { 00344 enum {DEBUG_LOC=false}; 00345 /* print flux array then quit */ 00346 if( DEBUG_LOC ) 00347 { 00348 for( i=0; i<rfield.nupper; ++i ) 00349 { 00350 fprintf(ioQQQ," consetintensityBUGGG\t%.2e\t%.2e\n" , 00351 rfield.anu[i] , rfield.flux[i]/rfield.widflx[i] ); 00352 } 00353 cdEXIT(EXIT_SUCCESS); 00354 } 00355 } 00356 00357 if( fntest > 0. ) 00358 { 00359 /* this test is not done in pdr conditions where NO H-ionizing radiation, 00360 * since fntest is zero*/ 00361 i = rfield.nupper; 00362 /* >>chng 05 aug 16, change ipeak to ipeak+3, to avoid possible one-off bugs 00363 * where continuum barely goes into H-ionizing radiation */ 00364 while( i > prt.ipeak+3 && 00365 rfield.flux[i-1]*rfield.anu2[i-1]/rfield.widflx[i-1] < (realnum)fntest ) 00366 { 00367 --i; 00368 } 00369 } 00370 else 00371 { 00372 /* when no H-ionizing radiation set to Lyman edge */ 00373 /* >>chng 05 aug 16, change ipeak to ipeak+3, to avoid possible one-off bugs 00374 * where continuum barely goes into H-ionizing radiation */ 00375 i = iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH1s]+3; 00376 } 00377 00378 /* 00379 * this line of code dates from 1979 and IOA Cambridge. removed july 19 95 00380 * I think it was the last line of the original Cambridge Fortran source 00381 nflux = MAX( ineon(1)+4 , i ) 00382 */ 00383 00384 /* >>chng 99 apr 28, reinstate the rfield.FluxFaint limit with nflux */ 00385 rfield.nflux = i; 00386 00387 /* trim down nflux, was set to rfield.nupper, the dimension of all vectors, in zero.c, 00388 * in ContCreatePointers was set to nupper, the number of cells needed to get up to the 00389 * high-energy limit of the code */ 00390 while( rfield.flux[rfield.nflux-1] < SMALLFLOAT && rfield.nflux > 1 ) 00391 { 00392 --rfield.nflux; 00393 } 00394 /* make sure we go high enough, avoid 1-off bugs as in draine field */ 00395 ++rfield.nflux; 00396 00397 if( rfield.nflux == 1 ) 00398 { 00399 fprintf( ioQQQ, " PROBLEM DISASTER This incident continuum appears to have no radiation.\n" ); 00400 fprintf( ioQQQ, " Sorry.\n" ); 00401 cdEXIT(EXIT_FAILURE); 00402 } 00403 00404 /* >>chng 04 oct 10, add this limit - arrays will malloc to nupper, but will add unit 00405 * continuum to [nflux] - this must be within array */ 00406 rfield.nflux = MIN2( rfield.nupper-1 , rfield.nflux ); 00407 00408 /* >>chng 05 mar 01, nfine should not extend beyond continuum 00409 * make sure fine opacity scale does not extend beyond continuum we will use */ 00410 rfield.nfine = rfield.nfine_malloc; 00411 while( rfield.nfine > 0 && rfield.fine_anu[rfield.nfine-1] > rfield.anu[rfield.nflux-1] ) 00412 { 00413 --rfield.nfine; 00414 } 00415 00416 /* check that continuum defined everywhere - look for zero's and comment if present */ 00417 continuum.lgCon0 = false; 00418 ip = 0; 00419 for( i=0; i < rfield.nflux; i++ ) 00420 { 00421 if( rfield.flux[i] == 0. ) 00422 { 00423 if( called.lgTalk && !continuum.lgCon0 ) 00424 { 00425 fprintf( ioQQQ, " NOTE Setcon: continuum has zero intensity starting at %11.4e Ryd.\n", 00426 rfield.anu[i] ); 00427 continuum.lgCon0 = true; 00428 } 00429 ++ip; 00430 } 00431 } 00432 00433 if( continuum.lgCon0 && called.lgTalk ) 00434 { 00435 fprintf( ioQQQ, 00436 "%6ld cells in the incident continuum have zero intensity. Problems???\n\n", 00437 ip ); 00438 } 00439 00440 /* check for devastating error in the continuum mesh or intensity */ 00441 lgCheckOK = true; 00442 for( i=1; i < rfield.nflux; i++ ) 00443 { 00444 if( rfield.flux[i] < 0. ) 00445 { 00446 fprintf( ioQQQ, 00447 " PROBLEM DISASTER Continuum has negative intensity at %.4e Ryd=%.2e %4.4s %4.4s\n", 00448 rfield.anu[i], rfield.flux[i], rfield.chLineLabel[i], rfield.chContLabel[i] ); 00449 lgCheckOK = false; 00450 } 00451 else if( rfield.anu[i] <= rfield.anu[i-1] ) 00452 { 00453 fprintf( ioQQQ, 00454 " PROBLEM DISASTER cont_setintensity - internal error - continuum energies not in increasing order: energies follow\n" ); 00455 fprintf( ioQQQ, 00456 "%ld %e %ld %e %ld %e\n", 00457 i -1 , rfield.anu[i-1], i, rfield.anu[i], i +1, rfield.anu[i+1] ); 00458 lgCheckOK = false; 00459 } 00460 } 00461 00462 /* either of the ways lgCheckOK would be set true would be a major internal error */ 00463 if( !lgCheckOK ) 00464 { 00465 ShowMe(); 00466 cdEXIT(EXIT_FAILURE); 00467 } 00468 00469 /* turn off recoil ionization if high energy < 190R */ 00470 if( rfield.anu[rfield.nflux-1] <= 190 ) 00471 { 00472 ionbal.lgCompRecoil = false; 00473 } 00474 00475 /* sum photons and energy, save mean */ 00476 00477 /* sum from low energy to Balmer edge */ 00478 sumcon(1,iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH2p]-1,&rfield.qrad,&prt.pradio,&p1); 00479 00480 /* sum over Balmer continuum */ 00481 sumcon(iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][2],iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH1s]-1,&rfield.qbal,&prt.pbal,&p1); 00482 00483 /* sum from Lyman edge to HeI edge */ 00484 sumcon(iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH1s], 00485 iso.ipIsoLevNIonCon[ipHE_LIKE][ipHELIUM][0]-1,&prt.q,&p,&p2); 00486 00487 /* sum from HeI to HeII edges */ 00488 sumcon(iso.ipIsoLevNIonCon[ipHE_LIKE][ipHELIUM][0], 00489 iso.ipIsoLevNIonCon[ipH_LIKE][1][ipH1s]-1,&rfield.qhe,&phe,&p3); 00490 00491 /* sum from Lyman edge to carbon k-shell */ 00492 sumcon(iso.ipIsoLevNIonCon[ipH_LIKE][1][ipH1s],opac.ipCKshell-1,&rfield.qheii,&pheii,&p4); 00493 00494 /* sum from c k-shell to gamma ray - where pairs start */ 00495 sumcon( opac.ipCKshell , rfield.ipEnerGammaRay-1 , &prt.qx, 00496 &prt.xpow , &p5); 00497 00498 /* complete sum up to high-energy limit */ 00499 sumcon(rfield.ipEnerGammaRay,rfield.nflux,&prt.qgam,&prt.GammaLumin, &p6); 00500 00501 /* find to estimate photoerosion timescale */ 00502 n = ipoint(7.35e5); 00503 sumcon(n,rfield.nflux,&qgn,&pgn,&p7); 00504 timesc.TimeErode = qgn; 00505 00506 /* find Compton temp */ 00507 tcompr = (p1 + p2 + p3 + p4 + p5 + p6)/(prt.pradio + prt.pbal + 00508 p + phe + pheii + prt.xpow + prt.GammaLumin); 00509 00510 tcomp = tcompr/(4.*6.272e-6); 00511 00512 if( trace.lgTrace ) 00513 { 00514 fprintf( ioQQQ, 00515 " mean photon energy=%10.3eR =%10.3eK low, high nu=%12.4e%12.4e\n", 00516 tcompr, tcomp, rfield.anu[0] - rfield.widflx[0]/2., rfield.anu[rfield.nflux-1] + 00517 rfield.widflx[rfield.nflux-1]/2. ); 00518 } 00519 00520 /* this is total power in ionizing radiation, per unit area */ 00521 prt.powion = p + phe + pheii + prt.xpow + prt.GammaLumin; 00522 00523 /* this is the total radiation intensity, erg cm-2 s-1 */ 00524 continuum.TotalLumin = prt.pradio + prt.powion + prt.pbal; 00525 00526 /* this is placed into the line stack on the first zone, then 00527 * reset to zero, to end up with luminosity in the emission lines array. 00528 * at end of iteration it is reset to TotalLumin */ 00529 continuum.totlsv = continuum.TotalLumin; 00530 00531 /* total H-ionizing photon number, */ 00532 rfield.qhtot = prt.q + rfield.qhe + rfield.qheii + prt.qx + prt.qgam; 00533 00534 /* ftotal photon number, all energies */ 00535 rfield.qtot = rfield.qhtot + rfield.qbal + rfield.qrad; 00536 00537 if( prt.powion <= 0. && called.lgTalk ) 00538 { 00539 rfield.lgHionRad = true; 00540 fprintf( ioQQQ, " NOTE There is no hydrogen-ionizing radiation.\n" ); 00541 fprintf( ioQQQ, " Was this intended?\n\n" ); 00542 /* check if any Balmer ionizing radiation since even metals will be 00543 * totally neutral */ 00544 if( prt.pbal <=0. && called.lgTalk ) 00545 { 00546 fprintf( ioQQQ, " NOTE There is no Balmer continuum radiation.<<<<\n" ); 00547 fprintf( ioQQQ, " Was this intended?\n\n" ); 00548 } 00549 } 00550 00551 else 00552 { 00553 rfield.lgHionRad = false; 00554 } 00555 00556 /* option to add energy deposition due to fast neutrons, 00557 * entered as fraction of total photon luminosity */ 00558 if( hextra.lgNeutrnHeatOn ) 00559 { 00560 /* hextra.totneu is erg cm-2 s-1 in neutrons 00561 * hextra.effneu - efficiency default is unity */ 00562 hextra.totneu = (realnum)(hextra.effneu*continuum.TotalLumin* 00563 pow((realnum)10.f,hextra.frcneu)); 00564 } 00565 else 00566 { 00567 hextra.totneu = (realnum)0.; 00568 } 00569 00570 /* temp correspond to energy density, printed in STARTR */ 00571 phycon.TEnerDen = pow(continuum.TotalLumin/SPEEDLIGHT/7.56464e-15,0.25); 00572 00573 /* sanity check for single blackbody, that energy density temperature 00574 * is smaller than black body temperature */ 00575 if( rfield.nspec==1 && 00576 strcmp( rfield.chSpType[rfield.nspec-1], "BLACK" )==0 ) 00577 { 00578 /* single black body, now confirm that TEnerDen is less than this temperature, 00579 * >>>chng 99 may 02, 00580 * in lte these are very very close, factor of 1.00001 allows for numerical 00581 * errors, and apparently slightly different atomic coef in different parts 00582 * of code. eventaully all mustuse physonst.h and agree exactly */ 00583 if( phycon.TEnerDen > 1.0001f*rfield.slope[rfield.nspec-1] ) 00584 { 00585 fprintf( ioQQQ, 00586 "\n WARNING: The energy density temperature (%g) is greater than the" 00587 " black body temperature (%g). This is unphysical.\n\n", 00588 phycon.TEnerDen , rfield.slope[rfield.nspec-1]); 00589 } 00590 } 00591 00592 /* incident continuum nu*f_nu at Hbeta and Ly-alpha */ 00593 continuum.cn4861 = (realnum)(ffun(0.1875)*HPLANCK*FR1RYD*0.1875*0.1875); 00594 continuum.cn1216 = (realnum)(ffun(0.75)*HPLANCK*FR1RYD*0.75*0.75); 00595 continuum.sv4861 = continuum.cn4861; 00596 continuum.sv1216 = continuum.cn1216; 00597 /* flux density in V, erg / s / cm2 / hz */ 00598 continuum.fluxv = (realnum)(ffun(0.1643)*HPLANCK*0.1643); 00599 continuum.fbeta = (realnum)(ffun(0.1875)*HPLANCK*0.1875*6.167e14); 00600 00601 /* flux density nu*Fnu = erg / s / cm2 00602 * EX1RYD is optional extinction factor at 1 ryd */ 00603 prt.fx1ryd = (realnum)(ffun(1.000)*HPLANCK*ex1ryd*FR1RYD); 00604 00605 /* check for plasma frequency - then zero out incident continuum 00606 * for energies below this 00607 * this is critical electron density, shielding of incident continuum 00608 * if electron density is greater than this */ 00609 ecrit = POW2(rfield.anu[0]/2.729e-12); 00610 00611 if( dense.gas_phase[ipHYDROGEN]*1.2 > ecrit ) 00612 { 00613 rfield.lgPlasNu = true; 00614 rfield.plsfrq = (realnum)(2.729e-12*sqrt(dense.gas_phase[ipHYDROGEN]*1.2)); 00615 rfield.plsfrqmax = rfield.plsfrq; 00616 rfield.ipPlasma = ipoint(rfield.plsfrq); 00617 00618 /* save max pointer too */ 00619 rfield.ipPlasmax = rfield.ipPlasma; 00620 00621 /* now loop over all affected energies, setting incident continuum 00622 * to zero there, and counting all as reflected */ 00623 /* >>chng 01 jul 14, from i < ipPlasma to ipPlasma-1 - 00624 * when ipPlasma is 1 plasma freq is not on energy scale */ 00625 for( i=0; i < rfield.ipPlasma-1; i++ ) 00626 { 00627 /* count as reflected incident continuum */ 00628 rfield.ConRefIncid[i] = rfield.flux[i]; 00629 /* set continuum to zero there */ 00630 rfield.flux_beam_const[i] = 0.; 00631 rfield.flux_beam_time[i] = 0.; 00632 rfield.flux_isotropic[i] = 0.; 00633 rfield.flux[i] = rfield.flux_beam_const[i] + rfield.flux_beam_time[i] + 00634 rfield.flux_isotropic[i]; 00635 } 00636 } 00637 else 00638 { 00639 rfield.lgPlasNu = false; 00640 /* >>chng 01 jul 14, from 0 to 1 - 1 is the first array element on the F scale, 00641 * ipoint would return this, so rest of code assumes ipPlasma is 1 plus correct index */ 00642 rfield.ipPlasma = 1; 00643 rfield.plsfrqmax = 0.; 00644 rfield.plsfrq = 0.; 00645 } 00646 00647 if( rfield.ipPlasma > 1 && called.lgTalk ) 00648 { 00649 fprintf( ioQQQ, 00650 " !The plasma frequency is %.2e Ryd. The incident continuum is set to 0 below this.\n", 00651 rfield.plsfrq ); 00652 } 00653 00654 rfield.occmax = 0.; 00655 rfield.tbrmax = 0.; 00656 for( i=0; i < rfield.nupper; i++ ) 00657 { 00658 /* set up occupation number array */ 00659 rfield.OccNumbIncidCont[i] = rfield.flux[i]*rfield.convoc[i]; 00660 if( rfield.OccNumbIncidCont[i] > rfield.occmax ) 00661 { 00662 rfield.occmax = rfield.OccNumbIncidCont[i]; 00663 rfield.occmnu = rfield.anu[i]; 00664 } 00665 /* following product is continuum brightness temperature */ 00666 if( rfield.OccNumbIncidCont[i]*TE1RYD*rfield.anu[i] > rfield.tbrmax ) 00667 { 00668 rfield.tbrmax = (realnum)(rfield.OccNumbIncidCont[i]*TE1RYD*rfield.anu[i]); 00669 rfield.tbrmnu = rfield.anu[i]; 00670 } 00671 /* save continuum for next iteration */ 00672 rfield.flux_total_incident[i] = rfield.flux[i]; 00673 rfield.flux_beam_const_save[i] = rfield.flux_beam_const[i]; 00674 rfield.flux_time_beam_save[i] = rfield.flux_beam_time[i]; 00675 rfield.flux_isotropic_save[i] = rfield.flux_isotropic[i]; 00676 /*fprintf(ioQQQ,"DEBUG type cont %li\t%.3e\t%.2e\t%.2e\t%.2e\t%.2e\n", 00677 i, rfield.anu[i], 00678 rfield.flux[i],rfield.flux_beam_const[i],rfield.flux_beam_time[i], 00679 rfield.flux_isotropic[i]); 00680 fflush(ioQQQ);*/ 00681 } 00682 00683 /* if continuum brightness temp is large, where does it fall below 00684 * 1e4k??? */ 00685 if( rfield.tbrmax > 1e4 ) 00686 { 00687 i = ipoint(rfield.tbrmnu)-1; 00688 while( i < rfield.nupper-1 && (rfield.OccNumbIncidCont[i]*TE1RYD* 00689 rfield.anu[i] > 1e4) ) 00690 { 00691 ++i; 00692 } 00693 rfield.tbr4nu = rfield.anu[i]; 00694 } 00695 else 00696 { 00697 rfield.tbr4nu = 0.; 00698 } 00699 00700 /* if continuum occ num is large, where does it fall below 1? */ 00701 if( rfield.occmax > 1. ) 00702 { 00703 i = ipoint(rfield.occmnu)-1; 00704 while( i < rfield.nupper && (rfield.OccNumbIncidCont[i] > 1.) ) 00705 { 00706 ++i; 00707 } 00708 rfield.occ1nu = rfield.anu[i]; 00709 } 00710 else 00711 { 00712 rfield.occ1nu = 0.; 00713 } 00714 00715 /* remember if incident radiation field is less than 10*Habing ISM */ 00716 /* >>chng 06 aug 01, change this test from continuum.TotalLumin to 00717 * energy in balmer and ionizing continuum, since this is the true habing field 00718 * and is the continuum that interacts with gas. When CMB set this 00719 * tests on total did not trigger due to cold blackbody, which has little 00720 * effect on gas, other than compton */ 00721 if( (prt.powion + prt.pbal) < 1.8e-2 ) 00722 { 00723 /* thermal.ConstTemp def is zero, set pos when constant temperature is set */ 00724 rfield.lgHabing = true; 00725 /* >>chng 06 aug 01 also print warning if substantially below Habing, this may be a mistake */ 00726 if( ((prt.powion + prt.pbal) < 1.8e-12) && 00727 /* this is test for not constant temperature */ 00728 (!thermal.lgTemperatureConstant) ) 00729 { 00730 fprintf( ioQQQ, "\n >>>\n" 00731 " >>> NOTE The incident continuum is surprisingly faint.\n" ); 00732 fprintf( ioQQQ, 00733 " >>> The total energy in the Balmer and Lyman continua is %.2e erg cm-2 s-1.\n" 00734 ,(prt.powion + prt.pbal)); 00735 fprintf( ioQQQ, " >>> This is many orders of magnitude fainter than the ISM galactic background.\n" ); 00736 fprintf( ioQQQ, " >>> This seems unphysical - please check that the continuum intensity has been properly set.\n" ); 00737 fprintf( ioQQQ, " >>> YOU MAY BE MAKING A BIG MISTAKE!!\n >>>\n\n\n\n" ); 00738 } 00739 } 00740 00741 /* fix ionization parameter (per hydrogen) at inner edge */ 00742 rfield.uh = (realnum)(rfield.qhtot/dense.gas_phase[ipHYDROGEN]/SPEEDLIGHT); 00743 rfield.uheii = (realnum)((rfield.qheii + prt.qx)/dense.gas_phase[ipHYDROGEN]/SPEEDLIGHT); 00744 if( rfield.uh > 1.e10 ) 00745 { 00746 fprintf( ioQQQ, "\n >>>\n" 00747 " NOTE The incident continuum is surprisingly intense.\n" ); 00748 fprintf( ioQQQ, 00749 " >>> The hydrogen ionization parameter is %.2e [dimensionless].\n" 00750 , rfield.uh ); 00751 fprintf( ioQQQ, " This is many orders of magnitude brighter than commonly seen.\n" ); 00752 fprintf( ioQQQ, " This seems unphysical - please check that the continuum intensity has been properly set.\n" ); 00753 fprintf( ioQQQ, " YOU MAY BE MAKING A BIG MISTAKE!!\n >>>\n\n\n\n" ); 00754 } 00755 00756 /* guess first temperature and neutral h density */ 00757 if( thermal.ConstTemp > 0. ) 00758 { 00759 phycon.te = thermal.ConstTemp; 00760 } 00761 else 00762 { 00763 if( rfield.uh > 0. ) 00764 { 00765 phycon.te = (20000.+log10(rfield.uh)*5000.); 00766 phycon.te = MAX2(8000. , phycon.te ); 00767 } 00768 else 00769 { 00770 phycon.te = 1000.; 00771 } 00772 } 00773 00774 /* this is an option to stop after printing header only */ 00775 if( noexec.lgNoExec ) 00776 return(0); 00777 00778 /* estimate secondary ionization rate - probably 0, but possible extra 00779 * SetCsupra set with "set csupra" */ 00780 /* >>>chng 99 apr 29, added cosmic ray ionization since this is used to get 00781 * helium ionization fraction, and was zero in pdr, so He turned off at start, 00782 * and never turned back on */ 00783 /* coef on cryden is from highen.c */ 00784 for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem ) 00785 { 00786 for( ion=0; ion<nelem+1; ++ion ) 00787 { 00788 secondaries.csupra[nelem][ion] = 00789 secondaries.SetCsupra + hextra.cryden*2e-9f; 00790 } 00791 } 00792 00793 /********************************************************************* 00794 * * 00795 * estimate hydrogen's level of ionization * 00796 * * 00797 *********************************************************************/ 00798 00799 /* create fake ionization balance, but will conserve number of hydrogens */ 00800 dense.xIonDense[ipHYDROGEN][0] = 0.; 00801 dense.xIonDense[ipHYDROGEN][1] = dense.gas_phase[ipHYDROGEN]; 00802 /* this must be zero since PresTotCurrent will do radiation pressure due to H */ 00803 StatesElem[ipH_LIKE][ipHYDROGEN][ipH1s].Pop = 0.; 00804 00805 /* "extra" electrons from command line, or assumed residual electrons */ 00806 double EdenExtraLocal = dense.EdenExtra + 00807 /* if we are in a molecular cloud the current logic could badly fail 00808 * do not let electron density fall below 1e-7 of H density */ 00809 1e-7*dense.gas_phase[ipHYDROGEN]; 00810 dense.eden = dense.xIonDense[ipHYDROGEN][1] + EdenExtraLocal; 00811 00812 /* hydrogen case B recombination coefficient */ 00813 HCaseBRecCoeff = (-9.9765209 + 0.158607055*phycon.telogn[0] + 0.30112749* 00814 phycon.telogn[1] - 0.063969007*phycon.telogn[2] + 0.0012691546* 00815 phycon.telogn[3])/(1. + 0.035055422*phycon.telogn[0] - 00816 0.037621619*phycon.telogn[1] + 0.0076175048*phycon.telogn[2] - 00817 0.00023432613*phycon.telogn[3]); 00818 HCaseBRecCoeff = pow(10.,HCaseBRecCoeff)/phycon.te; 00819 00820 double CollIoniz = t_ADfA::Inst().coll_ion(1,1,phycon.te); 00821 double OtherIonization = rfield.qhtot*2e-18 + secondaries.csupra[ipHYDROGEN][0]; 00822 double RatioIoniz = 00823 (CollIoniz*dense.eden+OtherIonization)/(HCaseBRecCoeff*dense.eden); 00824 if( RatioIoniz<1e-3 ) 00825 { 00826 /* very low ionization solution */ 00827 dense.xIonDense[ipHYDROGEN][1] = (realnum)( 00828 dense.gas_phase[ipHYDROGEN]*RatioIoniz); 00829 dense.xIonDense[ipHYDROGEN][0] = dense.gas_phase[ipHYDROGEN] - 00830 dense.xIonDense[ipHYDROGEN][1]; 00831 ASSERT( dense.xIonDense[ipHYDROGEN][0]>=0. && 00832 dense.xIonDense[ipHYDROGEN][0]<=dense.gas_phase[ipHYDROGEN]); 00833 ASSERT( dense.xIonDense[ipHYDROGEN][1]>=0. && 00834 dense.xIonDense[ipHYDROGEN][1]<dense.gas_phase[ipHYDROGEN]); 00835 } 00836 else if( RatioIoniz>1e3 ) 00837 { 00838 /* very high ionization solution */ 00839 dense.xIonDense[ipHYDROGEN][0] = (realnum)( 00840 dense.gas_phase[ipHYDROGEN]/RatioIoniz); 00841 dense.xIonDense[ipHYDROGEN][1] = dense.gas_phase[ipHYDROGEN] - 00842 dense.xIonDense[ipHYDROGEN][0]; 00843 ASSERT( dense.xIonDense[ipHYDROGEN][0]>=0. && 00844 dense.xIonDense[ipHYDROGEN][0]<dense.gas_phase[ipHYDROGEN]); 00845 ASSERT( dense.xIonDense[ipHYDROGEN][1]>=0. && 00846 dense.xIonDense[ipHYDROGEN][1]<=dense.gas_phase[ipHYDROGEN]); 00847 } 00848 else 00849 { 00850 /* intermediate ionization - solve quadratic */ 00851 double alpha = HCaseBRecCoeff + CollIoniz ; 00852 double beta = HCaseBRecCoeff*EdenExtraLocal + OtherIonization - 00853 dense.gas_phase[ipHYDROGEN]*CollIoniz; 00854 double gamma = EdenExtraLocal*CollIoniz - 00855 dense.gas_phase[ipHYDROGEN]*(OtherIonization+EdenExtraLocal*CollIoniz); 00856 00857 double discriminant = POW2(beta) - 4.*alpha*gamma; 00858 if( discriminant <0 ) 00859 { 00860 /* oops */ 00861 fprintf(ioQQQ," DISASTER PROBLEM cont_initensity found negative discriminant.\n"); 00862 TotalInsanity(); 00863 } 00864 00865 dense.xIonDense[ipHYDROGEN][1] = (realnum)((-beta + sqrt(discriminant))/(2.*alpha)); 00866 if( dense.xIonDense[ipHYDROGEN][1]> dense.gas_phase[ipHYDROGEN] ) 00867 { 00868 /* oops */ 00869 fprintf(ioQQQ," DISASTER PROBLEM cont_initensity found n(H+)>n(H).\n"); 00870 TotalInsanity(); 00871 } 00872 dense.xIonDense[ipHYDROGEN][0] = dense.gas_phase[ipHYDROGEN] - 00873 dense.xIonDense[ipHYDROGEN][1]; 00874 if( dense.xIonDense[ipHYDROGEN][0]<= 0 ) 00875 { 00876 /* oops */ 00877 fprintf(ioQQQ," DISASTER PROBLEM cont_initensity found n(H0)<0.\n"); 00878 TotalInsanity(); 00879 } 00880 } 00881 00882 dense.eden = dense.xIonDense[ipHYDROGEN][1] + (realnum)EdenExtraLocal; 00883 if( dense.eden <= SMALLFLOAT ) 00884 { 00885 /* no electrons! */ 00886 fprintf(ioQQQ,"\n PROBLEM DISASTER - this simulation has no source" 00887 " of ionization. The electron density is zero. Consider " 00888 "adding a source of ionization such as cosmic rays.\n\n"); 00889 cdEXIT(EXIT_FAILURE); 00890 } 00891 00892 if( dense.xIonDense[ipHYDROGEN][1] > 1e-30 ) 00893 { 00894 StatesElem[ipH_LIKE][ipHYDROGEN][ipH1s].Pop = dense.xIonDense[ipHYDROGEN][0]/dense.xIonDense[ipHYDROGEN][1]; 00895 } 00896 else 00897 { 00898 StatesElem[ipH_LIKE][ipHYDROGEN][ipH1s].Pop = 0.; 00899 } 00900 00901 /* now save estimates of whether induced recombination is going 00902 * to be important -this is a code pacesetter since GammaBn is slower 00903 * than GammaK */ 00904 hydro.lgHInducImp = false; 00905 for( i=ipH1s; i < iso.numLevels_max[ipH_LIKE][ipHYDROGEN]; i++ ) 00906 { 00907 if( rfield.OccNumbIncidCont[iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][i]-1] > 0.01 ) 00908 hydro.lgHInducImp = true; 00909 } 00910 00911 /******************************************************************* 00912 * * 00913 * estimate helium's level of ionization * 00914 * * 00915 *******************************************************************/ 00916 00917 /* only if helium is turned on */ 00918 if( dense.lgElmtOn[ipHELIUM] ) 00919 { 00920 /* next estimate level of helium singly ionized */ 00921 xIoniz = (realnum)t_ADfA::Inst().coll_ion(2,2,phycon.te); 00922 /* >>chng 05 jan 05, add cosmic rays */ 00923 xIoniz = (realnum)(xIoniz*dense.eden + rfield.qhe*1e-18 + secondaries.csupra[ipHELIUM][1] ); 00924 double RecTot = HCaseBRecCoeff*dense.eden; 00925 RatioIonizToRecomb = xIoniz/RecTot; 00926 00927 /* now estimate level of helium doubly ionized */ 00928 xIoniz = (realnum)t_ADfA::Inst().coll_ion(2,1,phycon.te); 00929 /* >>chng 05 jan 05, add cosmic rays */ 00930 xIoniz = (realnum)(xIoniz*dense.eden + rfield.qheii*1e-18 + ionbal.CosRayIonRate ); 00931 00932 /* rough charge dependence */ 00933 RecTot *= 4.; 00934 r3ov2 = xIoniz/RecTot; 00935 00936 /* now set level of helium ionization */ 00937 if( RatioIonizToRecomb > 0. ) 00938 { 00939 dense.xIonDense[ipHELIUM][1] = (realnum)(dense.gas_phase[ipHELIUM]/(1./RatioIonizToRecomb + 1. + r3ov2)); 00940 dense.xIonDense[ipHELIUM][0] = (realnum)(dense.xIonDense[ipHELIUM][1]/RatioIonizToRecomb); 00941 } 00942 else 00943 { 00944 /* no He ionizing radiation */ 00945 dense.xIonDense[ipHELIUM][1] = 0.; 00946 dense.xIonDense[ipHELIUM][0] = dense.gas_phase[ipHELIUM]; 00947 } 00948 00949 dense.xIonDense[ipHELIUM][2] = (realnum)(dense.xIonDense[ipHELIUM][1]*r3ov2); 00950 00951 if( dense.xIonDense[ipHELIUM][2] > 1e-30 ) 00952 { 00953 StatesElem[ipH_LIKE][1][ipH1s].Pop = dense.xIonDense[ipHELIUM][1]/dense.xIonDense[ipHELIUM][2]; 00954 } 00955 else 00956 { 00957 StatesElem[ipH_LIKE][1][ipH1s].Pop = 0.; 00958 } 00959 } 00960 else 00961 { 00962 /* case where helium is turned off */ 00963 dense.xIonDense[ipHELIUM][1] = 0.; 00964 dense.xIonDense[ipHELIUM][0] = 0.; 00965 dense.xIonDense[ipHELIUM][2] = 0.; 00966 StatesElem[ipH_LIKE][1][ipH1s].Pop = 0.; 00967 } 00968 00969 /* update electron density */ 00970 dense.eden = dense.xIonDense[ipHYDROGEN][1] + 00971 dense.xIonDense[ipHELIUM][1] + 2.*dense.xIonDense[ipHELIUM][2]; 00972 00973 /* fix number of stages of ionization */ 00974 for( nelem=ipHYDROGEN; nelem < LIMELM; nelem++ ) 00975 { 00976 if( dense.lgElmtOn[nelem] ) 00977 { 00978 dense.IonLow[nelem] = 0; 00979 /* 00980 * IonHigh[n] is the highest stage of ionization present 00981 * the IonHigh array index is on the C scale, so [0] is hydrogen 00982 * the value is also on the C scale, so element [nelem] can range 00983 * from 0 to nelem+1 00984 */ 00985 dense.IonHigh[nelem] = nelem + 1; 00986 /* >>chng 04 jan 13, add this test, caught by Orly Gnat */ 00987 /* check on actual zero abundances of lower stages - this will only 00988 * happen when ionization is set with element ionization command */ 00989 /* >>chng 04 jun 03, move this test here from conv ionopac do */ 00990 if( dense.lgSetIoniz[nelem] ) 00991 { 00992 while( dense.SetIoniz[nelem][dense.IonLow[nelem]] == 0. ) 00993 ++dense.IonLow[nelem]; 00994 while( dense.SetIoniz[nelem][dense.IonHigh[nelem]] == 0. ) 00995 --dense.IonHigh[nelem]; 00996 } 00997 } 00998 else 00999 { 01000 /* this element is turned off, make stages impossible */ 01001 dense.IonLow[nelem] = -1; 01002 dense.IonHigh[nelem] = -1; 01003 } 01004 } 01005 01006 /* estimate electrons from heavies, assuming each at least 01007 * 1 times ionized */ 01008 EdenHeav = 0.; 01009 realnum atomFrac = dense.xIonDense[ipHYDROGEN][0]/dense.gas_phase[ipHYDROGEN]; 01010 realnum firstIonFrac = dense.xIonDense[ipHYDROGEN][1]/dense.gas_phase[ipHYDROGEN]; 01011 for( nelem=ipLITHIUM; nelem < LIMELM; nelem++ ) 01012 { 01013 if( dense.lgElmtOn[nelem] ) 01014 { 01015 dense.xIonDense[nelem][0] = dense.gas_phase[nelem] * atomFrac; 01016 dense.xIonDense[nelem][1] = dense.gas_phase[nelem] * firstIonFrac; 01017 EdenHeav += dense.xIonDense[nelem][1]; 01018 } 01019 } 01020 01021 /* estimate of electron density */ 01022 dense.eden = 01023 dense.xIonDense[ipHYDROGEN][1] + dense.xIonDense[ipHELIUM][1] + 01024 2.*dense.xIonDense[ipHELIUM][2] + EdenHeav + dense.EdenExtra; 01025 /* >>chng 05 jan 05, insure positive eden */ 01026 dense.eden = MAX2( SMALLFLOAT , dense.eden ); 01027 01028 if( dense.EdenSet > 0. ) 01029 { 01030 dense.eden = dense.EdenSet; 01031 } 01032 01033 dense.EdenHCorr = dense.eden; 01034 01035 if( dense.eden < 0. ) 01036 { 01037 fprintf( ioQQQ, " PROBLEM DISASTER Negative electron density results in ContSetIntensity.\n" ); 01038 fprintf( ioQQQ, "%10.2e%10.2e%10.2e%10.2e%10.2e%10.2e\n", 01039 dense.eden, dense.xIonDense[ipHYDROGEN][1], dense.xIonDense[ipHELIUM][1], 01040 dense.xIonDense[ipHELIUM][2], dense.gas_phase[ipCARBON], dense.EdenExtra ); 01041 TotalInsanity(); 01042 } 01043 01044 dense.EdenTrue = dense.eden; 01045 01046 if( trace.lgTrace ) 01047 { 01048 fprintf( ioQQQ, 01049 " ContSetIntensity sets initial EDEN to %.4e, contributors H+=%.2e He+, ++= %.2e %.2e Heav %.2e extra %.2e\n", 01050 dense.eden , 01051 dense.xIonDense[ipHYDROGEN][1], 01052 dense.xIonDense[ipHELIUM][1], 01053 2.*dense.xIonDense[ipHELIUM][2], 01054 EdenHeav, 01055 dense.EdenExtra); 01056 } 01057 01058 /* next two just to make sure some values are set */ 01059 /* this sets values of pressure.PresTotlCurr */ 01060 /*PresTotCurrent();*/ 01061 01062 occ1 = (realnum)(prt.fx1ryd/HNU3C2/PI4/FR1RYD); 01063 01064 /* what is occupation number at 1 Ryd? */ 01065 if( occ1 > 1. ) 01066 { 01067 rfield.lgOcc1Hi = true; 01068 } 01069 else 01070 { 01071 rfield.lgOcc1Hi = false; 01072 } 01073 01074 if( trace.lgTrace && trace.lgConBug ) 01075 { 01076 /* print some useful pointers to ionization edges */ 01077 fprintf( ioQQQ, " H2,1=%5ld%5ld NX=%5ld IRC=%5ld\n", 01078 iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][2], 01079 iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH1s], 01080 opac.ipCKshell, 01081 ionbal.ipCompRecoil[ipHYDROGEN][0] ); 01082 01083 fprintf( ioQQQ, " CARBON" ); 01084 for( i=0; i < 6; i++ ) 01085 { 01086 fprintf( ioQQQ, "%5ld", Heavy.ipHeavy[ipCARBON][i] ); 01087 } 01088 fprintf( ioQQQ, "\n" ); 01089 01090 fprintf( ioQQQ, " OXY" ); 01091 for( i=0; i < 8; i++ ) 01092 { 01093 fprintf( ioQQQ, "%5ld", Heavy.ipHeavy[ipOXYGEN][i] ); 01094 } 01095 fprintf( ioQQQ, "%5ld%5ld%5ld\n", opac.ipo3exc[0], 01096 oxy.i2d, oxy.i2p ); 01097 01098 fprintf( ioQQQ, 01099 "\n\n PHOTONS PER CELL (NOT RYD)\n" ); 01100 fprintf( ioQQQ, 01101 "\n\n nu, flux, wid, occ \n" ); 01102 fprintf( ioQQQ, 01103 " " ); 01104 01105 for( i=0; i < rfield.nflux; i++ ) 01106 { 01107 fprintf( ioQQQ, "%4ld%10.2e%10.2e%10.2e%10.2e", i, 01108 rfield.anu[i], rfield.flux[i], rfield.widflx[i], 01109 rfield.OccNumbIncidCont[i] ); 01110 } 01111 fprintf( ioQQQ, " \n" ); 01112 } 01113 01114 /* zero out some continua related to the ots rates, 01115 * prototype and routine in RT_OTS_Update. This is done here since summed cont will 01116 * be set to rfield */ 01117 RT_OTS_Zero(); 01118 01119 /* >>chng 05 jul 01, do this 01120 * we are finished with these stored continuum arrays - free them */ 01121 for( rfield.ipspec=0; rfield.ipspec < rfield.nspec; rfield.ipspec++ ) 01122 { 01123 if( rfield.lgContMalloc[rfield.ipspec] ) 01124 { 01125 free(rfield.tNuRyd[rfield.ipspec] ); 01126 free(rfield.tslop[rfield.ipspec] ); 01127 free(rfield.tFluxLog[rfield.ipspec] ); 01128 rfield.lgContMalloc[rfield.ipspec] = false; 01129 } 01130 } 01131 01132 if( trace.lgTrace ) 01133 { 01134 fprintf( ioQQQ, " ContSetIntensity returns, nflux=%5ld anu(nflux)=%11.4e eden=%10.2e\n", 01135 rfield.nflux, rfield.anu[rfield.nflux-1], dense.eden ); 01136 } 01137 01138 return(0); 01139 } 01140 01141 /*sumcon sums L and Q for net incident continuum */ 01142 STATIC void sumcon(long int il, 01143 long int ih, 01144 realnum *q, 01145 realnum *p, 01146 realnum *panu) 01147 { 01148 long int i, 01149 iupper; /* used as upper limit to the sum */ 01150 01151 DEBUG_ENTRY( "sumcon()" ); 01152 01153 *q = 0.; 01154 *p = 0.; 01155 *panu = 0.; 01156 01157 /* soft continua may not go as high as the requested bin */ 01158 iupper = MIN2(rfield.nflux,ih); 01159 01160 /* n.b. - in F77 loop IS NOT executed when IUPPER < IL */ 01161 for( i=il-1; i < iupper; i++ ) 01162 { 01163 /* sum photon number */ 01164 *q += rfield.flux[i]; 01165 /* en1ryd is needed to stop overflow */ 01166 /* sum flux */ 01167 *p += (realnum)(rfield.flux[i]*(rfield.anu[i]*EN1RYD)); 01168 /* this sum needed for means */ 01169 *panu += (realnum)(rfield.flux[i]*(rfield.anu2[i]*EN1RYD)); 01170 } 01171 01172 return; 01173 } 01174 01175 /*ptrcer show continuum pointers in real time following drive pointers command */ 01176 STATIC void ptrcer(void) 01177 { 01178 char chCard[INPUT_LINE_LENGTH]; 01179 /* in case of checking everything, will write errors to this file */ 01180 FILE * ioERRORS=NULL; 01181 bool lgEOL; 01182 char chKey; 01183 long int i, 01184 ipnt, 01185 j; 01186 double pnt, 01187 t1, 01188 t2; 01189 01190 DEBUG_ENTRY( "ptrcer()" ); 01191 01192 fprintf( ioQQQ, " There are two ways to do this:\n"); 01193 fprintf( ioQQQ, " do you want me to test all the pointers (enter y)\n"); 01194 fprintf( ioQQQ, " or do you want to enter energies yourself? (enter n)\n" ); 01195 01196 if( read_whole_line( chCard , (int)sizeof(chCard) , ioStdin ) == NULL ) 01197 { 01198 fprintf( ioQQQ, " error getting input \n" ); 01199 cdEXIT(EXIT_FAILURE); 01200 } 01201 01202 /* this must be either y or n */ 01203 chKey = chCard[0]; 01204 01205 if( chKey == 'n' ) 01206 { 01207 /* this branch, enter energies by hand, and see what happens */ 01208 fprintf( ioQQQ, " Enter energy (Ryd); 0 to stop; negative is log.\n" ); 01209 pnt = 1.; 01210 while( pnt!=0. ) 01211 { 01212 if( read_whole_line( chCard , (int)sizeof(chCard) , ioStdin ) == NULL ) 01213 { 01214 fprintf( ioQQQ, " error getting input2 \n" ); 01215 cdEXIT(EXIT_FAILURE); 01216 } 01217 /* now get the number off the line */ 01218 i = 1; 01219 pnt = FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 01220 01221 /* bail if no number at all, or it is zero*/ 01222 if( lgEOL || pnt==0. ) 01223 { 01224 break; 01225 } 01226 01227 /* if number negative then interpret as log */ 01228 if( pnt < 0. ) 01229 { 01230 pnt = pow(10.,pnt); 01231 } 01232 01233 /* get pointer to call */ 01234 ipnt = ipoint(pnt); 01235 fprintf( ioQQQ, " Cell num%4ld center:%10.2e width:%10.2e low:%10.2e hi:%10.2e convoc:%10.2e\n", 01236 ipnt, rfield.anu[ipnt-1], rfield.widflx[ipnt-1], 01237 rfield.anu[ipnt-1] - rfield.widflx[ipnt-1]/2., 01238 rfield.anu[ipnt-1] + rfield.widflx[ipnt-1]/2., 01239 rfield.convoc[ipnt-1] ); 01240 } 01241 } 01242 01243 else if( chKey == 'y' ) 01244 { 01245 /* first check that ipoint will not crash due to out of range call*/ 01246 if( rfield.anu[0] - rfield.widflx[0]/2.*0.9 < continuum.filbnd[0] ) 01247 { 01248 fprintf( ioQQQ," ipoint would crash since lowest desired energy of %e ryd is below limit of %e\n", 01249 rfield.anu[0] - rfield.widflx[0]/2.*0.9 , continuum.filbnd[0] ); 01250 fprintf( ioQQQ," width of cell is %e\n",rfield.widflx[0]); 01251 cdEXIT(EXIT_FAILURE); 01252 } 01253 01254 else if( rfield.anu[rfield.nflux-1] + rfield.widflx[rfield.nflux-1]/2.*0.9 > 01255 continuum.filbnd[continuum.nrange] ) 01256 { 01257 fprintf( ioQQQ," ipoint would crash since highest desired energy of %e ryd is above limit of %e\n", 01258 rfield.anu[rfield.nflux-1] + rfield.widflx[rfield.nflux-1]/2.*0.9 , 01259 continuum.filbnd[continuum.nrange-1] ); 01260 fprintf( ioQQQ," width of cell is %e\n",rfield.widflx[rfield.nflux]); 01261 fprintf( ioQQQ," this, previous cells are %e %e\n", 01262 rfield.anu[rfield.nflux-1],rfield.anu[rfield.nflux-2]); 01263 cdEXIT(EXIT_FAILURE); 01264 } 01265 01266 /* this branch check everything, write errors to error file */ 01267 fprintf( ioQQQ, " errors output on errors.txt\n"); 01268 fprintf( ioQQQ, " IP(cor),IP(fount),nu lower, upper of found, desired cell.\n" ); 01269 01270 /* error file not open, set to null so we can check later */ 01271 ioERRORS = NULL; 01272 for( i=0; i < rfield.nflux-1; i++ ) 01273 { 01274 t1 = rfield.anu[i] - rfield.widflx[i]/2.*0.9; 01275 t2 = rfield.anu[i] + rfield.widflx[i]/2.*0.9; 01276 01277 j = ipoint(t1); 01278 if( j != i+1 ) 01279 { 01280 /* open file for errors if not already open */ 01281 if( ioERRORS == NULL ) 01282 { 01283 ioERRORS = open_data( "errors.txt", "w", AS_LOCAL_ONLY ); 01284 fprintf( ioQQQ," created errors.txt file with error summary\n"); 01285 } 01286 01287 fprintf( ioQQQ, " Pointers do not agree for lower bound of cell%4ld, %e\n", 01288 i, rfield.anu[i]); 01289 fprintf( ioERRORS, " Pointers do not agree for lower bound of cell%4ld, %e\n", 01290 i, rfield.anu[i] ); 01291 } 01292 01293 j = ipoint(t2); 01294 if( j != i+1 ) 01295 { 01296 /* open file for errors if not already open */ 01297 if( ioERRORS == NULL ) 01298 { 01299 ioERRORS = open_data( "errors.txt", "w", AS_LOCAL_ONLY ); 01300 fprintf( ioQQQ," created errors.txt file with error summary\n"); 01301 } 01302 fprintf( ioQQQ, " Pointers do not agree for upper bound of cell%4ld, %e\n", 01303 i , rfield.anu[i]); 01304 fprintf( ioERRORS, " Pointers do not agree for upper bound of cell%4ld, %e\n", 01305 i , rfield.anu[i]); 01306 } 01307 01308 } 01309 } 01310 01311 else 01312 { 01313 fprintf( ioQQQ, "I do not understand this key, sorry. %c\n", chKey ); 01314 cdEXIT(EXIT_FAILURE); 01315 } 01316 01317 if( ioERRORS!=NULL ) 01318 fclose( ioERRORS ); 01319 cdEXIT(EXIT_FAILURE); 01320 } 01321 01322 /*extin do extinction of incident continuum as set by extinguish command */ 01323 STATIC void extin(realnum *ex1ryd) 01324 { 01325 long int i, 01326 low; 01327 double absorb, 01328 factor; 01329 01330 DEBUG_ENTRY( "extin()" ); 01331 01332 /* modify input continuum by leaky absorber 01333 * power law fit to 01334 * >>refer XUV extinction Cruddace et al. 1974 ApJ 187, 497. */ 01335 if( extinc.excolm == 0. ) 01336 { 01337 *ex1ryd = 1.; 01338 } 01339 else 01340 { 01341 absorb = 1. - extinc.exleak; 01342 /*factor = extinc.excolm*6.22e-18;*/ 01343 /* >>chng 01 dec 19, use variable for the constant that gives extinction */ 01344 factor = extinc.excolm*extinc.cnst_col2optdepth; 01345 /* extinction at 1 and 4 Ryd */ 01346 *ex1ryd = (realnum)(extinc.exleak + absorb*sexp(factor)); 01347 01348 low = ipoint(extinc.exlow); 01349 for( i=low-1; i < rfield.nflux; i++ ) 01350 { 01351 realnum extfactor = (realnum)(extinc.exleak + absorb*sexp(factor* 01352 (pow(rfield.anu[i],extinc.cnst_power)))); 01353 rfield.flux_beam_const[i] *= extfactor; 01354 rfield.flux_beam_time[i] *= extfactor; 01355 rfield.flux_isotropic[i] *= extfactor; 01356 rfield.flux[i] = rfield.flux_beam_const[i] + rfield.flux_beam_time[i] + 01357 rfield.flux_isotropic[i]; 01358 } 01359 } 01360 return; 01361 } 01362 01363 /*conorm normalize continuum to proper intensity */ 01364 STATIC void conorm(void) 01365 { 01366 long int i , nd; 01367 double xLog_radius_inner, 01368 diff, 01369 f, 01370 flx1, 01371 flx2, 01372 pentrd, 01373 qentrd; 01374 01375 DEBUG_ENTRY( "conorm()" ); 01376 01377 xLog_radius_inner = log10(radius.Radius); 01378 01379 /* this is a sanity check, it can't happen */ 01380 for( i=0; i < rfield.nspec; i++ ) 01381 { 01382 if( strcmp(rfield.chRSpec[i],"UNKN") == 0 ) 01383 { 01384 fprintf( ioQQQ, " UNKN spectral normalization cannot happen.\n" ); 01385 fprintf( ioQQQ, " conorm punts.\n" ); 01386 cdEXIT(EXIT_FAILURE); 01387 } 01388 01389 else if( strcmp(rfield.chRSpec[i],"SQCM") != 0 && 01390 strcmp(rfield.chRSpec[i],"4 PI") != 0 ) 01391 { 01392 fprintf( ioQQQ, " chRSpec must be SQCM or 4 PI, and it was %4.4s. This cannot happen.\n", 01393 rfield.chRSpec[i] ); 01394 fprintf( ioQQQ, " conorm punts.\n" ); 01395 cdEXIT(EXIT_FAILURE); 01396 } 01397 01398 01399 /* this sanity check makes sure that atlas.mod or werner.mod grids 01400 * are for the current version of the code */ 01401 if( strcmp(rfield.chSpType[i],"VOLK ") == 0 ) 01402 { 01403 /* check that wavelength scale is actually defined outside here */ 01404 ASSERT( rfield.AnuOrg[rfield.nupper-1]>0. ); 01405 01406 diff = fabs(rfield.tNuRyd[i][rfield.nupper-1]-rfield.AnuOrg[rfield.nupper-1])/ 01407 rfield.AnuOrg[rfield.nupper-1]; 01408 01409 /* this was read from a binary file, so match should be precise */ 01410 if( diff > 10.*FLT_EPSILON ) 01411 { 01412 if( continuum.ResolutionScaleFactor == 1. ) 01413 { 01414 fprintf( ioQQQ, "%10.2e%10.2e\n", rfield.AnuOrg[rfield.nupper-1], 01415 rfield.tNuRyd[i][rfield.nupper-1] ); 01416 01417 fprintf( ioQQQ,"conorm: The energy grid of the stellar atmosphere file does not agree with the grid in this version of the code.\n" ); 01418 fprintf( ioQQQ,"A stellar atmosphere grid from an old version of the code is probably in place.\n" ); 01419 fprintf( ioQQQ,"A grid for the current version of Cloudy must be generated and used.\n" ); 01420 fprintf( ioQQQ,"This is done with the COMPILE STARS command.\n" ); 01421 fprintf( ioQQQ,"Sorry.\n" ); 01422 01423 cdEXIT(EXIT_FAILURE); 01424 } 01425 else 01426 { 01427 fprintf( ioQQQ,"\n\nThe continuum resolution has been chnaged with the SET CONTINUUM RESOLUTION command.\n" ); 01428 fprintf( ioQQQ,"The compiled stellar continua are not consistent with this changed continuum.\n" ); 01429 fprintf( ioQQQ,"Sorry.\n" ); 01430 cdEXIT(EXIT_FAILURE); 01431 } 01432 } 01433 } 01434 } 01435 01436 /* this sanity check is that the grains we have read in from opacity files agree 01437 * with the energy grid in this version of cloudy */ 01438 for( nd=0; nd < gv.nBin; nd++ ) 01439 { 01440 bool lgErrorFound = false; 01441 01442 /* these better agree */ 01443 if( gv.bin[nd]->NFPCheck != rfield.nupper ) 01444 { 01445 fprintf( ioQQQ, 01446 " conorm: expected:%ld found: %ld: number of frequency points in grains data do not match\n", 01447 rfield.nupper, gv.bin[nd]->NFPCheck ); 01448 lgErrorFound = true; 01449 } 01450 01451 /* check that wavelength scale is actually defined outside here */ 01452 ASSERT( rfield.AnuOrg[rfield.nupper-1]>0. ); 01453 01454 diff = fabs( gv.bin[nd]->EnergyCheck - rfield.AnuOrg[rfield.nupper-1] )/ 01455 rfield.AnuOrg[rfield.nupper-1]; 01456 01457 /* this was read from an ascii file, so we have to be more lenient 01458 * the last constant is determined by the number of decimal places in the .opc file */ 01459 if( diff > MAX2(10.*FLT_EPSILON,3.e-6f) ) 01460 { 01461 fprintf( ioQQQ, "%14.6e%14.6e: frequencies of last grid point do not match\n", 01462 rfield.AnuOrg[rfield.nupper-1], gv.bin[nd]->EnergyCheck ); 01463 lgErrorFound = true; 01464 } 01465 01466 if( lgErrorFound ) 01467 { 01468 if( continuum.ResolutionScaleFactor == 1. ) 01469 { 01470 fprintf( ioQQQ,"PROBLEM DISASTER conorm: The energy grid of the grain opacity file does not agree with the grid in this version of the code.\n" ); 01471 fprintf( ioQQQ,"A compiled grain grid from an old version of the code is probably in place.\n" ); 01472 fprintf( ioQQQ,"A grid for the current version of Cloudy must be generated and used.\n" ); 01473 fprintf( ioQQQ,"This is done with the COMPILE ALL GRAINS command.\n" ); 01474 fprintf( ioQQQ,"Sorry.\n" ); 01475 01476 cdEXIT(EXIT_FAILURE); 01477 } 01478 else 01479 { 01480 fprintf( ioQQQ,"\n\nPROBLEM DISASTER The continuum resolution has been chnaged with the SET CONTINUUM RESOLUTION command.\n" ); 01481 fprintf( ioQQQ,"The compiled grain opacities are not consistent with this changed continuum, so cannot be used.\n" ); 01482 fprintf( ioQQQ,"Sorry.\n" ); 01483 01484 cdEXIT(EXIT_FAILURE); 01485 } 01486 } 01487 } 01488 01489 /* default is is to predict line intensities, 01490 * but if any continuum specified as luminosity then would override this - 01491 * following two values are correct for intensities */ 01492 radius.pirsq = 0.; 01493 radius.lgPredLumin = false; 01494 01495 /* check whether ANY luminosities are present */ 01496 for( i=0; i < rfield.nspec; i++ ) 01497 { 01498 if( strcmp(rfield.chRSpec[i],"4 PI") == 0 ) 01499 { 01500 radius.pirsq = (realnum)(1.0992099 + 2.*xLog_radius_inner); 01501 radius.lgPredLumin = true; 01502 /* convert down to intensity */ 01503 rfield.totpow[i] -= radius.pirsq; 01504 01505 if( trace.lgTrace ) 01506 { 01507 fprintf( ioQQQ, 01508 " conorm converts continuum %ld from luminosity to intensity.\n", 01509 i ); 01510 } 01511 } 01512 } 01513 01514 /* if total luminosities are present, must have specified a starting radius */ 01515 if( radius.lgPredLumin && !radius.lgRadiusKnown ) 01516 { 01517 fprintf(ioQQQ,"PROBLEM DISASTER conorm: - A continuum source was specified as a luminosity, but the inner radius of the cloud was not set.\n"); 01518 fprintf(ioQQQ,"Please set an inner radius.\nSorry.\n"); 01519 cdEXIT(EXIT_FAILURE); 01520 } 01521 01522 /* convert ionization parameters to number of photons, called "q(h)" 01523 * at this stage q(h) and "PHI " are the same */ 01524 for( i=0; i < rfield.nspec; i++ ) 01525 { 01526 if( strcmp(rfield.chSpNorm[i],"IONI") == 0 ) 01527 { 01528 /* the log of the ionization parameter was stored here, this converts 01529 * it to the log of the number of photons per sq cm */ 01530 rfield.totpow[i] += log10(dense.gas_phase[ipHYDROGEN]) + log10(SPEEDLIGHT); 01531 strcpy( rfield.chSpNorm[i], "Q(H)" ); 01532 if( trace.lgTrace ) 01533 { 01534 fprintf( ioQQQ, 01535 " conorm converts continuum %ld from ionizat par to q(h).\n", 01536 i ); 01537 } 01538 } 01539 } 01540 01541 /* convert x-ray ionization parameter xi to intensity */ 01542 for( i=0; i < rfield.nspec; i++ ) 01543 { 01544 if( strcmp(rfield.chSpNorm[i],"IONX") == 0 ) 01545 { 01546 /* this converts it to an intensity */ 01547 rfield.totpow[i] += log10(dense.gas_phase[ipHYDROGEN]) - log10(PI4); 01548 strcpy( rfield.chSpNorm[i], "LUMI" ); 01549 if( trace.lgTrace ) 01550 { 01551 fprintf( ioQQQ, " conorm converts continuum%3ld from x-ray ionizat par to I.\n", 01552 i ); 01553 } 01554 } 01555 } 01556 01557 /* indicate whether we ended up with luminosity or intensity */ 01558 if( trace.lgTrace ) 01559 { 01560 if( radius.lgPredLumin ) 01561 { 01562 fprintf( ioQQQ, " Cloudy will predict lumin into 4pi\n" ); 01563 } 01564 else 01565 { 01566 fprintf( ioQQQ, " Cloudy will do surface flux for lumin\n" ); 01567 } 01568 } 01569 01570 /* if intensity per unit area is predicted then geometric 01571 * covering factor must be unity 01572 * variable can also be set elsewhere */ 01573 if( !radius.lgPredLumin ) 01574 { 01575 geometry.covgeo = 1.; 01576 } 01577 01578 /* main loop over all continuum shapes to find continuum normalization 01579 * for each one */ 01580 for( i=0; i < rfield.nspec; i++ ) 01581 { 01582 rfield.ipspec = i; 01583 01584 /* check that, if laser, bounds include laser energy */ 01585 if( strcmp(rfield.chSpType[rfield.ipspec],"LASER") == 0 ) 01586 { 01587 if( !( rfield.range[rfield.ipspec][0] < rfield.slope[rfield.ipspec] && 01588 rfield.range[rfield.ipspec][1] > rfield.slope[rfield.ipspec]) ) 01589 { 01590 fprintf(ioQQQ,"PROBLEM DISASTER, a continuum source is a laser at %f Ryd, but the intensity was specified over a range from %f to %f Ryd.\n", 01591 rfield.slope[rfield.ipspec], 01592 rfield.range[rfield.ipspec][0], 01593 rfield.range[rfield.ipspec][1]); 01594 fprintf(ioQQQ,"Please specify the continuum flux where the laser is active.\n"); 01595 cdEXIT(EXIT_FAILURE); 01596 } 01597 } 01598 01599 if( trace.lgTrace ) 01600 { 01601 long int jj; 01602 fprintf( ioQQQ, " conorm continuum number %ld is shape %s range is %.2e %.2e\n", 01603 i, 01604 rfield.chSpType[i], 01605 rfield.range[i][0], 01606 rfield.range[i][1] ); 01607 fprintf( ioQQQ, "the continuum points follow\n"); 01608 jj = 0; 01609 if( rfield.lgContMalloc[rfield.ipspec] ) 01610 { 01611 while( rfield.tNuRyd[rfield.ipspec][jj] != 0. && jj < 100 ) 01612 { 01613 fprintf( ioQQQ, "%li %e %e\n", 01614 jj , 01615 rfield.tNuRyd[rfield.ipspec][jj], 01616 rfield.tslop[rfield.ipspec][jj] ); 01617 ++jj; 01618 } 01619 } 01620 } 01621 01622 if( strcmp(rfield.chSpNorm[i],"RATI") == 0 ) 01623 { 01624 /* option to scale relative to previous continua 01625 * this must come first since otherwise may be too late 01626 * BUT ratio cannot be the first continuum source */ 01627 if( trace.lgTrace ) 01628 { 01629 fprintf( ioQQQ, " conorm this is ratio to 1st con\n" ); 01630 } 01631 01632 /* check that this is not the first continuum source, we must ratio */ 01633 if( i == 0 ) 01634 { 01635 fprintf( ioQQQ, " I cant form a ratio if continuum is first source.\n" ); 01636 cdEXIT(EXIT_FAILURE); 01637 } 01638 01639 /* first find photon flux and Q of prevous continuum */ 01640 rfield.ipspec -= 1; 01641 flx1 = ffun1(rfield.range[i][0])*rfield.spfac[rfield.ipspec]* 01642 rfield.range[i][0]; 01643 01644 /* check that previous continua were not zero where ratio is formed */ 01645 if( flx1 <= 0. ) 01646 { 01647 fprintf( ioQQQ, " Previous continua were zero where ratio is desired.\n" ); 01648 cdEXIT(EXIT_FAILURE); 01649 } 01650 01651 /* return pointer to previous (correct) value, find F, Q */ 01652 rfield.ipspec += 1; 01653 01654 /* we want a continuum totpow as powerful, flx is now desired flx */ 01655 flx1 *= rfield.totpow[i]; 01656 01657 /*. find flux of this new continuum at that point */ 01658 flx2 = ffun1(rfield.range[i][1])*rfield.range[i][1]; 01659 01660 /* this is ratio of desired to actual */ 01661 rfield.spfac[i] = flx1/flx2; 01662 if( trace.lgTrace ) 01663 { 01664 fprintf( ioQQQ, " conorm ratio will set scale fac to%10.3e at%10.2e Ryd.\n", 01665 rfield.totpow[i], rfield.range[i][0] ); 01666 } 01667 } 01668 01669 else if( strcmp(rfield.chSpNorm[i],"FLUX") == 0 ) 01670 { 01671 /* specify flux density 01672 * option to use arbitrary frequency or range */ 01673 f = ffun1(rfield.range[i][0]); 01674 01675 /* make sure this is positive, could be zero if out of range of table, 01676 * or for various forms of insanity */ 01677 if( f<=0. ) 01678 { 01679 fprintf( ioQQQ, "\n\n PROBLEM DISASTER\n The intensity of continuum source %ld is non-positive at the energy used to normalize it (%.3e Ryd). Something is seriously wrong.\n", 01680 i , rfield.range[i][0]); 01681 /* is this a table? if so, is ffun within its bounds? */ 01682 if( strcmp(rfield.chSpType[i],"INTER") == 0 ) 01683 fprintf( ioQQQ, " This continuum shape given by a table of points - check that the intensity is specified at an energy within the range of that table.\n"); 01684 fprintf( ioQQQ, " Also check that the numbers used to specify the shape and intensity do not under or overflow on this cpu.\n\n"); 01685 01686 cdEXIT(EXIT_FAILURE); 01687 } 01688 01689 /* now convert to log in continuum units we shall use */ 01690 f = MAX2(1e-37,f); 01691 f = log10(f) + log10(rfield.range[i][0]*EN1RYD/FR1RYD); 01692 01693 f = rfield.totpow[i] - f; 01694 /* >>chng 96 dec 31, added following test */ 01695 if( f > 35. ) 01696 { 01697 fprintf( ioQQQ, " PROBLEM DISASTER Continuum source %ld is too intense for this cpu - is it normalized correctly?\n", 01698 i ); 01699 cdEXIT(EXIT_FAILURE); 01700 } 01701 01702 rfield.spfac[i] = pow(10.,f); 01703 if( trace.lgTrace ) 01704 { 01705 fprintf( ioQQQ, " conorm will set log fnu to%10.3e at%10.2e Ryd. Factor=%11.4e\n", 01706 rfield.totpow[i], rfield.range[i][0], rfield.spfac[i] ); 01707 } 01708 } 01709 01710 else if( strcmp(rfield.chSpNorm[i],"Q(H)") == 0 || 01711 strcmp(rfield.chSpNorm[i],"PHI ") == 0 ) 01712 { 01713 /* some type of photon density entered */ 01714 if( trace.lgTrace ) 01715 { 01716 fprintf( ioQQQ, " conorm calling qintr range=%11.3e %11.3e desired val is %11.3e\n", 01717 rfield.range[i][0], 01718 rfield.range[i][1] , 01719 rfield.totpow[i]); 01720 } 01721 01722 /* the total number of photons over the specified range in 01723 * the arbitrary system of units that the code save the continuum shape */ 01724 qentrd = qintr(&rfield.range[i][0],&rfield.range[i][1]); 01725 /* this is the log of the scale factor that must multiply the 01726 * continuum shape to get the final set of numbers */ 01727 diff = rfield.totpow[i] - qentrd; 01728 01729 /* >>chng 03 mar 13, from diff < -25 to <-35, 01730 * tripped for very low U models used for H2 simulations */ 01731 /*if( diff < -25. || diff > 35. )*/ 01732 if( diff < -35. || diff > 35. ) 01733 { 01734 fprintf( ioQQQ, " PROBLEM DISASTER Continuum source specified is too extreme.\n" ); 01735 fprintf( ioQQQ, 01736 " The integral over the continuum shape gave (log) %.3e photons, and the command requested (log) %.3e.\n" , 01737 qentrd , rfield.totpow[i]); 01738 fprintf( ioQQQ, 01739 " The difference in the log is %.3e.\n" , 01740 diff ); 01741 if( diff>0. ) 01742 { 01743 fprintf( ioQQQ, " The continuum source is too bright.\n" ); 01744 } 01745 else 01746 { 01747 fprintf( ioQQQ, " The continuum source is too faint.\n" ); 01748 } 01749 /* explain what happened */ 01750 fprintf( ioQQQ, " The usual cause for this problem is an incorrect continuum intensity/luminosity or radius command.\n" ); 01751 fprintf( ioQQQ, " There were a total of %li continuum shape commands entered - the problem is with number %li.\n", 01752 rfield.nspec , i+1 ); 01753 cdEXIT(EXIT_FAILURE); 01754 } 01755 01756 else 01757 { 01758 rfield.spfac[i] = pow(10.,diff); 01759 } 01760 01761 if( trace.lgTrace ) 01762 { 01763 fprintf( ioQQQ, " conorm finds Q over range from%11.4e-%11.4e Ryd, integral= %10.4e Factor=%11.4e\n", 01764 rfield.range[i][0], 01765 rfield.range[i][1], 01766 qentrd , 01767 rfield.spfac[i] ); 01768 } 01769 } 01770 01771 else if( strcmp(rfield.chSpNorm[i],"LUMI") == 0 ) 01772 { 01773 /* luminosity entered, special since default is TOTAL lumin */ 01774 /*pintr integrates L for any continuum between two limits, used for normalization, 01775 * return units are log of ryd cm-2 s-1, last log conv to ergs */ 01776 pentrd = pintr(rfield.range[i][0],rfield.range[i][1]) + 01777 log10(EN1RYD); 01778 f = rfield.totpow[i] - pentrd; 01779 01780 /* >>chng 96 dec 31, added following test */ 01781 if( f > 35. ) 01782 { 01783 fprintf( ioQQQ, " PROBLEM DISASTER Continuum source %ld is too intense for this cpu - is it normalized correctly?\n", 01784 i ); 01785 cdEXIT(EXIT_FAILURE); 01786 } 01787 01788 rfield.spfac[i] = pow(10.,f); 01789 if( trace.lgTrace ) 01790 { 01791 fprintf( ioQQQ, " conorm finds luminosity range is %10.3e to %9.3e Ryd, factor is %11.4e\n", 01792 rfield.range[i][0], rfield.range[i][1], 01793 rfield.spfac[i] ); 01794 } 01795 } 01796 01797 else 01798 { 01799 fprintf( ioQQQ, "PROBLEM DISASTER What chSpNorm label is this? =%s=\n", rfield.chSpNorm[i]); 01800 TotalInsanity(); 01801 } 01802 01803 /* sec part after .or. added June 93 because sometimes spfac=0 01804 * got past first test */ 01805 if( 1./rfield.spfac[i] == 0. || rfield.spfac[i] == 0. ) 01806 { 01807 fprintf( ioQQQ, "PROBLEM DISASTER conorm finds infinite continuum scale factor.\n" ); 01808 fprintf( ioQQQ, "The continuum is too intense to compute with this cpu.\n" ); 01809 fprintf( ioQQQ, "Were the intensity and luminosity commands switched?\n" ); 01810 fprintf( ioQQQ, "Sorry, but I cannot go on.\n" ); 01811 cdEXIT(EXIT_FAILURE); 01812 } 01813 } 01814 01815 /* this is conversion factor for final units of line intensities or luminosities in printout, 01816 * will be intensities (==0) unless luminosity is to be printed, or flux at Earth 01817 * pirsq is the log of 4 pi r_in^2 */ 01818 radius.Conv2PrtInten = radius.pirsq; 01819 01820 /* >>chng 02 apr 25, add option for slit on aperture command */ 01821 if( geometry.iEmissPower == 1 ) 01822 { 01823 if( radius.lgPredLumin ) 01824 { 01825 /* factor should be divided by 2 r_in */ 01826 radius.Conv2PrtInten -= (log10(2.) + xLog_radius_inner); 01827 } 01828 else if( !radius.lgPredLumin ) 01829 { 01830 /* this is an error - slit requested but radius is not known */ 01831 fprintf( ioQQQ, "PROBLEM DISASTER conorm: Aperture slit specified, but not predicting luminosity.\n" ); 01832 fprintf( ioQQQ, "conorm: Please specify an inner radius to determine L.\nSorry\n" ); 01833 cdEXIT(EXIT_FAILURE); 01834 } 01835 } 01836 if( geometry.iEmissPower == 0 && radius.lgPredLumin) 01837 { 01838 /* leave Conv2PrtInten at zero if not predicting luminosity */ 01839 radius.Conv2PrtInten = log10(2.); 01840 } 01841 01842 /* this is option go give final absolute results as flux observed at Earth */ 01843 if( radius.distance > 0. && radius.lgRadiusKnown && prt.lgPrintFluxEarth ) 01844 { 01845 /*radius.Conv2PrtInten -= 2.*xLog_radius_inner - 2.*log10(radius.distance);*/ 01846 /* div (in log) by 4 pi dist^2 */ 01847 radius.Conv2PrtInten -= log10( 4.*PI*POW2(radius.distance) ); 01848 } 01849 01850 /* normally lines are into 4pi, this is option to do per sr or arcsec^2 */ 01851 if( prt.lgSurfaceBrightness ) 01852 { 01853 if( radius.pirsq!= 0. ) 01854 { 01855 /* make sure we are predicting line intensities, not luminosity */ 01856 fprintf( ioQQQ, " PROBLEM DISASTER Sorry, but both luminosity and surface brightness have been requested for lines.\n" ); 01857 fprintf( ioQQQ, " the PRINT LINE SURFACE BRIGHTNESS command can only be used when lines are predicted per unit cloud area.\n" ); 01858 cdEXIT(EXIT_FAILURE); 01859 } 01860 if( prt.lgSurfaceBrightness_SR ) 01861 { 01862 /* we want final units to be per sr */ 01863 radius.Conv2PrtInten -= log10( PI4 ); 01864 } 01865 else 01866 { 01867 /* we want final units to be per square arcsec, 01868 * first 4 pi converts to per sr, 01869 * there are 4pi sr in 360 deg, so term in square 01870 * goes from sr to square sec of arc */ 01871 radius.Conv2PrtInten -= log10( PI4 *POW2( (360./(PI2)*3600.) ) ); 01872 } 01873 } 01874 return; 01875 } 01876 01877 /*qintr integrates Q for any continuum between two limits, used for normalization */ 01878 STATIC double qintr(double *qenlo, 01879 double *qenhi) 01880 { 01881 long int i, 01882 ipHi, 01883 ipLo, 01884 j; 01885 double qintr_v, 01886 sum, 01887 wanu; 01888 01889 DEBUG_ENTRY( "qintr()" ); 01890 01891 /* returns LOG of number of photons over energy interval */ 01892 sum = 0.; 01893 /* >>chng 02 oct 27, do not use qg32, always use same method as 01894 * routine that does final set of continuum */ 01895 01896 /* this is copy of logic that occurs three times across code */ 01897 ipLo = ipoint(*qenlo); 01898 ipHi = ipoint(*qenhi); 01899 /* this is actual sum of photons within band */ 01900 for( i=ipLo-1; i < (ipHi - 1); i++ ) 01901 { 01902 /*sum += ffun1(rfield.anu[i])*rfield.widflx[i];*/ 01903 for( j=0; j < 4; j++ ) 01904 { 01905 wanu = rfield.anu[i] + rfield.widflx[i]*aweigh[j]; 01906 /* >>chng 02 jul 16, add test on continuum limits - 01907 * this was exceeded when resolution set very large */ 01908 wanu = MAX2( wanu , rfield.emm ); 01909 wanu = MIN2( wanu , rfield.egamry ); 01910 sum += fweigh[j]*ffun1(wanu)*rfield.widflx[i]; 01911 } 01912 } 01913 01914 if( sum <= 0. ) 01915 { 01916 fprintf( ioQQQ, " PROBLEM DISASTER Photon number sum in QINTR is %.3e\n", 01917 sum ); 01918 fprintf( ioQQQ, " This source has no ionizing radiation, and the number of ionizing photons was specified.\n" ); 01919 fprintf( ioQQQ, " This was continuum source number%3ld\n", 01920 rfield.ipspec ); 01921 fprintf( ioQQQ, " Sorry, but I cannot go on. ANU and FLUX arrays follow. Enjoy.\n" ); 01922 for( i=0; i < rfield.nupper; i++ ) 01923 { 01924 fprintf( ioQQQ, "%.2e\t%.2e\n", 01925 rfield.anu[i], 01926 rfield.flux[i] ); 01927 } 01928 cdEXIT(EXIT_FAILURE); 01929 } 01930 else 01931 { 01932 qintr_v = log10(sum); 01933 } 01934 return( qintr_v ); 01935 } 01936 01937 /*pintr integrates L for any continuum between two limits, used for normalization, 01938 * return units are log of ryd cm-2 s-1 */ 01939 STATIC double pintr(double penlo, 01940 double penhi) 01941 { 01942 long int i, 01943 j; 01944 double fsum, 01945 pintr_v, 01946 sum, 01947 wanu, 01948 wfun; 01949 long int ip1 , ip2; 01950 01951 DEBUG_ENTRY( "pintr()" ); 01952 01953 /* computes log of luminosity in radiation over some integral 01954 * answer is in Ryd per sec */ 01955 01956 sum = 0.; 01957 /* >>chng 02 oct 27, do not call qg32, do same type sum as 01958 * final integration */ 01959 /* laser is special since delta function, this is center of laser */ 01960 /* >>chng 01 jul 01, was +-21 cells, change to call to ipoint */ 01961 ip1 = ipoint( penlo ); 01962 01963 ip2 = ipoint( penhi ); 01964 01965 for( i=ip1-1; i < ip2-1; i++ ) 01966 { 01967 fsum = 0.; 01968 for( j=0; j < 4; j++ ) 01969 { 01970 wanu = rfield.anu[i] + rfield.widflx[i]*aweigh[j]; 01971 /*++iiii;fprintf(ioQQQ,"DEBUG iii %li %e \n",iiii, wanu );*/ 01972 wfun = fweigh[j]*ffun1(wanu)*wanu; 01973 fsum += wfun; 01974 } 01975 sum += fsum*rfield.widflx[i]; 01976 } 01977 01978 if( sum > 0. ) 01979 { 01980 pintr_v = log10(sum); 01981 } 01982 else 01983 { 01984 pintr_v = -38.; 01985 } 01986 01987 return( pintr_v ); 01988 }