00001
00002
00003
00004 #ifndef THIRDPARTY_H_
00005 #define THIRDPARTY_H_
00006
00007
00008
00009
00010
00011
00012 bool linfit(
00013 long n,
00014 const double x[],
00015 const double y[],
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
00059 void init_genrand(unsigned long s);
00060
00061
00062
00063
00064
00065 void init_by_array(unsigned long init_key[], int key_length);
00066
00067
00068 unsigned long genrand_int32(void);
00069
00070
00071 long genrand_int31(void);
00072
00073
00074
00075 double genrand_real1(void);
00076
00077
00078 double genrand_real2(void);
00079
00080
00081 double genrand_real3(void);
00082
00083
00084 double genrand_res53(void);
00085
00086
00087
00088
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
00153
00154
00155
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
00193
00194
00195
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[],
00241 const double y[],
00242 long n,
00243 double xval);
00244
00252 double linint(const double x[],
00253 const double y[],
00254 long n,
00255 double xval);
00256
00259 template<class T>
00260 inline long hunt_bisect(const T x[],
00261 long n,
00262 T xval)
00263 {
00264
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[],
00281 long n,
00282 T xval)
00283 {
00284
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
00300
00301
00302
00303
00304
00305
00306
00307
00308
00309
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 void humlik(int n, realnum x[], realnum y, realnum k[]);
00336
00338 static const unsigned int NMD5 = 32;
00339
00341 string MD5file(const char* fnam, access_scheme scheme=AS_DATA_ONLY);
00343 string MD5datafile(const char* fnam, access_scheme scheme=AS_DATA_ONLY);
00345 string MD5string(const string& str);
00346
00347 #endif