cloudy trunk
Loading...
Searching...
No Matches
thirdparty.h
Go to the documentation of this file.
1/* This file contains routines (perhaps in modified form) by third parties.
2 * Use and distribution of these works are determined by their respective copyrights. */
3
4#ifndef THIRDPARTY_H_
5#define THIRDPARTY_H_
6
7
8/*============================================================================*/
9
10/* these are the routines in thirdparty.cpp */
11
12bool linfit(
13 long n,
14 const double x[], /* x[n] */
15 const double y[], /* y[n] */
16 double &a,
17 double &siga,
18 double &b,
19 double &sigb
20);
21
24static const int NPRE_FACTORIAL = 171;
25
27double factorial(long n);
28
30double lfactorial(long n);
31
32complex<double> cdgamma(complex<double> x);
33
34double bessel_j0(double x);
35double bessel_y0(double x);
36double bessel_j1(double x);
37double bessel_y1(double x);
38double bessel_jn(int n, double x);
39double bessel_yn(int n, double x);
40
41double bessel_k0(double x);
42double bessel_k0_scaled(double x);
43double bessel_k1(double x);
44double bessel_k1_scaled(double x);
45double bessel_i0(double x);
46double bessel_i0_scaled(double x);
47double bessel_i1(double x);
48double bessel_i1_scaled(double x);
49
50double ellpk(double x);
51
56double expn(int n, double x);
57
59#ifndef HAVE_ERF
60double erf(double);
61#endif
62#ifndef HAVE_ERFC
63double erfc(double);
64#endif
66double erfce(double);
67
68/* initializes mt[N] with a seed */
69void init_genrand(unsigned long s);
70
71/* initialize by an array with array-length */
72/* init_key is the array for initializing keys */
73/* key_length is its length */
74/* slight change for C++, 2004/2/26 */
75void init_by_array(unsigned long init_key[], int key_length);
76
77/* generates a random number on [0,0xffffffff]-interval */
78unsigned long genrand_int32(void);
79
80/* generates a random number on [0,0x7fffffff]-interval */
81long genrand_int31(void);
82
83/* These real versions are due to Isaku Wada, 2002/01/09 added */
84/* generates a random number on [0,1]-real-interval */
85double genrand_real1(void);
86
87/* generates a random number on [0,1)-real-interval */
88double genrand_real2(void);
89
90/* generates a random number on (0,1)-real-interval */
91double genrand_real3(void);
92
93/* generates a random number on [0,1) with 53-bit resolution*/
94double genrand_res53(void);
95
96/*============================================================================*/
97
98/* these are the routines in thirdparty_interpolate.cpp */
99
100void spline_cubic_set( long n, const double t[], const double y[], double ypp[],
101 int ibcbeg, double ybcbeg, int ibcend, double ybcend );
102void spline_cubic_val( long n, const double t[], double tval, const double y[], const double ypp[],
103 double *yval, double *ypval, double *yppval );
104
117inline void spline(const double x[],
118 const double y[],
119 long int n,
120 double yp1,
121 double ypn,
122 double y2a[])
123{
124 int ibcbeg = ( yp1 > 0.99999e30 ) ? 2 : 1;
125 double ybcbeg = ( yp1 > 0.99999e30 ) ? 0. : yp1;
126 int ibcend = ( ypn > 0.99999e30 ) ? 2 : 1;
127 double ybcend = ( ypn > 0.99999e30 ) ? 0. : ypn;
128 spline_cubic_set( n, x, y, y2a, ibcbeg, ybcbeg, ibcend, ybcend );
129 return;
130}
131
140inline void splint(const double xa[],
141 const double ya[],
142 const double y2a[],
143 long int n,
144 double x,
145 double *y)
146{
147 spline_cubic_val( n, xa, x, ya, y2a, y, NULL, NULL );
148 return;
149}
150
151inline void spldrv(const double xa[],
152 const double ya[],
153 const double y2a[],
154 long int n,
155 double x,
156 double *y)
157{
158 spline_cubic_val( n, xa, x, ya, y2a, NULL, y, NULL );
159 return;
160}
161
162/* wrapper routine for splint that checks whether x-value is within bounds
163 * if the x-value is out of bounds, a flag will be raised and the function
164 * will be evaluated at the nearest boundary */
165/* >>chng 03 jan 15, added splint_safe, PvH */
166inline void splint_safe(const double xa[],
167 const double ya[],
168 const double y2a[],
169 long int n,
170 double x,
171 double *y,
172 bool *lgOutOfBounds)
173{
174 double xsafe;
175
176 const double lo_bound = MIN2(xa[0],xa[n-1]);
177 const double hi_bound = MAX2(xa[0],xa[n-1]);
178 const double SAFETY = MAX2(hi_bound-lo_bound,1.)*10.*DBL_EPSILON;
179
180 DEBUG_ENTRY( "splint_safe()" );
181
182 if( x < lo_bound-SAFETY )
183 {
184 xsafe = lo_bound;
185 *lgOutOfBounds = true;
186 }
187 else if( x > hi_bound+SAFETY )
188 {
189 xsafe = hi_bound;
190 *lgOutOfBounds = true;
191 }
192 else
193 {
194 xsafe = x;
195 *lgOutOfBounds = false;
196 }
197
198 splint(xa,ya,y2a,n,xsafe,y);
199 return;
200}
201
202/* wrapper routine for spldrv that checks whether x-value is within bounds
203 * if the x-value is out of bounds, a flag will be raised and the function
204 * will be evaluated at the nearest boundary */
205/* >>chng 03 jan 15, added spldrv_safe, PvH */
206inline void spldrv_safe(const double xa[],
207 const double ya[],
208 const double y2a[],
209 long int n,
210 double x,
211 double *y,
212 bool *lgOutOfBounds)
213{
214 double xsafe;
215
216 const double lo_bound = MIN2(xa[0],xa[n-1]);
217 const double hi_bound = MAX2(xa[0],xa[n-1]);
218 const double SAFETY = MAX2(fabs(lo_bound),fabs(hi_bound))*10.*DBL_EPSILON;
219
220 DEBUG_ENTRY( "spldrv_safe()" );
221
222 if( x < lo_bound-SAFETY )
223 {
224 xsafe = lo_bound;
225 *lgOutOfBounds = true;
226 }
227 else if( x > hi_bound+SAFETY )
228 {
229 xsafe = hi_bound;
230 *lgOutOfBounds = true;
231 }
232 else
233 {
234 xsafe = x;
235 *lgOutOfBounds = false;
236 }
237
238 spldrv(xa,ya,y2a,n,xsafe,y);
239 return;
240}
241
250double lagrange(const double x[], /* x[n] */
251 const double y[], /* y[n] */
252 long n,
253 double xval);
254
262double linint(const double x[], /* x[n] */
263 const double y[], /* y[n] */
264 long n,
265 double xval);
266
269template<class T>
270inline long hunt_bisect(const T x[], /* x[n] */
271 long n,
272 T xval)
273{
274 /* do bisection hunt */
275 long ilo = 0, ihi = n-1;
276 while( ihi-ilo > 1 )
277 {
278 long imid = (ilo+ihi)/2;
279 if( xval < x[imid] )
280 ihi = imid;
281 else
282 ilo = imid;
283 }
284 return ilo;
285}
286
289template<class T>
290inline long hunt_bisect_reverse(const T x[], /* x[n] */
291 long n,
292 T xval)
293{
294 /* do bisection hunt */
295 long ilo = 0, ihi = n-1;
296 while( ihi-ilo > 1 )
297 {
298 long imid = (ilo+ihi)/2;
299 if( xval <= x[imid] )
300 ilo = imid;
301 else
302 ihi = imid;
303 }
304 return ilo;
305}
306
307/*============================================================================*/
308
309/* these are the routines in thirdparty_lapack.cpp */
310
311/* there are wrappers for lapack linear algebra routines.
312 * there are two versions of the lapack routines - a fortran
313 * version that I converted to C with forc to use if nothing else is available
314 * (included in the Cloudy distribution),
315 * and an option to link into an external lapack library that may be optimized
316 * for your machine. By default the tralated C routines will be used.
317 * To use your machine's lapack library instead, define the macro
318 * LAPACK and link into your library. This is usually done with a command line
319 * switch "-DLAPACK" on the compiler command, and the linker option "-llapack"
320 */
321
330void getrf_wrapper(long M, long N, double *A, long lda, int32 *ipiv, int32 *info);
331
343void getrs_wrapper(char trans, long N, long nrhs, double *A, long lda, int32 *ipiv, double *B, long ldb, int32 *info);
344
345void humlik(int n, const realnum x[], realnum y, realnum k[]);
346
348
349// calculates y[i] = H(a,v[i]) as defined in Eq. 9-44 of Mihalas
350inline void VoigtH(realnum a, const realnum v[], realnum y[], int n)
351{
352 if( a <= 0.1f )
353 {
354 for( int i=0; i < n; ++i )
355 y[i] = FastVoigtH(a,v[i]);
356 }
357 else
358 {
359 humlik( n, v, a, y );
360 }
361}
362
363// calculates y[i] = U(a,v[i]) as defined in Eq. 9-45 of Mihalas
364inline void VoigtU(realnum a, const realnum v[], realnum y[], int n)
365{
366 VoigtH( a, v, y, n );
367 for( int i=0; i < n; ++i )
368 y[i] /= realnum(SQRTPI);
369}
370
371// VoigtH0(a) returns the value H(a,0) following Eq. 9-44 of Mihalas
372inline double VoigtH0(double a)
373{
374 return erfce(a);
375}
376
377// VoigtU0(a) returns the value U(a,0) following Eq. 9-45 of Mihalas
378inline double VoigtU0(double a)
379{
380 return erfce(a)/SQRTPI;
381}
382
384static const unsigned int NMD5 = 32;
385
387string MD5file(const char* fnam, access_scheme scheme=AS_DATA_ONLY);
389string MD5datafile(const char* fnam, access_scheme scheme=AS_DATA_ONLY);
391string MD5string(const string& str);
392
393#endif /* THIRDPARTY_H_ */
#define MIN2
Definition cddefines.h:761
float realnum
Definition cddefines.h:103
#define MAX2
Definition cddefines.h:782
#define DEBUG_ENTRY(funcname)
Definition cddefines.h:684
access_scheme
Definition cpu.h:207
@ AS_DATA_ONLY
Definition cpu.h:207
static const double SAFETY
UNUSED const double SQRTPI
Definition physconst.h:44
static const int M
static const int N
string MD5datafile(const char *fnam, access_scheme scheme=AS_DATA_ONLY)
double erfce(double)
double bessel_j1(double x)
double bessel_k1_scaled(double x)
double factorial(long n)
unsigned long genrand_int32(void)
double bessel_k0(double x)
double ellpk(double x)
double bessel_i1_scaled(double x)
double bessel_yn(int n, double x)
double erfc(double)
void spldrv(const double xa[], const double ya[], const double y2a[], long int n, double x, double *y)
Definition thirdparty.h:151
double bessel_y0(double x)
double bessel_k0_scaled(double x)
static const unsigned int NMD5
Definition thirdparty.h:384
double genrand_real2(void)
void humlik(int n, const realnum x[], realnum y, realnum k[])
void getrs_wrapper(char trans, long N, long nrhs, double *A, long lda, int32 *ipiv, double *B, long ldb, int32 *info)
double linint(const double x[], const double y[], long n, double xval)
void init_genrand(unsigned long s)
double bessel_jn(int n, double x)
double VoigtH0(double a)
Definition thirdparty.h:372
long genrand_int31(void)
long hunt_bisect_reverse(const T x[], long n, T xval)
Definition thirdparty.h:290
string MD5string(const string &str)
string MD5file(const char *fnam, access_scheme scheme=AS_DATA_ONLY)
double lfactorial(long n)
double lagrange(const double x[], const double y[], long n, double xval)
void spline_cubic_set(long n, const double t[], const double y[], double ypp[], int ibcbeg, double ybcbeg, int ibcend, double ybcend)
void spldrv_safe(const double xa[], const double ya[], const double y2a[], long int n, double x, double *y, bool *lgOutOfBounds)
Definition thirdparty.h:206
complex< double > cdgamma(complex< double > x)
void VoigtH(realnum a, const realnum v[], realnum y[], int n)
Definition thirdparty.h:350
long hunt_bisect(const T x[], long n, T xval)
Definition thirdparty.h:270
static const int NPRE_FACTORIAL
Definition thirdparty.h:24
void splint_safe(const double xa[], const double ya[], const double y2a[], long int n, double x, double *y, bool *lgOutOfBounds)
Definition thirdparty.h:166
void spline_cubic_val(long n, const double t[], double tval, const double y[], const double ypp[], double *yval, double *ypval, double *yppval)
double bessel_i1(double x)
double bessel_i0_scaled(double x)
realnum FastVoigtH(realnum a, realnum v)
double bessel_j0(double x)
double genrand_real3(void)
double VoigtU0(double a)
Definition thirdparty.h:378
double bessel_i0(double x)
void init_by_array(unsigned long init_key[], int key_length)
double genrand_res53(void)
double genrand_real1(void)
void spline(const double x[], const double y[], long int n, double yp1, double ypn, double y2a[])
Definition thirdparty.h:117
double bessel_k1(double x)
void getrf_wrapper(long M, long N, double *A, long lda, int32 *ipiv, int32 *info)
void splint(const double xa[], const double ya[], const double y2a[], long int n, double x, double *y)
Definition thirdparty.h:140
double expn(int n, double x)
double erf(double)
bool linfit(long n, const double x[], const double y[], double &a, double &siga, double &b, double &sigb)
double bessel_y1(double x)
void VoigtU(realnum a, const realnum v[], realnum y[], int n)
Definition thirdparty.h:364