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