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