00001
00002
00003
00004
00005
00006
00007
00008 bool linfit(
00009 long n,
00010 double x[],
00011 double y[],
00012 double &a,
00013 double &siga,
00014 double &b,
00015 double &sigb
00016 );
00017
00020 static const int NPRE_FACTORIAL = 171;
00021
00023 double factorial(long n);
00024
00026 double lfactorial(long n);
00027
00028 complex<double> cdgamma(complex<double> x);
00029
00030 double bessel_j0(double x);
00031 double bessel_y0(double x);
00032 double bessel_j1(double x);
00033 double bessel_y1(double x);
00034 double bessel_jn(int n, double x);
00035 double bessel_yn(int n, double x);
00036
00037 double bessel_k0(double x);
00038 double bessel_k0_scaled(double x);
00039 double bessel_k1(double x);
00040 double bessel_k1_scaled(double x);
00041 double bessel_i0(double x);
00042 double bessel_i0_scaled(double x);
00043 double bessel_i1(double x);
00044 double bessel_i1_scaled(double x);
00045
00046 double ellpk(double x);
00047
00052 double expn(int n, double x);
00053
00054
00055
00056
00057
00058 void spline_cubic_set( long n, const double t[], const double y[], double ypp[],
00059 int ibcbeg, double ybcbeg, int ibcend, double ybcend );
00060 void spline_cubic_val( long n, const double t[], double tval, const double y[], const double ypp[],
00061 double *yval, double *ypval, double *yppval );
00062
00063
00064
00065
00066
00067
00068
00069
00070 inline void spline(const double x[],
00071 const double y[],
00072 long int n,
00073 double yp1,
00074 double ypn,
00075 double y2a[])
00076 {
00077 int ibcbeg = ( yp1 > 0.99999e30 ) ? 2 : 1;
00078 double ybcbeg = ( yp1 > 0.99999e30 ) ? 0. : yp1;
00079 int ibcend = ( ypn > 0.99999e30 ) ? 2 : 1;
00080 double ybcend = ( ypn > 0.99999e30 ) ? 0. : ypn;
00081 spline_cubic_set( n, x, y, y2a, ibcbeg, ybcbeg, ibcend, ybcend );
00082 return;
00083 }
00084
00085 inline void splint(const double xa[],
00086 const double ya[],
00087 const double y2a[],
00088 long int n,
00089 double x,
00090 double *y)
00091 {
00092 spline_cubic_val( n, xa, x, ya, y2a, y, NULL, NULL );
00093 return;
00094 }
00095
00096 inline void spldrv(const double xa[],
00097 const double ya[],
00098 const double y2a[],
00099 long int n,
00100 double x,
00101 double *y)
00102 {
00103 spline_cubic_val( n, xa, x, ya, y2a, NULL, y, NULL );
00104 return;
00105 }
00106
00107
00108
00109
00110
00111 inline void splint_safe(const double xa[],
00112 const double ya[],
00113 const double y2a[],
00114 long int n,
00115 double x,
00116 double *y,
00117 bool *lgOutOfBounds)
00118 {
00119 double xsafe;
00120
00121 const double lo_bound = MIN2(xa[0],xa[n-1]);
00122 const double hi_bound = MAX2(xa[0],xa[n-1]);
00123 const double SAFETY = MAX2(hi_bound-lo_bound,1.)*10.*DBL_EPSILON;
00124
00125 DEBUG_ENTRY( "splint_safe()" );
00126
00127 if( x < lo_bound-SAFETY )
00128 {
00129 xsafe = lo_bound;
00130 *lgOutOfBounds = true;
00131 }
00132 else if( x > hi_bound+SAFETY )
00133 {
00134 xsafe = hi_bound;
00135 *lgOutOfBounds = true;
00136 }
00137 else
00138 {
00139 xsafe = x;
00140 *lgOutOfBounds = false;
00141 }
00142
00143 splint(xa,ya,y2a,n,xsafe,y);
00144
00145 DEBUG_EXIT( "splint_safe()" );
00146 return;
00147 }
00148
00149
00150
00151
00152
00153 inline void spldrv_safe(const double xa[],
00154 const double ya[],
00155 const double y2a[],
00156 long int n,
00157 double x,
00158 double *y,
00159 bool *lgOutOfBounds)
00160 {
00161 double xsafe;
00162
00163 const double lo_bound = MIN2(xa[0],xa[n-1]);
00164 const double hi_bound = MAX2(xa[0],xa[n-1]);
00165 const double SAFETY = MAX2(fabs(lo_bound),fabs(hi_bound))*10.*DBL_EPSILON;
00166
00167 DEBUG_ENTRY( "spldrv_safe()" );
00168
00169 if( x < lo_bound-SAFETY )
00170 {
00171 xsafe = lo_bound;
00172 *lgOutOfBounds = true;
00173 }
00174 else if( x > hi_bound+SAFETY )
00175 {
00176 xsafe = hi_bound;
00177 *lgOutOfBounds = true;
00178 }
00179 else
00180 {
00181 xsafe = x;
00182 *lgOutOfBounds = false;
00183 }
00184
00185 spldrv(xa,ya,y2a,n,xsafe,y);
00186
00187 DEBUG_EXIT( "spldrv_safe()" );
00188 return;
00189 }
00190
00199 double lagrange(const double x[],
00200 const double y[],
00201 long n,
00202 double xval);
00203
00211 double linint(const double x[],
00212 const double y[],
00213 long n,
00214 double xval);
00215
00218 inline long hunt_bisect(const double x[],
00219 long n,
00220 double xval)
00221 {
00222
00223 long ilo = 0, ihi = n-1;
00224 while( ihi-ilo > 1 )
00225 {
00226 long imid = (ilo+ihi)/2;
00227 if( xval < x[imid] )
00228 ihi = imid;
00229 else
00230 ilo = imid;
00231 }
00232 return ilo;
00233 }
00234
00237 inline long hunt_bisect_reverse(const double x[],
00238 long n,
00239 double xval)
00240 {
00241
00242 long ilo = 0, ihi = n-1;
00243 while( ihi-ilo > 1 )
00244 {
00245 long imid = (ilo+ihi)/2;
00246 if( xval <= x[imid] )
00247 ilo = imid;
00248 else
00249 ihi = imid;
00250 }
00251 return ilo;
00252 }
00253
00254
00255
00256
00257
00258
00259
00260
00261
00262
00263
00264
00265
00266
00267
00268
00277 void getrf_wrapper(long M, long N, double *A, long lda, int32 *ipiv, int32 *info);
00278
00290 void getrs_wrapper(char trans, long N, long nrhs, double *A, long lda, int32 *ipiv, double *B, long ldb, int32 *info);