00001
00002
00003
00004
00005
00006
00007
00008
00009
00010 #include "cddefines.h"
00011 #include "rfield.h"
00012 #include "doppvel.h"
00013 #include "dense.h"
00014 #include "opacity.h"
00015 #include "transition.h"
00016 #include "conv.h"
00017 #include "radius.h"
00018 #include "rt.h"
00019 #include "physconst.h"
00020 #include "cosmology.h"
00021 #include "thirdparty.h"
00022 #include "hydrogenic.h"
00023
00024
00025 STATIC void RT_line_pumping(
00026
00027 const TransitionProxy& t ,
00028
00029
00030
00031
00032 bool lgShield_this_zone,
00033 realnum DopplerWidth)
00034 {
00035 DEBUG_ENTRY( "RT_line_pumping()" );
00036
00037 ASSERT( t.ipCont() >= 1 );
00038
00039
00040
00041 if( !rfield.lgInducProcess )
00042 {
00043 t.Emis().pump() = 0.;
00044 }
00045 else if( conv.lgFirstSweepThisZone || t.Emis().iRedisFun() == ipLY_A )
00046 {
00047 double dTau;
00048 double shield_continuum;
00049
00050
00051 shield_continuum = RT_continuum_shield_fcn( t );
00052
00053
00054
00055
00056 t.Emis().pump() = t.Emis().Aul() * (*t.Hi()).g() / (*t.Lo()).g() * shield_continuum *(
00057 rfield.OccNumbIncidCont[t.ipCont()-1] + rfield.OccNumbContEmitOut[t.ipCont()-1] );
00058
00059 dTau = (t.Emis().PopOpc() * t.Emis().opacity() / DopplerWidth +
00060 opac.opacity_abs[t.ipCont()-1])* radius.drad_x_fillfac;
00061
00062 if( lgShield_this_zone && dTau > 1e-3 )
00063 {
00064
00065 t.Emis().pump() *= log(1. + dTau ) / dTau;
00066 }
00067
00068
00069
00070
00071
00072
00073
00074
00075
00076 if ( t.ipLo() == 0 && t.systemIs(iso_sp[ipH_LIKE][ipHYDROGEN].tr) )
00077 {
00078 t.Emis().pump() *= hydro.xLymanPumpingScaleFactor;
00079 }
00080
00081 }
00082
00083 return;
00084 }
00085
00086
00087 STATIC void RT_line_electron_scatter(
00088
00089 const TransitionProxy& t ,
00090 realnum DopplerWidth)
00091 {
00092
00093 DEBUG_ENTRY( "RT_line_electron_scatter()" );
00094
00095
00096
00097 if( !rt.lgElecScatEscape )
00098 {
00099 t.Emis().Pelec_esc() = 0.;
00100 return;
00101 }
00102
00103
00104
00105
00106 double opac_line = t.Emis().PopOpc() * t.Emis().opacity()/DopplerWidth;
00107
00108 if( opac_line > SMALLFLOAT )
00109 {
00110
00111 double opac_electron = dense.eden*6.65e-25;
00112
00113
00114 double opacity_ratio = opac_electron/(opac_electron+opac_line) * (1.-t.Emis().Pesc());
00115
00116 t.Emis().Pelec_esc() = (realnum)opacity_ratio * MAX2(0.f,1.f-t.Emis().Pesc()-t.Emis().Pdest());
00117 }
00118 else
00119
00120 t.Emis().Pelec_esc() *= 0.95f;
00121
00122 return;
00123 }
00124
00125
00126
00127 STATIC void RT_line_escape(
00128
00129 const TransitionProxy &t,
00130
00131 realnum pestrk,
00132 realnum DopplerWidth,
00133 bool lgGoodTau)
00134 {
00135 int nRedis = -1;
00136
00137 DEBUG_ENTRY( "RT_line_escape()" );
00138
00139
00140 if( t.Emis().TauIn() < -30. )
00141 {
00142 fprintf( ioQQQ, "PROBLEM RT_line_escape called with large negative "
00143 "optical depth, zone %.2f, setting lgAbort true.\n",
00144 fnzone );
00145 DumpLine(t);
00146
00147
00148 lgAbort = true;
00149 return;
00150 }
00151
00152 if( cosmology.lgDo )
00153 {
00154
00155 if( conv.lgFirstSweepThisZone && lgGoodTau )
00156 {
00157 realnum tau_Sobolev = t.Emis().TauIn();
00158
00159 if( tau_Sobolev < 1E-5 )
00160 {
00161 t.Emis().Pesc() = 1.;
00162 }
00163 else
00164 {
00165 t.Emis().Pesc() = ( 1.f - exp( -1.f * tau_Sobolev ) )/ tau_Sobolev;
00166 }
00167
00168
00169 t.Emis().FracInwd() = rt.fracin;
00170 }
00171 fixit();
00172 nRedis = ipDEST_K2;
00173 }
00174
00175
00176
00177
00178 else if( t.Emis().iRedisFun() == ipPRD )
00179 {
00180
00181 if( conv.lgFirstSweepThisZone && lgGoodTau )
00182 {
00183 t.Emis().Pesc() = (realnum)esc_PRD( t.Emis().TauIn(), t.Emis().TauTot(), t.Emis().damp() );
00184
00185
00186
00187 if( pestrk > 0.f && t.Emis().Pesc() < 1.f )
00188 t.Emis().Pesc() = min( 1.f, t.Emis().Pesc() + pestrk );
00189
00190
00191 t.Emis().FracInwd() = rt.fracin;
00192 }
00193 nRedis = ipDEST_INCOM;
00194 }
00195
00196
00197 else if( t.Emis().iRedisFun() == ipCRD )
00198 {
00199 if( conv.lgFirstSweepThisZone && lgGoodTau )
00200 {
00201
00202
00203
00204 t.Emis().Pesc() = (realnum)esc_CRDcore( t.Emis().TauIn(), t.Emis().TauTot() );
00205
00206 if( pestrk > 0.f && t.Emis().Pesc() < 1.f )
00207 t.Emis().Pesc() = min( 1.f, t.Emis().Pesc() + pestrk );
00208
00209
00210 t.Emis().FracInwd() = rt.fracin;
00211 }
00212 nRedis = ipDEST_K2;
00213 }
00214
00215
00216 else if( t.Emis().iRedisFun() == ipCRDW )
00217 {
00218
00219 if( conv.lgFirstSweepThisZone && lgGoodTau )
00220 {
00221 t.Emis().Pesc() = (realnum)esc_CRDwing( t.Emis().TauIn(), t.Emis().TauTot(), t.Emis().damp() );
00222
00223 if( pestrk > 0.f && t.Emis().Pesc() < 1.f )
00224 t.Emis().Pesc() = min( 1.f, t.Emis().Pesc() + pestrk );
00225
00226
00227 t.Emis().FracInwd() = rt.fracin;
00228 }
00229 nRedis = ipDEST_K2;
00230 }
00231
00232
00233 else if( t.Emis().iRedisFun() == ipLY_A )
00234 {
00235
00236
00237
00238
00239 if( lgGoodTau )
00240 {
00241 double dest , esin;
00242
00243
00244
00245 t.Emis().Pesc() = (realnum)RTesc_lya( &esin, &dest, t.Emis().PopOpc(), t, DopplerWidth );
00246
00247 if( pestrk > 0.f && t.Emis().Pesc() < 1.f )
00248 t.Emis().Pesc() = min( 1.f, t.Emis().Pesc() + pestrk );
00249
00250
00251 t.Emis().Pdest() = (realnum)dest;
00252
00253
00254 t.Emis().FracInwd() = rt.fracin;
00255 }
00256 }
00257 else
00258 {
00259 fprintf( ioQQQ, " RT_line_escape called with impossible redistribution function %d\n",
00260 t.Emis().iRedisFun());
00261 ShowMe();
00262 cdEXIT(EXIT_FAILURE);
00263 }
00264
00265
00266 if( lgGoodTau && t.Emis().iRedisFun() != ipLY_A && t.Emis().opacity() > 0. )
00267 {
00268 t.Emis().Pdest() = (realnum)RT_DestProb(
00269
00270 t.Emis().PopOpc(),
00271
00272 t.Emis().opacity(),
00273
00274
00275 t.ipCont(),
00276
00277 DopplerWidth,
00278
00279 t.Emis().Pesc(),
00280
00281 nRedis);
00282 }
00283
00284 return;
00285 }
00286
00287
00288 STATIC void RT_line_fine_opacity(
00289
00290 const TransitionProxy& t ,
00291 realnum DopplerWidth)
00292
00293 {
00294 DEBUG_ENTRY( "RT_line_fine_opacity()" );
00295
00296
00297 long int ipLineCenter = t.Emis().ipFine() + rfield.ipFineConVelShift;
00298
00299
00300
00301 if( !conv.lgLastSweepThisZone || ipLineCenter < 0 || t.Emis().PopOpc() < SMALLFLOAT ||
00302 ipLineCenter>rfield.nfine || !rfield.lgOpacityFine )
00303 {
00304 return;
00305 }
00306
00307
00308
00309
00310 realnum cells_wide_1x = DopplerWidth/rfield.fine_opac_velocity_width;
00311
00312
00313
00314 realnum opac_line = (realnum)t.Emis().PopOpc() * t.Emis().opacity() / DopplerWidth;
00315
00316
00317
00318 double dTauEffec = opac_line*radius.depth_x_fillfac;
00319 if( dTauEffec < SMALLFLOAT )
00320 return;
00321
00322
00323
00324 const bool doDamp = dTauEffec*t.Emis().damp()/9. > 0.1;
00325 long int nCells_core = (long)(cells_wide_1x*4.f + 1.5f);
00326
00327
00328 if( ipLineCenter - nCells_core < 1 )
00329 nCells_core = ipLineCenter - 1;
00330 if( ipLineCenter + nCells_core > rfield.nfine )
00331 nCells_core = ipLineCenter -rfield.nfine - 1;
00332
00333
00334 nCells_core = MAX2( 1 , nCells_core );
00335
00336 long int nCells_damp;
00337
00338 if( doDamp )
00339 {
00340
00341
00342
00343 realnum x = (realnum)sqrt( dTauEffec * t.Emis().damp()*100./SQRTPI ) * cells_wide_1x;
00344
00345
00346 if( x<LONG_MAX )
00347 {
00348 nCells_damp = (long)x;
00349
00350 if( ipLineCenter-nCells_damp < 1 )
00351 nCells_damp = ipLineCenter-1;
00352
00353 if( ipLineCenter+nCells_damp > rfield.nfine )
00354 nCells_damp = rfield.nfine - ipLineCenter-1;
00355 }
00356 else
00357 {
00358
00359
00360 nCells_damp = MIN2( rfield.nfine-ipLineCenter , ipLineCenter )-1;
00361 }
00362 }
00363 else
00364 {
00365 nCells_damp = nCells_core;
00366 }
00367
00368 static vector<realnum> xprofile, profile;
00369 xprofile.resize(nCells_damp);
00370 profile.resize(nCells_damp);
00371
00372 for( long int i=0; i<nCells_damp; ++i )
00373 {
00374
00375 xprofile[i] = (realnum) i/cells_wide_1x;
00376 }
00377
00378 VoigtH(t.Emis().damp(), &xprofile[0], &profile[0], nCells_damp);
00379
00380
00381 rfield.fine_opac_zone[ipLineCenter] += profile[0]*opac_line;
00382 for( long int i=1; i<nCells_damp; ++i )
00383 {
00384 rfield.fine_opac_zone[ipLineCenter+i] += profile[i]*opac_line;
00385 rfield.fine_opac_zone[ipLineCenter-i] += profile[i]*opac_line;
00386 }
00387
00388 return;
00389 }
00390
00391
00392 void RT_line_one(
00393
00394 const TransitionProxy &t,
00395
00396
00397
00398
00399 bool lgShield_this_zone,
00400
00401 realnum pestrk,
00402 realnum DopplerWidth )
00403 {
00404 DEBUG_ENTRY( "RT_line_one()" );
00405
00406
00407
00408 if( !rfield.lgDoLineTrans && (t.Emis().iRedisFun() != ipLY_A) )
00409 {
00410 return;
00411 }
00412
00413
00414 t.Emis().damp() = t.Emis().dampXvel() / DopplerWidth;
00415 ASSERT( t.Emis().damp() > 0. );
00416
00417
00418 if( (*t.Lo()).Pop()<=SMALLFLOAT )
00419 {
00420
00421 t.Emis().Pesc() = 1.f;
00422
00423
00424 t.Emis().FracInwd() = 0.5;
00425
00426
00427 t.Emis().pump() = 0.;
00428
00429
00430 t.Emis().Pdest() = 0.;
00431 t.Emis().Pelec_esc() = 0.;
00432
00433 return;
00434 }
00435
00436
00437
00438 enum {DEBUG_LOC=false};
00439 if( DEBUG_LOC )
00440 {
00441 static long int nTau[100];
00442 long n;
00443
00444 if( nzone==0 )
00445 {
00446 for(n=0; n<100; ++n )
00447 nTau[n] = 0;
00448 }
00449 if( (*t.Lo()).Pop()<=SMALLFLOAT )
00450 n = 0;
00451 else
00452 n = (long)log10( (*t.Lo()).Pop() )+37;
00453 n = MIN2( n , 99 );
00454 n = MAX2( n , 0 );
00455 ++nTau[n];
00456 if( nzone > 183 )
00457 {
00458 for(n=0; n<100; ++n )
00459 fprintf(ioQQQ,"%li\t%li\n", n , nTau[n] );
00460 cdEXIT(EXIT_SUCCESS);
00461 }
00462 }
00463
00464
00465 if( t.EnergyErg() / EN1RYD <= rfield.plsfrq )
00466 {
00467 t.Emis().Pesc() = SMALLFLOAT;
00468 t.Emis().Pdest() = SMALLFLOAT;
00469 t.Emis().Pelec_esc() = SMALLFLOAT;
00470 t.Emis().pump() = SMALLFLOAT;
00471 }
00472 else
00473 {
00474
00475
00476
00477
00478
00479
00480
00481 bool lgGoodTau = lgTauGood( t );
00482
00483
00484
00485
00486 if( conv.lgLastSweepThisZone )
00487 RT_line_fine_opacity( t , DopplerWidth );
00488 else
00489 {
00490 RT_line_escape( t, pestrk, DopplerWidth , lgGoodTau);
00491 RT_line_electron_scatter( t , DopplerWidth );
00492 RT_line_pumping( t , lgShield_this_zone , DopplerWidth );
00493 }
00494 }
00495
00496 return;
00497 }
00498