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-2022 by Gary J. Ferland and
2  * others. For conditions of distribution and use see copyright notice in license.txt */
3 #include "cdstd.h"
4 #include <UnitTest++.h>
5 #include "cddefines.h"
6 #include "thirdparty.h"
7 
8 namespace {
9 
10  // constants below were calculated with xmaxima 5.17.0 with fpprec set to 30
11 
12  TEST(TestFactorial)
13  {
14  CHECK( fp_equal( factorial(0), 1. ) );
15  CHECK( fp_equal( factorial(1), 1. ) );
16  CHECK( fp_equal( factorial(2), 2. ) );
17  CHECK( fp_equal( factorial(10), 3628800. ) );
18  CHECK( fp_equal( factorial(170), 7.257415615307998967396728211129316e306 ) );
19  // Validate against standard library function
20  CHECK( fp_equal( factorial(170), tgamma(171.), 40 ) );
21  }
22 
23  TEST(TestLogFactorial)
24  {
25  CHECK( fp_equal( lfactorial(0), 0. ) );
26  CHECK( fp_equal( lfactorial(1), 0. ) );
27  CHECK( fp_equal( lfactorial(2), 3.01029995663981195213738894725e-1, 10 ) );
28  CHECK( fp_equal( lfactorial(10), 6.55976303287679375117476123996e0, 10 ) );
29  CHECK( fp_equal( lfactorial(170), 3.06860781994828320192380273111e2, 10 ) );
30  CHECK( fp_equal( lfactorial(1700), 4.75547688279135071178908638597e3, 10 ) );
31  // Validate against standard library function
32  CHECK( fp_equal( lfactorial(1700), lgamma(1701.)/log(10.), 10 ) );
33  }
34 
35  TEST(TestCDGamma)
36  {
37  complex<double> y = cdgamma( complex<double>(1.,0.) );
38  CHECK( fp_equal( y.real(), 1., 10 ) && fp_equal( y.imag(), 0. ) );
39  y = cdgamma( complex<double>(2.,0.) );
40  CHECK( fp_equal( y.real(), 1., 10 ) && fp_equal( y.imag(), 0. ) );
41  y = cdgamma( complex<double>(11.,0.) );
42  CHECK( fp_equal( y.real(), 3628800., 50 ) && fp_equal( y.imag(), 0. ) );
43  y = cdgamma( complex<double>(-0.5,0.) );
44  CHECK( fp_equal( y.real(), -3.544907701811032054596334966682277e0, 10 ) &&
45  fp_equal( y.imag(), 0. ) );
46  y = cdgamma( complex<double>(0.,1.) );
47  CHECK( fp_equal( y.real(), -1.549498283018106851249551304838863e-1, 10 ) &&
48  fp_equal( y.imag(), -4.980156681183560427136911174621973e-1, 10 ) );
49  y = cdgamma( complex<double>(-1.,-2.) );
50  CHECK( fp_equal( y.real(), -3.23612885501927256868232032760969e-2, 30 ) &&
51  fp_equal( y.imag(), -1.122942423463261735043406872030743e-2, 30 ) );
52  }
53 
54  // constants below were taken from Abramowitz & Stegun, Handbook of Mathematical Functions
55 
56  TEST(TestBesselJ0)
57  {
58  CHECK( fp_equal( bessel_j0(0.), 1., 10 ) );
59  CHECK( fp_equal_tol( bessel_j0(1.), 0.765197686557967, 1.e-15 ) );
60  CHECK( fp_equal_tol( bessel_j0(2.9), -0.224311545791968, 1.e-15 ) );
61  CHECK( fp_equal_tol( bessel_j0(4.8), -0.240425327291183, 1.e-15 ) );
62  CHECK( fp_equal_tol( bessel_j0(8.3), 0.096006100895010, 1.e-15 ) );
63  CHECK( fp_equal_tol( bessel_j0(17.4), -0.118955856336348, 1.e-15 ) );
64  }
65 
66  TEST(TestBesselJ1)
67  {
68  CHECK( fp_equal( bessel_j1(0.), 0., 10 ) );
69  CHECK( fp_equal( bessel_j1(1.e-30), 5.e-31, 10 ) );
70  CHECK( fp_equal_tol( bessel_j1(1.e-15), 5.e-16, 1.e-27 ) );
71  CHECK( fp_equal_tol( bessel_j1(1.), 0.4400505857, 1.e-10 ) );
72  CHECK( fp_equal_tol( bessel_j1(2.9), 0.3754274818, 1.e-10 ) );
73  CHECK( fp_equal_tol( bessel_j1(4.8), -0.2984998581, 1.e-10 ) );
74  CHECK( fp_equal_tol( bessel_j1(8.3), 0.2657393020, 1.e-10 ) );
75  CHECK( fp_equal_tol( bessel_j1(17.4), -0.1532161760, 1.e-10 ) );
76  }
77 
78  TEST(TestBesselJn)
79  {
80  CHECK( fp_equal( bessel_jn(2,0.), 0., 10 ) );
81  CHECK( fp_equal( bessel_jn(2,1.e-30), 1.25e-61, 10 ) );
82  CHECK( fp_equal_tol( bessel_jn(2,2.e-15), 5.e-31, 1.e-42 ) );
83  CHECK( fp_equal_tol( bessel_jn(2,2.e-10), 5.e-21, 1.e-32 ) );
84  CHECK( fp_equal_tol( bessel_jn(2,0.099999999999), 0.0012489587, 1.e-10 ) );
85  CHECK( fp_equal_tol( bessel_jn(2,0.100000000001), 0.0012489587, 1.e-10 ) );
86  CHECK( fp_equal_tol( bessel_jn(2,1.), 0.1149034849, 1.e-10 ) );
87  CHECK( fp_equal_tol( bessel_jn(2,2.9), 0.4832270505, 1.e-10 ) );
88  CHECK( fp_equal_tol( bessel_jn(2,4.8), 0.1160503864, 1.e-10 ) );
89  CHECK( fp_equal_tol( bessel_jn(2,8.3), -0.0319725341, 1.e-10 ) );
90  CHECK( fp_equal_tol( bessel_jn(2,17.4), 0.1013448016, 1.e-10 ) );
91 
92  CHECK( fp_equal( bessel_jn(8,0.), 0., 10 ) );
93  CHECK( fp_equal( bessel_jn(8,2.e-20), 1.e-160/40320., 10 ) );
94  CHECK( fp_equal_tol( bessel_jn(8,2.e-15), 1.e-120/40320., 1.e-136 ) );
95  CHECK( fp_equal_tol( bessel_jn(8,1.), 9.4223e-8, 1.e-12 ) );
96  CHECK( fp_equal_tol( bessel_jn(8,2.8), 2.9367e-4, 1.e-8 ) );
97  CHECK( fp_equal_tol( bessel_jn(8,4.8), 1.4079e-2, 1.e-6 ) );
98  CHECK( fp_equal_tol( bessel_jn(8,8.2), 0.24257, 1.e-5 ) );
99 
100  CHECK( fp_equal_tol( bessel_jn(20,1.), 3.873503009e-25, 1.e-34 ) );
101  CHECK( fp_equal_tol( bessel_jn(50,1.), 2.906004948e-80, 1.e-89 ) );
102  CHECK( fp_equal_tol( bessel_jn(100,1.), 8.431828790e-189, 1.e-198 ) );
103  CHECK( fp_equal_tol( bessel_jn(100,100.), 9.636667330e-2, 1.e-11 ) );
104  }
105 
106  TEST(TestBesselY0)
107  {
108  CHECK( fp_equal_tol( bessel_y0(1.e-30), -44.0499402279, 1.e-10 ) );
109  CHECK( fp_equal_tol( bessel_y0(1.), 0.0882569642, 1.e-10 ) );
110  CHECK( fp_equal_tol( bessel_y0(2.9), 0.4079117692, 1.e-10 ) );
111  CHECK( fp_equal_tol( bessel_y0(4.8), -0.2723037945, 1.e-10 ) );
112  CHECK( fp_equal_tol( bessel_y0(8.3), 0.2595149638, 1.e-10 ) );
113  CHECK( fp_equal_tol( bessel_y0(17.4), -0.1497391883, 1.e-10 ) );
114  }
115 
116  TEST(TestBesselY1)
117  {
118  CHECK( fp_equal_tol( bessel_y1(1.e-30), -6.36619772368e29, 1.e18 ) );
119  CHECK( fp_equal_tol( bessel_y1(1.), -0.7812128213, 1.e-10 ) );
120  CHECK( fp_equal_tol( bessel_y1(2.9), 0.2959400546, 1.e-10 ) );
121  CHECK( fp_equal_tol( bessel_y1(4.8), 0.2135651673, 1.e-10 ) );
122  CHECK( fp_equal_tol( bessel_y1(8.3), -0.0805975035, 1.e-10 ) );
123  CHECK( fp_equal_tol( bessel_y1(17.4), 0.1147053859, 1.e-10 ) );
124  }
125 
126  TEST(TestBesselYn)
127  {
128  CHECK( fp_equal_tol( bessel_yn(2,1.e-30), -1.27323954474e60, 1.e49 ) );
129  CHECK( fp_equal_tol( bessel_yn(2,1.), -1.65068261, 1.e-8 ) );
130  CHECK( fp_equal_tol( bessel_yn(2,2.9), -0.20381518, 1.e-8 ) );
131  CHECK( fp_equal_tol( bessel_yn(2,4.8), 0.36128928, 1.e-8 ) );
132  CHECK( fp_equal_tol( bessel_yn(2,8.3), -0.27893605, 1.e-8 ) );
133  CHECK( fp_equal_tol( bessel_yn(2,17.4), 0.16292372, 1.e-8 ) );
134 
135  CHECK( fp_equal_tol( bessel_yn(8,2.e-20), -1.60428182637e163, 1.e152 ) );
136  CHECK( fp_equal_tol( bessel_yn(8,1.), -4.2567e5, 1.e1 ) );
137  CHECK( fp_equal_tol( bessel_yn(8,2.8), -1.4486e2, 1.e-2 ) );
138  CHECK( fp_equal_tol( bessel_yn(8,4.8), -3.5855, 1.e-4 ) );
139  CHECK( fp_equal_tol( bessel_yn(8,8.2), -0.35049, 1.e-5 ) );
140 
141  CHECK( fp_equal_tol( bessel_yn(20,1.), -4.113970315e22, 1.e13 ) );
142  CHECK( fp_equal_tol( bessel_yn(50,1.), -2.191142813e77, 1.e68 ) );
143  CHECK( fp_equal_tol( bessel_yn(100,1.), -3.775287810e185, 1.e176 ) );
144  CHECK( fp_equal_tol( bessel_yn(100,100.), -1.669214114e-1, 1.e-10 ) );
145  }
146 
147  // constants below were calculated using http://keisan.casio.com/exec/system/1180573473
148 
149  TEST(TestBesselI0)
150  {
151  CHECK( fp_equal( bessel_i0(0.), 1. ) );
152  CHECK( fp_equal( bessel_i0(1.), 1.266065877752008335598 ) );
153  CHECK( fp_equal( bessel_i0(2.9), 4.50274866132627436631 ) );
154  CHECK( fp_equal( bessel_i0(4.8), 22.79367799310579796012 ) );
155  CHECK( fp_equal( bessel_i0(8.3), 566.255055846491838162, 10 ) );
156  CHECK( fp_equal( bessel_i0(17.4), 3471961.712717317753354, 10 ) );
157 
158  CHECK( fp_equal( bessel_i0(-1.), 1.266065877752008335598 ) );
159  CHECK( fp_equal( bessel_i0(-17.4), 3471961.712717317753354, 10 ) );
160 
161  CHECK( fp_equal( bessel_i0_scaled(1.), exp(-1.)*bessel_i0(1.) ) );
162  CHECK( fp_equal( bessel_i0_scaled(17.4), exp(-17.4)*bessel_i0(17.4) ) );
163 
164  CHECK( fp_equal( bessel_i0_scaled(-1.), exp(-1.)*bessel_i0(-1.) ) );
165  CHECK( fp_equal( bessel_i0_scaled(-17.4), exp(-17.4)*bessel_i0(-17.4) ) );
166  }
167 
168  TEST(TestBesselI1)
169  {
170  CHECK( fp_equal( bessel_i1(0.), 0. ) );
171  CHECK( fp_equal( bessel_i1(1.), 0.5651591039924850272077 ) );
172  CHECK( fp_equal( bessel_i1(2.9), 3.6126072124369077367 ) );
173  CHECK( fp_equal( bessel_i1(4.8), 20.25283460023855998949 ) );
174  CHECK( fp_equal( bessel_i1(8.3), 530.959765452190423905, 7 ) );
175  CHECK( fp_equal( bessel_i1(17.4), 3370668.4094771099123, 10 ) );
176 
177  CHECK( fp_equal( bessel_i1(-1.), -0.5651591039924850272077 ) );
178  CHECK( fp_equal( bessel_i1(-17.4), -3370668.4094771099123, 10 ) );
179 
180  CHECK( fp_equal( bessel_i1_scaled(1.), exp(-1.)*bessel_i1(1.) ) );
181  CHECK( fp_equal( bessel_i1_scaled(17.4), exp(-17.4)*bessel_i1(17.4) ) );
182 
183  CHECK( fp_equal( bessel_i1_scaled(-1.), exp(-1.)*bessel_i1(-1.) ) );
184  CHECK( fp_equal( bessel_i1_scaled(-17.4), exp(-17.4)*bessel_i1(-17.4) ) );
185  }
186 
187  TEST(TestBesselI0I1)
188  {
189  double i0val, i1val;
190  bessel_i0_i1(0.3, &i0val, &i1val);
191  CHECK( fp_equal( i0val, 1.02262687935159699112 ) );
192  CHECK( fp_equal( i1val, 0.1516938400035927803287 ) );
193  bessel_i0_i1_scaled(0.3, &i0val, &i1val);
194  CHECK( fp_equal( i0val, 1.02262687935159699112*exp(-0.3) ) );
195  CHECK( fp_equal( i1val, 0.1516938400035927803287*exp(-0.3) ) );
196 
197  bessel_i0_i1(-0.3, &i0val, &i1val);
198  CHECK( fp_equal( i0val, 1.02262687935159699112 ) );
199  CHECK( fp_equal( i1val, -0.1516938400035927803287 ) );
200  bessel_i0_i1_scaled(-0.3, &i0val, &i1val);
201  CHECK( fp_equal( i0val, 1.02262687935159699112*exp(-0.3) ) );
202  CHECK( fp_equal( i1val, -0.1516938400035927803287*exp(-0.3) ) );
203 
204  bessel_i0_i1(22., &i0val, &i1val);
205  CHECK( fp_equal( i0val, 306692993.6403647456135 ) );
206  CHECK( fp_equal( i1val, 299639606.8773789797263 ) );
207  bessel_i0_i1_scaled(22., &i0val, &i1val);
208  CHECK( fp_equal( i0val, 306692993.6403647456135*exp(-22.) ) );
209  CHECK( fp_equal( i1val, 299639606.8773789797263*exp(-22.) ) );
210 
211  bessel_i0_i1(-22., &i0val, &i1val);
212  CHECK( fp_equal( i0val, 306692993.6403647456135 ) );
213  CHECK( fp_equal( i1val, -299639606.8773789797263 ) );
214  bessel_i0_i1_scaled(-22., &i0val, &i1val);
215  CHECK( fp_equal( i0val, 306692993.6403647456135*exp(-22.) ) );
216  CHECK( fp_equal( i1val, -299639606.8773789797263*exp(-22.) ) );
217  }
218 
219  TEST(TestBesselK0)
220  {
221  CHECK_THROW( bessel_k0(0.), domain_error );
222  CHECK_THROW( bessel_k0(-1.), domain_error );
223  CHECK( fp_equal( bessel_k0(2.e-30), 68.50033712491983765993 ) );
224  CHECK( fp_equal( bessel_k0(1.), 0.4210244382407083333356 ) );
225  CHECK( fp_equal( bessel_k0(2.9), 0.0390062345662234241006 ) );
226  CHECK( fp_equal( bessel_k0(4.8), 0.00459724631672465789944 ) );
227  CHECK( fp_equal( bessel_k0(8.3), 1.065830600438036175396e-4, 5 ) );
228  CHECK( fp_equal( bessel_k0(17.4), 8.279919509749720001621e-9, 10 ) );
229 
230  CHECK( fp_equal( bessel_k0_scaled(1.), exp(1.)*bessel_k0(1.) ) );
231  CHECK( fp_equal( bessel_k0_scaled(17.4), exp(17.4)*bessel_k0(17.4) ) );
232  }
233 
234  TEST(TestBesselK1)
235  {
236  CHECK_THROW( bessel_k1(0.), domain_error );
237  CHECK_THROW( bessel_k1(-1.), domain_error );
238  CHECK( fp_equal( bessel_k1(2.e-30), 5.e29 ) );
239  CHECK( fp_equal( bessel_k1(1.), 0.6019072301972345747375 ) );
240  CHECK( fp_equal( bessel_k1(2.9), 0.0452864232983614435612 ) );
241  CHECK( fp_equal( bessel_k1(4.8), 0.00505517644405629981585 ) );
242  CHECK( fp_equal( bessel_k1(8.3), 1.128300939464441907577e-4 ) );
243  CHECK( fp_equal( bessel_k1(17.4), 8.514610381504642118954e-9, 10 ) );
244 
245  CHECK( fp_equal( bessel_k1_scaled(1.), exp(1.)*bessel_k1(1.) ) );
246  CHECK( fp_equal( bessel_k1_scaled(17.4), exp(17.4)*bessel_k1(17.4) ) );
247  }
248 
249  TEST(TestBesselK0K1)
250  {
251  double k0val, k1val;
252  CHECK_THROW( bessel_k0_k1(0., &k0val, &k1val), domain_error );
253  CHECK_THROW( bessel_k0_k1(-1., &k0val, &k1val), domain_error );
254  bessel_k0_k1(0.3, &k0val, &k1val);
255  CHECK( fp_equal( k0val, 1.372460060544297376645 ) );
256  CHECK( fp_equal( k1val, 3.055992033457324978851 ) );
257  bessel_k0_k1_scaled(0.3, &k0val, &k1val);
258  CHECK( fp_equal( k0val, 1.372460060544297376645*exp(0.3) ) );
259  CHECK( fp_equal( k1val, 3.055992033457324978851*exp(0.3) ) );
260 
261  bessel_k0_k1(22., &k0val, &k1val);
262  CHECK( fp_equal( k0val, 7.412351614084865368946e-11 ) );
263  CHECK( fp_equal( k1val, 7.578981163485331089804e-11 ) );
264  bessel_k0_k1_scaled(22., &k0val, &k1val);
265  CHECK( fp_equal( k0val, 7.412351614084865368946e-11*exp(22.) ) );
266  CHECK( fp_equal( k1val, 7.578981163485331089804e-11*exp(22.) ) );
267  }
268 
269  // constants below were taken from Abramowitz & Stegun, Handbook of Mathematical Functions
270 
271  TEST(TestEllpk)
272  {
273  CHECK( fp_equal( ellpk(1.), PI/2., 10 ) );
274  CHECK( fp_equal_tol( ellpk(0.86), 1.630575548881754, 1.e-15 ) );
275  CHECK( fp_equal_tol( ellpk(0.56), 1.806327559107699, 1.e-15 ) );
276  CHECK( fp_equal_tol( ellpk(0.36), 1.995302777664729, 1.e-15 ) );
277  CHECK( fp_equal_tol( ellpk(0.13), 2.455338028321380, 1.e-15 ) );
278  CHECK( fp_equal_tol( ellpk(0.01), 3.695637362989875, 1.e-15 ) );
279  }
280 
281  // constants below were calculated with xmaxima 5.34.1 with fpprec set to 30
282 
283  TEST(TestIgam)
284  {
285  CHECK( fp_equal_tol( igam(0.3,0.2), 0.65750672426972171, 1.e-15 ) );
286  CHECK( fp_equal_tol( igam(0.46,0.8), 0.81316729720991554, 1.e-15 ) );
287  CHECK( fp_equal_tol( igam(0.67,2.2), 0.94307235949930435, 1.e-15 ) );
288  CHECK( fp_equal_tol( igam(1.2,4.5), 0.98302397308764477, 1.e-15 ) );
289  CHECK( fp_equal_tol( igam(1.,10.), 0.99995460007023752, 1.e-15 ) );
290  CHECK( fp_equal_tol( igam(0.3,1e-5), 0.03523536061556257, 1.e-15 ) );
291  }
292 
293  TEST(TestIgamc)
294  {
295  CHECK( fp_equal_tol( igamc(0.3,0.2), 1.-0.65750672426972171, 1.e-15 ) );
296  CHECK( fp_equal_tol( igamc(0.46,0.8), 1.-0.81316729720991554, 1.e-15 ) );
297  CHECK( fp_equal_tol( igamc(0.67,2.2), 1.-0.94307235949930435, 1.e-15 ) );
298  CHECK( fp_equal_tol( igamc(1.2,4.5), 1.-0.98302397308764477, 1.e-15 ) );
299  CHECK( fp_equal_tol( igamc(1.,10.), 1.-0.99995460007023752, 1.e-15 ) );
300  CHECK( fp_equal_tol( igamc(0.3,1e-5), 1.-0.03523536061556257, 1.e-15 ) );
301  }
302 
303  TEST(TestIgamc_scaled)
304  {
305  CHECK( fp_equal_tol( igamc_scaled(0.3,0.2), 0.41832223162827345394, 1.e-15 ) );
306  CHECK( fp_equal_tol( igamc_scaled(0.46,0.8), 0.41580382684020182212, 1.e-15 ) );
307  CHECK( fp_equal_tol( igamc_scaled(0.67,2.2), 0.51377272400971086105, 1.e-15 ) );
308  CHECK( fp_equal_tol( igamc_scaled(1.2,4.5), 1.52813324353067303945, 1.5e-15 ) );
309  CHECK( fp_equal_tol( igamc_scaled(1.,10.), 1., 1.e-15 ) );
310  CHECK( fp_equal_tol( igamc_scaled(0.3,1e5), 1.057055858519714446e-4 , 1.e-15 ) );
311  }
312 
313  // constants below were taken from Abramowitz & Stegun, Handbook of Mathematical Functions
314 
315  TEST(TestExpn)
316  {
317  CHECK( fp_equal_tol( expn(1,0.25), 0.25*0.9408157528 - log(0.25) - EULER, 1.e-10 ) );
318  CHECK( fp_equal_tol( expn(1,0.79), 0.316277004, 1.e-9 ) );
319  CHECK( fp_equal_tol( expn(1,1.97), 0.050976988, 1.e-9 ) );
320 
321  CHECK( fp_equal_tol( expn(2,0.25), 0.8643037 + log(0.25)/4., 1.e-7 ) );
322  CHECK( fp_equal_tol( expn(2,0.79), 0.2039860, 1.e-7 ) );
323  CHECK( fp_equal_tol( expn(2,1.97), 0.0390322, 1.e-7 ) );
324 
325  CHECK( fp_equal( expn(3,0.), 0.5, 10 ) );
326  CHECK( fp_equal_tol( expn(3,0.15), 0.3822761, 1.e-7 ) );
327  CHECK( fp_equal_tol( expn(3,0.79), 0.1463479, 1.e-7 ) );
328  CHECK( fp_equal_tol( expn(3,1.97), 0.0312817, 1.e-7 ) );
329 
330  CHECK( fp_equal( expn(4,0.), 1./3., 10 ) );
331  CHECK( fp_equal_tol( expn(4,0.15), 0.2677889, 1.e-7 ) );
332  CHECK( fp_equal_tol( expn(4,0.79), 0.1127433, 1.e-7 ) );
333  CHECK( fp_equal_tol( expn(4,1.97), 0.0259440, 1.e-7 ) );
334 
335  CHECK( fp_equal( expn(10,0.), 1./9., 10 ) );
336  CHECK( fp_equal_tol( expn(10,0.15), 0.0938786, 1.e-7 ) );
337  CHECK( fp_equal_tol( expn(10,0.79), 0.0459453, 1.e-7 ) );
338  CHECK( fp_equal_tol( expn(10,1.97), 0.0124964, 1.e-7 ) );
339  }
340 
341  TEST(TestE1)
342  {
343  CHECK( fp_equal_tol( e1(0.25), 0.25*0.9408157528 - log(0.25) - EULER, 2.e-7 ) );
344  CHECK( fp_equal_tol( e1(0.79), 0.316277004, 2.e-7 ) );
345  CHECK( fp_equal_tol( e1(1.97), 0.050976988, 1.4e-9 ) );
346 
347  CHECK( fp_equal_tol( e1_scaled(0.25), (0.25*0.9408157528 - log(0.25) - EULER)*exp(0.25), 2.6e-7 ) );
348  CHECK( fp_equal_tol( e1_scaled(0.79), 0.316277004*exp(0.79), 4.4e-7 ) );
349  CHECK( fp_equal_tol( e1_scaled(1.97), 0.050976988*exp(1.97), 1.e-8 ) );
350  }
351 
352  TEST(TestE2)
353  {
354  CHECK( fp_equal_tol( e2(0.25), 0.8643037 + log(0.25)/4., 1.e-7 ) );
355  CHECK( fp_equal_tol( e2(0.79), 0.2039860, 1.e-7 ) );
356  CHECK( fp_equal_tol( e2(1.97), 0.0390322, 1.e-7 ) );
357  }
358 
359  TEST(TestErf)
360  {
361  /* constants calculated with xmaxima 5.22.1 */
362  CHECK( fp_equal_tol( erf(0.), 0., 1.e-22 ) );
363  // erf(x) loses some precision around 1.e-10, but should still be plenty good...
364  CHECK( fp_equal_tol( erf(1.e-10), 1.1283791671081724525e-10, 3.e-21 ) );
365  CHECK( fp_equal_tol( erf(1.e-5), 1.1283791670579000307e-5, 1.e-17 ) );
366  CHECK( fp_equal_tol( erf(0.1), 1.1246291601828489259e-1, 1.e-13 ) );
367  CHECK( fp_equal_tol( erf(0.5), 5.2049987781304653768e-1, 1.e-12 ) );
368  CHECK( fp_equal_tol( erf(1.), 8.4270079294971486934e-1, 1.e-12 ) );
369  CHECK( fp_equal_tol( erf(2.), 9.9532226501895273416e-1, 1.e-12 ) );
370  CHECK( fp_equal_tol( erf(10.), 1.0, 1.e-12 ) );
371  CHECK( fp_equal( erf(-1.), -erf(1.) ) );
372  }
373 
374  TEST(TestErfc)
375  {
376  /* constants calculated with xmaxima 5.22.1 */
377  CHECK( fp_equal_tol( erfc(-1.), 1.8427007929497148693e0, 1.e-12 ) );
378  CHECK( fp_equal_tol( erfc(0.), 1., 1.e-12 ) );
379  CHECK( fp_equal_tol( erfc(1.e-5), 9.99988716208329421e-1, 1.e-12 ) );
380  CHECK( fp_equal_tol( erfc(0.1), 8.8753708398171510741e-1, 1.e-12 ) );
381  CHECK( fp_equal_tol( erfc(0.5), 4.7950012218695346232e-1, 1.e-12 ) );
382  CHECK( fp_equal_tol( erfc(1.), 1.5729920705028513066e-1, 1.e-13 ) );
383  CHECK( fp_equal_tol( erfc(2.), 4.6777349810472658396e-3, 1.e-15 ) );
384  CHECK( fp_equal_tol( erfc(10.), 2.088487583762544757e-45, 1.e-57 ) );
385  }
386 
387  TEST(TestErfce)
388  {
389  /* constants taken from Finn G.D. & Mugglestone D., 1965, MNRAS 129, 221 */
390  /* this is the voigt function H(a,0) normalized according to 9-44 of Mihalas */
391  CHECK( fp_equal_tol( erfce(0.), 1., 1.e-6 ) );
392  CHECK( fp_equal_tol( erfce(0.01), 9.88815e-1, 1.e-6 ) );
393  CHECK( fp_equal_tol( erfce(0.02), 9.77826e-1, 1.e-6 ) );
394  CHECK( fp_equal_tol( erfce(0.05), 9.45990e-1, 1.e-6 ) );
395  CHECK( fp_equal_tol( erfce(0.10), 8.96456e-1, 1.e-6 ) );
396  CHECK( fp_equal_tol( erfce(0.20), 8.09019e-1, 1.e-6 ) );
397  CHECK( fp_equal_tol( erfce(0.55), 5.90927e-1, 1.e-6 ) );
398  CHECK( fp_equal_tol( erfce(1.00), 4.27583e-1, 1.e-6 ) );
399  }
400 
401  TEST(TestVoigtH0)
402  {
403  /* constants taken from Finn G.D. & Mugglestone D., 1965, MNRAS 129, 221 */
404  /* this is the voigt function H(a,0) normalized according to 9-44 of Mihalas */
405  CHECK( fp_equal_tol( VoigtH0(0.), 1., 1.e-6 ) );
406  CHECK( fp_equal_tol( VoigtH0(0.01), 9.88815e-1, 1.e-6 ) );
407  CHECK( fp_equal_tol( VoigtH0(0.02), 9.77826e-1, 1.e-6 ) );
408  CHECK( fp_equal_tol( VoigtH0(0.05), 9.45990e-1, 1.e-6 ) );
409  CHECK( fp_equal_tol( VoigtH0(0.10), 8.96456e-1, 1.e-6 ) );
410  CHECK( fp_equal_tol( VoigtH0(0.20), 8.09019e-1, 1.e-6 ) );
411  CHECK( fp_equal_tol( VoigtH0(0.55), 5.90927e-1, 1.e-6 ) );
412  CHECK( fp_equal_tol( VoigtH0(1.00), 4.27583e-1, 1.e-6 ) );
413  }
414 
415  TEST(TestVoigtU0)
416  {
417  /* constants taken from Finn G.D. & Mugglestone D., 1965, MNRAS 129, 221 */
418  /* this is the voigt function U(a,0) normalized according to 9-45 of Mihalas */
419  CHECK( fp_equal_tol( VoigtU0(0.), 1./SQRTPI, 1.e-6 ) );
420  CHECK( fp_equal_tol( VoigtU0(0.01), 9.88815e-1/SQRTPI, 1.e-6 ) );
421  CHECK( fp_equal_tol( VoigtU0(0.02), 9.77826e-1/SQRTPI, 1.e-6 ) );
422  CHECK( fp_equal_tol( VoigtU0(0.05), 9.45990e-1/SQRTPI, 1.e-6 ) );
423  CHECK( fp_equal_tol( VoigtU0(0.10), 8.96456e-1/SQRTPI, 1.e-6 ) );
424  CHECK( fp_equal_tol( VoigtU0(0.20), 8.09019e-1/SQRTPI, 1.e-6 ) );
425  CHECK( fp_equal_tol( VoigtU0(0.55), 5.90927e-1/SQRTPI, 1.e-6 ) );
426  CHECK( fp_equal_tol( VoigtU0(1.00), 4.27583e-1/SQRTPI, 1.e-6 ) );
427  }
428 
429  TEST(TestGegen)
430  {
431  /* constants calculated with xmaxima 5.30.0 */
432  CHECK( fp_equal_tol( gegenbauer(0,1,1.), 1., 1.e-12 ) );
433  CHECK( fp_equal_tol( gegenbauer(0,1,0.5), 1., 1.e-12 ) );
434  CHECK( fp_equal_tol( gegenbauer(1,1,0.5), 1., 1.e-12 ) );
435  CHECK( fp_equal_tol( gegenbauer(2,2,0.55), 1.63, 1.e-12 ) );
436  CHECK( fp_equal_tol( gegenbauer(2,0.2,0.5), -0.08, 1.e-12 ) );
437  CHECK( fp_equal_tol( gegenbauer(2,3,0.553), 4.339416, 1.e-12 ) );
438  CHECK( fp_equal_tol( gegenbauer(3,3,0.553), 0.25699016, 1.e-12 ) );
439  CHECK( fp_equal_tol( gegenbauer(3,3,-0.553), -0.25699016, 1.e-12 ) );
440  // test case from Levrie & Piessens, http://nalag.cs.kuleuven.be/papers/paul/report74/report74.ps.gz, corrected using xmaxima, xmaxima quoted error is 3e-13
441  CHECK( fp_equal_tol( gegenbauer(40,0.5,0.999), 0.333404687422222, 1.e-9 ) );
442  // test case from http://nalag.cs.kuleuven.be/papers/paul/report74/report74.ps.gz, corrected using xmaxima, xmaxima quoted error is 3e-13
443  CHECK( fp_equal_tol( gegenbauer(40,0.5,0.999), 0.333404687422222, 1.e-9 ) );
444  {
445  UltraGen u(1,1.);
446  CHECK( fp_equal_tol( u.val(), 1., 1.e-12 ) );
447  }
448  {
449  UltraGen u(1,0.5);
450  CHECK( fp_equal_tol( u.val(), 1., 1.e-12 ) );
451  }
452  {
453  UltraGen u(2,0.5);
454  u.step();
455  CHECK( fp_equal_tol( u.val(), 1., 1.e-12 ) );
456  }
457  {
458  UltraGen u(4,0.55);
459  u.step(); u.step();
460  CHECK( fp_equal_tol( u.val(), 1.63, 1.e-12 ) );
461  }
462  {
463  UltraGen u(5,0.553);
464  u.step(); u.step();
465  CHECK( fp_equal_tol( u.val(), 4.339416, 1.e-12 ) );
466  }
467  {
468  UltraGen u(6,0.553);
469  u.step(); u.step(); u.step();
470  CHECK( fp_equal_tol( u.val(), 0.25699016, 1.e-12 ) );
471  }
472  CHECK( fp_equal_tol( gegenbauer(3,3,-0.553), -0.25699016, 1.e-12 ) );
473 
474  }
475 
476  TEST(TestSixj)
477  {
478  long a=1,b=2,c=3,s=a+b+c;
479  double cc = ipow(-1,s)/sqrt(double((2*b+1)*(2*c+1)));
480  CHECK( fp_equal_tol( sjs(2*a,2*b,2*c,0,2*c,2*b), cc, 1.e-15 ) );
481 
482  // test failure modes
483  // this set of j's fails only Triangle2( j1, j2, j3 )
484  CHECK( sjs( 0, 6, 2, 4, 2, 2 ) == 0. );
485  // this set of j's fails only Triangle2( j4, j5, j3 )
486  CHECK( sjs( 2, 2, 0, 0, 4, 2 ) == 0. );
487  // this set of j's fails only Triangle2( j1, j5, j6 )
488  CHECK( sjs( 2, 4, 6, 6, 2, 6 ) == 0. );
489  // this set of j's fails only Triangle2( j2, j4, j6 )
490  CHECK( sjs( 2, 6, 4, 0, 4, 2 ) == 0. );
491 
492  // test sjs() for small arguments by comparison with SixJFull()
493  for( int j1=0; j1 < 5; ++j1 )
494  for( int j2=0; j2 < 5; ++j2 )
495  for( int j3=0; j3 < 5; ++j3 )
496  for( int j4=0; j4 < 5; ++j4 )
497  for( int j5=0; j5 < 5; ++j5 )
498  for( int j6=0; j6 < 5; ++j6 )
499  {
500  double s1 = sjs( j1, j2, j3, j4, j5, j6 );
501  if( !Triangle2( j1, j2, j3 ) || !Triangle2( j4, j5, j3 ) ||
502  !Triangle2( j1, j5, j6 ) || !Triangle2( j2, j4, j6 ) )
503  {
504  CHECK( s1 == 0. );
505  }
506  else
507  {
508  double s2 = SixJFull( j1, j2, j3, j4, j5, j6 );
509  CHECK( fp_equal_tol( s1, s2, 1.e-15 ) );
510  }
511  }
512 
513  // test a few handpicked 6j symbols with large arguments
514  // results were taken from
515  // http://www-stone.ch.cam.ac.uk/cgi-bin/wigner.cgi
516  double s1 = sjs( 42, 38, 34, 46, 36, 32 );
517  double s2 = -146218147./1778945129108.*sqrt(20735./57.);
518  CHECK( fp_equal_tol( s1, s2, 1.e-15 ) );
519  s1 = sjs( 46, 44, 42, 40, 38, 36 );
520  s2 = -16185716849919./4093621176253.*sqrt(201./135605470.);
521  CHECK( fp_equal_tol( s1, s2, 1.e-15 ) );
522  s1 = sjs( 43, 39, 40, 12, 48, 49 );
523  s2 = 1111./11600540.*sqrt(2982639./595.);
524  CHECK( fp_equal_tol( s1, s2, 1.e-15 ) );
525 
526  // test sjs() for random arguments by comparison with SixJFull()
527  srand(245694237);
528  for( int i=0; i < 1000; ++i )
529  {
530  int j1, j2, j3, j4, j5, j6;
531  const int JMAX = 110;
532  do
533  {
534  j1 = rand()%JMAX;
535  j2 = rand()%JMAX;
536  j3 = rand()%JMAX;
537  j4 = rand()%JMAX;
538  j5 = rand()%JMAX;
539  j6 = rand()%JMAX;
540  }
541  while( !Triangle2( j1, j2, j3 ) || !Triangle2( j4, j5, j3 ) ||
542  !Triangle2( j1, j5, j6 ) || !Triangle2( j2, j4, j6 ) );
543 
544  s1 = sjs( j1, j2, j3, j4, j5, j6 );
545  s2 = SixJFull( j1, j2, j3, j4, j5, j6 );
546  // for very large arguments, cancellation errors start growing
547  // so we need to increase the tolerance to 1e-12 here...
548  CHECK( fp_equal_tol( s1, s2, 1.e-12 ) );
549  }
550  }
551  TEST(TestRecSixj)
552  {
553  const int NPT=256;
554  long a=1,b=2,c=3,s=a+b+c;
555  double cc = pow(-1,s)/sqrt((2*b+1)*(2*c+1));
556  double sixcof[NPT];
557  long ier;
558  double l1min, l1max, lmatch;
559  rec6j(sixcof, b, c, 0, c, b, &l1min, &l1max, &lmatch, NPT, &ier);
560  CHECK( fp_equal_tol( sixcof[a-(int)l1min], cc, 1.e-15 ) );
561 
562  // test failure modes
563  // this set of j's fails only Triangle2( j1, j2, j3 )
564  rec6j(sixcof, 3, 1, 2, 1, 1, &l1min, &l1max, &lmatch, NPT, &ier);
565  CHECK( 0 < l1min );
566  // this set of j's fails only Triangle2( j4, j5, j3 )
567  rec6j(sixcof, 1, 0, 0, 2, 1, &l1min, &l1max, &lmatch, NPT, &ier);
568  CHECK( ier == -1 );
569  // this set of j's fails only Triangle2( j1, j5, j6 )
570  rec6j(sixcof, 1, 3, 3, 1, 3, &l1min, &l1max, &lmatch, NPT, &ier);
571  CHECK( 1 < l1min );
572  // this set of j's fails only Triangle2( j2, j4, j6 )
573  rec6j(sixcof, 3, 2, 0, 2, 1, &l1min, &l1max, &lmatch, NPT, &ier);
574  CHECK( ier == -1 );
575 
576  // test rec6j() for small arguments by comparison with SixJFull()
577  for( int j2=0; j2 < 5; ++j2 )
578  for( int j3=0; j3 < 5; ++j3 )
579  for( int j4=0; j4 < 5; ++j4 )
580  for( int j5=0; j5 < 5; ++j5 )
581  for( int j6=0; j6 < 5; ++j6 )
582  {
583  rec6j(sixcof, 0.5*j2, 0.5*j3, 0.5*j4, 0.5*j5, 0.5*j6, &l1min, &l1max, &lmatch, NPT, &ier);
584  for( int j1=0; j1 < 5; ++j1 )
585  {
586  double pos = (0.5*j1-l1min);
587  if( !Triangle2( j1, j2, j3 ) || !Triangle2( j4, j5, j3 ) ||
588  !Triangle2( j1, j5, j6 ) || !Triangle2( j2, j4, j6 ) )
589  {
590  CHECK( ier == -1 || pos != int(pos) ||
591  0.5*j1 < l1min || 0.5*j1 > l1max );
592  }
593  else
594  {
595  CHECK( ier >= 0 );
596  double s1 = sixcof[int(pos)];
597  double s2 = SixJFull( j1, j2, j3, j4, j5, j6 );
598  CHECK( fp_equal_tol( s1, s2, 1.e-15 ) );
599  }
600  }
601  }
602 
603  // test a few handpicked 6j symbols with large arguments
604  // results were taken from
605  // http://www-stone.ch.cam.ac.uk/cgi-bin/wigner.cgi
606  rec6j(sixcof, 19, 17, 23, 18, 16, &l1min, &l1max, &lmatch, NPT, &ier);
607  double s1 = sixcof[int(21-l1min)];
608  double s2 = -146218147./1778945129108.*sqrt(20735./57.);
609  CHECK( fp_equal_tol( s1, s2, 1.e-15 ) );
610  rec6j(sixcof, 22, 21, 20, 19, 18, &l1min, &l1max, &lmatch, NPT, &ier);
611  s1 = sixcof[int(23-l1min)];
612  s2 = -16185716849919./4093621176253.*sqrt(201./135605470.);
613  CHECK( fp_equal_tol( s1, s2, 1.e-15 ) );
614  rec6j(sixcof, 19.5, 20, 6, 24, 24.5, &l1min, &l1max, &lmatch, NPT, &ier);
615  s1 = sixcof[int(21.5-l1min)];
616  s2 = 1111./11600540.*sqrt(2982639./595.);
617  CHECK( fp_equal_tol( s1, s2, 1.e-15 ) );
618 
619  // test rec6j() for random arguments by comparison with SixJFull()
620  srand(245694237);
621  for( int i=0; i < 1000; ++i )
622  {
623  int j1, j2, j3, j4, j5, j6;
624  const int JMAX = 110;
625  do
626  {
627  j1 = rand()%JMAX;
628  j2 = rand()%JMAX;
629  j3 = rand()%JMAX;
630  j4 = rand()%JMAX;
631  j5 = rand()%JMAX;
632  j6 = rand()%JMAX;
633  }
634  while( !Triangle2( j1, j2, j3 ) || !Triangle2( j4, j5, j3 ) ||
635  !Triangle2( j1, j5, j6 ) || !Triangle2( j2, j4, j6 ) );
636 
637  rec6j(sixcof, 0.5*j2, 0.5*j3, 0.5*j4, 0.5*j5, 0.5*j6,
638  &l1min, &l1max, &lmatch, NPT, &ier);
639  s1 = sixcof[int(0.5*j1-l1min)];
640  s2 = SixJFull( j1, j2, j3, j4, j5, j6 );
641  // for very large arguments, cancellation errors start growing
642  // for SixJFull, so we need to increase the tolerance to 1e-12 here...
643  CHECK( fp_equal_tol( s1, s2, 1.e-12 ) );
644  }
645 
646  // test rec6j for problematic arguments
647  // comparison values were calculated using a 512-bit version of sjs()
648  rec6j(sixcof, 85., 84., 89.5, 89.5, 89.5, &l1min, &l1max, &lmatch, NPT, &ier);
649  s1 = sixcof[int(86.-l1min)];
650  CHECK( fp_equal_tol( s1, 3.3988213869175631e-08, 4.e-18 ) );
651  // check endpoints of recursion
652  CHECK( fp_equal( l1min, 1. ) );
653  CHECK( fp_equal( l1max, 169. ) );
654  s1 = sixcof[int(1.-l1min)];
655  CHECK( fp_equal_tol( s1, 0.0035632864755963242, 4.e-18 ) );
656  s1 = sixcof[int(169.-l1min)];
657  CHECK( fp_equal_tol( s1, 1.785014472974611e-17, 2.e-31 ) );
658  }
659 
660  TEST(TestVoigtH)
661  {
662  // check that the Voigt profile returned by VoigtH() is properly normalized
663  const int NP = 200;
664  realnum v[NP], a, y[NP];
665  // for a > 0.1, VoigtH() calls humlik(), for smaller values it
666  // calls FastVoigtH().
667  //
668  // humlik() is set up for a relative precision of 1e-4 over its
669  // entire range, but looks to be more precise in practice (at
670  // least at v=0)
671  // FastVoigtH() is set up for a rel. precision of 2.5e-3 over its
672  // entire range, but should be better than 1e-4 for a < 0.0235
673  a = realnum(0.0002);
674  for( int i=0; i < 9; ++i )
675  {
676  // test both humlik() and FastVoigtH()
677  for( int i=0; i < NP; ++i )
678  v[i] = realnum(i)*max(a,1.f)/realnum(5.);
679  VoigtH( a, v, y, NP );
680  realnum integral = realnum(0.);
681  // We need the integral from -infinity to +infinity, which is simply
682  // two times the integral from 0 to +infinity. Hence we omit the
683  // division by 2 in the trapezoidal rule
684  for( int i=1; i < NP; ++i )
685  integral += (y[i]+y[i-1])*(v[i]-v[i-1]);
686  // add in the integral over the Lorentz wings assuming H(a,v) = c/v^2
687  integral += realnum(2.)*v[NP-1]*y[NP-1];
688  // VoigtH() calculates H(a,v), so integral should be sqrt(pi)
689  CHECK( fp_equal_tol( integral, realnum(SQRTPI), realnum(1.e-4) ) );
690  // also check the central value...
691  realnum h0 = realnum(VoigtH0(a));
692  CHECK( fp_equal_tol( y[0], h0, realnum(1.e-4)*h0 ) );
693  a *= realnum(10.);
694  }
695  // now do some spot checks
696  // constants taken from
697  // >>refer Zaghloul M.R. & Ali A.N. 2011, ACM Transactions on Mathematical Software, 38, 15
698  // available from http://arxiv.org/abs/1106.0151
699  v[0] = realnum(5.76);
700  VoigtH( realnum(1.e-20), v, y, 1 );
701  // this constant comes from page 21 ("Present algorithm").
702  CHECK( fp_equal_tol( y[0], realnum(3.900779639194698e-015), realnum(4.e-19) ) );
703  v[0] = realnum(6.3e-2);
704  v[1] = realnum(6.3e+0);
705  v[2] = realnum(6.3e+2);
706  VoigtH( realnum(1.e-20), v, y, 3 );
707  // these constants come from Table 2 of the same paper ("Present algorithm").
708  CHECK( fp_equal_tol( y[0], realnum(9.960388660702479e-001), realnum(1.e-4) ) );
709  CHECK( fp_equal_tol( y[1], realnum(5.792460778844116e-018), realnum(6.e-22) ) );
710  CHECK( fp_equal_tol( y[2], realnum(1.421495882582395e-026), realnum(1.4e-30) ) );
711  VoigtH( realnum(1.e-14), v, y, 3 );
712  CHECK( fp_equal_tol( y[0], realnum(9.960388660702366e-001), realnum(1.e-4) ) );
713  CHECK( fp_equal_tol( y[1], realnum(1.536857621303163e-016), realnum(1.5e-20) ) );
714  CHECK( fp_equal_tol( y[2], realnum(1.421495882582395e-020), realnum(1.4e-24) ) );
715  VoigtH( realnum(1.e-12), v, y, 3 );
716  CHECK( fp_equal_tol( y[0], realnum(9.960388660691284e-001), realnum(1.e-4) ) );
717  CHECK( fp_equal_tol( y[1], realnum(1.479513723737753e-014), realnum(1.5e-18) ) );
718  CHECK( fp_equal_tol( y[2], realnum(1.421495882582395e-018), realnum(1.4e-22) ) );
719  VoigtH( realnum(1.e-10), v, y, 3 );
720  CHECK( fp_equal_tol( y[0], realnum(9.960388659583033e-001), realnum(1.e-4) ) );
721  CHECK( fp_equal_tol( y[1], realnum(1.478940284762099e-012), realnum(1.5e-16) ) );
722  CHECK( fp_equal_tol( y[2], realnum(1.421495882582395e-016), realnum(1.4e-20) ) );
723  VoigtH( realnum(1.e-6), v, y, 3 );
724  CHECK( fp_equal_tol( y[0], realnum(9.960377466254801e-001), realnum(1.e-4) ) );
725  CHECK( fp_equal_tol( y[1], realnum(1.478934493028404e-008), realnum(1.5e-12) ) );
726  CHECK( fp_equal_tol( y[2], realnum(1.421495882582395e-012), realnum(1.4e-16) ) );
727  VoigtH( realnum(1.e-2), v, y, 3 );
728  CHECK( fp_equal_tol( y[0], realnum(9.849424862549039e-001), realnum(1.e-4) ) );
729  CHECK( fp_equal_tol( y[1], realnum(1.478930389133934e-004), realnum(1.5e-8) ) );
730  CHECK( fp_equal_tol( y[2], realnum(1.421495882224242e-008), realnum(1.4e-12) ) );
731  VoigtH( realnum(1.e+1), v, y, 3 );
732  CHECK( fp_equal_tol( y[0], realnum(5.613881832823886e-002), realnum(6.e-6) ) );
733  CHECK( fp_equal_tol( y[1], realnum(4.040671157393835e-002), realnum(4.e-6) ) );
734  CHECK( fp_equal_tol( y[2], realnum(1.421137820009847e-005), realnum(1.4e-9) ) );
735  VoigtH( realnum(1.2e+1), v, y, 3 );
736  CHECK( fp_equal_tol( y[0], realnum(4.685295149211637e-002), realnum(5.e-6) ) );
737  CHECK( fp_equal_tol( y[1], realnum(3.684277239564798e-002), realnum(4.e-6) ) );
738  CHECK( fp_equal_tol( y[2], realnum(1.705176395541707e-005), realnum(1.7e-9) ) );
739  VoigtH( realnum(1.5e+1), v, y, 3 );
740  CHECK( fp_equal_tol( y[0], realnum(3.752895161491574e-002), realnum(4.e-6) ) );
741  CHECK( fp_equal_tol( y[1], realnum(3.194834330452605e-002), realnum(3.e-6) ) );
742  CHECK( fp_equal_tol( y[2], realnum(2.131035743074598e-005), realnum(2.e-9) ) );
743  VoigtH( realnum(2.e+2), v, y, 3 );
744  CHECK( fp_equal_tol( y[0], realnum(2.820912377324508e-003), realnum(3.e-7) ) );
745  CHECK( fp_equal_tol( y[1], realnum(2.818116555672206e-003), realnum(3.e-7) ) );
746  CHECK( fp_equal_tol( y[2], realnum(2.582702147491469e-004), realnum(3.e-8) ) );
747  VoigtH( realnum(1.e+5), v, y, 3 );
748  CHECK( fp_equal_tol( y[0], realnum(5.641895835193228e-006), realnum(6.e-10) ) );
749  CHECK( fp_equal_tol( y[1], realnum(5.641895812802746e-006), realnum(6.e-10) ) );
750  CHECK( fp_equal_tol( y[2], realnum(5.641671917237128e-006), realnum(6.e-10) ) );
751  v[0] = realnum(1.e+0);
752  VoigtH( realnum(1.e-20), v, y, 1 );
753  CHECK( fp_equal_tol( y[0], realnum(3.678794411714423e-001), realnum(4.e-5) ) );
754  v[0] = realnum(5.5e+0);
755  VoigtH( realnum(1.e-14), v, y, 1 );
756  CHECK( fp_equal_tol( y[0], realnum(7.307386729528773e-014), realnum(7.e-18) ) );
757  v[0] = realnum(3.9e+4);
758  VoigtH( realnum(1.e+0), v, y, 1 );
759  CHECK( fp_equal_tol( y[0], realnum(3.709333226385423e-010), realnum(4.e-14) ) );
760  v[0] = realnum(1.e+0);
761  VoigtH( realnum(2.8e+4), v, y, 1 );
762  CHECK( fp_equal_tol( y[0], realnum(2.014962794529686e-005), realnum(2.e-9) ) );
763  }
764 
765  TEST(TestVoigtU)
766  {
767  // check that the Voigt profile returned by VoigtU() is properly normalized
768  const int NP = 200;
769  realnum v[NP], a, y[NP];
770  // for a > 0.1, VoigtU() calls humlik(), for smaller values it
771  // calls FastVoigtH() and divides by sqrt(pi).
772  //
773  // humlik() is set up for a relative precision of 1e-4 over its
774  // entire range, but looks to be more precise in practice (at
775  // least at v=0)
776  // FastVoigtH() is set up for a rel. precision of 2.5e-3 over its
777  // entire range, but should be better than 1e-4 for a < 0.0235
778  a = realnum(0.0002);
779  for( int i=0; i < 9; ++i )
780  {
781  // test both humlik() and FastVoigtH()
782  for( int i=0; i < NP; ++i )
783  v[i] = realnum(i)*max(a,1.f)/realnum(5.);
784  VoigtU( a, v, y, NP );
785  realnum integral = realnum(0.);
786  // We need the integral from -infinity to +infinity, which is simply
787  // two times the integral from 0 to +infinity. Hence we omit the
788  // division by 2 in the trapezoidal rule
789  for( int i=1; i < NP; ++i )
790  integral += (y[i]+y[i-1])*(v[i]-v[i-1]);
791  // add in the integral over the Lorentz wings assuming U(a,v) = c/v^2
792  integral += realnum(2.)*v[NP-1]*y[NP-1];
793  // VoigtU() calculates U(a,v), so integral should be 1
794  CHECK( fp_equal_tol( integral, realnum(1.), realnum(1.e-4) ) );
795  // also check the central value...
796  CHECK( fp_equal_tol( y[0], realnum(VoigtU0(a)), realnum(1.e-4) ) );
797  a *= realnum(10.);
798  }
799  }
800 
801  TEST(TestMD5string)
802  {
803  string test;
804  // md5sum of an empty file...
805  CHECK( MD5string( test ) == "d41d8cd98f00b204e9800998ecf8427e" );
806  CHECK( MD5string( test ).length() == NMD5 );
807  // check if padding is done correctly
808  // an extra block of padding needs to be added when length%64 == 56
809  test = "hhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhh";
810  CHECK( test.length() == 55 );
811  CHECK( MD5string( test ) == "426ec4ac35ad38d125f6efb39da03098" );
812  test = "hhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhh";
813  CHECK( test.length() == 56 );
814  CHECK( MD5string( test ) == "d03607b2c89adc0c4abf5a0f1d9e40c9" );
815  test = "hhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhh";
816  CHECK( test.length() == 57 );
817  CHECK( MD5string( test ) == "bac1b47748411cb6eee0cae3befb8377" );
818  string test64 = "hhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhh";
819  test = test64 + "hhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhh";
820  CHECK( test.length() == 64+55 );
821  CHECK( MD5string( test ) == "10d49aad1fc69976376fbe7c8c5ed118" );
822  test = test64 + "hhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhh";
823  CHECK( test.length() == 64+56 );
824  CHECK( MD5string( test ) == "61ec7da14576f3b585038c6d72cd5bd5" );
825  test = test64 + "hhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhh";
826  CHECK( test.length() == 64+57 );
827  CHECK( MD5string( test ) == "f17a0475a26d0930e2a35bb320c10e0d" );
828  // check that leading zeros are printed correctly
829  test = "ghijklmn";
830  CHECK( MD5string( test ) == "0256b9cea63bc1f97b8c5aea92c24a98" );
831  }
832 
833 }
bessel_i0_i1
void bessel_i0_i1(double x, double *i0val, double *i1val)
Definition: thirdparty.cpp:2976
bessel_k0_scaled
double bessel_k0_scaled(double x)
Definition: thirdparty.cpp:2723
bessel_i1_scaled
double bessel_i1_scaled(double x)
Definition: thirdparty.cpp:2955
e1
double e1(double x)
Definition: thirdparty.cpp:3258
lfactorial
double lfactorial(long n)
Definition: thirdparty.cpp:756
bessel_i1
double bessel_i1(double x)
Definition: thirdparty.cpp:2934
SixJFull
double SixJFull(long j1, long j2, long j3, long j4, long j5, long j6)
Definition: thirdparty.cpp:3558
Triangle2
bool Triangle2(long J1, long J2, long J3)
Definition: thirdparty.h:497
fp_equal_tol
bool fp_equal_tol(sys_float x, sys_float y, sys_float tol)
Definition: cddefines.h:904
sjs
double sjs(long j1, long j2, long j3, long l1, long l2, long l3)
Definition: thirdparty.cpp:3494
bessel_i0_i1_scaled
void bessel_i0_i1_scaled(double x, double *i0val, double *i1val)
Definition: thirdparty.cpp:3003
expn
double expn(int n, double x)
Definition: thirdparty.cpp:1756
cdgamma
complex< double > cdgamma(complex< double > x)
Definition: thirdparty.cpp:789
bessel_k0_k1
void bessel_k0_k1(double x, double *k0val, double *k1val)
Definition: thirdparty.cpp:2792
fp_equal
bool fp_equal(sys_float x, sys_float y, int n=3)
Definition: cddefines.h:854
VoigtU0
double VoigtU0(double a)
Definition: thirdparty.h:445
bessel_k1
double bessel_k1(double x)
Definition: thirdparty.cpp:2746
erfc
double erfc(double)
VoigtH0
double VoigtH0(double a)
Definition: thirdparty.h:439
bessel_i0_scaled
double bessel_i0_scaled(double x)
Definition: thirdparty.cpp:2913
bessel_i0
double bessel_i0(double x)
Definition: thirdparty.cpp:2892
bessel_jn
double bessel_jn(int n, double x)
Definition: thirdparty.cpp:1399
bessel_k0_k1_scaled
void bessel_k0_k1_scaled(double x, double *k0val, double *k1val)
Definition: thirdparty.cpp:2820
realnum
float realnum
Definition: cddefines.h:124
bessel_y0
double bessel_y0(double x)
Definition: thirdparty.cpp:1103
max
long max(int a, long b)
Definition: cddefines.h:817
bessel_j0
double bessel_j0(double x)
Definition: thirdparty.cpp:1065
bessel_j1
double bessel_j1(double x)
Definition: thirdparty.cpp:1296
igamc
double igamc(double a, double x)
Definition: thirdparty.cpp:2311
igam
double igam(double a, double x)
Definition: thirdparty.cpp:2356
ipow
long int ipow(long, long)
Definition: service.cpp:658
factorial
double factorial(long n)
Definition: thirdparty.cpp:356
bessel_yn
double bessel_yn(int n, double x)
Definition: thirdparty.cpp:1534
bessel_k1_scaled
double bessel_k1_scaled(double x)
Definition: thirdparty.cpp:2769
erfce
double erfce(double x)
Definition: thirdparty.cpp:2089
VoigtU
void VoigtU(realnum a, const realnum v[], realnum y[], int n)
Definition: thirdparty.h:431
bessel_y1
double bessel_y1(double x)
Definition: thirdparty.cpp:1323
NMD5
static const unsigned int NMD5
Definition: thirdparty.h:451
erf
double erf(double)
e1_scaled
double e1_scaled(double x)
Definition: thirdparty.cpp:3296
rec6j
void rec6j(double *sixcof, double l2, double l3, double l4, double l5, double l6, double *l1min, double *l1max, double *lmatch, long ndim, long *ier)
Definition: thirdparty.cpp:3608
ellpk
double ellpk(double x)
Definition: thirdparty.cpp:1676
igamc_scaled
double igamc_scaled(double a, double x)
Definition: thirdparty.cpp:2329
e2
double e2(double x)
Definition: thirdparty.cpp:3328
gegenbauer
double gegenbauer(long n, double al, double x)
Definition: thirdparty.cpp:3464
bessel_k0
double bessel_k0(double x)
Definition: thirdparty.cpp:2700
MD5string
string MD5string(const string &str)
Definition: thirdparty.cpp:4667
VoigtH
void VoigtH(realnum a, const realnum v[], realnum y[], int n)
Definition: thirdparty.h:418
Generated by   doxygen 1.8.5