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( const TransitionList::iterator& 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
00044 static int nCalled = 0;
00045
00046 double HIonPoten;
00047
00048 DEBUG_ENTRY( "iso_create()" );
00049
00050
00051 if( nCalled > 0 )
00052 {
00053 iso_zero();
00054 return;
00055 }
00056
00057
00058 ++nCalled;
00059
00060
00061 iso_ctrl.stat_ion[ipH_LIKE] = 1.f;
00062 iso_ctrl.stat_ion[ipHE_LIKE] = 2.f;
00063
00064
00065
00066 iso_allocate();
00067
00068
00069 iso_assign_quantum_numbers();
00070
00071
00072 (*TauDummy).Junk();
00073 (*TauDummy).AddHiState();
00074 (*TauDummy).AddLoState();
00075 (*TauDummy).AddLine2Stack();
00076
00077
00078
00079
00080 for( long ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
00081 {
00082
00083 for( long nelem=ipISO; nelem < LIMELM; nelem++ )
00084 {
00085
00086 if( nelem < 2 || dense.lgElmtOn[nelem] )
00087 {
00088
00089
00090
00091 HIonPoten = t_ADfA::Inst().ph1(0,0,nelem,0)/EVRYD* 0.9998787;
00092 ASSERT(HIonPoten > 0.);
00093
00094 double EnergyRydGround = 0.;
00095
00096 for( ipHi=0; ipHi < iso_sp[ipISO][nelem].numLevels_max; ipHi++ )
00097 {
00098 double EnergyWN, EnergyRyd;
00099
00100 if( ipISO == ipH_LIKE )
00101 {
00102 EnergyRyd = HIonPoten/POW2((double)N_(ipHi));
00103 }
00104 else if( ipISO == ipHE_LIKE )
00105 {
00106 EnergyRyd = helike_energy( nelem, ipHi ) * WAVNRYD;
00107 }
00108 else
00109 {
00110
00111 TotalInsanity();
00112 }
00113
00114
00115 ASSERT(EnergyRyd >= 0.);
00116
00117 iso_sp[ipISO][nelem].fb[ipHi].xIsoLevNIonRyd = EnergyRyd;
00118 if (ipHi == 0)
00119 EnergyRydGround = EnergyRyd;
00120 iso_sp[ipISO][nelem].st[ipHi].energy().set(EnergyRydGround-EnergyRyd);
00121
00122
00123 for( ipLo=0; ipLo < ipHi; ipLo++ )
00124 {
00125 EnergyWN = RYD_INF * (iso_sp[ipISO][nelem].fb[ipLo].xIsoLevNIonRyd -
00126 iso_sp[ipISO][nelem].fb[ipHi].xIsoLevNIonRyd);
00127
00128
00129
00130 if( EnergyWN==0 && ipISO==ipHE_LIKE )
00131 EnergyWN = 0.0001;
00132
00133 if( EnergyWN < 0. )
00134 EnergyWN = -1.0 * EnergyWN;
00135
00136
00137 iso_sp[ipISO][nelem].trans(ipHi,ipLo).EnergyWN() = (realnum)EnergyWN;
00138
00139 ASSERT(iso_sp[ipISO][nelem].trans(ipHi,ipLo).EnergyWN() >= 0.);
00140 ASSERT(iso_sp[ipISO][nelem].trans(ipHi,ipLo).EnergyErg() >= 0.);
00141 ASSERT(iso_sp[ipISO][nelem].trans(ipHi,ipLo).EnergyK() >= 0.);
00142
00144 if( N_(ipLo)==N_(ipHi) && ipISO==ipH_LIKE )
00145 {
00146 iso_sp[ipISO][nelem].trans(ipHi,ipLo).WLAng() = 0.;
00147 }
00148 else
00149 {
00150
00151 iso_sp[ipISO][nelem].trans(ipHi,ipLo).WLAng() =
00152 (realnum)(1.0e8/
00153 iso_sp[ipISO][nelem].trans(ipHi,ipLo).EnergyWN()/
00154 RefIndex( iso_sp[ipISO][nelem].trans(ipHi,ipLo).EnergyWN()));
00155 ASSERT(iso_sp[ipISO][nelem].trans(ipHi,ipLo).WLAng() > 0.);
00156 }
00157
00158 }
00159 }
00160
00161
00162 for( ipHi=2; ipHi < iso_ctrl.nLyman_malloc[ipISO]; ipHi++ )
00163 {
00164 FillExtraLymanLine( ExtraLymanLines[ipISO][nelem].begin()+ipExtraLymanLines[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( long 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( long nelem=ipISO; nelem < LIMELM; nelem++ )
00202 {
00203
00204 if( nelem < 2 || dense.lgElmtOn[nelem] )
00205 {
00206 for( ipLo=ipH1s; ipLo < (iso_sp[ipISO][nelem].numLevels_max - 1); ipLo++ )
00207 {
00208 for( ipHi=ipLo + 1; ipHi < iso_sp[ipISO][nelem].numLevels_max; 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_ctrl.SmallA )
00227 iso_sp[ipISO][nelem].trans(ipHi,ipLo).ipEmis() = -1;
00228 else
00229 iso_sp[ipISO][nelem].trans(ipHi,ipLo).AddLine2Stack();
00230
00231 iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().Aul() = Aul;
00232
00233 ASSERT(iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().Aul() > 0.);
00234
00235 if( ipLo == 0 && ipHi == iso_ctrl.nLyaLevel[ipISO] )
00236 {
00237 long redis = iso_ctrl.ipLyaRedist[ipISO];
00238
00239 if( ipISO==ipH_LIKE && nelem==ipHYDROGEN )
00240 redis = ipLY_A;
00241 iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().iRedisFun() = redis;
00242 }
00243 else if( ipLo == 0 )
00244 {
00245
00246
00247 iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().iRedisFun() = iso_ctrl.ipResoRedist[ipISO];
00248 }
00249 else
00250 {
00251
00252
00253 iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().iRedisFun() = iso_ctrl.ipSubRedist[ipISO];
00254 }
00255
00256 if( iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().Aul() <= iso_ctrl.SmallA ||
00257 iso_sp[ipISO][nelem].trans(ipHi,ipLo).EnergyWN() <= 0.)
00258 {
00259 iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().gf() = 0.;
00260 iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().opacity() = 0.;
00261 }
00262 else
00263 {
00264 iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().gf() =
00265 (realnum)(GetGF(iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().Aul(),
00266 iso_sp[ipISO][nelem].trans(ipHi,ipLo).EnergyWN(),
00267 iso_sp[ipISO][nelem].st[ipHi].g()));
00268 ASSERT(iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().gf() > 0.);
00269
00270
00271 iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().opacity() =
00272 (realnum)(abscf(iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().gf(),
00273 iso_sp[ipISO][nelem].trans(ipHi,ipLo).EnergyWN(),
00274 iso_sp[ipISO][nelem].st[ipLo].g()));
00275 ASSERT(iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().opacity() > 0.);
00276 }
00277 }
00278 }
00279 }
00280 }
00281 }
00282
00283
00284
00285
00286 if( iso_ctrl.lgFSM[ipHE_LIKE] )
00287 {
00288
00289 for( long nelem=ipHE_LIKE; nelem < LIMELM; nelem++ )
00290 {
00291
00292 if( nelem < 2 || dense.lgElmtOn[nelem] )
00293 {
00294 for( ipHi=ipHe2s3S; ipHi<iso_sp[ipHE_LIKE][nelem].numLevels_max; ipHi++ )
00295 {
00296 for( ipLo=ipHe1s1S; ipLo<ipHi; ipLo++ )
00297 {
00298 DoFSMixing( nelem, ipLo, ipHi );
00299 }
00300 }
00301 }
00302 }
00303 }
00304
00305
00306 iso_sp[ipH_LIKE][ipHELIUM].trans(ipH3s,ipH2s).WLAng() = 1640.f;
00307 iso_sp[ipH_LIKE][ipHELIUM].trans(ipH3s,ipH2p).WLAng() = 1640.f;
00308 if( iso_sp[ipH_LIKE][ipHELIUM].n_HighestResolved_max >=3 )
00309 {
00310 iso_sp[ipH_LIKE][ipHELIUM].trans(ipH3p,ipH2s).WLAng() = 1640.f;
00311 iso_sp[ipH_LIKE][ipHELIUM].trans(ipH3p,ipH2p).WLAng() = 1640.f;
00312 iso_sp[ipH_LIKE][ipHELIUM].trans(ipH3d,ipH2s).WLAng() = 1640.f;
00313 iso_sp[ipH_LIKE][ipHELIUM].trans(ipH3d,ipH2p).WLAng() = 1640.f;
00314 }
00315
00316
00317
00318
00319 for( long ipISO=ipH_LIKE; ipISO<NISO; ipISO++ )
00320 {
00321 for( long nelem=ipISO; nelem < LIMELM; nelem++ )
00322 {
00323
00324 if( nelem < 2 || dense.lgElmtOn[nelem] )
00325 {
00326
00327 iso_sp[ipISO][nelem].st[0].lifetime() = -FLT_MAX;
00328
00329 for( ipHi=1; ipHi < iso_sp[ipISO][nelem].numLevels_max; ipHi++ )
00330 {
00331 iso_sp[ipISO][nelem].st[ipHi].lifetime() = SMALLFLOAT;
00332
00333 for( ipLo=0; ipLo < ipHi; ipLo++ )
00334 {
00335 if( iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().Aul() <= iso_ctrl.SmallA )
00336 continue;
00337
00338 iso_sp[ipISO][nelem].st[ipHi].lifetime() += iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().Aul();
00339 }
00340
00341
00342 iso_sp[ipISO][nelem].st[ipHi].lifetime() = 1./iso_sp[ipISO][nelem].st[ipHi].lifetime();
00343
00344 for( ipLo=0; ipLo < ipHi; ipLo++ )
00345 {
00346 if( iso_sp[ipISO][nelem].trans(ipHi,ipLo).EnergyWN() <= 0. )
00347 continue;
00348
00349 if( iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().Aul() <= iso_ctrl.SmallA )
00350 continue;
00351
00352 iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().dampXvel() = (realnum)(
00353 (1.f/iso_sp[ipISO][nelem].st[ipHi].lifetime())/
00354 PI4/iso_sp[ipISO][nelem].trans(ipHi,ipLo).EnergyWN());
00355
00356 ASSERT(iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().dampXvel()> 0.);
00357 }
00358 }
00359 }
00360 }
00361 }
00362
00363
00364 iso_zero();
00365
00366
00367 for( long ipISO=ipH_LIKE; ipISO<NISO; ipISO++ )
00368 {
00369 for( long nelem = ipISO; nelem < LIMELM; nelem++ )
00370 {
00371
00372 if( nelem == ipISO || dense.lgElmtOn[nelem] )
00373 {
00374
00375 iso_cascade( ipISO, nelem);
00376 }
00377 }
00378 }
00379
00380 iso_satellite();
00381
00382 for( long nelem=ipHYDROGEN; nelem < LIMELM; ++nelem )
00383 iso_satellite_update( nelem );
00384
00385
00386
00387
00388 for( long ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
00389 {
00390 for( long nelem=ipISO; nelem < LIMELM; ++nelem )
00391 {
00392 if( nelem < 2 || dense.lgElmtOn[nelem] )
00393 {
00394 for( long ipHi= 1; ipHi < iso_sp[ipISO][nelem].numLevels_max; ++ipHi )
00395 {
00396 for( long ipLo=0; ipLo < ipHi; ++ipLo )
00397 {
00398 iso_sp[ipISO][nelem].ex[ipHi][ipLo].pestrk = 0.;
00399 iso_sp[ipISO][nelem].ex[ipHi][ipLo].pestrk_up = 0.;
00400 }
00401 }
00402 }
00403 }
00404 }
00405
00406 return;
00407 }
00408
00409
00410 STATIC void iso_zero(void)
00411 {
00412 DEBUG_ENTRY( "iso_zero()" );
00413
00414 hydro.HLineWidth = 0.;
00415
00416
00417
00418
00419 for( long ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
00420 {
00421 for( long nelem=ipISO; nelem < LIMELM; nelem++ )
00422 {
00423 if( nelem < 2 || dense.lgElmtOn[nelem] )
00424 {
00425 for( long ipHi=0; ipHi < iso_sp[ipISO][nelem].numLevels_max; ipHi++ )
00426 {
00427 iso_sp[ipISO][nelem].st[ipHi].Pop() = 0.;
00428 iso_sp[ipISO][nelem].fb[ipHi].Reset();
00429 }
00430 if (ipISO <= nelem)
00431 iso_sp[ipISO][nelem].st[0].Pop() =
00432 dense.xIonDense[nelem][nelem-ipISO];
00433 }
00434
00435 if( nelem < 2 )
00436 {
00437 iso_collapsed_bnl_set( ipISO, nelem );
00438
00439 iso_collapsed_Aul_update( ipISO, nelem );
00440 iso_collapsed_lifetimes_update( ipISO, nelem );
00441 }
00442 }
00443 }
00444
00445
00446
00447 iso_sp[ipH_LIKE][ipHYDROGEN].fb[0].ConOpacRatio = 1e-5;
00448 iso_sp[ipH_LIKE][ipHELIUM].fb[0].ConOpacRatio = 1e-5;
00449 iso_sp[ipHE_LIKE][ipHELIUM].fb[0].ConOpacRatio = 1e-5;
00450
00451 return;
00452 }
00453
00454 STATIC void iso_allocate(void)
00455 {
00456
00457 DEBUG_ENTRY( "iso_allocate()" );
00458
00459
00460 for( long ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
00461 {
00462 for( long nelem=ipISO; nelem < LIMELM; ++nelem )
00463 {
00464
00465 if( nelem < 2 || dense.lgElmtOn[nelem] )
00466 {
00467 t_iso_sp *sp = &iso_sp[ipISO][nelem];
00468 sp->numLevels_malloc = sp->numLevels_max;
00469
00470 ASSERT( sp->numLevels_max > 0 );
00471
00472 iso_ctrl.nLyman_malloc[ipISO] = iso_ctrl.nLyman[ipISO];
00473
00474 sp->CachedAs.reserve( MAX2(1, sp->nCollapsed_max) );
00475
00476 sp->ipTrans.reserve( sp->numLevels_max );
00477 sp->ex.reserve( sp->numLevels_max );
00478 sp->CascadeProb.reserve( sp->numLevels_max );
00479 sp->BranchRatio.reserve( sp->numLevels_max );
00480
00481 sp->fb.resize( sp->numLevels_max );
00482
00483 for( long i = 0; i < sp->nCollapsed_max; ++i )
00484 {
00485 sp->CachedAs.reserve( i, sp->numLevels_max - sp->nCollapsed_max );
00486 for( long i1 = 0; i1 < sp->numLevels_max - sp->nCollapsed_max; ++i1 )
00487 {
00488
00489 sp->CachedAs.reserve( i, i1, 2 );
00490 }
00491 }
00492
00493 sp->QuantumNumbers2Index.reserve( sp->n_HighestResolved_max + sp->nCollapsed_max + 1 );
00494
00495 for( long i = 1; i <= sp->n_HighestResolved_max + sp->nCollapsed_max; ++i )
00496 {
00497
00498 sp->QuantumNumbers2Index.reserve( i, i );
00499
00500 for( long i1 = 0; i1 < i; ++i1 )
00501 {
00502
00503 ASSERT( NISO == 2 );
00504
00505
00506 sp->QuantumNumbers2Index.reserve( i, i1, 4 );
00507 }
00508 }
00509
00510 for( long n=1; n < sp->numLevels_max; ++n )
00511 {
00512 sp->ipTrans.reserve( n, n );
00513 }
00514
00515 for( long n=0; n < sp->numLevels_max; ++n )
00516 {
00517 sp->ex.reserve( n, sp->numLevels_max );
00518 sp->CascadeProb.reserve( n, sp->numLevels_max );
00519 sp->BranchRatio.reserve( n, sp->numLevels_max );
00520 }
00521
00522 sp->ipTrans.alloc();
00523 sp->ex.alloc();
00524 sp->CascadeProb.alloc();
00525 sp->BranchRatio.alloc();
00526
00527 sp->CachedAs.alloc();
00528 sp->QuantumNumbers2Index.alloc();
00529 sp->bnl_effective.alloc( sp->QuantumNumbers2Index.clone() );
00530 sp->QuantumNumbers2Index.invalidate();
00531 sp->bnl_effective.invalidate();
00532 }
00533 }
00534 }
00535
00536 ipSatelliteLines.reserve( NISO );
00537 ipExtraLymanLines.reserve( NISO );
00538
00539 for( long ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
00540 {
00541 ipSatelliteLines.reserve( ipISO, LIMELM );
00542 ipExtraLymanLines.reserve( ipISO, LIMELM );
00543
00544 for( long nelem=ipISO; nelem < LIMELM; ++nelem )
00545 {
00546
00547 if( nelem < 2 || dense.lgElmtOn[nelem] )
00548 {
00549 ASSERT( iso_sp[ipISO][nelem].numLevels_max > 0 );
00550
00551 ipSatelliteLines.reserve( ipISO, nelem, iso_sp[ipISO][nelem].numLevels_max );
00552 ipExtraLymanLines.reserve( ipISO, nelem, iso_ctrl.nLyman_malloc[ipISO] );
00553 }
00554 }
00555 }
00556
00557 ipSatelliteLines.alloc();
00558 ipExtraLymanLines.alloc();
00559
00560 Transitions.resize(NISO);
00561 SatelliteLines.resize(NISO);
00562 ExtraLymanLines.resize(NISO);
00563 for( long ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
00564 {
00565 Transitions[ipISO].reserve(LIMELM);
00566 SatelliteLines[ipISO].reserve(LIMELM);
00567 ExtraLymanLines[ipISO].reserve(LIMELM);
00568 for( long nelem=0; nelem < ipISO; ++nelem )
00569 {
00570 Transitions[ipISO].push_back(
00571 TransitionList("Insanity",&AnonStates));
00572 SatelliteLines[ipISO].push_back(
00573 TransitionList("Insanity",&AnonStates));
00574 ExtraLymanLines[ipISO].push_back(
00575 TransitionList("Insanity",&AnonStates));
00576 }
00577 for( long nelem=ipISO; nelem < LIMELM; ++nelem )
00578 {
00579 if( nelem < 2 || dense.lgElmtOn[nelem] )
00580 {
00581 Transitions[ipISO].push_back(
00582 TransitionList("Isosequence",&iso_sp[ipISO][nelem].st));
00583 SatelliteLines[ipISO].push_back(
00584 TransitionList("SatelliteLines",&iso_sp[ipISO][nelem].st));
00585 ExtraLymanLines[ipISO].push_back(
00586 TransitionList("ExtraLymanLines",&iso_sp[ipISO][nelem].st));
00587 }
00588 else
00589 {
00590 Transitions[ipISO].push_back(
00591 TransitionList("Insanity",&AnonStates));
00592 SatelliteLines[ipISO].push_back(
00593 TransitionList("Insanity",&AnonStates));
00594 ExtraLymanLines[ipISO].push_back(
00595 TransitionList("Insanity",&AnonStates));
00596 }
00597 }
00598 for( long nelem=ipISO; nelem < LIMELM; ++nelem )
00599 {
00600
00601 if( nelem < 2 || dense.lgElmtOn[nelem] )
00602 {
00603 SatelliteLines[ipISO][nelem].resize( iso_sp[ipISO][nelem].numLevels_max );
00604 AllTransitions.push_back(SatelliteLines[ipISO][nelem]);
00605 unsigned int nLine = 0;
00606 for( long ipLo=0; ipLo<iso_sp[ipISO][nelem].numLevels_max; ipLo++ )
00607 {
00608
00609
00610 ipSatelliteLines[ipISO][nelem][ipLo] = nLine;
00611 SatelliteLines[ipISO][nelem][nLine].Junk();
00612 long ipHi = iso_sp[ipISO][nelem].numLevels_max;
00613 SatelliteLines[ipISO][nelem][nLine].setHi(ipHi);
00614 SatelliteLines[ipISO][nelem][nLine].setLo(ipLo);
00615 SatelliteLines[ipISO][nelem][nLine].AddLine2Stack();
00616 ++nLine;
00617 }
00618 ASSERT(SatelliteLines[ipISO][nelem].size() == nLine);
00619
00620
00621
00622 Transitions[ipISO][nelem].resize( iso_sp[ipISO][nelem].ipTrans.size() );
00623 AllTransitions.push_back(Transitions[ipISO][nelem]);
00624 unsigned int nTransition=0;
00625 for( long ipHi=1; ipHi<iso_sp[ipISO][nelem].numLevels_max; ipHi++ )
00626 {
00627 for( long ipLo=0; ipLo < ipHi; ipLo++ )
00628 {
00629
00630 iso_sp[ipISO][nelem].ipTrans[ipHi][ipLo] = nTransition;
00631 Transitions[ipISO][nelem][nTransition].Junk();
00632 Transitions[ipISO][nelem][nTransition].setHi(ipHi);
00633 Transitions[ipISO][nelem][nTransition].setLo(ipLo);
00634 ++nTransition;
00635 }
00636 }
00637 ASSERT(Transitions[ipISO][nelem].size() == nTransition);
00638 iso_sp[ipISO][nelem].tr = &Transitions[ipISO][nelem];
00639
00640
00641 AllTransitions.push_back(ExtraLymanLines[ipISO][nelem]);
00642 ExtraLymanLines[ipISO][nelem].resize(iso_ctrl.nLyman_malloc[ipISO]-2);
00643 ExtraLymanLines[ipISO][nelem].states() = &iso_sp[ipISO][nelem].st;
00644 unsigned int nExtraLyman = 0;
00645 for( long ipHi=2; ipHi < iso_ctrl.nLyman_malloc[ipISO]; ipHi++ )
00646 {
00647 ipExtraLymanLines[ipISO][nelem][ipHi] = nExtraLyman;
00648 ExtraLymanLines[ipISO][nelem][nExtraLyman].Junk();
00649 long ipHi_offset = iso_sp[ipISO][nelem].numLevels_max + 1 + ipHi - 2;
00650 ExtraLymanLines[ipISO][nelem][nExtraLyman].setHi(ipHi_offset);
00651
00652 ExtraLymanLines[ipISO][nelem][nExtraLyman].setLo(0);
00653 ExtraLymanLines[ipISO][nelem][nExtraLyman].AddLine2Stack();
00654 ++nExtraLyman;
00655 }
00656 ASSERT(ExtraLymanLines[ipISO][nelem].size() == nExtraLyman);
00657 }
00658 }
00659 }
00660
00661 return;
00662 }
00663
00664 STATIC void iso_assign_quantum_numbers(void)
00665 {
00666 long int
00667 ipLo,
00668 level,
00669 i,
00670 in,
00671 il,
00672 is,
00673 ij;
00674
00675 DEBUG_ENTRY( "iso_assign_quantum_numbers()" );
00676
00677 for( long nelem=ipHYDROGEN; nelem < LIMELM; nelem++ )
00678 {
00679 long ipISO = ipH_LIKE;
00680
00681 if( nelem == ipHELIUM || dense.lgElmtOn[nelem] )
00682 {
00683 i = 0;
00684
00685
00686 is = ipDOUBLET;
00687
00688
00689 for( in = 1L; in <= iso_sp[ipISO][nelem].n_HighestResolved_max; ++in )
00690 {
00691 for( il = 0L; il < in; ++il )
00692 {
00693 iso_sp[ipISO][nelem].st[i].n() = in;
00694 iso_sp[ipISO][nelem].st[i].S() = is;
00695 iso_sp[ipISO][nelem].st[i].l() = il;
00696 iso_sp[ipISO][nelem].st[i].j() = -1;
00697 iso_sp[ipISO][nelem].QuantumNumbers2Index[in][il][is] = i;
00698 ++i;
00699 }
00700 }
00701
00702 in = iso_sp[ipISO][nelem].n_HighestResolved_max + 1;
00703 for( level = i; level< iso_sp[ipISO][nelem].numLevels_max; ++level)
00704 {
00705 iso_sp[ipISO][nelem].st[level].n() = in;
00706 iso_sp[ipISO][nelem].st[level].S() = -LONG_MAX;
00707 iso_sp[ipISO][nelem].st[level].l() = -LONG_MAX;
00708 iso_sp[ipISO][nelem].st[level].j() = -1;
00709
00710 for( il = 0; il < in; ++il )
00711 {
00712 iso_sp[ipISO][nelem].QuantumNumbers2Index[in][il][is] = level;
00713 }
00714 ++in;
00715 }
00716 --in;
00717
00718
00719 ASSERT( i <= iso_sp[ipISO][nelem].numLevels_max );
00720
00721
00722 ASSERT( (in > 0) && (in < (iso_sp[ipISO][nelem].n_HighestResolved_max + iso_sp[ipISO][nelem].nCollapsed_max + 1) ) );
00723
00724
00725 for( in = 2L; in <= iso_sp[ipISO][nelem].n_HighestResolved_max + iso_sp[ipISO][nelem].nCollapsed_max; ++in )
00726 {
00727 for( il = 0L; il < in; ++il )
00728 {
00729 ASSERT( iso_sp[ipISO][nelem].st[ iso_sp[ipISO][nelem].QuantumNumbers2Index[in][il][is] ].n() == in );
00730 if( in <= iso_sp[ipISO][nelem].n_HighestResolved_max )
00731 {
00732
00733
00734 ASSERT( iso_sp[ipISO][nelem].st[ iso_sp[ipISO][nelem].QuantumNumbers2Index[in][il][is] ].l() == il );
00735 ASSERT( iso_sp[ipISO][nelem].st[ iso_sp[ipISO][nelem].QuantumNumbers2Index[in][il][is] ].S() == is );
00736 }
00737 }
00738 }
00739 }
00740 }
00741
00742
00743 for( long nelem=ipHELIUM; nelem < LIMELM; nelem++ )
00744 {
00745 long ipISO = ipHE_LIKE;
00746
00747 if( nelem == ipHELIUM || dense.lgElmtOn[nelem] )
00748 {
00749 i = 0;
00750
00751
00752 for( in = 1L; in <= iso_sp[ipISO][nelem].n_HighestResolved_max; ++in )
00753 {
00754 for( il = 0L; il < in; ++il )
00755 {
00756 for( is = 3L; is >= 1L; is -= 2 )
00757 {
00758
00759
00760
00761 if( (il == 1L) && (is == 1L) )
00762 continue;
00763
00764 if( (in == 1L) && (is == 3L) )
00765 continue;
00766
00767
00768 if( is == 1 )
00769 {
00770 iso_sp[ipISO][nelem].st[i].n() = in;
00771 iso_sp[ipISO][nelem].st[i].S() = is;
00772 iso_sp[ipISO][nelem].st[i].l() = il;
00773
00774 iso_sp[ipISO][nelem].st[i].j() = il;
00775 iso_sp[ipISO][nelem].QuantumNumbers2Index[in][il][is] = i;
00776 ++i;
00777 }
00778
00779 else if( (in == 2) && (il == 1) && (is == 3) )
00780 {
00781 ij = 0;
00782 do
00783 {
00784 iso_sp[ipISO][nelem].st[i].n() = in;
00785 iso_sp[ipISO][nelem].st[i].S() = is;
00786 iso_sp[ipISO][nelem].st[i].l() = il;
00787 iso_sp[ipISO][nelem].st[i].j() = ij;
00788 iso_sp[ipISO][nelem].QuantumNumbers2Index[in][il][is] = i;
00789 ++i;
00790 ++ij;
00791
00792 } while ( ij < 3 );
00793 }
00794 else
00795 {
00796 iso_sp[ipISO][nelem].st[i].n() = in;
00797 iso_sp[ipISO][nelem].st[i].S() = is;
00798 iso_sp[ipISO][nelem].st[i].l() = il;
00799 iso_sp[ipISO][nelem].st[i].j() = -1L;
00800 iso_sp[ipISO][nelem].QuantumNumbers2Index[in][il][is] = i;
00801 ++i;
00802 }
00803 }
00804 }
00805
00806 if( in > 1L )
00807 {
00808 iso_sp[ipISO][nelem].st[i].n() = in;
00809 iso_sp[ipISO][nelem].st[i].S() = 1L;
00810 iso_sp[ipISO][nelem].st[i].l() = 1L;
00811 iso_sp[ipISO][nelem].st[i].j() = 1L;
00812 iso_sp[ipISO][nelem].QuantumNumbers2Index[in][1][1] = i;
00813 ++i;
00814 }
00815 }
00816
00817 in = iso_sp[ipISO][nelem].n_HighestResolved_max + 1;
00818 for( level = i; level< iso_sp[ipISO][nelem].numLevels_max; ++level)
00819 {
00820 iso_sp[ipISO][nelem].st[level].n() = in;
00821 iso_sp[ipISO][nelem].st[level].S() = -LONG_MAX;
00822 iso_sp[ipISO][nelem].st[level].l() = -LONG_MAX;
00823 iso_sp[ipISO][nelem].st[level].j() = -1;
00824
00825 for( il = 0; il < in; ++il )
00826 {
00827 for( is = 1; is <= 3; is += 2 )
00828 {
00829 iso_sp[ipISO][nelem].QuantumNumbers2Index[in][il][is] = level;
00830 }
00831 }
00832 ++in;
00833 }
00834 --in;
00835
00836
00837 ASSERT( i <= iso_sp[ipISO][nelem].numLevels_max );
00838
00839
00840 ASSERT( (in > 0) && (in < (iso_sp[ipISO][nelem].n_HighestResolved_max + iso_sp[ipISO][nelem].nCollapsed_max + 1) ) );
00841
00842
00843 for( in = 2L; in <= iso_sp[ipISO][nelem].n_HighestResolved_max + iso_sp[ipISO][nelem].nCollapsed_max; ++in )
00844 {
00845 for( il = 0L; il < in; ++il )
00846 {
00847 for( is = 3L; is >= 1; is -= 2 )
00848 {
00849
00850 if( (in == 1L) && (is == 3L) )
00851 continue;
00852
00853 ASSERT( iso_sp[ipISO][nelem].st[ iso_sp[ipISO][nelem].QuantumNumbers2Index[in][il][is] ].n() == in );
00854 if( in <= iso_sp[ipISO][nelem].n_HighestResolved_max )
00855 {
00856
00857
00858 ASSERT( iso_sp[ipISO][nelem].st[ iso_sp[ipISO][nelem].QuantumNumbers2Index[in][il][is] ].l() == il );
00859 ASSERT( iso_sp[ipISO][nelem].st[ iso_sp[ipISO][nelem].QuantumNumbers2Index[in][il][is] ].S() == is );
00860 }
00861 }
00862 }
00863 }
00864 }
00865 }
00866
00867 for( long ipISO=ipH_LIKE; ipISO<NISO; ipISO++ )
00868 {
00869 for( long nelem=ipISO; nelem < LIMELM; nelem++ )
00870 {
00871
00872 if( nelem < 2 || dense.lgElmtOn[nelem] )
00873 {
00874 for( ipLo=ipH1s; ipLo < iso_sp[ipISO][nelem].numLevels_max; ipLo++ )
00875 {
00876 iso_sp[ipISO][nelem].st[ipLo].nelem() = (int)(nelem+1);
00877 iso_sp[ipISO][nelem].st[ipLo].IonStg() = (int)(nelem+1-ipISO);
00878
00879 if( iso_sp[ipISO][nelem].st[ipLo].j() >= 0 )
00880 {
00881 iso_sp[ipISO][nelem].st[ipLo].g() = 2.f*iso_sp[ipISO][nelem].st[ipLo].j()+1.f;
00882 }
00883 else if( iso_sp[ipISO][nelem].st[ipLo].l() >= 0 )
00884 {
00885 iso_sp[ipISO][nelem].st[ipLo].g() = (2.f*iso_sp[ipISO][nelem].st[ipLo].l()+1.f) *
00886 iso_sp[ipISO][nelem].st[ipLo].S();
00887 }
00888 else
00889 {
00890 if( ipISO == ipH_LIKE )
00891 iso_sp[ipISO][nelem].st[ipLo].g() = 2.f*(realnum)POW2( iso_sp[ipISO][nelem].st[ipLo].n() );
00892 else if( ipISO == ipHE_LIKE )
00893 iso_sp[ipISO][nelem].st[ipLo].g() = 4.f*(realnum)POW2( iso_sp[ipISO][nelem].st[ipLo].n() );
00894 else
00895 {
00896
00897 TotalInsanity();
00898 }
00899 }
00900 char chConfiguration[11] = " ";
00901 long nCharactersWritten = 0;
00902
00903 ASSERT( iso_sp[ipISO][nelem].st[ipLo].n() < 1000 );
00904
00905
00906 if( iso_sp[ipISO][nelem].st[ipLo].n() > iso_sp[ipISO][nelem].n_HighestResolved_max )
00907 {
00908 nCharactersWritten = sprintf( chConfiguration, "n=%3li",
00909 iso_sp[ipISO][nelem].st[ipLo].n() );
00910 }
00911 else if( iso_sp[ipISO][nelem].st[ipLo].j() > 0 )
00912 {
00913 nCharactersWritten = sprintf( chConfiguration, "%3li^%li%c_%li",
00914 iso_sp[ipISO][nelem].st[ipLo].n(),
00915 iso_sp[ipISO][nelem].st[ipLo].S(),
00916 chL[ MIN2( 20, iso_sp[ipISO][nelem].st[ipLo].l() ) ],
00917 iso_sp[ipISO][nelem].st[ipLo].j() );
00918 }
00919 else
00920 {
00921 nCharactersWritten = sprintf( chConfiguration, "%3li^%li%c",
00922 iso_sp[ipISO][nelem].st[ipLo].n(),
00923 iso_sp[ipISO][nelem].st[ipLo].S(),
00924 chL[ MIN2( 20, iso_sp[ipISO][nelem].st[ipLo].l()) ] );
00925 }
00926
00927 ASSERT( nCharactersWritten <= 10 );
00928 chConfiguration[10] = '\0';
00929
00930 strncpy( iso_sp[ipISO][nelem].st[ipLo].chConfig(), chConfiguration, 10 );
00931 }
00932 }
00933 }
00934 }
00935 return;
00936 }
00937
00938 #if defined(__ICC) && defined(__i386)
00939 #pragma optimization_level 1
00940 #endif
00941 STATIC void FillExtraLymanLine( const TransitionList::iterator& t, long ipISO, long nelem, long nHi )
00942 {
00943 double Enerwn, Aul;
00944
00945 DEBUG_ENTRY( "FillExtraLymanLine()" );
00946
00947
00948 (*(*t).Hi()).nelem() = (int)(nelem+1);
00949 (*(*t).Hi()).IonStg() = (int)(nelem+1-ipISO);
00950
00951 (*(*t).Hi()).n() = nHi;
00952
00953
00954 (*(*t).Hi()).g() = iso_sp[ipISO][nelem].st[iso_ctrl.nLyaLevel[ipISO]].g();
00955
00956
00957 Enerwn = iso_sp[ipISO][nelem].fb[0].xIsoLevNIonRyd * RYD_INF * ( 1. - 1./POW2((double)nHi) );
00958
00959
00960 (*t).EnergyWN() = (realnum)(Enerwn);
00961 (*t).WLAng() = (realnum)(1.0e8/ Enerwn/ RefIndex(Enerwn));
00962
00963 if( ipISO == ipH_LIKE )
00964 {
00965 Aul = H_Einstein_A( nHi, 1, 1, 0, nelem+1 );
00966 }
00967 else
00968 {
00969 if( nelem == ipHELIUM )
00970 {
00971
00972 Aul = (1.508e10) / pow((double)nHi,2.975);
00973 }
00974 else
00975 {
00976
00977
00978
00979
00980 Aul = 1.375E10 * pow((double)nelem, 3.9) / pow((double)nHi,3.1);
00981 }
00982 }
00983
00984 (*t).Emis().Aul() = (realnum)Aul;
00985
00986 (*(*t).Hi()).lifetime() = iso_state_lifetime( ipISO, nelem, nHi, 1 );
00987
00988 (*t).Emis().dampXvel() = (realnum)( 1.f / (*(*t).Hi()).lifetime() / PI4 / (*t).EnergyWN() );
00989
00990 (*t).Emis().iRedisFun() = iso_ctrl.ipResoRedist[ipISO];
00991
00992 (*t).Emis().gf() = (realnum)(GetGF((*t).Emis().Aul(), (*t).EnergyWN(), (*(*t).Hi()).g()));
00993
00994
00995 (*t).Emis().opacity() = (realnum)(abscf((*t).Emis().gf(), (*t).EnergyWN(), (*(*t).Lo()).g()));
00996
00997
00998 (*t).ipCont() = INT_MIN;
00999 (*t).Emis().ipFine() = INT_MIN;
01000
01001 {
01002
01003
01004
01005 enum {DEBUG_LOC=false};
01006 if( DEBUG_LOC )
01007 {
01008 fprintf(ioQQQ,"%li\t%li\t%.2e\t%.2e\n",
01009 nelem+1,
01010 nHi,
01011 (*t).Emis().Aul() ,
01012 (*t).Emis().opacity()
01013 );
01014 }
01015 }
01016 return;
01017 }
01018
01019
01020 double iso_state_lifetime( long ipISO, long nelem, long n, long l )
01021 {
01022
01023
01024 double tau, t0, eps2;
01025
01026 double m = ELECTRON_MASS;
01027
01028 double M = (double)dense.AtomicWeight[nelem] * ATOMIC_MASS_UNIT;
01029 double mu = (m*M)/(M+m);
01030 long z = 1;
01031 long Z = nelem + 1 - ipISO;
01032
01033 DEBUG_ENTRY( "iso_state_lifetime()" );
01034
01035
01036 ASSERT( l > 0 );
01037
01038 eps2 = 1. - ( l*l + l + 8./47. - (l+1.)/69./n ) / POW2( (double)n );
01039
01040 t0 = 3. * H_BAR * pow( (double)n, 5.) /
01041 ( 2. * POW4( (double)( z * Z ) ) * pow( FINE_STRUCTURE, 5. ) * mu * POW2( SPEEDLIGHT ) ) *
01042 POW2( (m + M)/(Z*m + z*M) );
01043
01044 tau = t0 * ( 1. - eps2 ) /
01045 ( 1. + 19./88.*( (1./eps2 - 1.) * log( 1. - eps2 ) + 1. -
01046 0.5 * eps2 - 0.025 * eps2 * eps2 ) );
01047
01048 if( ipISO == ipHE_LIKE )
01049 {
01050
01051 tau /= 3.;
01052
01053 tau *= 1.1722 * pow( (double)nelem, 0.1 );
01054 }
01055
01056
01057 ASSERT( ipISO <= ipHE_LIKE );
01058 ASSERT( tau > 0. );
01059
01060 return tau;
01061 }
01062
01063
01064 void iso_cascade( long ipISO, long nelem )
01065 {
01066
01067
01068 double *SumAPerN;
01069
01070 long int i, j, ipLo, ipHi;
01071
01072 DEBUG_ENTRY( "iso_cascade()" );
01073
01074 SumAPerN = ((double*)MALLOC( (size_t)(iso_sp[ipISO][nelem].n_HighestResolved_max + iso_sp[ipISO][nelem].nCollapsed_max + 1 )*sizeof(double )));
01075 memset( SumAPerN, 0, (iso_sp[ipISO][nelem].n_HighestResolved_max + iso_sp[ipISO][nelem].nCollapsed_max + 1 )*sizeof(double ) );
01076
01077
01078 iso_sp[ipISO][nelem].CascadeProb[0][0] = 1.;
01079 if( iso_ctrl.lgRandErrGen[ipISO] )
01080 {
01081 iso_sp[ipISO][nelem].fb[0].SigmaAtot = 0.;
01082 iso_sp[ipISO][nelem].ex[0][0].SigmaCascadeProb = 0.;
01083 }
01084
01085
01086
01087
01088 for( ipHi=1; ipHi<iso_sp[ipISO][nelem].numLevels_max; ipHi++ )
01089 {
01090 double SumAs = 0.;
01091
01097
01098 iso_sp[ipISO][nelem].CascadeProb[ipHi][ipHi] = 1.;
01099 iso_sp[ipISO][nelem].CascadeProb[ipHi][0] = 0.;
01100 iso_sp[ipISO][nelem].BranchRatio[ipHi][0] = 0.;
01101
01102 if( iso_ctrl.lgRandErrGen[ipISO] )
01103 {
01104 iso_sp[ipISO][nelem].fb[ipHi].SigmaAtot = 0.;
01105 iso_sp[ipISO][nelem].ex[ipHi][ipHi].SigmaCascadeProb = 0.;
01106 }
01107
01108 long ipLoStart = 0;
01109 if( opac.lgCaseB && L_(ipHi)==1 && (ipISO==ipH_LIKE || S_(ipHi)==1) )
01110 ipLoStart = 1;
01111
01112 for( ipLo=ipLoStart; ipLo<ipHi; ipLo++ )
01113 {
01114 SumAs += iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().Aul();
01115 }
01116
01117 for( ipLo=ipLoStart; ipLo<ipHi; ipLo++ )
01118 {
01119 if( iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().Aul() <= iso_ctrl.SmallA )
01120 {
01121 iso_sp[ipISO][nelem].CascadeProb[ipHi][ipLo] = 0.;
01122 iso_sp[ipISO][nelem].BranchRatio[ipHi][ipLo] = 0.;
01123 continue;
01124 }
01125
01126 iso_sp[ipISO][nelem].CascadeProb[ipHi][ipLo] = 0.;
01127 iso_sp[ipISO][nelem].BranchRatio[ipHi][ipLo] =
01128 iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().Aul() / SumAs;
01129
01130 ASSERT( iso_sp[ipISO][nelem].BranchRatio[ipHi][ipLo] <= 1.0000001 );
01131
01132 SumAPerN[N_(ipHi)] += iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().Aul();
01133
01134
01135
01136
01137 ASSERT( iso_sp[ipISO][nelem].trans(ipHi,ipLo).EnergyWN() > 0. ||
01138 iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().Aul() <= iso_ctrl.SmallA );
01139
01140 if( iso_ctrl.lgRandErrGen[ipISO] )
01141 {
01142 ASSERT( iso_sp[ipISO][nelem].ex[ipHi][ipLo].Error[IPRAD] >= 0. );
01143
01144 iso_sp[ipISO][nelem].fb[ipHi].SigmaAtot +=
01145 pow( iso_sp[ipISO][nelem].ex[ipHi][ipLo].Error[IPRAD] *
01146 (double)iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().Aul(), 2. );
01147 }
01148 }
01149
01150 if( iso_ctrl.lgRandErrGen[ipISO] )
01151 {
01152
01153 iso_sp[ipISO][nelem].fb[ipHi].SigmaAtot = sqrt( iso_sp[ipISO][nelem].fb[ipHi].SigmaAtot );
01154 }
01155
01156
01157 for( i=0; i<ipHi; i++ )
01158 {
01159 for( ipLo=0; ipLo<=i; ipLo++ )
01160 {
01161 iso_sp[ipISO][nelem].CascadeProb[ipHi][ipLo] += iso_sp[ipISO][nelem].BranchRatio[ipHi][i] * iso_sp[ipISO][nelem].CascadeProb[i][ipLo];
01162 }
01163 }
01164
01165 if( iso_ctrl.lgRandErrGen[ipISO] )
01166 {
01167 for( ipLo=0; ipLo<ipHi; ipLo++ )
01168 {
01169 double SigmaCul = 0.;
01170 for( i=ipLo; i<ipHi; i++ )
01171 {
01172 if( iso_sp[ipISO][nelem].trans(ipHi,i).Emis().Aul() > iso_ctrl.SmallA )
01173 {
01174
01175 double SigmaA = iso_sp[ipISO][nelem].ex[ipHi][i].Error[IPRAD] *
01176 iso_sp[ipISO][nelem].trans(ipHi,i).Emis().Aul();
01177 SigmaCul +=
01178 pow(SigmaA*iso_sp[ipISO][nelem].CascadeProb[i][ipLo]*iso_sp[ipISO][nelem].st[ipHi].lifetime(), 2.) +
01179 pow(iso_sp[ipISO][nelem].fb[ipHi].SigmaAtot*iso_sp[ipISO][nelem].BranchRatio[ipHi][i]*
01180 iso_sp[ipISO][nelem].CascadeProb[i][ipLo]*iso_sp[ipISO][nelem].st[ipHi].lifetime(), 2.) +
01181 pow(iso_sp[ipISO][nelem].ex[i][ipLo].SigmaCascadeProb*iso_sp[ipISO][nelem].BranchRatio[ipHi][i], 2.);
01182 }
01183 }
01184 SigmaCul = sqrt(SigmaCul);
01185 iso_sp[ipISO][nelem].ex[ipHi][ipLo].SigmaCascadeProb = SigmaCul;
01186 }
01187 }
01188 }
01189
01190
01191
01192
01193 {
01194 enum {DEBUG_LOC=false};
01195
01196 if( DEBUG_LOC && (nelem == ipHELIUM) && (ipISO==ipHE_LIKE) )
01197 {
01198
01199 long int hi_l,hi_s;
01200 double Bm;
01201
01202
01203
01204 hi_s = -100000;
01205 hi_l = -100000;
01206 ipLo = -100000;
01207
01208
01209
01210
01211
01212
01213
01214
01215
01216 ASSERT( hi_l != iso_sp[ipISO][nelem].st[ipLo].l() );
01217
01218 fprintf(ioQQQ,"Bm(n,%ld,%ld;%ld)\n",hi_l,hi_s,ipLo);
01219 fprintf(ioQQQ,"m\t2\t\t3\t\t4\t\t5\t\t6\n");
01220
01221 for( ipHi=ipHe2p3P2; ipHi<iso_sp[ipISO][nelem].numLevels_max-iso_sp[ipISO][nelem].nCollapsed_max; ipHi++ )
01222 {
01223
01224 if( (iso_sp[ipISO][nelem].st[ipHi].l() == 1) && (iso_sp[ipISO][nelem].st[ipHi].S() == 3) )
01225 {
01226 fprintf(ioQQQ,"\n%ld\t",iso_sp[ipISO][nelem].st[ipHi].n());
01227 j = 0;
01228 Bm = 0;
01229 for( i = ipLo; i<=ipHi; i++)
01230 {
01231 if( (iso_sp[ipISO][nelem].st[i].l() == hi_l) && (iso_sp[ipISO][nelem].st[i].S() == hi_s) )
01232 {
01233 if( (ipLo == ipHe2p3P0) && (i > ipHe2p3P2) )
01234 {
01235 Bm += iso_sp[ipISO][nelem].CascadeProb[ipHi][i] * ( iso_sp[ipISO][nelem].BranchRatio[i][ipHe2p3P0] +
01236 iso_sp[ipISO][nelem].BranchRatio[i][ipHe2p3P1] + iso_sp[ipISO][nelem].BranchRatio[i][ipHe2p3P2] );
01237 }
01238 else
01239 Bm += iso_sp[ipISO][nelem].CascadeProb[ipHi][i] * iso_sp[ipISO][nelem].BranchRatio[i][ipLo];
01240
01241 if( (i == ipHe2p3P0) || (i == ipHe2p3P1) || (i == ipHe2p3P2) )
01242 {
01243 j++;
01244 if(j == 3)
01245 {
01246 fprintf(ioQQQ,"%2.4e\t",Bm);
01247 Bm = 0;
01248 }
01249 }
01250 else
01251 {
01252 fprintf(ioQQQ,"%2.4e\t",Bm);
01253 Bm = 0;
01254 }
01255 }
01256 }
01257 }
01258 }
01259 fprintf(ioQQQ,"\n\n");
01260 }
01261 }
01262
01263
01264
01265
01266
01267 for( i=2; i < iso_sp[ipISO][nelem].n_HighestResolved_max; ++i)
01268 {
01269 ASSERT( (SumAPerN[i] > SumAPerN[i+1]) || opac.lgCaseB );
01270 }
01271
01272 {
01273 enum {DEBUG_LOC=false};
01274 if( DEBUG_LOC )
01275 {
01276 for( i = 2; i<= (iso_sp[ipISO][nelem].n_HighestResolved_max + iso_sp[ipISO][nelem].nCollapsed_max); ++i)
01277 {
01278 fprintf(ioQQQ,"n %ld\t lifetime %.4e\n", i, 1./SumAPerN[i]);
01279 }
01280 }
01281 }
01282
01283 free( SumAPerN );
01284
01285 return;
01286 }
01287
01289
01290
01291 STATIC void iso_satellite( void )
01292 {
01293 DEBUG_ENTRY( "iso_satellite()" );
01294
01295 for( long ipISO = ipHE_LIKE; ipISO < NISO; ipISO++ )
01296 {
01297 for( long nelem = ipISO; nelem < LIMELM; nelem++ )
01298 {
01299
01300 if( nelem == ipISO || dense.lgElmtOn[nelem] )
01301 {
01302 for( long i=0; i<iso_sp[ipISO][nelem].numLevels_max; i++ )
01303 {
01304 char chLab[5]=" ";
01305
01306 TransitionList::iterator tr = SatelliteLines[ipISO][nelem].begin()+ipSatelliteLines[ipISO][nelem][i];
01307 (*tr).Zero();
01308
01309
01310
01311
01312
01313 (*tr).WLAng() = (realnum)(RYDLAM/
01314 ((iso_sp[ipISO-1][nelem].fb[0].xIsoLevNIonRyd - iso_sp[ipISO-1][nelem].fb[1].xIsoLevNIonRyd) -
01315 (iso_sp[ipISO][nelem].fb[1].xIsoLevNIonRyd- iso_sp[ipISO][nelem].fb[i].xIsoLevNIonRyd)) );
01316
01317 (*tr).EnergyWN() = 1.e8f /
01318 (*tr).WLAng();
01319
01320
01321 sprintf( chLab, "%2s%2ld",elementnames.chElementSym[nelem], nelem+1-ipISO );
01322
01323 (*tr).Emis().iRedisFun() = ipCRDW;
01324
01325 (*(*tr).Hi()).nelem() = nelem + 1;
01326 (*(*tr).Hi()).IonStg() = nelem + 1 - ipISO;
01327 fixit();
01328 (*(*tr).Hi()).g() = 2.f;
01329 (*(*tr).Lo()).g() = iso_sp[ipISO][nelem].st[i].g();
01330 (*tr).Emis().PopOpc() =
01331 (*(*tr).Lo()).Pop();
01332
01333 (*tr).Emis().pump() = 0.;
01334
01335 }
01336 }
01337 }
01338 }
01339
01340 return;
01341 }
01342
01343 void iso_satellite_update( long nelem )
01344 {
01345 double ConBoltz, LTE_pop=SMALLFLOAT+FLT_EPSILON, factor1, ConvLTEPOP;
01346
01347 DEBUG_ENTRY( "iso_satellite_update()" );
01348
01349 for( long ipISO = ipHE_LIKE; ipISO < MIN2(NISO,nelem+1); ipISO++ )
01350 {
01351
01352 if( nelem == ipISO || dense.lgElmtOn[nelem] )
01353 {
01354 for( long i=0; i<iso_sp[ipISO][nelem].numLevels_max; i++ )
01355 {
01356 double dr_rate = iso_sp[ipISO][nelem].fb[i].DielecRecomb * iso_ctrl.lgDielRecom[ipISO];
01357
01358 TransitionList::iterator tr = SatelliteLines[ipISO][nelem].begin()+ipSatelliteLines[ipISO][nelem][i];
01359 (*tr).Emis().phots() =
01360 dr_rate * dense.eden * dense.xIonDense[nelem][nelem+1-ipISO];
01361
01362 (*tr).Emis().xIntensity() =
01363 (*tr).Emis().phots() *
01364 ERG1CM * (*tr).EnergyWN();
01365
01366
01367
01368
01369 factor1 = HION_LTE_POP*dense.AtomicWeight[nelem]/
01370 (dense.AtomicWeight[nelem]+ELECTRON_MASS/ATOMIC_MASS_UNIT);
01371
01372
01373 ConvLTEPOP = pow(factor1,1.5)/(2.*iso_ctrl.stat_ion[ipISO])/phycon.te32;
01374
01375
01376
01377
01378 ConBoltz = dsexp(iso_sp[ipISO-1][nelem].fb[1].xIsoLevNIonRyd/phycon.te_ryd);
01379
01380 if( ConBoltz >= SMALLDOUBLE )
01381 {
01382
01383
01384
01385
01386 LTE_pop = (*(*tr).Hi()).g() * ConBoltz * ConvLTEPOP;
01387 }
01388
01389 LTE_pop = max( LTE_pop, 1e-30f );
01390
01391
01392 (*tr).Emis().Aul() = (realnum)(dr_rate/LTE_pop);
01393 (*tr).Emis().Aul() =
01394 max( iso_ctrl.SmallA, (*tr).Emis().Aul() );
01395
01396 (*tr).Emis().gf() = (realnum)GetGF(
01397 (*tr).Emis().Aul(),
01398 (*tr).EnergyWN(),
01399 (*(*tr).Hi()).g());
01400
01401 (*tr).Emis().gf() =
01402 max( 1e-20f, (*tr).Emis().gf() );
01403
01404 (*(*tr).Hi()).Pop() = LTE_pop * dense.xIonDense[nelem][nelem+1-ipISO] * dense.eden;
01405
01406 (*tr).Emis().PopOpc() =
01407 (*(*tr).Lo()).Pop() -
01408 (*(*tr).Hi()).Pop() *
01409 (*(*tr).Lo()).g()/(*(*tr).Hi()).g();
01410
01411 (*tr).Emis().opacity() =
01412 (realnum)(abscf((*tr).Emis().gf(),
01413 (*tr).EnergyWN(),
01414 (*(*tr).Lo()).g()));
01415
01416
01417 double lifetime = 1e-10;
01418
01419 (*tr).Emis().dampXvel() = (realnum)(
01420 (1.f/lifetime)/PI4/(*tr).EnergyWN());
01421 }
01422 }
01423 }
01424
01425 return;
01426 }
01427
01428 long iso_get_total_num_levels( long ipISO, long nmaxResolved, long numCollapsed )
01429 {
01430 DEBUG_ENTRY( "iso_get_total_num_levels()" );
01431
01432 long tot_num_levels;
01433
01434
01435
01436
01437 if( ipISO == ipH_LIKE )
01438 {
01439 tot_num_levels = (long)( nmaxResolved * 0.5 *( nmaxResolved + 1 ) ) + numCollapsed;
01440 }
01441 else if( ipISO == ipHE_LIKE )
01442 {
01443 tot_num_levels = nmaxResolved*nmaxResolved + nmaxResolved + 1 + numCollapsed;
01444 }
01445 else
01446 TotalInsanity();
01447
01448 return tot_num_levels;
01449 }
01450
01451 void iso_update_num_levels( long ipISO, long nelem )
01452 {
01453 DEBUG_ENTRY( "iso_update_num_levels()" );
01454
01455
01456 ASSERT( iso_sp[ipISO][nelem].n_HighestResolved_max >= 3 );
01457
01458 iso_sp[ipISO][nelem].numLevels_max =
01459 iso_get_total_num_levels( ipISO, iso_sp[ipISO][nelem].n_HighestResolved_max, iso_sp[ipISO][nelem].nCollapsed_max );
01460
01461 if( iso_sp[ipISO][nelem].numLevels_max > iso_sp[ipISO][nelem].numLevels_malloc )
01462 {
01463 fprintf( ioQQQ, "The number of levels for ipISO %li, nelem %li, has been increased since the initial coreload.\n",
01464 ipISO, nelem );
01465 fprintf( ioQQQ, "This cannot be done.\n" );
01466 cdEXIT(EXIT_FAILURE);
01467 }
01468
01469
01470 iso_sp[ipISO][nelem].numLevels_local = iso_sp[ipISO][nelem].numLevels_max;
01471 iso_sp[ipISO][nelem].nCollapsed_local = iso_sp[ipISO][nelem].nCollapsed_max;
01472 iso_sp[ipISO][nelem].n_HighestResolved_local = iso_sp[ipISO][nelem].n_HighestResolved_max;
01473
01474
01475
01476
01477 max_num_levels = MAX2( max_num_levels, iso_sp[ipISO][nelem].numLevels_max);
01478
01479 return;
01480 }
01481
01482 void iso_collapsed_bnl_set( long ipISO, long nelem )
01483 {
01484
01485 DEBUG_ENTRY( "iso_collapsed_bnl_set()" );
01486
01487 double bnl_array[4][3][4][10] = {
01488 {
01489
01490 {
01491 {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},
01492 {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},
01493 {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},
01494 {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}
01495 },
01496 {
01497 {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},
01498 {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},
01499 {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},
01500 {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}
01501 },
01502 {
01503 {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},
01504 {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},
01505 {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},
01506 {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}
01507 }
01508 },
01509 {
01510
01511 {
01512 {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},
01513 {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},
01514 {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},
01515 {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}
01516 },
01517 {
01518 {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},
01519 {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},
01520 {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},
01521 {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}
01522 },
01523 {
01524 {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},
01525 {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},
01526 {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},
01527 {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}
01528 }
01529 },
01530 {
01531
01532 {
01533 {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},
01534 {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},
01535 {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},
01536 {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}
01537 },
01538 {
01539 {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},
01540 {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},
01541 {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},
01542 {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}
01543 },
01544 {
01545 {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},
01546 {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},
01547 {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},
01548 {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}
01549 }
01550 },
01551 {
01552
01553 {
01554 {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},
01555 {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},
01556 {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},
01557 {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}
01558 },
01559 {
01560 {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},
01561 {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},
01562 {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},
01563 {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}
01564 },
01565 {
01566 {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},
01567 {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},
01568 {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},
01569 {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}
01570 }
01571 }
01572 };
01573
01574 double temps[4] = {5000., 10000., 15000., 20000. };
01575 double log_dens[3] = {2., 4., 6.};
01576 long ipTe, ipDens;
01577
01578 ASSERT( nelem <= 1 );
01579
01580
01581 ipTe = hunt_bisect( temps, 4, phycon.te );
01582 ipDens = hunt_bisect( log_dens, 3, log10(dense.eden) );
01583
01584 ASSERT( (ipTe >=0) && (ipTe < 3) );
01585 ASSERT( (ipDens >=0) && (ipDens < 2) );
01586
01587 for( long nHi=iso_sp[ipISO][nelem].n_HighestResolved_max+1; nHi<=iso_sp[ipISO][nelem].n_HighestResolved_max+iso_sp[ipISO][nelem].nCollapsed_max; nHi++ )
01588 {
01589 for( long lHi=0; lHi<nHi; lHi++ )
01590 {
01591 for( long sHi=1; sHi<4; sHi++ )
01592 {
01593 if( ipISO == ipH_LIKE && sHi != 2 )
01594 continue;
01595 else if( ipISO == ipHE_LIKE && sHi != 1 && sHi != 3 )
01596 continue;
01597
01598 double bnl_at_lo_den, bnl_at_hi_den, bnl;
01599 double bnl_max, bnl_min, temp, dens;
01600
01601 long ipL = MIN2(9,lHi);
01602 long ip1;
01603
01604 if( nelem==ipHYDROGEN )
01605 ip1 = 0;
01606 else if( nelem==ipHELIUM )
01607 {
01608 if( ipISO==ipH_LIKE )
01609 ip1 = 1;
01610 else if( ipISO==ipHE_LIKE )
01611 {
01612 if( sHi==1 )
01613 ip1 = 2;
01614 else if( sHi==3 )
01615 ip1 = 3;
01616 else
01617 TotalInsanity();
01618 }
01619 else
01620 TotalInsanity();
01621 }
01622 else
01623 TotalInsanity();
01624
01625 temp = MAX2( temps[0], phycon.te );
01626 temp = MIN2( temps[3], temp );
01627
01628 dens = MAX2( log_dens[0], log10(dense.eden) );
01629 dens = MIN2( log_dens[2], dens );
01630
01631
01632
01633
01634
01635 if( temp < temps[0] && dens < log_dens[0] )
01636 bnl = bnl_array[ip1][0][0][ipL];
01637 else if( temp < temps[0] && dens >= log_dens[2] )
01638 bnl = bnl_array[ip1][2][0][ipL];
01639 else if( temp >= temps[3] && dens < log_dens[0] )
01640 bnl = bnl_array[ip1][0][3][ipL];
01641 else if( temp >= temps[3] && dens >= log_dens[2] )
01642 bnl = bnl_array[ip1][2][3][ipL];
01643 else
01644 {
01645 bnl_at_lo_den = ( temp - temps[ipTe]) / (temps[ipTe+1] - temps[ipTe]) *
01646 (bnl_array[ip1][ipDens][ipTe+1][ipL] - bnl_array[ip1][ipDens][ipTe][ipL]) + bnl_array[ip1][ipDens][ipTe][ipL];
01647
01648 bnl_at_hi_den = ( temp - temps[ipTe]) / (temps[ipTe+1] - temps[ipTe]) *
01649 (bnl_array[ip1][ipDens+1][ipTe+1][ipL] - bnl_array[ip1][ipDens+1][ipTe][ipL]) + bnl_array[ip1][ipDens+1][ipTe][ipL];
01650
01651 bnl = ( dens - log_dens[ipDens]) / (log_dens[ipDens+1] - log_dens[ipDens]) *
01652 (bnl_at_hi_den - bnl_at_lo_den) + bnl_at_lo_den;
01653 }
01654
01656 {
01657 bnl_max = MAX4( bnl_array[ip1][ipDens][ipTe+1][ipL], bnl_array[ip1][ipDens+1][ipTe+1][ipL],
01658 bnl_array[ip1][ipDens][ipTe][ipL], bnl_array[ip1][ipDens+1][ipTe][ipL] );
01659 ASSERT( bnl <= bnl_max );
01660
01661 bnl_min = MIN4( bnl_array[ip1][ipDens][ipTe+1][ipL], bnl_array[ip1][ipDens+1][ipTe+1][ipL],
01662 bnl_array[ip1][ipDens][ipTe][ipL], bnl_array[ip1][ipDens+1][ipTe][ipL] );
01663 ASSERT( bnl >= bnl_min );
01664 }
01665
01666 iso_sp[ipISO][nelem].bnl_effective[nHi][lHi][sHi] = bnl;
01667
01668 ASSERT( iso_sp[ipISO][nelem].bnl_effective[nHi][lHi][sHi] > 0. );
01669 }
01670 }
01671 }
01672
01673 return;
01674 }
01675
01676
01677 void iso_collapsed_bnl_print( long ipISO, long nelem )
01678 {
01679 DEBUG_ENTRY( "iso_collapsed_bnl_print()" );
01680
01681 for( long is = 1; is<=3; ++is)
01682 {
01683 if( ipISO == ipH_LIKE && is != 2 )
01684 continue;
01685 else if( ipISO == ipHE_LIKE && is != 1 && is != 3 )
01686 continue;
01687
01688 char chSpin[3][9]= {"singlets", "doublets", "triplets"};
01689
01690
01691 fprintf(ioQQQ," %s %s %s bnl\n",
01692 iso_ctrl.chISO[ipISO],
01693 elementnames.chElementSym[nelem],
01694 chSpin[is-1]);
01695
01696
01697 fprintf(ioQQQ," n\\l=> ");
01698 for( long i =0; i < iso_sp[ipISO][nelem].n_HighestResolved_max + iso_sp[ipISO][nelem].nCollapsed_max; ++i)
01699 {
01700 fprintf(ioQQQ,"%2ld ",i);
01701 }
01702 fprintf(ioQQQ,"\n");
01703
01704
01705 for( long in = 1; in <= iso_sp[ipISO][nelem].n_HighestResolved_max + iso_sp[ipISO][nelem].nCollapsed_max; ++in)
01706 {
01707 if( is==3 && in==1 )
01708 continue;
01709
01710 fprintf(ioQQQ," %2ld ",in);
01711
01712 for( long il = 0; il < in; ++il)
01713 {
01714 fprintf( ioQQQ, "%9.3e ", iso_sp[ipISO][nelem].bnl_effective[in][il][is] );
01715 }
01716 fprintf(ioQQQ,"\n");
01717 }
01718 }
01719
01720 return;
01721 }
01722
01723 void iso_collapsed_Aul_update( long ipISO, long nelem )
01724 {
01725 DEBUG_ENTRY( "iso_collapsed_Aul_update()" );
01726
01727 long ipFirstCollapsed = iso_sp[ipISO][nelem].numLevels_max - iso_sp[ipISO][nelem].nCollapsed_max;
01728
01729 for( long ipLo=0; ipLo<ipFirstCollapsed; ipLo++ )
01730 {
01731 long spin = iso_sp[ipISO][nelem].st[ipLo].S();
01732
01733
01734 for( long nHi=iso_sp[ipISO][nelem].n_HighestResolved_max+1; nHi<=iso_sp[ipISO][nelem].n_HighestResolved_max+iso_sp[ipISO][nelem].nCollapsed_max; nHi++ )
01735 {
01736 realnum Auls[2] = {
01737 iso_sp[ipISO][nelem].CachedAs[ nHi-iso_sp[ipISO][nelem].n_HighestResolved_max-1 ][ ipLo ][0],
01738 iso_sp[ipISO][nelem].CachedAs[ nHi-iso_sp[ipISO][nelem].n_HighestResolved_max-1 ][ ipLo ][1] };
01739
01740 realnum EffectiveAul =
01741 Auls[0]*spin*(2.f*(L_(ipLo)+1.f)+1.f)*(realnum)iso_sp[ipISO][nelem].bnl_effective[nHi][ L_(ipLo)+1 ][spin];
01742
01743
01744
01745 if( L_(ipLo) > 0 )
01746 {
01747 EffectiveAul +=
01748 Auls[1]*spin*(2.f*(L_(ipLo)-1.f)+1.f)*(realnum)iso_sp[ipISO][nelem].bnl_effective[nHi][ L_(ipLo)-1 ][spin];
01749 }
01750
01751 if( ipISO==ipH_LIKE )
01752 EffectiveAul /= (2.f*nHi*nHi);
01753 else if( ipISO==ipHE_LIKE )
01754 EffectiveAul /= (4.f*nHi*nHi);
01755 else
01756 TotalInsanity();
01757
01758 long ipHi = iso_sp[ipISO][nelem].QuantumNumbers2Index[nHi][ L_(ipLo)+1 ][spin];
01759
01760
01761 iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().Aul() = EffectiveAul;
01762
01763 ASSERT( iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().Aul() > 0. );
01764 }
01765 }
01766
01767 return;
01768 }
01769
01770 void iso_collapsed_lifetimes_update( long ipISO, long nelem )
01771 {
01772 DEBUG_ENTRY( "iso_collapsed_Aul_update()" );
01773
01774 for( long ipHi=iso_sp[ipISO][nelem].numLevels_max- iso_sp[ipISO][nelem].nCollapsed_max; ipHi < iso_sp[ipISO][nelem].numLevels_max; ipHi++ )
01775 {
01776 iso_sp[ipISO][nelem].st[ipHi].lifetime() = SMALLFLOAT;
01777
01778 for( long ipLo=0; ipLo < ipHi; ipLo++ )
01779 {
01780 if( iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().Aul() <= iso_ctrl.SmallA )
01781 continue;
01782
01783 iso_sp[ipISO][nelem].st[ipHi].lifetime() += iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().Aul();
01784 }
01785
01786
01787 iso_sp[ipISO][nelem].st[ipHi].lifetime() = 1./iso_sp[ipISO][nelem].st[ipHi].lifetime();
01788
01789 for( long ipLo=0; ipLo < ipHi; ipLo++ )
01790 {
01791 if( iso_sp[ipISO][nelem].trans(ipHi,ipLo).EnergyWN() <= 0. )
01792 continue;
01793
01794 if( iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().Aul() <= iso_ctrl.SmallA )
01795 continue;
01796
01797 iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().dampXvel() = (realnum)(
01798 (1.f/iso_sp[ipISO][nelem].st[ipHi].lifetime())/
01799 PI4/iso_sp[ipISO][nelem].trans(ipHi,ipLo).EnergyWN());
01800
01801 ASSERT(iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().dampXvel()> 0.);
01802 }
01803 }
01804
01805 return;
01806 }
01807
01808 #if 0
01809 STATIC void Prt_AGN_table( void )
01810 {
01811
01812 multi_arr<char,2> chLevel(max_num_levels,10);
01813
01814
01815 for( long ipLo=0; ipLo < iso_sp[ipISO][ipISO].numLevels_max-iso_sp[ipISO][ipISO].nCollapsed_max; ++ipLo )
01816 {
01817 long nelem = ipISO;
01818 sprintf( &chLevel.front(ipLo) , "%li %li%c", N_(ipLo), S_(ipLo), chL[MIN2(20,L_(ipLo))] );
01819 }
01820
01821
01822
01823 {
01824
01825 enum {AGN=false};
01826 if( AGN )
01827 {
01828 # define NTEMP 6
01829 double te[NTEMP]={6000.,8000.,10000.,15000.,20000.,25000. };
01830 double telog[NTEMP] ,
01831 cs ,
01832 ratecoef;
01833 long nelem = ipHELIUM;
01834 fprintf( ioQQQ,"trans");
01835 for( long i=0; i < NTEMP; ++i )
01836 {
01837 telog[i] = log10( te[i] );
01838 fprintf( ioQQQ,"\t%.3e",te[i]);
01839 }
01840 for( long i=0; i < NTEMP; ++i )
01841 {
01842 fprintf( ioQQQ,"\t%.3e",te[i]);
01843 }
01844 fprintf(ioQQQ,"\n");
01845
01846 for( long ipHi=ipHe2s3S; ipHi< iso_sp[ipHE_LIKE][ipHELIUM].numLevels_max; ++ipHi )
01847 {
01848 for( long ipLo=ipHe1s1S; ipLo < ipHi; ++ipLo )
01849 {
01850
01851
01852
01853 if( N_(ipHi) == N_(ipLo) )
01854 continue;
01855
01856
01857 fprintf( ioQQQ,"%s - %s",
01858 &chLevel.front(ipLo) , &chLevel.front(ipHi) );
01859
01860
01861 for( long i=0; i < NTEMP; ++i )
01862 {
01863 phycon.alogte = telog[i];
01864
01865 cs = HeCSInterp( nelem , ipHi , ipLo, ipELECTRON );
01866 fprintf(ioQQQ,"\t%.2e", cs );
01867 }
01868
01869
01870 for( long i=0; i < NTEMP; ++i )
01871 {
01872 phycon.alogte = telog[i];
01873 phycon.te = pow(10.,telog[i] );
01874 tfidle(false);
01875 cs = HeCSInterp( nelem , ipHi , ipLo, ipELECTRON );
01876
01877 ratecoef = cs/sqrt(phycon.te)*COLL_CONST/iso_sp[ipHE_LIKE][nelem].st[ipLo].g() *
01878 sexp( iso_sp[ipHE_LIKE][nelem].trans(ipHi,ipLo).EnergyK() / phycon.te );
01879 fprintf(ioQQQ,"\t%.2e", ratecoef );
01880 }
01881 fprintf(ioQQQ,"\n");
01882 }
01883 }
01884 cdEXIT(EXIT_FAILURE);
01885 }
01886 }
01887
01888 return;
01889 }
01890 #endif