#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