cloudy trunk
|
00001 /* This file contains routines (perhaps in modified form) by third parties. 00002 * Use and distribution of these works are determined by their respective copyrights. */ 00003 00004 #ifndef _THIRDPARTY_H_ 00005 #define _THIRDPARTY_H_ 00006 00007 00008 /*============================================================================*/ 00009 00010 /* these are the routines in thirdparty.cpp */ 00011 00012 bool linfit( 00013 long n, 00014 double x[], /* x[n] */ 00015 double y[], /* y[n] */ 00016 double &a, 00017 double &siga, 00018 double &b, 00019 double &sigb 00020 ); 00021 00024 static const int NPRE_FACTORIAL = 171; 00025 00027 double factorial(long n); 00028 00030 double lfactorial(long n); 00031 00032 complex<double> cdgamma(complex<double> x); 00033 00034 double bessel_j0(double x); 00035 double bessel_y0(double x); 00036 double bessel_j1(double x); 00037 double bessel_y1(double x); 00038 double bessel_jn(int n, double x); 00039 double bessel_yn(int n, double x); 00040 00041 double bessel_k0(double x); 00042 double bessel_k0_scaled(double x); 00043 double bessel_k1(double x); 00044 double bessel_k1_scaled(double x); 00045 double bessel_i0(double x); 00046 double bessel_i0_scaled(double x); 00047 double bessel_i1(double x); 00048 double bessel_i1_scaled(double x); 00049 00050 double ellpk(double x); 00051 00056 double expn(int n, double x); 00057 00058 /* initializes mt[N] with a seed */ 00059 void init_genrand(unsigned long s); 00060 00061 /* initialize by an array with array-length */ 00062 /* init_key is the array for initializing keys */ 00063 /* key_length is its length */ 00064 /* slight change for C++, 2004/2/26 */ 00065 void init_by_array(unsigned long init_key[], int key_length); 00066 00067 /* generates a random number on [0,0xffffffff]-interval */ 00068 unsigned long genrand_int32(void); 00069 00070 /* generates a random number on [0,0x7fffffff]-interval */ 00071 long genrand_int31(void); 00072 00073 /* These real versions are due to Isaku Wada, 2002/01/09 added */ 00074 /* generates a random number on [0,1]-real-interval */ 00075 double genrand_real1(void); 00076 00077 /* generates a random number on [0,1)-real-interval */ 00078 double genrand_real2(void); 00079 00080 /* generates a random number on (0,1)-real-interval */ 00081 double genrand_real3(void); 00082 00083 /* generates a random number on [0,1) with 53-bit resolution*/ 00084 double genrand_res53(void); 00085 00086 /*============================================================================*/ 00087 00088 /* these are the routines in thirdparty_interpolate.cpp */ 00089 00090 void spline_cubic_set( long n, const double t[], const double y[], double ypp[], 00091 int ibcbeg, double ybcbeg, int ibcend, double ybcend ); 00092 void spline_cubic_val( long n, const double t[], double tval, const double y[], const double ypp[], 00093 double *yval, double *ypval, double *yppval ); 00094 00107 inline void spline(const double x[], 00108 const double y[], 00109 long int n, 00110 double yp1, 00111 double ypn, 00112 double y2a[]) 00113 { 00114 int ibcbeg = ( yp1 > 0.99999e30 ) ? 2 : 1; 00115 double ybcbeg = ( yp1 > 0.99999e30 ) ? 0. : yp1; 00116 int ibcend = ( ypn > 0.99999e30 ) ? 2 : 1; 00117 double ybcend = ( ypn > 0.99999e30 ) ? 0. : ypn; 00118 spline_cubic_set( n, x, y, y2a, ibcbeg, ybcbeg, ibcend, ybcend ); 00119 return; 00120 } 00121 00130 inline void splint(const double xa[], 00131 const double ya[], 00132 const double y2a[], 00133 long int n, 00134 double x, 00135 double *y) 00136 { 00137 spline_cubic_val( n, xa, x, ya, y2a, y, NULL, NULL ); 00138 return; 00139 } 00140 00141 inline void spldrv(const double xa[], 00142 const double ya[], 00143 const double y2a[], 00144 long int n, 00145 double x, 00146 double *y) 00147 { 00148 spline_cubic_val( n, xa, x, ya, y2a, NULL, y, NULL ); 00149 return; 00150 } 00151 00152 /* wrapper routine for splint that checks whether x-value is within bounds 00153 * if the x-value is out of bounds, a flag will be raised and the function 00154 * will be evaluated at the nearest boundary */ 00155 /* >>chng 03 jan 15, added splint_safe, PvH */ 00156 inline void splint_safe(const double xa[], 00157 const double ya[], 00158 const double y2a[], 00159 long int n, 00160 double x, 00161 double *y, 00162 bool *lgOutOfBounds) 00163 { 00164 double xsafe; 00165 00166 const double lo_bound = MIN2(xa[0],xa[n-1]); 00167 const double hi_bound = MAX2(xa[0],xa[n-1]); 00168 const double SAFETY = MAX2(hi_bound-lo_bound,1.)*10.*DBL_EPSILON; 00169 00170 DEBUG_ENTRY( "splint_safe()" ); 00171 00172 if( x < lo_bound-SAFETY ) 00173 { 00174 xsafe = lo_bound; 00175 *lgOutOfBounds = true; 00176 } 00177 else if( x > hi_bound+SAFETY ) 00178 { 00179 xsafe = hi_bound; 00180 *lgOutOfBounds = true; 00181 } 00182 else 00183 { 00184 xsafe = x; 00185 *lgOutOfBounds = false; 00186 } 00187 00188 splint(xa,ya,y2a,n,xsafe,y); 00189 return; 00190 } 00191 00192 /* wrapper routine for spldrv that checks whether x-value is within bounds 00193 * if the x-value is out of bounds, a flag will be raised and the function 00194 * will be evaluated at the nearest boundary */ 00195 /* >>chng 03 jan 15, added spldrv_safe, PvH */ 00196 inline void spldrv_safe(const double xa[], 00197 const double ya[], 00198 const double y2a[], 00199 long int n, 00200 double x, 00201 double *y, 00202 bool *lgOutOfBounds) 00203 { 00204 double xsafe; 00205 00206 const double lo_bound = MIN2(xa[0],xa[n-1]); 00207 const double hi_bound = MAX2(xa[0],xa[n-1]); 00208 const double SAFETY = MAX2(fabs(lo_bound),fabs(hi_bound))*10.*DBL_EPSILON; 00209 00210 DEBUG_ENTRY( "spldrv_safe()" ); 00211 00212 if( x < lo_bound-SAFETY ) 00213 { 00214 xsafe = lo_bound; 00215 *lgOutOfBounds = true; 00216 } 00217 else if( x > hi_bound+SAFETY ) 00218 { 00219 xsafe = hi_bound; 00220 *lgOutOfBounds = true; 00221 } 00222 else 00223 { 00224 xsafe = x; 00225 *lgOutOfBounds = false; 00226 } 00227 00228 spldrv(xa,ya,y2a,n,xsafe,y); 00229 return; 00230 } 00231 00240 double lagrange(const double x[], /* x[n] */ 00241 const double y[], /* y[n] */ 00242 long n, 00243 double xval); 00244 00252 double linint(const double x[], /* x[n] */ 00253 const double y[], /* y[n] */ 00254 long n, 00255 double xval); 00256 00259 template<class T> 00260 inline long hunt_bisect(const T x[], /* x[n] */ 00261 long n, 00262 T xval) 00263 { 00264 /* do bisection hunt */ 00265 long ilo = 0, ihi = n-1; 00266 while( ihi-ilo > 1 ) 00267 { 00268 long imid = (ilo+ihi)/2; 00269 if( xval < x[imid] ) 00270 ihi = imid; 00271 else 00272 ilo = imid; 00273 } 00274 return ilo; 00275 } 00276 00279 template<class T> 00280 inline long hunt_bisect_reverse(const T x[], /* x[n] */ 00281 long n, 00282 T xval) 00283 { 00284 /* do bisection hunt */ 00285 long ilo = 0, ihi = n-1; 00286 while( ihi-ilo > 1 ) 00287 { 00288 long imid = (ilo+ihi)/2; 00289 if( xval <= x[imid] ) 00290 ilo = imid; 00291 else 00292 ihi = imid; 00293 } 00294 return ilo; 00295 } 00296 00297 /*============================================================================*/ 00298 00299 /* these are the routines in thirdparty_lapack.cpp */ 00300 00301 /* there are wrappers for lapack linear algebra routines. 00302 * there are two versions of the lapack routines - a fortran 00303 * version that I converted to C with forc to use if nothing else is available 00304 * (included in the Cloudy distribution), 00305 * and an option to link into an external lapack library that may be optimized 00306 * for your machine. By default the tralated C routines will be used. 00307 * To use your machine's lapack library instead, define the macro 00308 * LAPACK and link into your library. This is usually done with a command line 00309 * switch "-DLAPACK" on the compiler command, and the linker option "-llapack" 00310 */ 00311 00320 void getrf_wrapper(long M, long N, double *A, long lda, int32 *ipiv, int32 *info); 00321 00333 void getrs_wrapper(char trans, long N, long nrhs, double *A, long lda, int32 *ipiv, double *B, long ldb, int32 *info); 00334 00335 #endif /* _THIRDPARTY_H_ */