00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013 #include "cddefines.h"
00014 #include "physconst.h"
00015 #define SCALE 2.
00016 #include "dense.h"
00017 #include "conv.h"
00018 #include "rfield.h"
00019 #include "opacity.h"
00020 #include "lines_service.h"
00021 #include "taulines.h"
00022 #include "doppvel.h"
00023 #include "pressure.h"
00024 #include "wind.h"
00025 #include "rt.h"
00026
00027
00028
00029
00030 static double chnukt_ContTkt, chnukt_ctau;
00031
00032
00033 static double escmase(double tau);
00034
00035 static void RTesc_lya_1side(double taume,
00036 double beta,
00037 float *esc,
00038 float *dest,
00039
00040 long ipLine );
00041
00042 double esc_PRD_1side(double tau,
00043 double a)
00044 {
00045 double atau,
00046 b,
00047 escinc_v;
00048
00049 DEBUG_ENTRY( "esc_PRD_1side()" );
00050
00051
00052
00053
00054
00055
00056 # if 0
00057 if( strcmp(rt.chEscFunSubord,"SIMP") == 0 )
00058 {
00059
00060 escinc_v = 1./(1. + tau);
00061
00062 DEBUG_EXIT( "esc_PRD_1side()" );
00063 return( escinc_v );
00064 }
00065 # endif
00066
00067 if( tau < 0. )
00068 {
00069
00070 escinc_v = escmase(tau);
00071 }
00072 else if( tau < 10. )
00073 {
00074
00075 escinc_v = 1./(1. + 1.6*tau);
00076 }
00077 else
00078 {
00079
00080 atau = a*tau;
00081 if( atau > 1. )
00082 {
00083 b = 1.6 + (3.*pow(2.*a,-0.12))/(1. + atau);
00084 }
00085 else
00086 {
00087 b = 1.6 + (3.*pow(2.*a,-0.12))/(1. + 1./sqrt(atau));
00088 }
00089 b = MIN2(6.,b);
00090
00091 escinc_v = 1./(1. + b*tau);
00092 }
00093
00094 DEBUG_EXIT( "esc_PRD_1side()" );
00095 return( escinc_v );
00096 }
00097
00098
00099 double esc_CRDwing_1side(double tau,
00100 double a )
00101 {
00102 double esccom_v;
00103
00104 DEBUG_ENTRY( "esc_CRDwing_1side()" );
00105
00106
00107
00108
00109
00110
00111
00112
00113
00114
00115 esccom_v = esca0k2(tau);
00116 if( tau > 1e3 )
00117 {
00118 esccom_v += 0.333*sqrt(a/(SQRTPI*tau));
00119 }
00120
00121 DEBUG_EXIT( "esc_CRDwing_1side()" );
00122 return( esccom_v );
00123 }
00124
00125
00126
00127
00128 double RTesc_lya(
00129
00130 double *esin,
00131
00132 double *dest,
00133
00134 double abund,
00135
00136 long int nelem)
00137 {
00138 double beta,
00139 conopc,
00140 escla_v;
00141 float dstin,
00142 dstout;
00143
00144 DEBUG_ENTRY( "RTesc_lya()" );
00145
00146
00147
00148
00149
00150
00151
00152
00153
00154 ASSERT( nelem >= 0 && nelem < LIMELM );
00155
00156 if( EmisLines[ipH_LIKE][nelem][ipH2p][ipH1s].TauTot -
00157 EmisLines[ipH_LIKE][nelem][ipH2p][ipH1s].TauIn < 0. )
00158 {
00159
00160
00161 escla_v = EmisLines[ipH_LIKE][nelem][ipH2p][ipH1s].Pesc;
00162 rt.fracin = EmisLines[ipH_LIKE][nelem][ipH2p][ipH1s].FracInwd;
00163 *esin = rt.fracin;
00164 *dest = EmisLines[ipH_LIKE][nelem][ipH2p][ipH1s].Pdest;
00165
00166 DEBUG_EXIT( "RTesc_lya()" );
00167 return( escla_v );
00168 }
00169
00170
00171 conopc = opac.opacity_abs[EmisLines[ipH_LIKE][nelem][ipH2p][ipH1s].ipCont-1];
00172 if( abund > 0. )
00173 {
00174
00175 beta = conopc/(abund/SQRTPI*EmisLines[ipH_LIKE][nelem][ipH2p][ipH1s].opacity/
00176 DoppVel.doppler[nelem] + conopc);
00177 }
00178 else
00179 {
00180
00181 beta = 1e-10;
00182 }
00183
00184
00185 RTesc_lya_1side(
00186 EmisLines[ipH_LIKE][nelem][ipH2p][ipH1s].TauIn,
00187 beta,
00188 &rt.wayin,
00189 &dstin ,
00190
00191 EmisLines[ipH_LIKE][nelem][ipH2p][ipH1s].ipCont-1);
00192
00193 ASSERT( (rt.wayin <= 1.) && (rt.wayin >= 0.) && (dstin <= 1.) && (dstin >= 0.) );
00194
00195
00196 RTesc_lya_1side(MAX2(EmisLines[ipH_LIKE][nelem][ipH2p][ipH1s].TauTot/100.,
00197 EmisLines[ipH_LIKE][nelem][ipH2p][ipH1s].TauTot-EmisLines[ipH_LIKE][nelem][ipH2p][ipH1s].TauIn),
00198 beta,
00199 &rt.wayout,
00200 &dstout,
00201 EmisLines[ipH_LIKE][nelem][ipH2p][ipH1s].ipCont-1);
00202
00203 ASSERT( (rt.wayout <= 1.) && (rt.wayout >= 0.) && (dstout <= 1.) && (dstout >= 0.) );
00204
00205
00206 escla_v = (rt.wayin + rt.wayout)/2.;
00207
00208 *esin = rt.wayin;
00209
00210
00211 *dest = (dstin + dstout)/2.f;
00212
00213
00214
00215 *dest = (float)MIN2( *dest , 1.-escla_v );
00216
00217 *dest = (float)MAX2(0., *dest );
00218
00219
00220 rt.fracin = rt.wayin/(rt.wayin + rt.wayout);
00221 ASSERT( escla_v >=0. && *dest>=0. && *esin>=0. );
00222
00223 DEBUG_EXIT( "RTesc_lya()" );
00224 return( escla_v );
00225 }
00226
00227
00228 double esc_PRD(double tau,
00229 double tau_out,
00230 double damp )
00231 {
00232 double escgrd_v,
00233 tt;
00234
00235 DEBUG_ENTRY( "esc_PRD()" );
00236
00237
00238
00239 if( opac.lgTauOutOn )
00240 {
00241
00242 tt = tau_out - tau;
00243
00244
00245
00246
00247 if( tt < 0. )
00248 {
00249 tt = tau/SCALE;
00250 }
00251
00252 rt.wayin = (float)esc_PRD_1side(tau,damp);
00253 rt.wayout = (float)esc_PRD_1side(tt,damp);
00254 rt.fracin = rt.wayin/(rt.wayin + rt.wayout);
00255 escgrd_v = 0.5*(rt.wayin + rt.wayout);
00256 }
00257 else
00258 {
00259
00260 rt.fracin = 0.;
00261 rt.wayout = 1.;
00262 escgrd_v = esc_PRD_1side(tau,damp);
00263 rt.wayin = (float)escgrd_v;
00264 }
00265
00266 ASSERT( escgrd_v > 0. );
00267
00268 DEBUG_EXIT( "esc_PRD()" );
00269 return( escgrd_v );
00270 }
00271
00272
00273 double esc_CRDwing(double tau_in,
00274 double tau_out,
00275 double damp)
00276 {
00277 double escgrd_v,
00278 tt;
00279
00280 DEBUG_ENTRY( "esc_CRDwing()" );
00281
00282
00283
00284
00285 if( opac.lgTauOutOn )
00286 {
00287
00288
00289 if( tau_out <0 || tau_in < 0. )
00290 {
00291
00292 tt = MIN2( tau_out , tau_in );
00293 tau_in = tt;
00294 }
00295 else
00296 {
00297 tt = tau_out - tau_in;
00298
00299
00300
00301
00302 if( tt < 0. )
00303 {
00304 tt = tau_in/SCALE;
00305 }
00306 }
00307
00308 rt.wayin = (float)esc_CRDwing_1side(tau_in,damp);
00309 rt.wayout = (float)esc_CRDwing_1side(tt,damp);
00310 rt.fracin = rt.wayin/(rt.wayin + rt.wayout);
00311 escgrd_v = 0.5*(rt.wayin + rt.wayout);
00312 }
00313 else
00314 {
00315
00316 rt.fracin = 0.;
00317 rt.wayout = 1.;
00318 escgrd_v = esc_CRDwing_1side(tau_in,damp);
00319 rt.wayin = (float)escgrd_v;
00320 }
00321
00322 ASSERT( escgrd_v > 0. );
00323
00324 DEBUG_EXIT( "esc_CRDwing()" );
00325 return( escgrd_v );
00326 }
00327
00328
00329 double esc_CRDcore(double tau_in,
00330 double tau_out)
00331 {
00332 double escgrd_v,
00333 tt;
00334
00335 DEBUG_ENTRY( "esc_CRDcore()" );
00336
00337
00338
00339
00340 if( opac.lgTauOutOn )
00341 {
00342
00343
00344 if( tau_out <0 || tau_in < 0. )
00345 {
00346
00347 tt = MIN2( tau_out , tau_in );
00348 tau_in = tt;
00349 }
00350 else
00351 {
00352 tt = tau_out - tau_in;
00353
00354
00355
00356
00357 if( tt < 0. )
00358 {
00359 tt = tau_in/SCALE;
00360 }
00361 }
00362
00363 rt.wayin = (float)esca0k2(tau_in);
00364 rt.wayout = (float)esca0k2(tt);
00365 rt.fracin = rt.wayin/(rt.wayin + rt.wayout);
00366 escgrd_v = 0.5*(rt.wayin + rt.wayout);
00367 }
00368 else
00369 {
00370
00371 rt.fracin = 0.;
00372 rt.wayout = 1.;
00373 escgrd_v = esca0k2(tau_in);
00374 rt.wayin = (float)escgrd_v;
00375 }
00376
00377 ASSERT( escgrd_v > 0. );
00378
00379 DEBUG_EXIT( "esc_CRDcore()" );
00380 return( escgrd_v );
00381 }
00382
00383
00384 double esca0k2(double taume)
00385 {
00386 double arg,
00387 esca0k2_v,
00388 suma,
00389 sumb,
00390 sumc,
00391 sumd,
00392 tau;
00393 static double a[5]={1.00,-0.1117897,-0.1249099917,-9.136358767e-3,
00394 -3.370280896e-4};
00395 static double b[6]={1.00,0.1566124168,9.013261660e-3,1.908481163e-4,
00396 -1.547417750e-7,-6.657439727e-9};
00397 static double c[5]={1.000,19.15049608,100.7986843,129.5307533,-31.43372468};
00398 static double d[6]={1.00,19.68910391,110.2576321,169.4911399,-16.69969409,
00399 -36.664480000};
00400
00401 DEBUG_ENTRY( "esca0k2()" );
00402
00403
00404
00405
00406
00407
00408 tau = taume*SQRTPI;
00409
00410 if( tau < 0. )
00411 {
00412
00413 esca0k2_v = escmase(taume);
00414
00415 }
00416 else if( tau < 0.01 )
00417 {
00418 esca0k2_v = 1. - 2.*tau;
00419
00420 }
00421 else if( tau <= 11. )
00422 {
00423 suma = a[0] + tau*(a[1] + tau*(a[2] + tau*(a[3] + a[4]*tau)));
00424 sumb = b[0] + tau*(b[1] + tau*(b[2] + tau*(b[3] + tau*(b[4] +
00425 b[5]*tau))));
00426 esca0k2_v = tau/2.5066283*log(tau/SQRTPI) + suma/sumb;
00427
00428 }
00429 else
00430 {
00431
00432 arg = 1./log(tau/SQRTPI);
00433 sumc = c[0] + arg*(c[1] + arg*(c[2] + arg*(c[3] + c[4]*arg)));
00434 sumd = d[0] + arg*(d[1] + arg*(d[2] + arg*(d[3] + arg*(d[4] +
00435 d[5]*arg))));
00436 esca0k2_v = (sumc/sumd)/(2.*tau*sqrt(log(tau/SQRTPI)));
00437 }
00438
00439 DEBUG_EXIT( "esca0k2()" );
00440 return( esca0k2_v );
00441 }
00442
00443
00444 static void FindNeg( void )
00445 {
00446 long int i;
00447
00448 DEBUG_ENTRY( "FindNeg()" );
00449
00450
00451 for( i=1; i <= nLevel1; i++ )
00452 {
00453
00454 if( TauLines[i].TauIn < -1. )
00455 DumpLine(&TauLines[i]);
00456 }
00457
00458
00459 for( i=0; i < nWindLine; i++ )
00460 {
00461 if( TauLine2[i].IonStg < TauLine2[i].nelem+1-NISO )
00462 {
00463
00464 if( TauLine2[i].TauIn < -1. )
00465 DumpLine(&TauLine2[i]);
00466 }
00467 }
00468
00469
00470 for( i=0; i < nHFLines; i++ )
00471 {
00472
00473 if( HFLines[i].TauIn < -1. )
00474 DumpLine(&HFLines[i]);
00475 }
00476
00477
00478 for( i=0; i < nCORotate; i++ )
00479 {
00480
00481 if( C12O16Rotate[i].TauIn < -1. )
00482 DumpLine(&C12O16Rotate[i]);
00483
00484 if( C13O16Rotate[i].TauIn < -1. )
00485 DumpLine(&C13O16Rotate[i]);
00486 }
00487
00488 DEBUG_EXIT( "FindNeg()" );
00489 return;
00490 }
00491
00492 static double escmase(double tau)
00493 {
00494 double escmase_v;
00495
00496 DEBUG_ENTRY( "escmase()" );
00497
00498
00499 ASSERT( tau <= 0. );
00500
00501 if( tau > -0.1 )
00502 {
00503 escmase_v = 1. - tau*(0.5 + tau/6.);
00504 }
00505 else if( tau > -30. )
00506 {
00507 escmase_v = (1. - exp(-tau))/tau;
00508 }
00509 else
00510 {
00511 fprintf( ioQQQ, " DISASTER escmase called with 2big tau%10.2e\n",
00512 tau );
00513 fprintf( ioQQQ, " This is zone number%4ld\n", nzone );
00514 FindNeg();
00515 ShowMe();
00516 puts( "[Stop in escmase]" );
00517 cdEXIT(EXIT_FAILURE);
00518 }
00519
00520 ASSERT( escmase_v >= 1. );
00521
00522 DEBUG_EXIT( "escmase()" );
00523 return( escmase_v );
00524 }
00525
00526
00527
00528 static double cone2(double t);
00529
00530 double escConE2(double x)
00531 {
00532 double conesc_v;
00533
00534 DEBUG_ENTRY( "escConE2()" );
00535
00536 conesc_v = exp(-chnukt_ContTkt*(x-1.))/
00537 x*cone2(chnukt_ctau/x/x/x);
00538
00539 DEBUG_EXIT( "escConE2()" );
00540 return( conesc_v );
00541 }
00542
00543
00544 static double cone2(double t)
00545 {
00546 double cone2_v,
00547 remain,
00548 tln;
00549
00550 DEBUG_ENTRY( "cone2()" );
00551
00552 if( t < 80. )
00553 {
00554 tln = exp(-t);
00555 }
00556 else
00557 {
00558 cone2_v = 0.;
00559
00560 DEBUG_EXIT( "cone2()" );
00561 return( cone2_v );
00562 }
00563
00564
00565
00566
00567 if( t < 0.3 )
00568 {
00569 remain = (1.998069357 + t*(66.4037741 + t*107.2041376))/(1. +
00570 t*(37.4009646 + t*105.0388805));
00571
00572 }
00573 else if( t < 20. )
00574 {
00575 remain = (1.823707708 + t*2.395042899)/(1. + t*(2.488885899 -
00576 t*0.00430538));
00577
00578 }
00579 else
00580 {
00581 remain = 1.;
00582 }
00583
00584 cone2_v = remain*tln/(2. + t);
00585
00586 DEBUG_EXIT( "cone2()" );
00587 return( cone2_v );
00588 }
00589
00590
00591 static double conrec(double x)
00592 {
00593 double conrec_v;
00594
00595 DEBUG_ENTRY( "conrec()" );
00596
00597 conrec_v = exp(-chnukt_ContTkt*(x-1.))/x;
00598
00599 DEBUG_EXIT( "conrec()" );
00600 return( conrec_v );
00601 }
00602
00603
00604 double esccon(double tau,
00605 double hnukt)
00606 {
00607 double dinc,
00608 escpcn_v,
00609 sumesc,
00610 sumrec;
00611
00612 DEBUG_ENTRY( "esccon()" );
00613
00614
00615 if( tau < 0.01 )
00616 {
00617 escpcn_v = 1.;
00618
00619 DEBUG_EXIT( "esccon()" );
00620 return( escpcn_v );
00621 }
00622
00623 else if( hnukt > 1. && tau > 100. )
00624 {
00625 escpcn_v = 1e-20;
00626
00627 DEBUG_EXIT( "esccon()" );
00628 return( escpcn_v );
00629 }
00630
00631 chnukt_ContTkt = hnukt;
00632 chnukt_ctau = tau;
00633
00634 dinc = 10./hnukt;
00635 sumrec = qg32(1.,1.+dinc,conrec);
00636 sumesc = qg32(1.,1.+dinc,escConE2);
00637
00638 if( sumrec > 0. )
00639 {
00640 escpcn_v = sumesc/sumrec;
00641 }
00642 else
00643 {
00644 escpcn_v = 0.;
00645 }
00646
00647 DEBUG_EXIT( "esccon()" );
00648 return( escpcn_v );
00649 }
00650
00651
00652 static void RTesc_lya_1side(double taume,
00653 double beta,
00654 float *esc,
00655 float *dest,
00656
00657 long ipLine )
00658 {
00659 double esc0,
00660 fac,
00661 fac1,
00662 fac2,
00663 tau,
00664 taucon,
00665 taulog;
00666
00667
00668
00669
00670
00671 DEBUG_ENTRY( "RTesc_lya_1side()" );
00672
00673
00674 tau = taume*SQRTPI;
00675
00676
00677 esc0 = 1./((0.6451 + tau)*(0.47 + 1.08/(1. + 7.3e-6*tau)));
00678
00679 esc0 = MAX2(0.,esc0);
00680 esc0 = MIN2(1.,esc0);
00681
00682 if( tau > 0. )
00683 {
00684 taulog = log10(MIN2(1e8,tau));
00685 }
00686 else
00687 {
00688
00689
00690
00691 taulog = 0.;
00692 *dest = 0.;
00693 *esc = (float)esc0;
00694 }
00695
00696 if( beta > 0. )
00697 {
00698 taucon = MIN2(2.,beta*tau);
00699
00700 if( taucon > 1e-3 )
00701 {
00702 fac1 = -1.25 + 0.475*taulog;
00703 fac2 = -0.485 + 0.1615*taulog;
00704 fac = -fac1*taucon + fac2*POW2(taucon);
00705 fac = pow(10.,fac);
00706 fac = MIN2(1.,fac);
00707 }
00708 else
00709 {
00710 fac = 1.;
00711 }
00712
00713 *esc = (float)(esc0*fac);
00714
00715 *dest = (float)(beta/(0.30972 - MIN2(.28972,0.03541667*taulog)));
00716 }
00717
00718 else
00719 {
00720 *dest = 0.;
00721 *esc = (float)esc0;
00722 }
00723
00724 *dest = MIN2(*dest,1.f-*esc);
00725 *dest = MAX2(0.f,*dest);
00726
00727
00728
00729 *dest = (float)( (1. - opac.albedo[ipLine]) * *dest + opac.albedo[ipLine]*DEST0);
00730
00731 {
00732
00733 enum {BUG=false};
00734
00735 if( BUG )
00736 {
00737 fprintf(ioQQQ,"scatdest tau %.2e beta%.2e 1-al%.2e al%.2e dest%.2e \n",
00738 taume,
00739 beta,
00740 (1. - opac.albedo[ipLine]),
00741 opac.albedo[ipLine] ,
00742 *dest
00743 );
00744 }
00745 }
00746
00747 DEBUG_EXIT( "RTesc_lya_1side()" );
00748 return;
00749 }
00750
00751
00752 double RT_DestProb(
00753
00754 double abund,
00755
00756 double crsec,
00757
00758
00759 long int ipanu,
00760
00761 double widl,
00762
00763 double escp,
00764
00765 int nCore)
00766 {
00767
00768 double eovrlp_v;
00769
00770 double conopc,
00771 beta;
00772
00773
00774
00775
00776
00777
00778
00779
00780
00781
00782
00783
00784
00785
00786 DEBUG_ENTRY( "RT_DestProb()" );
00787
00788
00789
00790
00791
00792
00793
00794
00795 if( escp >= 1.0 || !conv.nTotalIoniz || ipanu >= rfield.nflux )
00796 {
00797 eovrlp_v = 0.;
00798
00799 DEBUG_EXIT( "RT_DestProb()" );
00800 return( eovrlp_v );
00801 }
00802
00803
00804 conopc = opac.opacity_abs[ipanu-1];
00805
00806 ASSERT( crsec > 0. );
00807
00808
00809 if( abund <= 0. || conopc <= 0. )
00810 {
00811
00812 eovrlp_v = 0.;
00813
00814 DEBUG_EXIT( "RT_DestProb()" );
00815 return( eovrlp_v );
00816 }
00817
00818
00819 beta = conopc/(abund*SQRTPI*crsec/widl + conopc);
00820
00821
00822
00823 if( nCore == ipDEST_INCOM )
00824 {
00825
00826
00827
00828
00829
00832 eovrlp_v = MIN2(1e-3,8.5*beta);
00833 }
00834 else if( nCore == ipDEST_K2 )
00835 {
00836
00837
00838 eovrlp_v = MIN2(1e-3,8.5*beta);
00839 }
00840 else if( nCore == ipDEST_SIMPL )
00841 {
00842
00843
00844
00845 eovrlp_v = MIN2(1e-3,8.5*beta);
00846 }
00847 else
00848 {
00849 fprintf( ioQQQ, " chCore of %i not understood by RT_DestProb.\n",
00850 nCore );
00851 puts( "[Stop in RT_DestProb]" );
00852 cdEXIT(EXIT_FAILURE);
00853 }
00854
00855
00856 eovrlp_v /= 1. + eovrlp_v;
00857
00858
00859 eovrlp_v *= 1. - escp;
00860
00861
00862 ASSERT( eovrlp_v >= 0. );
00863 ASSERT( eovrlp_v <= 1. );
00864
00865 {
00866
00867
00868 enum {DEBUG_LOC=false};
00869
00870 if( DEBUG_LOC )
00871 {
00872
00873 if( rfield.anu[ipanu-1]>0.73 && rfield.anu[ipanu-1]<0.76 &&
00874 abund==EmisLines[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].PopOpc )
00875
00876 {
00877 fprintf(ioQQQ,"%li RT_DestProb\t%g\n",
00878 nzone, eovrlp_v );
00879 }
00880 }
00881 }
00882
00883
00884
00885
00886 # if 0
00887
00888
00889
00890
00891
00892
00893
00894 eovrlp_v = POW2(1. - opac.albedo[ipanu-1]) * eovrlp_v +
00895 opac.albedo[ipanu-1]*0.;
00896 # endif
00897
00898 DEBUG_EXIT( "RT_DestProb()" );
00899
00900 return( eovrlp_v );
00901 }
00902
00903 #if 0
00904
00905
00906 double RT_LyaWidth(
00907 double tauin,
00908 double tauout,
00909 double a,
00910
00911
00912 double vth)
00913 {
00914 static bool lgTonSav;
00915 double aa,
00916 atau,
00917 b,
00918 r,
00919 t,
00920 tau,
00921 therm,
00922 widlin_v;
00923
00924 DEBUG_ENTRY( "RT_LyaWidth()" );
00925
00926
00927 lgTonSav = opac.lgTauOutOn;
00928 opac.lgTauOutOn = true;
00929
00930
00931
00932
00933
00934
00935
00936
00937
00938
00939 if( opac.lgTauOutOn )
00940 {
00941 tau = tauout/2.;
00942 t = MIN2(tauin,tauout-tauin);
00943 t = MAX2(tauin/100.,t);
00944 }
00945 else
00946 {
00947 tau = tauin;
00948 t = tauin;
00949 }
00950
00951
00952 therm = 5.3e16/dense.eden;
00953 if( tau > therm )
00954 {
00955 pressure.lgPradDen = true;
00956 tau = therm;
00957 }
00958
00959
00963 if( wind.windv <= 0 )
00964 {
00965
00966 atau = log(MAX2(0.0001,tau));
00967 if( tau <= 20. )
00968 {
00969 aa = 4.8 + 5.2*tau + (4.*tau - 1.)*atau;
00970 b = 6.5*tau - atau;
00971 }
00972 else
00973 {
00974 aa = 1. + 2.*atau/pow(1. + 0.3*a*tau,0.6667) + pow(6.5*
00975 a*tau,0.333);
00976 b = 1.6 + 1.5/(1. + 0.20*a*tau);
00977 }
00978
00979
00980 widlin_v = vth*0.8862*aa/b*(1. - esc_PRD_1side(t,a));
00981
00982 }
00983
00984 else
00985 {
00986
00987 r = a*tau/PI;
00988 if( r <= 1. )
00989 {
00990 widlin_v = vth*sqrt(log(MAX2(1.,tau))*.2821);
00991
00992 ASSERT( widlin_v >= 0. );
00993 DEBUG_EXIT( "RT_LyaWidth()" );
00994 return( widlin_v );
00995 }
00996
00997 widlin_v = 2.*fabs(wind.windv0);
00998 if( r*vth > fabs(widlin_v ) )
00999 {
01000 ASSERT( widlin_v >= 0. );
01001 DEBUG_EXIT( "RT_LyaWidth()" );
01002 return( widlin_v );
01003 }
01004
01005 widlin_v = vth*r*log(widlin_v/(r*vth));
01006 }
01007
01008
01009 opac.lgTauOutOn = lgTonSav;
01010
01011 ASSERT( widlin_v >= 0. );
01012
01013 DEBUG_EXIT( "RT_LyaWidth()" );
01014 return( widlin_v );
01015 }
01016 #endif
01017
01018
01019
01020 double RT_LineWidth(EmLine * t)
01021 {
01022 double RT_LineWidth_v,
01023 aa,
01024 atau,
01025 b,
01026 r,
01027 tau,
01028 tauin,
01029 tauout,
01030 therm,
01031 vth;
01032
01033 DEBUG_ENTRY( "RT_LineWidth()" );
01034
01035
01036
01037
01038
01039
01040
01041
01042
01043
01044 tauin = t->TauIn;
01045 vth = DoppVel.doppler[ t->nelem -1 ];
01046
01047 if( opac.lgTauOutOn )
01048 {
01049
01050
01051
01052
01053 tauout = t->TauTot-tauin;
01054
01055
01056 tau = MIN2( tauin , tauout );
01057 }
01058 else
01059 {
01060 tau = tauin;
01061 }
01062
01063
01064 therm = 5.3e16/dense.eden;
01065
01066 if( tau > therm )
01067 {
01068 pressure.lgPradDen = true;
01069 tau = therm;
01070 }
01071
01072
01073
01074
01075
01076
01077
01082 if( wind.windv <= 0. )
01083 {
01084
01085 # if 0
01086
01087
01088
01089 if( 0 && tau< 0.0001 )
01090 {
01091
01092 aa = 14.;
01093 b = 9.21;
01094
01095
01096 RT_LineWidth_v = vth*0.8862*aa/b*MAX2(0.,tau-opac.taumin);
01097 }
01098 # endif
01099
01100 if( (tau-opac.taumin)/100. < FLT_EPSILON )
01101 {
01102 RT_LineWidth_v = 0.;
01103 }
01104 else if( tau <= 20. )
01105 {
01106 atau = log(MAX2(0.0001,tau));
01107 aa = 4.8 + 5.2*tau + (4.*tau - 1.)*atau;
01108 b = 6.5*tau - atau;
01109
01110
01111 RT_LineWidth_v = vth*0.8862*aa/b*(1. - MIN2( (t->Pesc+ t->Pelec_esc + t->Pdest) , 1.) );
01112
01113 if( t->Pesc < 10. * FLT_EPSILON )
01114 RT_LineWidth_v = 0.;
01115 }
01116 else
01117 {
01118 ASSERT( t->damp*tau >= 0.);
01119 atau = log(MAX2(0.0001,tau));
01120 aa = 1. + 2.*atau/pow(1. + 0.3*t->damp*tau,0.6667) + pow(6.5*
01121 t->damp*tau,0.333);
01122 b = 1.6 + 1.5/(1. + 0.20*t->damp*tau);
01123
01124
01125 RT_LineWidth_v = vth*0.8862*aa/b*(1. - MIN2( (t->Pesc+ t->Pelec_esc + t->Pdest) , 1.) );
01126 }
01127
01128
01129 RT_LineWidth_v *= 2.;
01130
01131 }
01132 else
01133 {
01134
01135 r = t->damp*tau/PI;
01136 if( r <= 1. )
01137 {
01138 RT_LineWidth_v = vth*sqrt(log(MAX2(1.,tau))*.2821);
01139 }
01140 else
01141 {
01142 RT_LineWidth_v = 2.*fabs(wind.windv0);
01143 if( r*vth <= RT_LineWidth_v )
01144 {
01145 RT_LineWidth_v = vth*r*log(RT_LineWidth_v/(r*vth));
01146 }
01147 }
01148 }
01149
01150 ASSERT( RT_LineWidth_v >= 0. );
01151
01152 DEBUG_EXIT( "RT_LineWidth()" );
01153 return( RT_LineWidth_v );
01154 }
01155
01156
01157 double RT_DestHummer(double beta)
01158
01159 {
01160 double fhummr_v,
01161 x;
01162
01163 DEBUG_ENTRY( "RT_DestHummer()" );
01164
01165
01166
01167
01168
01169
01170
01171
01172 ASSERT( beta >= 0.);
01173 if( beta <= 0. )
01174 {
01175 fhummr_v = 0.;
01176 }
01177 else
01178 {
01179 x = log10(beta);
01180 if( x < -5.5 )
01181 {
01182 fhummr_v = 3.8363 - 0.56329*x;
01183 }
01184 else if( x < -3.5 )
01185 {
01186 fhummr_v = 2.79153 - 0.75325*x;
01187 }
01188 else if( x < -2. )
01189 {
01190 fhummr_v = 1.8446 - 1.0238*x;
01191 }
01192 else
01193 {
01194 fhummr_v = 0.72500 - 1.5836*x;
01195 }
01196 fhummr_v *= beta;
01197 }
01198
01199 DEBUG_EXIT( "RT_DestHummer()" );
01200 return( fhummr_v );
01201 }
01202