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
00059 #ifndef HAVE_ERF
00060 double erf(double);
00061 #endif
00062 #ifndef HAVE_ERFC
00063 double erfc(double);
00064 #endif
00065
00066 double erfce(double);
00067
00068
00069 void init_genrand(unsigned long s);
00070
00071
00072
00073
00074
00075 void init_by_array(unsigned long init_key[], int key_length);
00076
00077
00078 unsigned long genrand_int32(void);
00079
00080
00081 long genrand_int31(void);
00082
00083
00084
00085 double genrand_real1(void);
00086
00087
00088 double genrand_real2(void);
00089
00090
00091 double genrand_real3(void);
00092
00093
00094 double genrand_res53(void);
00095
00096
00097
00098
00099
00100 void spline_cubic_set( long n, const double t[], const double y[], double ypp[],
00101 int ibcbeg, double ybcbeg, int ibcend, double ybcend );
00102 void spline_cubic_val( long n, const double t[], double tval, const double y[], const double ypp[],
00103 double *yval, double *ypval, double *yppval );
00104
00117 inline void spline(const double x[],
00118 const double y[],
00119 long int n,
00120 double yp1,
00121 double ypn,
00122 double y2a[])
00123 {
00124 int ibcbeg = ( yp1 > 0.99999e30 ) ? 2 : 1;
00125 double ybcbeg = ( yp1 > 0.99999e30 ) ? 0. : yp1;
00126 int ibcend = ( ypn > 0.99999e30 ) ? 2 : 1;
00127 double ybcend = ( ypn > 0.99999e30 ) ? 0. : ypn;
00128 spline_cubic_set( n, x, y, y2a, ibcbeg, ybcbeg, ibcend, ybcend );
00129 return;
00130 }
00131
00140 inline void splint(const double xa[],
00141 const double ya[],
00142 const double y2a[],
00143 long int n,
00144 double x,
00145 double *y)
00146 {
00147 spline_cubic_val( n, xa, x, ya, y2a, y, NULL, NULL );
00148 return;
00149 }
00150
00151 inline void spldrv(const double xa[],
00152 const double ya[],
00153 const double y2a[],
00154 long int n,
00155 double x,
00156 double *y)
00157 {
00158 spline_cubic_val( n, xa, x, ya, y2a, NULL, y, NULL );
00159 return;
00160 }
00161
00162
00163
00164
00165
00166 inline void splint_safe(const double xa[],
00167 const double ya[],
00168 const double y2a[],
00169 long int n,
00170 double x,
00171 double *y,
00172 bool *lgOutOfBounds)
00173 {
00174 double xsafe;
00175
00176 const double lo_bound = MIN2(xa[0],xa[n-1]);
00177 const double hi_bound = MAX2(xa[0],xa[n-1]);
00178 const double SAFETY = MAX2(hi_bound-lo_bound,1.)*10.*DBL_EPSILON;
00179
00180 DEBUG_ENTRY( "splint_safe()" );
00181
00182 if( x < lo_bound-SAFETY )
00183 {
00184 xsafe = lo_bound;
00185 *lgOutOfBounds = true;
00186 }
00187 else if( x > hi_bound+SAFETY )
00188 {
00189 xsafe = hi_bound;
00190 *lgOutOfBounds = true;
00191 }
00192 else
00193 {
00194 xsafe = x;
00195 *lgOutOfBounds = false;
00196 }
00197
00198 splint(xa,ya,y2a,n,xsafe,y);
00199 return;
00200 }
00201
00202
00203
00204
00205
00206 inline void spldrv_safe(const double xa[],
00207 const double ya[],
00208 const double y2a[],
00209 long int n,
00210 double x,
00211 double *y,
00212 bool *lgOutOfBounds)
00213 {
00214 double xsafe;
00215
00216 const double lo_bound = MIN2(xa[0],xa[n-1]);
00217 const double hi_bound = MAX2(xa[0],xa[n-1]);
00218 const double SAFETY = MAX2(fabs(lo_bound),fabs(hi_bound))*10.*DBL_EPSILON;
00219
00220 DEBUG_ENTRY( "spldrv_safe()" );
00221
00222 if( x < lo_bound-SAFETY )
00223 {
00224 xsafe = lo_bound;
00225 *lgOutOfBounds = true;
00226 }
00227 else if( x > hi_bound+SAFETY )
00228 {
00229 xsafe = hi_bound;
00230 *lgOutOfBounds = true;
00231 }
00232 else
00233 {
00234 xsafe = x;
00235 *lgOutOfBounds = false;
00236 }
00237
00238 spldrv(xa,ya,y2a,n,xsafe,y);
00239 return;
00240 }
00241
00250 double lagrange(const double x[],
00251 const double y[],
00252 long n,
00253 double xval);
00254
00262 double linint(const double x[],
00263 const double y[],
00264 long n,
00265 double xval);
00266
00269 template<class T>
00270 inline long hunt_bisect(const T x[],
00271 long n,
00272 T xval)
00273 {
00274
00275 long ilo = 0, ihi = n-1;
00276 while( ihi-ilo > 1 )
00277 {
00278 long imid = (ilo+ihi)/2;
00279 if( xval < x[imid] )
00280 ihi = imid;
00281 else
00282 ilo = imid;
00283 }
00284 return ilo;
00285 }
00286
00289 template<class T>
00290 inline long hunt_bisect_reverse(const T x[],
00291 long n,
00292 T xval)
00293 {
00294
00295 long ilo = 0, ihi = n-1;
00296 while( ihi-ilo > 1 )
00297 {
00298 long imid = (ilo+ihi)/2;
00299 if( xval <= x[imid] )
00300 ilo = imid;
00301 else
00302 ihi = imid;
00303 }
00304 return ilo;
00305 }
00306
00307
00308
00309
00310
00311
00312
00313
00314
00315
00316
00317
00318
00319
00320
00321
00330 void getrf_wrapper(long M, long N, double *A, long lda, int32 *ipiv, int32 *info);
00331
00343 void getrs_wrapper(char trans, long N, long nrhs, double *A, long lda, int32 *ipiv, double *B, long ldb, int32 *info);
00344
00345 void humlik(int n, const realnum x[], realnum y, realnum k[]);
00346
00347 realnum FastVoigtH(realnum a, realnum v);
00348
00349
00350 inline void VoigtH(realnum a, const realnum v[], realnum y[], int n)
00351 {
00352 if( a <= 0.1f )
00353 {
00354 for( int i=0; i < n; ++i )
00355 y[i] = FastVoigtH(a,v[i]);
00356 }
00357 else
00358 {
00359 humlik( n, v, a, y );
00360 }
00361 }
00362
00363
00364 inline void VoigtU(realnum a, const realnum v[], realnum y[], int n)
00365 {
00366 VoigtH( a, v, y, n );
00367 for( int i=0; i < n; ++i )
00368 y[i] /= realnum(SQRTPI);
00369 }
00370
00371
00372 inline double VoigtH0(double a)
00373 {
00374 return erfce(a);
00375 }
00376
00377
00378 inline double VoigtU0(double a)
00379 {
00380 return erfce(a)/SQRTPI;
00381 }
00382
00384 static const unsigned int NMD5 = 32;
00385
00387 string MD5file(const char* fnam, access_scheme scheme=AS_DATA_ONLY);
00389 string MD5datafile(const char* fnam, access_scheme scheme=AS_DATA_ONLY);
00391 string MD5string(const string& str);
00392
00393 #endif