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