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