cloudy  trunk
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
TestThirdparty.cpp
Go to the documentation of this file.
1 /* This file is part of Cloudy and is copyright (C)1978-2013 by Gary J. Ferland and
2  * others. For conditions of distribution and use see copyright notice in license.txt */
3 #include <UnitTest++.h>
4 #include "cddefines.h"
5 #include "thirdparty.h"
6 
7 namespace {
8 
9  // constants below were calculated with xmaxima 5.17.0 with fpprec set to 30
10 
11  TEST(TestFactorial)
12  {
13  CHECK( fp_equal( factorial(0), 1. ) );
14  CHECK( fp_equal( factorial(1), 1. ) );
15  CHECK( fp_equal( factorial(2), 2. ) );
16  CHECK( fp_equal( factorial(10), 3628800. ) );
17  CHECK( fp_equal( factorial(170), 7.257415615307998967396728211129316e306 ) );
18  }
19 
20  TEST(TestLogFactorial)
21  {
22  CHECK( fp_equal( lfactorial(0), 0. ) );
23  CHECK( fp_equal( lfactorial(1), 0. ) );
24  CHECK( fp_equal( lfactorial(2), 3.01029995663981195213738894725e-1, 10 ) );
25  CHECK( fp_equal( lfactorial(10), 6.55976303287679375117476123996e0, 10 ) );
26  CHECK( fp_equal( lfactorial(170), 3.06860781994828320192380273111e2, 10 ) );
27  CHECK( fp_equal( lfactorial(1700), 4.75547688279135071178908638597e3, 10 ) );
28  }
29 
30  TEST(TestCDGamma)
31  {
32  complex<double> y = cdgamma( complex<double>(1.,0.) );
33  CHECK( fp_equal( y.real(), 1., 10 ) && fp_equal( y.imag(), 0. ) );
34  y = cdgamma( complex<double>(2.,0.) );
35  CHECK( fp_equal( y.real(), 1., 10 ) && fp_equal( y.imag(), 0. ) );
36  y = cdgamma( complex<double>(11.,0.) );
37  CHECK( fp_equal( y.real(), 3628800., 50 ) && fp_equal( y.imag(), 0. ) );
38  y = cdgamma( complex<double>(-0.5,0.) );
39  CHECK( fp_equal( y.real(), -3.544907701811032054596334966682277e0, 10 ) &&
40  fp_equal( y.imag(), 0. ) );
41  y = cdgamma( complex<double>(0.,1.) );
42  CHECK( fp_equal( y.real(), -1.549498283018106851249551304838863e-1, 10 ) &&
43  fp_equal( y.imag(), -4.980156681183560427136911174621973e-1, 10 ) );
44  y = cdgamma( complex<double>(-1.,-2.) );
45  CHECK( fp_equal( y.real(), -3.23612885501927256868232032760969e-2, 30 ) &&
46  fp_equal( y.imag(), -1.122942423463261735043406872030743e-2, 30 ) );
47  }
48 
49  // constants below were taken from Abramowitz & Stegun, Handbook of Mathematical Functions
50 
51  TEST(TestBesselJ0)
52  {
53  CHECK( fp_equal( bessel_j0(0.), 1., 10 ) );
54  CHECK( fp_equal_tol( bessel_j0(1.), 0.765197686557967, 1.e-15 ) );
55  CHECK( fp_equal_tol( bessel_j0(2.9), -0.224311545791968, 1.e-15 ) );
56  CHECK( fp_equal_tol( bessel_j0(4.8), -0.240425327291183, 1.e-15 ) );
57  CHECK( fp_equal_tol( bessel_j0(8.3), 0.096006100895010, 1.e-15 ) );
58  CHECK( fp_equal_tol( bessel_j0(17.4), -0.118955856336348, 1.e-15 ) );
59  }
60 
61  TEST(TestBesselJ1)
62  {
63  CHECK( fp_equal( bessel_j1(0.), 0., 10 ) );
64  CHECK( fp_equal( bessel_j1(1.e-30), 5.e-31, 10 ) );
65  CHECK( fp_equal_tol( bessel_j1(1.e-15), 5.e-16, 1.e-27 ) );
66  CHECK( fp_equal_tol( bessel_j1(1.), 0.4400505857, 1.e-10 ) );
67  CHECK( fp_equal_tol( bessel_j1(2.9), 0.3754274818, 1.e-10 ) );
68  CHECK( fp_equal_tol( bessel_j1(4.8), -0.2984998581, 1.e-10 ) );
69  CHECK( fp_equal_tol( bessel_j1(8.3), 0.2657393020, 1.e-10 ) );
70  CHECK( fp_equal_tol( bessel_j1(17.4), -0.1532161760, 1.e-10 ) );
71  }
72 
73  TEST(TestBesselJn)
74  {
75  CHECK( fp_equal( bessel_jn(2,0.), 0., 10 ) );
76  CHECK( fp_equal( bessel_jn(2,1.e-30), 1.25e-61, 10 ) );
77  CHECK( fp_equal_tol( bessel_jn(2,2.e-15), 5.e-31, 1.e-42 ) );
78  CHECK( fp_equal_tol( bessel_jn(2,2.e-10), 5.e-21, 1.e-32 ) );
79  CHECK( fp_equal_tol( bessel_jn(2,0.099999999999), 0.0012489587, 1.e-10 ) );
80  CHECK( fp_equal_tol( bessel_jn(2,0.100000000001), 0.0012489587, 1.e-10 ) );
81  CHECK( fp_equal_tol( bessel_jn(2,1.), 0.1149034849, 1.e-10 ) );
82  CHECK( fp_equal_tol( bessel_jn(2,2.9), 0.4832270505, 1.e-10 ) );
83  CHECK( fp_equal_tol( bessel_jn(2,4.8), 0.1160503864, 1.e-10 ) );
84  CHECK( fp_equal_tol( bessel_jn(2,8.3), -0.0319725341, 1.e-10 ) );
85  CHECK( fp_equal_tol( bessel_jn(2,17.4), 0.1013448016, 1.e-10 ) );
86 
87  CHECK( fp_equal( bessel_jn(8,0.), 0., 10 ) );
88  CHECK( fp_equal( bessel_jn(8,2.e-20), 1.e-160/40320., 10 ) );
89  CHECK( fp_equal_tol( bessel_jn(8,2.e-15), 1.e-120/40320., 1.e-136 ) );
90  CHECK( fp_equal_tol( bessel_jn(8,1.), 9.4223e-8, 1.e-12 ) );
91  CHECK( fp_equal_tol( bessel_jn(8,2.8), 2.9367e-4, 1.e-8 ) );
92  CHECK( fp_equal_tol( bessel_jn(8,4.8), 1.4079e-2, 1.e-6 ) );
93  CHECK( fp_equal_tol( bessel_jn(8,8.2), 0.24257, 1.e-5 ) );
94 
95  CHECK( fp_equal_tol( bessel_jn(20,1.), 3.873503009e-25, 1.e-34 ) );
96  CHECK( fp_equal_tol( bessel_jn(50,1.), 2.906004948e-80, 1.e-89 ) );
97  CHECK( fp_equal_tol( bessel_jn(100,1.), 8.431828790e-189, 1.e-198 ) );
98  CHECK( fp_equal_tol( bessel_jn(100,100.), 9.636667330e-2, 1.e-11 ) );
99  }
100 
101  TEST(TestBesselY0)
102  {
103  CHECK( fp_equal_tol( bessel_y0(1.e-30), -44.0499402279, 1.e-10 ) );
104  CHECK( fp_equal_tol( bessel_y0(1.), 0.0882569642, 1.e-10 ) );
105  CHECK( fp_equal_tol( bessel_y0(2.9), 0.4079117692, 1.e-10 ) );
106  CHECK( fp_equal_tol( bessel_y0(4.8), -0.2723037945, 1.e-10 ) );
107  CHECK( fp_equal_tol( bessel_y0(8.3), 0.2595149638, 1.e-10 ) );
108  CHECK( fp_equal_tol( bessel_y0(17.4), -0.1497391883, 1.e-10 ) );
109  }
110 
111  TEST(TestBesselY1)
112  {
113  CHECK( fp_equal_tol( bessel_y1(1.e-30), -6.36619772368e29, 1.e18 ) );
114  CHECK( fp_equal_tol( bessel_y1(1.), -0.7812128213, 1.e-10 ) );
115  CHECK( fp_equal_tol( bessel_y1(2.9), 0.2959400546, 1.e-10 ) );
116  CHECK( fp_equal_tol( bessel_y1(4.8), 0.2135651673, 1.e-10 ) );
117  CHECK( fp_equal_tol( bessel_y1(8.3), -0.0805975035, 1.e-10 ) );
118  CHECK( fp_equal_tol( bessel_y1(17.4), 0.1147053859, 1.e-10 ) );
119  }
120 
121  TEST(TestBesselYn)
122  {
123  CHECK( fp_equal_tol( bessel_yn(2,1.e-30), -1.27323954474e60, 1.e49 ) );
124  CHECK( fp_equal_tol( bessel_yn(2,1.), -1.65068261, 1.e-8 ) );
125  CHECK( fp_equal_tol( bessel_yn(2,2.9), -0.20381518, 1.e-8 ) );
126  CHECK( fp_equal_tol( bessel_yn(2,4.8), 0.36128928, 1.e-8 ) );
127  CHECK( fp_equal_tol( bessel_yn(2,8.3), -0.27893605, 1.e-8 ) );
128  CHECK( fp_equal_tol( bessel_yn(2,17.4), 0.16292372, 1.e-8 ) );
129 
130  CHECK( fp_equal_tol( bessel_yn(8,2.e-20), -1.60428182637e163, 1.e152 ) );
131  CHECK( fp_equal_tol( bessel_yn(8,1.), -4.2567e5, 1.e1 ) );
132  CHECK( fp_equal_tol( bessel_yn(8,2.8), -1.4486e2, 1.e-2 ) );
133  CHECK( fp_equal_tol( bessel_yn(8,4.8), -3.5855, 1.e-4 ) );
134  CHECK( fp_equal_tol( bessel_yn(8,8.2), -0.35049, 1.e-5 ) );
135 
136  CHECK( fp_equal_tol( bessel_yn(20,1.), -4.113970315e22, 1.e13 ) );
137  CHECK( fp_equal_tol( bessel_yn(50,1.), -2.191142813e77, 1.e68 ) );
138  CHECK( fp_equal_tol( bessel_yn(100,1.), -3.775287810e185, 1.e176 ) );
139  CHECK( fp_equal_tol( bessel_yn(100,100.), -1.669214114e-1, 1.e-10 ) );
140  }
141 
142  TEST(TestBesselI0)
143  {
144  CHECK( fp_equal( bessel_i0_scaled(0.), 1., 10 ) );
145  CHECK( fp_equal_tol( bessel_i0_scaled(1.), 0.4657596077, 2.e-10 ) );
146  CHECK( fp_equal_tol( bessel_i0_scaled(2.9), 0.2477557304, 1.e-10 ) );
147  CHECK( fp_equal_tol( bessel_i0_scaled(4.8), 0.1875862042, 1.e-10 ) );
148  CHECK( fp_equal_tol( bessel_i0_scaled(8.3), 0.1407239098, 1.e-10 ) );
149  CHECK( fp_equal_tol( bessel_i0_scaled(17.4), 0.0963498277, 1.e-10 ) );
150 
151  CHECK( fp_equal( bessel_i0_scaled(1.), exp(-1.)*bessel_i0(1.) ) );
152  CHECK( fp_equal( bessel_i0_scaled(17.4), exp(-17.4)*bessel_i0(17.4) ) );
153  }
154 
155  TEST(TestBesselI1)
156  {
157  CHECK( fp_equal( bessel_i1_scaled(0.), 0., 10 ) );
158  CHECK( fp_equal_tol( bessel_i1_scaled(1.), 0.2079104154, 1.e-10 ) );
159  CHECK( fp_equal_tol( bessel_i1_scaled(2.9), 0.1987772816, 1.e-10 ) );
160  CHECK( fp_equal_tol( bessel_i1_scaled(4.8), 0.1666757058, 1.e-10 ) );
161  CHECK( fp_equal_tol( bessel_i1_scaled(8.3), 0.1319524362, 2.e-10 ) );
162  CHECK( fp_equal_tol( bessel_i1_scaled(17.4), 0.0935388542, 1.e-10 ) );
163 
164  CHECK( fp_equal( bessel_i1_scaled(1.), exp(-1.)*bessel_i1(1.) ) );
165  CHECK( fp_equal( bessel_i1_scaled(17.4), exp(-17.4)*bessel_i1(17.4) ) );
166  }
167 
168  TEST(TestBesselK0)
169  {
170  CHECK( fp_equal_tol( bessel_k0(2.e-30), 68.5003371249, 1.e-10 ) );
171  CHECK( fp_equal_tol( bessel_k0_scaled(1.), 1.1444630797, 2.e-10 ) );
172  CHECK( fp_equal_tol( bessel_k0_scaled(2.9), 0.7089049774, 1.e-10 ) );
173  CHECK( fp_equal_tol( bessel_k0_scaled(4.8), 0.5586133194, 1.e-10 ) );
174  CHECK( fp_equal_tol( bessel_k0_scaled(8.3), 0.4288766329, 1.e-10 ) );
175  CHECK( fp_equal_tol( bessel_k0_scaled(17.4), 0.2983665276, 1.e-10 ) );
176 
177  CHECK( fp_equal( bessel_k0_scaled(1.), exp(1.)*bessel_k0(1.) ) );
178  CHECK( fp_equal( bessel_k0_scaled(17.4), exp(17.4)*bessel_k0(17.4) ) );
179  }
180 
181  TEST(TestBesselK1)
182  {
183  CHECK( fp_equal_tol( bessel_k1(2.e-30), 5e29, 1.e18 ) );
184  CHECK( fp_equal_tol( bessel_k1_scaled(1.), 1.6361534863, 1.e-10 ) );
185  CHECK( fp_equal_tol( bessel_k1_scaled(2.9), 0.8230420403, 1.e-10 ) );
186  CHECK( fp_equal_tol( bessel_k1_scaled(4.8), 0.6142566003, 1.e-10 ) );
187  CHECK( fp_equal_tol( bessel_k1_scaled(8.3), 0.4540139001, 2.e-10 ) );
188  CHECK( fp_equal_tol( bessel_k1_scaled(17.4), 0.3068236027, 1.e-10 ) );
189 
190  CHECK( fp_equal( bessel_k1_scaled(1.), exp(1.)*bessel_k1(1.) ) );
191  CHECK( fp_equal( bessel_k1_scaled(17.4), exp(17.4)*bessel_k1(17.4) ) );
192  }
193 
194  TEST(TestEllpk)
195  {
196  CHECK( fp_equal( ellpk(1.), PI/2., 10 ) );
197  CHECK( fp_equal_tol( ellpk(0.86), 1.630575548881754, 1.e-15 ) );
198  CHECK( fp_equal_tol( ellpk(0.56), 1.806327559107699, 1.e-15 ) );
199  CHECK( fp_equal_tol( ellpk(0.36), 1.995302777664729, 1.e-15 ) );
200  CHECK( fp_equal_tol( ellpk(0.13), 2.455338028321380, 1.e-15 ) );
201  CHECK( fp_equal_tol( ellpk(0.01), 3.695637362989875, 1.e-15 ) );
202  }
203 
204  TEST(TestExpn)
205  {
206  CHECK( fp_equal_tol( expn(1,0.25), 0.25*0.9408157528 - log(0.25) - EULER, 1.e-10 ) );
207  CHECK( fp_equal_tol( expn(1,0.79), 0.316277004, 1.e-9 ) );
208  CHECK( fp_equal_tol( expn(1,1.97), 0.050976988, 1.e-9 ) );
209 
210  CHECK( fp_equal_tol( expn(2,0.25), 0.8643037 + log(0.25)/4., 1.e-7 ) );
211  CHECK( fp_equal_tol( expn(2,0.79), 0.2039860, 1.e-7 ) );
212  CHECK( fp_equal_tol( expn(2,1.97), 0.0390322, 1.e-7 ) );
213 
214  CHECK( fp_equal( expn(3,0.), 0.5, 10 ) );
215  CHECK( fp_equal_tol( expn(3,0.15), 0.3822761, 1.e-7 ) );
216  CHECK( fp_equal_tol( expn(3,0.79), 0.1463479, 1.e-7 ) );
217  CHECK( fp_equal_tol( expn(3,1.97), 0.0312817, 1.e-7 ) );
218 
219  CHECK( fp_equal( expn(4,0.), 1./3., 10 ) );
220  CHECK( fp_equal_tol( expn(4,0.15), 0.2677889, 1.e-7 ) );
221  CHECK( fp_equal_tol( expn(4,0.79), 0.1127433, 1.e-7 ) );
222  CHECK( fp_equal_tol( expn(4,1.97), 0.0259440, 1.e-7 ) );
223 
224  CHECK( fp_equal( expn(10,0.), 1./9., 10 ) );
225  CHECK( fp_equal_tol( expn(10,0.15), 0.0938786, 1.e-7 ) );
226  CHECK( fp_equal_tol( expn(10,0.79), 0.0459453, 1.e-7 ) );
227  CHECK( fp_equal_tol( expn(10,1.97), 0.0124964, 1.e-7 ) );
228  }
229 
230  TEST(TestErf)
231  {
232  /* constants calculated with xmaxima 5.22.1 */
233  CHECK( fp_equal_tol( erf(0.), 0., 1.e-22 ) );
234  // erf(x) loses some precision around 1.e-10, but should still be plenty good...
235  CHECK( fp_equal_tol( erf(1.e-10), 1.1283791671081724525e-10, 3.e-21 ) );
236  CHECK( fp_equal_tol( erf(1.e-5), 1.1283791670579000307e-5, 1.e-17 ) );
237  CHECK( fp_equal_tol( erf(0.1), 1.1246291601828489259e-1, 1.e-13 ) );
238  CHECK( fp_equal_tol( erf(0.5), 5.2049987781304653768e-1, 1.e-12 ) );
239  CHECK( fp_equal_tol( erf(1.), 8.4270079294971486934e-1, 1.e-12 ) );
240  CHECK( fp_equal_tol( erf(2.), 9.9532226501895273416e-1, 1.e-12 ) );
241  CHECK( fp_equal_tol( erf(10.), 1.0, 1.e-12 ) );
242  CHECK( fp_equal( erf(-1.), -erf(1.) ) );
243  }
244 
245  TEST(TestErfc)
246  {
247  /* constants calculated with xmaxima 5.22.1 */
248  CHECK( fp_equal_tol( erfc(-1.), 1.8427007929497148693e0, 1.e-12 ) );
249  CHECK( fp_equal_tol( erfc(0.), 1., 1.e-12 ) );
250  CHECK( fp_equal_tol( erfc(1.e-5), 9.99988716208329421e-1, 1.e-12 ) );
251  CHECK( fp_equal_tol( erfc(0.1), 8.8753708398171510741e-1, 1.e-12 ) );
252  CHECK( fp_equal_tol( erfc(0.5), 4.7950012218695346232e-1, 1.e-12 ) );
253  CHECK( fp_equal_tol( erfc(1.), 1.5729920705028513066e-1, 1.e-13 ) );
254  CHECK( fp_equal_tol( erfc(2.), 4.6777349810472658396e-3, 1.e-15 ) );
255  CHECK( fp_equal_tol( erfc(10.), 2.088487583762544757e-45, 1.e-57 ) );
256  }
257 
258  TEST(TestErfce)
259  {
260  /* constants taken from Finn G.D. & Mugglestone D., 1965, MNRAS 129, 221 */
261  /* this is the voigt function H(a,0) normalized according to 9-44 of Mihalas */
262  CHECK( fp_equal_tol( erfce(0.), 1., 1.e-6 ) );
263  CHECK( fp_equal_tol( erfce(0.01), 9.88815e-1, 1.e-6 ) );
264  CHECK( fp_equal_tol( erfce(0.02), 9.77826e-1, 1.e-6 ) );
265  CHECK( fp_equal_tol( erfce(0.05), 9.45990e-1, 1.e-6 ) );
266  CHECK( fp_equal_tol( erfce(0.10), 8.96456e-1, 1.e-6 ) );
267  CHECK( fp_equal_tol( erfce(0.20), 8.09019e-1, 1.e-6 ) );
268  CHECK( fp_equal_tol( erfce(0.55), 5.90927e-1, 1.e-6 ) );
269  CHECK( fp_equal_tol( erfce(1.00), 4.27583e-1, 1.e-6 ) );
270  }
271 
272  TEST(TestVoigtH0)
273  {
274  /* constants taken from Finn G.D. & Mugglestone D., 1965, MNRAS 129, 221 */
275  /* this is the voigt function H(a,0) normalized according to 9-44 of Mihalas */
276  CHECK( fp_equal_tol( VoigtH0(0.), 1., 1.e-6 ) );
277  CHECK( fp_equal_tol( VoigtH0(0.01), 9.88815e-1, 1.e-6 ) );
278  CHECK( fp_equal_tol( VoigtH0(0.02), 9.77826e-1, 1.e-6 ) );
279  CHECK( fp_equal_tol( VoigtH0(0.05), 9.45990e-1, 1.e-6 ) );
280  CHECK( fp_equal_tol( VoigtH0(0.10), 8.96456e-1, 1.e-6 ) );
281  CHECK( fp_equal_tol( VoigtH0(0.20), 8.09019e-1, 1.e-6 ) );
282  CHECK( fp_equal_tol( VoigtH0(0.55), 5.90927e-1, 1.e-6 ) );
283  CHECK( fp_equal_tol( VoigtH0(1.00), 4.27583e-1, 1.e-6 ) );
284  }
285 
286  TEST(TestVoigtU0)
287  {
288  /* constants taken from Finn G.D. & Mugglestone D., 1965, MNRAS 129, 221 */
289  /* this is the voigt function U(a,0) normalized according to 9-45 of Mihalas */
290  CHECK( fp_equal_tol( VoigtU0(0.), 1./SQRTPI, 1.e-6 ) );
291  CHECK( fp_equal_tol( VoigtU0(0.01), 9.88815e-1/SQRTPI, 1.e-6 ) );
292  CHECK( fp_equal_tol( VoigtU0(0.02), 9.77826e-1/SQRTPI, 1.e-6 ) );
293  CHECK( fp_equal_tol( VoigtU0(0.05), 9.45990e-1/SQRTPI, 1.e-6 ) );
294  CHECK( fp_equal_tol( VoigtU0(0.10), 8.96456e-1/SQRTPI, 1.e-6 ) );
295  CHECK( fp_equal_tol( VoigtU0(0.20), 8.09019e-1/SQRTPI, 1.e-6 ) );
296  CHECK( fp_equal_tol( VoigtU0(0.55), 5.90927e-1/SQRTPI, 1.e-6 ) );
297  CHECK( fp_equal_tol( VoigtU0(1.00), 4.27583e-1/SQRTPI, 1.e-6 ) );
298  }
299 
300  TEST(TestVoigtH)
301  {
302  // check that the Voigt profile returned by VoigtH() is properly normalized
303  const int NP = 200;
304  realnum v[NP], a, y[NP];
305  // for a > 0.1, VoigtH() calls humlik(), for smaller values it
306  // calls FastVoigtH().
307  //
308  // humlik() is set up for a relative precision of 1e-4 over its
309  // entire range, but looks to be more precise in practice (at
310  // least at v=0)
311  // FastVoigtH() is set up for a rel. precision of 2.5e-3 over its
312  // entire range, but should be better than 1e-4 for a < 0.0235
313  a = realnum(0.0002);
314  for( int i=0; i < 9; ++i )
315  {
316  // test both humlik() and FastVoigtH()
317  for( int i=0; i < NP; ++i )
318  v[i] = realnum(i)*max(a,1.f)/realnum(5.);
319  VoigtH( a, v, y, NP );
320  realnum integral = realnum(0.);
321  // We need the integral from -infinity to +infinity, which is simply
322  // two times the integral from 0 to +infinity. Hence we omit the
323  // division by 2 in the trapezoidal rule
324  for( int i=1; i < NP; ++i )
325  integral += (y[i]+y[i-1])*(v[i]-v[i-1]);
326  // add in the integral over the Lorentz wings assuming H(a,v) = c/v^2
327  integral += realnum(2.)*v[NP-1]*y[NP-1];
328  // VoigtH() calculates H(a,v), so integral should be sqrt(pi)
329  CHECK( fp_equal_tol( integral, realnum(SQRTPI), realnum(1.e-4) ) );
330  // also check the central value...
331  realnum h0 = realnum(VoigtH0(a));
332  CHECK( fp_equal_tol( y[0], h0, realnum(1.e-4)*h0 ) );
333  a *= realnum(10.);
334  }
335  // now do some spot checks
336  // constants taken from
337  // >>refer Zaghloul M.R. & Ali A.N. 2011, ACM Transactions on Mathematical Software, 38, 15
338  // available from http://arxiv.org/abs/1106.0151
339  v[0] = realnum(5.76);
340  VoigtH( realnum(1.e-20), v, y, 1 );
341  // this constant comes from page 21 ("Present algorithm").
342  CHECK( fp_equal_tol( y[0], realnum(3.900779639194698e-015), realnum(4.e-19) ) );
343  v[0] = realnum(6.3e-2);
344  v[1] = realnum(6.3e+0);
345  v[2] = realnum(6.3e+2);
346  VoigtH( realnum(1.e-20), v, y, 3 );
347  // these constants come from Table 2 of the same paper ("Present algorithm").
348  CHECK( fp_equal_tol( y[0], realnum(9.960388660702479e-001), realnum(1.e-4) ) );
349  CHECK( fp_equal_tol( y[1], realnum(5.792460778844116e-018), realnum(6.e-22) ) );
350  CHECK( fp_equal_tol( y[2], realnum(1.421495882582395e-026), realnum(1.4e-30) ) );
351  VoigtH( realnum(1.e-14), v, y, 3 );
352  CHECK( fp_equal_tol( y[0], realnum(9.960388660702366e-001), realnum(1.e-4) ) );
353  CHECK( fp_equal_tol( y[1], realnum(1.536857621303163e-016), realnum(1.5e-20) ) );
354  CHECK( fp_equal_tol( y[2], realnum(1.421495882582395e-020), realnum(1.4e-24) ) );
355  VoigtH( realnum(1.e-12), v, y, 3 );
356  CHECK( fp_equal_tol( y[0], realnum(9.960388660691284e-001), realnum(1.e-4) ) );
357  CHECK( fp_equal_tol( y[1], realnum(1.479513723737753e-014), realnum(1.5e-18) ) );
358  CHECK( fp_equal_tol( y[2], realnum(1.421495882582395e-018), realnum(1.4e-22) ) );
359  VoigtH( realnum(1.e-10), v, y, 3 );
360  CHECK( fp_equal_tol( y[0], realnum(9.960388659583033e-001), realnum(1.e-4) ) );
361  CHECK( fp_equal_tol( y[1], realnum(1.478940284762099e-012), realnum(1.5e-16) ) );
362  CHECK( fp_equal_tol( y[2], realnum(1.421495882582395e-016), realnum(1.4e-20) ) );
363  VoigtH( realnum(1.e-6), v, y, 3 );
364  CHECK( fp_equal_tol( y[0], realnum(9.960377466254801e-001), realnum(1.e-4) ) );
365  CHECK( fp_equal_tol( y[1], realnum(1.478934493028404e-008), realnum(1.5e-12) ) );
366  CHECK( fp_equal_tol( y[2], realnum(1.421495882582395e-012), realnum(1.4e-16) ) );
367  VoigtH( realnum(1.e-2), v, y, 3 );
368  CHECK( fp_equal_tol( y[0], realnum(9.849424862549039e-001), realnum(1.e-4) ) );
369  CHECK( fp_equal_tol( y[1], realnum(1.478930389133934e-004), realnum(1.5e-8) ) );
370  CHECK( fp_equal_tol( y[2], realnum(1.421495882224242e-008), realnum(1.4e-12) ) );
371  VoigtH( realnum(1.e+1), v, y, 3 );
372  CHECK( fp_equal_tol( y[0], realnum(5.613881832823886e-002), realnum(6.e-6) ) );
373  CHECK( fp_equal_tol( y[1], realnum(4.040671157393835e-002), realnum(4.e-6) ) );
374  CHECK( fp_equal_tol( y[2], realnum(1.421137820009847e-005), realnum(1.4e-9) ) );
375  VoigtH( realnum(1.2e+1), v, y, 3 );
376  CHECK( fp_equal_tol( y[0], realnum(4.685295149211637e-002), realnum(5.e-6) ) );
377  CHECK( fp_equal_tol( y[1], realnum(3.684277239564798e-002), realnum(4.e-6) ) );
378  CHECK( fp_equal_tol( y[2], realnum(1.705176395541707e-005), realnum(1.7e-9) ) );
379  VoigtH( realnum(1.5e+1), v, y, 3 );
380  CHECK( fp_equal_tol( y[0], realnum(3.752895161491574e-002), realnum(4.e-6) ) );
381  CHECK( fp_equal_tol( y[1], realnum(3.194834330452605e-002), realnum(3.e-6) ) );
382  CHECK( fp_equal_tol( y[2], realnum(2.131035743074598e-005), realnum(2.e-9) ) );
383  VoigtH( realnum(2.e+2), v, y, 3 );
384  CHECK( fp_equal_tol( y[0], realnum(2.820912377324508e-003), realnum(3.e-7) ) );
385  CHECK( fp_equal_tol( y[1], realnum(2.818116555672206e-003), realnum(3.e-7) ) );
386  CHECK( fp_equal_tol( y[2], realnum(2.582702147491469e-004), realnum(3.e-8) ) );
387  VoigtH( realnum(1.e+5), v, y, 3 );
388  CHECK( fp_equal_tol( y[0], realnum(5.641895835193228e-006), realnum(6.e-10) ) );
389  CHECK( fp_equal_tol( y[1], realnum(5.641895812802746e-006), realnum(6.e-10) ) );
390  CHECK( fp_equal_tol( y[2], realnum(5.641671917237128e-006), realnum(6.e-10) ) );
391  v[0] = realnum(1.e+0);
392  VoigtH( realnum(1.e-20), v, y, 1 );
393  CHECK( fp_equal_tol( y[0], realnum(3.678794411714423e-001), realnum(4.e-5) ) );
394  v[0] = realnum(5.5e+0);
395  VoigtH( realnum(1.e-14), v, y, 1 );
396  CHECK( fp_equal_tol( y[0], realnum(7.307386729528773e-014), realnum(7.e-18) ) );
397  v[0] = realnum(3.9e+4);
398  VoigtH( realnum(1.e+0), v, y, 1 );
399  CHECK( fp_equal_tol( y[0], realnum(3.709333226385423e-010), realnum(4.e-14) ) );
400  v[0] = realnum(1.e+0);
401  VoigtH( realnum(2.8e+4), v, y, 1 );
402  CHECK( fp_equal_tol( y[0], realnum(2.014962794529686e-005), realnum(2.e-9) ) );
403  }
404 
405  TEST(TestVoigtU)
406  {
407  // check that the Voigt profile returned by VoigtU() is properly normalized
408  const int NP = 200;
409  realnum v[NP], a, y[NP];
410  // for a > 0.1, VoigtU() calls humlik(), for smaller values it
411  // calls FastVoigtH() and divides by sqrt(pi).
412  //
413  // humlik() is set up for a relative precision of 1e-4 over its
414  // entire range, but looks to be more precise in practice (at
415  // least at v=0)
416  // FastVoigtH() is set up for a rel. precision of 2.5e-3 over its
417  // entire range, but should be better than 1e-4 for a < 0.0235
418  a = realnum(0.0002);
419  for( int i=0; i < 9; ++i )
420  {
421  // test both humlik() and FastVoigtH()
422  for( int i=0; i < NP; ++i )
423  v[i] = realnum(i)*max(a,1.f)/realnum(5.);
424  VoigtU( a, v, y, NP );
425  realnum integral = realnum(0.);
426  // We need the integral from -infinity to +infinity, which is simply
427  // two times the integral from 0 to +infinity. Hence we omit the
428  // division by 2 in the trapezoidal rule
429  for( int i=1; i < NP; ++i )
430  integral += (y[i]+y[i-1])*(v[i]-v[i-1]);
431  // add in the integral over the Lorentz wings assuming U(a,v) = c/v^2
432  integral += realnum(2.)*v[NP-1]*y[NP-1];
433  // VoigtU() calculates U(a,v), so integral should be 1
434  CHECK( fp_equal_tol( integral, realnum(1.), realnum(1.e-4) ) );
435  // also check the central value...
436  CHECK( fp_equal_tol( y[0], realnum(VoigtU0(a)), realnum(1.e-4) ) );
437  a *= realnum(10.);
438  }
439  }
440 
441  TEST(TestMD5string)
442  {
443  string test;
444  // md5sum of an empty file...
445  CHECK( MD5string( test ) == "d41d8cd98f00b204e9800998ecf8427e" );
446  CHECK( MD5string( test ).length() == NMD5 );
447  // check if padding is done correctly
448  // an extra block of padding needs to be added when length%64 == 56
449  test = "hhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhh";
450  CHECK( test.length() == 55 );
451  CHECK( MD5string( test ) == "426ec4ac35ad38d125f6efb39da03098" );
452  test = "hhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhh";
453  CHECK( test.length() == 56 );
454  CHECK( MD5string( test ) == "d03607b2c89adc0c4abf5a0f1d9e40c9" );
455  test = "hhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhh";
456  CHECK( test.length() == 57 );
457  CHECK( MD5string( test ) == "bac1b47748411cb6eee0cae3befb8377" );
458  string test64 = "hhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhh";
459  test = test64 + "hhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhh";
460  CHECK( test.length() == 64+55 );
461  CHECK( MD5string( test ) == "10d49aad1fc69976376fbe7c8c5ed118" );
462  test = test64 + "hhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhh";
463  CHECK( test.length() == 64+56 );
464  CHECK( MD5string( test ) == "61ec7da14576f3b585038c6d72cd5bd5" );
465  test = test64 + "hhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhh";
466  CHECK( test.length() == 64+57 );
467  CHECK( MD5string( test ) == "f17a0475a26d0930e2a35bb320c10e0d" );
468  // check that leading zeros are printed correctly
469  test = "ghijklmn";
470  CHECK( MD5string( test ) == "0256b9cea63bc1f97b8c5aea92c24a98" );
471  }
472 
473 }
bessel_k0_scaled
double bessel_k0_scaled(double x)
Definition: thirdparty.cpp:1382
bessel_i1_scaled
double bessel_i1_scaled(double x)
Definition: thirdparty.cpp:1929
lfactorial
double lfactorial(long n)
Definition: thirdparty.cpp:399
bessel_i1
double bessel_i1(double x)
Definition: thirdparty.cpp:1908
fp_equal_tol
bool fp_equal_tol(sys_float x, sys_float y, sys_float tol)
Definition: cddefines.h:858
expn
double expn(int n, double x)
Definition: thirdparty.cpp:2121
cdgamma
complex< double > cdgamma(complex< double > x)
Definition: thirdparty.cpp:432
fp_equal
bool fp_equal(sys_float x, sys_float y, int n=3)
Definition: cddefines.h:816
VoigtU0
double VoigtU0(double a)
Definition: thirdparty.h:378
bessel_k1
double bessel_k1(double x)
Definition: thirdparty.cpp:1535
erfc
double erfc(double)
VoigtH0
double VoigtH0(double a)
Definition: thirdparty.h:372
bessel_i0_scaled
double bessel_i0_scaled(double x)
Definition: thirdparty.cpp:1743
bessel_i0
double bessel_i0(double x)
Definition: thirdparty.cpp:1726
bessel_jn
double bessel_jn(int n, double x)
Definition: thirdparty.cpp:1042
realnum
float realnum
Definition: cddefines.h:107
bessel_y0
double bessel_y0(double x)
Definition: thirdparty.cpp:746
max
long max(int a, long b)
Definition: cddefines.h:779
bessel_j0
double bessel_j0(double x)
Definition: thirdparty.cpp:708
bessel_j1
double bessel_j1(double x)
Definition: thirdparty.cpp:939
factorial
double factorial(long n)
Definition: thirdparty.cpp:356
bessel_yn
double bessel_yn(int n, double x)
Definition: thirdparty.cpp:1177
bessel_k1_scaled
double bessel_k1_scaled(double x)
Definition: thirdparty.cpp:1557
erfce
double erfce(double x)
Definition: thirdparty.cpp:2454
VoigtU
void VoigtU(realnum a, const realnum v[], realnum y[], int n)
Definition: thirdparty.h:364
SQRTPI
const double SQRTPI
Definition: physconst.h:44
EULER
const double EULER
Definition: physconst.h:26
bessel_y1
double bessel_y1(double x)
Definition: thirdparty.cpp:966
NMD5
static const unsigned int NMD5
Definition: thirdparty.h:384
erf
double erf(double)
ellpk
double ellpk(double x)
Definition: thirdparty.cpp:2041
bessel_k0
double bessel_k0(double x)
Definition: thirdparty.cpp:1359
e2
double e2(double t)
Definition: service.cpp:299
MD5string
string MD5string(const string &str)
Definition: thirdparty.cpp:3461
VoigtH
void VoigtH(realnum a, const realnum v[], realnum y[], int n)
Definition: thirdparty.h:350
PI
const double PI
Definition: physconst.h:29
Generated on Sat Oct 19 2013 10:43:14 for cloudy by   doxygen 1.8.5