thirdparty.h

Go to the documentation of this file.

/* This file contains routines (perhaps in modified form) by third parties.
 * Use and distribution of these works are determined by their respective copyrights. */
 
#ifndef THIRDPARTY_H_
#define THIRDPARTY_H_
 
 
/*============================================================================*/
 
/* these are the routines in thirdparty.cpp */
 
bool linfit(
        long n,
        const double x[], /* x[n] */
        const double y[], /* y[n] */
        double &a,
        double &siga,
        double &b,
        double &sigb
);
 
static const int NPRE_FACTORIAL = 171;
 
double factorial(long n);
 
double lfactorial(long n);
 
complex<double> cdgamma(complex<double> x);
 
double bessel_j0(double x);
double bessel_y0(double x);
double bessel_j1(double x);
double bessel_y1(double x);
double bessel_jn(int n, double x);
double bessel_yn(int n, double x);
 
double bessel_k0(double x);
double bessel_k0_scaled(double x);
double bessel_k1(double x);
double bessel_k1_scaled(double x);
double bessel_i0(double x);
double bessel_i0_scaled(double x);
double bessel_i1(double x);
double bessel_i1_scaled(double x);
 
double ellpk(double x);
 
double expn(int n, double x);
 
#ifndef HAVE_ERF
double erf(double);
#endif
#ifndef HAVE_ERFC
double erfc(double);
#endif
 
double erfce(double);
 
/* initializes mt[N] with a seed */
void init_genrand(unsigned long s);
 
/* initialize by an array with array-length */
/* init_key is the array for initializing keys */
/* key_length is its length */
/* slight change for C++, 2004/2/26 */
void init_by_array(unsigned long init_key[], int key_length);
 
/* generates a random number on [0,0xffffffff]-interval */
unsigned long genrand_int32(void);
 
/* generates a random number on [0,0x7fffffff]-interval */
long genrand_int31(void);
 
/* These real versions are due to Isaku Wada, 2002/01/09 added */
/* generates a random number on [0,1]-real-interval */
double genrand_real1(void);
 
/* generates a random number on [0,1)-real-interval */
double genrand_real2(void);
 
/* generates a random number on (0,1)-real-interval */
double genrand_real3(void);
 
/* generates a random number on [0,1) with 53-bit resolution*/
double genrand_res53(void);
 
/*============================================================================*/
 
/* these are the routines in thirdparty_interpolate.cpp */
 
void spline_cubic_set( long n, const double t[], const double y[], double ypp[],
                       int ibcbeg, double ybcbeg, int ibcend, double ybcend );
void spline_cubic_val( long n, const double t[], double tval, const double y[], const double ypp[],
                       double *yval, double *ypval, double *yppval );
 
inline void spline(const double x[], 
                   const double y[], 
                   long int n, 
                   double yp1, 
                   double ypn, 
                   double y2a[])
{
        int ibcbeg = ( yp1 > 0.99999e30 ) ? 2 : 1;
        double ybcbeg = ( yp1 > 0.99999e30 ) ? 0. : yp1;
        int ibcend = ( ypn > 0.99999e30 ) ? 2 : 1;
        double ybcend = ( ypn > 0.99999e30 ) ? 0. : ypn;
        spline_cubic_set( n, x, y, y2a, ibcbeg, ybcbeg, ibcend, ybcend );
        return;
}
 
inline void splint(const double xa[], 
                   const double ya[], 
                   const double y2a[], 
                   long int n, 
                   double x, 
                   double *y)
{
        spline_cubic_val( n, xa, x, ya, y2a, y, NULL, NULL );
        return;
}
 
inline void spldrv(const double xa[], 
                   const double ya[], 
                   const double y2a[], 
                   long int n, 
                   double x, 
                   double *y)
{
        spline_cubic_val( n, xa, x, ya, y2a, NULL, y, NULL );
        return;
}
 
/* wrapper routine for splint that checks whether x-value is within bounds
 * if the x-value is out of bounds, a flag will be raised and the function
 * will be evaluated at the nearest boundary */
/* >>chng 03 jan 15, added splint_safe, PvH */
inline void splint_safe(const double xa[], 
                        const double ya[], 
                        const double y2a[], 
                        long int n, 
                        double x, 
                        double *y,
                        bool *lgOutOfBounds)
{
        double xsafe;
 
        const double lo_bound = MIN2(xa[0],xa[n-1]);
        const double hi_bound = MAX2(xa[0],xa[n-1]);
        const double SAFETY = MAX2(hi_bound-lo_bound,1.)*10.*DBL_EPSILON;
 
        DEBUG_ENTRY( "splint_safe()" );
 
        if( x < lo_bound-SAFETY )
        {
                xsafe = lo_bound;
                *lgOutOfBounds = true;
        }
        else if( x > hi_bound+SAFETY )
        {
                xsafe = hi_bound;
                *lgOutOfBounds = true;
        }
        else
        {
                xsafe = x;
                *lgOutOfBounds = false;
        }
 
        splint(xa,ya,y2a,n,xsafe,y);
        return;
}
 
