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