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