cloudy  trunk
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
rt_continuum_shield_fcn.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 /*rt_continuum_shield_fcn computing continuum shielding due to single line */
4 /*conpmp local continuum pumping rate radiative transfer for all lines */
5 /* This file is part of Cloudy and is copyright (C)1978-2022 by Gary J. Ferland and
6  * others. For conditions of distribution and use see copyright notice in license.txt */
7 /*DrvContPump local continuum pumping rate radiative transfer for all lines */
8 /*con_pump_op routine used to get continuum pumping of lines
9  * used in DrvContPump in call to qg32 */
10 
11 #include "cddefines.h"
12 #include "rt.h"
13 #include "rt_escprob.h"
14 #include "transition.h"
15 #include "thirdparty.h"
16 #include "integrate.h"
17 
18 #include "cosmology.h"
19 
20 /*conpmp local continuum pumping rate radiative transfer for all lines */
21 STATIC double conpmp(double tau, double damp);
22 STATIC double conpmp_romb( double tau, double damp);
23 STATIC double growth_romb( double tau, double damp);
24 STATIC inline double fitted( double t );
25 
26 namespace {
27 
28 /*con_pump_op routine used to get continuum pumping of lines
29  * used in DrvContPump in call to qg32 */
30 class my_Integrand_con_pump_op
31 {
32 public:
33  /* variable used for inward optical depth for pumping */
34  realnum PumpTau;
35  /* damping constant used for pumping */
36  realnum damp;
37  my_Integrand_con_pump_op( realnum tau, realnum damp) :
38  PumpTau(tau), damp(damp)
39  {}
40 private:
41  // Do not implement
42  my_Integrand_con_pump_op( void );
43 public:
44  double operator() (double x) const
45  {
46  realnum v, rx = realnum(x);
47  VoigtH(damp,&rx,&v,1);
48  double opfun_v = sexp(PumpTau*v)*v;
49  return opfun_v;
50  }
51 };
52 
53 class con_pump_op_conv
54 {
55 public:
56  /* variable used for inward optical depth for pumping */
57  realnum PumpTau;
58  /* damping constant used for pumping */
59  realnum damp;
60  con_pump_op_conv( realnum tau, realnum damp) :
61  PumpTau(tau), damp(damp)
62  {}
63 private:
64  // Do not implement
65  con_pump_op_conv( void );
66 public:
67  double operator() (double y) const
68  {
69  // Substitute y = 1/(1+x) to map integral onto a finite domain
70  double x = 1./y - 1.;
71  realnum v, rx = realnum(x);
72  if (damp >= 0.0)
73  VoigtU(damp,&rx,&v,1);
74  else
75  v = 1.0/(PI*(1.0+x*x));
76  double opfun_v = sexp(PumpTau*v)*v;
77  return opfun_v/(y*y);
78  }
79 };
80 
81 class profile_pow
82 {
83 public:
84  /* power of profile function */
85  int npow;
86  /* damping constant used for pumping */
87  realnum damp;
88  profile_pow( int npow, realnum damp) :
89  npow(npow), damp(damp)
90  {}
91 private:
92  // Do not implement
93  profile_pow( void );
94 public:
95  double operator() (double y) const
96  {
97  // Substitute y = 1/(1+x) to map integral onto a finite domain
98  double x = 1./y - 1.;
99  realnum v, rx = realnum(x);
100  if (damp >= 0.)
101  {
102  VoigtU(damp,&rx,&v,1);
103  }
104  else
105  {
106  // Normalized Lorentz for negative a
107  v = 1./(PI*(1.0+x*x));
108  }
109  double opfun_v = powi(v,npow);
110  return opfun_v/(y*y);
111  }
112 };
113 
114 class growth
115 {
116 public:
117  /* variable used for inward optical depth for pumping */
118  realnum PumpTau;
119  /* damping constant used for pumping */
120  realnum damp;
121  growth( realnum tau, realnum damp) :
122  PumpTau(tau), damp(damp)
123  {}
124 private:
125  // Do not implement
126  growth( void );
127 public:
128  double operator() (double y) const
129  {
130  // Substitute y = 1/(1+x) to map integral onto a finite domain
131  double x = 1./y - 1.;
132  realnum v, rx = realnum(x);
133  if (damp >= 0.0)
134  VoigtU(damp,&rx,&v,1);
135  else
136  v = 1.0/(PI*(1.0+x*x));
137  double opfun_v = PumpTau*v;
138  if (opfun_v > 1e-6)
139  opfun_v = 1-sexp(opfun_v);
140  else
141  opfun_v *= 1-0.5*opfun_v;
142  return opfun_v/(y*y);
143  }
144 };
145 
146 }
147 
148 static double shieldFederman(double tau, double damp, bool lgBug);
149 
150 static double shieldRodgers(double tau, double damp);
151 static double growthRodgers(double tau, double damp);
152 
153 /*rt_continuum_shield_fcn computing continuum shielding due to single line */
155 {
156  DEBUG_ENTRY( "rt_continuum_shield_fcn_point()" );
157  double value = -1.;
158 
159  ASSERT( t.Emis().damp() > 0. );
160  if( cosmology.lgDo )
161  {
162  if( tau < 1e-5 )
163  value = (1. - tau/2.);
164  else
165  value = (1. - dsexp( tau ) )/ tau;
166  }
168  {
169  /* set continuum shielding pesc - shieding based on escape probability */
170  if( t.Emis().iRedisFun() == ipPRD )
171  {
172  value = esc_PRD_1side(tau,t.Emis().damp());
173  }
174  else if( t.Emis().iRedisFun() == ipCRD )
175  {
176  value = esca0k2(tau);
177  }
178  else if( t.Emis().iRedisFun() == ipCRDW )
179  {
180  value = esc_CRDwing_1side(tau,t.Emis().damp());
181  }
182  else if( t.Emis().iRedisFun() == ipLY_A )
183  {
184  value = esc_PRD_1side(tau,t.Emis().damp());
185  }
186  else
187  TotalInsanity();
188  }
190  {
191  // Last arg is whether to use buggy terms compatible
192  // with the implementation at 13.02 and previously, and
193  // as published by Federman et al. These have
194  // been shown to be incorrect, but may not be fixed in
195  // other codes. false means to use the corrected version.
196  value = shieldFederman(tau,t.Emis().damp(),false);
197  }
199  {
200  value = shieldFederman(tau,t.Emis().damp(),true);
201  }
203  {
204  /* set continuum shielding ferland */
205  value = conpmp( tau , t.Emis().damp() );
206  }
208  {
209  value = shieldRodgers(tau*SQRTPI,t.Emis().damp());
210  }
212  {
213  value = conpmp_romb(tau*SQRTPI,t.Emis().damp());
214  value = MIN2(1.,value);
215  }
216  else if( rt.nLineContShield == 0 )
217  {
218  /* set continuum shielding none */
219  value = 1.;
220  }
221  else
222  {
223  TotalInsanity();
224  }
225 
226  /* the returned pump shield function must be greater than zero,
227  * and less than 1 if a maser did not occur */
228  ASSERT( value>=0 && (value<=1.||tau<0.) );
229  return value;
230 }
231 
232 // Function to collect a range of options to handle finite optical
233 // depth effects for continuum shielding. s1 is the shielding factor
234 // evaluated at the front face, s2 is the shielding factor evaluated at
235 // the rear face.
236 STATIC double avg_shield(double s1, double s2, double dTau, double tau)
237 {
238  enum options { AVG_OLD, AVG_DEPTH, AVG_NO, AVG_BACK, AVG_ARITH,
239  AVG_GEOM, AVG_COUNT };
240  const enum options opt = AVG_DEPTH;
241  if (opt == AVG_OLD)
242  {
243  return s1*log(1.+dTau)/dTau;
244  }
245  else if (opt == AVG_DEPTH)
246  {
247  // Average (1+tau)/(1+tau+dTau) over zone...
248  double dTauRel = dTau/(1.+tau);
249  ASSERT( 1.+dTauRel >= 0. );
250  return s1*log(1.+dTauRel)/dTauRel;
251  }
252  else if (opt == AVG_NO)
253  {
254  return s1; // no shielding
255  }
256  else if (opt == AVG_BACK)
257  {
258  return s2; // implicit
259  }
260  else if (opt == AVG_ARITH)
261  {
262  return 0.5*(s1+s2); // Arithmetic average
263  }
264  else if (opt == AVG_GEOM)
265  {
266  return sqrt(s1*s2); // Geometric average
267  }
268  else if (opt == AVG_COUNT)
269  {
270  return (s1-s2)/dTau; // photon counting estimate [?]
271  }
272  TotalInsanity(); // invalid option
273 }
274 
275 STATIC double shieldFederman(double tau, double damp, bool lgBug)
276 {
277  DEBUG_ENTRY( "shieldFederman()" );
278  /* set continuum shielding Federman - this is the default */
279  double core, wings, value;
280 
281  /* these expressions implement the appendix of
282  * >>refer line shielding Federman, S.R., Glassgold, A.E., &
283  * >>refercon Kwan, J. 1979, ApJ, 227, 466 */
284  /* doppler core - equation A8 */
285  if( tau < 2. )
286  {
287  core = sexp( tau * 0.66666 );
288  }
289  else if( tau < 10. )
290  {
291  core = 0.638 * pow(tau,(realnum)-1.25f );
292  }
293  else if( tau < 100. )
294  {
295  core = 0.505 * pow(tau,(realnum)-1.15f );
296  }
297  else
298  {
299  core = 0.344 * pow(tau,(realnum)-1.0667f );
300  }
301 
302  /* do we add damping wings? */
303  wings = 0.;
304  if( damp>0. )
305  {
306  /* equation A6 */
307  double t1 = 3.02*pow(damp*1e3,-0.064 );
308  double tauwing = lgBug ? tau : tau*SQRTPI;
309  double u1 = sqrt(MAX2(tauwing,0.)*damp )/SDIV(t1);
310  wings = damp/SDIV(t1)/sqrt( 0.78540 + POW2(u1) );
311  /* add very large optical depth tail to converge this with respect
312  * to escape probabilities - if this function falls off more slowly
313  * than escape probability then upper level will become overpopulated.
314  * original paper was not intended for this regime */
315  if( lgBug && tau>1e7 )
316  wings *= pow( tau/1e7,-1.1 );
317  }
318  value = core + wings;
319  /* some x-ray lines have vastly large damping constants, greater than 1.
320  * in these cases the previous wings value does not work - approximation
321  * is for small damping constant - do not let pump efficiency exceed unity
322  * in this case */
323  if( tau>=0. )
324  value = MIN2(1., value );
325  return value;
326 }
327 
329  bool lgShield_this_zone, double dTau )
330 {
331  DEBUG_ENTRY( "rt_continuum_shield_fcn()" );
332 
333  double tau = t.Emis().TauCon();
334  double value = RT_continuum_shield_fcn_point(t,tau);
335 
336  if( lgShield_this_zone && dTau > 1e-3 )
337  {
338  if (0 && t.ipCont() == 3627)
339  fprintf(ioQQQ,"?shield %ld %15.8g %15.8g %15.8g %15.8g %s\n",
340  nzone,tau,dTau,1./(1. + dTau ),
341  RT_continuum_shield_fcn_point(t,tau+dTau)/value,
342  chLineLbl(t).c_str());
343  /* correction for line optical depth across zone */
344  value = avg_shield(value,RT_continuum_shield_fcn_point(t,tau+dTau),
345  dTau,tau);
346  }
347 
348  return value;
349 }
350 
351 /*conpmp local continuum pumping rate radiative transfer for all lines */
352 STATIC double conpmp_qg32( double tau, double damp)
353 {
354  DEBUG_ENTRY( "conpmp_qg32()" );
355 
356  my_Integrand_con_pump_op func(tau, damp);
358 
359  const double BREAK = 3.;
360 
361  double yinc1 = opfun.sum( 0., BREAK );
362  double yinc2 = opfun.sum( BREAK, 100. );
363 
364  double a0 = 0.5*SQRTPI;
365  return (yinc1 + yinc2)/a0;
366 }
367 
368 /*conpmp local continuum pumping rate radiative transfer for all lines */
369 STATIC double conpmp_romb( double tau, double damp)
370 {
371  DEBUG_ENTRY( "conpmp_romb()" );
372 
373  double tauint = tau;
374 
375  con_pump_op_conv func(tauint, damp);
376 
377  // Ignore underflowing values at small offsets x, i.e. y~=1
378  double top=1.0;
379  for (;;)
380  {
381  if (func(0.5*top) > 0.0)
382  break;
383  top *= 0.5;
384  }
385 
386  class integrate::Romberg<integrate::Midpoint<con_pump_op_conv> >
387  intl(integrate::Midpoint<con_pump_op_conv>(func,0.0,top));
388  intl.update(1e-10);
389  // 2.0 because only integrating RHS
390  return 2.0*intl.sum();
391 }
392 
393 //*conpmp local continuum pumping rate radiative transfer for all lines */
394 STATIC double growth_romb( double tau, double damp)
395 {
396  DEBUG_ENTRY( "growth_romb()" );
397 
398  double tauint = tau;
399 
400  growth func(tauint, damp);
401 
402  class integrate::Romberg<integrate::Midpoint<growth> >
403  intl(integrate::Midpoint<growth>(func,0.0,1.0));
404  intl.update(1e-10);
405  // 2.0 because only integrating RHS
406  return 2.0*intl.sum();
407 }
408 
409 /*conpmp local continuum pumping rate radiative transfer for all lines */
410 STATIC double conpmp( double tau, double damp)
411 {
412  DEBUG_ENTRY( "conpmp()" );
413  double conpmp_v;
414 
415  /* tau required is optical depth in center of next zone */
416  /* compute pumping probability */
417  if( tau <= 10+2.5*damp)
418  {
419  double tausc;
420  if (damp < 1e-2)
421  {
422  tausc = tau*(1.-1.5*damp);
423  }
424  else
425  {
426  tausc = tau*(1+damp)/(1.+2.5*damp*(1.+damp));
427  }
428  conpmp_v = fitted(tausc);
429  }
430  else if( tau > 1e6 )
431  {
432  /* this far in winds line opacity well below electron scattering
433  * so ignore for this problem */
434  conpmp_v = 0.;
435  }
436  else
437  {
438  conpmp_v = conpmp_qg32(tau, damp);
439  }
440  return conpmp_v;
441 }
442 
443 /* fit to results for tau less than 10 */
444 inline double fitted(double t)
445 {
446  return (0.98925439 + 0.084594094*t)/(1. + t*(0.64794212 + t*0.44743976));
447 }
448 
449 double DrvContPump(double tau, double damp)
450 {
451  DEBUG_ENTRY( "DrvContPump()" );
452  return conpmp(tau, damp);
453 }
454 
455 void DrivePump(double damp)
456 {
457  double taufac = SQRTPI;
458  double dampa = damp > 0.0 ? damp : 0.0;
459  fprintf(ioQQQ,"#tau damp QG32 conpmp Federman Romb\n");
460  for (long i=0; i<49; ++i)
461  {
462  double tau = 1e-4*exp10(i/4.);
463  fprintf(ioQQQ,"%13.8e %13.8e %13.8e %13.8e %13.8e %13.8e %13.8e\n",
464  tau, damp, conpmp_qg32(tau/taufac,dampa),
465  conpmp(tau/taufac,dampa),
466  shieldFederman(tau/taufac,dampa,true),conpmp_romb(tau,damp),
467  shieldRodgers(tau,damp)
468  );
469  }
470  fprintf(ioQQQ,"\n");
471 
472  if (0)
473  {
474  int NSAMP=10;
475  FILE* ioVALID = open_data("growth.dat","w");
476  fprintf(ioVALID,"#tau damp growth_romb growthRodgers\n");
477  for (long j=0; j<1+7*NSAMP; ++j)
478  {
479  double damp1 = 1e-5*exp10(j/double(NSAMP));
480  for (long i=0; i<1+12*NSAMP; ++i)
481  {
482  double tau = 1e-4*exp10(i/double(NSAMP));
483  double gromb = growth_romb(tau,damp1),
484  grodg = growthRodgers(tau,damp1),
485  err = safe_div(grodg,gromb,1.)-1.;
486  fprintf(ioVALID,"%13.8e %13.8e %13.8e %13.8e %13.8e\n",
487  tau, damp1,gromb, grodg, err );
488  }
489  }
490  fclose(ioVALID);
491  }
492 
493  if (0)
494  {
495  int NSAMP=10;
496  FILE* ioVALID = open_data("shield.dat","w");
497  fprintf(ioVALID,"#tau damp shield_romb shieldRodgers\n");
498  for (long j=0; j<1+7*NSAMP; ++j)
499  {
500  double damp1 = 1e-5*exp10(j/double(NSAMP));
501  for (long i=0; i<1+12*NSAMP; ++i)
502  {
503  double tau = 1e-4*exp10(i/double(NSAMP));
504  double gromb = conpmp_romb(tau,damp1),
505  grodg = shieldRodgers(tau,damp1), // shieldFederman(tau/taufac,damp1,false),
506  err = safe_div(grodg,gromb,1.)-1.;
507  fprintf(ioVALID,"%13.8e %13.8e %13.8e %13.8e %13.8e\n",
508  tau, damp1,gromb, grodg, err );
509  }
510  }
511  fclose(ioVALID);
512  }
513 
514  if (damp >= 0.0)
515  {
516  double tau = 1.0;
517  if (damp != 0.0 && damp < 0.5)
518  {
519  tau = 1.0/damp;
520  for (;;)
521  {
522  double tau0 = tau;
523  tau = 1.0/(damp*log(tau));
524  if (fabs(tau-tau0) < 1e-4*tau)
525  {
526  break;
527  }
528  }
529  }
530  for (long i=0; i<49; ++i)
531  {
532  double damp1 = 1e-6*damp*exp10(0.25*i);
533  fprintf(ioQQQ,"%13.8e %13.8e %13.8e %13.8e %13.8e %13.8e\n",
534  tau, damp1, conpmp_qg32(tau,damp1), conpmp(tau,damp1),
535  shieldFederman(tau,damp1,true),conpmp_romb(tau,damp1)
536  );
537  }
538  fprintf(ioQQQ,"\n");
539  }
540 
541  fprintf(ioQQQ,"Integrals of phi^n\n");
542  fprintf(ioQQQ,"%13.8e", damp);
543  for (int n=1; n<10; ++n)
544  {
545  class integrate::Romberg<integrate::Midpoint<profile_pow> >
546  intg(integrate::Midpoint<profile_pow>(profile_pow(n,damp),0.0,1.0));
547  intg.update(1e-10);
548  fprintf(ioQQQ," %13.8e",2.0*intg.sum());
549  }
550  fprintf(ioQQQ,"\n");
551 }
552 
553 // Rodgers & Williams JQSRT 14, 319 (1974)
554 // tau phi(x) = -k(v) m
555 // a -> y = alpha_L / alphaD
556 // tau -> S m / alphaD
557 // x -> (nu-nu0) / alphaD
558 // W = int 1-exp(-k m) dnu = alphaD int 1-exp(-phi tau) dx
559 // So can divide all W's by alphaD to normalize
560 // WM = [WL^2 + WD^2 - (WL WD/WW)^2]^{1/2}
561 // stands.
562 // WL = 2 pi a L ( tau / 2 pi a)
563 // WD = D ( tau / sqrtpi )
564 // WW = S m / alphaD = tau
565 
566 static void shieldRodgersLorentz(double tau, double damp,
567  double& w, double &dw)
568 {
569  // z = S m / 2 pi alphaL = tau / 2 sqrt(pi) y
570  double z = tau/(2 * PI * damp ); // Check scaling
571  if (z < 1.98)
572  {
573  static const double a[] = {9.99997697674e-1,
574  -4.99882252233e-1,
575  2.49005331874e-1,
576  -1.00956020660e-1,
577  3.13580341312e-2,
578  -6.50144170817e-3,
579  6.48325819427e-4};
580  w = z*(a[0]+z*(
581  a[1]+z*(
582  a[2]+z*(
583  a[3]+z*(
584  a[4]+z*(
585  a[5]+z*a[6]))))));
586  dw = a[0]+z*(
587  2*a[1]+z*(
588  3*a[2]+z*(
589  4*a[3]+z*(
590  5*a[4]+z*(
591  6*a[5]+z*7*a[6])))));
592  }
593  else
594  {
595  // Note that this is an expansion in z^-1, not z, a typo
596  // in the original paper
597  static const double b[] = {7.97885095473e-1,
598  -9.97555137671e-2,
599  -1.97623661059e-2,
600  1.05865022723e-2,
601  -1.32467496350e-1,
602  1.47947878333e-1};
603  double sz = sqrt(z),rz = 1./z;
604  w = sz*(b[0]+rz*(b[1]+rz*(b[2]+rz*(b[3]+rz*(b[4]+rz*b[5])))));
605  dw = (b[0]+rz*(
606  -1*b[1]+rz*(
607  -3*b[2]+rz*(
608  -5*b[3]+rz*(
609  -7*b[4]-rz*9*b[5])))))/
610  (2.*sz);
611  }
612  w *= 2.*PI*damp;
613  //dw *= damp;
614 }
615 
616 static void shieldRodgersDoppler(double tau,
617  double& w, double &dw)
618 {
619  // z = S m / pi^1/2 alphaD
620  double z = tau/SQRTPI; // Check scaling
621  if (z < 5)
622  {
623  static const double c[] = {9.99998291698e-1,
624  -3.53508187098e-1,
625  9.60267807976e-2,
626  -2.04969011013e-2,
627  3.43927368627e-3,
628  -4.27593051557e-4,
629  3.42209457833e-5,
630  -1.28380804108e-6};
631  w = z*SQRTPI*(c[0]+z*(
632  c[1]+z*(
633  c[2]+z*(
634  c[3]+z*(
635  c[4]+z*(
636  c[5]+z*(
637  c[6]+z*c[7])))))));
638  dw = (c[0]+z*(
639  2*c[1]+z*(
640  3*c[2]+z*(
641  4*c[3]+z*(
642  5*c[4]+z*(
643  6*c[5]+z*(
644  7*c[6]+z*8*c[7])))))));
645  }
646  else
647  {
648  static const double d[] = {1.99999898289e0,
649  5.77491987800e-1,
650  -5.05367549898e-1,
651  8.21896973657e-1,
652  -2.52226724530e0,
653  6.10070274810e0,
654  -8.51001627836e0,
655  4.65351167650e0};
656  double lz = log(z), slz=sqrt(lz), rlz=1./lz;
657  w = slz*(d[0]+rlz*(
658  d[1]+rlz*(
659  d[2]+rlz*(
660  d[3]+rlz*(
661  d[4]+rlz*(
662  d[5]+rlz*(
663  d[6]+rlz*d[7]))))))); //SQRTPI;
664  dw = (d[0]+rlz*(
665  -1*d[1]+rlz*(
666  -3*d[2]+rlz*(
667  -5*d[3]+rlz*(
668  -7*d[4]+rlz*(
669  -9*d[5]+rlz*(
670  -11*d[6]-rlz*13*d[7])))))))/
671  (2.*slz*z*SQRTPI);
672  }
673 }
674 
675 // Rodgers & Williams curve of growth
676 static double growthRodgers(double tau, double damp)
677 {
678  if (damp < 0)
679  {
680  double wl, dwl;
681  shieldRodgersLorentz(tau,1.0,wl,dwl);
682  return wl;
683  }
684 
685  double wd, dwd;
686  shieldRodgersDoppler(tau,wd,dwd);
687  if (damp == 0.0)
688  return wd;
689 
690  double wl, dwl;
691  shieldRodgersLorentz(tau,damp,wl,dwl);
692  double wls = wl*wl, wds = wd*wd, rtaus=1./(tau*tau);
693  double wv = sqrt(wls+wds-wls*wds*rtaus);
694  return wv;
695 }
696 
697 // Rodgers & Williams derived shielding function
698 static double shieldRodgers(double tau, double damp)
699 {
700  if (damp < 0)
701  {
702  double wl, dwl;
703  shieldRodgersLorentz(tau,1.0,wl,dwl);
704  return dwl;
705  }
706 
707  double wd, dwd;
708  shieldRodgersDoppler(tau,wd,dwd);
709  if (damp == 0.0)
710  return dwd;
711 
712  double wl, dwl;
713  shieldRodgersLorentz(tau,damp,wl,dwl);
714  double wls = wl*wl, wds = wd*wd, rtaus=1./(tau*tau);
715  double wv = sqrt(wls+wds-wls*wds*rtaus);
716  double dwv = (wl*dwl*(1.0-wds*rtaus)+
717  wd*dwd*(1.0-wls*rtaus)+
718  wls*wds*rtaus/tau)/wv;
719  return dwv;
720 }
static void shieldRodgersLorentz(double tau, double damp, double &w, double &dw)
FILE * open_data(const char *fname, const char *mode, access_scheme scheme)
Definition: cpu.cpp:765
int & iRedisFun() const
Definition: emission.h:450
double exp10(double x)
Definition: cddefines.h:1368
string chLineLbl(const TransitionProxy &t)
Definition: transition.h:599
NORETURN void TotalInsanity(void)
Definition: service.cpp:971
double DrvContPump(double tau, double damp)
STATIC double conpmp_qg32(double tau, double damp)
void DrivePump(double tau)
void update(double eps)
Definition: integrate.h:203
STATIC double RT_continuum_shield_fcn_point(const TransitionProxy &t, double tau)
STATIC double conpmp(double tau, double damp)
sys_float sexp(sys_float x)
Definition: service.cpp:999
static double growthRodgers(double tau, double damp)
FILE * ioQQQ
Definition: cddefines.cpp:7
long int nzone
Definition: cddefines.cpp:14
double sum() const
Definition: integrate.h:188
#define MIN2(a, b)
Definition: cddefines.h:803
bool lgDo
Definition: cosmology.h:44
double dsexp(double x)
Definition: service.cpp:1038
static double shieldFederman(double tau, double damp, bool lgBug)
STATIC double conpmp_romb(double tau, double damp)
#define POW2
Definition: cddefines.h:979
STATIC double fitted(double t)
#define STATIC
Definition: cddefines.h:118
EmissionList::reference Emis() const
Definition: transition.h:447
STATIC double growth_romb(double tau, double damp)
const int ipCRD
Definition: cddefines.h:342
long & ipCont() const
Definition: transition.h:489
float realnum
Definition: cddefines.h:124
const int ipLY_A
Definition: cddefines.h:346
static double shieldRodgers(double tau, double damp)
double powi(double, long int)
Definition: service.cpp:594
sys_float safe_div(sys_float x, sys_float y, sys_float res_0by0)
Definition: cddefines.h:1011
double sum(double min, double max) const
Definition: integrate.h:44
double RT_continuum_shield_fcn(const TransitionProxy &t, bool lgShieldThisZone, double dTau)
double esc_PRD_1side(double tau, double a)
Definition: rt_escprob.cpp:116
static void shieldRodgersDoppler(double tau, double &w, double &dw)
#define ASSERT(exp)
Definition: cddefines.h:613
realnum & TauCon() const
Definition: emission.h:510
void VoigtU(realnum a, const realnum v[], realnum y[], int n)
Definition: thirdparty.h:431
int nLineContShield
Definition: rt.h:190
#define DEBUG_ENTRY(funcname)
Definition: cddefines.h:723
t_cosmology cosmology
Definition: cosmology.cpp:8
#define MAX2(a, b)
Definition: cddefines.h:824
int fprintf(const Output &stream, const char *format,...)
Definition: service.cpp:1121
realnum & damp() const
Definition: emission.h:620
const int ipPRD
Definition: cddefines.h:340
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
static double a0[83]
const int ipCRDW
Definition: cddefines.h:344
void VoigtH(realnum a, const realnum v[], realnum y[], int n)
Definition: thirdparty.h:418
STATIC double avg_shield(double s1, double s2, double dTau, double tau)
t_rt rt
Definition: rt.cpp:5