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