30 const double SCALE = 2.;
31 double tt = tau_tot - tau_in;
112 if( strcmp(
rt.chEscFunSubord,
"SIMP") == 0 )
115 escinc_v = 1./(1. + tau);
131 b = 1.6 + (3.*pow(2.*a,-0.12))/(1. + atau);
135 double sqrtatau = sqrt(atau);
136 b = 1.6 + (3.*pow(2.*a,-0.12))*sqrtatau/(1. + sqrtatau);
140 escinc_v = 1./(1. + b*tau);
177 vector<double> vxmin, vxmax, fxmin, fxmax, step;
178 vxmin.push_back(xmin);
179 vxmax.push_back(xmax);
180 fxmin.push_back(func(xmin));
181 fxmax.push_back(func(xmax));
182 step.push_back(0.5*(xmax-xmin)*(fxmin[0]+fxmax[0]));
183 for (
long level = 0; level < 25; ++level)
185 size_t nstep = step.size();
186 double current = 0.0;
187 for (
size_t i=0; i<nstep; ++i)
189 for (
size_t i=0; i<nstep; ++i)
191 if (vxmax[i] > vxmin[i])
193 if (level <= 3 || step[i] > tol*current
194 || (i == nstep-1 && current <= 0.0) )
196 double xbar=0.5*(vxmin[i]+vxmax[i]);
197 vxmin.push_back(xbar);
198 fxmin.push_back(func(xbar));
199 vxmax.push_back(vxmax[i]);
200 fxmax.push_back(fxmax[i]);
201 long ilast = vxmax.size()-1;
203 fxmax[i] = fxmin[ilast];
204 step[i] = 0.5*(vxmax[i]-vxmin[i])*(fxmin[i]+fxmax[i]);
207 step.push_back(0.5*(vxmax[ilast]-vxmin[ilast])*
208 (fxmin[ilast]+fxmax[ilast]));
213 double current = 0.0;
214 for (
size_t i=0; i<step.size(); ++i)
241 double taustep = 2., taumin=1e-8,taumax=1e14;
243 for (tau = taumin; tau < taumax; tau *= taustep)
245 double esccom_v =
esca0k2(tau);
246 double sqrta = sqrt(a);
247 double pwing = tau > 0.0 ? sqrta/(sqrta+0.5*3.0*sqrt(
SQRTPI*tau)) :
249 double esctot = esccom_v*(1.0-pwing)+pwing;
252 fprintf(
ioQQQ,
"cfesc %15.8g %15.8g %15.8g %15.8g %15.8g %15.8g\n",
253 tau,a,esccom_v,esctot,2.0*intgral,2.0*intgralt);
257 double esccom_v =
esca0k2(tau);
262 double sqrta = sqrt(a);
263 double scal = a*(1.0+a+tau)/(
POW2(1.0+a)+a*tau);
264 double pwing = scal*((tau > 0.0) ? sqrta/sqrt(a+2.25*
SQRTPI*tau) : 1.0);
265 return esccom_v*(1.0-pwing)+pwing;
313 DopplerWidth + conopc);
330 ASSERT( (
rt.wayin <= 1.) && (
rt.wayin >= 0.) && (dstin <= 1.) && (dstin >= 0.) );
340 ASSERT( (
rt.wayout <= 1.) && (
rt.wayout >= 0.) && (dstout <= 1.) && (dstout >= 0.) );
343 escla_v = (
rt.wayin +
rt.wayout)/2.;
348 *dest = (dstin + dstout)/2.f;
357 rt.fracin =
rt.wayin/(
rt.wayin +
rt.wayout);
358 ASSERT( escla_v >=0. && *dest>=0. && *esin>=0. );
382 rt.fracin =
rt.wayin/(
rt.wayin +
rt.wayout);
383 escgrd_v = 0.5*(
rt.wayin +
rt.wayout);
415 if( tau_out <0 || tau_in < 0. )
418 tt =
MIN2( tau_out , tau_in );
428 rt.fracin =
rt.wayin/(
rt.wayin +
rt.wayout);
429 escgrd_v = 0.5*(
rt.wayin +
rt.wayout);
460 if( tau_out <0 || tau_in < 0. )
463 tt =
MIN2( tau_out , tau_in );
473 rt.fracin =
rt.wayin/(
rt.wayin +
rt.wayout);
474 escgrd_v = 0.5*(
rt.wayin +
rt.wayout);
499 static const double a[5]={1.00,-0.1117897,-0.1249099917,-9.136358767e-3,
501 static const double b[6]={1.00,0.1566124168,9.013261660e-3,1.908481163e-4,
502 -1.547417750e-7,-6.657439727e-9};
503 static const double c[5]={1.000,19.15049608,100.7986843,129.5307533,-31.43372468};
504 static const double d[6]={1.00,19.68910391,110.2576321,169.4911399,-16.69969409,
522 else if( tau < 0.01 )
524 esca0k2_v = 1. - 2.*tau;
527 else if( tau <= 11. )
529 suma = a[0] + tau*(a[1] + tau*(a[2] + tau*(a[3] + a[4]*tau)));
530 sumb = b[0] + tau*(b[1] + tau*(b[2] + tau*(b[3] + tau*(b[4] +
532 esca0k2_v = tau/2.5066283*log(tau/
SQRTPI) + suma/sumb;
539 sumc = c[0] + arg*(c[1] + arg*(c[2] + arg*(c[3] + c[4]*arg)));
540 sumd = d[0] + arg*(d[1] + arg*(d[2] + arg*(d[3] + arg*(d[4] +
542 esca0k2_v = (sumc/sumd)/(2.*tau*sqrt(log(tau/
SQRTPI)));
558 if(
TauLines[i].Emis().TauIn() < -1. )
564 for (
int ipSpecies=0; ipSpecies <
nSpecies; ++ipSpecies)
567 em !=
dBaseTrans[ipSpecies].Emis().end(); ++em)
569 if((*em).TauIn() < -1. )
580 if(
TauLine2[i].Emis().TauIn() < -1. )
589 if(
HFLines[i].Emis().TauIn() < -1. )
607 escmase_v = 1. - tau*(0.5 + tau/6.);
609 else if( tau > -30. )
611 escmase_v = (1. - exp(-tau))/tau;
615 fprintf(
ioQQQ,
" DISASTER escmase called with 2big tau%10.2e\n",
617 fprintf(
ioQQQ,
" This is zone number%4ld\n",
nzone );
623 ASSERT( escmase_v >= 1. );
645 else if( hnukt > 1. && tau > 100. )
661 sumrec = conrec.
sum(1.,1.+dinc,func_conrec);
662 sumesc = escConE2.
sum(1.,1.+dinc,func_escConE2);
666 escpcn_v = sumesc/sumrec;
701 esc0 = 1./((0.6451 + tau)*(0.47 + 1.08/(1. + 7.3e-6*tau)));
703 esc0 =
MAX2(0.,esc0);
704 esc0 =
MIN2(1.,esc0);
708 taulog = log10(
MIN2(1e8,tau));
722 taucon =
MIN2(2.,beta*tau);
726 fac1 = -1.25 + 0.475*taulog;
727 fac2 = -0.485 + 0.1615*taulog;
728 fac = -fac1*taucon + fac2*
POW2(taucon);
739 *dest = (
realnum)(beta/(0.30972 -
MIN2(.28972,0.03541667*taulog)));
748 *dest =
MIN2(*dest,1.f-*esc);
749 *dest =
MAX2(0.f,*dest);
761 fprintf(
ioQQQ,
"scatdest tau %.2e beta%.2e 1-al%.2e al%.2e dest%.2e \n",
817 if( escp >= 1.0 || !
conv.nTotalIoniz || ipanu >=
rfield.nflux )
824 conopc =
opac.opacity_abs[ipanu-1];
829 if(
abund <= 0. || conopc <= 0. )
850 eovrlp_v =
MIN2(1e-3,8.5*beta);
856 eovrlp_v =
MIN2(1e-3,8.5*beta);
863 eovrlp_v =
MIN2(1e-3,8.5*beta);
867 fprintf(
ioQQQ,
" chCore of %i not understood by RT_DestProb.\n",
873 eovrlp_v /= 1. + eovrlp_v;
876 eovrlp_v *= 1. - escp;
885 enum {DEBUG_LOC=
false};
889 if(
rfield.anu[ipanu-1]>0.73 &&
rfield.anu[ipanu-1]<0.76 &&
892 fprintf(
ioQQQ,
"%li RT_DestProb\t%g\n",
909 eovrlp_v =
POW2(1. -
opac.albedo[ipanu-1]) * eovrlp_v +
910 opac.albedo[ipanu-1]*0.;
920 double RT_LineWidth_v,
980 if( !
wind.lgBallistic() )
984 if( (tau-
opac.taumin)/100. < FLT_EPSILON )
988 else if( tau <= 20. )
993 aa = 4.8 + 5.2*tau + (4.*tau - 1.)*atau;
996 RT_LineWidth_v = vth*0.8862*aa/b*(1. -
MIN2( 1., escProb ) );
998 if( escProb >= 1. - 100. * FLT_EPSILON )
1004 atau = log(
MAX2(0.0001,tau));
1005 aa = 1. + 2.*atau/pow(1. + 0.3*t.
Emis().
damp()*tau,0.6667) + pow(6.5*
1007 b = 1.6 + 1.5/(1. + 0.20*t.
Emis().damp()*tau);
1008 RT_LineWidth_v = vth*0.8862*aa/b*(1. -
MIN2( 1. ,
1013 RT_LineWidth_v *= 2.;
1022 RT_LineWidth_v = vth*sqrt(log(
MAX2(1.,tau))*.2821);
1026 RT_LineWidth_v = 2.*fabs(
wind.windv0);
1027 if( r*vth <= RT_LineWidth_v )
1029 RT_LineWidth_v = vth*r*log(RT_LineWidth_v/(r*vth));
1034 ASSERT( RT_LineWidth_v >= 0. );
1035 return RT_LineWidth_v;
1064 fhummr_v = 3.8363 - 0.56329*x;
1068 fhummr_v = 2.79153 - 0.75325*x;
1072 fhummr_v = 1.8446 - 1.0238*x;
1076 fhummr_v = 0.72500 - 1.5836*x;
bool fp_equal(sys_float x, sys_float y, int n=3)
#define DEBUG_ENTRY(funcname)
EmissionProxy::iterator iterator
realnum & Pelec_esc() const
realnum & FracInwd() const
realnum & dampXvel() const
double sum(double min, double max, Integrand func)
EmissionList::reference Emis() const
realnum operator()(realnum y) const
k2DampArg(realnum damp, realnum tau)
double operator()(double x)
double operator()(double x)
t_iso_sp iso_sp[NISO][LIMELM]
UNUSED const double SQRTPI
double RT_DestHummer(double beta)
double tau_from_here(double tau_tot, double tau_in)
double esc_PRD_1side(double tau, double a)
double esc_CRDcore(double tau_in, double tau_out)
double trapezium(realnum xmin, realnum xmax, const F func)
double esca0k2(double taume)
double RT_LineWidth(const TransitionProxy &t, realnum DopplerWidth)
STATIC void FindNeg(void)
double RT_DestProb(double abund, double crsec, long int ipanu, double widl, double escp, int nCore)
STATIC double escmase(double tau)
double esc_CRDwing_1side(double tau, double a)
double esc_CRDwing(double tau_in, double tau_out, double damp)
double RTesc_lya(double *esin, double *dest, double abund, const TransitionProxy &t, realnum DopplerWidth)
double esc_PRD(double tau, double tau_out, double damp)
double esccon(double tau, double hnukt)
STATIC void RTesc_lya_1side(double taume, double beta, realnum *esc, realnum *dest, long ipLine)
TransitionList TauLine2("TauLine2", &AnonStates)
vector< TransitionList > dBaseTrans
TransitionList HFLines("HFLines", &AnonStates)
TransitionList TauLines("TauLines", &AnonStates)
double expn(int n, double x)
void VoigtU(realnum a, const realnum v[], realnum y[], int n)
void DumpLine(const TransitionProxy &t)