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_level solve for iso-sequence level populations */ 00004 #include "cddefines.h" 00005 #include "atmdat.h" 00006 #include "continuum.h" 00007 #include "conv.h" 00008 #include "dense.h" 00009 #include "dynamics.h" 00010 #include "elementnames.h" 00011 #include "grainvar.h" 00012 #include "he.h" 00013 #include "helike.h" 00014 #include "hydrogenic.h" 00015 #include "ionbal.h" 00016 #include "iso.h" 00017 #include "mole.h" 00018 #include "phycon.h" 00019 #include "physconst.h" 00020 #include "rfield.h" 00021 #include "secondaries.h" 00022 #include "taulines.h" 00023 #include "thirdparty.h" 00024 #include "trace.h" 00025 00026 /* solve for level populations */ 00027 void iso_level( const long int ipISO, const long int nelem) 00028 { 00029 long int ipHi, 00030 ipLo, 00031 i, 00032 level, 00033 level_error; 00034 00035 const long int numlevels_local = iso.numLevels_local[ipISO][nelem]; 00036 00037 double BigError; 00038 00039 int32 nerror; 00040 double HighestPColOut = 0., 00041 collider, 00042 sum_popn_ov_ion, 00043 TooSmall; 00044 bool lgNegPop=false; 00045 valarray<int32> ipiv(numlevels_local); 00046 /* this block of variables will be obtained and freed here */ 00047 valarray<double> 00048 creation(numlevels_local), 00049 error(numlevels_local), 00050 work(numlevels_local), 00051 Save_creation(numlevels_local); 00052 double source=0., 00053 sink=0.; 00054 valarray<double> PopPerN(iso.n_HighestResolved_local[ipISO][nelem]+1); 00055 00056 multi_arr<double,2,C_TYPE> z, SaveZ; 00057 00058 DEBUG_ENTRY( "iso_level()" ); 00059 00060 /* check that we were called with valid charge */ 00061 ASSERT( nelem >= ipISO ); 00062 ASSERT( nelem < LIMELM ); 00063 00064 /* now do the 2D array */ 00065 z.alloc(numlevels_local,numlevels_local); 00066 00067 /* >>chng 05 aug 17, must use real electron density for collisional ionization of H 00068 * since in Leiden v4 pdr with its artificial high temperature coll ion can be important 00069 * H on H is homonuclear and scaling laws for other elements does not apply 00070 * next two were multiplied by dense.EdenHCorr and changed to dense.eden for H, 00071 * EdenHCorr for rest */ 00072 if( ipISO == ipH_LIKE && nelem==ipHYDROGEN ) 00073 { 00074 /* special version for H0 onto H0 */ 00075 collider = dense.EdenHontoHCorr; 00076 } 00077 else 00078 { 00079 collider = dense.EdenHCorr; 00080 } 00081 00082 /* fill in recombination vector - values were set in iso_ionize_recombine.cpp */ 00083 for( level=0; level < numlevels_local; level++ ) 00084 { 00085 /* total recombination to level [s-1] */ 00086 creation[level] = iso.RateCont2Level[ipISO][nelem][level]; 00087 } 00088 00089 /* these two collision rates must be the same or we are in big trouble, 00090 * since used interchangably */ 00091 ASSERT( ionbal.CollIonRate_Ground[nelem][nelem-ipISO][0]< SMALLFLOAT || 00092 fabs( iso.ColIoniz[ipISO][nelem][0]* collider / 00093 SDIV(ionbal.CollIonRate_Ground[nelem][nelem-ipISO][0] ) - 1.) < 0.001 ); 00094 00095 00096 if( ipISO == ipH_LIKE ) 00097 { 00098 if( nelem==ipHYDROGEN ) 00099 TooSmall = 1e-28; 00100 else if( nelem==ipHELIUM ) 00101 TooSmall = 1e-18; 00102 else 00103 TooSmall = 2e-18; 00104 } 00105 else if( ipISO == ipHE_LIKE ) 00106 { 00107 /* >>chng 03 apr 30, different limit for He itself, 00108 * since so important in molecular regions */ 00109 if( nelem==ipHELIUM ) 00110 TooSmall = 1e-20; 00111 else 00112 TooSmall = 1e-10; 00113 } 00114 else 00115 TotalInsanity(); 00116 00117 /* which case atom to solve??? */ 00118 /* >>chng 03 apr 11, from ae-30 bail, to ae-35, will catch it in hydro 2 low */ 00119 if( ipISO == ipH_LIKE && 00120 ( iso.xIonSimple[ipISO][nelem] < 1e-35 ) ) 00121 { 00122 /* don't bother if no ionizing radiation */ 00123 strcpy( iso.chTypeAtomUsed[ipISO][nelem], "zero " ); 00124 if( trace.lgTrace && (nelem == trace.ipIsoTrace[ipISO]) ) 00125 { 00126 fprintf( ioQQQ, " iso_level iso=%2ld nelem=%2ld simple II/I=%10.2e so not doing equilibrium, doing %s.\n", 00127 ipISO, nelem, iso.xIonSimple[ipISO][nelem], iso.chTypeAtomUsed[ipISO][nelem] ); 00128 } 00129 00130 /* okay to leave as iso.numLevels_max */ 00131 for( i=0; i < iso.numLevels_max[ipISO][nelem]; i++ ) 00132 { 00133 iso.DepartCoef[ipISO][nelem][i] = 1.; 00134 StatesElem[ipISO][nelem][i].Pop = 0.; 00135 } 00136 00137 ionbal.RateRecomTot[nelem][nelem-ipISO] = 0.; 00138 iso.xIonSimple[ipISO][nelem] = 0.; 00139 iso.pop_ion_ov_neut[ipISO][nelem] = 0.; 00140 iso.qTot2S[ipISO][nelem] = 0.; 00141 } 00142 /* >>chng 99 nov 23, very dense model close to lte had very low simple 00143 * ionization fraction but moderate ionization since all pops in 00144 * excited states - added test for density 00145 * second change - for all cases use this routine if ionization is 00146 * very low indeed */ 00147 /* >>chng 00 aug 26, added matrix force option, 00148 * logic is to override default options, "none was set in zero */ 00149 00150 /* NB - this test is mostly bogus since chTypeAtom is "POPU" in zero */ 00151 00152 /* >>chng 02 mar 13, iso.chTypeAtomSet had been set to POPU in zero, look for 00153 * comment with this date. This had the effect of killing this following test. 00154 * This also meant that the code had been happily inverting these impossible 00155 * matrices for quite some time. This test changed to stringest one, in 00156 * light of this */ 00157 else if( ipISO == ipH_LIKE && 00158 ((strcmp( iso.chTypeAtomSet[ipISO] , "LOWT" )==0) || 00159 (iso.xIonSimple[ipISO][nelem] < TooSmall) ) ) 00160 { 00161 strcpy( iso.chTypeAtomUsed[ipISO][nelem], "Low T" ); 00162 if( trace.lgTrace && (nelem == trace.ipIsoTrace[ipISO]) ) 00163 { 00164 fprintf( ioQQQ, " iso_level iso=%2ld nelem=%2ld simple II/I=%10.2e so not doing equilibrium, doing %s\n", 00165 ipISO, nelem, iso.xIonSimple[ipISO][nelem], iso.chTypeAtomUsed[ipISO][nelem] ); 00166 } 00167 00168 /* the ionization ratio WILL be equal to iso.xIonSimple[ipISO][nelem] */ 00169 iso.pop_ion_ov_neut[ipISO][nelem] = iso.xIonSimple[ipISO][nelem]; 00170 00171 /* do simple cascade, should always work, 00172 * but not very well at high densities */ 00173 HydroT2Low( ipISO, nelem ); 00174 iso.qTot2S[ipISO][nelem] = 0.; 00175 } 00176 /* branch for very little ionization, use simple estimate but do not do level populations 00177 * never set to zero, use simple two-level system, since He+ important for chemistry */ 00178 else if( ipISO == ipHE_LIKE && 00179 ((strcmp( iso.chTypeAtomSet[ipISO] , "LOWT" )==0) || 00180 (iso.xIonSimple[ipISO][nelem] < TooSmall) ) ) 00181 { 00182 strcpy( iso.chTypeAtomUsed[ipISO][nelem], "Low T" ); 00183 /* total ionization will just be the ground state */ 00184 ionbal.RateIonizTot[nelem][nelem-ipISO] = iso.RateLevel2Cont[ipISO][nelem][0]; 00185 /* >>chng 01 nov 23, increase upper limit by 2 */ 00186 lgNegPop = false; 00187 /* >>chng 02 apr 10, had set to zero, 00188 * we will leave the level populations at zero but define an ion ratio, 00189 * since He+ abundance cannot go to zero in a molecular cloud, 00190 * - He+ charge transfer is the main CO destruction mechanism */ 00191 iso.pop_ion_ov_neut[ipISO][nelem] = iso.xIonSimple[ipISO][nelem]; 00192 StatesElem[ipISO][nelem][0].Pop = 1./SDIV(iso.pop_ion_ov_neut[ipISO][nelem]); 00193 for( long n=1; n < numlevels_local; n++ ) 00194 { 00195 StatesElem[ipISO][nelem][n].Pop = 0.; 00196 } 00197 iso.qTot2S[ipISO][nelem] = 0.; 00198 } 00199 else 00200 { 00201 ASSERT( NISO == 2 ); 00202 long SpinUsed[NISO] = {2,1}; 00203 long indexNmaxP = 00204 iso.QuantumNumbers2Index[ipISO][nelem][ iso.n_HighestResolved_local[ipISO][nelem] ][1][SpinUsed[ipISO]]; 00205 00206 strcpy( iso.chTypeAtomUsed[ipISO][nelem], "Popul" ); 00207 if( trace.lgTrace && (nelem == trace.ipIsoTrace[ipISO]) ) 00208 { 00209 fprintf( ioQQQ, " iso_level iso=%2ld nelem=%2ld doing regular matrix inversion, %s\n", 00210 ipISO, nelem, iso.chTypeAtomUsed[ipISO][nelem] ); 00211 } 00212 00213 multi_arr<quantumState,3>::const_iterator StElm = StatesElem.begin(ipISO,nelem); 00214 00215 /* master balance equation, use when significant population */ 00216 for( level=0; level < numlevels_local; level++ ) 00217 { 00218 /* all process depopulating level and placing into the continuum 00219 * this does NOT include grain charge transfer ionization, added below */ 00220 z[level][level] = iso.RateLevel2Cont[ipISO][nelem][level]; 00221 00222 if( level == 1 ) 00223 /* >>chng 05 dec 21, rm eden to make into rate coefficient */ 00224 iso.qTot2S[ipISO][nelem] = iso.ColIoniz[ipISO][nelem][level]; 00225 00226 multi_arr<transition,4>::const_iterator Trans = Transitions.begin(ipISO,nelem,level); 00227 md4ci Boltz = iso.Boltzmann.begin(ipISO,nelem,level); 00228 00229 /* all processes populating level from below */ 00230 for( ipLo=0; ipLo < level; ipLo++ ) 00231 { 00232 double coll_down = Trans[ipLo].Coll.ColUL * collider; 00233 00234 double RadDecay = MAX2( iso.SmallA, Trans[ipLo].Emis->Aul* 00235 (Trans[ipLo].Emis->Pesc + 00236 Trans[ipLo].Emis->Pelec_esc + 00237 Trans[ipLo].Emis->Pdest)* 00238 KILL_BELOW_PLASMA(Trans[ipLo].EnergyWN*WAVNRYD) ); 00239 00240 double pump = MAX2( iso.SmallA, Trans[ipLo].Emis->pump * 00241 KILL_BELOW_PLASMA(Trans[ipLo].EnergyWN*WAVNRYD) ); 00242 00243 if( iso.lgRandErrGen[ipISO] ) 00244 { 00245 coll_down *= iso.ErrorFactor[ipISO][nelem][level][ipLo][IPCOLLIS]; 00246 RadDecay *= iso.ErrorFactor[ipISO][nelem][level][ipLo][IPRAD]; 00247 pump *= iso.ErrorFactor[ipISO][nelem][level][ipLo][IPRAD]; 00248 } 00249 00250 double coll_up = coll_down * 00251 (double)StElm[level].g/ 00252 (double)StElm[ipLo].g* 00253 Boltz[ipLo]; 00254 00255 z[ipLo][ipLo] += coll_up + pump ; 00256 z[ipLo][level] = - ( coll_up + pump ); 00257 00258 double pump_down = pump * 00259 (double)StElm[ipLo].g/ 00260 (double)StElm[level].g; 00261 00262 z[level][level] += RadDecay + pump_down + coll_down; 00263 z[level][ipLo] = - (RadDecay + pump_down + coll_down); 00264 00265 if( level == indexNmaxP ) 00266 { 00267 HighestPColOut += coll_down; 00268 } 00269 if( ipLo == indexNmaxP ) 00270 { 00271 HighestPColOut += coll_up; 00272 } 00273 00274 /* find total collisions out of 2^3S to singlets. */ 00275 if( (level == 1) && (ipLo==0) ) 00276 { 00277 iso.qTot2S[ipISO][nelem] += coll_down/dense.EdenHCorr; 00278 } 00279 if( (ipLo == 1) && (ipISO==ipH_LIKE || StElm[level].S==1) ) 00280 { 00281 iso.qTot2S[ipISO][nelem] += coll_up/dense.EdenHCorr; 00282 } 00283 } 00284 } 00285 00286 if( ipISO == nelem ) 00287 { 00288 /* iso.lgCritDensLMix[ipISO] is a flag used to print warning if density is 00289 * too low for first collapsed level to be l-mixed. Check is if l-mixing collisions 00290 * out of highest resolved singlet P are greater than sum of transition probs out. */ 00291 if( HighestPColOut < 1./StatesElem[ipISO][nelem][indexNmaxP].lifetime ) 00292 { 00293 iso.lgCritDensLMix[ipISO] = false; 00294 } 00295 } 00296 00298 ASSERT( ipISO <= ipHE_LIKE ); 00299 00300 /* induced two photon emission - special because upward and downward are 00301 * not related by ratio of statistical weights */ 00302 /* iso.lgInd2nu_On is controlled with SET IND2 ON/OFF command */ 00303 z[1+ipISO][0] -= iso.TwoNu_induc_dn[ipISO][nelem]*iso.lgInd2nu_On; 00304 z[0][1+ipISO] -= iso.TwoNu_induc_up[ipISO][nelem]*iso.lgInd2nu_On; 00305 00306 /* rates out of 1s, and out of 2s */ 00307 z[0][0] += iso.TwoNu_induc_up[ipISO][nelem]*iso.lgInd2nu_On; 00308 z[1+ipISO][1+ipISO] += iso.TwoNu_induc_dn[ipISO][nelem]*iso.lgInd2nu_On; 00309 00310 /* grain charge transfer recombination and ionization to ALL other stages */ 00311 for( long ion=0; ion<=nelem+1; ++ion ) 00312 { 00313 if( ion!=nelem-ipISO ) 00314 { 00315 /* recombination must be multiplied by a ratio of densities to get proper rate. */ 00316 source += gv.GrainChTrRate[nelem][ion][nelem-ipISO] * 00317 dense.xIonDense[nelem][ion] ; 00318 /* take ionization out of every level. */ 00319 sink += gv.GrainChTrRate[nelem][nelem-ipISO][ion]; 00320 } 00321 } 00322 00323 /* >>chng 02 Sep 06 rjrw -- all elements have these terms */ 00324 /*>>>chng 02 oct 01, only include if lgAdvection is set */ 00325 if( dynamics.lgAdvection && dynamics.lgISO[ipISO]) 00326 { 00327 /* add in advection - these terms normally zero */ 00328 /* assume for present that all advection is into ground state */ 00329 source += dynamics.Source[nelem][nelem-ipISO]; 00330 /* >>chng 02 Sep 06 rjrw -- advective term not recombination */ 00331 /* can sink from all components (must do, for conservation) */ 00332 sink += dynamics.Rate; 00333 } 00334 00335 #if 0 00336 /* add in source and sink terms from molecular network. */ 00337 CodeReview(); /* Check... */ 00338 if( nelem == ipISO && conv.nTotalIoniz ) 00339 { 00340 /* these are the external source and sink terms */ 00341 /* source first */ 00342 source += mole.source[nelem][nelem-ipISO]; 00343 sink += mole.sink[nelem][nelem-ipISO]; 00344 } 00345 #endif 00346 00347 creation[0] += source/SDIV(dense.xIonDense[nelem][nelem+1-ipISO]); 00348 for( level=0; level < numlevels_local; level++ ) 00349 { 00350 z[level][level] += sink; 00351 } 00352 00353 /* >>chng 04 nov 30, atom XX-like collisions off turns this off too */ 00354 if( secondaries.Hx12[ipISO][nelem][iso.nLyaLevel[ipISO]] * iso.lgColl_excite[ipISO] > 0. ) 00355 { 00356 /* now add on supra thermal excitation */ 00357 for( level=1; level < numlevels_local; level++ ) 00358 { 00359 double RateUp , RateDown; 00360 00361 RateUp = secondaries.Hx12[ipISO][nelem][level]; 00362 RateDown = RateUp * (double)StatesElem[ipISO][nelem][ipH1s].g / 00363 (double)StatesElem[ipISO][nelem][level].g; 00364 00365 /* total rate out of lower level */ 00366 z[ipH1s][ipH1s] += RateUp; 00367 00368 /* rate from the upper level to ground */ 00369 z[level][ipH1s] -= RateDown; 00370 00371 /* rate from ground to upper level */ 00372 z[ipH1s][level] -= RateUp; 00373 00374 z[level][level] += RateDown; 00375 } 00376 } 00377 00378 /* =================================================================== 00379 * 00380 * at this point all matrix elements have been established 00381 * 00382 * ==================================================================== */ 00383 00384 /* save matrix, this allocates SaveZ */ 00385 SaveZ = z; 00386 00387 for( ipLo=0; ipLo < numlevels_local; ipLo++ ) 00388 Save_creation[ipLo] = creation[ipLo]; 00389 00390 if( trace.lgTrace && trace.lgIsoTraceFull[ipISO] && (nelem == trace.ipIsoTrace[ipISO]) ) 00391 { 00392 fprintf( ioQQQ, " pop level others => (iso_level)\n" ); 00393 for( long n=0; n < numlevels_local; n++ ) 00394 { 00395 fprintf( ioQQQ, " %s %s %2ld", iso.chISO[ipISO], elementnames.chElementNameShort[nelem], n ); 00396 for( long j=0; j < numlevels_local; j++ ) 00397 { 00398 fprintf( ioQQQ,"\t%.9e", z[j][n] ); 00399 } 00400 fprintf( ioQQQ, "\n" ); 00401 } 00402 fprintf(ioQQQ," recomb "); 00403 for( long n=0; n < numlevels_local; n++ ) 00404 { 00405 fprintf( ioQQQ,"\t%.9e", creation[n] ); 00406 } 00407 fprintf( ioQQQ, "\n" ); 00408 fprintf(ioQQQ," recomb ct %.2e co %.2e hectr %.2e hctr %.2e\n", 00409 atmdat.HeCharExcRecTotal, 00410 findspecies("CO")->hevmol , 00411 atmdat.HeCharExcRecTo[nelem][nelem-ipISO]*dense.xIonDense[ipHELIUM][0] , 00412 atmdat.HCharExcRecTo[nelem][nelem-ipISO]*dense.xIonDense[ipHYDROGEN][0] ); 00413 } 00414 00415 nerror = 0; 00416 00417 getrf_wrapper(numlevels_local,numlevels_local, 00418 z.data(),numlevels_local,&ipiv[0],&nerror); 00419 00420 getrs_wrapper('N',numlevels_local,1,z.data(),numlevels_local,&ipiv[0], 00421 &creation[0],numlevels_local,&nerror); 00422 00423 if( nerror != 0 ) 00424 { 00425 fprintf( ioQQQ, " iso_level: dgetrs finds singular or ill-conditioned matrix\n" ); 00426 cdEXIT(EXIT_FAILURE); 00427 } 00428 00429 /* check whether solution is valid */ 00430 /* >>chng 06 aug 28, both of these from numLevels_max to _local. */ 00431 for( level=ipH1s; level < numlevels_local; level++ ) 00432 { 00433 double BigSoln= 0.; 00434 error[level] = 0.; 00435 for( long n=ipH1s; n < numlevels_local; n++ ) 00436 { 00437 /* remember the largest value of the soln matrix to div by below */ 00438 if( fabs(SaveZ[n][level] ) > BigSoln ) 00439 BigSoln = fabs(SaveZ[n][level]); 00440 00441 error[level] += SaveZ[n][level]*creation[n]; 00442 } 00443 00444 if( BigSoln > 0. ) 00445 { 00446 error[level] = (error[level] - Save_creation[level])/ BigSoln; 00447 } 00448 else 00449 { 00450 error[level] = 0.; 00451 } 00452 } 00453 00454 /* remember largest residual in matrix inversion */ 00455 BigError = -1.; 00456 level_error = -1; 00457 /* >>chng 06 aug 28, from numLevels_max to _local. */ 00458 for( level=ipH1s; level < numlevels_local; level++ ) 00459 { 00460 double abserror; 00461 abserror = fabs( error[level]); 00462 /* this will be the largest residual in the matrix inversion */ 00463 if( abserror > BigError ) 00464 { 00465 BigError = abserror; 00466 level_error = level; 00467 } 00468 } 00469 00470 /* matrix inversion should be nearly as good as the accuracy of a double, 00471 * but demand that it is better than epsilon for a float */ 00472 if( BigError > FLT_EPSILON ) 00473 { 00474 if( !conv.lgSearch ) 00475 fprintf(ioQQQ," PROBLEM" ); 00476 00477 fprintf(ioQQQ, 00478 " iso_level zone %.2f - largest residual in iso=%li %s nelem=%li matrix inversion is %g " 00479 "level was %li Search?%c \n", 00480 fnzone, 00481 ipISO, 00482 elementnames.chElementName[nelem], 00483 nelem , 00484 BigError , 00485 level_error, 00486 TorF(conv.lgSearch) ); 00487 } 00488 00489 /* put departure coefficients and level populations into master array */ 00490 for( level=0; level < numlevels_local; level++ ) 00491 { 00492 StatesElem[ipISO][nelem][level].Pop = creation[level]; 00493 00494 /* check for negative level populations */ 00495 if( StatesElem[ipISO][nelem][level].Pop < 0. ) 00496 lgNegPop = true; 00497 00498 if( StatesElem[ipISO][nelem][level].Pop <= 0 ) 00499 { 00500 fprintf(ioQQQ," non-positive level pop for iso = %li, nelem = %li = %s, level=%li val=%.3e\n", 00501 ipISO, 00502 nelem , 00503 elementnames.chElementSym[nelem], 00504 level, 00505 StatesElem[ipISO][nelem][level].Pop ); 00506 } 00507 00508 if( iso.PopLTE[ipISO][nelem][level] > 0. ) 00509 { 00510 /* the LTE departure coefficients */ 00511 iso.DepartCoef[ipISO][nelem][level] = 00512 StatesElem[ipISO][nelem][level].Pop/ 00513 (iso.PopLTE[ipISO][nelem][level]* dense.eden); 00514 } 00515 else 00516 { 00517 iso.DepartCoef[ipISO][nelem][level] = 0.; 00518 } 00519 } 00520 00521 /* zero populations of unused levels. */ 00522 for( level=numlevels_local; level < iso.numLevels_max[ipISO][nelem]; level++ ) 00523 { 00524 StatesElem[ipISO][nelem][level].Pop = 0.; 00525 /* >>chng 06 jul 25, no need to zero this out, fix limit to 3-body heating elsewhere. */ 00526 /* iso.PopLTE[ipISO][nelem][level] = 0.; */ 00527 } 00528 } 00529 00530 /* all solvers end up here */ 00531 /* sum_popn_ov_ion will become the ratio of iso to parent 00532 * species, create sum of level pops per ion first */ 00533 sum_popn_ov_ion = 0.; 00534 00535 /* this is total ionization of this species referenced to the ground state */ 00536 ionbal.RateIonizTot[nelem][nelem-ipISO] = 0.; 00537 00538 for( level=0; level < numlevels_local; level++ ) 00539 { 00540 /* sum of all ionization processes from this atom to ion */ 00541 ionbal.RateIonizTot[nelem][nelem-ipISO] += 00542 StatesElem[ipISO][nelem][level].Pop * iso.RateLevel2Cont[ipISO][nelem][level]; 00543 00544 /* create sum of populations relative to ion */ 00545 sum_popn_ov_ion += StatesElem[ipISO][nelem][level].Pop; 00546 } 00547 00548 /* convert back to scaled from ground */ 00549 if( ( ionbal.RateIonizTot[nelem][nelem-ipISO]/MAX2(SMALLDOUBLE , sum_popn_ov_ion) ) > BIGDOUBLE ) 00550 { 00551 fprintf( ioQQQ, "DISASTER RateIonizTot for Z=%li, ion %li is larger than BIGDOUBLE. This is a big problem.", 00552 nelem+1, nelem-ipISO); 00553 cdEXIT(EXIT_FAILURE); 00554 } 00555 else 00556 ionbal.RateIonizTot[nelem][nelem-ipISO] /= MAX2(SMALLDOUBLE , sum_popn_ov_ion); 00557 00558 /* this is sum of all populations in model, div by ion 00559 * create ratio of ion pops to total atom */ 00560 if( sum_popn_ov_ion >= 1e-30 ) 00561 { 00562 iso.pop_ion_ov_neut[ipISO][nelem] = 1./sum_popn_ov_ion; 00563 } 00564 else 00565 { 00566 iso.pop_ion_ov_neut[ipISO][nelem] = 0.; 00567 00568 /* >>chng 03 aug 16, zero out parent density - had not been done */ 00569 dense.xIonDense[nelem][nelem+1-ipISO] = 0.; 00570 00571 /* reset pointer to one lower stage of ionization so this not 00572 * considered again, hydrogenic considered if IonHigh is nelem+2 */ 00573 dense.IonHigh[nelem] = nelem; 00574 ionbal.RateIonizTot[nelem][nelem-ipISO] = 0.; 00575 00576 /* now zero this out */ 00577 /* okay to leave as iso.numLevels_max */ 00578 for( ipHi=0; ipHi < iso.numLevels_max[ipISO][nelem]; ipHi++ ) 00579 { 00580 /* population of level relative to ion */ 00581 StatesElem[ipISO][nelem][ipHi].Pop = 0.; 00582 } 00583 } 00584 00585 ASSERT( ionbal.RateIonizTot[nelem][nelem-ipISO] >= 0. ); 00586 00587 /* check on the sum of the populations */ 00588 if( lgNegPop || iso.pop_ion_ov_neut[ipISO][nelem] < 0. ) 00589 { 00590 fprintf( ioQQQ, 00591 " DISASTER iso_level finds negative ion fraction for iso=%2ld nelem=%2ld "\ 00592 "%s using solver %s, IonFrac=%.3e simple=%.3e TE=%.3e ZONE=%4ld\n", 00593 ipISO, 00594 nelem, 00595 elementnames.chElementSym[nelem], 00596 iso.chTypeAtomUsed[ipISO][nelem], 00597 iso.pop_ion_ov_neut[ipISO][nelem], 00598 iso.xIonSimple[ipISO][nelem], 00599 phycon.te, 00600 nzone ); 00601 00602 fprintf( ioQQQ, " level pop are:\n" ); 00603 for( i=0; i < numlevels_local; i++ ) 00604 { 00605 fprintf( ioQQQ,PrintEfmt("%8.1e", StatesElem[ipISO][nelem][i].Pop )); 00606 if( (i!=0) && !(i%10) ) fprintf( ioQQQ,"\n" ); 00607 } 00608 fprintf( ioQQQ, "\n" ); 00609 ContNegative(); 00610 ShowMe(); 00611 cdEXIT(EXIT_FAILURE); 00612 } 00613 00614 if( ipISO == ipHE_LIKE && nelem==ipHELIUM && nzone>0 ) 00615 { 00616 /* find fraction of He0 destructions due to photoionization of 2^3S */ 00617 double ratio; 00618 double rateOutOf2TripS = StatesElem[ipISO][nelem][ipHe2s3S].Pop * iso.RateLevel2Cont[ipISO][nelem][ipHe2s3S]; 00619 if( rateOutOf2TripS > SMALLFLOAT ) 00620 { 00621 ratio = rateOutOf2TripS / 00622 ( rateOutOf2TripS + StatesElem[ipISO][nelem][ipHe1s1S].Pop*ionbal.RateIonizTot[nelem][nelem-ipISO] ); 00623 } 00624 else 00625 { 00626 ratio = 0.; 00627 } 00628 if( ratio > he.frac_he0dest_23S ) 00629 { 00630 /* remember zone where this happended and fraction, and frac due to photoionization */ 00631 he.nzone = nzone; 00632 he.frac_he0dest_23S = ratio; 00633 rateOutOf2TripS = StatesElem[ipISO][nelem][ipHe2s3S].Pop *iso.gamnc[ipISO][nelem][ipHe2s3S]; 00634 if( rateOutOf2TripS > SMALLFLOAT ) 00635 { 00636 he.frac_he0dest_23S_photo = rateOutOf2TripS / 00637 ( rateOutOf2TripS + StatesElem[ipISO][nelem][ipHe1s1S].Pop*ionbal.RateIonizTot[nelem][nelem-ipISO] ); 00638 } 00639 else 00640 { 00641 he.frac_he0dest_23S_photo = 0.; 00642 } 00643 } 00644 } 00645 00646 for( ipHi=1; ipHi<numlevels_local; ++ipHi ) 00647 { 00648 for( ipLo=0; ipLo<ipHi; ++ipLo ) 00649 { 00650 if( Transitions[ipISO][nelem][ipHi][ipLo].Emis->Aul <= iso.SmallA ) 00651 continue; 00652 00653 /* population of lower level rel to ion, corrected for stim em */ 00654 Transitions[ipISO][nelem][ipHi][ipLo].Emis->PopOpc = 00655 StatesElem[ipISO][nelem][ipLo].Pop - StatesElem[ipISO][nelem][ipHi].Pop* 00656 StatesElem[ipISO][nelem][ipLo].g/StatesElem[ipISO][nelem][ipHi].g; 00657 00658 /* don't allow masers from collapsed levels */ 00659 if( N_(ipHi) > iso.n_HighestResolved_local[ipISO][nelem] ) 00660 Transitions[ipISO][nelem][ipHi][ipLo].Emis->PopOpc = StatesElem[ipISO][nelem][ipLo].Pop; 00661 } 00662 } 00663 00664 return; 00665 }