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 /*iso_cool compute net cooling due to species in iso-sequences */ 00004 #include "cddefines.h" 00005 #include "physconst.h" 00006 #include "taulines.h" 00007 #include "hydrogenic.h" 00008 #include "elementnames.h" 00009 #include "phycon.h" 00010 #include "dense.h" 00011 #include "thermal.h" 00012 #include "cooling.h" 00013 #include "iso.h" 00014 00015 /* HP cc cannot compile this routine with any optimization */ 00016 #if defined(__HP_aCC) 00017 # pragma OPT_LEVEL 1 00018 #endif 00019 00020 // set to true to enable debug print of contributors to collisional ionization cooling 00021 const bool lgPrintIonizCooling = false; 00022 00023 void iso_cool( 00024 /* iso sequence, 0 for hydrogenic */ 00025 long int ipISO , 00026 /* nelem is charge -1, so 0 for H itself */ 00027 long int nelem) 00028 { 00029 long int ipHi, 00030 ipbig, 00031 ipLo, 00032 n; 00033 double RecCoolExtra, 00034 biggest = 0., 00035 dCdT_all, 00036 edenHCorr_IonAbund, 00037 edenIonAbund, 00038 CollisIonizatCoolingTotal, 00039 dCollisIonizatCoolingTotalDT, 00040 HeatExcited, 00041 heat_max, 00042 CollisIonizatCooling, 00043 CollisIonizatCoolingDT, 00044 hlone, 00045 thin, 00046 ThinCoolingCaseB, 00047 ThinCoolingSum, 00048 collider; 00049 00050 valarray<double> CollisIonizatCoolingArr, 00051 CollisIonizatCoolingDTArr, 00052 SavePhotoHeat, 00053 SaveInducCool, 00054 SaveRadRecCool; 00055 00056 long int nlo_heat_max , nhi_heat_max; 00057 00058 /* place to put string labels for iso lines */ 00059 char chLabel[NCOLNT_LAB_LEN+1]; 00060 00061 DEBUG_ENTRY( "iso_cool()" ); 00062 00063 /* validate the incoming data */ 00064 ASSERT( nelem >= ipISO ); 00065 ASSERT( ipISO <NISO ); 00066 ASSERT( nelem < LIMELM ); 00067 /* local number of levels may be less than malloced number if continuum 00068 * has been lowered due to high density */ 00069 ASSERT( iso.numLevels_local[ipISO][nelem] <= iso.numLevels_max[ipISO][nelem] ); 00070 00071 /* for zero abundance or if element has been turned off we need to 00072 * set some produced variables to zero */ 00073 if( dense.xIonDense[nelem][nelem+1-ipISO] <= 0. || 00074 /* >>chng 04 dec 07, add test for recombined species having zero abundance, 00075 * this occurs when set ion is in place so very artificial */ 00076 dense.xIonDense[nelem][nelem-ipISO]<=0. || 00077 !dense.lgElmtOn[nelem] ) 00078 { 00079 /* all global variables must be zeroed here or set below */ 00080 iso.coll_ion[ipISO][nelem] = 0.; 00081 iso.cLya_cool[ipISO][nelem] = 0.; 00082 iso.cLyrest_cool[ipISO][nelem] = 0.; 00083 iso.cBal_cool[ipISO][nelem] = 0.; 00084 iso.cRest_cool[ipISO][nelem] = 0.; 00085 iso.xLineTotCool[ipISO][nelem] = 0.; 00086 iso.RadRecCool[ipISO][nelem] = 0.; 00087 iso.FreeBnd_net_Cool_Rate[ipISO][nelem] = 0.; 00088 iso.dLTot[ipISO][nelem] = 0.; 00089 iso.RecomInducCool_Rate[ipISO][nelem] = 0.; 00090 return; 00091 } 00092 /* >>chng 05 aug 17, must use real electron density for collisional ionization of H 00093 * since in Leiden v4 pdr with its artificial high temperature coll ion can be important 00094 * H on H is homonuclear and scaling laws for other elements does not apply 00095 * next two were multiplied by dense.EdenHCorr and changed to dense.eden for H, 00096 * EdenHCorr for rest */ 00097 if( nelem==ipHYDROGEN ) 00098 { 00099 /* special version for H0 onto H0 */ 00100 collider = dense.EdenHontoHCorr; 00101 } 00102 else 00103 { 00104 collider = dense.EdenHCorr; 00105 } 00106 00107 /* create some space, these go to numLevels_local instead of _max 00108 * since continuum may have been lowered by density */ 00109 if( lgPrintIonizCooling ) 00110 { 00111 CollisIonizatCoolingArr.resize( iso.numLevels_local[ipISO][nelem] ); 00112 CollisIonizatCoolingDTArr.resize( iso.numLevels_local[ipISO][nelem] ); 00113 } 00114 SavePhotoHeat.resize( iso.numLevels_local[ipISO][nelem] ); 00115 SaveInducCool.resize( iso.numLevels_local[ipISO][nelem] ); 00116 SaveRadRecCool.resize( iso.numLevels_local[ipISO][nelem] ); 00117 00118 /*********************************************************************** 00119 * * 00120 * collisional ionization cooling, less three-body recombination heat * 00121 * * 00122 ***********************************************************************/ 00123 00124 /* will be net collisional ionization cooling, units erg/cm^3/s */ 00125 CollisIonizatCoolingTotal = 0.; 00126 dCollisIonizatCoolingTotalDT = 0.; 00127 00128 /* collisional ionization cooling minus three body heating 00129 * depending on how topoff is done, highest level can have large population 00130 * and its coupling to continuum can be large, at various times code 00131 * had to ignore effects of very highest level, but starting in mid 00132 * 20006 all levels have been included 00133 * 2008 apr 18, do not include highest - when only 2 collapsed levels 00134 * are used several density 13 BLR sims have serious convergence 00135 * problems */ 00136 for( n=0; n < iso.numLevels_local[ipISO][nelem]-1; ++n ) 00137 { 00138 /* total collisional ionization cooling less three body heating */ 00139 CollisIonizatCooling = 00140 iso.xIsoLevNIonRyd[ipISO][nelem][n]*iso.ColIoniz[ipISO][nelem][n]*collider* 00141 (StatesElem[ipISO][nelem][n].Pop -iso.PopLTE[ipISO][nelem][n]*dense.eden); 00142 CollisIonizatCoolingTotal += CollisIonizatCooling; 00143 00144 /* the derivative of the cooling 00145 * need extra factor of temp of 1 Ryd since div by square of temp in Ryd */ 00146 CollisIonizatCoolingDT = CollisIonizatCooling* 00147 (iso.xIsoLevNIonRyd[ipISO][nelem][n]/TE1RYD/POW2(phycon.te_ryd)); 00148 00149 dCollisIonizatCoolingTotalDT += CollisIonizatCoolingDT; 00150 // save values for debug printout 00151 if( lgPrintIonizCooling ) 00152 { 00153 CollisIonizatCoolingArr[n] = CollisIonizatCooling; 00154 CollisIonizatCoolingDTArr[n] = CollisIonizatCoolingDT; 00155 } 00156 } 00157 00158 /* save net collisional ionization cooling less H-3 body heating 00159 * for inclusion in printout 00160 * convert to physical units, need to convert Ryd to ergs, 00161 * and bring to density per vol not per ion */ 00162 iso.coll_ion[ipISO][nelem] = CollisIonizatCoolingTotal * EN1RYD* 00163 dense.xIonDense[nelem][nelem+1-ipISO]; 00164 00165 /* derivative wrt temp 00166 * convert to physical units, need to convert Ryd to ergs, 00167 * and bring to density per vol not per ion */ 00168 double dcool = dCollisIonizatCoolingTotalDT * EN1RYD * 00169 dense.xIonDense[nelem][nelem+1-ipISO]; 00170 fixit(); 00173 dcool = iso.coll_ion[ipISO][nelem] / phycon.te; 00174 00175 /* create a meaningful label */ 00176 sprintf(chLabel , "IScion%2s%2s" , elementnames.chElementSym[ipISO] , 00177 elementnames.chElementSym[nelem] ); 00178 /* dump the coolant onto the cooling stack */ 00179 CoolAdd(chLabel,(realnum)nelem,MAX2(0.,iso.coll_ion[ipISO][nelem])); 00180 00181 thermal.dCooldT += dcool; 00182 00183 /* heating[0][3] is all three body heating, opposite of collisional 00184 * ionization cooling, 00185 * would be unusual for this to be non-zero since would require excited 00186 * state departure coefficients to be greater than unity */ 00187 thermal.heating[0][3] += MAX2(0.,-iso.coll_ion[ipISO][nelem]); 00188 00189 /* debug block printing individual contributors to collisional ionization cooling */ 00190 if( lgPrintIonizCooling && nelem==1 && ipISO==1 ) 00191 { 00192 fprintf(ioQQQ,"DEBUG coll ioniz cool contributors:"); 00193 for( n=0; n < iso.numLevels_local[ipISO][nelem]; n++ ) 00194 { 00195 if( CollisIonizatCoolingArr[n] / SDIV( CollisIonizatCoolingTotal ) > 0.1 ) 00196 fprintf(ioQQQ," %2li %.1e", 00197 n, 00198 CollisIonizatCoolingArr[n]/ SDIV( CollisIonizatCoolingTotal ) ); 00199 } 00200 fprintf(ioQQQ,"\n"); 00201 fprintf(ioQQQ,"DEBUG coll ioniz derivcontributors:"); 00202 for( n=0; n < iso.numLevels_local[ipISO][nelem]; n++ ) 00203 { 00204 if( CollisIonizatCoolingDTArr[n] / SDIV( dCollisIonizatCoolingTotalDT ) > 0.1 ) 00205 fprintf(ioQQQ," %2li %.1e", 00206 n, 00207 CollisIonizatCoolingDTArr[n]/ SDIV( dCollisIonizatCoolingTotalDT ) ); 00208 } 00209 fprintf(ioQQQ,"\n"); 00210 } 00211 00212 /*********************************************************************** 00213 * * 00214 * hydrogen recombination free-bound free bound cooling * 00215 * * 00216 ***********************************************************************/ 00217 00218 /* this is the product of the ion abundance times the electron density, 00219 * will multiply level pops which are stored relative to ion 00220 * EdenHCorr is used in level pops, so should be used here too */ 00221 edenIonAbund = dense.eden*dense.xIonDense[nelem][nelem+1-ipISO]; 00222 00223 /* this is the product of the ion abundance times the electron density with 00224 * a small correction for the presence of neutrals, which should be used in 00225 * reactions that involve collisions between states, but NOT radiative recombination 00226 * but should be used in collisional ionization / three body recombination */ 00227 edenHCorr_IonAbund = collider*dense.xIonDense[nelem][nelem+1-ipISO]; 00228 00229 /* now do case b sum to compare with exact value below */ 00230 iso.RadRecCool[ipISO][nelem] = 0.; 00231 ThinCoolingSum = 0.; 00232 00233 if( ipISO == ipH_LIKE ) 00234 { 00235 /* do ground with special approximate fits to Ferland et al. '92 */ 00236 thin = HydroRecCool( 00237 /* n is the prin quantum number on the physical scale */ 00238 1 , 00239 /* nelem is the charge on the C scale, 0 is hydrogen */ 00240 nelem); 00241 } 00242 else 00243 { 00244 /* this is the cooling before correction for optical depths */ 00245 thin = iso.RadRecomb[ipISO][nelem][0][ipRecRad]* 00246 /* arg is the scaled temperature, T * n^2 / Z^2, 00247 * n is principal quantum number, Z is charge, 1 for H */ 00248 HCoolRatio( 00249 phycon.te * POW2( (double)StatesElem[ipISO][nelem][0].n / (double)(nelem+1-ipISO) ))* 00250 /* convert results to energy per unit vol */ 00251 phycon.te * BOLTZMANN; 00252 } 00253 /* the cooling, corrected for optical depth */ 00254 SaveRadRecCool[0] = iso.RadRecomb[ipISO][nelem][0][ipRecNetEsc] * thin; 00255 /* this is now total free-bound cooling */ 00256 iso.RadRecCool[ipISO][nelem] += SaveRadRecCool[0] * edenIonAbund; 00257 00258 /* radiative recombination cooling for all excited states */ 00259 for( n=1; n < iso.numLevels_local[ipISO][nelem]; n++ ) 00260 { 00261 /* this is the cooling before correction for optical depths */ 00262 thin = iso.RadRecomb[ipISO][nelem][n][ipRecRad]* 00263 /* arg is the scaled temperature, T * n^2 / Z^2, 00264 * n is principal quantum number, Z is charge, 1 for H */ 00265 HCoolRatio( 00266 phycon.te * POW2( (double)StatesElem[ipISO][nelem][n].n / (double)(nelem+1-ipISO) ))* 00267 /* convert results to energy per unit vol */ 00268 phycon.te * BOLTZMANN; 00269 00270 /* the cooling, corrected for optical depth */ 00271 SaveRadRecCool[n] = iso.RadRecomb[ipISO][nelem][n][ipRecNetEsc] * thin * edenIonAbund; 00272 /* this is now total free-bound cooling */ 00273 iso.RadRecCool[ipISO][nelem] += SaveRadRecCool[n]; 00274 00275 /* keep track of case b sum for topoff below */ 00276 ThinCoolingSum += thin; 00277 } 00278 { 00279 /* debug block for state specific recombination cooling */ 00280 enum {DEBUG_LOC=false}; 00281 if( DEBUG_LOC ) 00282 { 00283 if( nelem==ipISO ) 00284 { 00285 /* >>chng 06 aug 17, should go to numLevels_local instead of _max. */ 00286 for( n=0; n < (iso.numLevels_local[ipISO][nelem] - 1); n++ ) 00287 { 00288 fprintf(ioQQQ,"\t%.2f",SaveRadRecCool[n]/ThinCoolingSum); 00289 } 00290 fprintf(ioQQQ,"\n"); 00291 } 00292 } 00293 } 00294 00295 /* Case b sum of optically thin radiative recombination cooling. 00296 * add any remainder to the sum from above - high precision is needed 00297 * to get STE result to converge close to equilibrium - only done for 00298 * H-like ions where exact result is known */ 00299 if( ipISO == ipH_LIKE ) 00300 { 00301 /* these expressions are only valid for hydrogenic sequence */ 00302 if( nelem == 0 ) 00303 { 00304 /*expansion for hydrogen itself */ 00305 ThinCoolingCaseB = (-25.859117 + 00306 0.16229407*phycon.telogn[0] + 00307 0.34912863*phycon.telogn[1] - 00308 0.10615964*phycon.telogn[2])/(1. + 00309 0.050866793*phycon.telogn[0] - 00310 0.014118924*phycon.telogn[1] + 00311 0.0044980897*phycon.telogn[2] + 00312 6.0969594e-5*phycon.telogn[3]); 00313 } 00314 else 00315 { 00316 /* same expansion but for hydrogen ions */ 00317 ThinCoolingCaseB = (-25.859117 + 00318 0.16229407*(phycon.telogn[0]-phycon.sqlogz[nelem-ipISO]) + 00319 0.34912863*POW2(phycon.telogn[0]-phycon.sqlogz[nelem-ipISO]) - 00320 0.10615964*powi( (phycon.telogn[0]-phycon.sqlogz[nelem-ipISO]),3) )/(1. + 00321 0.050866793*(phycon.telogn[0]-phycon.sqlogz[nelem-ipISO]) - 00322 0.014118924*POW2(phycon.telogn[0]-phycon.sqlogz[nelem-ipISO]) + 00323 0.0044980897*powi( (phycon.telogn[0]-phycon.sqlogz[nelem-ipISO]),3) + 00324 6.0969594e-5*powi( (phycon.telogn[0]-phycon.sqlogz[nelem-ipISO]),4) ); 00325 } 00326 00327 /* now convert to linear cooling coefficient */ 00328 ThinCoolingCaseB = POW3(1.+nelem-ipISO)*pow(10.,ThinCoolingCaseB)/(phycon.te/POW2(1.+nelem-ipISO) ); 00329 00330 /* this is the error, expect positive since do not include infinite number of levels */ 00331 RecCoolExtra = ThinCoolingCaseB - ThinCoolingSum; 00332 } 00333 else 00334 { 00335 ThinCoolingCaseB = 0.; 00336 RecCoolExtra = 0.; 00337 } 00338 00339 /* don't let the extra be negative - should be positive if H-like, negative for 00340 * he-like only due to real difference in recombination coefficients */ 00341 RecCoolExtra = MAX2(0., RecCoolExtra ); 00342 00343 /* add error onto total - this is significant for approach to STE */ 00344 iso.RadRecCool[ipISO][nelem] += RecCoolExtra* edenIonAbund *iso.RadRecomb[ipISO][nelem][iso.numLevels_local[ipISO][nelem]-1][ipRecNetEsc]; 00345 00346 /*********************************************************************** 00347 * 00348 * heating by photoionization of 00349 * excited states of all species 00350 * 00351 ***********************************************************************/ 00352 00353 /* photoionization of excited levels */ 00354 HeatExcited = 0.; 00355 ipbig = -1000; 00356 for( n=1; n < (iso.numLevels_local[ipISO][nelem] - 1); ++n ) 00357 { 00358 SavePhotoHeat[n] = dense.xIonDense[nelem][nelem+1-ipISO]*StatesElem[ipISO][nelem][n].Pop* 00359 iso.PhotoHeat[ipISO][nelem][n]; 00360 HeatExcited += SavePhotoHeat[n]; 00361 if( SavePhotoHeat[n] > biggest ) 00362 { 00363 biggest = SavePhotoHeat[n]; 00364 ipbig = (int)n; 00365 } 00366 } 00367 { 00368 /* debug block for heating due to photo of each n */ 00369 enum {DEBUG_LOC=false}; 00370 if( DEBUG_LOC && ipISO==0 && nelem==0 && nzone > 700) 00371 { 00372 /* this was not done above */ 00373 SavePhotoHeat[ipH1s] = dense.xIonDense[nelem][nelem+1-ipISO]*StatesElem[ipISO][nelem][ipH1s].Pop* 00374 iso.PhotoHeat[ipISO][nelem][ipH1s]; 00375 fprintf(ioQQQ,"ipISO = %li nelem=%li ipbig=%li biggest=%g HeatExcited=%.2e ctot=%.2e\n", 00376 ipISO , nelem, 00377 ipbig , 00378 biggest, 00379 HeatExcited , 00380 thermal.ctot); 00381 /* >>chng 06 aug 17, should go to numLevels_local instead of _max. */ 00382 for(n=ipH1s; n< (iso.numLevels_local[ipISO][nelem] - 1); ++n ) 00383 { 00384 fprintf(ioQQQ,"DEBUG phot heat%2li\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\n", 00385 n, 00386 SavePhotoHeat[n]/HeatExcited, 00387 dense.xIonDense[nelem][nelem+1-ipISO], 00388 StatesElem[ipISO][nelem][n].Pop, 00389 iso.PhotoHeat[ipISO][nelem][n], 00390 iso.gamnc[ipISO][nelem][n]); 00391 } 00392 } 00393 } 00394 00395 /* FreeBnd_net_Cool_Rate is net cooling due to recombination 00396 * RadRecCool is total radiative recombination cooling sum to all levels, 00397 * with n>=2 photoionization heating subtracted */ 00398 iso.FreeBnd_net_Cool_Rate[ipISO][nelem] = iso.RadRecCool[ipISO][nelem] - HeatExcited; 00399 /*fprintf(ioQQQ,"DEBUG Hn2\t%.3e\t%.3e\n", 00400 -iso.RadRecCool[ipISO][nelem]/SDIV(thermal.htot), 00401 HeatExcited/SDIV(thermal.htot));*/ 00402 00403 /* heating[0][1] is all excited state photoionization heating from ALL 00404 * species, this is set to zero in CoolEvaluate before loop where this 00405 * is called. */ 00406 thermal.heating[0][1] += MAX2(0.,-iso.FreeBnd_net_Cool_Rate[ipISO][nelem]); 00407 00408 /* net free-bound cooling minus free-free heating */ 00409 /* create a meaningful label */ 00410 sprintf(chLabel , "ISrcol%2s%2s" , elementnames.chElementSym[ipISO] , 00411 elementnames.chElementSym[nelem]); 00412 CoolAdd(chLabel, (realnum)nelem, MAX2(0.,iso.FreeBnd_net_Cool_Rate[ipISO][nelem])); 00413 00414 /* if rec coef goes as T^-.8, but times T, so propto t^.2 */ 00415 thermal.dCooldT += 0.2*iso.FreeBnd_net_Cool_Rate[ipISO][nelem]*phycon.teinv; 00416 00417 /*********************************************************************** 00418 * * 00419 * induced recombination cooling * 00420 * * 00421 ***********************************************************************/ 00422 00423 iso.RecomInducCool_Rate[ipISO][nelem] = 0.; 00424 /* >>chng 06 aug 17, should go to numLevels_local instead of _max. */ 00425 for( n=0; n < (iso.numLevels_local[ipISO][nelem] - 1); ++n ) 00426 { 00427 /* >>chng 02 jan 22, removed cinduc, replace with RecomInducCool */ 00428 SaveInducCool[n] = iso.RecomInducCool_Coef[ipISO][nelem][n]*iso.PopLTE[ipISO][nelem][n]*edenIonAbund; 00429 iso.RecomInducCool_Rate[ipISO][nelem] += SaveInducCool[n]; 00430 thermal.dCooldT += SaveInducCool[n]* 00431 (iso.xIsoLevNIonRyd[ipISO][nelem][n]/phycon.te_ryd - 1.5)*phycon.teinv; 00432 } 00433 00434 { 00435 /* print rec cool, induced rec cool, photo heat */ 00436 enum {DEBUG_LOC=false}; 00437 if( DEBUG_LOC && ipISO==0 && nelem==5 ) 00438 { 00439 fprintf(ioQQQ," ipISO=%li nelem=%li ctot = %.2e\n", 00440 ipISO, 00441 nelem, 00442 thermal.ctot); 00443 fprintf(ioQQQ,"sum\t%.2e\t%.2e\t%.2e\n", 00444 HeatExcited, 00445 iso.RadRecCool[ipISO][nelem], 00446 iso.RecomInducCool_Rate[ipISO][nelem]); 00447 fprintf(ioQQQ,"sum\tp ht\tr cl\ti cl\n"); 00448 00449 /* >>chng 06 aug 17, should go to numLevels_local instead of _max. */ 00450 for(n=0; n< (iso.numLevels_local[ipISO][nelem] - 1); ++n ) 00451 { 00452 fprintf(ioQQQ,"%li\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e \n", 00453 n, 00454 SavePhotoHeat[n], 00455 SaveRadRecCool[n], 00456 SaveInducCool[n] , 00457 iso.RecomInducCool_Coef[ipISO][nelem][n], 00458 iso.PopLTE[ipISO][nelem][n], 00459 iso.RecomInducRate[ipISO][nelem][n]); 00460 } 00461 fprintf(ioQQQ," \n"); 00462 } 00463 } 00464 /* create a meaningful label - induced rec cooling */ 00465 sprintf(chLabel , "ISicol%2s%2s" , elementnames.chElementSym[ipISO] , 00466 elementnames.chElementSym[nelem]); 00467 /* induced rec cooling */ 00468 CoolAdd(chLabel,(realnum)nelem,iso.RecomInducCool_Rate[ipISO][nelem]); 00469 00470 /* find total collisional energy exchange due to bound-bound */ 00471 iso.xLineTotCool[ipISO][nelem] = 0.; 00472 dCdT_all = 0.; 00473 heat_max = 0.; 00474 nlo_heat_max = -1; 00475 nhi_heat_max = -1; 00476 00477 /* loop does not include highest levels - their population may be 00478 * affected by topoff */ 00479 for( ipLo=0; ipLo < iso.numLevels_local[ipISO][nelem]-2; ipLo++ ) 00480 { 00481 for( ipHi=ipLo + 1; ipHi < iso.numLevels_local[ipISO][nelem]-1; ipHi++ ) 00482 { 00483 /* collisional energy exchange between ipHi and ipLo - net cool */ 00484 hlone = 00485 Transitions[ipISO][nelem][ipHi][ipLo].Coll.ColUL* 00486 (StatesElem[ipISO][nelem][ipLo].Pop* 00487 iso.Boltzmann[ipISO][nelem][ipHi][ipLo]* 00488 StatesElem[ipISO][nelem][ipHi].g/StatesElem[ipISO][nelem][ipLo].g - 00489 StatesElem[ipISO][nelem][ipHi].Pop)*edenHCorr_IonAbund* 00490 Transitions[ipISO][nelem][ipHi][ipLo].EnergyErg; 00491 00492 iso.xLineTotCool[ipISO][nelem] += hlone; 00493 00494 /* next get derivative */ 00495 if( hlone > 0. ) 00496 { 00497 /* usual case, this line was a net coolant - for derivative 00498 * taking the exponential from ground gave better 00499 * representation of effects */ 00500 dCdT_all += hlone* 00501 (Transitions[ipISO][nelem][ipHi][ipH1s].EnergyK*thermal.tsq1 - thermal.halfte); 00502 } 00503 else 00504 { 00505 /* this line heats the gas, remember which one it was, 00506 * the following three vars never are used, but could be for 00507 * debugging */ 00508 if( hlone < heat_max ) 00509 { 00510 heat_max = hlone; 00511 nlo_heat_max = ipLo; 00512 nhi_heat_max = ipHi; 00513 } 00514 dCdT_all -= hlone*thermal.halfte; 00515 } 00516 } 00517 } 00518 { 00519 /* debug block announcing which line was most important */ 00520 enum {DEBUG_LOC=false}; 00521 if( DEBUG_LOC ) 00522 { 00523 if( nelem==ipISO ) 00524 fprintf(ioQQQ,"%li %li %.2e\n", nlo_heat_max, nhi_heat_max, heat_max ); 00525 } 00526 } 00527 00528 /* Lyman line collisional heating/cooling */ 00529 /* Lya itself */ 00530 iso.cLya_cool[ipISO][nelem] = 0.; 00531 /* lines higher than Lya */ 00532 iso.cLyrest_cool[ipISO][nelem] = 0.; 00533 00534 for( ipHi=1; ipHi < iso.numLevels_local[ipISO][nelem]; ipHi++ ) 00535 { 00536 hlone = Transitions[ipISO][nelem][ipHi][ipH1s].Coll.ColUL* 00537 (StatesElem[ipISO][nelem][0].Pop*iso.Boltzmann[ipISO][nelem][ipHi][0]* 00538 StatesElem[ipISO][nelem][ipHi].g/StatesElem[ipISO][nelem][0].g - 00539 StatesElem[ipISO][nelem][ipHi].Pop)* edenHCorr_IonAbund* 00540 Transitions[ipISO][nelem][ipHi][0].EnergyErg; 00541 00542 if( ipHi == iso.nLyaLevel[ipISO] ) 00543 { 00544 iso.cLya_cool[ipISO][nelem] = hlone; 00545 00546 } 00547 else 00548 { 00549 /* sum energy in higher lyman lines */ 00550 iso.cLyrest_cool[ipISO][nelem] += hlone; 00551 } 00552 } 00553 00554 /* Balmer line collisional heating/cooling and derivative 00555 * only used for print out, incorrect if not H */ 00556 iso.cBal_cool[ipISO][nelem] = 0.; 00557 for( ipHi=3; ipHi < iso.numLevels_local[ipISO][nelem]; ipHi++ ) 00558 { 00559 /* single line cooling */ 00560 ipLo = ipH2s; 00561 hlone = 00562 Transitions[ipISO][nelem][ipHi][ipLo].Coll.ColUL*( 00563 StatesElem[ipISO][nelem][ipLo].Pop* 00564 iso.Boltzmann[ipISO][nelem][ipHi][ipLo]* 00565 StatesElem[ipISO][nelem][ipHi].g/StatesElem[ipISO][nelem][ipLo].g - 00566 StatesElem[ipISO][nelem][ipHi].Pop)*edenHCorr_IonAbund* 00567 Transitions[ipISO][nelem][ipHi][ipLo].EnergyErg; 00568 00569 ipLo = ipH2p; 00570 hlone += 00571 Transitions[ipISO][nelem][ipHi][ipLo].Coll.ColUL*( 00572 StatesElem[ipISO][nelem][ipLo].Pop* 00573 iso.Boltzmann[ipISO][nelem][ipHi][ipLo]* 00574 StatesElem[ipISO][nelem][ipHi].g/StatesElem[ipISO][nelem][ipLo].g - 00575 StatesElem[ipISO][nelem][ipHi].Pop)*edenHCorr_IonAbund* 00576 Transitions[ipISO][nelem][ipHi][ipLo].EnergyErg; 00577 00578 /* this is only used to add to line array */ 00579 iso.cBal_cool[ipISO][nelem] += hlone; 00580 } 00581 00582 /* all hydrogen lines from Paschen on up, but not Balmer 00583 * only used for printout, incorrect if not H */ 00584 iso.cRest_cool[ipISO][nelem] = 0.; 00585 for( ipLo=3; ipLo < iso.numLevels_local[ipISO][nelem]-1; ipLo++ ) 00586 { 00587 for( ipHi=ipLo + 1; ipHi < iso.numLevels_local[ipISO][nelem]; ipHi++ ) 00588 { 00589 hlone = 00590 Transitions[ipISO][nelem][ipHi][ipLo].Coll.ColUL*( 00591 StatesElem[ipISO][nelem][ipLo].Pop* 00592 iso.Boltzmann[ipISO][nelem][ipHi][ipLo]* 00593 StatesElem[ipISO][nelem][ipHi].g/StatesElem[ipISO][nelem][ipLo].g - 00594 StatesElem[ipISO][nelem][ipHi].Pop)*edenHCorr_IonAbund* 00595 Transitions[ipISO][nelem][ipHi][ipLo].EnergyErg; 00596 00597 iso.cRest_cool[ipISO][nelem] += hlone; 00598 } 00599 } 00600 00601 /* add total line heating or cooling into stacks, derivatives */ 00602 /* line energy exchange can be either heating or coolant 00603 * must add this to total heating or cooling, as appropriate */ 00604 /* create a meaningful label */ 00605 sprintf(chLabel , "ISclin%2s%2s" , elementnames.chElementSym[ipISO] , 00606 elementnames.chElementSym[nelem]); 00607 if( iso.xLineTotCool[ipISO][nelem] > 0. ) 00608 { 00609 /* species is a net coolant label gives iso sequence, "wavelength" gives element */ 00610 CoolAdd(chLabel,(realnum)nelem,iso.xLineTotCool[ipISO][nelem]); 00611 thermal.dCooldT += dCdT_all; 00612 iso.dLTot[ipISO][nelem] = 0.; 00613 } 00614 else 00615 { 00616 /* species is a net heat source, thermal.heating[0][23]was set to 0 in CoolEvaluate*/ 00617 thermal.heating[0][23] -= iso.xLineTotCool[ipISO][nelem]; 00618 CoolAdd(chLabel,(realnum)nelem,0.); 00619 iso.dLTot[ipISO][nelem] = -dCdT_all; 00620 } 00621 00622 { 00623 /* debug print for understanding contributors to HLineTotCool */ 00624 enum {DEBUG_LOC=false}; 00625 if( DEBUG_LOC ) 00626 { 00627 if( nelem == 0 ) 00628 { 00629 fprintf(ioQQQ,"%.2e la %.2f restly %.2f barest %.2f hrest %.2f\n", 00630 iso.xLineTotCool[ipISO][nelem] , 00631 iso.cLya_cool[ipISO][nelem]/iso.xLineTotCool[ipISO][nelem] , 00632 iso.cLyrest_cool[ipISO][nelem]/iso.xLineTotCool[ipISO][nelem] , 00633 iso.cBal_cool[ipISO][nelem]/iso.xLineTotCool[ipISO][nelem] , 00634 iso.cRest_cool[ipISO][nelem]/iso.xLineTotCool[ipISO][nelem] ); 00635 } 00636 } 00637 } 00638 { 00639 /* this is an option to print out each rate affecting each level in strict ste, 00640 * this is a check that the rates are indeed in detailed balance */ 00641 enum {DEBUG_LOC=false}; 00642 enum {LTEPOP=true}; 00643 /* special debug print for gas near STE */ 00644 if( DEBUG_LOC && (nelem==1 || nelem==0) ) 00645 { 00646 /* LTEPOP is option to assume STE for rates */ 00647 if( LTEPOP ) 00648 { 00649 /* >>chng 06 aug 17, should go to numLevels_local instead of _max. */ 00650 for(n=ipH1s; n<iso.numLevels_local[ipISO][nelem]-1; ++n ) 00651 { 00652 fprintf(ioQQQ,"%li\t%li\t%g\t%g\t%g\t%g\tT=\t%g\t%g\t%g\t%g\n", nelem,n, 00653 iso.gamnc[ipISO][nelem][n] *iso.PopLTE[ipISO][nelem][n], 00654 /* induced recom, intergral times hlte */ 00655 /*iso.RadRecomb[ipISO][nelem][n][ipRecRad]+iso.rinduc[ipISO][nelem][n] ,*/ 00656 /* >>chng 02 jan 22, remove rinduc, replace with RecomInducRate */ 00657 iso.RadRecomb[ipISO][nelem][n][ipRecRad]+ 00658 iso.RecomInducRate[ipISO][nelem][n]*iso.PopLTE[ipISO][nelem][n] , 00659 /* induced rec */ 00660 iso.RecomInducRate[ipISO][nelem][n]*iso.PopLTE[ipISO][nelem][n] , 00661 /* spontaneous recombination */ 00662 iso.RadRecomb[ipISO][nelem][n][ipRecRad] , 00663 /* heating, followed by two processes that must balance it */ 00664 iso.PhotoHeat[ipISO][nelem][n]*iso.PopLTE[ipISO][nelem][n], 00665 iso.RecomInducCool_Coef[ipISO][nelem][n]*iso.PopLTE[ipISO][nelem][n]+SaveRadRecCool[n] , 00666 /* induced rec cooling, integral times hlte */ 00667 iso.RecomInducCool_Coef[ipISO][nelem][n]*iso.PopLTE[ipISO][nelem][n] , 00668 /* free-bound cooling per unit vol */ 00669 SaveRadRecCool[n] ); 00670 } 00671 } 00672 else 00673 { 00674 /* >>chng 06 aug 17, should go to numLevels_local instead of _max. */ 00675 for(n=ipH1s; n<iso.numLevels_local[ipISO][nelem]-1; ++n ) 00676 { 00677 fprintf(ioQQQ,"%li\t%li\t%g\t%g\t%g\t%g\tT=\t%g\t%g\t%g\t%g\n", nelem,n, 00678 iso.gamnc[ipISO][nelem][n] *dense.xIonDense[nelem][nelem+1-ipISO]*StatesElem[ipISO][nelem][n].Pop, 00679 /* induced recom, intergral times hlte */ 00680 iso.RadRecomb[ipISO][nelem][n][ipRecRad]*edenIonAbund+ 00681 iso.RecomInducRate[ipISO][nelem][n]*iso.PopLTE[ipISO][nelem][n] *edenIonAbund , 00682 iso.RadRecomb[ipISO][nelem][n][ipRecRad]*edenIonAbund , 00683 iso.RecomInducRate[ipISO][nelem][n]*iso.PopLTE[ipISO][nelem][n] *edenIonAbund , 00684 /* heating, followed by two processes that must balance it */ 00685 SavePhotoHeat[n], 00686 SaveInducCool[n]+SaveRadRecCool[n]*edenIonAbund , 00687 /* induced rec cooling, integral times hlte */ 00688 SaveInducCool[n] , 00689 /* free-bound cooling per unit vol */ 00690 SaveRadRecCool[n]*edenIonAbund ); 00691 } 00692 } 00693 } 00694 } 00695 return; 00696 } 00697 #if defined(__HP_aCC) 00698 #pragma OPTIMIZE OFF 00699 #pragma OPTIMIZE ON 00700 #endif