cloudy  trunk
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
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 
12 bool linfit(
13  long n,
14  double x[], /* x[n] */
15  double y[], /* y[n] */
16  double &a,
17  double &siga,
18  double &b,
19  double &sigb
20 );
21 
24 static const int NPRE_FACTORIAL = 171;
25 
27 double factorial(long n);
28 
30 double lfactorial(long n);
31 
32 complex<double> cdgamma(complex<double> x);
33 
34 double bessel_j0(double x);
35 double bessel_y0(double x);
36 double bessel_j1(double x);
37 double bessel_y1(double x);
38 double bessel_jn(int n, double x);
39 double bessel_yn(int n, double x);
40 
41 double bessel_k0(double x);
42 double bessel_k0_scaled(double x);
43 double bessel_k1(double x);
44 double bessel_k1_scaled(double x);
45 double bessel_i0(double x);
46 double bessel_i0_scaled(double x);
47 double bessel_i1(double x);
48 double bessel_i1_scaled(double x);
49 
50 double ellpk(double x);
51 
56 double expn(int n, double x);
57 
58 /* initializes mt[N] with a seed */
59 void init_genrand(unsigned long s);
60 
61 /* initialize by an array with array-length */
62 /* init_key is the array for initializing keys */
63 /* key_length is its length */
64 /* slight change for C++, 2004/2/26 */
65 void init_by_array(unsigned long init_key[], int key_length);
66 
67 /* generates a random number on [0,0xffffffff]-interval */
68 unsigned long genrand_int32(void);
69 
70 /* generates a random number on [0,0x7fffffff]-interval */
71 long genrand_int31(void);
72 
73 /* These real versions are due to Isaku Wada, 2002/01/09 added */
74 /* generates a random number on [0,1]-real-interval */
75 double genrand_real1(void);
76 
77 /* generates a random number on [0,1)-real-interval */
78 double genrand_real2(void);
79 
80 /* generates a random number on (0,1)-real-interval */
81 double genrand_real3(void);
82 
83 /* generates a random number on [0,1) with 53-bit resolution*/
84 double genrand_res53(void);
85 
86 /*============================================================================*/
87 
88 /* these are the routines in thirdparty_interpolate.cpp */
89 
90 void spline_cubic_set( long n, const double t[], const double y[], double ypp[],
91  int ibcbeg, double ybcbeg, int ibcend, double ybcend );
92 void spline_cubic_val( long n, const double t[], double tval, const double y[], const double ypp[],
93  double *yval, double *ypval, double *yppval );
94 
107 inline void spline(const double x[],
108  const double y[],
109  long int n,
110  double yp1,
111  double ypn,
112  double y2a[])
113 {
114  int ibcbeg = ( yp1 > 0.99999e30 ) ? 2 : 1;
115  double ybcbeg = ( yp1 > 0.99999e30 ) ? 0. : yp1;
116  int ibcend = ( ypn > 0.99999e30 ) ? 2 : 1;
117  double ybcend = ( ypn > 0.99999e30 ) ? 0. : ypn;
118  spline_cubic_set( n, x, y, y2a, ibcbeg, ybcbeg, ibcend, ybcend );
119  return;
120 }
121 
130 inline void splint(const double xa[],
131  const double ya[],
132  const double y2a[],
133  long int n,
134  double x,
135  double *y)
136 {
137  spline_cubic_val( n, xa, x, ya, y2a, y, NULL, NULL );
138  return;
139 }
140 
141 inline void spldrv(const double xa[],
142  const double ya[],
143  const double y2a[],
144  long int n,
145  double x,
146  double *y)
147 {
148  spline_cubic_val( n, xa, x, ya, y2a, NULL, y, NULL );
149  return;
150 }
151 
152 /* wrapper routine for splint that checks whether x-value is within bounds
153  * if the x-value is out of bounds, a flag will be raised and the function
154  * will be evaluated at the nearest boundary */
155 /* >>chng 03 jan 15, added splint_safe, PvH */
156 inline void splint_safe(const double xa[],
157  const double ya[],
158  const double y2a[],
159  long int n,
160  double x,
161  double *y,
162  bool *lgOutOfBounds)
163 {
164  double xsafe;
165 
166  const double lo_bound = MIN2(xa[0],xa[n-1]);
167  const double hi_bound = MAX2(xa[0],xa[n-1]);
168  const double SAFETY = MAX2(hi_bound-lo_bound,1.)*10.*DBL_EPSILON;
169 
170  DEBUG_ENTRY( "splint_safe()" );
171 
172  if( x < lo_bound-SAFETY )
173  {
174  xsafe = lo_bound;
175  *lgOutOfBounds = true;
176  }
177  else if( x > hi_bound+SAFETY )
178  {
179  xsafe = hi_bound;
180  *lgOutOfBounds = true;
181  }
182  else
183  {
184  xsafe = x;
185  *lgOutOfBounds = false;
186  }
187 
188  splint(xa,ya,y2a,n,xsafe,y);
189  return;
190 }
191 
192 /* wrapper routine for spldrv that checks whether x-value is within bounds
193  * if the x-value is out of bounds, a flag will be raised and the function
194  * will be evaluated at the nearest boundary */
195 /* >>chng 03 jan 15, added spldrv_safe, PvH */
196 inline void spldrv_safe(const double xa[],
197  const double ya[],
198  const double y2a[],
199  long int n,
200  double x,
201  double *y,
202  bool *lgOutOfBounds)
203 {
204  double xsafe;
205 
206  const double lo_bound = MIN2(xa[0],xa[n-1]);
207  const double hi_bound = MAX2(xa[0],xa[n-1]);
208  const double SAFETY = MAX2(fabs(lo_bound),fabs(hi_bound))*10.*DBL_EPSILON;
209 
210  DEBUG_ENTRY( "spldrv_safe()" );
211 
212  if( x < lo_bound-SAFETY )
213  {
214  xsafe = lo_bound;
215  *lgOutOfBounds = true;
216  }
217  else if( x > hi_bound+SAFETY )
218  {
219  xsafe = hi_bound;
220  *lgOutOfBounds = true;
221  }
222  else
223  {
224  xsafe = x;
225  *lgOutOfBounds = false;
226  }
227 
228  spldrv(xa,ya,y2a,n,xsafe,y);
229  return;
230 }
231 
240 double lagrange(const double x[], /* x[n] */
241  const double y[], /* y[n] */
242  long n,
243  double xval);
244 
252 double linint(const double x[], /* x[n] */
253  const double y[], /* y[n] */
254  long n,
255  double xval);
256 
259 template<class T>
260 inline long hunt_bisect(const T x[], /* x[n] */
261  long n,
262  T xval)
263 {
264  /* do bisection hunt */
265  long ilo = 0, ihi = n-1;
266  while( ihi-ilo > 1 )
267  {
268  long imid = (ilo+ihi)/2;
269  if( xval < x[imid] )
270  ihi = imid;
271  else
272  ilo = imid;
273  }
274  return ilo;
275 }
276 
279 template<class T>
280 inline long hunt_bisect_reverse(const T x[], /* x[n] */
281  long n,
282  T xval)
283 {
284  /* do bisection hunt */
285  long ilo = 0, ihi = n-1;
286  while( ihi-ilo > 1 )
287  {
288  long imid = (ilo+ihi)/2;
289  if( xval <= x[imid] )
290  ilo = imid;
291  else
292  ihi = imid;
293  }
294  return ilo;
295 }
296 
297 /*============================================================================*/
298 
299 /* these are the routines in thirdparty_lapack.cpp */
300 
301 /* there are wrappers for lapack linear algebra routines.
302  * there are two versions of the lapack routines - a fortran
303  * version that I converted to C with forc to use if nothing else is available
304  * (included in the Cloudy distribution),
305  * and an option to link into an external lapack library that may be optimized
306  * for your machine. By default the tralated C routines will be used.
307  * To use your machine's lapack library instead, define the macro
308  * LAPACK and link into your library. This is usually done with a command line
309  * switch "-DLAPACK" on the compiler command, and the linker option "-llapack"
310  */
311 
320 void getrf_wrapper(long M, long N, double *A, long lda, int32 *ipiv, int32 *info);
321 
333 void getrs_wrapper(char trans, long N, long nrhs, double *A, long lda, int32 *ipiv, double *B, long ldb, int32 *info);
334 
335 #endif /* _THIRDPARTY_H_ */

Generated for cloudy by doxygen 1.8.4