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
00023 STATIC void getVoigt(long int nCells_core, long int nCells_damp,
00024 realnum *xprofile, realnum *profile, realnum damp);
00025
00026
00027 STATIC void RT_line_pumping(
00028
00029 transition* t ,
00030
00031
00032
00033
00034 bool lgShield_this_zone,
00035 realnum DopplerWidth)
00036 {
00037 DEBUG_ENTRY( "RT_line_pumping()" );
00038
00039 ASSERT( t->ipCont >= 1 );
00040
00041
00042
00043 if( !rfield.lgInducProcess )
00044 {
00045 t->Emis->pump = 0.;
00046 }
00047 else if( conv.lgFirstSweepThisZone || t->Emis->iRedisFun == ipLY_A )
00048 {
00049 double dTau;
00050 double shield_continuum;
00051
00052
00053 shield_continuum = RT_continuum_shield_fcn( t );
00054
00055
00056
00057
00058 t->Emis->pump = t->Emis->Aul * t->Hi->g / t->Lo->g * shield_continuum *(
00059 rfield.OccNumbIncidCont[t->ipCont-1] + rfield.OccNumbContEmitOut[t->ipCont-1] );
00060
00061
00062
00063
00064 dTau = (t->Emis->PopOpc * t->Emis->opacity / DopplerWidth +
00065 opac.opacity_abs[t->ipCont-1])* radius.drad_x_fillfac_mean;
00066
00067
00068
00069 dTau = MIN2(1., dTau);
00070
00071 if( lgShield_this_zone && dTau > 1e-3 && dTau < 1e10 )
00072 {
00073
00074 t->Emis->pump *= log(1. + dTau ) / dTau;
00075 }
00076
00077
00078 }
00079
00080 return;
00081 }
00082
00083
00084 STATIC void RT_line_electron_scatter(
00085
00086 transition* t ,
00087 realnum DopplerWidth)
00088 {
00089
00090 DEBUG_ENTRY( "RT_line_electron_scatter()" );
00091
00092
00093
00094 if( !rt.lgElecScatEscape )
00095 {
00096 t->Emis->Pelec_esc = 0.;
00097 return;
00098 }
00099
00100
00101
00102
00103 double opac_line = t->Emis->PopOpc * t->Emis->opacity/DopplerWidth;
00104
00105 if( opac_line > SMALLFLOAT )
00106 {
00107
00108 double opac_electron = dense.eden*6.65e-25;
00109
00110
00111 double opacity_ratio = opac_electron/(opac_electron+opac_line) * (1.-t->Emis->Pesc);
00112
00113 t->Emis->Pelec_esc = (realnum)opacity_ratio * MAX2(0.f,1.f-t->Emis->Pesc-t->Emis->Pdest);
00114 }
00115 else
00116
00117 t->Emis->Pelec_esc *= 0.95f;
00118
00119 return;
00120 }
00121
00122
00123
00124 STATIC void RT_line_escape(
00125
00126 transition* t,
00127
00128 realnum pestrk,
00129 realnum DopplerWidth,
00130 bool lgGoodTau)
00131 {
00132 int nRedis = -1;
00133
00134 DEBUG_ENTRY( "RT_line_escape()" );
00135
00136
00137 if( t->Emis->TauIn < -30. )
00138 {
00139 fprintf( ioQQQ, "PROBLEM RT_line_escape called with large negative "
00140 "optical depth, zone %.2f, setting lgAbort true.\n",
00141 fnzone );
00142 DumpLine(t);
00143
00144
00145 lgAbort = true;
00146 return;
00147 }
00148
00149
00150 long int nelem = t->Hi->nelem-1;
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, nelem, 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 transition* 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 getVoigt(nCells_core, nCells_damp, &xprofile[0], &profile[0],
00379 t->Emis->damp);
00380
00381
00382 rfield.fine_opac_zone[ipLineCenter] += profile[0]*opac_line;
00383 for( long int i=1; i<nCells_damp; ++i )
00384 {
00385 rfield.fine_opac_zone[ipLineCenter+i] += profile[i]*opac_line;
00386 rfield.fine_opac_zone[ipLineCenter-i] += profile[i]*opac_line;
00387 }
00388
00389 return;
00390 }
00391
00392 STATIC void getVoigt(long int nCells_core, long int nCells_damp,
00393 realnum *xprofile, realnum *profile, realnum damp)
00394 {
00395 if (1)
00396 {
00397 profile[0] = 1.;
00398 for( long int i=1; i<nCells_core; ++i )
00399 {
00400 realnum xsq = POW2(xprofile[i]);
00401
00402
00403
00404 profile[i] = sexp(xsq) + damp / MAX2(1.f,((realnum)SQRTPI*xsq));
00405
00406 }
00407 ASSERT( nCells_core>0 );
00408 for( long int i=nCells_core; i<nCells_damp; ++i )
00409 {
00410
00411 realnum xsq = POW2(xprofile[i]);
00412 profile[i] = damp/((realnum)SQRTPI*xsq);
00413 }
00414 }
00415 else
00416 {
00417 humlik(nCells_damp,xprofile,damp,profile);
00418 }
00419
00420 }
00421
00422
00423 void RT_line_one(
00424
00425 transition* t,
00426
00427
00428
00429
00430 bool lgShield_this_zone,
00431
00432 realnum pestrk,
00433 realnum DopplerWidth )
00434 {
00435 DEBUG_ENTRY( "RT_line_one()" );
00436
00437
00438
00439 if( !rfield.lgDoLineTrans && (t->Emis->iRedisFun != ipLY_A) )
00440 {
00441 return;
00442 }
00443
00444
00445 t->Emis->damp = t->Emis->dampXvel / DopplerWidth;
00446 ASSERT( t->Emis->damp > 0. );
00447
00448
00449 if( t->Lo->Pop<=SMALLFLOAT )
00450 {
00451
00452 t->Emis->Pesc = 1.f;
00453
00454
00455 t->Emis->FracInwd = 0.5;
00456
00457
00458 t->Emis->pump = 0.;
00459
00460
00461 t->Emis->Pdest = 0.;
00462 t->Emis->Pelec_esc = 0.;
00463
00464 return;
00465 }
00466
00467
00468
00469 enum {DEBUG_LOC=false};
00470 if( DEBUG_LOC )
00471 {
00472 static long int nTau[100];
00473 long n;
00474
00475 if( nzone==0 )
00476 {
00477 for(n=0; n<100; ++n )
00478 nTau[n] = 0;
00479 }
00480 if( t->Lo->Pop<=SMALLFLOAT )
00481 n = 0;
00482 else
00483 n = (long)log10( t->Lo->Pop )+37;
00484 n = MIN2( n , 99 );
00485 n = MAX2( n , 0 );
00486 ++nTau[n];
00487 if( nzone > 183 )
00488 {
00489 for(n=0; n<100; ++n )
00490 fprintf(ioQQQ,"%li\t%li\n", n , nTau[n] );
00491 cdEXIT(EXIT_SUCCESS);
00492 }
00493 }
00494
00495
00496 if( t->EnergyErg / EN1RYD <= rfield.plsfrq )
00497 {
00498 t->Emis->Pesc = SMALLFLOAT;
00499 t->Emis->Pdest = SMALLFLOAT;
00500 t->Emis->Pelec_esc = SMALLFLOAT;
00501 t->Emis->pump = SMALLFLOAT;
00502 }
00503 else
00504 {
00505
00506
00507
00508
00509
00510
00511
00512 bool lgGoodTau = lgTauGood( t );
00513
00514
00515
00516
00517 if( conv.lgLastSweepThisZone )
00518 RT_line_fine_opacity( t , DopplerWidth );
00519 else
00520 {
00521 RT_line_escape( t, pestrk, DopplerWidth , lgGoodTau);
00522 RT_line_electron_scatter( t , DopplerWidth );
00523 RT_line_pumping( t , lgShield_this_zone , DopplerWidth );
00524 }
00525 }
00526
00527 return;
00528 }
00529