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_create create data for hydrogen and helium, 1 per coreload, called by ContCreatePointers 00004 * in turn called after commands parsed */ 00005 /*iso_zero zero data for hydrogen and helium */ 00006 #include "cddefines.h" 00007 #include "atmdat.h" 00008 #include "dense.h" 00009 #include "elementnames.h" 00010 #include "helike.h" 00011 #include "helike_einsta.h" 00012 #include "hydro_bauman.h" 00013 #include "hydrogenic.h" 00014 #include "hydroeinsta.h" 00015 #include "iso.h" 00016 #include "lines_service.h" 00017 #include "opacity.h" 00018 #include "phycon.h" 00019 #include "physconst.h" 00020 #include "secondaries.h" 00021 #include "taulines.h" 00022 #include "thirdparty.h" 00023 00024 /*iso_zero zero data for hydrogen and helium */ 00025 STATIC void iso_zero(void); 00026 00027 /* allocate memory for iso sequence structures */ 00028 STATIC void iso_allocate(void); 00029 00030 /* define levels of iso sequences and assign quantum numbers to those levels */ 00031 STATIC void iso_assign_quantum_numbers(void); 00032 00033 STATIC void FillExtraLymanLine( transition *t, long ipISO, long nelem, long nHi ); 00034 00035 /* calculate radiative lifetime of an individual iso state */ 00036 STATIC double iso_state_lifetime( long ipISO, long nelem, long n, long l ); 00037 00038 STATIC void iso_satellite( void ); 00039 00040 char chL[21]={'S','P','D','F','G','H','I','K','L','M','N','O','Q','R','T','U','V','W','X','Y','Z'}; 00041 00042 void iso_create(void) 00043 { 00044 long int ipHi, 00045 ipLo, 00046 ipISO, 00047 nelem; 00048 00049 static int nCalled = 0; 00050 00051 double HIonPoten; 00052 00053 DEBUG_ENTRY( "iso_create()" ); 00054 00055 /* > 1 if not first call, then just zero arrays out */ 00056 if( nCalled > 0 ) 00057 { 00058 iso_zero(); 00059 return; 00060 } 00061 00062 /* this is first call, increment the nCalled counterso never do this again */ 00063 ++nCalled; 00064 00065 /* these are the statistical weights of the ions */ 00066 iso.stat_ion[ipH_LIKE] = 1.f; 00067 iso.stat_ion[ipHE_LIKE] = 2.f; 00068 00069 /* this routine allocates all the memory 00070 * needed for the iso sequence structures */ 00071 iso_allocate(); 00072 00073 /* loop over iso sequences and assign quantum numbers to all levels */ 00074 iso_assign_quantum_numbers(); 00075 00076 /* this is a dummy line, junk it too. */ 00077 TransitionJunk( &TauDummy ); 00078 TauDummy.Hi = AddState2Stack(); 00079 TauDummy.Lo = AddState2Stack(); 00080 TauDummy.Emis = AddLine2Stack( true ); 00081 00082 /********************************************/ 00083 /********** Line and level energies ********/ 00084 /********************************************/ 00085 for( ipISO=ipH_LIKE; ipISO<NISO; ++ipISO ) 00086 { 00087 /* main hydrogenic arrays, fill with sane values */ 00088 for( nelem=ipISO; nelem < LIMELM; nelem++ ) 00089 { 00090 /* must always do helium even if turned off */ 00091 if( nelem < 2 || dense.lgElmtOn[nelem] ) 00092 { 00093 /* Dima's array has ionization potentials in eV, but not on same 00094 * scale as cloudy itself*/ 00095 /* extra factor accounts for this */ 00096 HIonPoten = t_ADfA::Inst().ph1(0,0,nelem,0)/EVRYD* 0.9998787; 00097 ASSERT(HIonPoten > 0.); 00098 00099 /* go from ground to the highest level */ 00100 for( ipHi=0; ipHi < iso.numLevels_max[ipISO][nelem]; ipHi++ ) 00101 { 00102 double EnergyWN, EnergyRyd; 00103 00104 if( ipISO == ipH_LIKE ) 00105 { 00106 EnergyRyd = HIonPoten/POW2((double)N_(ipHi)); 00107 } 00108 else if( ipISO == ipHE_LIKE ) 00109 { 00110 EnergyRyd = helike_energy( nelem, ipHi ) * WAVNRYD; 00111 } 00112 else 00113 { 00114 /* Other iso sequences don't exist yet. */ 00115 TotalInsanity(); 00116 } 00117 00118 /* >>chng 02 feb 09, change test to >= 0 since we now use 0 for 2s-2p */ 00119 ASSERT(EnergyRyd >= 0.); 00120 00121 iso.xIsoLevNIonRyd[ipISO][nelem][ipHi] = EnergyRyd; 00122 00123 /* now loop from ground to level below ipHi */ 00124 for( ipLo=0; ipLo < ipHi; ipLo++ ) 00125 { 00126 EnergyWN = RYD_INF * (iso.xIsoLevNIonRyd[ipISO][nelem][ipLo] - 00127 iso.xIsoLevNIonRyd[ipISO][nelem][ipHi]); 00128 00129 /* This is the minimum line energy we will allow. */ 00130 /* \todo 2 wire this to lowest energy of code. */ 00131 if( EnergyWN==0 && ipISO==ipHE_LIKE ) 00132 EnergyWN = 0.0001; 00133 00134 if( EnergyWN < 0. ) 00135 EnergyWN = -1.0 * EnergyWN; 00136 00137 /* transition energy in various units: */ 00138 Transitions[ipISO][nelem][ipHi][ipLo].EnergyWN = (realnum)EnergyWN; 00139 Transitions[ipISO][nelem][ipHi][ipLo].EnergyErg = (realnum)(EnergyWN*WAVNRYD*EN1RYD); 00140 Transitions[ipISO][nelem][ipHi][ipLo].EnergyK = (realnum)(EnergyWN*WAVNRYD*TE1RYD ); 00141 00142 ASSERT(Transitions[ipISO][nelem][ipHi][ipLo].EnergyWN >= 0.); 00143 ASSERT(Transitions[ipISO][nelem][ipHi][ipLo].EnergyErg >= 0.); 00144 ASSERT(Transitions[ipISO][nelem][ipHi][ipLo].EnergyK >= 0.); 00145 00147 if( N_(ipLo)==N_(ipHi) && ipISO==ipH_LIKE ) 00148 { 00149 Transitions[ipISO][nelem][ipHi][ipLo].WLAng = 0.; 00150 } 00151 else 00152 { 00153 /* make following an air wavelength */ 00154 Transitions[ipISO][nelem][ipHi][ipLo].WLAng = 00155 (realnum)(1.0e8/ 00156 Transitions[ipISO][nelem][ipHi][ipLo].EnergyWN/ 00157 RefIndex( Transitions[ipISO][nelem][ipHi][ipLo].EnergyWN)); 00158 ASSERT(Transitions[ipISO][nelem][ipHi][ipLo].WLAng > 0.); 00159 } 00160 00161 } 00162 } 00163 00164 /* fill the extra Lyman lines */ 00165 for( ipHi=2; ipHi < iso.nLyman[ipISO]; ipHi++ ) 00166 { 00167 FillExtraLymanLine( &ExtraLymanLines[ipISO][nelem][ipHi], ipISO, nelem, ipHi ); 00168 } 00169 } 00170 } 00171 } 00172 00173 /***************************************************************/ 00174 /***** Set up recombination tables for later interpolation *****/ 00175 /***************************************************************/ 00176 /* NB - the above is all we need if we are compiling recombination tables. */ 00177 iso_recomb_malloc(); 00178 iso_recomb_setup( ipH_LIKE ); 00179 iso_recomb_setup( ipHE_LIKE ); 00180 iso_recomb_auxiliary_free(); 00181 00182 /* set up helium collision tables */ 00183 HeCollidSetup(); 00184 00185 /***********************************************************************************/ 00186 /********** Transition Probabilities, Redistribution Functions, Opacitites ********/ 00187 /***********************************************************************************/ 00188 for( ipISO=ipH_LIKE; ipISO<NISO; ++ipISO ) 00189 { 00190 if( ipISO == ipH_LIKE ) 00191 { 00192 /* do nothing here */ 00193 } 00194 else if( ipISO == ipHE_LIKE ) 00195 { 00196 /* This routine reads in transition probabilities from a file. */ 00197 HelikeTransProbSetup(); 00198 } 00199 else 00200 { 00201 TotalInsanity(); 00202 } 00203 00204 for( nelem=ipISO; nelem < LIMELM; nelem++ ) 00205 { 00206 /* must always do helium even if turned off */ 00207 if( nelem < 2 || dense.lgElmtOn[nelem] ) 00208 { 00209 for( ipLo=ipH1s; ipLo < (iso.numLevels_max[ipISO][nelem] - 1); ipLo++ ) 00210 { 00211 for( ipHi=ipLo + 1; ipHi < iso.numLevels_max[ipISO][nelem]; ipHi++ ) 00212 { 00213 realnum Aul; 00214 00215 /* transition prob, EinstA uses current H atom indices */ 00216 if( ipISO == ipH_LIKE ) 00217 { 00218 Aul = hydro_transprob( nelem, ipHi, ipLo ); 00219 } 00220 else if( ipISO == ipHE_LIKE ) 00221 { 00222 Aul = helike_transprob(nelem, ipHi, ipLo); 00223 } 00224 else 00225 { 00226 TotalInsanity(); 00227 } 00228 00229 if( Aul <= iso.SmallA ) 00230 Transitions[ipISO][nelem][ipHi][ipLo].Emis = AddLine2Stack( false ); 00231 else 00232 Transitions[ipISO][nelem][ipHi][ipLo].Emis = AddLine2Stack( true ); 00233 00234 Transitions[ipISO][nelem][ipHi][ipLo].Emis->Aul = Aul; 00235 00236 ASSERT(Transitions[ipISO][nelem][ipHi][ipLo].Emis->Aul > 0.); 00237 00238 if( ipLo == 0 && ipHi == iso.nLyaLevel[ipISO] ) 00239 { 00240 /* this is La, special redistribution, default is ipLY_A */ 00241 Transitions[ipISO][nelem][ipHi][ipLo].Emis->iRedisFun = iso.ipLyaRedist[ipISO]; 00242 } 00243 else if( ipLo == 0 ) 00244 { 00245 /* these are rest of Lyman lines, 00246 * complete redistribution, doppler core only, K2 core, default ipCRD */ 00247 Transitions[ipISO][nelem][ipHi][ipLo].Emis->iRedisFun = iso.ipResoRedist[ipISO]; 00248 } 00249 else 00250 { 00251 /* all lines coming from excited states, default is complete 00252 * redis with wings, ipCRDW*/ 00253 Transitions[ipISO][nelem][ipHi][ipLo].Emis->iRedisFun = iso.ipSubRedist[ipISO]; 00254 } 00255 00256 if( Transitions[ipISO][nelem][ipHi][ipLo].Emis->Aul <= iso.SmallA || 00257 Transitions[ipISO][nelem][ipHi][ipLo].EnergyWN <= 0.) 00258 { 00259 Transitions[ipISO][nelem][ipHi][ipLo].Emis->gf = 0.; 00260 Transitions[ipISO][nelem][ipHi][ipLo].Emis->opacity = 0.; 00261 } 00262 else 00263 { 00264 Transitions[ipISO][nelem][ipHi][ipLo].Emis->gf = 00265 (realnum)(GetGF(Transitions[ipISO][nelem][ipHi][ipLo].Emis->Aul, 00266 Transitions[ipISO][nelem][ipHi][ipLo].EnergyWN, 00267 Transitions[ipISO][nelem][ipHi][ipLo].Hi->g)); 00268 ASSERT(Transitions[ipISO][nelem][ipHi][ipLo].Emis->gf > 0.); 00269 00270 /* derive the abs coef, call to function is gf, wl (A), g_low */ 00271 Transitions[ipISO][nelem][ipHi][ipLo].Emis->opacity = 00272 (realnum)(abscf(Transitions[ipISO][nelem][ipHi][ipLo].Emis->gf, 00273 Transitions[ipISO][nelem][ipHi][ipLo].EnergyWN, 00274 Transitions[ipISO][nelem][ipHi][ipLo].Lo->g)); 00275 ASSERT(Transitions[ipISO][nelem][ipHi][ipLo].Emis->opacity > 0.); 00276 } 00277 } 00278 } 00279 } 00280 } 00281 } 00282 00283 /************************************************/ 00284 /********** Fine Structure Mixing - FSM ********/ 00285 /************************************************/ 00286 if( iso.lgFSM[ipHE_LIKE] ) 00287 { 00288 /* set some special optical depth values */ 00289 for( nelem=ipISO; nelem < LIMELM; nelem++ ) 00290 { 00291 /* must always do helium even if turned off */ 00292 if( nelem < 2 || dense.lgElmtOn[nelem] ) 00293 { 00294 for( ipHi=ipHe2s3S; ipHi<iso.numLevels_max[ipHE_LIKE][nelem]; ipHi++ ) 00295 { 00296 for( ipLo=ipHe1s1S; ipLo<ipHi; ipLo++ ) 00297 { 00298 DoFSMixing( nelem, ipLo, ipHi ); 00299 } 00300 } 00301 } 00302 } 00303 } 00304 00305 /* following comes out very slightly off, correct here */ 00306 Transitions[ipH_LIKE][ipHELIUM][ipH3s][ipH2s].WLAng = 1640.f; 00307 Transitions[ipH_LIKE][ipHELIUM][ipH3s][ipH2p].WLAng = 1640.f; 00308 if( iso.n_HighestResolved_max[ipH_LIKE][ipHELIUM] >=3 ) 00309 { 00310 Transitions[ipH_LIKE][ipHELIUM][ipH3p][ipH2s].WLAng = 1640.f; 00311 Transitions[ipH_LIKE][ipHELIUM][ipH3p][ipH2p].WLAng = 1640.f; 00312 Transitions[ipH_LIKE][ipHELIUM][ipH3d][ipH2s].WLAng = 1640.f; 00313 Transitions[ipH_LIKE][ipHELIUM][ipH3d][ipH2p].WLAng = 1640.f; 00314 } 00315 00316 /****************************************************/ 00317 /********** lifetimes and damping constants ********/ 00318 /****************************************************/ 00319 for( ipISO=ipH_LIKE; ipISO<NISO; ipISO++ ) 00320 { 00321 for( nelem=ipISO; nelem < LIMELM; nelem++ ) 00322 { 00323 /* define these for H and He always */ 00324 if( nelem < 2 || dense.lgElmtOn[nelem] ) 00325 { 00326 /* these are not defined and must never be used */ 00327 StatesElem[ipISO][nelem][0].lifetime = -FLT_MAX; 00328 00329 for( ipHi=1; ipHi < iso.numLevels_max[ipISO][nelem]; ipHi++ ) 00330 { 00331 StatesElem[ipISO][nelem][ipHi].lifetime = SMALLFLOAT; 00332 00333 for( ipLo=0; ipLo < ipHi; ipLo++ ) 00334 { 00335 if( Transitions[ipISO][nelem][ipHi][ipLo].Emis->Aul <= iso.SmallA ) 00336 continue; 00337 00338 StatesElem[ipISO][nelem][ipHi].lifetime += Transitions[ipISO][nelem][ipHi][ipLo].Emis->Aul; 00339 } 00340 00341 /* sum of A's was just stuffed, now invert for lifetime. */ 00342 StatesElem[ipISO][nelem][ipHi].lifetime = 1./StatesElem[ipISO][nelem][ipHi].lifetime; 00343 00344 for( ipLo=0; ipLo < ipHi; ipLo++ ) 00345 { 00346 if( Transitions[ipISO][nelem][ipHi][ipLo].EnergyWN <= 0. ) 00347 continue; 00348 00349 if( Transitions[ipISO][nelem][ipHi][ipLo].Emis->Aul <= iso.SmallA ) 00350 continue; 00351 00352 Transitions[ipISO][nelem][ipHi][ipLo].Emis->dampXvel = (realnum)( 00353 (1.f/StatesElem[ipISO][nelem][ipHi].lifetime)/ 00354 PI4/Transitions[ipISO][nelem][ipHi][ipLo].EnergyWN); 00355 00356 ASSERT(Transitions[ipISO][nelem][ipHi][ipLo].Emis->dampXvel> 0.); 00357 } 00358 } 00359 } 00360 } 00361 } 00362 00363 /********************************************/ 00364 /********** Fix some 2-photon stuff ********/ 00365 /********************************************/ 00366 for( ipISO=ipH_LIKE; ipISO<NISO; ipISO++ ) 00367 { 00368 /* set some special optical depth values */ 00369 for( nelem=ipISO; nelem < LIMELM; nelem++ ) 00370 { 00371 /* must always do helium even if turned off */ 00372 if( nelem < 2 || dense.lgElmtOn[nelem] ) 00373 { 00374 /* Must think out two-photon treatment for other sequences. */ 00375 ASSERT( ipISO <= ipHE_LIKE ); 00376 00377 /* total optical depth in 1s - 2s */ 00378 Transitions[ipISO][nelem][1+ipISO][0].Emis->TauTot = 0.; 00379 00380 /* opacity in two-photon transition is incorrect - actually a continuum, 00381 * so divide by typical width*/ 00382 Transitions[ipISO][nelem][1+ipISO][0].Emis->opacity /= 1e4f; 00383 00384 /* wavelength of 0 is sentinel for two-photon emission */ 00385 Transitions[ipISO][nelem][1+ipISO][0].WLAng = 0.; 00386 } 00387 } 00388 } 00389 00390 00391 /* zero out some line information */ 00392 iso_zero(); 00393 00394 /* loop over iso sequences */ 00395 for( ipISO=ipH_LIKE; ipISO<NISO; ipISO++ ) 00396 { 00397 for( nelem = ipISO; nelem < LIMELM; nelem++ ) 00398 { 00399 /* must always do helium even if turned off */ 00400 if( nelem == ipISO || dense.lgElmtOn[nelem] ) 00401 { 00402 /* calculate cascade probabilities, branching ratios, and associated errors. */ 00403 iso_cascade( ipISO, nelem); 00404 } 00405 } 00406 } 00407 00408 iso_satellite(); 00409 00410 iso_satellite_update(); 00411 00412 /***************************************/ 00413 /********** Stark Broadening **********/ 00414 /***************************************/ 00415 /* fill in iso.strkar array */ 00416 for( ipISO=ipH_LIKE; ipISO<NISO; ++ipISO ) 00417 { 00418 for( nelem=ipISO; nelem < LIMELM; nelem++ ) 00419 { 00420 if( nelem < 2 || dense.lgElmtOn[nelem] ) 00421 { 00422 for( ipLo=ipH1s; ipLo < (iso.numLevels_max[ipISO][nelem] - 1); ipLo++ ) 00423 { 00424 for( ipHi= ipLo + 1; ipHi < iso.numLevels_max[ipISO][nelem]; ipHi++ ) 00425 { 00426 long nHi = StatesElem[ipISO][nelem][ipHi].n; 00427 long nLo = StatesElem[ipISO][nelem][ipLo].n; 00428 00429 iso.strkar[ipISO][nelem][ipLo][ipHi] = (realnum)pow((realnum)( nLo * nHi ),(realnum)1.2f); 00430 iso.pestrk[ipISO][nelem][ipLo][ipHi] = 0.; 00431 } 00432 } 00433 } 00434 } 00435 } 00436 00437 return; 00438 } 00439 00440 /* ============================================================================== */ 00441 STATIC void iso_zero(void) 00442 { 00443 long int i, 00444 ipHi, 00445 ipISO, 00446 ipLo, 00447 nelem; 00448 00449 DEBUG_ENTRY( "iso_zero()" ); 00450 00451 fixit(); /* this routine appears to be completely unnecessary because all of 00452 * this is done in RT_tau_init shortly later. */ 00453 00454 hydro.HLineWidth = 0.; 00455 00456 /****************************************************/ 00457 /********** initialize some variables **********/ 00458 /****************************************************/ 00459 for( ipISO=ipH_LIKE; ipISO<NISO; ++ipISO ) 00460 { 00461 for( nelem=ipISO; nelem < LIMELM; nelem++ ) 00462 { 00463 if( nelem < 2 || dense.lgElmtOn[nelem] ) 00464 { 00465 for( ipHi=0; ipHi < iso.numLevels_max[ipISO][nelem]; ipHi++ ) 00466 { 00467 StatesElem[ipISO][nelem][ipHi].Pop = 0.; 00468 00469 iso.PopLTE[ipISO][nelem][ipHi] = 0.; 00470 iso.ColIoniz[ipISO][nelem][ipHi] = 0.; 00471 iso.gamnc[ipISO][nelem][ipHi] = -DBL_MAX; 00472 iso.RecomInducRate[ipISO][nelem][ipHi] = -DBL_MAX; 00473 iso.DepartCoef[ipISO][nelem][ipHi] = -DBL_MAX; 00474 iso.RateLevel2Cont[ipISO][nelem][ipHi] = 0.; 00475 iso.RateCont2Level[ipISO][nelem][ipHi] = 0.; 00476 iso.ConOpacRatio[ipISO][nelem][ipHi] = 1.f; 00477 iso.RadRecomb[ipISO][nelem][ipHi][ipRecRad] = 0.; 00478 iso.RadRecomb[ipISO][nelem][ipHi][ipRecNetEsc] = 1.; 00479 iso.RadRecomb[ipISO][nelem][ipHi][ipRecEsc] = 1.; 00480 iso.DielecRecomb[ipISO][nelem][ipHi] = 0.; 00481 } 00482 } 00483 } 00484 } 00485 00486 /* ground state of H and He is different since totally determine 00487 * their own opacities */ 00488 iso.ConOpacRatio[ipH_LIKE][ipHYDROGEN][0] = 1e-5f; 00489 iso.ConOpacRatio[ipH_LIKE][ipHELIUM][0] = 1e-5f; 00490 iso.ConOpacRatio[ipHE_LIKE][ipHELIUM][0] = 1e-5f; 00491 00492 00493 /* main isoelectronic arrays, fill with sane values */ 00494 for( ipISO=ipH_LIKE; ipISO<NISO; ++ipISO ) 00495 { 00496 for( nelem=ipISO; nelem < LIMELM; nelem++ ) 00497 { 00498 /* >>chng 05 dec 30, move nelem.numLevel to numLevels_local and int here */ 00499 iso.numLevels_local[ipISO][nelem] = iso.numLevels_max[ipISO][nelem]; 00500 00501 /* must always do helium even if turned off */ 00502 if( nelem < 2 || dense.lgElmtOn[nelem] ) 00503 { 00504 for( ipHi=1; ipHi < iso.numLevels_max[ipISO][nelem]; ipHi++ ) 00505 { 00506 TransitionZero( &SatelliteLines[ipISO][nelem][ipHi] ); 00507 00508 for( ipLo=0; ipLo < ipHi; ipLo++ ) 00509 { 00510 TransitionZero( &Transitions[ipISO][nelem][ipHi][ipLo] ); 00511 } 00512 } 00513 00514 for( ipHi=2; ipHi < iso.nLyman[ipISO]; ipHi++ ) 00515 { 00516 TransitionZero( &ExtraLymanLines[ipISO][nelem][ipHi] ); 00517 } 00518 } 00519 } 00520 } 00521 00522 /* initialize heavy element line array */ 00523 for( i=0; i <= nLevel1; i++ ) 00524 { 00525 TransitionZero( &TauLines[i] ); 00526 } 00527 00528 /* this is a dummy optical depth array for non-existant lines 00529 * when this goes over to struc, make sure all are set to zero here since 00530 * init in grid may depend on it */ 00531 TransitionZero( &TauDummy ); 00532 00533 for( i=0; i < nUTA; i++ ) 00534 { 00535 /* heat is special for this array - it is heat per pump */ 00536 double hsave = UTALines[i].Coll.heat; 00537 TransitionZero( &UTALines[i] ); 00538 UTALines[i].Coll.heat = hsave; 00539 } 00540 00541 for( i=0; i < nWindLine; i++ ) 00542 { 00543 TransitionZero( &TauLine2[i] ); 00544 } 00545 00546 for( i=0; i < nHFLines; i++ ) 00547 { 00548 TransitionZero( &HFLines[i] ); 00549 } 00550 00551 /* database lines - DB lines */ 00552 for( i=0; i <linesAdded2; i++) 00553 { 00554 TransitionZero( atmolEmis[i].tran ); 00555 } 00556 00557 for( i=0; i < nCORotate; i++ ) 00558 { 00559 TransitionZero( &C12O16Rotate[i] ); 00560 TransitionZero( &C13O16Rotate[i] ); 00561 } 00562 return; 00563 } 00564 00565 STATIC void iso_allocate(void) 00566 { 00567 long int 00568 ipISO, 00569 nelem, 00570 n, 00571 i, 00572 i1, 00573 j, 00574 ipHi, 00575 ipLo; 00576 00577 DEBUG_ENTRY( "iso_allocate()" ); 00578 00579 iso.strkar.reserve( NISO ); 00580 iso.ipOpac.reserve( NISO ); 00581 iso.RadRecomb.reserve( NISO ); 00582 iso.DielecRecombVsTemp.reserve( NISO ); 00583 iso.Boltzmann.reserve( NISO ); 00584 iso.QuantumNumbers2Index.reserve( NISO ); 00585 iso.BranchRatio.reserve( NISO ); 00586 iso.CascadeProb.reserve( NISO ); 00587 iso.CachedAs.reserve( NISO ); 00588 00589 /* the hydrogen and helium like iso-sequences */ 00590 for( ipISO=ipH_LIKE; ipISO<NISO; ++ipISO ) 00591 { 00592 iso.strkar.reserve( ipISO, LIMELM ); 00593 iso.ipOpac.reserve( ipISO, LIMELM ); 00594 iso.RadRecomb.reserve( ipISO, LIMELM ); 00595 iso.DielecRecombVsTemp.reserve( ipISO, LIMELM ); 00596 iso.Boltzmann.reserve( ipISO, LIMELM ); 00597 iso.QuantumNumbers2Index.reserve( ipISO, LIMELM ); 00598 iso.BranchRatio.reserve( ipISO, LIMELM ); 00599 iso.CascadeProb.reserve( ipISO, LIMELM ); 00600 iso.CachedAs.reserve( ipISO, LIMELM ); 00601 00602 for( nelem=ipISO; nelem < LIMELM; ++nelem ) 00603 { 00604 /* only grab core for elements that are turned on */ 00605 if( nelem < 2 || dense.lgElmtOn[nelem] ) 00606 { 00607 iso.numLevels_malloc[ipISO][nelem] = iso.numLevels_max[ipISO][nelem]; 00608 00609 /*fprintf(ioQQQ,"assert failll %li\t%li\t%li\n", ipISO, nelem , iso.numLevels_max[ipISO][nelem] );*/ 00610 ASSERT( iso.numLevels_max[ipISO][nelem] > 0 ); 00611 00612 iso.nLyman_malloc[ipISO] = iso.nLyman[ipISO]; 00613 00614 iso.strkar.reserve( ipISO, nelem, iso.numLevels_max[ipISO][nelem] ); 00615 iso.ipOpac.reserve( ipISO, nelem, iso.numLevels_max[ipISO][nelem] ); 00616 iso.RadRecomb.reserve( ipISO, nelem, iso.numLevels_max[ipISO][nelem] ); 00617 iso.DielecRecombVsTemp.reserve( ipISO, nelem, iso.numLevels_max[ipISO][nelem] ); 00618 iso.Boltzmann.reserve( ipISO, nelem, iso.numLevels_max[ipISO][nelem] ); 00619 00620 iso.CachedAs.reserve( ipISO, nelem, MAX2(1, iso.nCollapsed_max[ipISO][nelem]) ); 00621 00622 for( i = 0; i < iso.nCollapsed_max[ipISO][nelem]; ++i ) 00623 { 00624 iso.CachedAs.reserve( ipISO, nelem, i, iso.numLevels_max[ipISO][nelem] - iso.nCollapsed_max[ipISO][nelem] ); 00625 for( i1 = 0; i1 < iso.numLevels_max[ipISO][nelem] - iso.nCollapsed_max[ipISO][nelem]; ++i1 ) 00626 { 00627 /* allocate two spaces delta L +/- 1 */ 00628 iso.CachedAs.reserve( ipISO, nelem, i, i1, 2 ); 00629 } 00630 } 00631 00632 iso.QuantumNumbers2Index.reserve( ipISO, nelem, iso.n_HighestResolved_max[ipISO][nelem] + iso.nCollapsed_max[ipISO][nelem] + 1 ); 00633 00634 00635 for( i = 1; i <= (iso.n_HighestResolved_max[ipISO][nelem] + iso.nCollapsed_max[ipISO][nelem]); ++i ) 00636 { 00637 /* Allocate proper number of angular momentum quantum number. */ 00638 iso.QuantumNumbers2Index.reserve( ipISO, nelem, i, i ); 00639 00640 for( i1 = 0; i1 < i; ++i1 ) 00641 { 00642 /* This may have to change for other iso sequences. */ 00643 ASSERT( NISO == 2 ); 00644 /* Allocate 4 spaces for multiplicity. H-like will be accessed with "2" for doublet, 00645 * He-like will be accessed via "1" for singlet or "3" for triplet. "0" will not be used. */ 00646 iso.QuantumNumbers2Index.reserve( ipISO, nelem, i, i1, 4 ); 00647 } 00648 } 00649 00650 for( n=0; n < iso.numLevels_max[ipISO][nelem]; ++n ) 00651 { 00652 iso.RadRecomb.reserve( ipISO, nelem, n, 3 ); 00653 00654 /* sec to last dim is upper level n, 00655 * last dim of this array is lower level, will go from 0 to n-1 */ 00656 if( n > 0 ) 00657 iso.Boltzmann.reserve( ipISO, nelem, n, n ); 00658 00659 /* malloc space for number of temperature points in Badnell grids */ 00660 iso.DielecRecombVsTemp.reserve( ipISO, nelem, n, NUM_DR_TEMPS ); 00661 } 00662 00663 /* In this array are stored the C values described in Robbins 68. */ 00664 iso.CascadeProb.reserve( ipISO, nelem, iso.numLevels_max[ipISO][nelem] ); 00665 iso.BranchRatio.reserve( ipISO, nelem, iso.numLevels_max[ipISO][nelem]+1 ); 00666 00667 for( i = 1; i < iso.numLevels_max[ipISO][nelem]; i++ ) 00668 iso.BranchRatio.reserve( ipISO, nelem, i, i+1 ); 00669 00670 /* reserve final dimension of cascade probabilities. */ 00671 for( i = 0; i < iso.numLevels_max[ipISO][nelem]; ++i ) 00672 { 00673 iso.strkar.reserve( ipISO, nelem, i, iso.numLevels_max[ipISO][nelem] ); 00674 iso.CascadeProb.reserve( ipISO, nelem, i, i+1 ); 00675 } 00676 } 00677 } 00678 } 00679 00680 iso.strkar.alloc(); 00681 iso.pestrk.alloc( iso.strkar.clone() ); 00682 iso.ipOpac.alloc(); 00683 secondaries.Hx12.alloc( iso.ipOpac.clone() ); 00684 iso.ipIsoLevNIonCon.alloc( iso.ipOpac.clone() ); 00685 iso.xIsoLevNIonRyd.alloc( iso.ipOpac.clone() ); 00686 iso.DielecRecomb.alloc( iso.ipOpac.clone() ); 00687 iso.DepartCoef.alloc( iso.ipOpac.clone() ); 00688 iso.RateLevel2Cont.alloc( iso.ipOpac.clone() ); 00689 iso.RateCont2Level.alloc( iso.ipOpac.clone() ); 00690 iso.ConOpacRatio.alloc( iso.ipOpac.clone() ); 00691 iso.gamnc.alloc( iso.ipOpac.clone() ); 00692 iso.RecomInducRate.alloc( iso.ipOpac.clone() ); 00693 iso.RecomInducCool_Coef.alloc( iso.ipOpac.clone() ); 00694 iso.PhotoHeat.alloc( iso.ipOpac.clone() ); 00695 iso.PopLTE.alloc( iso.ipOpac.clone() ); 00696 iso.ColIoniz.alloc( iso.ipOpac.clone() ); 00697 iso.RadEffec.alloc( iso.ipOpac.clone() ); 00698 iso.RadRecomb.alloc(); 00699 iso.Boltzmann.alloc(); 00700 iso.DielecRecombVsTemp.alloc(); 00701 iso.CachedAs.alloc(); 00702 iso.QuantumNumbers2Index.alloc(); 00703 iso.BranchRatio.alloc(); 00704 iso.CascadeProb.alloc(); 00705 00706 iso.bnl_effective.alloc( iso.QuantumNumbers2Index.clone() ); 00707 00708 iso.DielecRecombVsTemp.zero(); 00709 iso.QuantumNumbers2Index.invalidate(); 00710 iso.bnl_effective.invalidate(); 00711 iso.BranchRatio.invalidate(); 00712 iso.CascadeProb.invalidate(); 00713 00714 StatesElem.reserve( NISO ); 00715 Transitions.reserve( NISO ); 00716 ExtraLymanLines.reserve( NISO ); 00717 00718 for( ipISO=ipH_LIKE; ipISO<NISO; ++ipISO ) 00719 { 00720 StatesElem.reserve( ipISO, LIMELM ); 00721 Transitions.reserve( ipISO, LIMELM ); 00722 ExtraLymanLines.reserve( ipISO, LIMELM ); 00723 00724 for( nelem=ipISO; nelem < LIMELM; ++nelem ) 00725 { 00726 /* only grab core for elements that are turned on */ 00727 if( nelem < 2 || dense.lgElmtOn[nelem] ) 00728 { 00729 ASSERT( iso.numLevels_max[ipISO][nelem] > 0 ); 00730 00731 StatesElem.reserve( ipISO, nelem, iso.numLevels_max[ipISO][nelem] ); 00732 Transitions.reserve( ipISO, nelem, iso.numLevels_max[ipISO][nelem] ); 00733 ExtraLymanLines.reserve( ipISO, nelem, iso.nLyman[ipISO] ); 00734 00735 for( ipHi=1; ipHi<iso.numLevels_max[ipISO][nelem]; ipHi++ ) 00736 Transitions.reserve( ipISO, nelem, ipHi, ipHi ); 00737 } 00738 } 00739 } 00740 00741 StatesElem.alloc(); 00742 SatelliteLines.alloc( StatesElem.clone() ); 00743 Transitions.alloc(); 00744 ExtraLymanLines.alloc(); 00745 00746 for( ipISO=ipH_LIKE; ipISO<NISO; ++ipISO ) 00747 { 00748 for( nelem=ipISO; nelem < LIMELM; ++nelem ) 00749 { 00750 /* only grab core for elements that are turned on */ 00751 if( nelem < 2 || dense.lgElmtOn[nelem] ) 00752 { 00753 for( ipLo=0; ipLo<iso.numLevels_max[ipISO][nelem]; ipLo++ ) 00754 { 00755 /* Upper level is continuum, use a generic state 00756 * lower level is the same as the index. */ 00757 TransitionJunk( &SatelliteLines[ipISO][nelem][ipLo] ); 00758 SatelliteLines[ipISO][nelem][ipLo].Hi = AddState2Stack(); 00759 SatelliteLines[ipISO][nelem][ipLo].Lo = &StatesElem[ipISO][nelem][ipLo]; 00760 SatelliteLines[ipISO][nelem][ipLo].Emis = AddLine2Stack( true ); 00761 } 00762 00763 for( ipHi=1; ipHi<iso.numLevels_max[ipISO][nelem]; ipHi++ ) 00764 { 00765 for( ipLo=0; ipLo < ipHi; ipLo++ ) 00766 { 00767 /* set ENTIRE array to impossible values, in case of bad pointer */ 00768 TransitionJunk( &Transitions[ipISO][nelem][ipHi][ipLo] ); 00769 Transitions[ipISO][nelem][ipHi][ipLo].Hi = &StatesElem[ipISO][nelem][ipHi]; 00770 Transitions[ipISO][nelem][ipHi][ipLo].Lo = &StatesElem[ipISO][nelem][ipLo]; 00771 } 00772 } 00773 00774 /* junk the extra Lyman lines */ 00775 for( ipHi=2; ipHi < iso.nLyman[ipISO]; ipHi++ ) 00776 { 00777 TransitionJunk( &ExtraLymanLines[ipISO][nelem][ipHi] ); 00778 ExtraLymanLines[ipISO][nelem][ipHi].Hi = AddState2Stack(); 00779 /* lower level is just ground state of the ion */ 00780 ExtraLymanLines[ipISO][nelem][ipHi].Lo = &StatesElem[ipISO][nelem][0]; 00781 ExtraLymanLines[ipISO][nelem][ipHi].Emis = AddLine2Stack( true ); 00782 } 00783 } 00784 } 00785 } 00786 00787 /* malloc space for random error generation, if appropriate flags set. */ 00788 for( ipISO=ipH_LIKE; ipISO<NISO; ++ipISO ) 00789 { 00790 static bool lgNotSetYet = true; 00791 00792 if( iso.lgRandErrGen[ipISO] && lgNotSetYet ) 00793 { 00794 iso.Error.reserve( NISO ); 00795 iso.SigmaCascadeProb.reserve( NISO ); 00796 lgNotSetYet = false; 00797 } 00798 00799 if( iso.lgRandErrGen[ipISO] ) 00800 { 00801 ASSERT( !lgNotSetYet ); 00802 iso.Error.reserve( ipISO, LIMELM ); 00803 iso.SigmaCascadeProb.reserve( ipISO, LIMELM ); 00804 00805 for( nelem=ipISO; nelem < LIMELM; ++nelem ) 00806 { 00807 if( nelem == ipHELIUM || dense.lgElmtOn[nelem] ) 00808 { 00809 /* Here we add one slot for rates involving the continuum. */ 00810 iso.Error.reserve( ipISO, nelem, iso.numLevels_max[ipISO][nelem]+1 ); 00811 00812 for( i = 1; i< iso.numLevels_max[ipISO][nelem] + 1; i++ ) 00813 { 00814 iso.Error.reserve( ipISO, nelem, i, i+1 ); 00815 00816 for( j = 0; j<i+1; j++ ) 00817 iso.Error.reserve( ipISO, nelem, i, j, 3 ); 00818 } 00819 00820 /* Uncertainty in cascade probability and effective recombination, 00821 * only when error generation is turned on. */ 00822 iso.SigmaCascadeProb.reserve( ipISO, nelem, iso.numLevels_max[ipISO][nelem] ); 00823 for( i=0; i < iso.numLevels_max[ipISO][nelem]; i++ ) 00824 /* error in cascade probabilities, for all levels */ 00825 iso.SigmaCascadeProb.reserve( ipISO, nelem, i, i+1 ); 00826 } 00827 } 00828 } 00829 } 00830 00831 iso.Error.alloc(); 00832 iso.ErrorFactor.alloc( iso.Error.clone() ); 00833 iso.SigmaAtot.alloc( iso.ipOpac.clone() ); 00834 iso.SigmaRadEffec.alloc( iso.RadEffec.clone() ); 00835 iso.SigmaCascadeProb.alloc(); 00836 00837 iso.Error.invalidate(); 00838 iso.ErrorFactor.invalidate(); 00839 00840 return; 00841 } 00842 00843 STATIC void iso_assign_quantum_numbers(void) 00844 { 00845 long int 00846 ipISO, 00847 ipLo, 00848 level, 00849 nelem, 00850 i, 00851 in, 00852 il, 00853 is, 00854 ij; 00855 00856 DEBUG_ENTRY( "iso_assign_quantum_numbers()" ); 00857 00858 ipISO = ipH_LIKE; 00859 for( nelem=ipISO; nelem < LIMELM; nelem++ ) 00860 { 00861 /* only check elements that are turned on */ 00862 if( nelem == ipHELIUM || dense.lgElmtOn[nelem] ) 00863 { 00864 i = 0; 00865 00866 /* 2 for doublet */ 00867 is = 2; 00868 00869 /* fill in StatesElem and iso.QuantumNumbers2Index arrays. */ 00870 /* this loop is over quantum number n */ 00871 for( in = 1L; in <= iso.n_HighestResolved_max[ipISO][nelem]; ++in ) 00872 { 00873 for( il = 0L; il < in; ++il ) 00874 { 00875 StatesElem[ipISO][nelem][i].n = in; 00876 StatesElem[ipISO][nelem][i].S = is; 00877 StatesElem[ipISO][nelem][i].l = il; 00878 StatesElem[ipISO][nelem][i].j = -1; 00879 iso.QuantumNumbers2Index[ipISO][nelem][in][il][is] = i; 00880 ++i; 00881 } 00882 } 00883 /* now do the collapsed levels */ 00884 in = iso.n_HighestResolved_max[ipISO][nelem] + 1; 00885 for( level = i; level< iso.numLevels_max[ipISO][nelem]; ++level) 00886 { 00887 StatesElem[ipISO][nelem][level].n = in; 00888 StatesElem[ipISO][nelem][level].S = -LONG_MAX; 00889 StatesElem[ipISO][nelem][level].l = -LONG_MAX; 00890 StatesElem[ipISO][nelem][level].j = -1; 00891 /* Point every l to same index for collapsed levels. */ 00892 for( il = 0; il < in; ++il ) 00893 { 00894 iso.QuantumNumbers2Index[ipISO][nelem][in][il][is] = level; 00895 } 00896 ++in; 00897 } 00898 --in; 00899 00900 /* confirm that we did not overrun the array */ 00901 ASSERT( i <= iso.numLevels_max[ipISO][nelem] ); 00902 00903 /* confirm that n is positive and not greater than the max n. */ 00904 ASSERT( (in > 0) && (in < (iso.n_HighestResolved_max[ipISO][nelem] + iso.nCollapsed_max[ipISO][nelem] + 1) ) ); 00905 00906 /* Verify StatesElem[ipISO] and iso.QuantumNumbers2Index[ipISO] agree in all cases */ 00907 for( in = 2L; in <= iso.n_HighestResolved_max[ipISO][nelem] + iso.nCollapsed_max[ipISO][nelem]; ++in ) 00908 { 00909 for( il = 0L; il < in; ++il ) 00910 { 00911 ASSERT( StatesElem[ipISO][nelem][ iso.QuantumNumbers2Index[ipISO][nelem][in][il][is] ].n == in ); 00912 if( in <= iso.n_HighestResolved_max[ipISO][nelem] ) 00913 { 00914 /* Must only check these for resolved levels... 00915 * collapsed levels in StatesElem[ipISO] have pointers for l and s that will blow if used. */ 00916 ASSERT( StatesElem[ipISO][nelem][ iso.QuantumNumbers2Index[ipISO][nelem][in][il][is] ].l == il ); 00917 ASSERT( StatesElem[ipISO][nelem][ iso.QuantumNumbers2Index[ipISO][nelem][in][il][is] ].S == is ); 00918 } 00919 } 00920 } 00921 } 00922 } 00923 00924 /* then do he-like */ 00925 ipISO = ipHE_LIKE; 00926 for( nelem=ipHELIUM; nelem < LIMELM; nelem++ ) 00927 { 00928 /* only check elements that are turned on */ 00929 if( nelem == ipHELIUM || dense.lgElmtOn[nelem] ) 00930 { 00931 i = 0; 00932 00933 /* fill in StatesElem and iso.QuantumNumbers2Index arrays. */ 00934 /* this loop is over quantum number n */ 00935 for( in = 1L; in <= iso.n_HighestResolved_max[ipISO][nelem]; ++in ) 00936 { 00937 for( il = 0L; il < in; ++il ) 00938 { 00939 for( is = 3L; is >= 1L; is -= 2 ) 00940 { 00941 /* All levels except singlet P follow the ordering scheme: */ 00942 /* lower l's have lower energy */ 00943 /* triplets have lower energy */ 00944 if( (il == 1L) && (is == 1L) ) 00945 continue; 00946 /* n = 1 has no triplet, of course. */ 00947 if( (in == 1L) && (is == 3L) ) 00948 continue; 00949 00950 /* singlets */ 00951 if( is == 1 ) 00952 { 00953 StatesElem[ipISO][nelem][i].n = in; 00954 StatesElem[ipISO][nelem][i].S = is; 00955 StatesElem[ipISO][nelem][i].l = il; 00956 /* this is not a typo, J=L for singlets. */ 00957 StatesElem[ipISO][nelem][i].j = il; 00958 iso.QuantumNumbers2Index[ipISO][nelem][in][il][is] = i; 00959 ++i; 00960 } 00961 /* 2 triplet P is j-resolved */ 00962 else if( (in == 2) && (il == 1) && (is == 3) ) 00963 { 00964 ij = 0; 00965 do 00966 { 00967 StatesElem[ipISO][nelem][i].n = in; 00968 StatesElem[ipISO][nelem][i].S = is; 00969 StatesElem[ipISO][nelem][i].l = il; 00970 StatesElem[ipISO][nelem][i].j = ij; 00971 iso.QuantumNumbers2Index[ipISO][nelem][in][il][is] = i; 00972 ++i; 00973 ++ij; 00974 /* repeat this for the separate j-levels within 2^3P. */ 00975 } while ( ij < 3 ); 00976 } 00977 else 00978 { 00979 StatesElem[ipISO][nelem][i].n = in; 00980 StatesElem[ipISO][nelem][i].S = is; 00981 StatesElem[ipISO][nelem][i].l = il; 00982 StatesElem[ipISO][nelem][i].j = -1L; 00983 iso.QuantumNumbers2Index[ipISO][nelem][in][il][is] = i; 00984 ++i; 00985 } 00986 } 00987 } 00988 /* Insert singlet P at the end of every sequence for a given n. */ 00989 if( in > 1L ) 00990 { 00991 StatesElem[ipISO][nelem][i].n = in; 00992 StatesElem[ipISO][nelem][i].S = 1L; 00993 StatesElem[ipISO][nelem][i].l = 1L; 00994 StatesElem[ipISO][nelem][i].j = 1L; 00995 iso.QuantumNumbers2Index[ipISO][nelem][in][1][1] = i; 00996 ++i; 00997 } 00998 } 00999 /* now do the collapsed levels */ 01000 in = iso.n_HighestResolved_max[ipISO][nelem] + 1; 01001 for( level = i; level< iso.numLevels_max[ipISO][nelem]; ++level) 01002 { 01003 StatesElem[ipISO][nelem][level].n = in; 01004 StatesElem[ipISO][nelem][level].S = -LONG_MAX; 01005 StatesElem[ipISO][nelem][level].l = -LONG_MAX; 01006 StatesElem[ipISO][nelem][level].j = -1; 01007 /* Point every l and s to same index for collapsed levels. */ 01008 for( il = 0; il < in; ++il ) 01009 { 01010 for( is = 1; is <= 3; is += 2 ) 01011 { 01012 iso.QuantumNumbers2Index[ipISO][nelem][in][il][is] = level; 01013 } 01014 } 01015 ++in; 01016 } 01017 --in; 01018 01019 /* confirm that we did not overrun the array */ 01020 ASSERT( i <= iso.numLevels_max[ipISO][nelem] ); 01021 01022 /* confirm that n is positive and not greater than the max n. */ 01023 ASSERT( (in > 0) && (in < (iso.n_HighestResolved_max[ipISO][nelem] + iso.nCollapsed_max[ipISO][nelem] + 1) ) ); 01024 01025 /* Verify StatesElem[ipISO] and iso.QuantumNumbers2Index[ipISO] agree in all cases */ 01026 for( in = 2L; in <= iso.n_HighestResolved_max[ipISO][nelem] + iso.nCollapsed_max[ipISO][nelem]; ++in ) 01027 { 01028 for( il = 0L; il < in; ++il ) 01029 { 01030 for( is = 3L; is >= 1; is -= 2 ) 01031 { 01032 /* Ground state is not triplicate. */ 01033 if( (in == 1L) && (is == 3L) ) 01034 continue; 01035 01036 ASSERT( StatesElem[ipISO][nelem][ iso.QuantumNumbers2Index[ipISO][nelem][in][il][is] ].n == in ); 01037 if( in <= iso.n_HighestResolved_max[ipISO][nelem] ) 01038 { 01039 /* Must only check these for resolved levels... 01040 * collapsed levels in StatesElem[ipISO] have pointers for l and s that will blow if used. */ 01041 ASSERT( StatesElem[ipISO][nelem][ iso.QuantumNumbers2Index[ipISO][nelem][in][il][is] ].l == il ); 01042 ASSERT( StatesElem[ipISO][nelem][ iso.QuantumNumbers2Index[ipISO][nelem][in][il][is] ].S == is ); 01043 } 01044 } 01045 } 01046 } 01047 } 01048 } 01049 01050 for( ipISO=ipH_LIKE; ipISO<NISO; ipISO++ ) 01051 { 01052 for( nelem=ipISO; nelem < LIMELM; nelem++ ) 01053 { 01054 /* must always do helium even if turned off */ 01055 if( nelem < 2 || dense.lgElmtOn[nelem] ) 01056 { 01057 for( ipLo=ipH1s; ipLo < iso.numLevels_max[ipISO][nelem]; ipLo++ ) 01058 { 01059 StatesElem[ipISO][nelem][ipLo].nelem = (int)(nelem+1); 01060 StatesElem[ipISO][nelem][ipLo].IonStg = (int)(nelem+1-ipISO); 01061 01062 if( StatesElem[ipISO][nelem][ipLo].j >= 0 ) 01063 { 01064 StatesElem[ipISO][nelem][ipLo].g = 2.f*StatesElem[ipISO][nelem][ipLo].j+1.f; 01065 } 01066 else if( StatesElem[ipISO][nelem][ipLo].l >= 0 ) 01067 { 01068 StatesElem[ipISO][nelem][ipLo].g = (2.f*StatesElem[ipISO][nelem][ipLo].l+1.f) * 01069 StatesElem[ipISO][nelem][ipLo].S; 01070 } 01071 else 01072 { 01073 if( ipISO == ipH_LIKE ) 01074 StatesElem[ipISO][nelem][ipLo].g = 2.f*(realnum)POW2( StatesElem[ipISO][nelem][ipLo].n ); 01075 else if( ipISO == ipHE_LIKE ) 01076 StatesElem[ipISO][nelem][ipLo].g = 4.f*(realnum)POW2( StatesElem[ipISO][nelem][ipLo].n ); 01077 else 01078 { 01079 /* replace this with correct thing if more sequences are added. */ 01080 TotalInsanity(); 01081 } 01082 } 01083 char chConfiguration[11] = " "; 01084 long nCharactersWritten = 0; 01085 01086 /* include j only if defined. */ 01087 if( StatesElem[ipISO][nelem][ipLo].n > iso.n_HighestResolved_max[ipISO][nelem] ) 01088 { 01089 nCharactersWritten = sprintf( chConfiguration, "n=%3li", 01090 StatesElem[ipISO][nelem][ipLo].n ); 01091 } 01092 else if( StatesElem[ipISO][nelem][ipLo].j > 0 ) 01093 { 01094 nCharactersWritten = sprintf( chConfiguration, "%3li^%2li%c_%li", 01095 StatesElem[ipISO][nelem][ipLo].n, 01096 StatesElem[ipISO][nelem][ipLo].S, 01097 chL[ MIN2( 20, StatesElem[ipISO][nelem][ipLo].l ) ], 01098 StatesElem[ipISO][nelem][ipLo].j ); 01099 } 01100 else 01101 { 01102 nCharactersWritten = sprintf( chConfiguration, "%3li^%2li%c", 01103 StatesElem[ipISO][nelem][ipLo].n, 01104 StatesElem[ipISO][nelem][ipLo].S, 01105 chL[ MIN2( 20, StatesElem[ipISO][nelem][ipLo].l) ] ); 01106 } 01107 01108 ASSERT( nCharactersWritten <= 10 ); 01109 chConfiguration[10] = '\0'; 01110 01111 strncpy( StatesElem[ipISO][nelem][ipLo].chConfig, chConfiguration, 10 ); 01112 } 01113 } 01114 } 01115 } 01116 return; 01117 } 01118 01119 #if defined(__ICC) && defined(__i386) 01120 #pragma optimization_level 1 01121 #endif 01122 STATIC void FillExtraLymanLine( transition *t, long ipISO, long nelem, long nHi ) 01123 { 01124 double Enerwn, Aul; 01125 01126 DEBUG_ENTRY( "FillExtraLymanLine()" ); 01127 01128 /* atomic number or charge and stage: */ 01129 t->Hi->nelem = (int)(nelem+1); 01130 t->Hi->IonStg = (int)(nelem+1-ipISO); 01131 01132 /* statistical weight is same as statistical weight of corresponding LyA. */ 01133 t->Hi->g = StatesElem[ipISO][nelem][iso.nLyaLevel[ipISO]].g; 01134 01135 /* energies */ 01136 Enerwn = iso.xIsoLevNIonRyd[ipISO][nelem][0] * RYD_INF * ( 1. - 1./POW2((double)nHi) ); 01137 01138 /* transition energy in various units:*/ 01139 t->EnergyWN = (realnum)(Enerwn); 01140 t->EnergyErg = (realnum)(Enerwn * WAVNRYD * EN1RYD); 01141 t->EnergyK = (realnum)(Enerwn * WAVNRYD * TE1RYD); 01142 t->WLAng = (realnum)(1.0e8/ Enerwn/ RefIndex(Enerwn)); 01143 01144 if( ipISO == ipH_LIKE ) 01145 { 01146 Aul = H_Einstein_A( nHi, 1, 1, 0, nelem+1 ); 01147 } 01148 else 01149 { 01150 if( nelem == ipHELIUM ) 01151 { 01152 /* A simple fit for the calculation of Helium lyman Aul's. */ 01153 Aul = (1.508e10) / pow((double)nHi,2.975); 01154 } 01155 else 01156 { 01157 /* Fit to values given in 01158 * >>refer He-like As Johnson, W.R., Savukov, I.M., Safronova, U.I., & 01159 * >>refercon Dalgarno, A., 2002, ApJS 141, 543J */ 01160 /* originally astro.ph. 0201454 */ 01161 Aul = 1.375E10 * pow((double)nelem, 3.9) / pow((double)nHi,3.1); 01162 } 01163 } 01164 01165 t->Emis->Aul = (realnum)Aul; 01166 01167 t->Hi->lifetime = iso_state_lifetime( ipISO, nelem, nHi, 1 ); 01168 01169 t->Emis->dampXvel = (realnum)( 1.f / t->Hi->lifetime / PI4 / t->EnergyWN ); 01170 01171 t->Emis->iRedisFun = iso.ipResoRedist[ipISO]; 01172 01173 t->Emis->gf = (realnum)(GetGF(t->Emis->Aul, t->EnergyWN, t->Hi->g)); 01174 01175 /* derive the abs coef, call to function is Emis.gf, wl (A), g_low */ 01176 t->Emis->opacity = (realnum)(abscf(t->Emis->gf, t->EnergyWN, t->Lo->g)); 01177 01178 /* create array indices that will blow up */ 01179 t->ipCont = INT_MIN; 01180 t->Emis->ipFine = INT_MIN; 01181 01182 { 01183 /* option to print particulars of some line when called 01184 * a prettier print statement is near where chSpin is defined below 01185 * search for "pretty print" */ 01186 enum {DEBUG_LOC=false}; 01187 if( DEBUG_LOC ) 01188 { 01189 fprintf(ioQQQ,"%li\t%li\t%.2e\t%.2e\n", 01190 nelem+1, 01191 nHi, 01192 t->Emis->Aul , 01193 t->Emis->opacity 01194 ); 01195 } 01196 } 01197 return; 01198 } 01199 01200 /* calculate radiative lifetime of an individual iso state */ 01201 STATIC double iso_state_lifetime( long ipISO, long nelem, long n, long l ) 01202 { 01203 /* >>refer hydro lifetimes Horbatsch, M. W., Horbatsch, M. and Hessels, E. A. 2005, JPhysB, 38, 1765 */ 01204 01205 double tau, t0, eps2; 01206 /* mass of electron */ 01207 double m = ELECTRON_MASS; 01208 /* nuclear mass */ 01209 double M = (double)dense.AtomicWeight[nelem] * ATOMIC_MASS_UNIT; 01210 double mu = (m*M)/(M+m); 01211 long z = 1; 01212 long Z = nelem + 1 - ipISO; 01213 01214 DEBUG_ENTRY( "iso_state_lifetime()" ); 01215 01216 eps2 = 1. - ( l*l + l + 8./47. - (l+1.)/69./n ) / POW2( (double)n ); 01217 01218 t0 = 3. * H_BAR * pow( (double)n, 5.) / 01219 ( 2. * POW4( (double)( z * Z ) ) * pow( FINE_STRUCTURE, 5. ) * mu * POW2( SPEEDLIGHT ) ) * 01220 POW2( (m + M)/(Z*m + z*M) ); 01221 01222 tau = t0 * ( 1. - eps2 ) / 01223 ( 1. + 19./88.*( (1./eps2 - 1.) * log( 1. - eps2 ) + 1. - 01224 0.5 * eps2 - 0.025 * eps2 * eps2 ) ); 01225 01226 if( ipISO == ipHE_LIKE ) 01227 { 01228 /* iso_state_lifetime is not spin specific, must exclude helike triplet here. */ 01229 tau /= 3.; 01230 /* this is also necessary to correct the helike lifetimes */ 01231 tau *= 1.1722 * pow( (double)nelem, 0.1 ); 01232 } 01233 01234 /* would probably need a new lifetime algorithm for any other iso sequences. */ 01235 ASSERT( ipISO <= ipHE_LIKE ); 01236 ASSERT( tau > 0. ); 01237 01238 return tau; 01239 } 01240 01241 /* calculate cascade probabilities, branching ratios, and associated errors. */ 01242 void iso_cascade( long ipISO, long nelem ) 01243 { 01244 /* The sum of all A's coming out of a given n, 01245 * Below we assert a monotonic trend. */ 01246 double *SumAPerN; 01247 01248 long int i, j, ipLo, ipHi; 01249 01250 DEBUG_ENTRY( "iso_cascade()" ); 01251 01252 SumAPerN = ((double*)MALLOC( (size_t)(iso.n_HighestResolved_max[ipISO][nelem] + iso.nCollapsed_max[ipISO][nelem] + 1 )*sizeof(double ))); 01253 memset( SumAPerN, 0, (iso.n_HighestResolved_max[ipISO][nelem] + iso.nCollapsed_max[ipISO][nelem] + 1 )*sizeof(double ) ); 01254 01255 /* Initialize some ground state stuff, easier here than in loops. */ 01256 iso.CascadeProb[ipISO][nelem][0][0] = 1.; 01257 if( iso.lgRandErrGen[ipISO] ) 01258 { 01259 iso.SigmaAtot[ipISO][nelem][0] = 0.; 01260 iso.SigmaCascadeProb[ipISO][nelem][0][0] = 0.; 01261 } 01262 01263 /***************************************************************************/ 01264 /****** Cascade probabilities, Branching ratios, and associated errors *****/ 01265 /***************************************************************************/ 01266 for( ipHi=1; ipHi<iso.numLevels_max[ipISO][nelem]; ipHi++ ) 01267 { 01268 double SumAs = 0.; 01269 01275 /* initialize variables. */ 01276 iso.CascadeProb[ipISO][nelem][ipHi][ipHi] = 1.; 01277 iso.CascadeProb[ipISO][nelem][ipHi][0] = 0.; 01278 iso.BranchRatio[ipISO][nelem][ipHi][0] = 0.; 01279 01280 if( iso.lgRandErrGen[ipISO] ) 01281 { 01282 iso.SigmaAtot[ipISO][nelem][ipHi] = 0.; 01283 iso.SigmaCascadeProb[ipISO][nelem][ipHi][ipHi] = 0.; 01284 } 01285 01286 long ipLoStart = 0; 01287 if( opac.lgCaseB && L_(ipHi)==1 && (ipISO==ipH_LIKE || S_(ipHi)==1) ) 01288 ipLoStart = 1; 01289 01290 for( ipLo=ipLoStart; ipLo<ipHi; ipLo++ ) 01291 { 01292 SumAs += Transitions[ipISO][nelem][ipHi][ipLo].Emis->Aul; 01293 } 01294 01295 for( ipLo=ipLoStart; ipLo<ipHi; ipLo++ ) 01296 { 01297 if( Transitions[ipISO][nelem][ipHi][ipLo].Emis->Aul <= iso.SmallA ) 01298 { 01299 iso.CascadeProb[ipISO][nelem][ipHi][ipLo] = 0.; 01300 iso.BranchRatio[ipISO][nelem][ipHi][ipLo] = 0.; 01301 continue; 01302 } 01303 01304 iso.CascadeProb[ipISO][nelem][ipHi][ipLo] = 0.; 01305 iso.BranchRatio[ipISO][nelem][ipHi][ipLo] = 01306 Transitions[ipISO][nelem][ipHi][ipLo].Emis->Aul / SumAs; 01307 01308 ASSERT( iso.BranchRatio[ipISO][nelem][ipHi][ipLo] <= 1.0000001 ); 01309 01310 SumAPerN[N_(ipHi)] += Transitions[ipISO][nelem][ipHi][ipLo].Emis->Aul; 01311 01312 /* there are some negative energy transitions, where the order 01313 * has changed, but these are not optically allowed, these are 01314 * same n, different L, forbidden transitions */ 01315 ASSERT( Transitions[ipISO][nelem][ipHi][ipLo].EnergyWN > 0. || 01316 Transitions[ipISO][nelem][ipHi][ipLo].Emis->Aul <= iso.SmallA ); 01317 01318 if( iso.lgRandErrGen[ipISO] ) 01319 { 01320 ASSERT( iso.Error[ipISO][nelem][ipHi][ipLo][IPRAD] >= 0. ); 01321 /* Uncertainties in A's are added in quadrature, square root is taken below. */ 01322 iso.SigmaAtot[ipISO][nelem][ipHi] += 01323 pow( iso.Error[ipISO][nelem][ipHi][ipLo][IPRAD] * 01324 (double)Transitions[ipISO][nelem][ipHi][ipLo].Emis->Aul, 2. ); 01325 } 01326 } 01327 01328 if( iso.lgRandErrGen[ipISO] ) 01329 { 01330 /* Uncertainties in A's are added in quadrature above, square root taken here. */ 01331 iso.SigmaAtot[ipISO][nelem][ipHi] = sqrt( iso.SigmaAtot[ipISO][nelem][ipHi] ); 01332 } 01333 01334 /* cascade probabilities */ 01335 for( ipLo=0; ipLo<ipHi; ipLo++ ) 01336 { 01337 double SigmaA, SigmaCul = 0.; 01338 01339 for( i=ipLo; i<ipHi; i++ ) 01340 { 01341 iso.CascadeProb[ipISO][nelem][ipHi][ipLo] += iso.BranchRatio[ipISO][nelem][ipHi][i] * iso.CascadeProb[ipISO][nelem][i][ipLo]; 01342 if( iso.lgRandErrGen[ipISO] ) 01343 { 01344 if( Transitions[ipISO][nelem][ipHi][i].Emis->Aul > iso.SmallA ) 01345 { 01346 /* Uncertainties in A's and cascade probabilities */ 01347 SigmaA = iso.Error[ipISO][nelem][ipHi][i][IPRAD] * 01348 Transitions[ipISO][nelem][ipHi][i].Emis->Aul; 01349 SigmaCul += 01350 pow(SigmaA*iso.CascadeProb[ipISO][nelem][i][ipLo]*StatesElem[ipISO][nelem][ipHi].lifetime, 2.) + 01351 pow(iso.SigmaAtot[ipISO][nelem][ipHi]*iso.BranchRatio[ipISO][nelem][ipHi][i]*iso.CascadeProb[ipISO][nelem][i][ipLo]*StatesElem[ipISO][nelem][ipHi].lifetime, 2.) + 01352 pow(iso.SigmaCascadeProb[ipISO][nelem][i][ipLo]*iso.BranchRatio[ipISO][nelem][ipHi][i], 2.); 01353 } 01354 } 01355 } 01356 if( iso.lgRandErrGen[ipISO] ) 01357 { 01358 SigmaCul = sqrt(SigmaCul); 01359 iso.SigmaCascadeProb[ipISO][nelem][ipHi][ipLo] = SigmaCul; 01360 } 01361 } 01362 } 01363 01364 /************************************************************************/ 01365 /*** Allowed decay conversion probabilities. See Robbins68b, Table 1. ***/ 01366 /************************************************************************/ 01367 { 01368 enum {DEBUG_LOC=false}; 01369 01370 if( DEBUG_LOC && (nelem == ipHELIUM) && (ipISO==ipHE_LIKE) ) 01371 { 01372 /* To output Bm(n,l; ipLo), set ipLo, hi_l, and hi_s accordingly. */ 01373 long int hi_l,hi_s; 01374 double Bm; 01375 01376 /* these must be set for following output to make sense 01377 * as is, a dangerous bit of code - set NaN for safety */ 01378 hi_s = -100000; 01379 hi_l = -100000; 01380 ipLo = -100000; 01381 /* tripS to 2^3P */ 01382 //hi_l = 0, hi_s = 3, ipLo = ipHe2p3P0; 01383 01384 /* tripD to 2^3P */ 01385 //hi_l = 2, hi_s = 3, ipLo = ipHe2p3P0; 01386 01387 /* tripP to 2^3S */ 01388 //hi_l = 1, hi_s = 3, ipLo = ipHe2s3S; 01389 01390 ASSERT( hi_l != StatesElem[ipISO][nelem][ipLo].l ); 01391 01392 fprintf(ioQQQ,"Bm(n,%ld,%ld;%ld)\n",hi_l,hi_s,ipLo); 01393 fprintf(ioQQQ,"m\t2\t\t3\t\t4\t\t5\t\t6\n"); 01394 01395 for( ipHi=ipHe2p3P2; ipHi<iso.numLevels_max[ipISO][nelem]-iso.nCollapsed_max[ipISO][nelem]; ipHi++ ) 01396 { 01397 /* Pick out excitations from metastable 2tripS to ntripP. */ 01398 if( (StatesElem[ipISO][nelem][ipHi].l == 1) && (StatesElem[ipISO][nelem][ipHi].S == 3) ) 01399 { 01400 fprintf(ioQQQ,"\n%ld\t",StatesElem[ipISO][nelem][ipHi].n); 01401 j = 0; 01402 Bm = 0; 01403 for( i = ipLo; i<=ipHi; i++) 01404 { 01405 if( (StatesElem[ipISO][nelem][i].l == hi_l) && (StatesElem[ipISO][nelem][i].S == hi_s) ) 01406 { 01407 if( (ipLo == ipHe2p3P0) && (i > ipHe2p3P2) ) 01408 { 01409 Bm += iso.CascadeProb[ipISO][nelem][ipHi][i] * ( iso.BranchRatio[ipISO][nelem][i][ipHe2p3P0] + 01410 iso.BranchRatio[ipISO][nelem][i][ipHe2p3P1] + iso.BranchRatio[ipISO][nelem][i][ipHe2p3P2] ); 01411 } 01412 else 01413 Bm += iso.CascadeProb[ipISO][nelem][ipHi][i] * iso.BranchRatio[ipISO][nelem][i][ipLo]; 01414 01415 if( (i == ipHe2p3P0) || (i == ipHe2p3P1) || (i == ipHe2p3P2) ) 01416 { 01417 j++; 01418 if(j == 3) 01419 { 01420 fprintf(ioQQQ,"%2.4e\t",Bm); 01421 Bm = 0; 01422 } 01423 } 01424 else 01425 { 01426 fprintf(ioQQQ,"%2.4e\t",Bm); 01427 Bm = 0; 01428 } 01429 } 01430 } 01431 } 01432 } 01433 fprintf(ioQQQ,"\n\n"); 01434 } 01435 } 01436 01437 /******************************************************/ 01438 /*** Lifetimes should increase monotonically with ***/ 01439 /*** increasing n...Make sure the A's decrease. ***/ 01440 /******************************************************/ 01441 for( i=2; i < iso.n_HighestResolved_max[ipISO][nelem] + iso.nCollapsed_max[ipISO][nelem] - 1; ++i) 01442 { 01443 /* exclude the first collapsed level, since different treatments can cause discontinuity */ 01444 ASSERT( (SumAPerN[i] > SumAPerN[i+1]) || opac.lgCaseB || (i==iso.n_HighestResolved_max[ipISO][nelem]+1) ); 01445 } 01446 01447 { 01448 enum {DEBUG_LOC=false}; 01449 if( DEBUG_LOC /* && (ipISO == ipH_LIKE) && (nelem == ipHYDROGEN) */) 01450 { 01451 for( i = 2; i<= (iso.n_HighestResolved_max[ipISO][nelem] + iso.nCollapsed_max[ipISO][nelem]); ++i) 01452 { 01453 fprintf(ioQQQ,"n %ld\t lifetime %.4e\n", i, 1./SumAPerN[i]); 01454 } 01455 } 01456 } 01457 01458 free( SumAPerN ); 01459 01460 return; 01461 } 01462 01464 /* For double-ionization discussions, see Lindsay, Rejoub, & Stebbings 2002 */ 01465 /* Also read Itza-Ortiz, Godunov, Wang, and McGuire 2001. */ 01466 STATIC void iso_satellite( void ) 01467 { 01468 long i, ipISO, nelem; 01469 01470 DEBUG_ENTRY( "iso_satellite()" ); 01471 01472 for( ipISO = ipHE_LIKE; ipISO < NISO; ipISO++ ) 01473 { 01474 for( nelem = ipISO; nelem < LIMELM; nelem++ ) 01475 { 01476 if( !iso.lgDielRecom[ipISO] ) 01477 continue; 01478 01479 /* must always do helium even if turned off */ 01480 if( nelem == ipISO || dense.lgElmtOn[nelem] ) 01481 { 01482 for( i=0; i<iso.numLevels_max[ipISO][nelem]; i++ ) 01483 { 01484 char chLab[5]=" "; 01485 01486 TransitionZero( &SatelliteLines[ipISO][nelem][i] ); 01487 01488 /* Make approximation that all levels have energy of H-like 2s level */ 01489 /* Lines to 1s2s have roughly energy of parent Ly-alpha. So lines to 1snL will have an energy 01490 * smaller by the difference between nL and 2s energies. Therefore, the following has 01491 * energy of parent Ly-alpha MINUS the difference between daughter level and daughter n=2 level. */ 01492 SatelliteLines[ipISO][nelem][i].WLAng = (realnum)(RYDLAM/ 01493 ((iso.xIsoLevNIonRyd[ipISO-1][nelem][0] - iso.xIsoLevNIonRyd[ipISO-1][nelem][1]) - 01494 (iso.xIsoLevNIonRyd[ipISO][nelem][1]- iso.xIsoLevNIonRyd[ipISO][nelem][i])) ); 01495 01496 SatelliteLines[ipISO][nelem][i].EnergyWN = 1.e8f / 01497 SatelliteLines[ipISO][nelem][i].WLAng; 01498 01499 SatelliteLines[ipISO][nelem][i].EnergyErg = (realnum)ERG1CM * 01500 SatelliteLines[ipISO][nelem][i].EnergyWN; 01501 01502 SatelliteLines[ipISO][nelem][i].EnergyK = (realnum)T1CM * 01503 SatelliteLines[ipISO][nelem][i].EnergyWN; 01504 01505 /* generate label for this ion */ 01506 sprintf( chLab, "%2s%2ld",elementnames.chElementSym[nelem], nelem+1-ipISO ); 01507 01508 SatelliteLines[ipISO][nelem][i].Emis->iRedisFun = ipCRDW; 01509 /* this is not the usual nelem, is it atomic not C scale. */ 01510 SatelliteLines[ipISO][nelem][i].Hi->nelem = nelem + 1; 01511 SatelliteLines[ipISO][nelem][i].Hi->IonStg = nelem + 1 - ipISO; 01512 fixit(); /* what should the stat weight of the upper level be? For now say 2. */ 01513 SatelliteLines[ipISO][nelem][i].Hi->g = 2.f; 01514 SatelliteLines[ipISO][nelem][i].Lo->g = StatesElem[ipISO][nelem][i].g; 01515 SatelliteLines[ipISO][nelem][i].Emis->PopOpc = 01516 SatelliteLines[ipISO][nelem][i].Lo->Pop; 01517 } 01518 } 01519 } 01520 } 01521 01522 return; 01523 } 01524 01525 void iso_satellite_update( void ) 01526 { 01527 double ConBoltz, LTE_pop=SMALLFLOAT+FLT_EPSILON, factor1, ConvLTEPOP; 01528 long i, ipISO, nelem; 01529 double dr_rate; 01530 01531 DEBUG_ENTRY( "iso_satellite()" ); 01532 01533 for( ipISO = ipHE_LIKE; ipISO < NISO; ipISO++ ) 01534 { 01535 for( nelem = ipISO; nelem < LIMELM; nelem++ ) 01536 { 01537 if( !iso.lgDielRecom[ipISO] ) 01538 continue; 01539 01540 /* must always do helium even if turned off */ 01541 if( nelem == ipISO || dense.lgElmtOn[nelem] ) 01542 { 01543 for( i=0; i<iso.numLevels_max[ipISO][nelem]; i++ ) 01544 { 01545 dr_rate = MAX2( iso.SmallA, iso.DielecRecomb[ipISO][nelem][i] ); 01546 01547 SatelliteLines[ipISO][nelem][i].Emis->phots = 01548 dr_rate * dense.eden * dense.xIonDense[nelem][nelem+1-ipISO]; 01549 01550 SatelliteLines[ipISO][nelem][i].Emis->xIntensity = 01551 SatelliteLines[ipISO][nelem][i].Emis->phots * 01552 ERG1CM * SatelliteLines[ipISO][nelem][i].EnergyWN; 01553 01554 /* We set line intensity above using a rate, but here we need a transition probability. 01555 * We can obtain this by dividing dr_rate by the population of the autoionizing level. 01556 * We assume this level is in statistical equilibrium. */ 01557 factor1 = HION_LTE_POP*dense.AtomicWeight[nelem]/ 01558 (dense.AtomicWeight[nelem]+ELECTRON_MASS/ATOMIC_MASS_UNIT); 01559 01560 /* term in () is stat weight of electron * ion */ 01561 ConvLTEPOP = pow(factor1,1.5)/(2.*iso.stat_ion[ipISO])/phycon.te32; 01562 01563 /* This Boltzmann factor is exp( +ioniz energy / Te ). For simplicity, we make 01564 * the fair approximation that all of the autoionizing levels have an energy 01565 * equal to the parents n=2. */ 01566 ConBoltz = dsexp(iso.xIsoLevNIonRyd[ipISO-1][nelem][1]/phycon.te_ryd); 01567 01568 if( ConBoltz >= SMALLDOUBLE ) 01569 { 01570 /* The energy used to calculate ConBoltz above 01571 * should be negative since this is above the continuum, but 01572 * to be safe we calculate ConBoltz with a positive energy above 01573 * and multiply by it here instead of dividing. */ 01574 LTE_pop = SatelliteLines[ipISO][nelem][i].Hi->g * ConBoltz * ConvLTEPOP; 01575 } 01576 01577 LTE_pop = max( LTE_pop, 1e-30f ); 01578 01579 /* Now the transtion probability is simply dr_rate/LTE_pop. */ 01580 SatelliteLines[ipISO][nelem][i].Emis->Aul = (realnum)(dr_rate/LTE_pop); 01581 SatelliteLines[ipISO][nelem][i].Emis->Aul = 01582 max( iso.SmallA, SatelliteLines[ipISO][nelem][i].Emis->Aul ); 01583 01584 SatelliteLines[ipISO][nelem][i].Emis->gf = (realnum)GetGF( 01585 SatelliteLines[ipISO][nelem][i].Emis->Aul, 01586 SatelliteLines[ipISO][nelem][i].EnergyWN, 01587 SatelliteLines[ipISO][nelem][i].Hi->g); 01588 01589 SatelliteLines[ipISO][nelem][i].Emis->gf = 01590 max( 1e-20f, SatelliteLines[ipISO][nelem][i].Emis->gf ); 01591 01592 SatelliteLines[ipISO][nelem][i].Emis->PopOpc = 01593 SatelliteLines[ipISO][nelem][i].Lo->Pop - 01594 LTE_pop * SatelliteLines[ipISO][nelem][i].Lo->g/SatelliteLines[ipISO][nelem][i].Hi->g; 01595 01596 SatelliteLines[ipISO][nelem][i].Emis->opacity = 01597 (realnum)(abscf(SatelliteLines[ipISO][nelem][i].Emis->gf, 01598 SatelliteLines[ipISO][nelem][i].EnergyWN, 01599 SatelliteLines[ipISO][nelem][i].Lo->g)); 01600 01601 /* a typical transition probability is of order 1e10 s-1 */ 01602 double lifetime = 1e-10; 01603 01604 SatelliteLines[ipISO][nelem][i].Emis->dampXvel = (realnum)( 01605 (1.f/lifetime)/PI4/SatelliteLines[ipISO][nelem][i].EnergyWN); 01606 } 01607 } 01608 } 01609 } 01610 01611 return; 01612 } 01613 01614 long iso_get_total_num_levels( long ipISO, long nmaxResolved, long numCollapsed ) 01615 { 01616 DEBUG_ENTRY( "iso_get_total_num_levels()" ); 01617 01618 long tot_num_levels; 01619 01620 /* return the number of levels up to and including nmaxResolved PLUS 01621 * the number (numCollapsed) of collapsed n-levels */ 01622 01623 if( ipISO == ipH_LIKE ) 01624 { 01625 tot_num_levels = (long)( nmaxResolved * 0.5 *( nmaxResolved + 1 ) ) + numCollapsed; 01626 } 01627 else if( ipISO == ipHE_LIKE ) 01628 { 01629 tot_num_levels = nmaxResolved*nmaxResolved + nmaxResolved + 1 + numCollapsed; 01630 } 01631 else 01632 TotalInsanity(); 01633 01634 return tot_num_levels; 01635 } 01636 01637 void iso_update_num_levels( long ipISO, long nelem ) 01638 { 01639 DEBUG_ENTRY( "iso_update_num_levels()" ); 01640 01641 /* This is the minimum resolved nmax. */ 01642 ASSERT( iso.n_HighestResolved_max[ipISO][nelem] >= 3 ); 01643 01644 iso.numLevels_max[ipISO][nelem] = 01645 iso_get_total_num_levels( ipISO, iso.n_HighestResolved_max[ipISO][nelem], iso.nCollapsed_max[ipISO][nelem] ); 01646 01647 if( iso.numLevels_max[ipISO][nelem] > iso.numLevels_malloc[ipISO][nelem] ) 01648 { 01649 fprintf( ioQQQ, "The number of levels for ipISO %li, nelem %li, has been increased since the initial coreload.\n", 01650 ipISO, nelem ); 01651 fprintf( ioQQQ, "This cannot be done.\n" ); 01652 cdEXIT(EXIT_FAILURE); 01653 } 01654 01655 /* set local copies to the max values */ 01656 iso.numLevels_local[ipISO][nelem] = iso.numLevels_max[ipISO][nelem]; 01657 iso.nCollapsed_local[ipISO][nelem] = iso.nCollapsed_max[ipISO][nelem]; 01658 iso.n_HighestResolved_local[ipISO][nelem] = iso.n_HighestResolved_max[ipISO][nelem]; 01659 01660 /* only print results from resolved levels. */ 01661 iso.numPrintLevels[ipISO][nelem] = iso.numLevels_max[ipISO][nelem] - iso.nCollapsed_max[ipISO][nelem]; 01662 01663 /* find the largest number of levels in any element in all iso sequences 01664 * we will allocate one matrix for ionization solver, and just use a piece of that memory 01665 * for smaller models. */ 01666 max_num_levels = MAX2( max_num_levels, iso.numLevels_max[ipISO][nelem]); 01667 01668 return; 01669 } 01670 01671 void iso_collapsed_bnl_set( long ipISO, long nelem ) 01672 { 01673 01674 DEBUG_ENTRY( "iso_collapsed_bnl_set()" ); 01675 01676 double bnl_array[4][3][4][10] = { 01677 { 01678 /* H */ 01679 { 01680 {6.13E-01, 2.56E-01, 1.51E-01, 2.74E-01, 3.98E-01, 4.98E-01, 5.71E-01, 6.33E-01, 7.28E-01, 9.59E-01}, 01681 {1.31E+00, 5.17E-01, 2.76E-01, 4.47E-01, 5.87E-01, 6.82E-01, 7.44E-01, 8.05E-01, 9.30E-01, 1.27E+00}, 01682 {1.94E+00, 7.32E-01, 3.63E-01, 5.48E-01, 6.83E-01, 7.66E-01, 8.19E-01, 8.80E-01, 1.02E+00, 1.43E+00}, 01683 {2.53E+00, 9.15E-01, 4.28E-01, 6.16E-01, 7.42E-01, 8.13E-01, 8.60E-01, 9.22E-01, 1.08E+00, 1.56E+00} 01684 }, 01685 { 01686 {5.63E-01, 2.65E-01, 1.55E-01, 2.76E-01, 3.91E-01, 4.75E-01, 5.24E-01, 5.45E-01, 5.51E-01, 5.53E-01}, 01687 {1.21E+00, 5.30E-01, 2.81E-01, 4.48E-01, 5.80E-01, 6.62E-01, 7.05E-01, 7.24E-01, 7.36E-01, 7.46E-01}, 01688 {1.81E+00, 7.46E-01, 3.68E-01, 5.49E-01, 6.78E-01, 7.51E-01, 7.88E-01, 8.09E-01, 8.26E-01, 8.43E-01}, 01689 {2.38E+00, 9.27E-01, 4.33E-01, 6.17E-01, 7.38E-01, 8.05E-01, 8.40E-01, 8.65E-01, 8.92E-01, 9.22E-01} 01690 }, 01691 { 01692 {2.97E-01, 2.76E-01, 2.41E-01, 3.04E-01, 3.66E-01, 4.10E-01, 4.35E-01, 4.48E-01, 4.52E-01, 4.53E-01}, 01693 {5.63E-01, 5.04E-01, 3.92E-01, 4.67E-01, 5.39E-01, 5.85E-01, 6.10E-01, 6.20E-01, 6.23E-01, 6.23E-01}, 01694 {7.93E-01, 6.90E-01, 4.94E-01, 5.65E-01, 6.36E-01, 6.79E-01, 7.00E-01, 7.09E-01, 7.11E-01, 7.11E-01}, 01695 {1.04E+00, 8.66E-01, 5.62E-01, 6.31E-01, 7.01E-01, 7.43E-01, 7.63E-01, 7.71E-01, 7.73E-01, 7.73E-01} 01696 } 01697 }, 01698 { 01699 /* He+ */ 01700 { 01701 {6.70E-02, 2.93E-02, 1.94E-02, 4.20E-02, 7.40E-02, 1.12E-01, 1.51E-01, 1.86E-01, 2.26E-01, 3.84E-01}, 01702 {2.39E-01, 1.03E-01, 6.52E-02, 1.31E-01, 2.11E-01, 2.91E-01, 3.61E-01, 4.17E-01, 4.85E-01, 8.00E-01}, 01703 {4.26E-01, 1.80E-01, 1.10E-01, 2.09E-01, 3.18E-01, 4.15E-01, 4.93E-01, 5.54E-01, 6.34E-01, 1.04E+00}, 01704 {6.11E-01, 2.55E-01, 1.51E-01, 2.74E-01, 3.99E-01, 5.02E-01, 5.80E-01, 6.41E-01, 7.30E-01, 1.21E+00} 01705 }, 01706 { 01707 {6.79E-02, 3.00E-02, 2.00E-02, 4.30E-02, 7.48E-02, 1.11E-01, 1.44E-01, 1.70E-01, 1.87E-01, 1.96E-01}, 01708 {2.40E-01, 1.04E-01, 6.62E-02, 1.32E-01, 2.11E-01, 2.87E-01, 3.51E-01, 3.98E-01, 4.32E-01, 4.58E-01}, 01709 {4.26E-01, 1.81E-01, 1.11E-01, 2.10E-01, 3.17E-01, 4.12E-01, 4.89E-01, 5.53E-01, 6.14E-01, 6.84E-01}, 01710 {6.12E-01, 2.55E-01, 1.51E-01, 2.73E-01, 3.97E-01, 4.98E-01, 5.77E-01, 6.51E-01, 7.82E-01, 1.18E+00} 01711 }, 01712 { 01713 {4.98E-02, 3.47E-02, 2.31E-02, 4.54E-02, 7.14E-02, 9.37E-02, 1.08E-01, 1.13E-01, 1.13E-01, 1.11E-01}, 01714 {1.75E-01, 1.16E-01, 7.36E-02, 1.36E-01, 2.01E-01, 2.50E-01, 2.76E-01, 2.84E-01, 2.81E-01, 2.77E-01}, 01715 {3.38E-01, 1.97E-01, 1.18E-01, 2.13E-01, 3.06E-01, 3.72E-01, 4.06E-01, 4.15E-01, 4.10E-01, 4.04E-01}, 01716 {6.01E-01, 2.60E-01, 1.53E-01, 2.76E-01, 3.95E-01, 4.87E-01, 5.45E-01, 5.76E-01, 5.93E-01, 6.05E-01} 01717 } 01718 }, 01719 { 01720 /* He singlets */ 01721 { 01722 {1.77E-01, 3.59E-01, 1.54E-01, 2.75E-01, 3.98E-01, 4.94E-01, 5.51E-01, 5.68E-01, 5.46E-01, 4.97E-01}, 01723 {4.09E-01, 7.23E-01, 2.83E-01, 4.48E-01, 5.89E-01, 6.78E-01, 7.22E-01, 7.30E-01, 7.07E-01, 6.65E-01}, 01724 {6.40E-01, 1.02E+00, 3.74E-01, 5.49E-01, 6.85E-01, 7.63E-01, 7.98E-01, 8.03E-01, 7.84E-01, 7.53E-01}, 01725 {8.70E-01, 1.28E+00, 4.42E-01, 6.17E-01, 7.44E-01, 8.13E-01, 8.42E-01, 8.46E-01, 8.34E-01, 8.13E-01} 01726 }, 01727 { 01728 {1.78E-01, 3.62E-01, 1.55E-01, 2.73E-01, 3.91E-01, 4.73E-01, 5.10E-01, 5.04E-01, 4.70E-01, 4.32E-01}, 01729 {4.08E-01, 7.26E-01, 2.83E-01, 4.45E-01, 5.79E-01, 6.54E-01, 6.78E-01, 6.64E-01, 6.30E-01, 5.98E-01}, 01730 {6.37E-01, 1.03E+00, 3.73E-01, 5.46E-01, 6.75E-01, 7.40E-01, 7.57E-01, 7.43E-01, 7.15E-01, 6.92E-01}, 01731 {8.65E-01, 1.28E+00, 4.41E-01, 6.14E-01, 7.35E-01, 7.92E-01, 8.05E-01, 7.95E-01, 7.74E-01, 7.59E-01} 01732 }, 01733 { 01734 {2.07E-01, 3.73E-01, 1.73E-01, 2.85E-01, 4.03E-01, 4.76E-01, 5.06E-01, 5.03E-01, 4.84E-01, 4.63E-01}, 01735 {4.32E-01, 7.13E-01, 3.06E-01, 4.54E-01, 5.81E-01, 6.44E-01, 6.59E-01, 6.49E-01, 6.28E-01, 6.11E-01}, 01736 {6.40E-01, 9.85E-01, 3.98E-01, 5.53E-01, 6.74E-01, 7.27E-01, 7.36E-01, 7.26E-01, 7.10E-01, 6.98E-01}, 01737 {8.38E-01, 1.21E+00, 4.67E-01, 6.20E-01, 7.34E-01, 7.79E-01, 7.87E-01, 7.79E-01, 7.69E-01, 7.63E-01} 01738 } 01739 }, 01740 { 01741 /* He triplets */ 01742 { 01743 {9.31E-02, 3.96E-01, 1.36E-01, 2.74E-01, 3.99E-01, 4.95E-01, 5.52E-01, 5.70E-01, 5.48E-01, 4.96E-01}, 01744 {2.25E-01, 8.46E-01, 2.49E-01, 4.46E-01, 5.89E-01, 6.79E-01, 7.23E-01, 7.31E-01, 7.08E-01, 6.64E-01}, 01745 {3.59E-01, 1.24E+00, 3.30E-01, 5.47E-01, 6.85E-01, 7.63E-01, 7.98E-01, 8.04E-01, 7.85E-01, 7.53E-01}, 01746 {4.93E-01, 1.60E+00, 3.91E-01, 6.15E-01, 7.44E-01, 8.13E-01, 8.42E-01, 8.47E-01, 8.35E-01, 8.12E-01} 01747 }, 01748 { 01749 {9.32E-02, 3.99E-01, 1.35E-01, 2.72E-01, 3.91E-01, 4.75E-01, 5.14E-01, 5.09E-01, 4.73E-01, 4.31E-01}, 01750 {2.25E-01, 8.49E-01, 2.49E-01, 4.44E-01, 5.79E-01, 6.56E-01, 6.81E-01, 6.68E-01, 6.31E-01, 5.96E-01}, 01751 {3.58E-01, 1.25E+00, 3.29E-01, 5.44E-01, 6.76E-01, 7.42E-01, 7.60E-01, 7.46E-01, 7.16E-01, 6.91E-01}, 01752 {4.92E-01, 1.60E+00, 3.90E-01, 6.12E-01, 7.36E-01, 7.93E-01, 8.07E-01, 7.97E-01, 7.74E-01, 7.58E-01} 01753 }, 01754 { 01755 {1.13E-01, 4.21E-01, 1.47E-01, 2.83E-01, 4.04E-01, 4.80E-01, 5.13E-01, 5.12E-01, 4.93E-01, 4.71E-01}, 01756 {2.52E-01, 8.56E-01, 2.61E-01, 4.50E-01, 5.82E-01, 6.48E-01, 6.66E-01, 6.56E-01, 6.35E-01, 6.16E-01}, 01757 {3.85E-01, 1.23E+00, 3.41E-01, 5.49E-01, 6.75E-01, 7.30E-01, 7.41E-01, 7.31E-01, 7.15E-01, 7.02E-01}, 01758 {5.14E-01, 1.56E+00, 4.01E-01, 6.15E-01, 7.34E-01, 7.82E-01, 7.90E-01, 7.83E-01, 7.72E-01, 7.65E-01} 01759 } 01760 } 01761 }; 01762 01763 double temps[4] = {5000., 10000., 15000., 20000. }; 01764 double log_dens[3] = {2., 4., 6.}; 01765 long ipTe, ipDens; 01766 01767 ASSERT( nelem <= 1 ); 01768 01769 /* find temperature in tabulated values. */ 01770 ipTe = hunt_bisect( temps, 4, phycon.te ); 01771 ipDens = hunt_bisect( log_dens, 3, log10(dense.eden) ); 01772 01773 ASSERT( (ipTe >=0) && (ipTe < 3) ); 01774 ASSERT( (ipDens >=0) && (ipDens < 2) ); 01775 01776 for( long nHi=iso.n_HighestResolved_max[ipISO][nelem]+1; nHi<=iso.n_HighestResolved_max[ipISO][nelem]+iso.nCollapsed_max[ipISO][nelem]; nHi++ ) 01777 { 01778 for( long lHi=0; lHi<nHi; lHi++ ) 01779 { 01780 for( long sHi=1; sHi<4; sHi++ ) 01781 { 01782 if( ipISO == ipH_LIKE && sHi != 2 ) 01783 continue; 01784 else if( ipISO == ipHE_LIKE && sHi != 1 && sHi != 3 ) 01785 continue; 01786 01787 double bnl_at_lo_den, bnl_at_hi_den, bnl; 01788 double bnl_max, bnl_min, temp, dens; 01789 01790 long ipL = MIN2(9,lHi); 01791 long ip1; 01792 01793 if( nelem==ipHYDROGEN ) 01794 ip1 = 0; 01795 else if( nelem==ipHELIUM ) 01796 { 01797 if( ipISO==ipH_LIKE ) 01798 ip1 = 1; 01799 else if( ipISO==ipHE_LIKE ) 01800 { 01801 if( sHi==1 ) 01802 ip1 = 2; 01803 else if( sHi==3 ) 01804 ip1 = 3; 01805 else 01806 TotalInsanity(); 01807 } 01808 else 01809 TotalInsanity(); 01810 } 01811 else 01812 TotalInsanity(); 01813 01814 temp = MAX2( temps[0], phycon.te ); 01815 temp = MIN2( temps[3], temp ); 01816 01817 dens = MAX2( log_dens[0], log10(dense.eden) ); 01818 dens = MIN2( log_dens[2], dens ); 01819 01820 /* Calculate the answer...must interpolate on two variables. 01821 * First interpolate on T, at both the lower and upper densities. 01822 * Then interpolate between these results for the right density. */ 01823 01824 if( temp < temps[0] && dens < log_dens[0] ) 01825 bnl = bnl_array[ip1][0][0][ipL]; 01826 else if( temp < temps[0] && dens >= log_dens[2] ) 01827 bnl = bnl_array[ip1][2][0][ipL]; 01828 else if( temp >= temps[3] && dens < log_dens[0] ) 01829 bnl = bnl_array[ip1][0][3][ipL]; 01830 else if( temp >= temps[3] && dens >= log_dens[2] ) 01831 bnl = bnl_array[ip1][2][3][ipL]; 01832 else 01833 { 01834 bnl_at_lo_den = ( temp - temps[ipTe]) / (temps[ipTe+1] - temps[ipTe]) * 01835 (bnl_array[ip1][ipDens][ipTe+1][ipL] - bnl_array[ip1][ipDens][ipTe][ipL]) + bnl_array[ip1][ipDens][ipTe][ipL]; 01836 01837 bnl_at_hi_den = ( temp - temps[ipTe]) / (temps[ipTe+1] - temps[ipTe]) * 01838 (bnl_array[ip1][ipDens+1][ipTe+1][ipL] - bnl_array[ip1][ipDens+1][ipTe][ipL]) + bnl_array[ip1][ipDens+1][ipTe][ipL]; 01839 01840 bnl = ( dens - log_dens[ipDens]) / (log_dens[ipDens+1] - log_dens[ipDens]) * 01841 (bnl_at_hi_den - bnl_at_lo_den) + bnl_at_lo_den; 01842 } 01843 01845 { 01846 bnl_max = MAX4( bnl_array[ip1][ipDens][ipTe+1][ipL], bnl_array[ip1][ipDens+1][ipTe+1][ipL], 01847 bnl_array[ip1][ipDens][ipTe][ipL], bnl_array[ip1][ipDens+1][ipTe][ipL] ); 01848 ASSERT( bnl <= bnl_max ); 01849 01850 bnl_min = MIN4( bnl_array[ip1][ipDens][ipTe+1][ipL], bnl_array[ip1][ipDens+1][ipTe+1][ipL], 01851 bnl_array[ip1][ipDens][ipTe][ipL], bnl_array[ip1][ipDens+1][ipTe][ipL] ); 01852 ASSERT( bnl >= bnl_min ); 01853 } 01854 01855 iso.bnl_effective[ipISO][nelem][nHi][lHi][sHi] = bnl; 01856 01857 ASSERT( iso.bnl_effective[ipISO][nelem][nHi][lHi][sHi] > 0. ); 01858 } 01859 } 01860 } 01861 01862 return; 01863 } 01864 01865 01866 void iso_collapsed_bnl_print( long ipISO, long nelem ) 01867 { 01868 DEBUG_ENTRY( "iso_collapsed_bnl_print()" ); 01869 01870 for( long is = 1; is<=3; ++is) 01871 { 01872 if( ipISO == ipH_LIKE && is != 2 ) 01873 continue; 01874 else if( ipISO == ipHE_LIKE && is != 1 && is != 3 ) 01875 continue; 01876 01877 char chSpin[3][9]= {"singlets", "doublets", "triplets"}; 01878 01879 /* give element number and spin */ 01880 fprintf(ioQQQ," %s %s %s bnl\n", 01881 iso.chISO[ipISO], 01882 elementnames.chElementSym[nelem], 01883 chSpin[is-1]); 01884 01885 /* header with the l states */ 01886 fprintf(ioQQQ," n\\l=> "); 01887 for( long i =0; i < iso.n_HighestResolved_max[ipISO][nelem] + iso.nCollapsed_max[ipISO][nelem]; ++i) 01888 { 01889 fprintf(ioQQQ,"%2ld ",i); 01890 } 01891 fprintf(ioQQQ,"\n"); 01892 01893 /* loop over prin quant numbers, one per line, with l across */ 01894 for( long in = 1; in <= iso.n_HighestResolved_max[ipISO][nelem] + iso.nCollapsed_max[ipISO][nelem]; ++in) 01895 { 01896 if( is==3 && in==1 ) 01897 continue; 01898 01899 fprintf(ioQQQ," %2ld ",in); 01900 01901 for( long il = 0; il < in; ++il) 01902 { 01903 fprintf( ioQQQ, "%9.3e ", iso.bnl_effective[ipISO][nelem][in][il][is] ); 01904 } 01905 fprintf(ioQQQ,"\n"); 01906 } 01907 } 01908 01909 return; 01910 } 01911 01912 void iso_collapsed_Aul_update( long ipISO, long nelem ) 01913 { 01914 DEBUG_ENTRY( "iso_collapsed_Aul_update()" ); 01915 01916 long ipFirstCollapsed = iso.numLevels_max[ipISO][nelem] - iso.nCollapsed_max[ipISO][nelem]; 01917 01918 for( long ipLo=0; ipLo<ipFirstCollapsed; ipLo++ ) 01919 { 01920 long spin = StatesElem[ipISO][nelem][ipLo].S; 01921 01922 /* calculate effective Aul's from collapsed levels */ 01923 for( long nHi=iso.n_HighestResolved_max[ipISO][nelem]+1; nHi<=iso.n_HighestResolved_max[ipISO][nelem]+iso.nCollapsed_max[ipISO][nelem]; nHi++ ) 01924 { 01925 realnum Auls[2] = { 01926 iso.CachedAs[ipISO][nelem][ nHi-iso.n_HighestResolved_max[ipISO][nelem]-1 ][ ipLo ][0], 01927 iso.CachedAs[ipISO][nelem][ nHi-iso.n_HighestResolved_max[ipISO][nelem]-1 ][ ipLo ][1] }; 01928 01929 realnum EffectiveAul = 01930 Auls[0]*spin*(2.f*(L_(ipLo)+1.f)+1.f)*(realnum)iso.bnl_effective[ipISO][nelem][nHi][ L_(ipLo)+1 ][spin]; 01931 01932 /* this is for n,L-1 -> n',L 01933 * make sure L-1 exists. */ 01934 if( L_(ipLo) > 0 ) 01935 { 01936 EffectiveAul += 01937 Auls[1]*spin*(2.f*(L_(ipLo)-1.f)+1.f)*(realnum)iso.bnl_effective[ipISO][nelem][nHi][ L_(ipLo)-1 ][spin]; 01938 } 01939 01940 if( ipISO==ipH_LIKE ) 01941 EffectiveAul /= (2.f*nHi*nHi); 01942 else if( ipISO==ipHE_LIKE ) 01943 EffectiveAul /= (4.f*nHi*nHi); 01944 else 01945 TotalInsanity(); 01946 01947 long ipHi = iso.QuantumNumbers2Index[ipISO][nelem][nHi][ L_(ipLo)+1 ][spin]; 01948 01949 /* FINALLY, put the effective A in the proper Emis structure. */ 01950 Transitions[ipISO][nelem][ipHi][ipLo].Emis->Aul = EffectiveAul; 01951 01952 ASSERT( Transitions[ipISO][nelem][ipHi][ipLo].Emis->Aul > 0. ); 01953 } 01954 } 01955 01956 return; 01957 } 01958 01959 void iso_collapsed_lifetimes_update( long ipISO, long nelem ) 01960 { 01961 DEBUG_ENTRY( "iso_collapsed_Aul_update()" ); 01962 01963 for( long ipHi=iso.numLevels_max[ipISO][nelem]- iso.nCollapsed_max[ipISO][nelem]; ipHi < iso.numLevels_max[ipISO][nelem]; ipHi++ ) 01964 { 01965 StatesElem[ipISO][nelem][ipHi].lifetime = SMALLFLOAT; 01966 01967 for( long ipLo=0; ipLo < ipHi; ipLo++ ) 01968 { 01969 if( Transitions[ipISO][nelem][ipHi][ipLo].Emis->Aul <= iso.SmallA ) 01970 continue; 01971 01972 StatesElem[ipISO][nelem][ipHi].lifetime += Transitions[ipISO][nelem][ipHi][ipLo].Emis->Aul; 01973 } 01974 01975 /* sum of A's was just stuffed, now invert for lifetime. */ 01976 StatesElem[ipISO][nelem][ipHi].lifetime = 1./StatesElem[ipISO][nelem][ipHi].lifetime; 01977 01978 for( long ipLo=0; ipLo < ipHi; ipLo++ ) 01979 { 01980 if( Transitions[ipISO][nelem][ipHi][ipLo].EnergyWN <= 0. ) 01981 continue; 01982 01983 if( Transitions[ipISO][nelem][ipHi][ipLo].Emis->Aul <= iso.SmallA ) 01984 continue; 01985 01986 Transitions[ipISO][nelem][ipHi][ipLo].Emis->dampXvel = (realnum)( 01987 (1.f/StatesElem[ipISO][nelem][ipHi].lifetime)/ 01988 PI4/Transitions[ipISO][nelem][ipHi][ipLo].EnergyWN); 01989 01990 ASSERT(Transitions[ipISO][nelem][ipHi][ipLo].Emis->dampXvel> 0.); 01991 } 01992 } 01993 01994 return; 01995 } 01996 01997 #if 0 01998 STATIC void Prt_AGN_table( void ) 01999 { 02000 /* the designation of the levels, chLevel[n][string] */ 02001 multi_arr<char,2> chLevel(max_num_levels,10); 02002 02003 /* create spectroscopic designation of labels */ 02004 for( long ipLo=0; ipLo < iso.numLevels_max[ipISO][ipISO]-iso.nCollapsed_max[ipISO][ipISO]; ++ipLo ) 02005 { 02006 long nelem = ipISO; 02007 sprintf( &chLevel.front(ipLo) , "%li %li%c", N_(ipLo), S_(ipLo), chL[MIN2(20,L_(ipLo))] ); 02008 } 02009 02010 /* option to print cs data for AGN */ 02011 /* create spectroscopic designation of labels */ 02012 { 02013 /* option to print particulars of some line when called */ 02014 enum {AGN=false}; 02015 if( AGN ) 02016 { 02017 # define NTEMP 6 02018 double te[NTEMP]={6000.,8000.,10000.,15000.,20000.,25000. }; 02019 double telog[NTEMP] , 02020 cs , 02021 ratecoef; 02022 long nelem = ipHELIUM; 02023 fprintf( ioQQQ,"trans"); 02024 for( long i=0; i < NTEMP; ++i ) 02025 { 02026 telog[i] = log10( te[i] ); 02027 fprintf( ioQQQ,"\t%.3e",te[i]); 02028 } 02029 for( long i=0; i < NTEMP; ++i ) 02030 { 02031 fprintf( ioQQQ,"\t%.3e",te[i]); 02032 } 02033 fprintf(ioQQQ,"\n"); 02034 02035 for( long ipHi=ipHe2s3S; ipHi< iso.numLevels_max[ipHE_LIKE][ipHELIUM]; ++ipHi ) 02036 { 02037 for( long ipLo=ipHe1s1S; ipLo < ipHi; ++ipLo ) 02038 { 02039 02040 /* deltaN = 0 transitions may be wrong because 02041 * COLL_CONST below is only correct for electron colliders */ 02042 if( N_(ipHi) == N_(ipLo) ) 02043 continue; 02044 02045 /* print the designations of the lower and upper levels */ 02046 fprintf( ioQQQ,"%s - %s", 02047 &chLevel.front(ipLo) , &chLevel.front(ipHi) ); 02048 02049 /* print the interpolated collision strengths */ 02050 for( long i=0; i < NTEMP; ++i ) 02051 { 02052 phycon.alogte = telog[i]; 02053 /* print cs */ 02054 cs = HeCSInterp( nelem , ipHi , ipLo, ipELECTRON ); 02055 fprintf(ioQQQ,"\t%.2e", cs ); 02056 } 02057 02058 /* print the rate coefficients */ 02059 for( long i=0; i < NTEMP; ++i ) 02060 { 02061 phycon.alogte = telog[i]; 02062 phycon.te = pow(10.,telog[i] ); 02063 tfidle(false); 02064 cs = HeCSInterp( nelem , ipHi , ipLo, ipELECTRON ); 02065 /* collisional deexcitation rate */ 02066 ratecoef = cs/sqrt(phycon.te)*COLL_CONST/StatesElem[ipHE_LIKE][nelem][ipLo].g * 02067 sexp( Transitions[ipHE_LIKE][nelem][ipHi][ipLo].EnergyK / phycon.te ); 02068 fprintf(ioQQQ,"\t%.2e", ratecoef ); 02069 } 02070 fprintf(ioQQQ,"\n"); 02071 } 02072 } 02073 cdEXIT(EXIT_FAILURE); 02074 } 02075 } 02076 02077 return; 02078 } 02079 #endif