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_solve main routine to call iso_level and determine iso level balances */ 00004 /* HydroRenorm - renormalize H so that it agrees with the chemistry */ 00005 /* AGN_He1_CS routine to punch table needed for AGN3 - collision strengths of HeI */ 00006 #include "cddefines.h" 00007 #include "atmdat.h" 00008 #include "conv.h" 00009 #include "dense.h" 00010 #include "opacity.h" 00011 #include "elementnames.h" 00012 #include "h2.h" 00013 #include "helike.h" 00014 #include "helike_cs.h" 00015 #include "hmi.h" 00016 #include "hydrogenic.h" 00017 #include "ionbal.h" 00018 #include "iso.h" 00019 #include "phycon.h" 00020 #include "secondaries.h" 00021 #include "taulines.h" 00022 #include "thermal.h" 00023 #include "trace.h" 00024 00025 /* iso_solve main routine to call iso_level and determine iso level balances */ 00026 void iso_solve(long ipISO) 00027 { 00028 long int ipLo, 00029 ipHi, 00030 nelem, 00031 mol, 00032 lowsav=-1, 00033 ihisav=-1; 00034 double coltot, 00035 gamtot, 00036 sum, 00037 error; 00038 static bool lgFinitePop[LIMELM]; 00039 static bool lgMustInit[NISO]={true, true}; 00040 bool lgH_chem_conv; 00041 int loop_H_chem; 00042 double solomon_assumed; 00043 double *PumpSave=NULL; 00044 00045 DEBUG_ENTRY( "iso_solve()" ); 00046 00047 if( lgMustInit[ipISO] ) 00048 { 00049 for( nelem=ipHELIUM; nelem<LIMELM; ++nelem ) 00050 lgFinitePop[nelem] = true; 00051 } 00052 00053 /* 00054 * This is an option to account for intrinsic absorption or emission line the lyman 00055 * lines. This is done without changing the incident coarse continuum. 00056 * Check if continuum pumping of H Lyman lines to be multiplied by scale factor 00057 * hydro.xLymanPumpingScaleFactor is set with atom h-like Lyman pumping scale command 00058 * Lya pump rate is always updated - this is simplest way to finese intrinsic absorption 00059 * or emission 00060 */ 00061 if( ipISO==ipH_LIKE && hydro.xLymanPumpingScaleFactor!=1.f ) 00062 { 00063 /* >> chng 06 aug 31, from numLevels_max to _local */ 00064 PumpSave = (double *)MALLOC( (unsigned)iso.numLevels_local[ipISO][ipHYDROGEN]*sizeof(double) ); 00065 ipLo = 0; 00066 /* do not include the highest levels since pumping to them could create 00067 * artificial ionizations */ 00068 /* >> chng 06 aug 31, from numLevels_max to _local */ 00069 for( ipHi=2; ipHi < iso.numLevels_local[ipH_LIKE][ipHYDROGEN]; ++ipHi ) 00070 { 00071 PumpSave[ipHi] = Transitions[ipISO][ipHYDROGEN][ipHi][ipLo].Emis->pump; 00072 Transitions[ipISO][ipHYDROGEN][ipHi][ipLo].Emis->pump *= hydro.xLymanPumpingScaleFactor; 00073 } 00074 } 00075 00076 if( ipISO == ipHE_LIKE ) 00077 { 00078 /* save state of he-like low and high ionization, since must not change here, 00079 * since there is a parallel helium atom that must decrement he */ 00081 fixit(); /* why does this routine not control helium? Shouldn't it? */ 00082 // can just remove these entirely?, only applies to helium itself anyway. 00083 lowsav = dense.IonLow[ipHELIUM]; 00084 ihisav = dense.IonHigh[ipHELIUM]; 00085 } 00086 00087 for( nelem=ipISO; nelem < LIMELM; nelem++ ) 00088 { 00089 /* do not consider elements that have been turned off */ 00090 if( dense.lgElmtOn[nelem] ) 00091 { 00092 /* note that nelem scale is totally on c not physical scale, so 0 is h */ 00093 /* evaluate balance if ionization reaches this high */ 00094 if( (dense.IonHigh[nelem] >= nelem + 1 - ipISO) && 00095 (dense.IonLow[nelem] <= nelem + 1 - ipISO) ) 00096 { 00097 /* truncate atom if physical conditions limit the maximum principal quantum number of a 00098 * bound electron to a number less than the malloc'd size */ 00099 iso_continuum_lower( ipISO, nelem ); 00100 00101 /* evaluate collisional rates */ 00102 iso_collide( ipISO, nelem ); 00103 00104 /* evaluate photoionization rates */ 00105 iso_photo( ipISO , nelem ); 00106 00107 fixit(); /* the order of error calculation is screwed up... rationalize this */ 00108 /* Generate Gaussian errors if turned on. */ 00109 if( iso.lgRandErrGen[ipISO] && nzone==0 && !iso.lgErrGenDone[ipISO][nelem] ) 00110 { 00111 iso_error_generation(ipISO, nelem ); 00112 } 00113 00114 /* evaluate recombination rates */ 00115 iso_radiative_recomb( ipISO , nelem ); 00116 00117 if( opac.lgRedoStatic ) 00118 { 00119 if( nelem<=ipHELIUM ) 00120 { 00121 iso_collapsed_bnl_set( ipISO, nelem ); 00122 00123 //iso_collapsed_bnl_print( ipISO, nelem ); 00124 00125 iso_collapsed_Aul_update( ipISO, nelem ); 00126 00127 iso_collapsed_lifetimes_update( ipISO, nelem ); 00128 00129 iso_cascade( ipISO, nelem ); 00130 00131 iso_radiative_recomb_effective( ipISO, nelem ); 00132 } 00133 else 00134 { 00135 iso_cascade( ipISO, nelem ); 00136 00137 iso_radiative_recomb_effective( ipISO, nelem ); 00138 } 00139 } 00140 00141 /* evaluate state specific creation and destruction processes, 00142 * also define iso.xIonSimple */ 00143 iso_ionize_recombine( ipISO , nelem ); 00144 00145 /* solve for the level populations */ 00146 iso_level( ipISO , nelem ); 00147 00148 /* this contains a bunch of trace statements and HydroRenorm. 00149 * Will eventually come over to iso treatment. */ 00150 if( ipISO == ipH_LIKE ) 00151 HydroLevel(nelem); 00152 00153 /* say that we have set the populations */ 00154 lgFinitePop[nelem] = true; 00155 00156 } 00157 else if( lgFinitePop[nelem] ) 00158 { 00159 /* this branch, pops were set previously, but now are zero, 00160 * so must zero them out this time */ 00161 lgFinitePop[nelem] = false; 00162 00163 iso.pop_ion_ov_neut[ipISO][nelem] = 0.; 00164 iso.xIonSimple[ipISO][nelem] = 0.; 00165 00166 /* zero it out since no population*/ 00167 StatesElem[ipISO][nelem][0].Pop = 0.; 00168 for( ipHi=1; ipHi < iso.numLevels_max[ipISO][nelem]; ipHi++ ) 00169 { 00170 StatesElem[ipISO][nelem][ipHi].Pop = 0.; 00171 for( ipLo=0; ipLo < ipHi; ipLo++ ) 00172 { 00173 if( Transitions[ipISO][nelem][ipHi][ipLo].Emis->Aul <= iso.SmallA ) 00174 continue; 00175 00176 /* population of lower level rel to ion, corrected for stim em */ 00177 Transitions[ipISO][nelem][ipHi][ipLo].Emis->PopOpc = 0.; 00178 } 00179 } 00180 } 00181 /* option to force ionization */ 00182 if( dense.lgSetIoniz[nelem] ) 00183 { 00184 dense.xIonDense[nelem][nelem+1-ipISO] = dense.SetIoniz[nelem][nelem+1-ipISO]*dense.gas_phase[nelem]; 00185 dense.xIonDense[nelem][nelem-ipISO] = dense.SetIoniz[nelem][nelem-ipISO]*dense.gas_phase[nelem]; 00186 if( dense.SetIoniz[nelem][nelem+1-ipISO] > SMALLFLOAT ) 00187 { 00188 StatesElem[ipISO][nelem][ipH1s].Pop = dense.SetIoniz[nelem][nelem-ipISO] / dense.SetIoniz[nelem][nelem+1-ipISO]; 00189 } 00190 else 00191 { 00192 StatesElem[ipISO][nelem][ipH1s].Pop = 0.; 00193 } 00194 ASSERT( Transitions[ipISO][nelem][iso.nLyaLevel[ipISO]][0].Lo->Pop == StatesElem[ipISO][nelem][ipH1s].Pop ); 00195 } 00196 } 00197 } 00198 00199 /* setting this flag false means that we have been through this loop at least once already. */ 00200 lgMustInit[ipISO] = false; 00201 00202 if( ipISO == ipHE_LIKE ) 00203 { 00204 dense.IonLow[ipHELIUM] = lowsav; 00205 dense.IonHigh[ipHELIUM] = ihisav; 00206 } 00207 00208 if( ipISO==ipH_LIKE ) 00209 { 00210 /* ============================================================================== */ 00211 /* rest is for hydrogen only */ 00212 00213 /* this block appears redundant and could be removed? */ 00214 /* >> 02 nov 21 rjrw -- xIonDense is used in hmole but only set in ion_solver, 00215 * so we need this at least for the first iteration. */ 00216 00217 /* do molecular balance 00218 * hmovh1 will be ratio of molecular to atomic hydrogen 00219 * HIonFrac is fraction of H that is ionized, ratio of ion to atom */ 00220 00221 lgH_chem_conv = false; 00222 loop_H_chem = 0; 00223 while( loop_H_chem < 5 && !lgH_chem_conv ) 00224 { 00225 /* save solomon rate that was assumed in the above - want to make sure that assumed 00226 * solomon rate is stable, and does not change after call to large H2 molecule */ 00227 solomon_assumed = hmi.H2_Solomon_dissoc_rate_used_H2g/SDIV(hmi.H2_rate_destroy); 00228 00229 /* now do chem, this will reset hmi.H2_Solomon_dissoc_rate_used */ 00230 hmole(); 00231 /* now do level populations for H2 */ 00232 H2_LevelPops(); 00233 /* >>chng 05 feb 12 - add logic checking that Solomon rate from big molecule is equal to 00234 * rate used in H chem */ 00235 lgH_chem_conv = true; 00236 /* check whether solomon rate has changed */ 00237 if( h2.lgH2ON && hmi.lgBigH2_evaluated && hmi.lgH2_Chemistry_BigH2 ) 00238 { 00239 /* >>chng 05 mar 11, go from "total" to H2g for solomon rate */ 00240 if( fabs( solomon_assumed - hmi.H2_Solomon_dissoc_rate_BigH2_H2g/SDIV(hmi.H2_rate_destroy) ) > 00241 conv.EdenErrorAllowed/5.) 00242 { 00243 lgH_chem_conv = false; 00244 } 00245 } 00246 ++loop_H_chem; 00247 } 00248 00249 { 00250 /*@-redef@*/ 00251 /* often the H- route is the most efficient formation mechanism for H2, 00252 * will be through rate called ratach 00253 * this debug print statement is to trace h2 oscillations */ 00254 enum {DEBUG_LOC=false}; 00255 /*@+redef@*/ 00256 if(DEBUG_LOC ) 00257 { 00258 fprintf(ioQQQ,"DEBUG \t%.2f\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\n", 00259 fnzone, 00260 hmi.H2_total , 00261 hmi.Hmolec[ipMH2g], 00262 dense.xIonDense[ipHYDROGEN][0], 00263 dense.xIonDense[ipHYDROGEN][1], 00264 hmi.H2_Solomon_dissoc_rate_used_H2g, 00265 hmi.H2_Solomon_dissoc_rate_BD96_H2g, 00266 hmi.H2_Solomon_dissoc_rate_TH85_H2g); 00267 } 00268 } 00269 /* >>chng 01 may 09, add option to force abundance, with element name ioniz cmmnd */ 00270 if( dense.lgSetIoniz[ipHYDROGEN] ) 00271 { 00272 dense.xIonDense[ipHYDROGEN][1] = dense.SetIoniz[ipHYDROGEN][1]*dense.gas_phase[ipHYDROGEN]; 00273 dense.xIonDense[ipHYDROGEN][0] = dense.SetIoniz[ipHYDROGEN][0]*dense.gas_phase[ipHYDROGEN]; 00274 00275 /* initialize ground state pop too */ 00276 StatesElem[ipISO][ipHYDROGEN][ipH1s].Pop = dense.xIonDense[ipHYDROGEN][0]/dense.xIonDense[ipHYDROGEN][1]; 00277 } 00278 else 00279 { 00280 /* 00281 * >> chng 03 jan 15 rjrw:- terms are now in ion_solver, to allow for 00282 * molecular sources and sinks of H and H+. ion_solver renormalizes 00283 * to keep the total H abundance correct -- only the molecular 00284 * network is allowed to change this. 00285 */ 00286 ion_solver( ipHYDROGEN , false ); 00287 } 00288 00289 /* confirm that species still add up correctly */ 00290 /* this exposes a "leak" that occurs somewhere, almost certainly in hmole 00291 * if DEBUG_LOC is set true in the following there will be comments printing 00292 * due to a loss of some protons */ 00293 00294 sum = dense.xIonDense[ipHYDROGEN][0] + dense.xIonDense[ipHYDROGEN][1]; 00295 for(mol=0;mol<N_H_MOLEC;mol++) 00296 { 00297 sum += hmi.Hmolec[mol]*hmi.nProton[mol]; 00298 } 00299 /* do not double count H0 and H+ */ 00300 sum -= hmi.Hmolec[ipMH]+hmi.Hmolec[ipMHp]; 00301 00302 /* >>chng 06 jun 30, remove H in CO network, as much as a few percent can be 00303 * in the form of OH, H2O, etc 00304 * >>chng 06 jul 03, undo this change so function is as before - need to 00305 * think about co protons */ 00306 error = ( dense.gas_phase[ipHYDROGEN] - sum )/dense.gas_phase[ipHYDROGEN]; 00307 /*error = ( dense.gas_phase[ipHYDROGEN] - sum - dense.H_sum_in_CO)/dense.gas_phase[ipHYDROGEN];*/ 00308 { 00309 /*@-redef@*/ 00310 /* often the H- route is the most efficient formation mechanism for H2, 00311 * will be through rate called ratach 00312 * this debug print statement is to trace h2 oscillations */ 00313 enum {DEBUG_LOC=false}; 00314 /*@+redef@*/ 00315 if(DEBUG_LOC && (fabs(error) > 1e-4) ) 00316 fprintf(ioQQQ,"PROBLEM hydrogenic zone %li hden %.4e, sum %.4e (h-s)/h %.3e \n", nzone, dense.gas_phase[ipHYDROGEN] , sum , 00317 error ); 00318 } 00319 00320 fixit(); /* this is called in HydroLevel above, is it needed in both places? */ 00321 /* >>hcng 05 mar 24, 00322 * renormalize the populations and emission of H atom to agree with chemistry */ 00323 HydroRenorm(); 00324 00325 /* this is test whether collisions are important first define ratio of exits 00326 * from ground due to coll, rel to everthing else, then flag is large */ 00327 if( iso.PopLTE[ipH_LIKE][ipHYDROGEN][ipH1s] > 1e-30 ) 00328 { 00329 coltot = 00330 Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].Coll.ColUL* 00331 iso.PopLTE[ipH_LIKE][ipHYDROGEN][ipH2p]/iso.PopLTE[ipH_LIKE][ipHYDROGEN][ipH1s]* 00332 dense.eden*(iso.gamnc[ipH_LIKE][ipHYDROGEN][ipH2p] + iso.ColIoniz[ipH_LIKE][ipHYDROGEN][ipH2p] + 00333 Transitions[ipH_LIKE][ipHYDROGEN][ipH3s][ipH2p].Coll.ColUL*iso.PopLTE[ipH_LIKE][ipHYDROGEN][ipH3s]/iso.PopLTE[ipH_LIKE][ipHYDROGEN][ipH2p] + 00334 Transitions[ipH_LIKE][ipHYDROGEN][ipH3p][ipH2p].Coll.ColUL*iso.PopLTE[ipH_LIKE][ipHYDROGEN][ipH3p]/iso.PopLTE[ipH_LIKE][ipHYDROGEN][ipH2p] + 00335 Transitions[ipH_LIKE][ipHYDROGEN][ipH3d][ipH2p].Coll.ColUL*iso.PopLTE[ipH_LIKE][ipHYDROGEN][ipH3d]/iso.PopLTE[ipH_LIKE][ipHYDROGEN][ipH2p] )/ 00336 (Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].Coll.ColUL*dense.eden + 00337 Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].Emis->Aul* 00338 (Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].Emis->Pesc + 00339 Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].Emis->Pelec_esc + 00340 Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].Emis->Pdest) ); 00341 00342 /* caution that collisions are important (and so only small changes in 00343 * temperature should happen) if more than 20% of the total */ 00344 if( coltot > 0.2 ) 00345 { 00346 hydro.lgHColionImp = true; 00347 } 00348 } 00349 else 00350 { 00351 hydro.lgHColionImp = false; 00352 } 00353 00354 /* remember the ratio of pops of 2p to 1s for possible printout in prtComment 00355 * and to obtain Lya excitation temps. the pop of ground is not defined if 00356 * NO ionization at all since these pops are relative to ion */ 00357 /* >>chng 99 jun 03, added MAX2 to protect against totally neutral gas */ 00358 if( StatesElem[ipH_LIKE][ipHYDROGEN][ipH2p].Pop/MAX2(SMALLDOUBLE,StatesElem[ipH_LIKE][ipHYDROGEN][ipH1s].Pop) > 0.1 && 00359 StatesElem[ipH_LIKE][ipHYDROGEN][ipH1s].Pop > SMALLDOUBLE ) 00360 { 00361 hydro.lgHiPop2 = true; 00362 hydro.pop2mx = (realnum)MAX2(StatesElem[ipH_LIKE][ipHYDROGEN][ipH2p].Pop/StatesElem[ipH_LIKE][ipHYDROGEN][ipH1s].Pop, 00363 hydro.pop2mx); 00364 } 00365 00366 gamtot = iso.gamnc[ipH_LIKE][ipHYDROGEN][ipH1s] + secondaries.csupra[ipHYDROGEN][0]; 00367 00368 coltot = iso.ColIoniz[ipH_LIKE][ipHYDROGEN][ipH1s] + 00369 Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].Coll.ColUL*4.* 00370 iso.Boltzmann[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s]; 00371 00372 /* if ground state destruction rate is significant, recall different dest procceses */ 00373 if( iso.RateLevel2Cont[ipH_LIKE][ipHYDROGEN][ipH1s] > SMALLFLOAT ) 00374 { 00375 hydro.H_ion_frac_photo = 00376 (realnum)(iso.gamnc[ipH_LIKE][ipHYDROGEN][ipH1s]/iso.RateLevel2Cont[ipH_LIKE][ipHYDROGEN][ipH1s] ); 00377 00378 /* fraction of ionizations of H from ground, due to thermal collisions */ 00379 hydro.H_ion_frac_collis = 00380 (realnum)(iso.ColIoniz[ipH_LIKE][ipHYDROGEN][ipH1s]*dense.eden/iso.RateLevel2Cont[ipH_LIKE][ipHYDROGEN][ipH1s]); 00381 00382 /* this flag is used in ConvBase to decide whether we 00383 * really need to converge the secondary ionization rates */ 00384 secondaries.sec2total = 00385 (realnum)(secondaries.csupra[ipHYDROGEN][0] / iso.RateLevel2Cont[ipH_LIKE][ipHYDROGEN][ipH1s]); 00386 00387 /* frac of ionizations due to ct */ 00388 atmdat.HIonFrac = atmdat.HCharExcIonTotal / iso.RateLevel2Cont[ipH_LIKE][ipHYDROGEN][ipH1s]; 00389 } 00390 else 00391 { 00392 hydro.H_ion_frac_collis = 0.; 00393 hydro.H_ion_frac_photo = 0.; 00394 secondaries.sec2total = 0.; 00395 atmdat.HIonFrac = 0.; 00396 } 00397 00398 if( trace.lgTrace ) 00399 { 00400 fprintf( ioQQQ, " Hydrogenic return %.2f ",fnzone); 00401 fprintf(ioQQQ,"H0:%.3e ", dense.xIonDense[ipHYDROGEN][0]); 00402 fprintf(ioQQQ,"H+:%.3e ", dense.xIonDense[ipHYDROGEN][1]); 00403 fprintf(ioQQQ,"H2:%.3e ", hmi.H2_total); 00404 fprintf(ioQQQ,"H-:%.3e ", hmi.Hmolec[ipMHm]); 00405 fprintf(ioQQQ,"ne:%.3e ", dense.eden); 00406 fprintf( ioQQQ, " REC, COL, GAMT= "); 00407 /* recomb rate coef, cm^3 s-1 */ 00408 fprintf(ioQQQ,"%.2e ", iso.RadRec_effec[ipH_LIKE][ipHYDROGEN] ); 00409 fprintf(ioQQQ,"%.2e ", coltot); 00410 fprintf(ioQQQ,"%.2e ", gamtot); 00411 fprintf( ioQQQ, " CSUP="); 00412 PrintE82( ioQQQ, secondaries.csupra[ipHYDROGEN][0]); 00413 fprintf( ioQQQ, "\n"); 00414 } 00415 00416 if( hydro.xLymanPumpingScaleFactor!=1.f ) 00417 { 00418 ipLo = 0; 00419 /* do not include the highest levels since pumping to them could create 00420 * artificial ionizaitons */ 00421 /* >> chng 06 aug 31, from numLevels_max to _local */ 00422 for( ipHi=2; ipHi < iso.numLevels_local[ipH_LIKE][ipHYDROGEN]; ++ipHi ) 00423 { 00424 Transitions[ipH_LIKE][ipHYDROGEN][ipHi][ipLo].Emis->pump = PumpSave[ipHi]; 00425 } 00426 free(PumpSave); 00427 } 00428 } 00429 return; 00430 } 00431 00432 /* HydroRenorm - renormalize H so that it agrees with the chemistry */ 00433 void HydroRenorm( void ) 00434 { 00435 double sum_atom_iso , renorm; 00436 long int level, 00437 nelem, 00438 ipHi , 00439 ipLo; 00440 00441 DEBUG_ENTRY( "HydroRenorm()" ); 00442 00443 /*>>chng 04 mar 23, add this renorm */ 00444 /* renormalize the state specific populations, so that they 00445 * add up to the results that came from ion_solver 00446 * units at first is sum div by H+ density, since Pop2Ion defn this way */ 00447 sum_atom_iso = 0.; 00448 /* >> chng 06 aug 31, from numLevels_max to _local */ 00449 for( level=ipH1s; level < iso.numLevels_local[ipH_LIKE][ipHYDROGEN]; level++ ) 00450 { 00451 sum_atom_iso += StatesElem[ipH_LIKE][ipHYDROGEN][level].Pop; 00452 } 00453 /* now convert to cm-3 - this is total population in iso solved model */ 00454 sum_atom_iso *= dense.xIonDense[ipHYDROGEN][1]; 00455 00456 /* >>chng 04 may 25, humunculus sum_atom_iso is zero */ 00457 if( sum_atom_iso > SMALLFLOAT ) 00458 { 00459 renorm = dense.xIonDense[ipHYDROGEN][0] / sum_atom_iso; 00460 } 00461 else 00462 { 00463 renorm = 0.; 00464 } 00465 ASSERT( renorm < BIGFLOAT ); 00466 /*fprintf(ioQQQ,"DEBUG renorm\t%.2f\t%.3e\n",fnzone, renorm);*/ 00467 /* renormalize populations from iso-model atom so that they agree with ion solver & chemistry */ 00468 StatesElem[ipH_LIKE][ipHYDROGEN][ipH1s].Pop *= renorm; 00469 /*fprintf(ioQQQ,"DEBUG h \t%.3e hydrogenic renorm %.3e\n", 00470 StatesElem[ipH_LIKE][ipHYDROGEN][ipH1s].Pop , 00471 StatesElem[ipH_LIKE][ipHYDROGEN][ipH1s].Pop/renorm );*/ 00472 00473 nelem = ipHYDROGEN; 00474 /* >> chng 06 aug 31, from numLevels_max to _local */ 00475 for( ipHi=ipH2s; ipHi < iso.numLevels_local[ipH_LIKE][nelem]; ipHi++ ) 00476 { 00477 StatesElem[ipH_LIKE][ipHYDROGEN][ipHi].Pop *= renorm; 00478 00479 for( ipLo=ipH1s; ipLo < ipHi; ipLo++ ) 00480 { 00481 if( Transitions[ipH_LIKE][nelem][ipHi][ipLo].Emis->Aul <= iso.SmallA ) 00482 continue; 00483 00484 /* population of lower level rel to ion, corrected for stim em */ 00485 Transitions[ipH_LIKE][nelem][ipHi][ipLo].Emis->PopOpc *= renorm; 00486 } 00487 } 00488 00489 return; 00490 } 00491 00492 /*iso_prt_pops print out iso sequence populations or departure coefficients */ 00493 void iso_prt_pops( long ipISO, long nelem, bool lgPrtDeparCoef ) 00494 { 00495 long int in, il, is, i, ipLo, nResolved, ipFirstCollapsed=LONG_MIN; 00496 char chPrtType[2][12]={"populations","departure"}; 00497 /* first dimension is multiplicity */ 00498 char chSpin[3][9]= {"singlets", "doublets", "triplets"}; 00499 00500 #define ITEM_TO_PRINT(A_) ( lgPrtDeparCoef ? iso.DepartCoef[ipISO][nelem][A_] : StatesElem[ipISO][nelem][A_].Pop ) 00501 00502 DEBUG_ENTRY( "iso_prt_pops()" ); 00503 00504 ASSERT( ipISO < NISO ); 00505 00506 for( is = 1; is<=3; ++is) 00507 { 00508 if( ipISO == ipH_LIKE && is != 2 ) 00509 continue; 00510 else if( ipISO == ipHE_LIKE && is != 1 && is != 3 ) 00511 continue; 00512 00513 ipFirstCollapsed= iso.numLevels_local[ipISO][nelem]-iso.nCollapsed_local[ipISO][nelem]; 00514 nResolved = StatesElem[ipISO][nelem][ipFirstCollapsed-1].n; 00515 ASSERT( nResolved == iso.n_HighestResolved_local[ipISO][nelem] ); 00516 ASSERT(nResolved > 0 ); 00517 00518 /* give element number and spin */ 00519 fprintf(ioQQQ," %s %s %s %s\n", 00520 iso.chISO[ipISO], 00521 elementnames.chElementSym[nelem], 00522 chSpin[is-1], 00523 chPrtType[lgPrtDeparCoef]); 00524 00525 /* header with the l states */ 00526 fprintf(ioQQQ," n\\l=> "); 00527 for( i =0; i < nResolved; ++i) 00528 { 00529 fprintf(ioQQQ,"%2ld ",i); 00530 } 00531 fprintf(ioQQQ,"\n"); 00532 00533 /* loop over prin quant numbers, one per line, with l across */ 00534 for( in = 1; in <= nResolved; ++in) 00535 { 00536 if( is==3 && in==1 ) 00537 continue; 00538 00539 fprintf(ioQQQ," %2ld ",in); 00540 00541 for( il = 0; il < in; ++il) 00542 { 00543 if( ipISO==ipHE_LIKE && (in==2) && (il==1) && (is==3) ) 00544 { 00545 fprintf( ioQQQ, "%9.3e ", ITEM_TO_PRINT(ipHe2p3P0) ); 00546 fprintf( ioQQQ, "%9.3e ", ITEM_TO_PRINT(ipHe2p3P1) ); 00547 fprintf( ioQQQ, "%9.3e ", ITEM_TO_PRINT(ipHe2p3P2) ); 00548 } 00549 else 00550 { 00551 ipLo = iso.QuantumNumbers2Index[ipISO][nelem][in][il][is]; 00552 fprintf( ioQQQ, "%9.3e ", ITEM_TO_PRINT(ipLo) ); 00553 } 00554 } 00555 fprintf(ioQQQ,"\n"); 00556 } 00557 } 00558 /* above loop was over spin, now do collapsed levels, no spin or ang momen */ 00559 for( il = ipFirstCollapsed; il < iso.numLevels_local[ipISO][nelem]; ++il) 00560 { 00561 in = StatesElem[ipISO][nelem][il].n; 00562 /* prin quan number of collapsed levels */ 00563 fprintf(ioQQQ," %2ld ",in); 00564 fprintf( ioQQQ, "%9.3e ", ITEM_TO_PRINT(il) ); 00565 fprintf(ioQQQ,"\n"); 00566 } 00567 return; 00568 } 00569 00570 /* routine to punch table needed for AGN3 - collision strengths of HeI */ 00571 void AGN_He1_CS( FILE *ioPun ) 00572 { 00573 00574 long int i; 00575 00576 /* list of temperatures where cs will be printed */ 00577 const int NTE = 5; 00578 double TeList[NTE] = {6000.,10000.,15000.,20000.,25000.}; 00579 double TempSave; 00580 00581 DEBUG_ENTRY( "AGN_He1_CS()" ); 00582 00583 /* put on a header */ 00584 fprintf(ioPun, "Te\t2 3s 33s\n"); 00585 00586 /* Restore the original temp when this routine done. */ 00587 TempSave = phycon.te; 00588 00589 for( i=0; i<NTE; ++i ) 00590 { 00591 TempChange(TeList[i] , false); 00592 00593 fprintf(ioPun , "%.0f\t", 00594 TeList[i] ); 00595 fprintf(ioPun , "%.2f\t", 00596 HeCSInterp( 1 , ipHe3s3S , ipHe2s3S, ipELECTRON ) ); 00597 fprintf(ioPun , "%.2f\t", 00598 HeCSInterp( 1 , ipHe3p3P , ipHe2s3S, ipELECTRON ) ); 00599 fprintf(ioPun , "%.2f\t", 00600 HeCSInterp( 1 , ipHe3d3D , ipHe2s3S, ipELECTRON ) ); 00601 fprintf(ioPun , "%.3f\t", 00602 HeCSInterp( 1 , ipHe3d1D , ipHe2s3S, ipELECTRON ) ); 00603 /*fprintf(ioPun , "%.1f\t%.1f\t%.1f\n", */ 00604 fprintf(ioPun , "%.1f\n", 00605 HeCSInterp( 1 , ipHe2p3P0 , ipHe2s3S, ipELECTRON ) + 00606 HeCSInterp( 1 , ipHe2p3P1 , ipHe2s3S, ipELECTRON ) + 00607 HeCSInterp( 1 , ipHe2p3P2 , ipHe2s3S, ipELECTRON )); 00608 } 00609 00610 /* no need to force update since didn't do above */ 00611 TempChange(TempSave , false); 00612 return; 00613 }