00001
00002
00003
00004
00005
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
00025 STATIC void iso_zero(void);
00026
00027
00028 STATIC void iso_allocate(void);
00029
00030
00031 STATIC void iso_assign_quantum_numbers(void);
00032
00033 STATIC void FillExtraLymanLine( transition *t, long ipISO, long nelem, long nHi );
00034
00035 STATIC void iso_satellite( void );
00036
00037 char chL[21]={'S','P','D','F','G','H','I','K','L','M','N','O','Q','R','T','U','V','W','X','Y','Z'};
00038
00039 void iso_create(void)
00040 {
00041 long int ipHi,
00042 ipLo,
00043 ipISO,
00044 nelem;
00045
00046 static int nCalled = 0;
00047
00048 double HIonPoten;
00049
00050 DEBUG_ENTRY( "iso_create()" );
00051
00052
00053 if( nCalled > 0 )
00054 {
00055 iso_zero();
00056 return;
00057 }
00058
00059
00060 ++nCalled;
00061
00062
00063 iso.stat_ion[ipH_LIKE] = 1.f;
00064 iso.stat_ion[ipHE_LIKE] = 2.f;
00065
00066
00067
00068 iso_allocate();
00069
00070
00071 iso_assign_quantum_numbers();
00072
00073
00074 TauDummy.Junk();
00075 TauDummy.Hi = AddState2Stack();
00076 TauDummy.Lo = AddState2Stack();
00077 TauDummy.Emis = AddLine2Stack( true );
00078
00079
00080
00081
00082 for( ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
00083 {
00084
00085 for( nelem=ipISO; nelem < LIMELM; nelem++ )
00086 {
00087
00088 if( nelem < 2 || dense.lgElmtOn[nelem] )
00089 {
00090
00091
00092
00093 HIonPoten = t_ADfA::Inst().ph1(0,0,nelem,0)/EVRYD* 0.9998787;
00094 ASSERT(HIonPoten > 0.);
00095
00096
00097 for( ipHi=0; ipHi < iso.numLevels_max[ipISO][nelem]; ipHi++ )
00098 {
00099 double EnergyWN, EnergyRyd;
00100
00101 if( ipISO == ipH_LIKE )
00102 {
00103 EnergyRyd = HIonPoten/POW2((double)N_(ipHi));
00104 }
00105 else if( ipISO == ipHE_LIKE )
00106 {
00107 EnergyRyd = helike_energy( nelem, ipHi ) * WAVNRYD;
00108 }
00109 else
00110 {
00111
00112 TotalInsanity();
00113 }
00114
00115
00116 ASSERT(EnergyRyd >= 0.);
00117
00118 iso.xIsoLevNIonRyd[ipISO][nelem][ipHi] = EnergyRyd;
00119
00120
00121 for( ipLo=0; ipLo < ipHi; ipLo++ )
00122 {
00123 EnergyWN = RYD_INF * (iso.xIsoLevNIonRyd[ipISO][nelem][ipLo] -
00124 iso.xIsoLevNIonRyd[ipISO][nelem][ipHi]);
00125
00126
00127
00128 if( EnergyWN==0 && ipISO==ipHE_LIKE )
00129 EnergyWN = 0.0001;
00130
00131 if( EnergyWN < 0. )
00132 EnergyWN = -1.0 * EnergyWN;
00133
00134
00135 Transitions[ipISO][nelem][ipHi][ipLo].EnergyWN = (realnum)EnergyWN;
00136 Transitions[ipISO][nelem][ipHi][ipLo].EnergyErg = (realnum)(EnergyWN*WAVNRYD*EN1RYD);
00137 Transitions[ipISO][nelem][ipHi][ipLo].EnergyK = (realnum)(EnergyWN*WAVNRYD*TE1RYD );
00138
00139 ASSERT(Transitions[ipISO][nelem][ipHi][ipLo].EnergyWN >= 0.);
00140 ASSERT(Transitions[ipISO][nelem][ipHi][ipLo].EnergyErg >= 0.);
00141 ASSERT(Transitions[ipISO][nelem][ipHi][ipLo].EnergyK >= 0.);
00142
00144 if( N_(ipLo)==N_(ipHi) && ipISO==ipH_LIKE )
00145 {
00146 Transitions[ipISO][nelem][ipHi][ipLo].WLAng = 0.;
00147 }
00148 else
00149 {
00150
00151 Transitions[ipISO][nelem][ipHi][ipLo].WLAng =
00152 (realnum)(1.0e8/
00153 Transitions[ipISO][nelem][ipHi][ipLo].EnergyWN/
00154 RefIndex( Transitions[ipISO][nelem][ipHi][ipLo].EnergyWN));
00155 ASSERT(Transitions[ipISO][nelem][ipHi][ipLo].WLAng > 0.);
00156 }
00157
00158 }
00159 }
00160
00161
00162 for( ipHi=2; ipHi < iso.nLyman_malloc[ipISO]; ipHi++ )
00163 {
00164 FillExtraLymanLine( &ExtraLymanLines[ipISO][nelem][ipHi], ipISO, nelem, ipHi );
00165 }
00166 }
00167 }
00168 }
00169
00170
00171
00172
00173
00174 iso_recomb_malloc();
00175 iso_recomb_setup( ipH_LIKE );
00176 iso_recomb_setup( ipHE_LIKE );
00177 iso_recomb_auxiliary_free();
00178
00179
00180 HeCollidSetup();
00181
00182
00183
00184
00185 for( ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
00186 {
00187 if( ipISO == ipH_LIKE )
00188 {
00189
00190 }
00191 else if( ipISO == ipHE_LIKE )
00192 {
00193
00194 HelikeTransProbSetup();
00195 }
00196 else
00197 {
00198 TotalInsanity();
00199 }
00200
00201 for( nelem=ipISO; nelem < LIMELM; nelem++ )
00202 {
00203
00204 if( nelem < 2 || dense.lgElmtOn[nelem] )
00205 {
00206 for( ipLo=ipH1s; ipLo < (iso.numLevels_max[ipISO][nelem] - 1); ipLo++ )
00207 {
00208 for( ipHi=ipLo + 1; ipHi < iso.numLevels_max[ipISO][nelem]; ipHi++ )
00209 {
00210 realnum Aul;
00211
00212
00213 if( ipISO == ipH_LIKE )
00214 {
00215 Aul = hydro_transprob( nelem, ipHi, ipLo );
00216 }
00217 else if( ipISO == ipHE_LIKE )
00218 {
00219 Aul = helike_transprob(nelem, ipHi, ipLo);
00220 }
00221 else
00222 {
00223 TotalInsanity();
00224 }
00225
00226 if( Aul <= iso.SmallA )
00227 Transitions[ipISO][nelem][ipHi][ipLo].Emis = AddLine2Stack( false );
00228 else
00229 Transitions[ipISO][nelem][ipHi][ipLo].Emis = AddLine2Stack( true );
00230
00231 Transitions[ipISO][nelem][ipHi][ipLo].Emis->Aul = Aul;
00232
00233 ASSERT(Transitions[ipISO][nelem][ipHi][ipLo].Emis->Aul > 0.);
00234
00235 if( ipLo == 0 && ipHi == iso.nLyaLevel[ipISO] )
00236 {
00237
00238 Transitions[ipISO][nelem][ipHi][ipLo].Emis->iRedisFun = iso.ipLyaRedist[ipISO];
00239 }
00240 else if( ipLo == 0 )
00241 {
00242
00243
00244 Transitions[ipISO][nelem][ipHi][ipLo].Emis->iRedisFun = iso.ipResoRedist[ipISO];
00245 }
00246 else
00247 {
00248
00249
00250 Transitions[ipISO][nelem][ipHi][ipLo].Emis->iRedisFun = iso.ipSubRedist[ipISO];
00251 }
00252
00253 if( Transitions[ipISO][nelem][ipHi][ipLo].Emis->Aul <= iso.SmallA ||
00254 Transitions[ipISO][nelem][ipHi][ipLo].EnergyWN <= 0.)
00255 {
00256 Transitions[ipISO][nelem][ipHi][ipLo].Emis->gf = 0.;
00257 Transitions[ipISO][nelem][ipHi][ipLo].Emis->opacity = 0.;
00258 }
00259 else
00260 {
00261 Transitions[ipISO][nelem][ipHi][ipLo].Emis->gf =
00262 (realnum)(GetGF(Transitions[ipISO][nelem][ipHi][ipLo].Emis->Aul,
00263 Transitions[ipISO][nelem][ipHi][ipLo].EnergyWN,
00264 Transitions[ipISO][nelem][ipHi][ipLo].Hi->g));
00265 ASSERT(Transitions[ipISO][nelem][ipHi][ipLo].Emis->gf > 0.);
00266
00267
00268 Transitions[ipISO][nelem][ipHi][ipLo].Emis->opacity =
00269 (realnum)(abscf(Transitions[ipISO][nelem][ipHi][ipLo].Emis->gf,
00270 Transitions[ipISO][nelem][ipHi][ipLo].EnergyWN,
00271 Transitions[ipISO][nelem][ipHi][ipLo].Lo->g));
00272 ASSERT(Transitions[ipISO][nelem][ipHi][ipLo].Emis->opacity > 0.);
00273 }
00274 }
00275 }
00276 }
00277 }
00278 }
00279
00280
00281
00282
00283 if( iso.lgFSM[ipHE_LIKE] )
00284 {
00285
00286 for( nelem=ipISO; nelem < LIMELM; nelem++ )
00287 {
00288
00289 if( nelem < 2 || dense.lgElmtOn[nelem] )
00290 {
00291 for( ipHi=ipHe2s3S; ipHi<iso.numLevels_max[ipHE_LIKE][nelem]; ipHi++ )
00292 {
00293 for( ipLo=ipHe1s1S; ipLo<ipHi; ipLo++ )
00294 {
00295 DoFSMixing( nelem, ipLo, ipHi );
00296 }
00297 }
00298 }
00299 }
00300 }
00301
00302
00303 Transitions[ipH_LIKE][ipHELIUM][ipH3s][ipH2s].WLAng = 1640.f;
00304 Transitions[ipH_LIKE][ipHELIUM][ipH3s][ipH2p].WLAng = 1640.f;
00305 if( iso.n_HighestResolved_max[ipH_LIKE][ipHELIUM] >=3 )
00306 {
00307 Transitions[ipH_LIKE][ipHELIUM][ipH3p][ipH2s].WLAng = 1640.f;
00308 Transitions[ipH_LIKE][ipHELIUM][ipH3p][ipH2p].WLAng = 1640.f;
00309 Transitions[ipH_LIKE][ipHELIUM][ipH3d][ipH2s].WLAng = 1640.f;
00310 Transitions[ipH_LIKE][ipHELIUM][ipH3d][ipH2p].WLAng = 1640.f;
00311 }
00312
00313
00314
00315
00316 for( ipISO=ipH_LIKE; ipISO<NISO; ipISO++ )
00317 {
00318 for( nelem=ipISO; nelem < LIMELM; nelem++ )
00319 {
00320
00321 if( nelem < 2 || dense.lgElmtOn[nelem] )
00322 {
00323
00324 StatesElemNEW[nelem][nelem-ipISO][0].lifetime = -FLT_MAX;
00325
00326 for( ipHi=1; ipHi < iso.numLevels_max[ipISO][nelem]; ipHi++ )
00327 {
00328 StatesElemNEW[nelem][nelem-ipISO][ipHi].lifetime = SMALLFLOAT;
00329
00330 for( ipLo=0; ipLo < ipHi; ipLo++ )
00331 {
00332 if( Transitions[ipISO][nelem][ipHi][ipLo].Emis->Aul <= iso.SmallA )
00333 continue;
00334
00335 StatesElemNEW[nelem][nelem-ipISO][ipHi].lifetime += Transitions[ipISO][nelem][ipHi][ipLo].Emis->Aul;
00336 }
00337
00338
00339 StatesElemNEW[nelem][nelem-ipISO][ipHi].lifetime = 1./StatesElemNEW[nelem][nelem-ipISO][ipHi].lifetime;
00340
00341 for( ipLo=0; ipLo < ipHi; ipLo++ )
00342 {
00343 if( Transitions[ipISO][nelem][ipHi][ipLo].EnergyWN <= 0. )
00344 continue;
00345
00346 if( Transitions[ipISO][nelem][ipHi][ipLo].Emis->Aul <= iso.SmallA )
00347 continue;
00348
00349 Transitions[ipISO][nelem][ipHi][ipLo].Emis->dampXvel = (realnum)(
00350 (1.f/StatesElemNEW[nelem][nelem-ipISO][ipHi].lifetime)/
00351 PI4/Transitions[ipISO][nelem][ipHi][ipLo].EnergyWN);
00352
00353 ASSERT(Transitions[ipISO][nelem][ipHi][ipLo].Emis->dampXvel> 0.);
00354 }
00355 }
00356 }
00357 }
00358 }
00359
00360
00361
00362
00363 for( ipISO=ipH_LIKE; ipISO<NISO; ipISO++ )
00364 {
00365
00366 for( nelem=ipISO; nelem < LIMELM; nelem++ )
00367 {
00368
00369 if( nelem < 2 || dense.lgElmtOn[nelem] )
00370 {
00371
00372 ASSERT( ipISO <= ipHE_LIKE );
00373
00374
00375 Transitions[ipISO][nelem][1+ipISO][0].Emis->TauTot = 0.;
00376
00377
00378
00379 Transitions[ipISO][nelem][1+ipISO][0].Emis->opacity /= 1e4f;
00380
00381
00382 Transitions[ipISO][nelem][1+ipISO][0].WLAng = 0.;
00383 }
00384 }
00385 }
00386
00387
00388
00389 iso_zero();
00390
00391
00392 for( ipISO=ipH_LIKE; ipISO<NISO; ipISO++ )
00393 {
00394 for( nelem = ipISO; nelem < LIMELM; nelem++ )
00395 {
00396
00397 if( nelem == ipISO || dense.lgElmtOn[nelem] )
00398 {
00399
00400 iso_cascade( ipISO, nelem);
00401 }
00402 }
00403 }
00404
00405 iso_satellite();
00406
00407 iso_satellite_update();
00408
00409
00410
00411
00412
00413 for( ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
00414 {
00415 for( nelem=ipISO; nelem < LIMELM; nelem++ )
00416 {
00417 if( nelem < 2 || dense.lgElmtOn[nelem] )
00418 {
00419 for( ipLo=ipH1s; ipLo < (iso.numLevels_max[ipISO][nelem] - 1); ipLo++ )
00420 {
00421 for( ipHi= ipLo + 1; ipHi < iso.numLevels_max[ipISO][nelem]; ipHi++ )
00422 {
00423 long nHi = StatesElemNEW[nelem][nelem-ipISO][ipHi].n;
00424 long nLo = StatesElemNEW[nelem][nelem-ipISO][ipLo].n;
00425
00426 iso.strkar[ipISO][nelem][ipLo][ipHi] = (realnum)pow((realnum)( nLo * nHi ),(realnum)1.2f);
00427 iso.pestrk[ipISO][nelem][ipLo][ipHi] = 0.;
00428 }
00429 }
00430 }
00431 }
00432 }
00433
00434 return;
00435 }
00436
00437
00438 STATIC void iso_zero(void)
00439 {
00440 long int i,
00441 ipHi,
00442 ipISO,
00443 ipLo,
00444 nelem;
00445
00446 DEBUG_ENTRY( "iso_zero()" );
00447
00448 fixit();
00449
00450
00451 hydro.HLineWidth = 0.;
00452
00453
00454
00455
00456 for( ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
00457 {
00458 for( nelem=ipISO; nelem < LIMELM; nelem++ )
00459 {
00460 if( nelem < 2 || dense.lgElmtOn[nelem] )
00461 {
00462 for( ipHi=0; ipHi < iso.numLevels_max[ipISO][nelem]; ipHi++ )
00463 {
00464 StatesElemNEW[nelem][nelem-ipISO][ipHi].Pop = 0.;
00465
00466 iso.PopLTE[ipISO][nelem][ipHi] = 0.;
00467 iso.ColIoniz[ipISO][nelem][ipHi] = 0.;
00468 iso.gamnc[ipISO][nelem][ipHi] = -DBL_MAX;
00469 iso.RecomInducRate[ipISO][nelem][ipHi] = -DBL_MAX;
00470 iso.DepartCoef[ipISO][nelem][ipHi] = -DBL_MAX;
00471 iso.RateLevel2Cont[ipISO][nelem][ipHi] = 0.;
00472 iso.RateCont2Level[ipISO][nelem][ipHi] = 0.;
00473 iso.ConOpacRatio[ipISO][nelem][ipHi] = 1.;
00474 iso.RadRecCon[ipISO][nelem][ipHi] = 0.;
00475 iso.RadRecomb[ipISO][nelem][ipHi][ipRecRad] = 0.;
00476 iso.RadRecomb[ipISO][nelem][ipHi][ipRecNetEsc] = 1.;
00477 iso.RadRecomb[ipISO][nelem][ipHi][ipRecEsc] = 1.;
00478 iso.DielecRecomb[ipISO][nelem][ipHi] = 0.;
00479 }
00480 }
00481
00482 if( nelem < 2 )
00483 {
00484 iso_collapsed_bnl_set( ipISO, nelem );
00485
00486 iso_collapsed_Aul_update( ipISO, nelem );
00487 iso_collapsed_lifetimes_update( ipISO, nelem );
00488 }
00489 }
00490 }
00491
00492
00493
00494 iso.ConOpacRatio[ipH_LIKE][ipHYDROGEN][0] = 1e-5;
00495 iso.ConOpacRatio[ipH_LIKE][ipHELIUM][0] = 1e-5;
00496 iso.ConOpacRatio[ipHE_LIKE][ipHELIUM][0] = 1e-5;
00497
00498
00499
00500 for( ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
00501 {
00502 for( nelem=ipISO; nelem < LIMELM; nelem++ )
00503 {
00504
00505 iso.numLevels_local[ipISO][nelem] = iso.numLevels_max[ipISO][nelem];
00506
00507
00508 if( nelem < 2 || dense.lgElmtOn[nelem] )
00509 {
00510 SatelliteLines[ipISO][nelem][0].Zero();
00511 for( ipHi=1; ipHi < iso.numLevels_max[ipISO][nelem]; ipHi++ )
00512 {
00513 SatelliteLines[ipISO][nelem][ipHi].Zero();
00514
00515 for( ipLo=0; ipLo < ipHi; ipLo++ )
00516 {
00517 Transitions[ipISO][nelem][ipHi][ipLo].Zero();
00518 }
00519 }
00520
00521 for( ipHi=2; ipHi < iso.nLyman[ipISO]; ipHi++ )
00522 {
00523 ExtraLymanLines[ipISO][nelem][ipHi].Zero();
00524 }
00525 }
00526 }
00527 }
00528
00529
00530 for( i=0; i <= nLevel1; i++ )
00531 {
00532 TauLines[i].Zero();
00533 }
00534
00535
00536
00537
00538 TauDummy.Zero();
00539
00540 for( i=0; i < nUTA; i++ )
00541 {
00542
00543 double hsave = UTALines[i].Coll.heat;
00544 UTALines[i].Zero();
00545 UTALines[i].Coll.heat = hsave;
00546 }
00547
00548 for( i=0; i < nWindLine; i++ )
00549 {
00550 TauLine2[i].Zero();
00551 }
00552
00553 for( i=0; i < nHFLines; i++ )
00554 {
00555 HFLines[i].Zero();
00556 }
00557
00558
00559 for( i=0; i <linesAdded2; i++)
00560 {
00561 dBaseLines[i].tran->Zero();
00562 }
00563
00564 return;
00565 }
00566
00567 STATIC void iso_allocate(void)
00568 {
00569
00570 DEBUG_ENTRY( "iso_allocate()" );
00571
00572 iso.strkar.reserve( NISO );
00573 iso.ipOpac.reserve( NISO );
00574 iso.RadRecomb.reserve( NISO );
00575 iso.RadRecCon.reserve( NISO );
00576 iso.DielecRecombVsTemp.reserve( NISO );
00577 iso.Boltzmann.reserve( NISO );
00578 iso.QuantumNumbers2Index.reserve( NISO );
00579 iso.BranchRatio.reserve( NISO );
00580 iso.CascadeProb.reserve( NISO );
00581 iso.CachedAs.reserve( NISO );
00582
00583
00584 for( long ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
00585 {
00586 iso.strkar.reserve( ipISO, LIMELM );
00587 iso.ipOpac.reserve( ipISO, LIMELM );
00588 iso.RadRecomb.reserve( ipISO, LIMELM );
00589 iso.RadRecCon.reserve( ipISO, LIMELM );
00590 iso.DielecRecombVsTemp.reserve( ipISO, LIMELM );
00591 iso.Boltzmann.reserve( ipISO, LIMELM );
00592 iso.QuantumNumbers2Index.reserve( ipISO, LIMELM );
00593 iso.BranchRatio.reserve( ipISO, LIMELM );
00594 iso.CascadeProb.reserve( ipISO, LIMELM );
00595 iso.CachedAs.reserve( ipISO, LIMELM );
00596
00597 for( long nelem=ipISO; nelem < LIMELM; ++nelem )
00598 {
00599
00600 if( nelem < 2 || dense.lgElmtOn[nelem] )
00601 {
00602 iso.numLevels_malloc[ipISO][nelem] = iso.numLevels_max[ipISO][nelem];
00603
00604
00605 ASSERT( iso.numLevels_max[ipISO][nelem] > 0 );
00606
00607 iso.nLyman_malloc[ipISO] = iso.nLyman[ipISO];
00608
00609 iso.strkar.reserve( ipISO, nelem, iso.numLevels_max[ipISO][nelem] );
00610 iso.ipOpac.reserve( ipISO, nelem, iso.numLevels_max[ipISO][nelem] );
00611 iso.RadRecomb.reserve( ipISO, nelem, iso.numLevels_max[ipISO][nelem] );
00612 iso.RadRecCon.reserve( ipISO, nelem, iso.numLevels_max[ipISO][nelem] );
00613 iso.DielecRecombVsTemp.reserve( ipISO, nelem, iso.numLevels_max[ipISO][nelem] );
00614 iso.Boltzmann.reserve( ipISO, nelem, iso.numLevels_max[ipISO][nelem] );
00615
00616 iso.CachedAs.reserve( ipISO, nelem, MAX2(1, iso.nCollapsed_max[ipISO][nelem]) );
00617
00618 for( long i = 0; i < iso.nCollapsed_max[ipISO][nelem]; ++i )
00619 {
00620 iso.CachedAs.reserve( ipISO, nelem, i, iso.numLevels_max[ipISO][nelem] - iso.nCollapsed_max[ipISO][nelem] );
00621 for( long i1 = 0; i1 < iso.numLevels_max[ipISO][nelem] - iso.nCollapsed_max[ipISO][nelem]; ++i1 )
00622 {
00623
00624 iso.CachedAs.reserve( ipISO, nelem, i, i1, 2 );
00625 }
00626 }
00627
00628 iso.QuantumNumbers2Index.reserve( ipISO, nelem, iso.n_HighestResolved_max[ipISO][nelem] + iso.nCollapsed_max[ipISO][nelem] + 1 );
00629
00630
00631 for( long i = 1; i <= (iso.n_HighestResolved_max[ipISO][nelem] + iso.nCollapsed_max[ipISO][nelem]); ++i )
00632 {
00633
00634 iso.QuantumNumbers2Index.reserve( ipISO, nelem, i, i );
00635
00636 for( long i1 = 0; i1 < i; ++i1 )
00637 {
00638
00639 ASSERT( NISO == 2 );
00640
00641
00642 iso.QuantumNumbers2Index.reserve( ipISO, nelem, i, i1, 4 );
00643 }
00644 }
00645
00646 for( long n=0; n < iso.numLevels_max[ipISO][nelem]; ++n )
00647 {
00648 iso.RadRecomb.reserve( ipISO, nelem, n, 3 );
00649
00650
00651
00652 if( n > 0 )
00653 iso.Boltzmann.reserve( ipISO, nelem, n, n );
00654
00655
00656 iso.DielecRecombVsTemp.reserve( ipISO, nelem, n, NUM_DR_TEMPS );
00657 }
00658
00659
00660 iso.CascadeProb.reserve( ipISO, nelem, iso.numLevels_max[ipISO][nelem] );
00661 iso.BranchRatio.reserve( ipISO, nelem, iso.numLevels_max[ipISO][nelem]+1 );
00662
00663 for( long i = 1; i < iso.numLevels_max[ipISO][nelem]; i++ )
00664 iso.BranchRatio.reserve( ipISO, nelem, i, i+1 );
00665
00666
00667 for( long i = 0; i < iso.numLevels_max[ipISO][nelem]; ++i )
00668 {
00669 iso.strkar.reserve( ipISO, nelem, i, iso.numLevels_max[ipISO][nelem] );
00670 iso.CascadeProb.reserve( ipISO, nelem, i, i+1 );
00671 }
00672 }
00673 }
00674 }
00675
00676 iso.strkar.alloc();
00677 iso.pestrk.alloc( iso.strkar.clone() );
00678 iso.ipOpac.alloc();
00679 secondaries.Hx12.alloc( iso.ipOpac.clone() );
00680 iso.ipIsoLevNIonCon.alloc( iso.ipOpac.clone() );
00681 iso.xIsoLevNIonRyd.alloc( iso.ipOpac.clone() );
00682 iso.DielecRecomb.alloc( iso.ipOpac.clone() );
00683 iso.DepartCoef.alloc( iso.ipOpac.clone() );
00684 iso.RateLevel2Cont.alloc( iso.ipOpac.clone() );
00685 iso.RateCont2Level.alloc( iso.ipOpac.clone() );
00686 iso.ConOpacRatio.alloc( iso.ipOpac.clone() );
00687 iso.gamnc.alloc( iso.ipOpac.clone() );
00688 iso.RecomInducRate.alloc( iso.ipOpac.clone() );
00689 iso.RecomInducCool_Coef.alloc( iso.ipOpac.clone() );
00690 iso.PhotoHeat.alloc( iso.ipOpac.clone() );
00691 iso.PopLTE.alloc( iso.ipOpac.clone() );
00692 iso.ColIoniz.alloc( iso.ipOpac.clone() );
00693 iso.RadEffec.alloc( iso.ipOpac.clone() );
00694 iso.RadRecCon.alloc();
00695 iso.RadRecomb.alloc();
00696 iso.Boltzmann.alloc();
00697 iso.DielecRecombVsTemp.alloc();
00698 iso.CachedAs.alloc();
00699 iso.QuantumNumbers2Index.alloc();
00700 iso.BranchRatio.alloc();
00701 iso.CascadeProb.alloc();
00702
00703 iso.bnl_effective.alloc( iso.QuantumNumbers2Index.clone() );
00704
00705 iso.DielecRecombVsTemp.zero();
00706 iso.QuantumNumbers2Index.invalidate();
00707 iso.bnl_effective.invalidate();
00708 iso.BranchRatio.invalidate();
00709 iso.CascadeProb.invalidate();
00710
00711 SatelliteLines.reserve( NISO );
00712 Transitions.reserve( NISO );
00713 ExtraLymanLines.reserve( NISO );
00714
00715 for( long ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
00716 {
00717 SatelliteLines.reserve( ipISO, LIMELM );
00718 Transitions.reserve( ipISO, LIMELM );
00719 ExtraLymanLines.reserve( ipISO, LIMELM );
00720
00721 for( long nelem=ipISO; nelem < LIMELM; ++nelem )
00722 {
00723
00724 if( nelem < 2 || dense.lgElmtOn[nelem] )
00725 {
00726 ASSERT( iso.numLevels_max[ipISO][nelem] > 0 );
00727
00728 SatelliteLines.reserve( ipISO, nelem, iso.numLevels_max[ipISO][nelem] );
00729 Transitions.reserve( ipISO, nelem, iso.numLevels_max[ipISO][nelem] );
00730 ExtraLymanLines.reserve( ipISO, nelem, iso.nLyman_malloc[ipISO] );
00731
00732 for( long ipHi=1; ipHi<iso.numLevels_max[ipISO][nelem]; ipHi++ )
00733 Transitions.reserve( ipISO, nelem, ipHi, ipHi );
00734 }
00735 }
00736 }
00737
00738 SatelliteLines.alloc();
00739 Transitions.alloc();
00740 ExtraLymanLines.alloc();
00741
00742 for( long ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
00743 {
00744 for( long nelem=ipISO; nelem < LIMELM; ++nelem )
00745 {
00746
00747 if( nelem < 2 || dense.lgElmtOn[nelem] )
00748 {
00749 for( long ipLo=0; ipLo<iso.numLevels_max[ipISO][nelem]; ipLo++ )
00750 {
00751
00752
00753 SatelliteLines[ipISO][nelem][ipLo].Junk();
00754 SatelliteLines[ipISO][nelem][ipLo].Hi = AddState2Stack();
00755 SatelliteLines[ipISO][nelem][ipLo].Lo = &StatesElemNEW[nelem][nelem-ipISO][ipLo];
00756 SatelliteLines[ipISO][nelem][ipLo].Emis = AddLine2Stack( true );
00757 }
00758
00759 for( long ipHi=1; ipHi<iso.numLevels_max[ipISO][nelem]; ipHi++ )
00760 {
00761 for( long ipLo=0; ipLo < ipHi; ipLo++ )
00762 {
00763
00764 Transitions[ipISO][nelem][ipHi][ipLo].Junk();
00765 Transitions[ipISO][nelem][ipHi][ipLo].Hi = &StatesElemNEW[nelem][nelem-ipISO][ipHi];
00766 Transitions[ipISO][nelem][ipHi][ipLo].Lo = &StatesElemNEW[nelem][nelem-ipISO][ipLo];
00767 }
00768 }
00769
00770
00771 for( long ipHi=2; ipHi < iso.nLyman_malloc[ipISO]; ipHi++ )
00772 {
00773 ExtraLymanLines[ipISO][nelem][ipHi].Junk();
00774 ExtraLymanLines[ipISO][nelem][ipHi].Hi = AddState2Stack();
00775
00776 ExtraLymanLines[ipISO][nelem][ipHi].Lo = &StatesElemNEW[nelem][nelem-ipISO][0];
00777 ExtraLymanLines[ipISO][nelem][ipHi].Emis = AddLine2Stack( true );
00778 }
00779 }
00780 }
00781 }
00782
00783
00784 for( long ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
00785 {
00786 static bool lgNotSetYet = true;
00787
00788 if( iso.lgRandErrGen[ipISO] && lgNotSetYet )
00789 {
00790 iso.Error.reserve( NISO );
00791 iso.SigmaCascadeProb.reserve( NISO );
00792 lgNotSetYet = false;
00793 }
00794
00795 if( iso.lgRandErrGen[ipISO] )
00796 {
00797 ASSERT( !lgNotSetYet );
00798 iso.Error.reserve( ipISO, LIMELM );
00799 iso.SigmaCascadeProb.reserve( ipISO, LIMELM );
00800
00801 for( long nelem=ipISO; nelem < LIMELM; ++nelem )
00802 {
00803 if( nelem == ipHELIUM || dense.lgElmtOn[nelem] )
00804 {
00805
00806 iso.Error.reserve( ipISO, nelem, iso.numLevels_max[ipISO][nelem]+1 );
00807
00808 for( long i = 1; i< iso.numLevels_max[ipISO][nelem] + 1; i++ )
00809 {
00810 iso.Error.reserve( ipISO, nelem, i, i+1 );
00811
00812 for( long j = 0; j<i+1; j++ )
00813 iso.Error.reserve( ipISO, nelem, i, j, 3 );
00814 }
00815
00816
00817
00818 iso.SigmaCascadeProb.reserve( ipISO, nelem, iso.numLevels_max[ipISO][nelem] );
00819 for( long i=0; i < iso.numLevels_max[ipISO][nelem]; i++ )
00820
00821 iso.SigmaCascadeProb.reserve( ipISO, nelem, i, i+1 );
00822 }
00823 }
00824 iso.Error.alloc();
00825 iso.ErrorFactor.alloc( iso.Error.clone() );
00826 iso.SigmaAtot.alloc( iso.ipOpac.clone() );
00827 iso.SigmaRadEffec.alloc( iso.RadEffec.clone() );
00828 iso.SigmaCascadeProb.alloc();
00829
00830 iso.Error.invalidate();
00831 iso.ErrorFactor.invalidate();
00832 }
00833 }
00834
00835 return;
00836 }
00837
00838 STATIC void iso_assign_quantum_numbers(void)
00839 {
00840 long int
00841 ipISO,
00842 ipLo,
00843 level,
00844 nelem,
00845 i,
00846 in,
00847 il,
00848 is,
00849 ij;
00850
00851 DEBUG_ENTRY( "iso_assign_quantum_numbers()" );
00852
00853 ipISO = ipH_LIKE;
00854 for( nelem=ipISO; nelem < LIMELM; nelem++ )
00855 {
00856
00857 if( nelem == ipHELIUM || dense.lgElmtOn[nelem] )
00858 {
00859 i = 0;
00860
00861
00862 is = 2;
00863
00864
00865
00866 for( in = 1L; in <= iso.n_HighestResolved_max[ipISO][nelem]; ++in )
00867 {
00868 for( il = 0L; il < in; ++il )
00869 {
00870 StatesElemNEW[nelem][nelem-ipISO][i].n = in;
00871 StatesElemNEW[nelem][nelem-ipISO][i].S = is;
00872 StatesElemNEW[nelem][nelem-ipISO][i].l = il;
00873 StatesElemNEW[nelem][nelem-ipISO][i].j = -1;
00874 iso.QuantumNumbers2Index[ipISO][nelem][in][il][is] = i;
00875 ++i;
00876 }
00877 }
00878
00879 in = iso.n_HighestResolved_max[ipISO][nelem] + 1;
00880 for( level = i; level< iso.numLevels_max[ipISO][nelem]; ++level)
00881 {
00882 StatesElemNEW[nelem][nelem-ipISO][level].n = in;
00883 StatesElemNEW[nelem][nelem-ipISO][level].S = -LONG_MAX;
00884 StatesElemNEW[nelem][nelem-ipISO][level].l = -LONG_MAX;
00885 StatesElemNEW[nelem][nelem-ipISO][level].j = -1;
00886
00887 for( il = 0; il < in; ++il )
00888 {
00889 iso.QuantumNumbers2Index[ipISO][nelem][in][il][is] = level;
00890 }
00891 ++in;
00892 }
00893 --in;
00894
00895
00896 ASSERT( i <= iso.numLevels_max[ipISO][nelem] );
00897
00898
00899 ASSERT( (in > 0) && (in < (iso.n_HighestResolved_max[ipISO][nelem] + iso.nCollapsed_max[ipISO][nelem] + 1) ) );
00900
00901
00902 for( in = 2L; in <= iso.n_HighestResolved_max[ipISO][nelem] + iso.nCollapsed_max[ipISO][nelem]; ++in )
00903 {
00904 for( il = 0L; il < in; ++il )
00905 {
00906 ASSERT( StatesElemNEW[nelem][nelem-ipISO][ iso.QuantumNumbers2Index[ipISO][nelem][in][il][is] ].n == in );
00907 if( in <= iso.n_HighestResolved_max[ipISO][nelem] )
00908 {
00909
00910
00911 ASSERT( StatesElemNEW[nelem][nelem-ipISO][ iso.QuantumNumbers2Index[ipISO][nelem][in][il][is] ].l == il );
00912 ASSERT( StatesElemNEW[nelem][nelem-ipISO][ iso.QuantumNumbers2Index[ipISO][nelem][in][il][is] ].S == is );
00913 }
00914 }
00915 }
00916 }
00917 }
00918
00919
00920 ipISO = ipHE_LIKE;
00921 for( nelem=ipHELIUM; nelem < LIMELM; nelem++ )
00922 {
00923
00924 if( nelem == ipHELIUM || dense.lgElmtOn[nelem] )
00925 {
00926 i = 0;
00927
00928
00929
00930 for( in = 1L; in <= iso.n_HighestResolved_max[ipISO][nelem]; ++in )
00931 {
00932 for( il = 0L; il < in; ++il )
00933 {
00934 for( is = 3L; is >= 1L; is -= 2 )
00935 {
00936
00937
00938
00939 if( (il == 1L) && (is == 1L) )
00940 continue;
00941
00942 if( (in == 1L) && (is == 3L) )
00943 continue;
00944
00945
00946 if( is == 1 )
00947 {
00948 StatesElemNEW[nelem][nelem-ipISO][i].n = in;
00949 StatesElemNEW[nelem][nelem-ipISO][i].S = is;
00950 StatesElemNEW[nelem][nelem-ipISO][i].l = il;
00951
00952 StatesElemNEW[nelem][nelem-ipISO][i].j = il;
00953 iso.QuantumNumbers2Index[ipISO][nelem][in][il][is] = i;
00954 ++i;
00955 }
00956
00957 else if( (in == 2) && (il == 1) && (is == 3) )
00958 {
00959 ij = 0;
00960 do
00961 {
00962 StatesElemNEW[nelem][nelem-ipISO][i].n = in;
00963 StatesElemNEW[nelem][nelem-ipISO][i].S = is;
00964 StatesElemNEW[nelem][nelem-ipISO][i].l = il;
00965 StatesElemNEW[nelem][nelem-ipISO][i].j = ij;
00966 iso.QuantumNumbers2Index[ipISO][nelem][in][il][is] = i;
00967 ++i;
00968 ++ij;
00969
00970 } while ( ij < 3 );
00971 }
00972 else
00973 {
00974 StatesElemNEW[nelem][nelem-ipISO][i].n = in;
00975 StatesElemNEW[nelem][nelem-ipISO][i].S = is;
00976 StatesElemNEW[nelem][nelem-ipISO][i].l = il;
00977 StatesElemNEW[nelem][nelem-ipISO][i].j = -1L;
00978 iso.QuantumNumbers2Index[ipISO][nelem][in][il][is] = i;
00979 ++i;
00980 }
00981 }
00982 }
00983
00984 if( in > 1L )
00985 {
00986 StatesElemNEW[nelem][nelem-ipISO][i].n = in;
00987 StatesElemNEW[nelem][nelem-ipISO][i].S = 1L;
00988 StatesElemNEW[nelem][nelem-ipISO][i].l = 1L;
00989 StatesElemNEW[nelem][nelem-ipISO][i].j = 1L;
00990 iso.QuantumNumbers2Index[ipISO][nelem][in][1][1] = i;
00991 ++i;
00992 }
00993 }
00994
00995 in = iso.n_HighestResolved_max[ipISO][nelem] + 1;
00996 for( level = i; level< iso.numLevels_max[ipISO][nelem]; ++level)
00997 {
00998 StatesElemNEW[nelem][nelem-ipISO][level].n = in;
00999 StatesElemNEW[nelem][nelem-ipISO][level].S = -LONG_MAX;
01000 StatesElemNEW[nelem][nelem-ipISO][level].l = -LONG_MAX;
01001 StatesElemNEW[nelem][nelem-ipISO][level].j = -1;
01002
01003 for( il = 0; il < in; ++il )
01004 {
01005 for( is = 1; is <= 3; is += 2 )
01006 {
01007 iso.QuantumNumbers2Index[ipISO][nelem][in][il][is] = level;
01008 }
01009 }
01010 ++in;
01011 }
01012 --in;
01013
01014
01015 ASSERT( i <= iso.numLevels_max[ipISO][nelem] );
01016
01017
01018 ASSERT( (in > 0) && (in < (iso.n_HighestResolved_max[ipISO][nelem] + iso.nCollapsed_max[ipISO][nelem] + 1) ) );
01019
01020
01021 for( in = 2L; in <= iso.n_HighestResolved_max[ipISO][nelem] + iso.nCollapsed_max[ipISO][nelem]; ++in )
01022 {
01023 for( il = 0L; il < in; ++il )
01024 {
01025 for( is = 3L; is >= 1; is -= 2 )
01026 {
01027
01028 if( (in == 1L) && (is == 3L) )
01029 continue;
01030
01031 ASSERT( StatesElemNEW[nelem][nelem-ipISO][ iso.QuantumNumbers2Index[ipISO][nelem][in][il][is] ].n == in );
01032 if( in <= iso.n_HighestResolved_max[ipISO][nelem] )
01033 {
01034
01035
01036 ASSERT( StatesElemNEW[nelem][nelem-ipISO][ iso.QuantumNumbers2Index[ipISO][nelem][in][il][is] ].l == il );
01037 ASSERT( StatesElemNEW[nelem][nelem-ipISO][ iso.QuantumNumbers2Index[ipISO][nelem][in][il][is] ].S == is );
01038 }
01039 }
01040 }
01041 }
01042 }
01043 }
01044
01045 for( ipISO=ipH_LIKE; ipISO<NISO; ipISO++ )
01046 {
01047 for( nelem=ipISO; nelem < LIMELM; nelem++ )
01048 {
01049
01050 if( nelem < 2 || dense.lgElmtOn[nelem] )
01051 {
01052 for( ipLo=ipH1s; ipLo < iso.numLevels_max[ipISO][nelem]; ipLo++ )
01053 {
01054 StatesElemNEW[nelem][nelem-ipISO][ipLo].nelem = (int)(nelem+1);
01055 StatesElemNEW[nelem][nelem-ipISO][ipLo].IonStg = (int)(nelem+1-ipISO);
01056
01057 if( StatesElemNEW[nelem][nelem-ipISO][ipLo].j >= 0 )
01058 {
01059 StatesElemNEW[nelem][nelem-ipISO][ipLo].g = 2.f*StatesElemNEW[nelem][nelem-ipISO][ipLo].j+1.f;
01060 }
01061 else if( StatesElemNEW[nelem][nelem-ipISO][ipLo].l >= 0 )
01062 {
01063 StatesElemNEW[nelem][nelem-ipISO][ipLo].g = (2.f*StatesElemNEW[nelem][nelem-ipISO][ipLo].l+1.f) *
01064 StatesElemNEW[nelem][nelem-ipISO][ipLo].S;
01065 }
01066 else
01067 {
01068 if( ipISO == ipH_LIKE )
01069 StatesElemNEW[nelem][nelem-ipISO][ipLo].g = 2.f*(realnum)POW2( StatesElemNEW[nelem][nelem-ipISO][ipLo].n );
01070 else if( ipISO == ipHE_LIKE )
01071 StatesElemNEW[nelem][nelem-ipISO][ipLo].g = 4.f*(realnum)POW2( StatesElemNEW[nelem][nelem-ipISO][ipLo].n );
01072 else
01073 {
01074
01075 TotalInsanity();
01076 }
01077 }
01078 char chConfiguration[11] = " ";
01079 long nCharactersWritten = 0;
01080
01081 ASSERT( StatesElemNEW[nelem][nelem-ipISO][ipLo].n < 1000 );
01082
01083
01084 if( StatesElemNEW[nelem][nelem-ipISO][ipLo].n > iso.n_HighestResolved_max[ipISO][nelem] )
01085 {
01086 nCharactersWritten = sprintf( chConfiguration, "n=%3li",
01087 StatesElemNEW[nelem][nelem-ipISO][ipLo].n );
01088 }
01089 else if( StatesElemNEW[nelem][nelem-ipISO][ipLo].j > 0 )
01090 {
01091 nCharactersWritten = sprintf( chConfiguration, "%3li^%li%c_%li",
01092 StatesElemNEW[nelem][nelem-ipISO][ipLo].n,
01093 StatesElemNEW[nelem][nelem-ipISO][ipLo].S,
01094 chL[ MIN2( 20, StatesElemNEW[nelem][nelem-ipISO][ipLo].l ) ],
01095 StatesElemNEW[nelem][nelem-ipISO][ipLo].j );
01096 }
01097 else
01098 {
01099 nCharactersWritten = sprintf( chConfiguration, "%3li^%li%c",
01100 StatesElemNEW[nelem][nelem-ipISO][ipLo].n,
01101 StatesElemNEW[nelem][nelem-ipISO][ipLo].S,
01102 chL[ MIN2( 20, StatesElemNEW[nelem][nelem-ipISO][ipLo].l) ] );
01103 }
01104
01105 ASSERT( nCharactersWritten <= 10 );
01106 chConfiguration[10] = '\0';
01107
01108 strncpy( StatesElemNEW[nelem][nelem-ipISO][ipLo].chConfig, chConfiguration, 10 );
01109 }
01110 }
01111 }
01112 }
01113 return;
01114 }
01115
01116 #if defined(__ICC) && defined(__i386)
01117 #pragma optimization_level 1
01118 #endif
01119 STATIC void FillExtraLymanLine( transition *t, long ipISO, long nelem, long nHi )
01120 {
01121 double Enerwn, Aul;
01122
01123 DEBUG_ENTRY( "FillExtraLymanLine()" );
01124
01125
01126 t->Hi->nelem = (int)(nelem+1);
01127 t->Hi->IonStg = (int)(nelem+1-ipISO);
01128
01129
01130 t->Hi->g = StatesElemNEW[nelem][nelem-ipISO][iso.nLyaLevel[ipISO]].g;
01131
01132
01133 Enerwn = iso.xIsoLevNIonRyd[ipISO][nelem][0] * RYD_INF * ( 1. - 1./POW2((double)nHi) );
01134
01135
01136 t->EnergyWN = (realnum)(Enerwn);
01137 t->EnergyErg = (realnum)(Enerwn * WAVNRYD * EN1RYD);
01138 t->EnergyK = (realnum)(Enerwn * WAVNRYD * TE1RYD);
01139 t->WLAng = (realnum)(1.0e8/ Enerwn/ RefIndex(Enerwn));
01140
01141 if( ipISO == ipH_LIKE )
01142 {
01143 Aul = H_Einstein_A( nHi, 1, 1, 0, nelem+1 );
01144 }
01145 else
01146 {
01147 if( nelem == ipHELIUM )
01148 {
01149
01150 Aul = (1.508e10) / pow((double)nHi,2.975);
01151 }
01152 else
01153 {
01154
01155
01156
01157
01158 Aul = 1.375E10 * pow((double)nelem, 3.9) / pow((double)nHi,3.1);
01159 }
01160 }
01161
01162 t->Emis->Aul = (realnum)Aul;
01163
01164 t->Hi->lifetime = iso_state_lifetime( ipISO, nelem, nHi, 1 );
01165
01166 t->Emis->dampXvel = (realnum)( 1.f / t->Hi->lifetime / PI4 / t->EnergyWN );
01167
01168 t->Emis->iRedisFun = iso.ipResoRedist[ipISO];
01169
01170 t->Emis->gf = (realnum)(GetGF(t->Emis->Aul, t->EnergyWN, t->Hi->g));
01171
01172
01173 t->Emis->opacity = (realnum)(abscf(t->Emis->gf, t->EnergyWN, t->Lo->g));
01174
01175
01176 t->ipCont = INT_MIN;
01177 t->Emis->ipFine = INT_MIN;
01178
01179 {
01180
01181
01182
01183 enum {DEBUG_LOC=false};
01184 if( DEBUG_LOC )
01185 {
01186 fprintf(ioQQQ,"%li\t%li\t%.2e\t%.2e\n",
01187 nelem+1,
01188 nHi,
01189 t->Emis->Aul ,
01190 t->Emis->opacity
01191 );
01192 }
01193 }
01194 return;
01195 }
01196
01197
01198 double iso_state_lifetime( long ipISO, long nelem, long n, long l )
01199 {
01200
01201
01202 double tau, t0, eps2;
01203
01204 double m = ELECTRON_MASS;
01205
01206 double M = (double)dense.AtomicWeight[nelem] * ATOMIC_MASS_UNIT;
01207 double mu = (m*M)/(M+m);
01208 long z = 1;
01209 long Z = nelem + 1 - ipISO;
01210
01211 DEBUG_ENTRY( "iso_state_lifetime()" );
01212
01213
01214 ASSERT( l > 0 );
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
01229 tau /= 3.;
01230
01231 tau *= 1.1722 * pow( (double)nelem, 0.1 );
01232 }
01233
01234
01235 ASSERT( ipISO <= ipHE_LIKE );
01236 ASSERT( tau > 0. );
01237
01238 return tau;
01239 }
01240
01241
01242 void iso_cascade( long ipISO, long nelem )
01243 {
01244
01245
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
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
01265
01266 for( ipHi=1; ipHi<iso.numLevels_max[ipISO][nelem]; ipHi++ )
01267 {
01268 double SumAs = 0.;
01269
01275
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
01313
01314
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
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
01331 iso.SigmaAtot[ipISO][nelem][ipHi] = sqrt( iso.SigmaAtot[ipISO][nelem][ipHi] );
01332 }
01333
01334
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
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]*StatesElemNEW[nelem][nelem-ipISO][ipHi].lifetime, 2.) +
01351 pow(iso.SigmaAtot[ipISO][nelem][ipHi]*iso.BranchRatio[ipISO][nelem][ipHi][i]*iso.CascadeProb[ipISO][nelem][i][ipLo]*StatesElemNEW[nelem][nelem-ipISO][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
01366
01367 {
01368 enum {DEBUG_LOC=false};
01369
01370 if( DEBUG_LOC && (nelem == ipHELIUM) && (ipISO==ipHE_LIKE) )
01371 {
01372
01373 long int hi_l,hi_s;
01374 double Bm;
01375
01376
01377
01378 hi_s = -100000;
01379 hi_l = -100000;
01380 ipLo = -100000;
01381
01382
01383
01384
01385
01386
01387
01388
01389
01390 ASSERT( hi_l != StatesElemNEW[nelem][nelem-ipISO][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
01398 if( (StatesElemNEW[nelem][nelem-ipISO][ipHi].l == 1) && (StatesElemNEW[nelem][nelem-ipISO][ipHi].S == 3) )
01399 {
01400 fprintf(ioQQQ,"\n%ld\t",StatesElemNEW[nelem][nelem-ipISO][ipHi].n);
01401 j = 0;
01402 Bm = 0;
01403 for( i = ipLo; i<=ipHi; i++)
01404 {
01405 if( (StatesElemNEW[nelem][nelem-ipISO][i].l == hi_l) && (StatesElemNEW[nelem][nelem-ipISO][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
01439
01440
01441 for( i=2; i < iso.n_HighestResolved_max[ipISO][nelem]; ++i)
01442 {
01443 ASSERT( (SumAPerN[i] > SumAPerN[i+1]) || opac.lgCaseB );
01444 }
01445
01446 {
01447 enum {DEBUG_LOC=false};
01448 if( DEBUG_LOC )
01449 {
01450 for( i = 2; i<= (iso.n_HighestResolved_max[ipISO][nelem] + iso.nCollapsed_max[ipISO][nelem]); ++i)
01451 {
01452 fprintf(ioQQQ,"n %ld\t lifetime %.4e\n", i, 1./SumAPerN[i]);
01453 }
01454 }
01455 }
01456
01457 free( SumAPerN );
01458
01459 return;
01460 }
01461
01463
01464
01465 STATIC void iso_satellite( void )
01466 {
01467 long i, ipISO, nelem;
01468
01469 DEBUG_ENTRY( "iso_satellite()" );
01470
01471 for( ipISO = ipHE_LIKE; ipISO < NISO; ipISO++ )
01472 {
01473 for( nelem = ipISO; nelem < LIMELM; nelem++ )
01474 {
01475
01476 if( nelem == ipISO || dense.lgElmtOn[nelem] )
01477 {
01478 for( i=0; i<iso.numLevels_max[ipISO][nelem]; i++ )
01479 {
01480 char chLab[5]=" ";
01481
01482 SatelliteLines[ipISO][nelem][i].Zero();
01483
01484
01485
01486
01487
01488 SatelliteLines[ipISO][nelem][i].WLAng = (realnum)(RYDLAM/
01489 ((iso.xIsoLevNIonRyd[ipISO-1][nelem][0] - iso.xIsoLevNIonRyd[ipISO-1][nelem][1]) -
01490 (iso.xIsoLevNIonRyd[ipISO][nelem][1]- iso.xIsoLevNIonRyd[ipISO][nelem][i])) );
01491
01492 SatelliteLines[ipISO][nelem][i].EnergyWN = 1.e8f /
01493 SatelliteLines[ipISO][nelem][i].WLAng;
01494
01495 SatelliteLines[ipISO][nelem][i].EnergyErg = (realnum)ERG1CM *
01496 SatelliteLines[ipISO][nelem][i].EnergyWN;
01497
01498 SatelliteLines[ipISO][nelem][i].EnergyK = (realnum)T1CM *
01499 SatelliteLines[ipISO][nelem][i].EnergyWN;
01500
01501
01502 sprintf( chLab, "%2s%2ld",elementnames.chElementSym[nelem], nelem+1-ipISO );
01503
01504 SatelliteLines[ipISO][nelem][i].Emis->iRedisFun = ipCRDW;
01505
01506 SatelliteLines[ipISO][nelem][i].Hi->nelem = nelem + 1;
01507 SatelliteLines[ipISO][nelem][i].Hi->IonStg = nelem + 1 - ipISO;
01508 fixit();
01509 SatelliteLines[ipISO][nelem][i].Hi->g = 2.f;
01510 SatelliteLines[ipISO][nelem][i].Lo->g = StatesElemNEW[nelem][nelem-ipISO][i].g;
01511 SatelliteLines[ipISO][nelem][i].Emis->PopOpc =
01512 SatelliteLines[ipISO][nelem][i].Lo->Pop;
01513
01514 SatelliteLines[ipISO][nelem][i].Emis->pump = 0.;
01515
01516 }
01517 }
01518 }
01519 }
01520
01521 return;
01522 }
01523
01524 void iso_satellite_update( void )
01525 {
01526 double ConBoltz, LTE_pop=SMALLFLOAT+FLT_EPSILON, factor1, ConvLTEPOP;
01527 long i, ipISO, nelem;
01528 double dr_rate;
01529
01530 DEBUG_ENTRY( "iso_satellite()" );
01531
01532 for( ipISO = ipHE_LIKE; ipISO < NISO; ipISO++ )
01533 {
01534 for( nelem = ipISO; nelem < LIMELM; nelem++ )
01535 {
01536
01537 if( nelem == ipISO || dense.lgElmtOn[nelem] )
01538 {
01539 for( i=0; i<iso.numLevels_max[ipISO][nelem]; i++ )
01540 {
01541 dr_rate = MAX2( iso.SmallA, iso.DielecRecomb[ipISO][nelem][i] * iso.lgDielRecom[ipISO] );
01542
01543 SatelliteLines[ipISO][nelem][i].Emis->phots =
01544 dr_rate * dense.eden * dense.xIonDense[nelem][nelem+1-ipISO];
01545
01546 SatelliteLines[ipISO][nelem][i].Emis->xIntensity =
01547 SatelliteLines[ipISO][nelem][i].Emis->phots *
01548 ERG1CM * SatelliteLines[ipISO][nelem][i].EnergyWN;
01549
01550
01551
01552
01553 factor1 = HION_LTE_POP*dense.AtomicWeight[nelem]/
01554 (dense.AtomicWeight[nelem]+ELECTRON_MASS/ATOMIC_MASS_UNIT);
01555
01556
01557 ConvLTEPOP = pow(factor1,1.5)/(2.*iso.stat_ion[ipISO])/phycon.te32;
01558
01559
01560
01561
01562 ConBoltz = dsexp(iso.xIsoLevNIonRyd[ipISO-1][nelem][1]/phycon.te_ryd);
01563
01564 if( ConBoltz >= SMALLDOUBLE )
01565 {
01566
01567
01568
01569
01570 LTE_pop = SatelliteLines[ipISO][nelem][i].Hi->g * ConBoltz * ConvLTEPOP;
01571 }
01572
01573 LTE_pop = max( LTE_pop, 1e-30f );
01574
01575
01576 SatelliteLines[ipISO][nelem][i].Emis->Aul = (realnum)(dr_rate/LTE_pop);
01577 SatelliteLines[ipISO][nelem][i].Emis->Aul =
01578 max( iso.SmallA, SatelliteLines[ipISO][nelem][i].Emis->Aul );
01579
01580 SatelliteLines[ipISO][nelem][i].Emis->gf = (realnum)GetGF(
01581 SatelliteLines[ipISO][nelem][i].Emis->Aul,
01582 SatelliteLines[ipISO][nelem][i].EnergyWN,
01583 SatelliteLines[ipISO][nelem][i].Hi->g);
01584
01585 SatelliteLines[ipISO][nelem][i].Emis->gf =
01586 max( 1e-20f, SatelliteLines[ipISO][nelem][i].Emis->gf );
01587
01588 SatelliteLines[ipISO][nelem][i].Hi->Pop = LTE_pop * dense.xIonDense[nelem][nelem+1-ipISO] * dense.eden;
01589
01590 SatelliteLines[ipISO][nelem][i].Emis->PopOpc =
01591 SatelliteLines[ipISO][nelem][i].Lo->Pop -
01592 SatelliteLines[ipISO][nelem][i].Hi->Pop *
01593 SatelliteLines[ipISO][nelem][i].Lo->g/SatelliteLines[ipISO][nelem][i].Hi->g;
01594
01595 SatelliteLines[ipISO][nelem][i].Emis->opacity =
01596 (realnum)(abscf(SatelliteLines[ipISO][nelem][i].Emis->gf,
01597 SatelliteLines[ipISO][nelem][i].EnergyWN,
01598 SatelliteLines[ipISO][nelem][i].Lo->g));
01599
01600
01601 double lifetime = 1e-10;
01602
01603 SatelliteLines[ipISO][nelem][i].Emis->dampXvel = (realnum)(
01604 (1.f/lifetime)/PI4/SatelliteLines[ipISO][nelem][i].EnergyWN);
01605 }
01606 }
01607 }
01608 }
01609
01610 return;
01611 }
01612
01613 long iso_get_total_num_levels( long ipISO, long nmaxResolved, long numCollapsed )
01614 {
01615 DEBUG_ENTRY( "iso_get_total_num_levels()" );
01616
01617 long tot_num_levels;
01618
01619
01620
01621
01622 if( ipISO == ipH_LIKE )
01623 {
01624 tot_num_levels = (long)( nmaxResolved * 0.5 *( nmaxResolved + 1 ) ) + numCollapsed;
01625 }
01626 else if( ipISO == ipHE_LIKE )
01627 {
01628 tot_num_levels = nmaxResolved*nmaxResolved + nmaxResolved + 1 + numCollapsed;
01629 }
01630 else
01631 TotalInsanity();
01632
01633 return tot_num_levels;
01634 }
01635
01636 void iso_update_num_levels( long ipISO, long nelem )
01637 {
01638 DEBUG_ENTRY( "iso_update_num_levels()" );
01639
01640
01641 ASSERT( iso.n_HighestResolved_max[ipISO][nelem] >= 3 );
01642
01643 iso.numLevels_max[ipISO][nelem] =
01644 iso_get_total_num_levels( ipISO, iso.n_HighestResolved_max[ipISO][nelem], iso.nCollapsed_max[ipISO][nelem] );
01645
01646 if( iso.numLevels_max[ipISO][nelem] > iso.numLevels_malloc[ipISO][nelem] )
01647 {
01648 fprintf( ioQQQ, "The number of levels for ipISO %li, nelem %li, has been increased since the initial coreload.\n",
01649 ipISO, nelem );
01650 fprintf( ioQQQ, "This cannot be done.\n" );
01651 cdEXIT(EXIT_FAILURE);
01652 }
01653
01654
01655 iso.numLevels_local[ipISO][nelem] = iso.numLevels_max[ipISO][nelem];
01656 iso.nCollapsed_local[ipISO][nelem] = iso.nCollapsed_max[ipISO][nelem];
01657 iso.n_HighestResolved_local[ipISO][nelem] = iso.n_HighestResolved_max[ipISO][nelem];
01658
01659
01660
01661
01662 max_num_levels = MAX2( max_num_levels, iso.numLevels_max[ipISO][nelem]);
01663
01664 return;
01665 }
01666
01667 void iso_collapsed_bnl_set( long ipISO, long nelem )
01668 {
01669
01670 DEBUG_ENTRY( "iso_collapsed_bnl_set()" );
01671
01672 double bnl_array[4][3][4][10] = {
01673 {
01674
01675 {
01676 {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},
01677 {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},
01678 {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},
01679 {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}
01680 },
01681 {
01682 {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},
01683 {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},
01684 {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},
01685 {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}
01686 },
01687 {
01688 {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},
01689 {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},
01690 {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},
01691 {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}
01692 }
01693 },
01694 {
01695
01696 {
01697 {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},
01698 {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},
01699 {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},
01700 {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}
01701 },
01702 {
01703 {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},
01704 {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},
01705 {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},
01706 {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}
01707 },
01708 {
01709 {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},
01710 {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},
01711 {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},
01712 {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}
01713 }
01714 },
01715 {
01716
01717 {
01718 {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},
01719 {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},
01720 {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},
01721 {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}
01722 },
01723 {
01724 {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},
01725 {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},
01726 {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},
01727 {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}
01728 },
01729 {
01730 {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},
01731 {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},
01732 {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},
01733 {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}
01734 }
01735 },
01736 {
01737
01738 {
01739 {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},
01740 {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},
01741 {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},
01742 {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}
01743 },
01744 {
01745 {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},
01746 {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},
01747 {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},
01748 {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}
01749 },
01750 {
01751 {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},
01752 {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},
01753 {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},
01754 {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}
01755 }
01756 }
01757 };
01758
01759 double temps[4] = {5000., 10000., 15000., 20000. };
01760 double log_dens[3] = {2., 4., 6.};
01761 long ipTe, ipDens;
01762
01763 ASSERT( nelem <= 1 );
01764
01765
01766 ipTe = hunt_bisect( temps, 4, phycon.te );
01767 ipDens = hunt_bisect( log_dens, 3, log10(dense.eden) );
01768
01769 ASSERT( (ipTe >=0) && (ipTe < 3) );
01770 ASSERT( (ipDens >=0) && (ipDens < 2) );
01771
01772 for( long nHi=iso.n_HighestResolved_max[ipISO][nelem]+1; nHi<=iso.n_HighestResolved_max[ipISO][nelem]+iso.nCollapsed_max[ipISO][nelem]; nHi++ )
01773 {
01774 for( long lHi=0; lHi<nHi; lHi++ )
01775 {
01776 for( long sHi=1; sHi<4; sHi++ )
01777 {
01778 if( ipISO == ipH_LIKE && sHi != 2 )
01779 continue;
01780 else if( ipISO == ipHE_LIKE && sHi != 1 && sHi != 3 )
01781 continue;
01782
01783 double bnl_at_lo_den, bnl_at_hi_den, bnl;
01784 double bnl_max, bnl_min, temp, dens;
01785
01786 long ipL = MIN2(9,lHi);
01787 long ip1;
01788
01789 if( nelem==ipHYDROGEN )
01790 ip1 = 0;
01791 else if( nelem==ipHELIUM )
01792 {
01793 if( ipISO==ipH_LIKE )
01794 ip1 = 1;
01795 else if( ipISO==ipHE_LIKE )
01796 {
01797 if( sHi==1 )
01798 ip1 = 2;
01799 else if( sHi==3 )
01800 ip1 = 3;
01801 else
01802 TotalInsanity();
01803 }
01804 else
01805 TotalInsanity();
01806 }
01807 else
01808 TotalInsanity();
01809
01810 temp = MAX2( temps[0], phycon.te );
01811 temp = MIN2( temps[3], temp );
01812
01813 dens = MAX2( log_dens[0], log10(dense.eden) );
01814 dens = MIN2( log_dens[2], dens );
01815
01816
01817
01818
01819
01820 if( temp < temps[0] && dens < log_dens[0] )
01821 bnl = bnl_array[ip1][0][0][ipL];
01822 else if( temp < temps[0] && dens >= log_dens[2] )
01823 bnl = bnl_array[ip1][2][0][ipL];
01824 else if( temp >= temps[3] && dens < log_dens[0] )
01825 bnl = bnl_array[ip1][0][3][ipL];
01826 else if( temp >= temps[3] && dens >= log_dens[2] )
01827 bnl = bnl_array[ip1][2][3][ipL];
01828 else
01829 {
01830 bnl_at_lo_den = ( temp - temps[ipTe]) / (temps[ipTe+1] - temps[ipTe]) *
01831 (bnl_array[ip1][ipDens][ipTe+1][ipL] - bnl_array[ip1][ipDens][ipTe][ipL]) + bnl_array[ip1][ipDens][ipTe][ipL];
01832
01833 bnl_at_hi_den = ( temp - temps[ipTe]) / (temps[ipTe+1] - temps[ipTe]) *
01834 (bnl_array[ip1][ipDens+1][ipTe+1][ipL] - bnl_array[ip1][ipDens+1][ipTe][ipL]) + bnl_array[ip1][ipDens+1][ipTe][ipL];
01835
01836 bnl = ( dens - log_dens[ipDens]) / (log_dens[ipDens+1] - log_dens[ipDens]) *
01837 (bnl_at_hi_den - bnl_at_lo_den) + bnl_at_lo_den;
01838 }
01839
01841 {
01842 bnl_max = MAX4( bnl_array[ip1][ipDens][ipTe+1][ipL], bnl_array[ip1][ipDens+1][ipTe+1][ipL],
01843 bnl_array[ip1][ipDens][ipTe][ipL], bnl_array[ip1][ipDens+1][ipTe][ipL] );
01844 ASSERT( bnl <= bnl_max );
01845
01846 bnl_min = MIN4( 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_min );
01849 }
01850
01851 iso.bnl_effective[ipISO][nelem][nHi][lHi][sHi] = bnl;
01852
01853 ASSERT( iso.bnl_effective[ipISO][nelem][nHi][lHi][sHi] > 0. );
01854 }
01855 }
01856 }
01857
01858 return;
01859 }
01860
01861
01862 void iso_collapsed_bnl_print( long ipISO, long nelem )
01863 {
01864 DEBUG_ENTRY( "iso_collapsed_bnl_print()" );
01865
01866 for( long is = 1; is<=3; ++is)
01867 {
01868 if( ipISO == ipH_LIKE && is != 2 )
01869 continue;
01870 else if( ipISO == ipHE_LIKE && is != 1 && is != 3 )
01871 continue;
01872
01873 char chSpin[3][9]= {"singlets", "doublets", "triplets"};
01874
01875
01876 fprintf(ioQQQ," %s %s %s bnl\n",
01877 iso.chISO[ipISO],
01878 elementnames.chElementSym[nelem],
01879 chSpin[is-1]);
01880
01881
01882 fprintf(ioQQQ," n\\l=> ");
01883 for( long i =0; i < iso.n_HighestResolved_max[ipISO][nelem] + iso.nCollapsed_max[ipISO][nelem]; ++i)
01884 {
01885 fprintf(ioQQQ,"%2ld ",i);
01886 }
01887 fprintf(ioQQQ,"\n");
01888
01889
01890 for( long in = 1; in <= iso.n_HighestResolved_max[ipISO][nelem] + iso.nCollapsed_max[ipISO][nelem]; ++in)
01891 {
01892 if( is==3 && in==1 )
01893 continue;
01894
01895 fprintf(ioQQQ," %2ld ",in);
01896
01897 for( long il = 0; il < in; ++il)
01898 {
01899 fprintf( ioQQQ, "%9.3e ", iso.bnl_effective[ipISO][nelem][in][il][is] );
01900 }
01901 fprintf(ioQQQ,"\n");
01902 }
01903 }
01904
01905 return;
01906 }
01907
01908 void iso_collapsed_Aul_update( long ipISO, long nelem )
01909 {
01910 DEBUG_ENTRY( "iso_collapsed_Aul_update()" );
01911
01912 long ipFirstCollapsed = iso.numLevels_max[ipISO][nelem] - iso.nCollapsed_max[ipISO][nelem];
01913
01914 for( long ipLo=0; ipLo<ipFirstCollapsed; ipLo++ )
01915 {
01916 long spin = StatesElemNEW[nelem][nelem-ipISO][ipLo].S;
01917
01918
01919 for( long nHi=iso.n_HighestResolved_max[ipISO][nelem]+1; nHi<=iso.n_HighestResolved_max[ipISO][nelem]+iso.nCollapsed_max[ipISO][nelem]; nHi++ )
01920 {
01921 realnum Auls[2] = {
01922 iso.CachedAs[ipISO][nelem][ nHi-iso.n_HighestResolved_max[ipISO][nelem]-1 ][ ipLo ][0],
01923 iso.CachedAs[ipISO][nelem][ nHi-iso.n_HighestResolved_max[ipISO][nelem]-1 ][ ipLo ][1] };
01924
01925 realnum EffectiveAul =
01926 Auls[0]*spin*(2.f*(L_(ipLo)+1.f)+1.f)*(realnum)iso.bnl_effective[ipISO][nelem][nHi][ L_(ipLo)+1 ][spin];
01927
01928
01929
01930 if( L_(ipLo) > 0 )
01931 {
01932 EffectiveAul +=
01933 Auls[1]*spin*(2.f*(L_(ipLo)-1.f)+1.f)*(realnum)iso.bnl_effective[ipISO][nelem][nHi][ L_(ipLo)-1 ][spin];
01934 }
01935
01936 if( ipISO==ipH_LIKE )
01937 EffectiveAul /= (2.f*nHi*nHi);
01938 else if( ipISO==ipHE_LIKE )
01939 EffectiveAul /= (4.f*nHi*nHi);
01940 else
01941 TotalInsanity();
01942
01943 long ipHi = iso.QuantumNumbers2Index[ipISO][nelem][nHi][ L_(ipLo)+1 ][spin];
01944
01945
01946 Transitions[ipISO][nelem][ipHi][ipLo].Emis->Aul = EffectiveAul;
01947
01948 ASSERT( Transitions[ipISO][nelem][ipHi][ipLo].Emis->Aul > 0. );
01949 }
01950 }
01951
01952 return;
01953 }
01954
01955 void iso_collapsed_lifetimes_update( long ipISO, long nelem )
01956 {
01957 DEBUG_ENTRY( "iso_collapsed_Aul_update()" );
01958
01959 for( long ipHi=iso.numLevels_max[ipISO][nelem]- iso.nCollapsed_max[ipISO][nelem]; ipHi < iso.numLevels_max[ipISO][nelem]; ipHi++ )
01960 {
01961 StatesElemNEW[nelem][nelem-ipISO][ipHi].lifetime = SMALLFLOAT;
01962
01963 for( long ipLo=0; ipLo < ipHi; ipLo++ )
01964 {
01965 if( Transitions[ipISO][nelem][ipHi][ipLo].Emis->Aul <= iso.SmallA )
01966 continue;
01967
01968 StatesElemNEW[nelem][nelem-ipISO][ipHi].lifetime += Transitions[ipISO][nelem][ipHi][ipLo].Emis->Aul;
01969 }
01970
01971
01972 StatesElemNEW[nelem][nelem-ipISO][ipHi].lifetime = 1./StatesElemNEW[nelem][nelem-ipISO][ipHi].lifetime;
01973
01974 for( long ipLo=0; ipLo < ipHi; ipLo++ )
01975 {
01976 if( Transitions[ipISO][nelem][ipHi][ipLo].EnergyWN <= 0. )
01977 continue;
01978
01979 if( Transitions[ipISO][nelem][ipHi][ipLo].Emis->Aul <= iso.SmallA )
01980 continue;
01981
01982 Transitions[ipISO][nelem][ipHi][ipLo].Emis->dampXvel = (realnum)(
01983 (1.f/StatesElemNEW[nelem][nelem-ipISO][ipHi].lifetime)/
01984 PI4/Transitions[ipISO][nelem][ipHi][ipLo].EnergyWN);
01985
01986 ASSERT(Transitions[ipISO][nelem][ipHi][ipLo].Emis->dampXvel> 0.);
01987 }
01988 }
01989
01990 return;
01991 }
01992
01993 #if 0
01994 STATIC void Prt_AGN_table( void )
01995 {
01996
01997 multi_arr<char,2> chLevel(max_num_levels,10);
01998
01999
02000 for( long ipLo=0; ipLo < iso.numLevels_max[ipISO][ipISO]-iso.nCollapsed_max[ipISO][ipISO]; ++ipLo )
02001 {
02002 long nelem = ipISO;
02003 sprintf( &chLevel.front(ipLo) , "%li %li%c", N_(ipLo), S_(ipLo), chL[MIN2(20,L_(ipLo))] );
02004 }
02005
02006
02007
02008 {
02009
02010 enum {AGN=false};
02011 if( AGN )
02012 {
02013 # define NTEMP 6
02014 double te[NTEMP]={6000.,8000.,10000.,15000.,20000.,25000. };
02015 double telog[NTEMP] ,
02016 cs ,
02017 ratecoef;
02018 long nelem = ipHELIUM;
02019 fprintf( ioQQQ,"trans");
02020 for( long i=0; i < NTEMP; ++i )
02021 {
02022 telog[i] = log10( te[i] );
02023 fprintf( ioQQQ,"\t%.3e",te[i]);
02024 }
02025 for( long i=0; i < NTEMP; ++i )
02026 {
02027 fprintf( ioQQQ,"\t%.3e",te[i]);
02028 }
02029 fprintf(ioQQQ,"\n");
02030
02031 for( long ipHi=ipHe2s3S; ipHi< iso.numLevels_max[ipHE_LIKE][ipHELIUM]; ++ipHi )
02032 {
02033 for( long ipLo=ipHe1s1S; ipLo < ipHi; ++ipLo )
02034 {
02035
02036
02037
02038 if( N_(ipHi) == N_(ipLo) )
02039 continue;
02040
02041
02042 fprintf( ioQQQ,"%s - %s",
02043 &chLevel.front(ipLo) , &chLevel.front(ipHi) );
02044
02045
02046 for( long i=0; i < NTEMP; ++i )
02047 {
02048 phycon.alogte = telog[i];
02049
02050 cs = HeCSInterp( nelem , ipHi , ipLo, ipELECTRON );
02051 fprintf(ioQQQ,"\t%.2e", cs );
02052 }
02053
02054
02055 for( long i=0; i < NTEMP; ++i )
02056 {
02057 phycon.alogte = telog[i];
02058 phycon.te = pow(10.,telog[i] );
02059 tfidle(false);
02060 cs = HeCSInterp( nelem , ipHi , ipLo, ipELECTRON );
02061
02062 ratecoef = cs/sqrt(phycon.te)*COLL_CONST/StatesElemNEW[nelem][nelem-ipHE_LIKE][ipLo].g *
02063 sexp( Transitions[ipHE_LIKE][nelem][ipHi][ipLo].EnergyK / phycon.te );
02064 fprintf(ioQQQ,"\t%.2e", ratecoef );
02065 }
02066 fprintf(ioQQQ,"\n");
02067 }
02068 }
02069 cdEXIT(EXIT_FAILURE);
02070 }
02071 }
02072
02073 return;
02074 }
02075 #endif