00001
00002
00003
00004
00005
00006
00007 #include "cddefines.h"
00008 #include "physconst.h"
00009 #include "taulines.h"
00010 #include "grains.h"
00011 #include "grainvar.h"
00012 #include "iso.h"
00013 #include "dense.h"
00014 #include "opacity.h"
00015 #include "trace.h"
00016 #include "coolheavy.h"
00017 #include "rfield.h"
00018 #include "phycon.h"
00019 #include "hmi.h"
00020 #include "radius.h"
00021 #include "atmdat.h"
00022 #include "heavy.h"
00023 #include "atomfeii.h"
00024 #include "lines_service.h"
00025 #include "h2.h"
00026 #include "ipoint.h"
00027 #include "rt.h"
00028 #include "mole.h"
00029 #include "conv.h"
00030
00031 #if defined (__ICC) && defined(__ia64) && __INTEL_COMPILER < 910
00032 #pragma optimization_level 0
00033 #endif
00034 void RT_diffuse(void)
00035 {
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062 long int ip=-100000,
00063 ipla=-100000,
00064 limit=-100000,
00065 nu=-10000;
00066
00067 double EdenAbund,
00068 difflya,
00069 fac,
00070 factor,
00071 gamma,
00072 gion,
00073 gn,
00074 photon;
00075
00076 DEBUG_ENTRY( "RT_diffuse()" );
00077
00078
00079
00080 ASSERT(rfield.nflux < rfield.nupper );
00081
00082
00083
00084 memset(rfield.DiffuseEscape , 0 , (unsigned)rfield.nupper*sizeof(realnum) );
00085 memset(rfield.ConEmitLocal[nzone] , 0 , (unsigned)rfield.nupper*sizeof(realnum) );
00086 memset(rfield.ConOTS_local_photons , 0 , (unsigned)rfield.nupper*sizeof(realnum) );
00087 memset(rfield.TotDiff2Pht , 0 , (unsigned)rfield.nupper*sizeof(realnum) );
00088 memset(rfield.DiffuseLineEmission , 0 , (unsigned)rfield.nupper*sizeof(realnum) );
00089
00090
00091
00092 if( lgAbort )
00093 {
00094
00095 return;
00096 }
00097
00098
00099
00100 for( long ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
00101 {
00102
00103 for( long nelem=ipISO; nelem < LIMELM; nelem++ )
00104 {
00105
00106 RT_iso_integrate_RRC( ipISO, nelem, true );
00107
00108
00109 EdenAbund = dense.eden*dense.xIonDense[nelem][nelem+1-ipISO];
00110
00111
00112
00113 if( dense.IonHigh[nelem] >= nelem+1-ipISO )
00114 {
00115 t_iso_sp* sp = &iso_sp[ipISO][nelem];
00116
00117
00118 for( long ipHi=1; ipHi < sp->numLevels_local; ipHi++ )
00119 {
00120 for( long ipLo=0; ipLo < ipHi; ipLo++ )
00121 {
00122
00123 if( sp->trans(ipHi,ipLo).ipCont() < 1 )
00124 continue;
00125
00126
00127
00128 sp->trans(ipHi,ipLo).Emis().phots() =
00129 sp->trans(ipHi,ipLo).Emis().Aul()*
00130 sp->st[ipHi].Pop()*
00131 sp->trans(ipHi,ipLo).Emis().Pesc();
00132
00133
00134
00135
00136
00137
00138
00139 const bool lgDoChecks = false;
00140 sp->trans(ipHi,ipLo).outline(1.0, lgDoChecks );
00141 }
00142 }
00143
00144
00145
00146
00147
00148 ASSERT( ipISO <= ipHE_LIKE );
00149
00150
00151 for( vector<two_photon>::iterator tnu = sp->TwoNu.begin(); tnu != sp->TwoNu.end(); ++tnu )
00152 {
00153 CalcTwoPhotonEmission( *tnu, rfield.lgInducProcess && iso_ctrl.lgInd2nu_On );
00154
00155 for( nu=0; nu < tnu->ipTwoPhoE; nu++ )
00156 {
00157
00158 rfield.TotDiff2Pht[nu] += tnu->local_emis[nu];
00159
00160
00161 rfield.ConEmitLocal[nzone][nu] += tnu->local_emis[nu];
00162
00163
00164
00165 rfield.DiffuseEscape[nu] += tnu->local_emis[nu] * opac.ExpmTau[nu];
00166
00167
00168
00169 rfield.ConOTS_local_photons[nu] += tnu->local_emis[nu] * (1.f - opac.ExpmTau[nu]);
00170 }
00171 enum {DEBUG_LOC=false};
00172 if( DEBUG_LOC )
00173 {
00174 fprintf( ioQQQ, "Two-photon emission coefficients - ipISO, nelem = %2li, %2li\n", ipISO, nelem );
00175 PrtTwoPhotonEmissCoef( *tnu, EdenAbund );
00176 }
00177 }
00178 }
00179 }
00180 }
00181
00182
00183 for( long nelem=NISO; nelem < LIMELM; nelem++ )
00184 {
00185
00186
00187 for( long ion=dense.IonLow[nelem]; ion < nelem-NISO+1; ion++ )
00188 {
00189 if( dense.xIonDense[nelem][ion+1] > 0. )
00190 {
00191 long int ns, nshell,igRec , igIon,
00192 iplow , iphi , ipop;
00193
00194 ip = Heavy.ipHeavy[nelem][ion]-1;
00195 ASSERT( ip >= 0 );
00196
00197
00198
00199
00200
00201
00202
00203
00204
00205 if( ip >= rfield.nflux )
00206 continue;
00207
00208
00209 atmdat_outer_shell( nelem+1 , nelem+1-ion , &nshell, &igRec , &igIon );
00210 gn = (double)igRec;
00211 gion = (double)igIon;
00212
00213
00214 ns = Heavy.nsShells[nelem][ion]-1;
00215 ASSERT( ns == (nshell-1) );
00216
00217
00218 iplow = opac.ipElement[nelem][ion][ns][0]-1;
00219 iphi = opac.ipElement[nelem][ion][ns][1];
00220 iphi = MIN2( iphi , rfield.nflux );
00221 ipop = opac.ipElement[nelem][ion][ns][2];
00222
00223
00224 ipop = ipop - iplow;
00225
00226 EdenAbund = dense.eden*dense.xIonDense[nelem][ion+1];
00227 gamma = 0.5*MILNE_CONST*gn/gion/phycon.te/phycon.sqrte;
00228
00229
00230 Heavy.RadRecCon[nelem][ion] = 0;
00231 if( rfield.ContBoltz[iplow] > SMALLFLOAT )
00232 {
00233 for( nu=iplow; nu < iphi; ++nu )
00234 {
00235 photon = gamma*rfield.ContBoltz[nu]/rfield.ContBoltz[iplow]*
00236 rfield.widflx[nu]*opac.OpacStack[nu+ipop]*rfield.anu2[nu];
00237
00242 rfield.ConEmitLocal[nzone][nu] += (realnum)photon*EdenAbund;
00243 rfield.DiffuseEscape[nu] += (realnum)photon*EdenAbund*opac.ExpmTau[nu];
00244
00245
00246 Heavy.RadRecCon[nelem][ion] += rfield.anu[nu] *
00247 emergent_line( photon*EdenAbund/2. , photon*EdenAbund/2. ,
00248
00249 nu+1 );
00250 }
00251 }
00252
00253 Heavy.RadRecCon[nelem][ion] *= EN1RYD;
00254
00255
00256 ipla = Heavy.ipLyHeavy[nelem][ion]-1;
00257 ASSERT( ipla >= 0 );
00258
00259 difflya = Heavy.xLyaHeavy[nelem][ion]*dense.xIonDense[nelem][ion+1];
00260 rfield.DiffuseLineEmission[ipla] += (realnum)difflya;
00261
00262
00263 rfield.outlin_noplot[ipla] += (realnum)(difflya*radius.dVolOutwrd*opac.tmn[ipla]*opac.ExpmTau[ipla]);
00264
00265
00266 ipla = Heavy.ipBalHeavy[nelem][ion]-1;
00267 ASSERT( ipla >= 0 );
00268
00269 difflya = Heavy.xLyaHeavy[nelem][ion]*dense.xIonDense[nelem][ion+1];
00270 rfield.outlin_noplot[ipla] += (realnum)(difflya*radius.dVolOutwrd*opac.tmn[ipla]*opac.ExpmTau[ipla]);
00271 }
00272 }
00273 }
00274
00275
00276 limit = MIN2( rfield.ipMaxBolt , rfield.nflux );
00277 for( nu=0; nu < limit; nu++ )
00278 {
00279 double TotBremsAllIons = 0., BremsThisIon;
00280
00281
00282
00283
00284 TotBremsAllIons += rfield.gff[1][nu] * opac.OpacStack[nu-1+opac.iphmra] * iso_sp[ipH_LIKE][ipHYDROGEN].st[ipH1s].Pop();
00285
00286
00287
00288 for( long nelem=ipHYDROGEN; nelem < LIMELM; nelem++ )
00289 {
00290
00291
00292 for( long ion=MAX2(1,dense.IonLow[nelem]); ion<=dense.IonHigh[nelem]; ++ion )
00293 {
00294
00295 BremsThisIon = POW2( (realnum)ion )*dense.xIonDense[nelem][ion]*rfield.gff[ion][nu];
00296 TotBremsAllIons += BremsThisIon;
00297 }
00298 }
00299
00300
00301 for( long ipMol = 0; ipMol<mole_global.num_calc; ipMol++ )
00302 {
00303 if( !mole_global.list[ipMol]->isMonatomic() && mole_global.list[ipMol]->charge > 0 && mole_global.list[ipMol]->parentLabel.empty()
00304
00305
00306 )
00307 {
00308
00309 long ion = mole_global.list[ipMol]->charge;
00310 BremsThisIon = POW2( (double)ion )*mole.species[ipMol].den*rfield.gff[ion][nu];
00311 TotBremsAllIons += BremsThisIon;
00312 }
00313 }
00314
00316
00317 TotBremsAllIons *= dense.eden*1.032e-11*rfield.widflx[nu]*rfield.ContBoltz[nu]/rfield.anu[nu]/phycon.sqrte *
00318 CoolHeavy.lgFreeOn;
00319 ASSERT( TotBremsAllIons >= 0.);
00320
00321
00322
00323
00324
00325
00326
00327 if( nu >= rfield.ipEnergyBremsThin )
00328 {
00329
00330
00331
00332
00333
00334
00335 rfield.ConEmitLocal[nzone][nu] += (realnum)TotBremsAllIons;
00336
00337
00338
00339
00340
00341
00342
00343 rfield.DiffuseEscape[nu] += (realnum)TotBremsAllIons;
00344 }
00345 }
00346
00347
00348
00349 if( gv.lgDustOn() && gv.lgGrainPhysicsOn )
00350 {
00351
00352
00353 GrainMakeDiffuse();
00354
00355 for( nu=0; nu < rfield.nflux; nu++ )
00356 {
00357 rfield.ConEmitLocal[nzone][nu] += gv.GrainEmission[nu];
00358 rfield.DiffuseEscape[nu] += gv.GrainEmission[nu];
00359 }
00360 }
00361
00362
00363 fac = dense.eden*(double)dense.xIonDense[ipHYDROGEN][0];
00364 gn = 1.;
00365 gion = 2.;
00366 gamma = 0.5*MILNE_CONST*gn/gion/phycon.te/phycon.sqrte;
00367
00368 limit = MIN2(iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH1s].ipIsoLevNIonCon-1,rfield.nflux);
00369
00370 if( rfield.ContBoltz[hmi.iphmin-1] > 0. )
00371 {
00372 for( nu=hmi.iphmin-1; nu < limit; nu++ )
00373 {
00374
00375
00376 factor = gamma*rfield.ContBoltz[nu]/rfield.ContBoltz[hmi.iphmin-1]*rfield.widflx[nu]*
00377 opac.OpacStack[nu-hmi.iphmin+opac.iphmop]*
00378 rfield.anu2[nu]*fac;
00379 rfield.ConEmitLocal[nzone][nu] += (realnum)factor;
00380 rfield.DiffuseEscape[nu] += (realnum)factor;
00381 }
00382 }
00383 else
00384 {
00385 for( nu=hmi.iphmin-1; nu < limit; nu++ )
00386 {
00387 double arg = MAX2(0.,TE1RYD*(rfield.anu[nu]-0.05544)/phycon.te);
00388
00389 if( arg > SEXP_LIMIT )
00390 break;
00391
00392
00393 factor = gamma*exp(-arg)*rfield.widflx[nu]*
00394 opac.OpacStack[nu-hmi.iphmin+opac.iphmop]*
00395 rfield.anu2[nu]*fac;
00396 rfield.ConEmitLocal[nzone][nu] += (realnum)factor;
00397 rfield.DiffuseEscape[nu] += (realnum)factor;
00398 }
00399 }
00400
00401
00402 for( long i=1; i <= nLevel1; i++ )
00403 TauLines[i].outline_resonance();
00404
00405 for( long ipISO = ipHE_LIKE; ipISO < NISO; ipISO++ )
00406 {
00407 for( long nelem = ipISO; nelem < LIMELM; nelem++ )
00408 {
00409 if( dense.lgElmtOn[nelem] && iso_ctrl.lgDielRecom[ipISO] )
00410 {
00411 for( long i=0; i<iso_sp[ipISO][nelem].numLevels_local; i++ )
00412 {
00413 const TransitionList::iterator& tr = SatelliteLines[ipISO][nelem].begin()+ipSatelliteLines[ipISO][nelem][i];
00414 (*tr).Emis().phots() =
00415 (*tr).Emis().Aul()*
00416 (*(*tr).Hi()).Pop()*
00417 ((*tr).Emis().Pesc()+
00418 (*tr).Emis().Pelec_esc());
00419
00420 (*tr).outline_resonance();
00421 }
00422 }
00423 }
00424 }
00425
00426
00427 for( long i=0; i < nWindLine; i++ )
00428 {
00429
00430
00431 if( (*TauLine2[i].Hi()).IonStg() < (*TauLine2[i].Hi()).nelem()+1-NISO )
00432 {
00433 {
00434 enum {DEBUG_LOC=false};
00435 if( DEBUG_LOC && i==4821 )
00436 {
00437
00438 fprintf(ioQQQ,"DEBUG dump lev2 line %li\n", i );
00439 DumpLine( TauLine2[i] );
00440 fprintf(ioQQQ,"DEBUG dump %.3e %.3e %.3e\n",
00441 rfield.outlin[0][TauLine2[i].ipCont()-1],
00442 TauLine2[i].Emis().phots()*TauLine2[i].Emis().FracInwd()*radius.BeamInOut*opac.tmn[i]*TauLine2[i].Emis().ColOvTot(),
00443 TauLine2[i].Emis().phots()*(1. - TauLine2[i].Emis().FracInwd())*radius.BeamOutOut* TauLine2[i].Emis().ColOvTot() );
00444 }
00445 }
00446 TauLine2[i].outline_resonance();
00447
00448
00449 }
00450 }
00451
00452
00453 for( long i=0; i < nHFLines; i++ )
00454 {
00455 HFLines[i].outline_resonance();
00456 }
00457
00458
00459 for( long ipSpecies=0; ipSpecies<nSpecies; ipSpecies++ )
00460 {
00461 if( dBaseSpecies[ipSpecies].lgActive )
00462 {
00463 for (TransitionList::iterator tr=dBaseTrans[ipSpecies].begin();
00464 tr != dBaseTrans[ipSpecies].end(); ++tr)
00465 {
00466 int ipHi = (*tr).ipHi();
00467 if (ipHi >= dBaseSpecies[ipSpecies].numLevels_local || (*tr).ipCont() <= 0)
00468 continue;
00469 (*tr).outline_resonance();
00470 }
00471 }
00472 }
00473
00474
00475 for( diatom_iter diatom = diatoms.begin(); diatom != diatoms.end(); ++diatom )
00476 (*diatom)->H2_RT_diffuse();
00477
00478
00479 FeII_RT_Out();
00483 if( trace.lgTrace )
00484 fprintf( ioQQQ, " RT_diffuse returns.\n" );
00485
00486
00487 for( nu=0; nu < rfield.ipPlasma-1; nu++ )
00488 {
00489 rfield.flux_beam_const[nu] = 0.;
00490 rfield.flux_beam_time[nu] = 0.;
00491 rfield.flux_isotropic[nu] = 0.;
00492 rfield.flux[0][nu] = 0.;
00493 rfield.ConEmitLocal[nzone][nu] = 0.;
00494 rfield.otscon[nu] = 0.;
00495 rfield.otslin[nu] = 0.;
00496 rfield.outlin[0][nu] = 0.;
00497 rfield.outlin_noplot[nu] = 0.;
00498 rfield.reflin[0][nu] = 0.;
00499 rfield.TotDiff2Pht[nu] = 0.;
00500 rfield.ConOTS_local_photons[nu] = 0.;
00501 rfield.ConInterOut[nu] = 0.;
00502 }
00503
00504
00505 for( nu=0; nu < rfield.nflux; nu++ )
00506 {
00507
00508
00509 rfield.OccNumbDiffCont[nu] =
00510 rfield.ConEmitLocal[nzone][nu]*rfield.convoc[nu];
00511
00512
00513 rfield.ConSourceFcnLocal[nzone][nu] =
00514
00515 (realnum)safe_div( (double)rfield.ConEmitLocal[nzone][nu],
00516
00517 opac.opacity_abs[nu] );
00518
00519
00520 ASSERT( rfield.flux_beam_const[nu] >= 0.);
00521 ASSERT( rfield.flux_beam_time[nu] >= 0.);
00522 ASSERT( rfield.flux_isotropic[nu] >= 0.);
00523 ASSERT( rfield.flux[0][nu] >= 0.);
00524 ASSERT( rfield.ConEmitLocal[nzone][nu] >= 0.);
00525 ASSERT( rfield.otscon[nu] >= 0.);
00526 ASSERT( rfield.otslin[nu] >= 0.);
00527 ASSERT( rfield.outlin[0][nu] >= 0.);
00528 ASSERT( rfield.outlin_noplot[nu] >= 0.);
00529 ASSERT( rfield.reflin[0][nu] >= 0.);
00530 ASSERT( rfield.TotDiff2Pht[nu] >= 0.);
00531 ASSERT( rfield.ConOTS_local_photons[nu] >= 0.);
00532 ASSERT( rfield.ConInterOut[nu] >= 0.);
00533 }
00534
00535
00536 if( rfield.lgKillOutLine )
00537 {
00538 for( nu=0; nu < rfield.nflux; nu++ )
00539 {
00540 rfield.outlin[0][nu] = 0.;
00541 rfield.outlin_noplot[nu] = 0.;
00542 }
00543 }
00544
00545
00546 if( rfield.lgKillOutCont )
00547 {
00548 for( nu=0; nu < rfield.nflux; nu++ )
00549 {
00550 rfield.ConInterOut[nu] = 0.;
00551 }
00552 }
00553 return;
00554 }
00555
00556 void RT_iso_integrate_RRC( const long ipISO, const long nelem, const bool lgUpdateContinuum )
00557 {
00558 DEBUG_ENTRY( "RT_iso_integrate_RRC()" );
00559
00560
00561 static double TeUsed[NISO][LIMELM]={
00562 {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0},
00563 {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0} };
00564
00565 if( !lgUpdateContinuum && fp_equal( phycon.te, TeUsed[ipISO][nelem] ) && conv.nTotalIoniz )
00566 return;
00567
00568 ASSERT( nelem >= ipISO );
00569 ASSERT( nelem < LIMELM );
00570
00571
00572 double SumCaseB = 0.;
00573
00574
00575 double EdenAbund = dense.eden*dense.xIonDense[nelem][nelem+1-ipISO];
00576
00577
00578
00579 if( dense.IonHigh[nelem] >= nelem+1-ipISO )
00580 {
00581 t_iso_sp* sp = &iso_sp[ipISO][nelem];
00582
00583
00584
00585 long ipHi = rfield.nflux;
00586
00587 for( long n=0; n < sp->numLevels_local; n++ )
00588 {
00589 double Sum1level = 0.;
00590 sp->fb[n].RadRecCon = 0.;
00591 sp->fb[n].RadRecCoolCoef = 0.;
00592
00593
00594 double gamma = 0.5*MILNE_CONST*sp->st[n].g()/iso_ctrl.stat_ion[ipISO]/phycon.te/phycon.sqrte;
00595
00596
00597
00598
00599 for( long nu=sp->fb[n].ipIsoLevNIonCon-1; nu < ipHi; nu++ )
00600 {
00601
00602
00603 double dwid = 0.2;
00604
00605
00606
00607 double arg = (rfield.anu[nu]-sp->fb[n].xIsoLevNIonRyd+
00608 rfield.widflx[nu]*dwid)/phycon.te_ryd;
00609 arg = MAX2(0.,arg);
00610
00611
00612 if( arg > SEXP_LIMIT )
00613 break;
00614
00615
00616 double photon = gamma*exp(-arg)*rfield.widflx[nu]*
00617 opac.OpacStack[ nu-sp->fb[n].ipIsoLevNIonCon + sp->fb[n].ipOpac ] *
00618 rfield.anu2[nu];
00619
00620 Sum1level += photon;
00621
00622 fixit();
00623
00624
00625 if( lgUpdateContinuum )
00626 {
00627
00628 rfield.ConEmitLocal[nzone][nu] += (realnum)(photon*EdenAbund);
00629
00630
00631
00632 rfield.DiffuseEscape[nu] +=
00633 (realnum)(photon*EdenAbund*sp->fb[n].RadRecomb[ipRecEsc]);
00634 }
00635
00636
00637 sp->fb[n].RadRecCon += rfield.anu[nu] *
00638 emergent_line( photon*EdenAbund/2. , photon*EdenAbund/2. ,
00639
00640 nu+1 );
00641
00642 double energyAboveThresh = rfield.anu[nu] - sp->fb[n].xIsoLevNIonRyd;
00643 energyAboveThresh = MAX2( 0., energyAboveThresh );
00644 sp->fb[n].RadRecCoolCoef += energyAboveThresh * photon *
00645 sp->fb[n].RadRecomb[ipRecNetEsc];
00646 }
00647
00648
00649 sp->fb[n].RadRecCon *= EN1RYD;
00650 sp->fb[n].RadRecCoolCoef *= EN1RYD;
00651
00652 if( n > 0 )
00653 {
00654
00655 SumCaseB += Sum1level;
00656 }
00657 }
00658
00659
00660 for( long n=sp->numLevels_local;n<sp->numLevels_max; n++ )
00661 {
00662 sp->fb[n].RadRecCon = 0.;
00663 sp->fb[n].RadRecCoolCoef = 0.;
00664 }
00665
00666
00667 sp->CaseBCheck = MAX2(sp->CaseBCheck,
00668 (realnum)(SumCaseB/sp->RadRec_caseB));
00669 }
00670
00671 TeUsed[ipISO][nelem] = phycon.te;
00672
00673 return;
00674 }
00675