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_

#include "physconst.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;

/* largest value of fct function corresponding to factorial in six j sjs calculation */
static const int MXDSF = 340;

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 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);

double igam (double a, double x);
double igamc (double a, double x);
double igamc_scaled(double a, double x);

double bessel_k0(double x);
double bessel_k0_scaled(double x);
double bessel_k1(double x);
double bessel_k1_scaled(double x);
void bessel_k0_k1(double x, double* k0val, double* k1val);
void bessel_k0_k1_scaled(double x, double* k0val, double* k1val);

double bessel_i0(double x);
double bessel_i0_scaled(double x);
double bessel_i1(double x);
double bessel_i1_scaled(double x);
void bessel_i0_i1(double x, double* k0val, double* k1val);
void bessel_i0_i1_scaled(double x, double* k0val, double* k1val);

#ifndef HAVE_SINCOS
// this is a GNU extension
inline void sincos(double x, double* s, double* c)
{
        *s = sin(x);
        *c = cos(x);
}
#endif

double e1(double x);

double e1_scaled(double x);

double e2(double t);

/* 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(const vector<T>& x,
                        T xval)
{
        return hunt_bisect(get_ptr(x), x.size(), xval);
}

template<class T>
inline long hunt_bisect(const vector<T>& x,
                        size_t n,
                        T xval)
{
        ASSERT( n <= x.size() );
        return hunt_bisect(get_ptr(x), n, xval);
}

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;
}

template<class T>
inline long hunt_bisect_reverse(const vector<T>& x,
                                T xval)
{
        return hunt_bisect_reverse(get_ptr(x), x.size(), xval);
}

template<class T>
inline long hunt_bisect_reverse(const vector<T>& x,
                                size_t n,
                                T xval)
{
        ASSERT( n <= x.size() );
        return hunt_bisect_reverse(get_ptr(x), n, xval);
}

/*============================================================================*/

/* 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);
void FastVoigtH(realnum a, const realnum v[], realnum y[], size_t n);

// 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 )
        {
                FastVoigtH( a, v, y, n );
        }
        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_DEFAULT);
string MD5datafile(const char* fnam, access_scheme scheme=AS_DEFAULT);
string MD5datastream(fstream& ioFile);
string MD5string(const string& str);

void test_expn();
void chbfit(double a, double b, vector<double>& c, double (*func)(double));

void svd(const int nm, const int m, const int n, double *a, 
                  double *w, bool matu, double *u, bool matv, 
                  double *v, int *ierr, double *rv1);

// Evaluate the Gegenbauer (aka ultraspherical) polynomial C_n^(alpha)(x)
double gegenbauer(long n, double al, double x);

// Diagonal generating function for ultraspherical polynomials,
// following Badnell analysis.  Integer argument to constructor is
// sum of raised and lowered index to the ultraspherical function.
class UltraGen
{
        long   m_a, m_n;
        double m_x, m_c1, m_c2;
public:
UltraGen(long a, double x) : m_a(a), m_n(0), m_x(x), m_c1(0.), m_c2(1.) 
        {}
        void step()
        {
                double c1 = m_c1;
                --m_a;
                m_c1 = 2.*m_a*(m_c2-m_x*c1)/(m_n+2*m_a);
                m_c2 = 2.*m_a*(m_x*m_c2-c1)/(m_n+1);
                ++m_n;
        }
        double val() const
        {
                return m_c2;
        }
};

// this is the triangular inequality for angular momenta in 2*J format
inline bool Triangle2(long J1, long J2, long J3)
{
        return ( J1 >= 0 && J2 >= 0 && J3 >= 0 &&
                 J1 >= abs(J2-J3) && J1 <= J2+J3 &&
                 J1%2 == (J2+J3)%2 );
}

// 6j Wigner evaluation, original routine Fortran Nigel Badnell 6j for Autostructure
double sjs(long j1, long j2, long j3, long l1, long l2, long l3);
// version of the above with less range in j, is exported for testing purposes
double SixJFull(long j1, long j2, long j3, long j4, long j5, long j6);
// Implementation utlizing recursion relation for vector of sixj values,
// should be robust to high j
void rec6j(double *sixcof, double l2, double l3, 
                          double l4, double l5, double l6, double *l1min,
                          double *l1max, double *lmatch, long ndim, long *ier);


#endif /* THIRDPARTY_H_ */