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 /*HydroLevel calls iso_level to solve for ionization balance 00004 * and level populations of model hydrogen atom */ 00005 /*PrtHydroTrace1 print trace info for hydrogen-like species */ 00006 #include "cddefines.h" 00007 #include "taulines.h" 00008 #include "iso.h" 00009 #include "dense.h" 00010 #include "secondaries.h" 00011 #include "trace.h" 00012 #include "phycon.h" 00013 #include "ionbal.h" 00014 #include "hydrogenic.h" 00015 00016 /*PrtHydroTrace1a print trace info for hydrogen-like species */ 00017 STATIC void PrtHydroTrace1a(long nelem ) 00018 { 00019 double colfrc, 00020 phtfrc, 00021 secfrc, 00022 collider; 00023 00024 DEBUG_ENTRY( "PrtHydroTrace1a()" ); 00025 00026 /* >>chng 05 aug 17, must use real electron density for collisional ionization of H 00027 * since in Leiden v4 pdr with its artificial high temperature coll ion can be important 00028 * H on H is homonuclear and scaling laws for other elements does not apply 00029 * next two were multiplied by dense.EdenHCorr and changed to dense.eden for H, 00030 * EdenHCorr for rest */ 00031 if( nelem==ipHYDROGEN ) 00032 { 00033 /* special version for H0 onto H0 */ 00034 collider = dense.EdenHontoHCorr; 00035 } 00036 else 00037 { 00038 collider = dense.EdenHCorr; 00039 } 00040 00041 /* identify how atom is ionized for full trace */ 00042 if( iso.xIonSimple[ipH_LIKE][nelem] > 0. ) 00043 { 00044 /* fraction of ionization due to photoionization */ 00045 phtfrc = iso.gamnc[ipH_LIKE][nelem][ipH1s]/((dense.eden*(iso.RadRec_effec[ipH_LIKE][nelem] + 00046 ionbal.CotaRate[nelem]) )* 00047 iso.xIonSimple[ipH_LIKE][nelem]); 00048 00049 /* fraction of ionization due to collisional ionization 00050 * >>chng 05 aug 17 from EdenHCorr to collider */ 00051 colfrc = (iso.ColIoniz[ipH_LIKE][nelem][ipH1s]*collider )/ 00052 ((dense.eden*(iso.RadRec_effec[ipH_LIKE][nelem] + 00053 ionbal.CotaRate[0]) )* 00054 iso.xIonSimple[ipH_LIKE][nelem]); 00055 00056 /* fraction of ionization due to secondary ionization */ 00057 secfrc = secondaries.csupra[nelem][nelem]/((dense.eden*(iso.RadRec_effec[ipH_LIKE][nelem] + 00058 ionbal.CotaRate[0]) )* 00059 iso.xIonSimple[ipH_LIKE][nelem]); 00060 } 00061 else 00062 { 00063 phtfrc = 0.; 00064 colfrc = 0.; 00065 secfrc = 0.; 00066 } 00067 00068 fprintf( ioQQQ, " HydroLevel Z=%2ld called, simple II/I=",nelem); 00069 PrintE93( ioQQQ, iso.xIonSimple[ipH_LIKE][nelem]); 00070 fprintf( ioQQQ," PhotFrc:"); 00071 PrintE93( ioQQQ,phtfrc); 00072 fprintf(ioQQQ," ColFrc:"); 00073 PrintE93( ioQQQ,colfrc); 00074 fprintf( ioQQQ," SecFrc"); 00075 PrintE93( ioQQQ, secfrc); 00076 fprintf( ioQQQ," Te:"); 00077 PrintE93( ioQQQ,phycon.te); 00078 fprintf( ioQQQ," eden:"); 00079 PrintE93( ioQQQ,dense.eden); 00080 fprintf( ioQQQ,"\n"); 00081 return; 00082 } 00083 00084 /*PrtHydroTrace1 print trace info for hydrogen-like species */ 00085 STATIC void PrtHydroTrace1(long nelem ) 00086 { 00087 long int ipHi , ipLo , i; 00088 00089 DEBUG_ENTRY( "PrtHydroTrace1()" ); 00090 00091 fprintf( ioQQQ, 00092 " HydroLevel%3ld finds arrays, with optical depths defined? %li induced 2ph=%12.3e\n", 00093 nelem, iteration, Transitions[ipH_LIKE][nelem][ipH2s][ipH1s].Emis->pump ); 00094 /* 06 aug 28, from numLevels_max to _local. */ 00095 for( ipHi=ipH2s; ipHi < iso.numLevels_local[ipH_LIKE][nelem]; ipHi++ ) 00096 { 00097 fprintf( ioQQQ, "up:%2ld", ipHi ); 00098 fprintf( ioQQQ, "lo" ); 00099 for( ipLo=ipH1s; ipLo < ipHi; ipLo++ ) 00100 { 00101 fprintf( ioQQQ, "%9ld", ipLo ); 00102 } 00103 fprintf( ioQQQ, "\n" ); 00104 00105 fprintf( ioQQQ, "%3ld", ipHi ); 00106 fprintf( ioQQQ, " A*esc" ); 00107 for( ipLo=ipH1s; ipLo < ipHi; ipLo++ ) 00108 { 00109 fprintf( ioQQQ,PrintEfmt("%9.2e", Transitions[ipH_LIKE][nelem][ipHi][ipLo].Emis->Aul* 00110 Transitions[ipH_LIKE][nelem][ipHi][ipLo].Emis->Pesc )); 00111 } 00112 fprintf( ioQQQ, "\n" ); 00113 00114 fprintf( ioQQQ, "%3ld", ipHi ); 00115 fprintf( ioQQQ, " A*ees" ); 00116 for( ipLo=ipH1s; ipLo < ipHi; ipLo++ ) 00117 { 00118 fprintf( ioQQQ,PrintEfmt("%9.2e", Transitions[ipH_LIKE][nelem][ipHi][ipLo].Emis->Aul* 00119 Transitions[ipH_LIKE][nelem][ipHi][ipLo].Emis->Pelec_esc )); 00120 } 00121 fprintf( ioQQQ, "\n" ); 00122 00123 fprintf( ioQQQ, "%3ld", ipHi ); 00124 fprintf( ioQQQ, " tauin" ); 00125 for( ipLo=ipH1s; ipLo < ipHi; ipLo++ ) 00126 { 00127 fprintf( ioQQQ,PrintEfmt("%9.2e", Transitions[ipH_LIKE][nelem][ipHi][ipLo].Emis->TauIn )); 00128 } 00129 fprintf( ioQQQ, "\n" ); 00130 00131 fprintf( ioQQQ, "%3ld", ipHi ); 00132 fprintf( ioQQQ, " t tot" ); 00133 for( ipLo=ipH1s; ipLo < ipHi; ipLo++ ) 00134 { 00135 fprintf( ioQQQ,PrintEfmt("%9.2e", Transitions[ipH_LIKE][nelem][ipHi][ipLo].Emis->TauTot )); 00136 } 00137 fprintf( ioQQQ, "\n" ); 00138 00139 fprintf( ioQQQ, "%3ld", ipHi ); 00140 fprintf( ioQQQ, " Esc " ); 00141 for( ipLo=ipH1s; ipLo < ipHi; ipLo++ ) 00142 { 00143 fprintf( ioQQQ,PrintEfmt("%9.2e", Transitions[ipH_LIKE][nelem][ipHi][ipLo].Emis->Pesc )); 00144 } 00145 fprintf( ioQQQ, "\n" ); 00146 00147 fprintf( ioQQQ, "%3ld", ipHi ); 00148 fprintf( ioQQQ, " Eesc " ); 00149 for( ipLo=ipH1s; ipLo < ipHi; ipLo++ ) 00150 { 00151 fprintf( ioQQQ,PrintEfmt("%9.2e", Transitions[ipH_LIKE][nelem][ipHi][ipLo].Emis->Pelec_esc )); 00152 } 00153 fprintf( ioQQQ, "\n" ); 00154 00155 fprintf( ioQQQ, "%3ld", ipHi ); 00156 fprintf( ioQQQ, " Dest " ); 00157 for( ipLo=ipH1s; ipLo < ipHi; ipLo++ ) 00158 { 00159 fprintf( ioQQQ,PrintEfmt("%9.2e", Transitions[ipH_LIKE][nelem][ipHi][ipLo].Emis->Pdest) ); 00160 } 00161 fprintf( ioQQQ, "\n" ); 00162 00163 fprintf( ioQQQ, "%3ld", ipHi ); 00164 fprintf( ioQQQ, " A*dst" ); 00165 for( ipLo=ipH1s; ipLo < ipHi; ipLo++ ) 00166 { 00167 fprintf( ioQQQ,PrintEfmt("%9.2e", Transitions[ipH_LIKE][nelem][ipHi][ipLo].Emis->Aul* 00168 Transitions[ipH_LIKE][nelem][ipHi][ipLo].Emis->Pdest )); 00169 } 00170 fprintf( ioQQQ, "\n" ); 00171 00172 fprintf( ioQQQ, "%3ld", ipHi ); 00173 fprintf( ioQQQ, " StrkE" ); 00174 for( ipLo=ipH1s; ipLo < ipHi; ipLo++ ) 00175 { 00176 fprintf( ioQQQ,PrintEfmt("%9.2e", iso.pestrk[ipH_LIKE][nelem][ipLo][ipHi] )); 00177 } 00178 fprintf( ioQQQ, "\n" ); 00179 00180 fprintf( ioQQQ, "%3ld", ipHi ); 00181 fprintf( ioQQQ, " B(ul)" ); 00182 for( ipLo=ipH1s; ipLo < ipHi; ipLo++ ) 00183 { 00184 fprintf( ioQQQ,PrintEfmt("%9.2e", Transitions[ipH_LIKE][nelem][ipHi][ipLo].Emis->pump* 00185 StatesElem[ipH_LIKE][nelem][ipLo].g/StatesElem[ipH_LIKE][nelem][ipHi].g )); 00186 } 00187 fprintf( ioQQQ, "\n" ); 00188 00189 fprintf( ioQQQ, "%3ld", ipHi ); 00190 fprintf( ioQQQ, " tcont" ); 00191 for( ipLo=ipH1s; ipLo < ipHi; ipLo++ ) 00192 { 00193 fprintf( ioQQQ,PrintEfmt("%9.2e", Transitions[ipH_LIKE][nelem][ipHi][ipLo].Emis->TauCon )); 00194 } 00195 fprintf( ioQQQ, "\n" ); 00196 00197 fprintf( ioQQQ, "%3ld", ipHi ); 00198 fprintf( ioQQQ, " C(ul)" ); 00199 for( ipLo=ipH1s; ipLo < ipHi; ipLo++ ) 00200 { 00201 fprintf( ioQQQ,PrintEfmt("%9.2e", 00202 Transitions[ipH_LIKE][nelem][ipHi][ipLo].Coll.ColUL*dense.eden )); 00203 } 00204 fprintf( ioQQQ, "\n" ); 00205 00206 if( ipHi == 2 ) 00207 { 00208 fprintf( ioQQQ, " FeIIo"); 00209 fprintf( ioQQQ,PrintEfmt("%9.2e", 00210 hydro.dstfe2lya* Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].Emis->Aul )); 00211 fprintf( ioQQQ, "\n"); 00212 } 00213 } 00214 00215 fprintf( ioQQQ, " " ); 00216 /* 06 aug 28, from numLevels_max to _local. */ 00217 for( i=1; i < iso.numLevels_local[ipH_LIKE][nelem]; i++ ) 00218 { 00219 fprintf( ioQQQ, "%9ld", i ); 00220 } 00221 fprintf( ioQQQ, "\n" ); 00222 return; 00223 } 00224 00225 /*HydroLevel calls iso_level to solve for ionization balance 00226 * and level populations of model hydrogen atom */ 00227 void HydroLevel(long int nelem) 00228 { 00229 long int i; 00230 double collider; 00231 00232 int ipISO = ipH_LIKE; 00233 00234 DEBUG_ENTRY( "HydroLevel()" ); 00235 00236 /* check that we were called with valid charge */ 00237 ASSERT( nelem >= 0); 00238 ASSERT( nelem < LIMELM ); 00239 00240 /* >>chng 05 aug 17, must use real electron density for collisional ionization of H 00241 * since in Leiden v4 pdr with its artificial high temperature coll ion can be important 00242 * H on H is homonuclear and scaling laws for other elements does not apply 00243 * next two were multiplied by dense.EdenHCorr and changed to dense.eden for H, 00244 * EdenHCorr for rest */ 00245 if( ipISO == ipH_LIKE && nelem==ipHYDROGEN ) 00246 { 00247 /* special version for H0 onto H0 */ 00248 collider = dense.EdenHontoHCorr; 00249 } 00250 else 00251 { 00252 collider = dense.EdenHCorr; 00253 } 00254 00255 /* option to print some rates */ 00256 if( (trace.lgTrace && trace.lgIsoTraceFull[ipISO]) && (nelem == trace.ipIsoTrace[ipISO]) ) 00257 PrtHydroTrace1(nelem); 00258 00259 if( trace.lgHBug && trace.lgTrace ) 00260 PrtHydroTrace1a(nelem); 00261 00262 /* this is main trace h-like printout */ 00263 if( (trace.lgIsoTraceFull[ipISO] && trace.lgTrace) && (nelem == trace.ipIsoTrace[ipISO]) ) 00264 { 00265 fprintf( ioQQQ, " HLEV HGAMNC" ); 00266 PrintE93( ioQQQ, iso.gamnc[ipISO][nelem][ipH1s] ); 00267 /* 06 aug 28, from numLevels_max to _local. */ 00268 for( i=ipH2s; i < iso.numLevels_local[ipISO][nelem]; i++ ) 00269 { 00270 fprintf(ioQQQ,PrintEfmt("%9.2e", iso.gamnc[ipISO][nelem][i] )); 00271 } 00272 fprintf( ioQQQ, "\n" ); 00273 00274 fprintf( ioQQQ, " HLEV TOTCAP" ); 00275 /* 06 aug 28, from numLevels_max to _local. */ 00276 for( i=1; i < iso.numLevels_local[ipISO][nelem]; i++ ) 00277 { 00278 fprintf(ioQQQ,PrintEfmt("%9.2e", iso.RateCont2Level[ipISO][nelem][i]/dense.eden )); 00279 } 00280 fprintf( ioQQQ," tot"); 00281 fprintf( ioQQQ,PrintEfmt("%10.2e", ionbal.RateRecomTot[nelem][nelem-ipISO]/dense.eden ) ); 00282 fprintf( ioQQQ, "\n" ); 00283 00284 fprintf( ioQQQ, " HLEV IND Rc" ); 00285 /* 06 aug 28, from numLevels_max to _local. */ 00286 for( i=ipH1s; i < iso.numLevels_local[ipISO][nelem]; i++ ) 00287 { 00288 fprintf(ioQQQ,PrintEfmt("%9.2e", iso.RecomInducRate[ipISO][nelem][i]*iso.PopLTE[ipISO][nelem][i] )); 00289 } 00290 fprintf( ioQQQ, "\n" ); 00291 00292 /* incuded recombination rate coefficients */ 00293 fprintf( ioQQQ, " IND Rc LTE " ); 00294 /* 06 aug 28, from numLevels_max to _local. */ 00295 for( i=ipH1s; i < iso.numLevels_local[ipISO][nelem]; i++ ) 00296 { 00297 fprintf(ioQQQ,PrintEfmt("%9.2e", 00298 iso.gamnc[ipISO][nelem][i]*iso.PopLTE[ipISO][nelem][i] )); 00299 } 00300 fprintf( ioQQQ, "\n" ); 00301 00302 /* LTE level populations */ 00303 fprintf( ioQQQ, " HLEV HLTE" ); 00304 /* 06 aug 28, from numLevels_max to _local. */ 00305 for( i=ipH1s; i < iso.numLevels_local[ipISO][nelem]; i++ ) 00306 { 00307 fprintf(ioQQQ,PrintEfmt("%9.2e", iso.PopLTE[ipISO][nelem][i] )); 00308 } 00309 fprintf( ioQQQ, "\n" ); 00310 00311 /* fraction of total ionization due to collisions from given level */ 00312 fprintf( ioQQQ, " HLEVfr cion" ); 00313 /* 06 aug 28, from numLevels_max to _local. */ 00314 for( i=ipH1s; i < iso.numLevels_local[ipISO][nelem]; i++ ) 00315 { 00316 fprintf(ioQQQ,PrintEfmt("%9.2e", 00317 iso.ColIoniz[ipISO][nelem][i]* 00318 StatesElem[ipISO][nelem][i].Pop*collider/MAX2(1e-30,iso.RateLevel2Cont[ipISO][nelem][i]) ) ); 00319 } 00320 fprintf( ioQQQ, "\n" ); 00321 00322 /* fraction of total ionization due to photoionization from given level */ 00323 if( ionbal.RateRecomTot[nelem][nelem]> 0. ) 00324 { 00325 fprintf( ioQQQ, " HLEVfrPhIon" ); 00326 /* 06 aug 28, from numLevels_max to _local. */ 00327 for( i=ipH1s; i < iso.numLevels_local[ipISO][nelem]; i++ ) 00328 { 00329 fprintf(ioQQQ,PrintEfmt("%9.2e", 00330 iso.gamnc[ipISO][nelem][i]*StatesElem[ipISO][ipHYDROGEN][i].Pop/MAX2(1e-30,iso.RateLevel2Cont[ipISO][nelem][i]) ) ); 00331 } 00332 fprintf( ioQQQ, "\n" ); 00333 } 00334 00335 fprintf( ioQQQ, " HLEV HN" ); 00336 /* 06 aug 28, from numLevels_max to _local. */ 00337 for( i=ipH1s; i < iso.numLevels_local[ipISO][nelem]; i++ ) 00338 { 00339 fprintf(ioQQQ,PrintEfmt("%9.2e", StatesElem[ipISO][nelem][i].Pop )); 00340 } 00341 fprintf( ioQQQ, "\n" ); 00342 00343 fprintf( ioQQQ, " HLEV b(n)" ); 00344 /* 06 aug 28, from numLevels_max to _local. */ 00345 for( i=ipH1s; i < iso.numLevels_local[ipISO][nelem]; i++ ) 00346 { 00347 fprintf(ioQQQ,PrintEfmt("%9.2e", iso.DepartCoef[ipISO][nelem][i] )); 00348 } 00349 fprintf( ioQQQ, "\n" ); 00350 00351 fprintf( ioQQQ, " HLEV X12tot"); 00352 fprintf(ioQQQ,PrintEfmt("%9.2e", secondaries.x12tot )); 00353 fprintf( ioQQQ," Grn dest:"); 00354 fprintf(ioQQQ,PrintEfmt("%9.2e", 00355 ionbal.RateIonizTot[nelem][nelem] )); 00356 fprintf(ioQQQ, "\n"); 00357 } 00358 00359 00360 /* >>chng 05 mar 24, 00361 * renormalize the populations and emission of H atom to agree with chemistry 00362 * these were not being kept parallel with chemistry, and caused large changes in O+ 00363 * abundance when finally done */ 00364 fixit(); /* this probably does not need to be called here. */ 00365 if( nelem == ipHYDROGEN ) 00366 HydroRenorm(); 00367 00368 if( trace.lgTrace ) 00369 { 00370 /* iso.RecomTotal[nelem] is gross rec coef, computed here while filling in matrix 00371 * elements, all physical processes included. 00372 * RadRec_effec is total effective radiative only */ 00373 fprintf( ioQQQ, " HydroLevel%3ld return %s te=", 00374 nelem, 00375 iso.chTypeAtomUsed[ipISO][nelem] ); 00376 PrintE93( ioQQQ,phycon.te); 00377 fprintf( ioQQQ," ion/atom=%.4e",iso.pop_ion_ov_neut[ipISO][nelem]); 00378 00379 fprintf( ioQQQ," simple=%.4e",iso.xIonSimple[ipISO][nelem]); 00380 00381 fprintf( ioQQQ," b1=%.2e",iso.DepartCoef[ipISO][nelem][ipH1s]); 00382 00383 fprintf( ioQQQ," ion rate=%.4e",ionbal.RateIonizTot[nelem][nelem-ipISO]); 00384 00385 fprintf( ioQQQ," TotRec=%.4e",ionbal.RateRecomTot[nelem][nelem-ipISO]); 00386 00387 fprintf( ioQQQ," RadRec=%.4e",iso.RadRec_effec[ipISO][nelem]); 00388 fprintf( ioQQQ, "\n"); 00389 } 00390 return; 00391 }