00001
00002
00003
00004
00005
00006
00007
00008
00009 #include "cddefines.h"
00010 #define TAULIM 1e8
00011 #include "taulines.h"
00012 #include "doppvel.h"
00013 #include "iso.h"
00014 #include "h2.h"
00015 #include "lines_service.h"
00016 #include "rfield.h"
00017 #include "dense.h"
00018 #include "opacity.h"
00019 #include "thermal.h"
00020 #include "geometry.h"
00021 #include "stopcalc.h"
00022 #include "ipoint.h"
00023 #include "abund.h"
00024 #include "conv.h"
00025 #include "atomfeii.h"
00026 #include "rt.h"
00027 #include "trace.h"
00028
00029 void RT_tau_init(void)
00030 {
00031 long int i,
00032 nelem,
00033 ipISO,
00034 ipHi,
00035 ipLo,
00036 nHi;
00037 long lgHit;
00038
00039 double AbunRatio,
00040 balc,
00041 coleff;
00042
00043 realnum f, BalmerTauOn;
00044
00045 DEBUG_ENTRY( "RT_tau_init()" );
00046
00047 ASSERT( dense.eden > 0. );
00048
00049
00050 for( ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
00051 {
00052 for( nelem=ipISO; nelem < LIMELM; nelem++ )
00053 {
00054 if( dense.lgElmtOn[nelem] )
00055 {
00056
00057 for( ipLo=0; ipLo < iso.numLevels_max[ipISO][nelem]; ipLo++ )
00058 {
00059 if( iso.lgDielRecom[ipISO] )
00060 SatelliteLines[ipISO][nelem][ipLo].Zero();
00061 }
00062
00063 for( ipHi=1; ipHi < iso.numLevels_max[ipISO][nelem]; ipHi++ )
00064 {
00065 for( ipLo=0; ipLo < ipHi; ipLo++ )
00066 {
00067 Transitions[ipISO][nelem][ipHi][ipLo].Zero();
00068 }
00069 }
00070 for( ipHi=2; ipHi <iso.nLyman[ipISO]; ipHi++ )
00071 {
00072 ExtraLymanLines[ipISO][nelem][ipHi].Zero();
00073 }
00074 }
00075 }
00076 }
00077
00078
00079 for( i=1; i <= nLevel1; i++ )
00080 {
00081 TauLines[i].Zero();
00082 }
00083
00084
00085
00086 TauDummy.Zero();
00087
00088
00089 for( i=0; i < nWindLine; i++ )
00090 {
00091 if( TauLine2[i].Hi->IonStg < TauLine2[i].Hi->nelem+1-NISO )
00092 {
00093
00094 TauLine2[i].Zero();
00095 }
00096 }
00097
00098
00099 for( i=0; i < nUTA; i++ )
00100 {
00101
00102
00103
00104 double hsave = UTALines[i].Coll.heat;
00105 UTALines[i].Zero();
00106 UTALines[i].Coll.heat = hsave;
00107 }
00108
00109
00110 for( i=0; i < nHFLines; i++ )
00111 {
00112 HFLines[i].Zero();
00113 }
00114
00115
00116 for( i=0; i < linesAdded2; i++)
00117 {
00118 dBaseLines[i].tran->Zero();
00119 }
00120
00121
00122 H2_LineZero();
00123
00124
00125 FeII_LineZero();
00126
00127
00128
00129
00130
00131
00132
00133 if( StopCalc.taunu > 0. )
00134 {
00135 StopCalc.iptnu = ipoint(StopCalc.taunu);
00136 StopCalc.iptnu = MIN2(StopCalc.iptnu,rfield.nupper-1);
00137 }
00138 else
00139 {
00140 StopCalc.iptnu = rfield.nupper;
00141
00142
00143 }
00144
00145
00146
00147 ASSERT( StopCalc.taunu == 0. || StopCalc.iptnu >= 0 );
00148
00149
00150
00151 if( StopCalc.taunu > 0. )
00152 {
00153
00154 if( StopCalc.iptnu >= iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH1s] )
00155 {
00156
00157 for( i=0; i < (iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH1s] - 1); i++ )
00158 {
00159
00160 opac.TauAbsGeo[1][i] = opac.taumin;
00161 opac.TauScatGeo[1][i] = opac.taumin;
00162 opac.TauTotalGeo[1][i] = opac.taumin;
00163 }
00164
00165 for( i=iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH1s]-1; i < rfield.nupper; i++ )
00166 {
00167
00168 opac.TauAbsGeo[1][i] = StopCalc.tauend;
00169 opac.TauScatGeo[1][i] = opac.taumin;
00170 opac.TauTotalGeo[1][i] = opac.TauAbsGeo[1][i] + opac.TauScatGeo[1][i];
00171 }
00172 }
00173
00174 else
00175 {
00176
00177 for( i=0; i < (iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH1s] - 1); i++ )
00178 {
00179 opac.TauAbsGeo[1][i] = StopCalc.tauend;
00180 opac.TauScatGeo[1][i] = StopCalc.tauend;
00181 opac.TauTotalGeo[1][i] = StopCalc.tauend;
00182 }
00183
00184 for( i=iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH1s]-1; i < rfield.nupper; i++ )
00185 {
00186 opac.TauAbsGeo[1][i] = (realnum)(TAULIM*pow(rfield.anu[i],(realnum)-2.43f));
00187 opac.TauScatGeo[1][i] = opac.taumin;
00188 opac.TauTotalGeo[1][i] = opac.TauAbsGeo[1][i] + opac.TauScatGeo[1][i];
00189 }
00190
00191 }
00192 }
00193
00194 else
00195 {
00196
00197 for( i=0; i < (iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH1s] - 1); i++ )
00198 {
00199 opac.TauAbsGeo[1][i] = opac.taumin;
00200 opac.TauScatGeo[1][i] = opac.taumin;
00201 opac.TauTotalGeo[1][i] = opac.taumin;
00202 }
00203
00204 for( i=iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH1s]-1; i < rfield.nupper; i++ )
00205 {
00206 opac.TauAbsGeo[1][i] = (realnum)(TAULIM*pow(rfield.anu[i],(realnum)-2.43f));
00207 opac.TauScatGeo[1][i] = opac.taumin;
00208 opac.TauTotalGeo[1][i] = opac.TauAbsGeo[1][i] + opac.TauScatGeo[1][i];
00209 }
00210 }
00211
00212
00213
00214 if( geometry.lgSphere || opac.lgCaseB )
00215 {
00216 for( i=0; i < (iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH1s] - 1); i++ )
00217 {
00218 opac.TauAbsGeo[0][i] = opac.taumin;
00219 opac.TauAbsGeo[1][i] = opac.taumin*2.f;
00220 opac.TauScatGeo[0][i] = opac.taumin;
00221 opac.TauScatGeo[1][i] = opac.taumin*2.f;
00222 opac.TauTotalGeo[0][i] = 2.f*opac.taumin;
00223 opac.TauTotalGeo[1][i] = 4.f*opac.taumin;
00224 }
00225
00226 for( i=iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH1s]-1; i < rfield.nupper; i++ )
00227 {
00228 opac.TauAbsGeo[0][i] = opac.TauAbsGeo[1][i];
00229 opac.TauAbsGeo[1][i] *= 2.;
00230 opac.TauScatGeo[0][i] = opac.TauScatGeo[1][i];
00231 opac.TauScatGeo[1][i] *= 2.;
00232 opac.TauTotalGeo[0][i] = opac.TauTotalGeo[1][i];
00233 opac.TauTotalGeo[1][i] *= 2.;
00234 }
00235
00236 if( StopCalc.taunu > 0. )
00237 {
00238
00239 StopCalc.tauend *= 2.;
00240 }
00241 }
00242
00243
00244
00245 if( !thermal.lgTemperatureConstant )
00246 {
00247 double TeNew;
00248
00249
00250
00251 TeNew = 2e4;
00252
00253
00254
00255
00256 if( dense.gas_phase[ipHYDROGEN] >= 1e9 )
00257 {
00258
00259 TeNew = 7000.;
00260 }
00261 else if( dense.gas_phase[ipHYDROGEN] <= 1e5 )
00262 {
00263
00264 TeNew = 100.;
00265 }
00266 else
00267 {
00268
00269 TeNew = 0.5012 * pow( (double)dense.gas_phase[ipHYDROGEN], 0.46 );
00270 }
00271
00272
00273
00274 TempChange( TeNew );
00275 }
00276
00277
00278 dense.xNucleiTotal = 0.;
00279 for( nelem=0; nelem < LIMELM; nelem++ )
00280 {
00281 if( dense.lgElmtOn[nelem] )
00282 {
00283 dense.xNucleiTotal += dense.gas_phase[nelem];
00284 }
00285 }
00286
00287
00288 for( nelem=0; nelem < LIMELM; nelem++ )
00289 {
00290 if( dense.lgElmtOn[nelem] )
00291 {
00292
00293 AbunRatio = dense.gas_phase[nelem]/dense.xNucleiTotal;
00294 for( ipLo=ipH1s; ipLo < (iso.numLevels_max[ipH_LIKE][nelem] - 1); ipLo++ )
00295 {
00296 for( ipHi=ipLo + 1; ipHi < iso.numLevels_max[ipH_LIKE][nelem]; ipHi++ )
00297 {
00298
00299
00300 Transitions[ipH_LIKE][nelem][ipHi][ipLo].Emis->TauIn = opac.taumin;
00301 Transitions[ipH_LIKE][nelem][ipHi][ipLo].Emis->TauTot = 2.f*opac.taumin;
00302 }
00303 }
00304
00305
00306
00307 Transitions[ipH_LIKE][nelem][ipH2p][ipH1s].Emis->TauIn = opac.tlamin;
00308
00309
00310 f = opac.tlamin/Transitions[ipH_LIKE][nelem][ipH2p][ipH1s].Emis->opacity;
00311 fixit();
00312
00313 for( nHi=3; nHi<=iso.n_HighestResolved_max[ipH_LIKE][nelem]; nHi++ )
00314 {
00315 ipHi = iso.QuantumNumbers2Index[ipH_LIKE][nelem][nHi][1][2];
00316 Transitions[ipH_LIKE][nelem][ipHi][ipH1s].Emis->TauIn = MAX2( opac.taumin,
00317 f*Transitions[ipH_LIKE][nelem][ipHi][ipH1s].Emis->opacity );
00318 }
00319 for( ipHi=iso.numLevels_max[ipH_LIKE][nelem] - iso.nCollapsed_max[ipH_LIKE][nelem]; ipHi < iso.numLevels_max[ipH_LIKE][nelem]; ipHi++ )
00320 {
00321 Transitions[ipH_LIKE][nelem][ipHi][ipH1s].Emis->TauIn = MAX2( opac.taumin,
00322 f*Transitions[ipH_LIKE][nelem][ipHi][ipH1s].Emis->opacity );
00323 }
00324
00325
00326
00327
00328 if( opac.lgCaseB )
00329 {
00330
00331 Transitions[ipH_LIKE][nelem][ipH2p][ipH1s].Emis->TauTot =
00332 (realnum)(2.*Transitions[ipH_LIKE][nelem][ipH2p][ipH1s].Emis->TauIn);
00333
00334 BalmerTauOn = 0.f;
00335 }
00336
00337 else
00338 {
00339
00340 if( StopCalc.colnut < 6e29 )
00341 {
00342
00343 Transitions[ipH_LIKE][nelem][ipH2p][ipH1s].Emis->TauTot = (realnum)(StopCalc.colnut*
00344 rt.DoubleTau*Transitions[ipH_LIKE][nelem][ipH2p][ipH1s].Emis->opacity/
00345 GetDopplerWidth(dense.AtomicWeight[nelem])*AbunRatio);
00346 ASSERT( Transitions[ipH_LIKE][nelem][ipH2p][ipH1s].Emis->TauTot > 0. );
00347 }
00348
00349
00350
00351 else if( StopCalc.taunu < 3. && StopCalc.taunu >= 0.99 )
00352 {
00353
00354 Transitions[ipH_LIKE][nelem][ipH2p][ipH1s].Emis->TauTot = (realnum)(StopCalc.tauend*
00355 1.2e4*1.28e6/GetDopplerWidth(dense.AtomicWeight[nelem])*rt.DoubleTau*AbunRatio);
00356 ASSERT( Transitions[ipH_LIKE][nelem][ipH2p][ipH1s].Emis->TauTot > 0. );
00357 }
00358 else if( StopCalc.HColStop < 6e29 )
00359 {
00360
00361 coleff = StopCalc.HColStop -
00362 MIN2(MIN2(rfield.qhtot/dense.eden,1e24)/2.6e-13,0.6*StopCalc.HColStop);
00363
00364 Transitions[ipH_LIKE][nelem][ipH2p][ipH1s].Emis->TauTot = (realnum)(coleff*
00365 7.6e-14*AbunRatio);
00366 ASSERT( Transitions[ipH_LIKE][nelem][ipH2p][ipH1s].Emis->TauTot > 0. );
00367 }
00368 else
00369 {
00370
00371 Transitions[ipH_LIKE][nelem][ipH2p][ipH1s].Emis->TauTot = (realnum)(1e20*
00372 AbunRatio);
00373 ASSERT( Transitions[ipH_LIKE][nelem][ipH2p][ipH1s].Emis->TauTot > 0. );
00374 }
00375
00376 BalmerTauOn = 1.f;
00377 }
00378
00379
00380 if( Transitions[ipH_LIKE][nelem][ipH2p][ipH1s].Emis->TauTot > 1. )
00381 {
00382 rt.TAddHLya = (realnum)MIN2(1.,Transitions[ipH_LIKE][nelem][ipH2p][ipH1s].Emis->TauTot/
00383 1e4);
00384 Transitions[ipH_LIKE][nelem][ipH2p][ipH1s].Emis->TauIn += rt.TAddHLya;
00385 }
00386
00387 else
00388 {
00389 rt.TAddHLya = opac.tlamin;
00390 }
00391
00392
00393 f = Transitions[ipH_LIKE][nelem][ipH2p][ipH1s].Emis->TauTot/
00394 Transitions[ipH_LIKE][nelem][ipH2p][ipH1s].Emis->opacity;
00395
00396 ipISO = ipH_LIKE;
00397 ASSERT( ipISO<NISO && nelem < LIMELM );
00398 for( nHi=3; nHi<=iso.n_HighestResolved_max[ipISO][nelem]; nHi++ )
00399 {
00400 ipHi = iso.QuantumNumbers2Index[ipISO][nelem][nHi][1][2];
00401 Transitions[ipH_LIKE][nelem][ipHi][ipH1s].Emis->TauTot = MAX2( opac.taumin,
00402 Transitions[ipH_LIKE][nelem][ipHi][ipH1s].Emis->opacity * f );
00403
00404
00405
00406
00407 Transitions[ipH_LIKE][nelem][ipHi][ipH1s].Emis->TauIn += rt.TAddHLya*
00408 Transitions[ipH_LIKE][nelem][ipHi][ipH1s].Emis->opacity/
00409 Transitions[ipH_LIKE][nelem][ipH2p][ipH1s].Emis->opacity;
00410
00411 Transitions[ipH_LIKE][nelem][ipHi][ipH1s].Emis->TauIn = MAX2(
00412 opac.taumin, Transitions[ipH_LIKE][nelem][ipHi][ipH1s].Emis->TauIn );
00413 }
00414 for( ipHi=iso.numLevels_max[ipH_LIKE][nelem] - iso.nCollapsed_max[ipH_LIKE][nelem]; ipHi < iso.numLevels_max[ipH_LIKE][nelem]; ipHi++ )
00415 {
00416
00417 Transitions[ipH_LIKE][nelem][ipHi][ipH1s].Emis->TauTot = MAX2( opac.taumin,
00418 Transitions[ipH_LIKE][nelem][ipHi][ipH1s].Emis->opacity * f );
00419
00420
00421
00422
00423 Transitions[ipH_LIKE][nelem][ipHi][ipH1s].Emis->TauIn += rt.TAddHLya*
00424 Transitions[ipH_LIKE][nelem][ipHi][ipH1s].Emis->opacity/
00425 Transitions[ipH_LIKE][nelem][ipH2p][ipH1s].Emis->opacity;
00426
00427 Transitions[ipH_LIKE][nelem][ipHi][ipH1s].Emis->TauIn = MAX2(
00428 opac.taumin, Transitions[ipH_LIKE][nelem][ipHi][ipH1s].Emis->TauIn );
00429 }
00430
00431
00432
00433
00434 if( StopCalc.taunu > 0.24 && StopCalc.taunu < 0.7 )
00435 {
00436 Transitions[ipH_LIKE][nelem][ipH3p][ipH2s].Emis->TauTot = (realnum)(StopCalc.tauend*
00437 3.7e4*BalmerTauOn*AbunRatio + 1e-20);
00438
00439 Transitions[ipH_LIKE][nelem][ipH3s][ipH2p].Emis->TauTot = (realnum)(StopCalc.tauend*
00440 3.7e4*BalmerTauOn*AbunRatio + 1e-20);
00441
00442 Transitions[ipH_LIKE][nelem][ipH3d][ipH2p].Emis->TauTot = (realnum)(StopCalc.tauend*
00443 3.7e4*BalmerTauOn*AbunRatio + 1e-20);
00444 }
00445
00446 else
00447 {
00448
00449 balc = rfield.qhtot*2.1e-19*BalmerTauOn*AbunRatio + 1e-20;
00450
00451 Transitions[ipH_LIKE][nelem][ipH3p][ipH2s].Emis->TauTot =
00452 (realnum)(MIN2(2e5, balc*3.7e4*BalmerTauOn+1e-20));
00453
00454 Transitions[ipH_LIKE][nelem][ipH3s][ipH2p].Emis->TauTot =
00455 (realnum)(MIN2(2e5, balc*3.7e4*BalmerTauOn+1e-20));
00456
00457 Transitions[ipH_LIKE][nelem][ipH3d][ipH2p].Emis->TauTot =
00458 (realnum)(MIN2(2e5, balc*3.7e4*BalmerTauOn+1e-20));
00459
00460 ASSERT( Transitions[ipH_LIKE][nelem][ipH3p][ipH2s].Emis->TauTot > 0.);
00461 ASSERT( Transitions[ipH_LIKE][nelem][ipH3s][ipH2p].Emis->TauTot > 0.);
00462 ASSERT( Transitions[ipH_LIKE][nelem][ipH3d][ipH2p].Emis->TauTot > 0.);
00463 }
00464
00465
00466 Transitions[ipH_LIKE][nelem][ipH2s][ipH1s].Emis->TauTot = 2.f*opac.taumin;
00467 Transitions[ipH_LIKE][nelem][ipH2s][ipH1s].Emis->TauIn = opac.taumin;
00468
00469
00470
00471 f = Transitions[ipH_LIKE][nelem][ipH3s][ipH2p].Emis->TauTot/
00472 Transitions[ipH_LIKE][nelem][ipH3s][ipH2p].Emis->opacity* BalmerTauOn;
00473
00474 ipLo = ipH2s;
00475 for( ipHi=ipLo + 1; ipHi < iso.numLevels_max[ipH_LIKE][nelem]; ipHi++ )
00476 {
00477 if( Transitions[ipH_LIKE][nelem][ipHi][ipLo].ipCont <= 0 )
00478 continue;
00479
00480 Transitions[ipH_LIKE][nelem][ipHi][ipLo].Emis->TauTot = opac.taumin +
00481 f* Transitions[ipH_LIKE][nelem][ipHi][ipLo].Emis->opacity;
00482 ASSERT(Transitions[ipH_LIKE][nelem][ipHi][ipLo].Emis->TauTot > 0.);
00483 }
00484
00485
00486 for( ipLo=ipH2p; ipLo < (iso.numLevels_max[ipH_LIKE][nelem] - 1); ipLo++ )
00487 {
00488 long ipISO = ipH_LIKE, ipNS, ipNPlusOneP;
00489
00490
00491 ipNS = iso.QuantumNumbers2Index[ipH_LIKE][nelem][ N_(ipLo) ][0][2];
00492 if( ( N_(ipLo) + 1 ) <=
00493 (iso.n_HighestResolved_max[ipH_LIKE][nelem] + iso.nCollapsed_max[ipH_LIKE][nelem] ) )
00494 {
00495 ipNPlusOneP = iso.QuantumNumbers2Index[ipH_LIKE][nelem][ N_(ipLo)+1 ][1][2];
00496 }
00497 else
00498 {
00499 ipNPlusOneP = ipNS + 1;
00500 }
00501
00502 #if 1
00503 ipNS = ipLo;
00504 ipNPlusOneP = ipNS + 1;
00505 #endif
00506
00507 if( Transitions[ipH_LIKE][nelem][ipNPlusOneP][ipNS].ipCont <= 0 )
00508 {
00509 f = SMALLFLOAT;
00510 }
00511 else
00512 {
00513 f = Transitions[ipH_LIKE][nelem][ipNPlusOneP][ipNS].Emis->TauTot/
00514 Transitions[ipH_LIKE][nelem][ipNPlusOneP][ipNS].Emis->opacity*
00515 BalmerTauOn;
00516 }
00517
00518 for( ipHi=ipLo + 1; ipHi < iso.numLevels_max[ipH_LIKE][nelem]; ipHi++ )
00519 {
00520 if( Transitions[ipH_LIKE][nelem][ipHi][ipLo].ipCont <= 0 )
00521 continue;
00522
00523 Transitions[ipH_LIKE][nelem][ipHi][ipLo].Emis->TauTot = opac.taumin +
00524 f* Transitions[ipH_LIKE][nelem][ipHi][ipLo].Emis->opacity;
00525 ASSERT(Transitions[ipH_LIKE][nelem][ipHi][ipLo].Emis->TauTot > 0.);
00526 }
00527 }
00528
00529
00530 for( ipLo=ipH1s; ipLo < (iso.numLevels_max[ipH_LIKE][nelem] - 1); ipLo++ )
00531 {
00532 for( ipHi=ipLo + 1; ipHi < iso.numLevels_max[ipH_LIKE][nelem]; ipHi++ )
00533 {
00534 if( Transitions[ipH_LIKE][nelem][ipHi][ipLo].ipCont <= 0 )
00535 continue;
00536
00537
00538
00539
00540 Transitions[ipH_LIKE][nelem][ipHi][ipLo].Emis->TauCon =
00541 Transitions[ipH_LIKE][nelem][ipHi][ipLo].Emis->TauIn;
00542
00543
00544 Transitions[ipH_LIKE][nelem][ipHi][ipLo].Emis->TauIn =
00545 MIN2( Transitions[ipH_LIKE][nelem][ipHi][ipLo].Emis->TauIn ,
00546 Transitions[ipH_LIKE][nelem][ipHi][ipLo].Emis->TauTot/2.f );
00547 ASSERT(Transitions[ipH_LIKE][nelem][ipHi][ipLo].Emis->TauIn > 0.f);
00548
00549
00550 Transitions[ipH_LIKE][nelem][ipHi][ipLo].Emis->FracInwd = 0.5;
00551 }
00552 }
00553 }
00554 }
00555
00556
00557 for( nelem=ipHELIUM; nelem < LIMELM; nelem++ )
00558 {
00559 if( dense.lgElmtOn[nelem] )
00560 {
00561 for( ipLo=0; ipLo < (iso.numLevels_max[ipHE_LIKE][nelem] - 1); ipLo++ )
00562 {
00563 for( ipHi=ipLo + 1; ipHi < iso.numLevels_max[ipHE_LIKE][nelem]; ipHi++ )
00564 {
00565
00566
00567 Transitions[ipHE_LIKE][nelem][ipHi][ipLo].Emis->TauIn = opac.taumin;
00568 Transitions[ipHE_LIKE][nelem][ipHi][ipLo].Emis->TauTot = 2.f*opac.taumin;
00569 }
00570 }
00571 }
00572 }
00573
00574
00575 if( opac.lgCaseB )
00576 {
00577 for( nelem=1; nelem < LIMELM; nelem++ )
00578 {
00579 if( dense.lgElmtOn[nelem] )
00580 {
00581 double Aprev;
00582 realnum ratio;
00583
00584 Transitions[ipHE_LIKE][nelem][ipHe2p1P][ipHe1s1S].Emis->TauIn = opac.tlamin;
00585
00586 Transitions[ipHE_LIKE][nelem][ipHe2p1P][ipHe1s1S].Emis->TauCon =
00587 Transitions[ipHE_LIKE][nelem][ipHe2p1P][ipHe1s1S].Emis->TauIn;
00588
00589 Transitions[ipHE_LIKE][nelem][ipHe2p1P][ipHe1s1S].Emis->TauTot =
00590 2.f*Transitions[ipHE_LIKE][nelem][ipHe2p1P][ipHe1s1S].Emis->TauIn;
00591
00592 ratio = opac.tlamin/Transitions[ipHE_LIKE][nelem][ipHe2p1P][ipHe1s1S].Emis->opacity;
00593
00594
00595
00596 Aprev = Transitions[ipHE_LIKE][nelem][ipHe2p1P][ipHe1s1S].Emis->Aul;
00597
00598 i = ipHe2p1P+1;
00599
00600
00601 for( i=ipHe2p1P+1; i < iso.numLevels_max[ipHE_LIKE][nelem]; i++ )
00602 {
00603 if( Transitions[ipHE_LIKE][nelem][i][ipHe1s1S].ipCont <= 0 )
00604 continue;
00605
00606
00607
00608 if( Transitions[ipHE_LIKE][nelem][i][ipHe1s1S].Emis->Aul> Aprev/10. ||
00609 StatesElemNEW[nelem][nelem-ipHE_LIKE][i].S < 0 )
00610 {
00611 Aprev = Transitions[ipHE_LIKE][nelem][i][ipHe1s1S].Emis->Aul;
00612
00613 Transitions[ipHE_LIKE][nelem][i][ipHe1s1S].Emis->TauIn =
00614 ratio*Transitions[ipHE_LIKE][nelem][i][ipHe1s1S].Emis->opacity;
00615
00616 Transitions[ipHE_LIKE][nelem][i][ipHe1s1S].Emis->TauCon =
00617 Transitions[ipHE_LIKE][nelem][i][ipHe1s1S].Emis->TauIn;
00618 Transitions[ipHE_LIKE][nelem][i][ipHe1s1S].Emis->TauTot =
00619 2.f*Transitions[ipHE_LIKE][nelem][i][ipHe1s1S].Emis->TauIn;
00620 }
00621 }
00622 }
00623 }
00624 }
00625
00626
00627 lgHit = false;
00628 for( nelem=0; nelem < LIMELM; nelem++ )
00629 {
00630 if( dense.lgElmtOn[nelem] )
00631 {
00632 for( ipLo=ipH1s; ipLo < (iso.numLevels_max[ipH_LIKE][nelem] - 1); ipLo++ )
00633 {
00634 for( ipHi=ipLo + 1; ipHi < iso.numLevels_max[ipH_LIKE][nelem]; ipHi++ )
00635 {
00636 if( Transitions[ipH_LIKE][nelem][ipHi][ipLo].ipCont <= 0 )
00637 continue;
00638
00639 if( Transitions[ipH_LIKE][nelem][ipHi][ipLo].Emis->TauTot <
00640 0.9f*Transitions[ipH_LIKE][nelem][ipHi][ipLo].Emis->TauIn )
00641 {
00642 fprintf(ioQQQ,
00643 "RT_tau_init insanity for h line, Z=%li lo=%li hi=%li tot=%g in=%g \n",
00644 nelem , ipLo, ipHi , Transitions[ipH_LIKE][nelem][ipHi][ipLo].Emis->TauTot ,
00645 Transitions[ipH_LIKE][nelem][ipHi][ipLo].Emis->TauIn );
00646 lgHit = true;
00647 }
00648 }
00649 }
00650 }
00651 }
00652 if( lgHit )
00653 {
00654 fprintf( ioQQQ," stopping due to insanity in RT_tau_init\n");
00655 ShowMe();
00656 cdEXIT(EXIT_FAILURE);
00657 }
00658
00659
00660
00661 rt.tauxry = opac.TauAbsGeo[0][rt.ipxry-1];
00662
00663
00664
00665 conv.lgOscilOTS = false;
00666
00667
00668
00669
00670 conv.lgFirstSweepThisZone = true;
00671
00672
00673
00674
00675 if( 0 && trace.lgTrace )
00676 {
00677
00678 ipISO = ipH_LIKE;
00679 if( trace.lgIsoTraceFull[ipHE_LIKE] )
00680 ipISO = ipHE_LIKE;
00681
00682 if( trace.lgIsoTraceFull[ipISO] )
00683 {
00684 fprintf( ioQQQ, "\n\n up Transitions[ipISO][hi][lo].TauTot array as set in RT_tau_init ipZTrace (nelem)= %ld\n",
00685 trace.ipIsoTrace[ipH_LIKE] );
00686 for( ipHi=2; ipHi < iso.numLevels_max[ipISO][trace.ipIsoTrace[ipISO]]; ipHi++ )
00687 {
00688 fprintf( ioQQQ, " %3ld", ipHi );
00689 for( ipLo=ipH1s; ipLo < ipHi; ipLo++ )
00690 {
00691 if( Transitions[ipISO][trace.ipIsoTrace[ipISO]][ipHi][ipLo].Emis->Aul <= iso.SmallA )
00692 continue;
00693
00694 fprintf( ioQQQ,PrintEfmt("%9.2e",
00695 Transitions[ipISO][trace.ipIsoTrace[ipISO]][ipHi][ipLo].Emis->TauTot ));
00696 }
00697 fprintf( ioQQQ, "\n" );
00698 }
00699
00700 fprintf( ioQQQ, "\n\n TauIn array\n" );
00701 for( ipHi=1; ipHi < iso.numLevels_max[ipISO][trace.ipIsoTrace[ipISO]]; ipHi++ )
00702 {
00703 fprintf( ioQQQ, " %3ld", ipHi );
00704 for( ipLo=ipH1s; ipLo < ipHi; ipLo++ )
00705 {
00706 if( Transitions[ipISO][trace.ipIsoTrace[ipISO]][ipHi][ipLo].Emis->Aul <= iso.SmallA )
00707 continue;
00708
00709 fprintf( ioQQQ,PrintEfmt("%9.2e",
00710 Transitions[ipISO][trace.ipIsoTrace[ipISO]][ipHi][ipLo].Emis->TauIn ));
00711 }
00712 fprintf( ioQQQ, "\n" );
00713 }
00714
00715 fprintf( ioQQQ, "\n\n Aul As array\n" );
00716 for( ipHi=1; ipHi < iso.numLevels_max[ipISO][trace.ipIsoTrace[ipISO]]; ipHi++ )
00717 {
00718 fprintf( ioQQQ, " %3ld", ipHi );
00719 for( ipLo=ipH1s; ipLo < ipHi; ipLo++ )
00720 {
00721 if( Transitions[ipISO][trace.ipIsoTrace[ipISO]][ipHi][ipLo].Emis->Aul <= iso.SmallA )
00722 continue;
00723
00724 fprintf( ioQQQ,PrintEfmt("%9.2e",
00725 Transitions[ipISO][trace.ipIsoTrace[ipISO]][ipHi][ipLo].Emis->Aul) );
00726 }
00727 fprintf( ioQQQ, "\n" );
00728 }
00729
00730 fprintf( ioQQQ, "\n\n Aul*esc array\n" );
00731 for( ipHi=1; ipHi < iso.numLevels_max[ipISO][trace.ipIsoTrace[ipISO]]; ipHi++ )
00732 {
00733 fprintf( ioQQQ, " %3ld", ipHi );
00734 for( ipLo=ipH1s; ipLo < ipHi; ipLo++ )
00735 {
00736 if( Transitions[ipISO][trace.ipIsoTrace[ipISO]][ipHi][ipLo].Emis->Aul <= iso.SmallA )
00737 continue;
00738
00739 fprintf( ioQQQ,PrintEfmt("%9.2e",
00740 Transitions[ipISO][trace.ipIsoTrace[ipISO]][ipHi][ipLo].Emis->Aul*
00741 (Transitions[ipISO][trace.ipIsoTrace[ipISO]][ipHi][ipLo].Emis->Pdest +
00742 Transitions[ipISO][trace.ipIsoTrace[ipISO]][ipHi][ipLo].Emis->Pelec_esc +
00743 Transitions[ipISO][trace.ipIsoTrace[ipISO]][ipHi][ipLo].Emis->Pesc) ));
00744 }
00745 fprintf( ioQQQ, "\n" );
00746 }
00747
00748 fprintf( ioQQQ, "\n\n H opac array\n" );
00749 for( ipHi=1; ipHi < iso.numLevels_max[ipISO][trace.ipIsoTrace[ipISO]]; ipHi++ )
00750 {
00751 fprintf( ioQQQ, " %3ld", ipHi );
00752 for( ipLo=ipH1s; ipLo < ipHi; ipLo++ )
00753 {
00754 if( Transitions[ipISO][trace.ipIsoTrace[ipISO]][ipHi][ipLo].Emis->Aul <= iso.SmallA )
00755 continue;
00756
00757 fprintf( ioQQQ,PrintEfmt("%9.2e",
00758 Transitions[ipISO][trace.ipIsoTrace[ipISO]][ipHi][ipLo].Emis->opacity ));
00759 }
00760 fprintf( ioQQQ, "\n" );
00761 }
00762 }
00763
00764 else
00765 {
00766 fprintf( ioQQQ, " RT_tau_init called.\n" );
00767 }
00768 }
00769
00770 ASSERT( Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].Emis->TauIn> 0. );
00771 ASSERT( Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].Emis->PopOpc>= 0. );
00772 return;
00773 }