#ifndef _CONSTS_H_
#define _CONSTS_H_
const double GEOCO = 0.993277; /* geographical to geocentral coordinate */
const double HALF_PI = 1.5707963267949;
const double TWO_PI = 6.2831853071796;
const double FOUR_PI = 12.56637061435916;
const double PI = 3.1415926535898;
const double Y00 = 0.28209479177388; /* sqrt(1/fpi) */
const double Y10 = 0.48860251190292; /* sqrt(3/fpi) */
const double RADIAN = 0.017453293; /* transform degrees into radians */
const double REPRAD = 57.295779513; /* transform radians into degrees */
const double RADIUS = 6371.0; /* radius of the earth in km */
const double RCMB = 3479.95776; /* radius of the earth's CMB */
const double uRCMB = 0.546218452; /* unified radius of the earth's CMB */
const double RMOHO = 6346.62891; /* radius of the earth's Moho */
const double uRMOHO = 0.996174684; /* radius of the earth's Moho */
const double RICMB = 1221.49072; /* radius of the earth's ICMB */
const double R670 = 5701.0; /* radius of 670km discontinuity */
const double uR670 = 0.894835976; /* radius of 670km discontinuity */
const int MAXSTR = 133; /* Maximum number of characters of string */
#endif
#ifndef _NEWTYPE_H
#define _NEWTYPE_H
typedef enum /* a Boolean type, compatible with C++ definition */
{
FALSE, TRUE
} Boolean; /* C++ definition: bool with two values: false, true. */
typedef enum {CHEB_UPPER=14, CHEB_LOWWER, CHEB_WHOLE, B_SPLINE, NO_RADFUN} RADFUN_T;
typedef enum {CHEB_WHOLE_MODEL = 8, SPH_SPLINE_MODEL, BSP_WHOLE_MODEL,
BSP_SPLIT_MODEL, CHEB_SPLIT_MODEL, XV_SPH_SPLINE_MODEL} MODEL3D_T;
#endif
#ifndef _B_SPLINE_H
#define _B_SPLINE_H
class BSpline {
protected:
double h; // average bell size
const static double sixth; // 1/6
public:
BSpline (double h0);
~BSpline() {}
double fdelta(double delta) const;
double f_df_delta(double delta, double *df) const;
double dfddelta(double delta) const;
double f_df_ddf_delta(double delta, double *df, double *ddf) const;
double dfdddelta(double delta) const;
};
// Radial B-spline class.
class RadialBSpline {
protected:
double *xnode; // Radial knot value of B-Spline
int N; // number of radial knots
Boolean increase; // depth or radius. If TRUE, depth; if FALSE, radius.
public:
RadialBSpline(double *xnode0, int N0) : xnode(xnode0), N(N0)
{
if (xnode[N] < xnode[1])
increase = FALSE;
else
increase = TRUE;
}
~RadialBSpline() {}
int BSplKer(double *B, double x) const;
int BSplKer1(double **coeff, double x) const;
int DBBSplKer(double *B, double *DB, double x) const;
int DDBBSplKer(double *B, double *DB, double *DDB, double x) const;
};
// This class assumes that the radial basis function is B-spline.
class RadialBSplineModel {
protected:
int n_rad; // number of radial knots
double *co_rad; // Radial knot value of B-Spline
RadialBSpline *radialBSpline; // radial B-Spline object
public:
RadialBSplineModel (void) {}
RadialBSplineModel (double radialKnot[], int numKnots);
~RadialBSplineModel() { delete radialBSpline; }
};
// This class assumes that the radial basis function is B-spline.
class RadialBSplineSplitModel {
protected:
int n_radUM; // number of radial knots for UM
double *co_radUM; // Radial knot value of B-Spline
RadialBSpline *radialBSplineUM; // radial B-Spline object for UM
int n_radLM; // number of radial knots for LM
double *co_radLM; // Radial knot value of B-Spline
RadialBSpline *radialBSplineLM; // radial B-Spline object
public:
RadialBSplineSplitModel (void) {}
RadialBSplineSplitModel (double radialKnotUM[], int numKnotsUM,
double radialKnotLM[], int numKnotsLM);
~RadialBSplineSplitModel() { delete radialBSplineUM; delete radialBSplineLM; }
};
#endif
#ifndef _MODEL_3D_H_
#define _MODEL_3D_H_
void radfun(RADFUN_T iopt, double *func, double r, int kmax);
void Ylm(double theta, int lmax, double **p);
void dradfun(RADFUN_T iopt, double *func, double *dfunc, double r, int kmax);
void dYlm(double theta, int lmax, double **p, double **dp);
void ddradfun(RADFUN_T iopt, double *func, double *dfunc,
double *ddfunc, double r, int kmax);
void ddYlm(double theta, int lmax, double **p, double **dp, double **ddp);
MODEL3D_T GetModelType(void);
class MODEL_3D {
protected:
MODEL3D_T type; // type of 3-D model
RADFUN_T radfun_t; // radial function type
char name[MAXSTR]; // model name
char *description; // model description
// Euler angles used for rotation in order that the source
// and the receiver are both in the equator
double alpha;
double beta;
double gamma;
virtual void CreateModel3D() = 0;
public:
MODEL_3D() {}
virtual ~MODEL_3D() {}
MODEL_3D(char *model3D_name, MODEL3D_T type, double alpha = 0.0,
double beta = 0.0, double gamma = 0.0,
char* description = NULL);
MODEL3D_T GetType(void) { return type; }
virtual void read_model3D() = 0;
virtual void rotate(double alp, double bet, double gam)
{ alpha = alp; beta = bet; gamma = gam; }
/* Function dv2vonly--
* radius (km)
* xlat,xlong: co-latitude, co-longitude in radians
* output: dV = velocity anomaly divided by average velocity (PREM)
*/
virtual void dv2vonly(double xr, double xlat,
double xlong, double *dV) const = 0;
virtual void dv2v(double xr, double xlat, double xlong,
double *dv, double *dvdx, double *dvdtheta,
double *dvdphi) const = 0;
virtual void alldv2v(double xr, double xlat, double xlong,
double *dv, double *dvdx, double *dvdtheta,
double *dvdphi, double *d2vd2x, double *d2vd2theta,
double *d2vdxdtheta, double *d2vdxdphi,
double *d2vdthetadphi) const = 0;
virtual void free_model() = 0;
// For some kind of 3-D models (such as BSP_SPLINE_T),
// this function is not implemented.
// Judge if the model is expressed in spherical harmonics.
Boolean IsSpheric() {
return (Boolean) (type == CHEB_WHOLE_MODEL || type == BSP_WHOLE_MODEL
|| type == BSP_SPLIT_MODEL);
}
};
// This class assumes that the horizontal basis function is
// spherical harmonics.
class SPHERIC_MODEL_3D : public MODEL_3D {
protected:
double **Alm; // coefficients Alm after radial summation
double **Blm; // coefficients Blm after radial summation
void CreateAB(double xr); // Allocate memory for Alm and Blm
void FreeAB(double xr = 0.5 * (uRMOHO + uR670)); // free memory of Alm and Blm
public:
SPHERIC_MODEL_3D() {}
virtual ~SPHERIC_MODEL_3D() {}
SPHERIC_MODEL_3D(char *model3D_name, MODEL3D_T typ, double alp = 0.0,
double bet = 0.0, double gam = 0.0,
char* descrip = NULL)
: MODEL_3D (model3D_name, typ, alp, bet, gam, descrip) { }
virtual int GetLmax(double xr = 0.5 * (uRMOHO + uR670)) const = 0;
virtual void SumRad(double xr) = 0;
// Given two spherical harmonics models, calculate their
// correlation value at the given radius (normalized)
friend double *Correlation(SPHERIC_MODEL_3D *m1, SPHERIC_MODEL_3D *m2,
double xr);
// Given two spherical harmonics models, calculate their
// correlation value at the given radius (normalized)
friend double *Correlation(SPHERIC_MODEL_3D *m1, SPHERIC_MODEL_3D *m2,
double xr1, double xr2);
};
class CHEB_WHOLE_T : public SPHERIC_MODEL_3D {
protected:
int kmax; // highest radial order
int Lmax; // highest degree of horizonal scales
double ***a; // coefficients Almk
double ***b; // coefficients Blmk
double ***a0; // coefficients Almk without rotation
double ***b0; // coefficients Blmk without rotation
virtual void CreateModel3D(void);
public:
CHEB_WHOLE_T() {}
CHEB_WHOLE_T(char *model3D_name, MODEL3D_T typ = CHEB_WHOLE_MODEL,
double alp = 0.0, double bet = 0.0, double gam = 0.0,
char* descrip = NULL)
: SPHERIC_MODEL_3D (model3D_name, typ, alp, bet, gam, descrip)
{ read_model3D(); }
~CHEB_WHOLE_T();
int get_kmax() { return kmax; }
int GetLmax(double xr = 0.5*(uRMOHO+uR670)) const { return Lmax; }
void read_model3D();
void rotate(double alp, double bet, double gam);
void dv2vonly(double xr, double xlat,
double xlong, double *dV) const;
void dv2v(double xr, double xlat, double xlong,
double *dv, double *dvdx, double *dvdtheta,
double *dvdphi) const;
void alldv2v(double xr, double xlat, double xlong,
double *dv, double *dvdx, double *dvdtheta,
double *dvdphi, double *d2vd2x, double *d2vd2theta,
double *d2vdxdtheta, double *d2vdxdphi,
double *d2vdthetadphi) const;
void free_model();
void SumRad(double xr);
};
class CHEB_SPLIT_T : public SPHERIC_MODEL_3D {
protected:
int kmaxUM; // highest radial order in the upper mantle
int kmaxLM; // highest radial order in the lower mantle
int LmaxUM; // highest degree of horizonal scales in the upper mantle
int LmaxLM; // highest degree of horizonal scales in the lower mantle
double ***aUM; // coefficients Almk in the upper mantle
double ***bUM; // coefficients Blmk in the upper mantle
double ***a0UM; // coefficients Almk without rotation in the upper mantle
double ***b0UM; // coefficients Blmk without rotation in the upper mantle
double ***aLM; // coefficients Almk in the lower mantle
double ***bLM; // coefficients Blmk in the lower mantle
double ***a0LM; // coefficients Almk without rotation in the lower mantle
double ***b0LM; // coefficients Blmk without rotation in the lower mantle
virtual void CreateModel3D(void);
public:
CHEB_SPLIT_T() {}
CHEB_SPLIT_T(char *model3D_name, MODEL3D_T typ = CHEB_SPLIT_MODEL,
double alp = 0.0, double bet = 0.0, double gam = 0.0,
char* descrip = NULL)
: SPHERIC_MODEL_3D (model3D_name, typ, alp, bet, gam, descrip)
{ read_model3D(); }
~CHEB_SPLIT_T();
int getKmaxUM() { return kmaxUM; }
int getKmaxLM() { return kmaxLM; }
int GetLmax(double xr = 0.5 * (uRMOHO+uR670)) const
{
if (xr >= uR670)
return LmaxUM;
else
return LmaxLM;
}
int GetLmaxUM() const { return LmaxUM; }
int GetLmaxLM() const { return LmaxLM; }
void read_model3D();
void rotate(double alp, double bet, double gam);
void dv2vonly(double xr, double xlat,
double xlong, double *dV) const;
void dv2v(double xr, double xlat, double xlong,
double *dv, double *dvdx, double *dvdtheta,
double *dvdphi) const;
void alldv2v(double xr, double xlat, double xlong,
double *dv, double *dvdx, double *dvdtheta,
double *dvdphi, double *d2vd2x, double *d2vd2theta,
double *d2vdxdtheta, double *d2vdxdphi,
double *d2vdthetadphi) const;
void free_model();
void SumRad(double xr);
};
// Mixed whole mantle model: radial function is continuous for the
// whole mantle.
class BSP_WHOLE_T : public RadialBSplineModel, public SPHERIC_MODEL_3D {
protected:
int Lmax; // highest degree of horizonal scales
double ***a; // coefficients Almk
double ***b; // coefficients Blmk
double ***a0; // coefficients Almk without rotation
double ***b0; // coefficients Blmk without rotation
double ***cr; // real coefficients of Cklm (complex coefficients)
double ***ci; // imaginary coefficients of Cklm
virtual void CreateModel3D(void);
public:
BSP_WHOLE_T() {}
BSP_WHOLE_T(int lmax, char RadKnotFile[], char ModelName[]);
BSP_WHOLE_T(char *model3D_name, MODEL3D_T typ = BSP_WHOLE_MODEL,
double alp = 0.0, double bet = 0.0, double gam = 0.0,
char* descrip = NULL)
: SPHERIC_MODEL_3D (model3D_name, typ, alp, bet, gam, descrip)
{
read_model3D();
// initialize RadialBSplineModel
// Other elements co_rad, n_rad are set by read_model3D.
radialBSpline = new RadialBSpline(co_rad, n_rad);
}
~BSP_WHOLE_T();
void ConvCoeff(double x[]);
int GetLmax(double xr = 0.5 * (uRMOHO+uR670)) const { return Lmax; }
int get_n_rad() { return n_rad; }
void convert_to_complex(); // convert a and b into cr and ci
void read_model3D();
void rotate(double alp, double bet, double gam);
void dv2vonly(double xr, double xlat,
double xlong, double *dV) const;
void dv2v(double xr, double xlat, double xlong,
double *dv, double *dvdx, double *dvdtheta,
double *dvdphi) const;
void alldv2v(double xr, double xlat, double xlong,
double *dv, double *dvdx, double *dvdtheta,
double *dvdphi, double *d2vd2x, double *d2vd2theta,
double *d2vdxdtheta, double *d2vdxdphi,
double *d2vdthetadphi) const;
void free_model();
void SumRad(double xr);
};
// Mixed Model: radial direction: B-spline;
// Horizontal direction: spherical harmonics.
class BSP_SPLIT_T : public RadialBSplineSplitModel, public SPHERIC_MODEL_3D {
protected:
int LmaxUM; // highest degree of horizonal scales in the upper mantle
int LmaxLM; // highest degree of horizonal scales in the lower mantle
double ***aUM; // coefficients Almk for the upper mantle
double ***bUM; // coefficients Blmk for the upper mantle
double ***aLM; // coefficients Almk for the lower mantle
double ***bLM; // coefficients Blmk for the lower mantle
double ***a0UM; // coefficients Almk without rotation (upper mantle)
double ***b0UM; // coefficients Blmk without rotation (upper mantle)
double ***a0LM; // coefficients Almk without rotation (lower mantle)
double ***b0LM; // coefficients Blmk without rotation (lower mantle)
double ***crUM; // real coefficients of Cklm (complex coefficients)
double ***ciUM; // imaginary coefficients of Cklm (upper mantle)
double ***crLM; // real coefficients of Cklm (complex coefficients)
double ***ciLM; // imaginary coefficients of Cklm (lower mantle)
virtual void CreateModel3D(void);
public:
BSP_SPLIT_T() {}
BSP_SPLIT_T(int lmax, char RadKnotFile[], char ModelName[]);
BSP_SPLIT_T(int lmaxUM, char RadKnotFile[], char ModelFile[], int lmaxLM);
BSP_SPLIT_T(char *model3D_name, MODEL3D_T typ = BSP_SPLIT_MODEL,
double alp = 0.0, double bet = 0.0, double gam = 0.0,
char* descrip = NULL);
~BSP_SPLIT_T();
void ConvCoeff(double x[]);
int GetLmaxUM() const { return LmaxUM; }
int GetLmax(double xr = 0.5 * (uRMOHO+uR670)) const
{
if (xr >= uR670)
return LmaxUM;
else
return LmaxLM;
}
int GetLmaxLM() const { return LmaxLM; }
int get_n_radUM() { return n_radUM; }
int get_n_radLM() { return n_radLM; }
void convert_to_complex(); // convert a and b into cr and ci
void read_model3D();
// The input file of the 3-D model has no information about its
// radial knot value. This function is going to set them.
// Used by t14.
void read_NO_RADIUS_model3D();
void rotate(double alp, double bet, double gam);
void dv2vonly(double xr, double xlat,
double xlong, double *dV) const;
void dv2v(double xr, double xlat, double xlong,
double *dv, double *dvdx, double *dvdtheta,
double *dvdphi) const;
void alldv2v(double xr, double xlat, double xlong,
double *dv, double *dvdx, double *dvdtheta,
double *dvdphi, double *d2vd2x, double *d2vd2theta,
double *d2vdxdtheta, double *d2vdxdphi,
double *d2vdthetadphi) const;
void free_model();
void SumRad(double xr);
};
// Both radial direction and horizontal direction use B-spline.
class SPH_SPLINE_T : public RadialBSplineModel, public MODEL_3D {
protected:
int n_node; // number of nodes in horizontal direction
double **co_node; // coordinates of horizontal nodes
double **coeff; // coefficients
double avg_dist; // average distance
int **aj_node; // adjacent nodes
virtual void CreateModel3D(void);
public:
SPH_SPLINE_T() {}
~SPH_SPLINE_T() { free_model(); }
SPH_SPLINE_T(char *model3D_name, MODEL3D_T typ = SPH_SPLINE_MODEL,
double alp = 0.0, double bet = 0.0, double gam = 0.0,
char* descrip = NULL)
: MODEL_3D (model3D_name, typ, alp, bet, gam, descrip)
{
read_model3D();
radialBSpline = new RadialBSpline(co_rad, n_rad);
}
void read_model3D();
void dv2vonly(double xr, double xlat,
double xlong, double *dV) const;
void dv2v(double xr, double xlat, double xlong,
double *dv, double *dvdx, double *dvdtheta,
double *dvdphi) const;
void alldv2v(double xr, double xlat, double xlong,
double *dv, double *dvdx, double *dvdtheta,
double *dvdphi, double *d2vd2x, double *d2vd2theta,
double *d2vdxdtheta, double *d2vdxdphi,
double *d2vdthetadphi) const;
double SumSphBspline(double *ar, double theta, double phi) const;
double d_SumSphBspline(double *ar, double *ar_d, double theta,
double phi, double *dv_r, double *dv_dth,
double *dv_dphi) const;
double dd_SumSphBspline(double *ar, double *ar_d,
double *ar_dd, double theta, double phi,
double *dv_dx, double *dv_dth, double *dv_dphi,
double *ddv_d2x, double *ddv_dthth,
double *ddv_dxdth,
double *ddv_dxdph, double *ddv_dthdph) const;
int search_node(double theta, double phi) const;
void free_model();
};
void convert_coef(int nl, double **a, double **b, double **c_r, double **c_i);
void rotate_velocity(double alpha, double beta, double gamma, int nl,
double **a, double **b, double **ap, double **bp);
int BSplKer(double *B, double x, double *xnode, int N);
int DBBSplKer(double *B, double *DB, double x, double *xnode, int N);
int DDBBSplKer(double *B, double *DB, double *DDB,
double x, double *xnode, int N);
void SetNode(double xnode[]);
void SetSplitNode(double xnodeUM[], double xnodeLM[]);
int th_to_index(int n_strip, int th_index, int *total_phi);
// Serves as a virtual constructor for MODEL_3D
MODEL_3D *MakeModel3D(char *model3D_name, MODEL3D_T type,
double alpha = 0.0, double beta = 0.0,
double gamma = 0.0, char *descrip = "3D model");
// Serves as a virtual constructor for SPHERIC_MODEL_3D
SPHERIC_MODEL_3D *MakeSphericModel3D(char *model3D_name, MODEL3D_T type,
double alpha = 0.0, double beta = 0.0,
double gamma = 0.0);
double *Correlation(SPHERIC_MODEL_3D *m1, SPHERIC_MODEL_3D *m2,
double xr);
double *Correlation(SPHERIC_MODEL_3D *m1, SPHERIC_MODEL_3D *m2,
double xr1, double xr2);
#endif
#ifndef _GEN_UTIL_H_
#define _GEN_UTIL_H_
double lat2co(double latitude);
double Geocen2theta(double geocen_latitude);
double long2phi(double longitude);
void euler(double scol, double slon, double rcol, double rlon,
double *delta, double *alpha, double *beta, double *gamma);
void rot(double alpha, double beta, double gamma,
double *theta, double *phi, double *zeta);
void rot(double alpha, double beta, double gamma,
double *theta, double *phi);
void reduce(double *pth, double *pph);
void choldc(double **a, int n, double p[]);
void cholsl(double **a, int n, double p[], double b[], double x[]);
int binary_search(double x, double v[], int n);
void XlnHpi(double Xln[], int L);
void rotylm(double alpha, double beta, double gamma,
int s, int mp, int t, double *bmp);
void pnab(int n, int nalpha, int nbeta, double anorm, double *x,
double *pjac);
void fact(int n1, int n2, double *a);
void exitLog(char *logMessage);
void polcoe(double x[], double y[], int n, double cof[]);
#endif
#ifndef _NR_UTILS_H_
#define _NR_UTILS_H_
static float sqrarg;
#define SQR(a) ((sqrarg=(a)) == 0.0 ? 0.0 : sqrarg*sqrarg)
static double dsqrarg;
#define DSQR(a) ((dsqrarg=(a)) == 0.0 ? 0.0 : dsqrarg*dsqrarg)
static double dmaxarg1,dmaxarg2;
#define DMAX(a,b) (dmaxarg1=(a),dmaxarg2=(b),(dmaxarg1) > (dmaxarg2) ?\
(dmaxarg1) : (dmaxarg2))
static double dminarg1,dminarg2;
#define DMIN(a,b) (dminarg1=(a),dminarg2=(b),(dminarg1) < (dminarg2) ?\
(dminarg1) : (dminarg2))
static float maxarg1,maxarg2;
#define FMAX(a,b) (maxarg1=(a),maxarg2=(b),(maxarg1) > (maxarg2) ?\
(maxarg1) : (maxarg2))
static float minarg1,minarg2;
#define FMIN(a,b) (minarg1=(a),minarg2=(b),(minarg1) < (minarg2) ?\
(minarg1) : (minarg2))
static long lmaxarg1,lmaxarg2;
#define LMAX(a,b) (lmaxarg1=(a),lmaxarg2=(b),(lmaxarg1) > (lmaxarg2) ?\
(lmaxarg1) : (lmaxarg2))
static long lminarg1,lminarg2;
#define LMIN(a,b) (lminarg1=(a),lminarg2=(b),(lminarg1) < (lminarg2) ?\
(lminarg1) : (lminarg2))
static int imaxarg1,imaxarg2;
#define IMAX(a,b) (imaxarg1=(a),imaxarg2=(b),(imaxarg1) > (imaxarg2) ?\
(imaxarg1) : (imaxarg2))
static int iminarg1,iminarg2;
#define IMIN(a,b) (iminarg1=(a),iminarg2=(b),(iminarg1) < (iminarg2) ?\
(iminarg1) : (iminarg2))
#define SIGN(a,b) ((b) >= 0.0 ? fabs(a) : -fabs(a))
void nrerror(char error_text[]);
void error(char error_text[]);
void error(char error_text[], int error_value);
void error(char error_text[], double error_value);
double *vector(long nl, long nh);
int *ivector(long nl, long nh);
unsigned char *cvector(long nl, long nh);
unsigned long *lvector(long nl, long nh);
float *fvector(long nl, long nh);
double **matrix(long nrl, long nrh, long ncl, long nch);
float **fmatrix(long nrl, long nrh, long ncl, long nch);
int **imatrix(long nrl, long nrh, long ncl, long nch);
double **submatrix(double **a, long oldrl, long oldrh, long oldcl, long oldch,
long newrl, long newcl);
double **convert_matrix(double *a, long nrl, long nrh, long ncl, long nch);
double ***f3tensor(long nrl, long nrh, long ncl, long nch, long ndl, long ndh);
void free_vector(double *v, long nl, long nh);
void free_vector(int *v, long nl, long nh);
void free_ivector(int *v, long nl, long nh);
void free_vector(unsigned char *v, long nl, long nh);
void free_vector(unsigned long *v, long nl, long nh);
void free_vector(double *v, long nl, long nh);
void free_vector(float *v, int nl, int nh);
void free_vector(float *v, long nl, long nh);
void free_matrix(double **m, long nrl, long nrh, long ncl, long nch);
void free_matrix(float **m, long nrl, long nrh, long ncl, long nch);
void free_matrix(int **m, long nrl, long nrh, long ncl, long nch);
void free_submatrix(double **b, long nrl, long nrh, long ncl, long nch);
void free_convert_matrix(double **b, long nrl, long nrh, long ncl, long nch);
void free_f3tensor(double ***t, long nrl, long nrh, long ncl, long nch,
long ndl, long ndh);
#endif /* _NR_UTILS_H_ */
#ifndef _MACROS_H_
#define _MACROS_H_
/* test if integer i is even */
inline int EVEN (int i) { return ((i-2*((i)/2)==0) ? 1:0); }
/* test if integer i is odd */
inline int ODD(int i) { return (i & 1); }
/* test sign */
inline int ISIGN(int i) { return ((i<0) ? -1:1); }
/* calculate distance between two points */
inline double ARC_DIST(double costh, double sinth,
double phi1, double theta2, double phi2)
{
return (acos(costh*cos(theta2)+sinth*sin(theta2)*cos(phi1-phi2)));
}
#endif