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 /* GammaBn evaluate photoionization rate for single shell with induced recomb */ 00004 /* GammaBnPL evaluate photoionization rate for single shell with induced recomb */ 00005 /* GammaPrt special version of gamma function to print strong contributors */ 00006 /* GammaK evaluate photoionization rate for single shell */ 00007 /* GammaPL evaluate photoionization rate for power law photo cross section */ 00008 /* GammaPrtRate print photo rates for all shells of a ion and element */ 00009 /*GammaPrtShells for the element nelem and ion, print total photo rate, subshells, 00010 * and call GamaPrt for important subshells */ 00011 #include "cddefines.h" 00012 #include "physconst.h" 00013 #include "iso.h" 00014 #include "thermal.h" 00015 #include "secondaries.h" 00016 #include "opacity.h" 00017 #include "rfield.h" 00018 #include "ionbal.h" 00019 #include "atmdat.h" 00020 #include "heavy.h" 00021 #include "gammas.h" 00022 00023 /* 00024 * these are the routines that evaluate the photoionization rates, gammas, 00025 * throughout cloudy. a considerable amount of time is spent in these routines, 00026 * so they should be compiled at the highest possible efficientcy. 00027 */ 00028 00029 /*>>chng 99 apr 16, ConInterOut had not been included in the threshold point for any 00030 * of these routines. added it. also moved thresholds above loop for a few */ 00031 00032 double GammaBn( 00033 /* index for low energy */ 00034 long int ipLoEnr, 00035 /* index for high energy */ 00036 long int ipHiEnr, 00037 /* offset within the opacity stack */ 00038 long int ipOpac, 00039 /* threshold (Ryd) for arbitrary photoionization */ 00040 double thresh, 00041 /* induced rec rate */ 00042 double *ainduc, 00043 /* rec cooling */ 00044 double *rcool) 00045 { 00046 long int i, 00047 ilo, 00048 iup, 00049 limit; 00050 double GamHi, 00051 bnfun_v, 00052 emin, 00053 g, 00054 phisig, 00055 RateInducRec, 00056 RateInducRecCool, 00057 prod; 00058 00059 DEBUG_ENTRY( "GammaBn()" ); 00060 00061 /********************************************************************** 00062 * 00063 * special version of GAMFUN for arbitrary threshold, and induc rec 00064 * compute ainduc, rate of inducted recombinations, rcool, induc rec cool 00065 * 00066 **********************************************************************/ 00067 00068 /* special version of GAMFUN for arbitrary threshold, and induc rec 00069 * compute ainduc, rate of inducted recombinations, rcool, induc rec cool */ 00070 if( ipLoEnr >= rfield.nflux || ipLoEnr >= ipHiEnr ) 00071 { 00072 bnfun_v = 0.; 00073 thermal.HeatNet = 0.; 00074 thermal.HeatLowEnr = 0.; 00075 thermal.HeatHiEnr = 0.; 00076 *ainduc = 0.; 00077 *rcool = 0.; 00078 return( bnfun_v ); 00079 } 00080 00081 ASSERT( ipLoEnr >= 0 && ipHiEnr >= 0 ); 00082 00083 /* >>chng 96 oct 25, first photo elec energy may have been negative 00084 * possible for first any point to be below threshold due to 00085 * finite resolution of continuum mesh */ 00086 emin = MAX2(thresh,rfield.anu[ipLoEnr-1]); 00087 /* >>chng 99 jun 08 heating systematically too small, change to correct 00088 * threshold, and protect first cell */ 00089 emin = thresh; 00090 00091 thermal.HeatNet = 0.; 00092 g = 0.; 00093 RateInducRec = 0.; 00094 RateInducRecCool = 0.; 00095 00096 /* do first point without otscon, which may have diffuse cont, 00097 * extra minus one after ipOpac is correct since not present in i = */ 00098 i = ipLoEnr; 00099 /* >>chng 99 apr 16, add ConInterOut */ 00100 /* >>chng 01 jul 04, add *rfield.lgOutOnly logic */ 00101 phisig = (rfield.flux[i-1] + rfield.otslin[i-1] + 00102 rfield.ConInterOut[i-1]*rfield.lgOutOnly)* 00103 opac.OpacStack[i-ipLoEnr+ipOpac-1]; 00104 00105 g += phisig; 00106 thermal.HeatNet += phisig*rfield.anu[i-1]; 00107 00108 /* integral part of induced rec rate */ 00109 prod = phisig*rfield.ContBoltz[i-1]; 00110 RateInducRec += prod; 00111 /* incuded rec cooling */ 00112 RateInducRecCool += prod*(rfield.anu[i-1] - emin); 00113 00114 iup = MIN2(ipHiEnr,rfield.nflux); 00115 limit = MIN2(iup,secondaries.ipSecIon-1); 00116 00117 /* these is no -1 after the ipLoEnr in the following, due to c off by one counting */ 00118 for( i=ipLoEnr; i < limit; i++ ) 00119 { 00120 /* SummedCon defined in RT_OTS_Update, 00121 * includes flux, otscon, otslin, ConInterOut, outlin */ 00122 phisig = rfield.SummedCon[i]*opac.OpacStack[i-ipLoEnr+ipOpac]; 00123 00124 g += phisig; 00125 thermal.HeatNet += phisig*rfield.anu[i]; 00126 00127 /* phi-sigma times exp(-hnu/kT) */ 00128 prod = phisig*rfield.ContBoltz[i]; 00129 /* induced recombination rate */ 00130 RateInducRec += prod; 00131 /* incuded rec cooling */ 00132 RateInducRecCool += prod*(rfield.anu[i] - emin); 00133 } 00134 00135 /* heating in Rydbergs, no secondary ionization */ 00136 /* >>chng 02 mar 31, added MAX2 due to roundoff error 00137 * creating slightly negative number, caught downstream as insanity */ 00138 thermal.HeatNet = MAX2(0.,thermal.HeatNet - g*thresh); 00139 00140 /* we will save this heat - the part that does not secondary ionize */ 00141 thermal.HeatLowEnr = thermal.HeatNet; 00142 00143 /* now do part where secondary ionization is possible */ 00144 thermal.HeatHiEnr = 0.; 00145 GamHi = 0.; 00146 00147 /* make sure we don't count twice when ipSecIon = ipLoEnr */ 00148 ilo = MAX2(ipLoEnr+1,secondaries.ipSecIon); 00149 for( i=ilo-1; i < iup; i++ ) 00150 { 00151 phisig = rfield.SummedCon[i]*opac.OpacStack[i-ipLoEnr+ipOpac]; 00152 00153 GamHi += phisig; 00154 thermal.HeatHiEnr += phisig*rfield.anu[i]; 00155 00156 /* integral part of induced rec rate */ 00157 prod = phisig*rfield.ContBoltz[i]; 00158 RateInducRec += prod; 00159 /* incuded rec cooling */ 00160 RateInducRecCool += prod*(rfield.anu[i] - emin); 00161 } 00162 00163 /* this is total photo rate, low and high energy parts */ 00164 bnfun_v = g + GamHi; 00165 00166 /* heating in Rydbergs, uncorrected for secondary ionization */ 00167 thermal.HeatHiEnr = thermal.HeatHiEnr - GamHi*thresh; 00168 00169 /* add corrected high energy heat to total */ 00170 thermal.HeatNet += thermal.HeatHiEnr*secondaries.HeatEfficPrimary; 00171 00172 /* final product is heating in erg */ 00173 thermal.HeatNet *= EN1RYD; 00174 thermal.HeatHiEnr *= EN1RYD; 00175 thermal.HeatLowEnr *= EN1RYD; 00176 00177 /* this is an option to turn off induced processes */ 00178 if( rfield.lgInducProcess ) 00179 { 00180 *rcool = RateInducRecCool*EN1RYD; 00181 *ainduc = RateInducRec; 00182 } 00183 else 00184 { 00185 /* the "no induced" command was given */ 00186 *rcool = 0.; 00187 *ainduc = 0.; 00188 } 00189 00190 ASSERT( bnfun_v >= 0. ); 00191 ASSERT( thermal.HeatNet>= 0. ); 00192 return( bnfun_v ); 00193 } 00194 00195 /*GammaPrtShells for the element nelem and ion, print total photo rate, subshells, 00196 * and call GamaPrt for important subshells */ 00197 void GammaPrtShells( long nelem , long ion ) 00198 { 00199 double sum; 00200 long int ns; 00201 00202 DEBUG_ENTRY( "GammaPrtShells()" ); 00203 00204 fprintf(ioQQQ," GammaPrtShells nz\t%.2f \t%.2li %.2li ", 00205 fnzone , 00206 nelem, 00207 ion 00208 ); 00209 sum = 0; 00210 for( ns=0; ns < Heavy.nsShells[nelem][ion]; ++ns ) 00211 { 00212 sum += ionbal.PhotoRate_Shell[nelem][ion][ns][0]; 00213 } 00214 fprintf(ioQQQ,"\ttot\t%.2e", sum); 00215 /* loop over all shells for this ion */ 00216 for( ns=0; ns < Heavy.nsShells[nelem][ion]; ++ns ) 00217 { 00218 fprintf(ioQQQ,"\t%.2e", ionbal.PhotoRate_Shell[nelem][ion][ns][0] ); 00219 # if 0 00220 /* always reevaluate the outer shell, and all shells if lgRedoStatic is set */ 00221 if( (ns==(Heavy.nsShells[nelem][ion]-1) || opac.lgRedoStatic) ) 00222 { 00223 /* option to redo the rates only on occasion */ 00224 iplow = opac.ipElement[nelem][ion][ns][0]; 00225 iphi = opac.ipElement[nelem][ion][ns][1]; 00226 ipop = opac.ipElement[nelem][ion][ns][2]; 00227 00228 /* compute the photoionization rate, ionbal.lgPhotoIoniz_On is 1, set 0 00229 * with "no photoionization" command */ 00230 ionbal.PhotoRate_Shell[nelem][ion][ns][0] = 00231 GammaK(iplow,iphi, 00232 ipop,t_yield::Inst().elec_eject_frac(nelem,ion,ns,0))*ionbal.lgPhotoIoniz_On; 00233 } 00234 # endif 00235 } 00236 fprintf(ioQQQ,"\n"); 00237 return; 00238 } 00239 00240 /********************************************************************** 00241 * 00242 * GammaPrt special version of gamma function to print strong contributors 00243 * this only prints info - does not update heating rates properly since 00244 * this is only a diagnostic 00245 * 00246 **********************************************************************/ 00247 00248 void GammaPrt( 00249 /* first three par are identical to GammaK */ 00250 long int ipLoEnr, 00251 long int ipHiEnr, 00252 long int ipOpac, 00253 /* io unit we will write to */ 00254 FILE * ioFILE, 00255 /* total photo rate from previous call to GammaK */ 00256 double total, 00257 /* we will print contributors that are more than this rate */ 00258 double threshold) 00259 { 00260 long int i, 00261 j, 00262 k; 00263 double flxcor, 00264 phisig; 00265 00266 DEBUG_ENTRY( "GammaPrt()" ); 00267 00268 /* special special version of GAMFUN to punch step-by-step results */ 00269 /* >>chng 99 apr 02, test for lower greater than higher (when shell 00270 * does not exist) added. This caused incorrect photo rates for 00271 * non-existant shells */ 00272 if( ipLoEnr >= rfield.nflux || ipLoEnr >= ipHiEnr ) 00273 { 00274 return; 00275 } 00276 00277 fprintf( ioFILE, " GammaPrt %.2f from ",fnzone); 00278 fprintf( ioFILE,PrintEfmt("%9.2e",rfield.anu[ipLoEnr-1])); 00279 fprintf( ioFILE, " to "); 00280 fprintf( ioFILE,PrintEfmt("%9.2e",rfield.anu[ipHiEnr-1])); 00281 fprintf( ioFILE, "R rates >"); 00282 fprintf( ioFILE,PrintEfmt("%9.2e",threshold)); 00283 fprintf( ioFILE, " of total="); 00284 fprintf( ioFILE,PrintEfmt("%9.2e",total)); 00285 fprintf( ioFILE, " (frac inc, otslin, otscon, ConInterOut, outlin ConOTS_local_OTS_rate ) chL, C\n"); 00286 00287 if( threshold <= 0. || total <= 0. ) 00288 { 00289 return; 00290 } 00291 00292 k = ipLoEnr; 00293 j = MIN2(ipHiEnr,rfield.nflux); 00294 00295 /* do theshold special, do not pick up otscon */ 00296 i = k-1; 00297 00298 /* >>chng 01 jul 04, add *rfield.lgOutOnly logic */ 00299 phisig = (rfield.flux[i] + rfield.otslin[i]+ rfield.ConInterOut[i]*rfield.lgOutOnly)* 00300 opac.OpacStack[i-k+ipOpac]; 00301 if( phisig > threshold || phisig < 0.) 00302 { 00303 flxcor = rfield.flux[i] + rfield.otslin[i] + rfield.ConInterOut[i]*rfield.lgOutOnly; 00304 00305 /* this really is array index on C scale */ 00306 fprintf( ioFILE, "[%5ld]" , i ); 00307 fprintf( ioFILE, PrintEfmt("%9.2e",rfield.anu[i])); 00308 fprintf( ioFILE, PrintEfmt("%9.2e",phisig/total)); 00309 fprintf( ioFILE, "%5.2f%5.2f%5.2f%5.2f%5.2f%5.2f %4.4s %4.4s %.2e \n", 00310 rfield.flux[i]/SDIV(flxcor), 00311 rfield.otslin[i]/SDIV(flxcor), 00312 /* otscon will appear below, but is not counted here, so do not print it (deceiving) */ 00313 0./SDIV(flxcor), 00314 rfield.ConInterOut[i]*rfield.lgOutOnly/SDIV(flxcor), 00315 (rfield.outlin[i]+rfield.outlin_noplot[i])/SDIV(flxcor), 00316 rfield.ConOTS_local_OTS_rate[i]/SDIV(flxcor), 00317 rfield.chLineLabel[i], 00318 rfield.chContLabel[i], 00319 opac.OpacStack[i-k+ipOpac]); 00320 } 00321 00322 for( i=k; i < j; i++ ) 00323 { 00324 phisig = rfield.SummedCon[i]*opac.OpacStack[i-k+ipOpac]; 00325 if( phisig > threshold || phisig < 0.) 00326 { 00327 /* >>chng 01 jul 04, add *rfield.lgOutOnly logic */ 00328 flxcor = rfield.flux[i] + rfield.otslin[i] + rfield.otscon[i] + 00329 rfield.outlin[i] + rfield.outlin_noplot[i] +rfield.ConInterOut[i]*rfield.lgOutOnly; 00330 00331 fprintf( ioFILE, "%5ld", i ); 00332 fprintf(ioFILE,PrintEfmt("%9.2e",rfield.anu[i])); 00333 fprintf(ioFILE,PrintEfmt("%9.2e",phisig/total)); 00334 fprintf( ioFILE, "%5.2f%5.2f%5.2f%5.2f%5.2f%5.2f %4.4s %4.4s %.2e \n", 00335 rfield.flux[i]/SDIV(flxcor), 00336 rfield.otslin[i]/SDIV(flxcor), 00337 rfield.otscon[i]/SDIV(flxcor), 00338 rfield.ConInterOut[i]*rfield.lgOutOnly/SDIV(flxcor), 00339 (rfield.outlin[i] + rfield.outlin_noplot[i])/SDIV(flxcor), 00340 rfield.ConOTS_local_OTS_rate[i]/SDIV(flxcor), 00341 rfield.chLineLabel[i], 00342 rfield.chContLabel[i], 00343 opac.OpacStack[i-k+ipOpac]); 00344 } 00345 } 00346 return; 00347 } 00348 00349 /*GammaK evaluate photoionization rate for single shell */ 00350 00351 /* this routine is a major pace setter for the code 00352 * carefully check anything put into the following do loop */ 00353 00354 double GammaK( 00355 /* ipLoEnr and ipHiEnr are pointers within frequency array to upper and 00356 * lower bounds of evaluation */ 00357 long int ipLoEnr, 00358 long int ipHiEnr, 00359 /* ipOpac is offset to map onto OPSV*/ 00360 long int ipOpac, 00361 /* yield1 is fraction of ionizations that emit 1 electron only, 00362 * only used for energy balance */ 00363 double yield1) 00364 { 00365 long int i, 00366 ilo, 00367 ipOffset, 00368 iup, 00369 limit; 00370 double GamHi, 00371 eauger; 00372 00373 double gamk_v , 00374 phisig; 00375 00376 DEBUG_ENTRY( "GammaK()" ); 00377 00378 /* evaluate photoioinzation rate and photoelectric heating 00379 * returns photoionization rate as GAMK, heating is H */ 00380 00381 /* >>chng 99 apr 02, test for lower greater than higher (when shell 00382 * does not exist) added. This caused incorrect photo rates for 00383 * non-existant shells */ 00384 if( ipLoEnr >= rfield.nflux || ipLoEnr >= ipHiEnr) 00385 { 00386 gamk_v = 0.; 00387 thermal.HeatNet = 0.; 00388 thermal.HeatHiEnr = 0.; 00389 thermal.HeatLowEnr = 0.; 00390 return( gamk_v ); 00391 } 00392 00393 iup = MIN2(ipHiEnr,rfield.nflux); 00394 00395 /* anu(i) is threshold, assume each auger electron has this energy 00396 * less threshold energy, IP lost in initial photoionizaiton 00397 * yield1 is the fraction of photos that emit 1 electron */ 00398 eauger = rfield.anu[ipLoEnr-1]*yield1; 00399 00400 /* low energies where secondary ionization cannot occur 00401 * will do threshold point, ipLoEnr, later */ 00402 gamk_v = 0.; 00403 thermal.HeatNet = 0.; 00404 00405 /* set up total offset for pointer addition not in loop */ 00406 ipOffset = ipOpac - ipLoEnr; 00407 00408 /* >>>chng 99 apr 16, this had followed the loop below, moved here to 00409 * be like rest of gamma functions */ 00410 /* first do the threshold point 00411 * do not include otscon, which may contain diffuse continuum */ 00412 /* >>chng 01 jul 04, add *rfield.lgOutOnly logic */ 00413 phisig = (rfield.flux[ipLoEnr-1] + rfield.otslin[ipLoEnr-1]+ 00414 rfield.ConInterOut[ipLoEnr-1]*rfield.lgOutOnly) * opac.OpacStack[ipOpac-1]; 00415 gamk_v += phisig; 00416 thermal.HeatNet += phisig*rfield.anu[ipLoEnr-1]; 00417 00418 /* this loop only executed for energies than cannot sec ioniz */ 00419 limit = MIN2(iup,secondaries.ipSecIon-1); 00420 for( i=ipLoEnr; i < limit; i++ ) 00421 { 00422 phisig = rfield.SummedCon[i]*opac.OpacStack[i+ipOffset]; 00423 gamk_v += phisig; 00424 thermal.HeatNet += phisig*rfield.anu[i]; 00425 } 00426 00427 /* correct heating for work function */ 00428 /* >>chng 02 apr 10, from first to sec, due to neg heat in blister.in */ 00429 /*thermal.HeatNet += -gamk_v*eauger;*/ 00430 ASSERT( thermal.HeatNet >= 0. ); 00431 thermal.HeatNet = MAX2(0.,thermal.HeatNet - gamk_v*eauger); 00432 00433 /* this is the total low-energy heating, in ryd, that cannot secondary ionize */ 00434 thermal.HeatLowEnr = thermal.HeatNet; 00435 00436 /* now do part where secondary ionization possible */ 00437 thermal.HeatHiEnr = 0.; 00438 GamHi = 0.; 00439 /* make sure we don't count twice when ipSecIon = ipLoEnr */ 00440 ilo = MAX2(ipLoEnr+1,secondaries.ipSecIon); 00441 for( i=ilo-1; i < iup; i++ ) 00442 { 00443 /* produce of flux and cross section */ 00444 phisig = rfield.SummedCon[i]*opac.OpacStack[i+ipOffset]; 00445 GamHi += phisig; 00446 thermal.HeatHiEnr += phisig*rfield.anu[i]; 00447 } 00448 00449 /* add on the high energy part */ 00450 gamk_v += GamHi; 00451 /* correct for work function */ 00452 thermal.HeatHiEnr -= GamHi*eauger; 00453 00454 /* net heating include high energy heat, with correction for 00455 * secondary ionization */ 00456 thermal.HeatNet += thermal.HeatHiEnr*secondaries.HeatEfficPrimary; 00457 00458 /* final product is heating in erg */ 00459 thermal.HeatNet *= EN1RYD; 00460 thermal.HeatLowEnr *= EN1RYD; 00461 thermal.HeatHiEnr *= EN1RYD; 00462 00463 ASSERT( gamk_v >= 0. ); 00464 ASSERT( thermal.HeatNet>= 0. ); 00465 return( gamk_v ); 00466 } 00467 00468 #if 0 00469 //this is never used 00470 /*GammaPL evaluate photoionization rate for power law photoionization cross section */ 00471 00472 double GammaPL( 00473 /* this is prin quantum number for level */ 00474 long int n, 00475 /* nelem = 0 for H */ 00476 long int nelem) 00477 { 00478 long int i, 00479 ilo, 00480 iup, 00481 limit, 00482 ipLoEnr, 00483 ipHiEnr; 00484 realnum GamHi, 00485 csThresh, 00486 gamPL_v, 00487 phisig; 00488 00489 DEBUG_ENTRY( "GammaPL()" ); 00490 00491 /* evaluate photoioinzation rate and photoelectric heating 00492 * returns photoionization rate as GammaPL, heating is H 00493 * n and nelem and principal quantum number and charge 00494 */ 00495 00496 /* validate the input */ 00497 ASSERT( nelem >= 0); 00498 ASSERT( nelem < LIMELM); 00499 ASSERT( n >= 0); 00500 00501 /* get pointer to photo threshold */ 00502 ipLoEnr = iso.ipIsoLevNIonCon[ipH_LIKE][nelem][n]; 00503 ASSERT( ipLoEnr>0 ); 00504 00505 if( ipLoEnr >= rfield.nflux ) 00506 { 00507 gamPL_v = 0.; 00508 thermal.HeatNet = 0.; 00509 thermal.HeatHiEnr = 0.; 00510 thermal.HeatLowEnr = 0.; 00511 return( gamPL_v ); 00512 } 00513 00514 ipHiEnr = iso.ipIsoLevNIonCon[ipH_LIKE][nelem][ipH2s]; 00515 iup = MIN2(ipHiEnr,rfield.nflux); 00516 00517 csThresh = t_ADfA::Inst().sth(n)*rfield.anu3[ipLoEnr-1]/ POW2(nelem+1.f); 00518 00519 /* low energies where secondary ionization cannot occur 00520 * will do threshold point, ipLoEnr, later */ 00521 gamPL_v = 0.; 00522 thermal.HeatNet = 0.; 00523 00524 /* first do the threshold point 00525 * do not include otscon, which may contain diffuse continuum */ 00526 /* >>chng 01 jul 04, add *rfield.lgOutOnly logic */ 00527 phisig = (rfield.flux[ipLoEnr-1] + rfield.otslin[ipLoEnr-1] + 00528 rfield.ConInterOut[ipLoEnr-1]*rfield.lgOutOnly ) /rfield.anu3[ipLoEnr-1]; 00529 gamPL_v += phisig; 00530 00531 /* this loop only executed for energies than cannot sec ioniz */ 00532 limit = MIN2(iup,secondaries.ipSecIon-1); 00533 for( i=ipLoEnr; i < limit; i++ ) 00534 { 00535 phisig = rfield.SummedCon[i]/rfield.anu3[i]; 00536 gamPL_v += phisig; 00537 thermal.HeatNet += phisig*rfield.anu[i]; 00538 } 00539 00540 /* convert photo rate into proper units with CS at threshold */ 00541 gamPL_v *= csThresh; 00542 00543 /* correct heating for CS at threshold and work function */ 00544 thermal.HeatNet *= csThresh; 00545 thermal.HeatNet += -gamPL_v*rfield.anu[ipLoEnr-1]; 00546 thermal.HeatLowEnr = thermal.HeatNet; 00547 00548 /* now do part where secondary ionization possible */ 00549 thermal.HeatHiEnr = 0.; 00550 GamHi = 0.; 00551 00552 /* make sure we don't count twice when ipSecIon = ipLoEnr */ 00553 ilo = MAX2(ipLoEnr+1,secondaries.ipSecIon); 00554 for( i=ilo-1; i < iup; i++ ) 00555 { 00556 /* produce of flux and cross section */ 00557 phisig = rfield.SummedCon[i]/rfield.anu3[i]; 00558 GamHi += phisig; 00559 thermal.HeatHiEnr += phisig*rfield.anu[i]; 00560 } 00561 00562 /* add on the high energy part */ 00563 gamPL_v += GamHi*csThresh; 00564 00565 /* correct heating for work function and convert to true cross section */ 00566 thermal.HeatHiEnr = (thermal.HeatHiEnr - GamHi*rfield.anu[ipLoEnr-1])*csThresh; 00567 00568 /* net heating includes high energy heating corrected for secondary ionizations */ 00569 thermal.HeatNet += thermal.HeatHiEnr*secondaries.HeatEfficPrimary; 00570 00571 /* final product is heating in erg */ 00572 thermal.HeatNet *= EN1RYD; 00573 thermal.HeatHiEnr *= EN1RYD; 00574 thermal.HeatLowEnr *= EN1RYD; 00575 00576 ASSERT( gamPL_v>= 0.); 00577 return( gamPL_v ); 00578 } 00579 #endif 00580 00581 /*GammaBnPL evaluate hydrogenic photoionization rate for single level with induced recomb */ 00582 00583 double GammaBnPL(long int n, 00584 long int nelem, /* 0 for H, etc */ 00585 double *ainduc, 00586 double *rcool) 00587 { 00588 long int i, 00589 ilo, 00590 iup, 00591 limit, 00592 ipLoEnr, 00593 ipHiEnr; 00594 double GamHi, 00595 bnfunPL_v, 00596 csThresh, 00597 emin, 00598 g, 00599 phisig, 00600 RateInducRec, 00601 RateInducRecCool, 00602 prod, 00603 thresh; 00604 00605 DEBUG_ENTRY( "GammaBnPL()" ); 00606 00607 /* validate the incoming charge */ 00608 ASSERT( nelem >= 0); 00609 ASSERT( nelem < LIMELM ); 00610 00611 /* special version of GAMFUN for arbitrary threshold, and induc rec 00612 * compute ainduc, rate of inducted recombinations, rcool, induc rec cool */ 00613 00614 /* these are the upper and lower energy bounds, which we know from iso-sequence */ 00615 ipLoEnr = iso.ipIsoLevNIonCon[ipH_LIKE][nelem][n]; 00616 ipHiEnr = iso.ipIsoLevNIonCon[ipH_LIKE][nelem][ipH1s]; 00617 00618 ilo = ipLoEnr; 00619 iup = MIN2(ipHiEnr,rfield.nflux); 00620 if( ipLoEnr >= rfield.nflux || ipLoEnr >= ipHiEnr ) 00621 { 00622 bnfunPL_v = 0.; 00623 thermal.HeatNet = 0.; 00624 thermal.HeatHiEnr = 0.; 00625 thermal.HeatLowEnr = 0.; 00626 *ainduc = 0.; 00627 *rcool = 0.; 00628 return( bnfunPL_v ); 00629 } 00630 00631 /* >>chng 96 oct 25, first photo elec energy may have been negative 00632 * possible for first any point to be below threshold due to 00633 * finite resolution of continuum mesh */ 00634 thresh = iso.xIsoLevNIonRyd[ipH_LIKE][nelem][n]; 00635 emin = MAX2(thresh,rfield.anu[ipLoEnr-1]); 00636 /* >>chng 99 jun 08 heating systematically too small, change to correct 00637 * threshold, and protect first cell */ 00638 emin = thresh; 00639 00640 thermal.HeatNet = 0.; 00641 g = 0.; 00642 RateInducRec = 0.; 00643 RateInducRecCool = 0.; 00644 00645 /*csThresh = t_ADfA::Inst().sth(n)*POW3(thresh) / POW2(nelem+1.f);*/ 00646 csThresh = t_ADfA::Inst().sth(n)*rfield.anu3[ipLoEnr-1] / POW2(nelem+1.f); 00647 00648 /* do first point without otscon, which may have diffuse cont */ 00649 /* >>chng 01 jul 04, add *rfield.lgOutOnly logic */ 00650 phisig = (rfield.flux[ipLoEnr-1] + rfield.otslin[ipLoEnr-1]+ 00651 rfield.ConInterOut[ipLoEnr-1]*rfield.lgOutOnly)/ 00652 rfield.anu3[ipLoEnr-1]; 00653 00654 g += phisig; 00655 thermal.HeatNet += phisig*rfield.anu[ipLoEnr-1]; 00656 00657 /* integral part of induced rec rate */ 00658 prod = phisig*rfield.ContBoltz[ipLoEnr-1]; 00659 RateInducRec += prod; 00660 /* induced rec cooling */ 00661 RateInducRecCool += prod*(rfield.anu[ipLoEnr-1] - emin); 00662 00663 limit = MIN2(iup,secondaries.ipSecIon-1); 00664 for( i=ipLoEnr; i < limit; i++ ) 00665 { 00666 phisig = rfield.SummedCon[i]/rfield.anu3[i]; 00667 00668 g += phisig; 00669 thermal.HeatNet += phisig*rfield.anu[i]; 00670 00671 /* integral part of induced rec rate */ 00672 prod = phisig*rfield.ContBoltz[i]; 00673 RateInducRec += prod; 00674 /* incuded rec cooling */ 00675 RateInducRecCool += prod*(rfield.anu[i] - emin); 00676 } 00677 00678 /* heating in Rydbergs, no secondary ionization */ 00679 thermal.HeatNet = MAX2(0.,(thermal.HeatNet - g*thresh)*csThresh); 00680 thermal.HeatLowEnr = thermal.HeatNet; 00681 00682 /* now do part where secondary ionization is possible */ 00683 thermal.HeatHiEnr = 0.; 00684 GamHi = 0.; 00685 /* make sure we don't count twice when ipSecIon = ipLoEnr */ 00686 ilo = MAX2(ipLoEnr+1,secondaries.ipSecIon); 00687 for( i=ilo-1; i < iup; i++ ) 00688 { 00689 phisig = rfield.SummedCon[i]/rfield.anu3[i]; 00690 GamHi += phisig; 00691 thermal.HeatHiEnr += phisig*rfield.anu[i]; 00692 00693 /* integral part of induced rec rate */ 00694 prod = phisig*rfield.ContBoltz[i]; 00695 RateInducRec += prod; 00696 /* incuded rec cooling */ 00697 RateInducRecCool += prod*(rfield.anu[i] - emin); 00698 } 00699 RateInducRec *= csThresh; 00700 RateInducRecCool *= csThresh; 00701 00702 /* this is total photo rate, low and high energy parts */ 00703 bnfunPL_v = (g + GamHi)*csThresh; 00704 00705 /* heating in Rydbergs, uncorrected for secondary ionization */ 00706 thermal.HeatHiEnr = (thermal.HeatHiEnr - GamHi*thresh)*csThresh; 00707 00708 /* net heating includes high energy heat corrected for secondary ionization */ 00709 thermal.HeatNet += thermal.HeatHiEnr*secondaries.HeatEfficPrimary; 00710 00711 /* final product is heating in erg */ 00712 thermal.HeatNet *= EN1RYD; 00713 thermal.HeatHiEnr *= EN1RYD; 00714 thermal.HeatLowEnr *= EN1RYD; 00715 00716 /* this is an option to turn off induced processes */ 00717 if( rfield.lgInducProcess ) 00718 { 00719 *rcool = RateInducRecCool*EN1RYD; 00720 *ainduc = RateInducRec; 00721 } 00722 else 00723 { 00724 /* the "no induced" command was given */ 00725 *rcool = 0.; 00726 *ainduc = 0.; 00727 } 00728 00729 ASSERT( bnfunPL_v>= 0.); 00730 return( bnfunPL_v ); 00731 } 00732 00733 /* GammaPrtRate will print resulting rates for ion and element */ 00734 void GammaPrtRate( 00735 /* io unit we will write to */ 00736 FILE * ioFILE, 00737 /* stage of ionization on C scale, 0 for atom */ 00738 long int ion , 00739 /* 0 for H, etc */ 00740 long int nelem, 00741 /* true - then print photo sources for each shell */ 00742 bool lgPRT ) 00743 { 00744 long int nshell , ns; 00745 00746 DEBUG_ENTRY( "GammaPrtRate()" ); 00747 00748 /* number of shells for this species */ 00749 nshell = Heavy.nsShells[nelem][ion]; 00750 00751 /* now print subshell photo rate */ 00752 fprintf(ioFILE , "GammaPrtRate: %li %li",ion , nelem ); 00753 for( ns=nshell-1; ns>=0; --ns ) 00754 { 00755 fprintf(ioFILE , " %.2e" , ionbal.PhotoRate_Shell[nelem][ion][ns][0]); 00756 00757 /* option to print individual contributors to each shell */ 00758 if( lgPRT ) 00759 { 00760 fprintf(ioFILE , "\n"); 00761 GammaPrt( 00762 /* first three par are identical to GammaK */ 00763 opac.ipElement[nelem][ion][ns][0], 00764 opac.ipElement[nelem][ion][ns][1], 00765 opac.ipElement[nelem][ion][ns][2], 00766 /* io unit we will write to */ 00767 ioFILE, 00768 /* total photo rate from previous call to GammaK */ 00769 ionbal.PhotoRate_Shell[nelem][ion][ns][0], 00770 /* we will print contributors that are more than this rate */ 00771 ionbal.PhotoRate_Shell[nelem][ion][ns][0]*0.05); 00772 } 00773 } 00774 fprintf(ioFILE , "\n"); 00775 return; 00776 }