Cantera  2.1.2
1 #include "cantera/base/ct_defs.h"
2 #include <math.h>
3 #include <iostream>
4
5 #ifdef USE_BOOST_MATH
6 #include <boost/math/special_functions/erf.hpp>
7 using boost::math::erf;
8 #endif
9
11
12 namespace Cantera
13 {
14
15 LorentzianProfile::LorentzianProfile(doublereal gamma)
16 {
17  m_hwhm = gamma;
18  m_hwhm2 = m_hwhm*m_hwhm;
19 }
20
21 /**
22  * The Lorentzian profile for collision-broadened lines.
23  *
24  *\f[
25  * \frac{1}{\pi} \frac{\gamma}{ (\Delta\nu)^2 + \gamma^2}
26  *\f]
27  * Units: 1/wavenumber (or cm).
28  */
29 doublereal LorentzianProfile::profile(doublereal deltaFreq)
30 {
31  return (1.0/Cantera::Pi) *m_hwhm/(deltaFreq*deltaFreq + m_hwhm2);
32 }
33
34 /**
35  *
36  * The cumulative profile, given by
37  * \f[
38  * \frac{1}{\pi} \tan^{-1}\left(\frac{\Delta\nu}{gamma}\right) + 0.5
39  * \f]
40  */
41 doublereal LorentzianProfile::cumulative(doublereal deltaFreq)
42 {
43  return (1.0/Pi) * atan(deltaFreq/m_hwhm) + 0.5;
44 }
45
46 doublereal LorentzianProfile::width()
47 {
48  return 2.0*m_hwhm;
49 }
50
52 {
53  m_sigma = sigma;
54  m_sigma2 = m_sigma*m_sigma;
55 }
56
57 doublereal GaussianProfile::profile(doublereal deltaFreq)
58 {
59  //cout << "entered Gaussian::profile" << endl;
60  //cout << "deltaFreq = " << deltaFreq << endl;
61  //cout << "m_sigma = " << m_sigma << endl;
62  return 1.0/(m_sigma * Cantera::SqrtTwo *Cantera::SqrtPi) *
63  exp(-deltaFreq*deltaFreq/(2.0*m_sigma2));
64 }
65
66 doublereal GaussianProfile::cumulative(doublereal deltaFreq)
67 {
68  return 0.5*(1.0 + erf(deltaFreq/(m_sigma*SqrtTwo)));
69 }
70
71 doublereal GaussianProfile::width()
72 {
73  return 2.0*m_sigma*sqrt(log(4.0));
74 }
75
76
77
78 /**
79  * @param sigma The standard deviation of the Gaussian
80  * @param gamma The half-width of the Lorentzian.
81  */
82 Voigt::Voigt(doublereal sigma, doublereal gamma)
83 {
84  m_sigma = sigma;
85  m_sigma2 = m_sigma*m_sigma;
86  m_gamma_lor = gamma;
87  m_sigsqrt2 = SqrtTwo*m_sigma;
88  m_gamma = gamma/m_sigsqrt2;
89  m_eps = 1.0e-20;
90 }
91
92 void Voigt::testv()
93 {
94  m_gamma = 1.0e1;
95  std::cout << F(1.0) << std::endl;
96  m_gamma = 0.5;
97  std::cout << F(1.0) << std::endl;
98  m_gamma = 0.0001;
99  std::cout << F(10.0) << std::endl;
100 }
101 /**
102  * This method evaluates the function
103  * \f[
104  * F(x, y) = \frac{y}{\pi}\int_{-\infty}^{+\infty} \frac{e^{-z^2}}
105  * {(x - z)^2 + y^2} dz
106  * \f]
107  * The algorithm used to cmpute this function is described in the
108  * reference below. @see F. G. Lether and P. R. Wenston, "The
109  * numerical computation of the %Voigt function by a corrected
110  * midpoint quadrature rule for \f$(-\infty, \infty) \f$. Journal
111  * of Computational and Applied Mathematics}, 34 (1):75--92, 1991.
112  */
113 doublereal Voigt::F(doublereal x)
114 {
115
116  if (x < 0.0) {
117  x = -x;
118  }
119  double y = m_gamma;
120
121  double c3 = log(Pi*m_eps/2.0);
122  double tau = sqrt(-log(y) - c3);
123  double b = (tau + x)/y;
124  double t = b*y;
125  double f1, f2, f3;
126  const double c0 = 2.0/(Pi*exp(0.0));
127  const double c1 = 1.0/SqrtTwo;
128  const double c2 = 2.0/SqrtPi;
129
130  if (y > c0/m_eps) {
131  return 0.0;
132  }
133  double f0, ef0;
134  while (1 > 0) {
135  f0 = Pi*Pi/(t*t);
136  ef0 = exp(-f0);
137  f1 = c2*y*ef0;
138  f2 = fabs(y*y - Pi*Pi/(t*t));
139  f3 = 1.0 - ef0*ef0;
140  t *= c1;
141  if (f1/(f2*f3) < 0.5*m_eps) {
142  break;
143  }
144  }
145  double h = t/y;
146  int N = int(0.5 + b/h);
147  double S = 0.0;
148  double u = h/2;
149  for (int i = 0; i < N; i++) {
150  S += (1.0 + exp(-4.0*x*y*u))*exp(-pow(y*u-x,2))/(u*u+1.0);
151  u += h;
152  }
153  double Q = h*S/Pi;
154  double C = 0.0;
155  if (y*y < Pi/h) {
156  C = 2.0*exp(y*y - x*x)*cos(2*x*y)/(1.0 + exp(2*Pi/h));
157  } else {
158  return 0.0;
159  }
160  return Q + C;
161 }
162
163 /**
164  * Voigt profile.
165  *
166  * Not sure that constant is right.
167  */
168 doublereal Voigt::profile(doublereal deltaFreq)
169 {
170  const double ff = 1.0/(m_sigsqrt2*SqrtPi);
171  return ff*F(deltaFreq/m_sigsqrt2);
172 }
173
174 }
Voigt(doublereal sigma, doublereal gamma)
Constructor.
virtual doublereal profile(doublereal deltaFreq)
Voigt profile.
doublereal F(doublereal x)
This method evaluates the function The algorithm used to cmpute this function is described in the re...
virtual doublereal cumulative(doublereal deltaFreq)
The cumulative profile, given by .
This file contains definitions of terms that are used in internal routines and are unlikely to need m...
const doublereal Pi
Pi.
Definition: ct_defs.h:51
virtual doublereal profile(doublereal deltaFreq)
The line shape profile as a function of distance from line center .
virtual doublereal cumulative(doublereal deltaFreq)
The cumulative profile, defined as .