cloudy  trunk
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
mc_escape.cpp
Go to the documentation of this file.
1 /* This file is part of Cloudy and is copyright (C)1978-2022 by Gary J. Ferland and
2  * others. For conditions of distribution and use see copyright notice in license.txt */
3 
4 #include "cddefines.h"
5 #include "mc_escape.h"
6 
7 #include "integrate.h"
8 #include "rt_escprob.h"
9 #include "thirdparty.h"
10 
11 namespace
12 {
13  static void mc_table(bool lgHeader, double tau, double a, double pesc1,
14  double nbar, double lbar, double taul, double pdest)
15  {
16  if (lgHeader)
17  {
18  // 1234567890---1234567890---1234567890---1234567890--
19  fprintf(ioQQQ,
20  "# tau1/2 pesc1 <N> <L>"
21  " <taul>"
22  " pdest esca0k2 esctot\n");
23  }
24  else
25  {
26  fprintf(ioQQQ,
27  "%12.5e %12.5e %12.5e %12.5e %12.5e %12.5e %12.5e %12.5e\n",
28  tau, pesc1, nbar, lbar, taul, pdest,esca0k2(tau/SQRTPI),
29  esc_CRDwing_1side(tau/SQRTPI,a));
30  }
31  fflush(ioQQQ);
32  }
33 
34  static double phigauss(double nu)
35  {
36  return exp(-nu*nu)/SQRTPI;
37  }
38 
39  // Implement equation (157) of Avrett & Loeser 1966
40  // 1966SAOSR.201.....A for comparison with escape probability with
41  // damping. Uses transformation x -> 1/y-1 to map domain of
42  // integral to (0,1]
43  class k2DampArg
44  {
45  realnum damp, tau;
46  public:
47  k2DampArg(realnum damp, realnum tau) : damp(damp), tau(tau) {}
48 
49  realnum operator()(realnum y) const
50  {
51  if (y <= 0.0)
52  return 0.0;
53  realnum x = 1./y-1.;
54  realnum phi;
55  VoigtU(damp,&x,&phi,1);
56  if (phi <= 0.)
57  {
58  return 0.;
59  }
60  else
61  {
62  return phi*expn(2,phi*tau)/POW2(y);
63  }
64  }
65  };
66 
67  class overlap
68  {
69  realnum damp, beta;
70  public:
71  overlap(realnum damp, realnum beta) : damp(damp), beta(beta) {}
72 
73  realnum operator()(realnum x) const
74  {
75  realnum phi;
76  if (0)
77  {
78  phi = (fabs(x)<0.5);
79  }
80  else
81  {
82  VoigtU(damp,&x,&phi,1);
83  }
84  return beta*phi/(beta+phi);
85  }
86  };
87 }
88 
89 
90 // Simple complete redistribution Monte Carlo
91 // tau_in is half-depth in mean optical depth scale
92 void mc_escape(double tau_in, double a, double beta)
93 {
94  // Parameters
95  // beta = k_c/k_L -- Hummer & Kunasz 1980
96 
97  fprintf(ioQQQ, "# a = %.4e, beta = %.4e\n", a, beta);
98  if (0)
99  {
100  // Confirm distributions
101  const long NSAMP = 1e7;
102  double totnu1 = 0.0, totnu2 = 0.0, totphi = 0.0;
103  for (long i=0; i<NSAMP; ++i)
104  {
105  double nu = RandGauss(0.,1./sqrt(2.));
106  totnu1 += nu*nu;
107  double nu2 = (10.*(i+0.5))/NSAMP;
108  double phi = phigauss(nu2);
109  totphi += phi;
110  totnu2 += nu2*nu2*phi;
111  }
112  // All these numbers should be 0.5 -- 1 and 3 are mean square of
113  // distributions, 2 is integral over half of profile
114  fprintf(ioQQQ,"Check %g %g %g\n",totnu1/NSAMP,totphi/(0.1*NSAMP),
115  totnu2/totphi);
116  }
117 
118  realnum width = 10.;
119  long npt = 1000;
120  vector<realnum> v(npt), y(npt), cv(npt+1);
121  if (a != 0.0)
122  {
123  // Construct cumulative pdf for generation of random deviates
124 
125  // Other distributions of ordinates will be better...
126  for (long i = 0; i<npt; ++i)
127  {
128  v[i] = width*(i+0.5)/npt;
129  }
130  VoigtU(a,&v[0],&y[0],npt);
131  cv[0] = 0.;
132  for (long i = 0; i<npt; ++i)
133  {
134  cv[i+1] = cv[i] + (width/npt)*y[i];
135  }
136  for (long i = 0; i<npt; ++i)
137  {
138  cv[i+1] /= cv[npt];
139  }
140  }
141 
142  const long NPART = 10000;
143  mc_table(true,0.,a,0.,0.,0.,0.,0.);
144  for (double tau0=0.01; tau0<tau_in; tau0 *= 1.1)
145  {
146  double nfirst = 0.0, ntot = 0.0, ltot=0.0, taulast=0.0,
147  ndest = 0.0;
148  for (long i=0; i<NPART; ++i)
149  {
150  double tau = tau0;
151  double nscat = 0;
152  for (;;)
153  {
154  // Save last scattering point
155  double tauprev = tau;
156  // CRD frequency, angle
157  double nu,phix;
158  if (a == 0.0)
159  {
160  nu = RandGauss(0.,1./sqrt(2.));
161  phix = phigauss(nu); // Line profile function
162  }
163  else
164  {
165  realnum offx = 2.0*genrand_real1()-1.0;
166  int sign = 1;
167  if (offx < 0.0)
168  {
169  sign = -1;
170  offx = -offx;
171  }
172  long ii = hunt_bisect(&cv[0],npt,offx)+1;
173  ASSERT (ii > 0 && ii < npt);
174  realnum frac = (offx-cv[ii])/SDIV(cv[ii+1]-cv[ii]);
175  nu = sign*((1.0-frac)*v[ii]+frac*v[ii+1]);
176  phix = (1.0-frac)*y[ii]+frac*y[ii+1];
177  }
178  double costheta = 2.0*genrand_real1()-1.0;
179 
180  // Random tau step at nu
181  double dtaunu = -log(genrand_real3());
182  // Convert tau step from nu to mean optical depth scale
183  // by dividing by relative asorption at nu, including
184  // continuous absorption
185  double dlength = dtaunu / SDIV(beta+phix);
186  // Convert to mean optical depth step perpendicular to slab
187  // cf HK equation 2.5, dtau ~ mu/(beta+phix)
188  double dtauperp = dlength * costheta;
189  // Calculate length of step within slab
190  if (tau-dtauperp < 0.0)
191  ltot += tau/costheta;
192  else if (tau-dtauperp > 2.0*tau0)
193  ltot += (2.0*tau0-tau)/costheta;
194  else
195  ltot += dlength;
196  // Calculate new position, reflecting across midplane
197  tau -= dtauperp;
198  if (tau > tau0)
199  tau = 2.0*tau0-tau;
200  if (tau < 0.0)
201  {
202  // Has escaped
203  taulast += tauprev;
204  if (nscat == 0)
205  ++nfirst;
206  break;
207  }
208  if (beta > 0.0 && genrand_real1() < beta/(beta+phix) )
209  {
210  // beta = k_c/k_L
211  // Destruction fraction is k_c/(k_c+k_L*phi)
212  ++ndest;
213  break;
214  }
215  // Has not escaped
216  ++nscat;
217  }
218  ntot += nscat;
219  }
220  // Print half-depth of slab, escape probability on midplane
221  // emission, mean number of scatterings and mean path length
222  // within slab
223  mc_table(false,tau0,a,nfirst/NPART,ntot/NPART,ltot/NPART,
224  taulast/(NPART-ndest),ndest/NPART);
225  }
226 
227  fprintf(ioQQQ,"\n\n#Avrett tau a esca0k2 esctot \n");
228 
229  // try to compare the formulae from esc_CRDwing_1side with
230  // Avrett & Loeser exact expression
231  double taustep = 2., taumin=1e-8,taumax=1e14;
232  for (double tau = taumin; tau < taumax; tau *= taustep)
233  {
234  double esccom_v = esca0k2(tau);
235  double esctot = esc_CRDwing_1side(tau,a);
236  k2DampArg k2fun(a,SQRTPI*tau);
237  double rmax1 = 1.0,rmax;
238  for(int idiv=0;idiv<200;++idiv)
239  {
240  rmax = rmax1;
241  rmax1 *= 0.75;
242  if (k2fun(rmax1) != 0.0)
243  {
244  break;
245  }
246  }
247  Integrator<k2DampArg,Gaussian32> IntDamp( k2fun );
248  double intgral = IntDamp.sum(0.0,rmax);
249  class integrate::Romberg<integrate::Midpoint<k2DampArg> >
250  k2r(integrate::Midpoint<k2DampArg>(k2fun,0.0,rmax));
251  k2r.update(1e-10);
252  fprintf(ioQQQ,"%15.8g %15.8g %15.8g %15.8g %15.8g %15.8g\n",
253  tau,a,esccom_v,esctot,2.0*intgral,2.0*k2r.sum());
254  }
255 
256  fprintf(ioQQQ,"\n\n#betaFbeta at a \n");
257  double betastep = 2., betamin=1e-6,betamax=1e3;
258  for (double beta = betamin; beta < betamax; beta *= betastep)
259  {
260  overlap ov(a,beta);
261  class integrate::Romberg<integrate::Midpoint<overlap> >
262  ovr(integrate::Midpoint<overlap>(ov,-100.0,100.0));
263  ovr.update(1e-10);
264  fprintf(ioQQQ,"%15.8g %15.8g %15.8g %15.8g\n",
265  beta,a,beta/(1.0+beta),ovr.sum());
266  }
267 }
double genrand_real1()
void update(double eps)
Definition: integrate.h:203
T sign(T x, T y)
Definition: cddefines.h:842
FILE * ioQQQ
Definition: cddefines.cpp:7
double expn(int n, double x)
long hunt_bisect(const T x[], long n, T xval)
Definition: thirdparty.h:303
void mc_escape(double tau_in, double a, double beta)
Definition: mc_escape.cpp:92
#define POW2
Definition: cddefines.h:979
float realnum
Definition: cddefines.h:124
double sum(double min, double max) const
Definition: integrate.h:44
double RandGauss(double xMean, double s)
Definition: service.cpp:1696
#define ASSERT(exp)
Definition: cddefines.h:613
void VoigtU(realnum a, const realnum v[], realnum y[], int n)
Definition: thirdparty.h:431
int fprintf(const Output &stream, const char *format,...)
Definition: service.cpp:1121
sys_float SDIV(sys_float x)
Definition: cddefines.h:1002
double esca0k2(double taume)
Definition: rt_escprob.cpp:424
double esc_CRDwing_1side(double tau, double a)
Definition: rt_escprob.cpp:165
double genrand_real3()
double frac(double d)