/* wrapper routine for spldrv that checks whether x-value is within bounds
 * if the x-value is out of bounds, a flag will be raised and the function
 * will be evaluated at the nearest boundary */
/* >>chng 03 jan 15, added spldrv_safe, PvH */
inline void spldrv_safe(const double xa[],
                        const double ya[],
                        const double y2a[],
                        long int n,
                        double x,
                        double *y,
                        bool *lgOutOfBounds)
{
        double xsafe;
 
        const double lo_bound = MIN2(xa[0],xa[n-1]);
        const double hi_bound = MAX2(xa[0],xa[n-1]);
        const double SAFETY = MAX2(fabs(lo_bound),fabs(hi_bound))*10.*DBL_EPSILON;
 
        DEBUG_ENTRY( "spldrv_safe()" );
 
        if( x < lo_bound-SAFETY )
        {
                xsafe = lo_bound;
                *lgOutOfBounds = true;
        }
        else if( x > hi_bound+SAFETY )
        {
                xsafe = hi_bound;
                *lgOutOfBounds = true;
        }
        else
        {
                xsafe = x;
                *lgOutOfBounds = false;
        }
 
        spldrv(xa,ya,y2a,n,xsafe,y);
        return;
}
 
double lagrange(const double x[], /* x[n] */
                const double y[], /* y[n] */
                long n,
                double xval);
 
double linint(const double x[], /* x[n] */
              const double y[], /* y[n] */
              long n,
              double xval);
 
template<class T>
inline long hunt_bisect(const T x[], /* x[n] */
                        long n,
                        T xval)
{
        /* do bisection hunt */
        long ilo = 0, ihi = n-1;
        while( ihi-ilo > 1 )
        {
                long imid = (ilo+ihi)/2;
                if( xval < x[imid] )
                        ihi = imid;
                else
                        ilo = imid;
        }
        return ilo;
}
 
template<class T>
inline long hunt_bisect_reverse(const T x[], /* x[n] */
                                long n,
                                T xval)
{
        /* do bisection hunt */
        long ilo = 0, ihi = n-1;
        while( ihi-ilo > 1 )
        {
                long imid = (ilo+ihi)/2;
                if( xval <= x[imid] )
                        ilo = imid;
                else
                        ihi = imid;
        }
        return ilo;
}
 
/*============================================================================*/
 
/* these are the routines in thirdparty_lapack.cpp */
 
/* there are wrappers for lapack linear algebra routines.
 * there are two versions of the lapack routines - a fortran
 * version that I converted to C with forc to use if nothing else is available
 * (included in the Cloudy distribution),
 * and an option to link into an external lapack library that may be optimized
 * for your machine.  By default the tralated C routines will be used.
 * To use your machine's lapack library instead, define the macro
 * LAPACK and link into your library.  This is usually done with a command line
 * switch "-DLAPACK" on the compiler command, and the linker option "-llapack"
 */
 
void getrf_wrapper(long M, long N, double *A, long lda, int32 *ipiv, int32 *info);
 
void getrs_wrapper(char trans, long N, long nrhs, double *A, long lda, int32 *ipiv, double *B, long ldb, int32 *info);
 
void humlik(int n, const realnum x[], realnum y, realnum k[]);
 
realnum FastVoigtH(realnum a, realnum v);
 
// calculates y[i] = H(a,v[i]) as defined in Eq. 9-44 of Mihalas
inline void VoigtH(realnum a, const realnum v[], realnum y[], int n)
{
        if( a <= 0.1f )
        {
                for( int i=0; i < n; ++i )
                        y[i] = FastVoigtH(a,v[i]);
        }
        else
        {
                humlik( n, v, a, y );
        }
}
 
// calculates y[i] = U(a,v[i]) as defined in Eq. 9-45 of Mihalas
inline void VoigtU(realnum a, const realnum v[], realnum y[], int n)
{
        VoigtH( a, v, y, n );
        for( int i=0; i < n; ++i )
                y[i] /= realnum(SQRTPI);
}
 
// VoigtH0(a) returns the value H(a,0) following Eq. 9-44 of Mihalas
inline double VoigtH0(double a)
{
        return erfce(a);
}
 
// VoigtU0(a) returns the value U(a,0) following Eq. 9-45 of Mihalas
inline double VoigtU0(double a)
{
        return erfce(a)/SQRTPI;
}
 
static const unsigned int NMD5 = 32;
 
string MD5file(const char* fnam, access_scheme scheme=AS_DATA_ONLY);
string MD5datafile(const char* fnam, access_scheme scheme=AS_DATA_ONLY);
string MD5string(const string& str);
 
#endif /* THIRDPARTY_H_ */