/* function model3D-- * read 3-D velocity perturbation coefficients * rotate the coefficients so that the source and receiver * is on the Equatorial Plane. * It should return the absolute velocity perturbation and * its first, second derivatives with respect to r, theta. */ #include #include #include #include #include #include #include "model3Dlib.h" /* Constructor-- * open model file, and read 3-D model, and return a 3-D model. * All other derived class will call this contructor. */ MODEL_3D::MODEL_3D(char *model3D_name, MODEL3D_T typ, double alp, double bet, double gam, char * descrip) { strcpy(name, model3D_name); type = typ; alpha = alp; beta = bet; gamma = gam; if (descrip) { description = new char[strlen(descrip)+1]; strcpy(description, descrip); } else { description = "No description available."; } // It is illegal to call a pure virtual function in a constructor // A constructor always calls the local copy of the virtual // function, which is a pure virtual function in our case. // read_model3D(); // illegal } BSP_SPLIT_T::BSP_SPLIT_T(char *model3D_name, MODEL3D_T typ, double alp, double bet, double gam, char* descrip) : SPHERIC_MODEL_3D (model3D_name, typ, alp, bet, gam, descrip) { read_model3D(); // initialize RadialBSplineSplitModel // Other elements co_rad, n_rad are set by read_model3D. radialBSplineUM = new RadialBSpline(co_radUM, n_radUM); radialBSplineLM = new RadialBSpline(co_radLM, n_radLM); } // Read the radial value from file, allow memory for model // coefficients. Model coefficients are still undetermined. BSP_WHOLE_T::BSP_WHOLE_T(int lmax, char RadKnotFile[], char ModelFile[]) { strcpy(name, ModelFile); ifstream fin(RadKnotFile); if (!fin) { cerr << "Cannot open: " << RadKnotFile << endl; char msg[MAXSTR]; sprintf(msg, "Cannot open: %s\n", RadKnotFile); exitLog(msg); } fin >> n_rad; co_rad = vector(-1, n_rad + 2); if (!co_rad) { cerr << "BSP_WHOLE_T: no memory for radial coordinates." << endl; exitLog("BSP_WHOLE_T: no memory for radial coordinates.\n"); } for (int i = 1; i <= n_rad; i++) { fin >> co_rad[i]; } co_rad[-1] = co_rad[0] = co_rad[1]; co_rad[n_rad+2] = co_rad[n_rad+1] = co_rad[n_rad]; Lmax = lmax; a = f3tensor(0, Lmax, 0, Lmax, 0, n_rad + 1); if (a == NULL) { cerr << "BSP_WHOLE_T: no memory for 3D coefficients" << endl; exitLog("BSP_WHOLE_T: no memory for 3D coefficients\n"); } b = f3tensor(0, Lmax, 0, Lmax, 0, n_rad + 1); if (b == NULL) { cerr << "BSP_WHOLE_T: no memory for 3D coefficients" << endl; exitLog("BSP_WHOLE_T: no memory for 3D coefficients\n"); } radialBSpline = new RadialBSpline(co_rad, n_rad); int n_radp1 = n_rad + 1; // set unneccessary coefficients to zero for (int L = 0; L <= Lmax; L++) { for (int n = 0; n <= n_rad; n++) { b[L][0][n] = 0.0; b0[L][0][n] = 0.0; } for (int m = 0; m <= L; m++) { a[L][m][0] = 0.0; b[L][m][0] = 0.0; a[L][m][n_radp1] = 0.0; b[L][m][n_radp1] = 0.0; a0[L][m][0] = 0.0; b0[L][m][0] = 0.0; a0[L][m][n_radp1] = 0.0; b0[L][m][n_radp1] = 0.0; } } return; } // Read the radial value from file, allow memory for model // coefficients. Model coefficients are still undetermined. BSP_SPLIT_T::BSP_SPLIT_T(int lmaxUM, char RadKnotFile[], char ModelFile[], int lmaxLM) { strcpy(name, ModelFile); ifstream fin(RadKnotFile); if (!fin) { cerr << "Cannot open: " << RadKnotFile << endl; char msg[MAXSTR]; sprintf(msg, "Cannot open: %s\n", RadKnotFile); exitLog(msg); } fin >> n_radUM >> n_radLM; if (n_radLM < 0 || n_radLM >= 1000 || n_radUM < 0 || n_radUM >= 1000) { cerr << "Split model expected 2 radial knot values" << endl; exitLog("Split model expected 2 radial knot values\n"); } co_radUM = vector(-1, n_radUM + 2); co_radLM = vector(-1, n_radLM + 2); if (co_radUM == NULL || co_radLM == NULL) { cerr << "BSP_SPLIT_T: no memory for radial coordinates." << endl; exitLog("BSP_SPLIT_T: no memory for radial coordinates.\n"); } for (int i = 1; i <= n_radUM; i++) { fin >> co_radUM[i]; } co_radUM[0] = co_radUM[1]; co_radUM[-1] = co_radUM[0]; co_radUM[n_radUM+1] = co_radUM[n_radUM]; co_radUM[n_radUM+2] = co_radUM[n_radUM+1]; co_radLM[0] = co_radLM[1]; co_radLM[-1] = co_radLM[1]; co_radLM[n_radLM+1] = co_radLM[n_radLM]; co_radLM[n_radLM+2] = co_radLM[n_radLM]; LmaxUM = lmaxUM; LmaxLM = lmaxLM; aUM = f3tensor(0, LmaxUM, 0, LmaxUM, 0, n_radUM+1); aLM = f3tensor(0, LmaxLM, 0, LmaxLM, 0, n_radLM+1); if (aUM == NULL || aLM == NULL) { fprintf(stderr, "BSP_SPLIT_T: no memory for 3D coefficients\n"); exitLog("BSP_SPLIT_T: no memory for 3D coefficients\n"); } a0UM = f3tensor(0, LmaxUM, 0, LmaxUM, 0, n_radUM+1); a0LM = f3tensor(0, LmaxLM, 0, LmaxLM, 0, n_radLM+1); if (a0UM == NULL || a0LM == NULL) { fprintf(stderr, "BSP_SPLIT_T: no memory for 3D coefficients\n"); exitLog("BSP_SPLIT_T: no memory for 3D coefficients\n"); } bUM = f3tensor(0, LmaxUM, 0, LmaxUM, 0, n_radUM+1); bLM = f3tensor(0, LmaxLM, 0, LmaxLM, 0, n_radLM+1); if (bUM == NULL || bLM == NULL) { fprintf(stderr, "BSP_SPLIT_T: no memory for 3D coefficients\n"); exitLog("BSP_SPLIT_T: no memory for 3D coefficients\n"); } b0UM = f3tensor(0, LmaxUM, 0, LmaxUM, 0, n_radUM+1); b0LM = f3tensor(0, LmaxLM, 0, LmaxLM, 0, n_radLM+1); if (b0UM == NULL || b0LM == NULL) { fprintf(stderr, "BSP_SPLIT_T: no memory for 3D coefficients\n"); exitLog("BSP_SPLIT_T: no memory for 3D coefficients\n"); } radialBSplineUM = new RadialBSpline(co_radUM, n_radUM); radialBSplineLM = new RadialBSpline(co_radLM, n_radLM); int n_radp1 = n_radUM + 1; // set unneccessary coefficients to zero int L; for (L = 0; L <= LmaxUM; L++) { for (int m = 0; m <= L; m++) { aUM[L][m][0] = 0.0; bUM[L][m][0] = 0.0; aUM[L][m][n_radp1] = 0.0; bUM[L][m][n_radp1] = 0.0; a0UM[L][m][0] = 0.0; b0UM[L][m][0] = 0.0; a0UM[L][m][n_radp1] = 0.0; b0UM[L][m][n_radp1] = 0.0; } } n_radp1 = n_radLM + 1; for (L = 0; L <= LmaxLM; L++) { for (int m = 0; m <= L; m++) { aLM[L][m][0] = 0.0; bLM[L][m][0] = 0.0; aLM[L][m][n_radp1] = 0.0; bLM[L][m][n_radp1] = 0.0; a0LM[L][m][0] = 0.0; b0LM[L][m][0] = 0.0; a0LM[L][m][n_radp1] = 0.0; b0LM[L][m][n_radp1] = 0.0; } } return; } // Read the radial value from file, allow memory for model coefficients // Model coefficients are still undetermined. // This constructor assume that LmaxUM = LmaxLM BSP_SPLIT_T::BSP_SPLIT_T(int lmax, char RadKnotFile[], char ModelFile[]) { strcpy(name, ModelFile); ifstream fin(RadKnotFile); if (!fin) { cerr << "Cannot open: " << RadKnotFile << endl; char msg[MAXSTR]; sprintf(msg, "Cannot open: %s\n", RadKnotFile); exitLog(msg); } fin >> n_radUM >> n_radLM; if (n_radLM < 0 || n_radLM >= 1000 || n_radUM < 0 || n_radUM >= 1000) { cerr << "Split model expected 2 radial knot values" << endl; exitLog("Split model expected 2 radial knot values\n"); } co_radUM = vector(-1, n_radUM + 2); co_radLM = vector(-1, n_radLM + 2); if (co_radUM == NULL || co_radLM == NULL) { cerr << "BSP_SPLIT_T: no memory for radial coordinates." << endl; exitLog("BSP_SPLIT_T: no memory for radial coordinates.\n"); } int i; for (i = 1; i <= n_radUM; i++) { fin >> co_radUM[i]; if (co_radUM[i] > (RADIUS-RCMB) || co_radUM[i] < (RADIUS-RMOHO)) { fprintf(stderr, "Wrong with your 3-D model coefficients file\n"); fprintf(stderr, "Radial knot value co_radUM[%d] = %f\n", i, co_radUM[i]); char msg[MAXSTR]; sprintf(msg, "Wrong with your 3-D model coefficients file, co_radUM[%d] = %f\n", i, co_radUM[i]); exitLog(msg); } } for (i = 1; i <= n_radLM; i++) { fin >> co_radLM[i]; if (co_radLM[i] > (RADIUS-RCMB) || co_radLM[i] < (RADIUS-R670)) { fprintf(stderr, "Wrong with your 3-D model coefficients file\n"); fprintf(stderr, "Radial knot value co_radLM[%d] = %f\n", i, co_radLM[i]); char msg[MAXSTR]; sprintf(msg, "Wrong with your 3-D model coefficients file, co_radLM[%d] = %f\n", i, co_radLM[i]); exitLog(msg); } } if (co_radLM[1] != (RADIUS-R670)) { fprintf(stderr, "Wrong with your 3-D model coefficients file\n"); fprintf(stderr, "Radial knot value R670 = %f\n", co_radLM[1]); char msg[MAXSTR]; sprintf(msg, "Wrong with your 3-D model coefficients file, R670 = %f\n", co_radLM[1]); exitLog(msg); } co_radUM[0] = co_radUM[1]; co_radUM[-1] = co_radUM[0]; co_radUM[n_radUM+1] = co_radUM[n_radUM]; co_radUM[n_radUM+2] = co_radUM[n_radUM+1]; co_radLM[0] = co_radLM[1]; co_radLM[-1] = co_radLM[1]; co_radLM[n_radLM+1] = co_radLM[n_radLM]; co_radLM[n_radLM+2] = co_radLM[n_radLM]; LmaxUM = LmaxLM = lmax; aUM = f3tensor(0, LmaxUM, 0, LmaxUM, 0, n_radUM+1); aLM = f3tensor(0, LmaxLM, 0, LmaxLM, 0, n_radLM+1); if (aUM == NULL || aLM == NULL) { fprintf(stderr, "BSP_SPLIT_T: no memory for 3D coefficients\n"); exitLog("BSP_SPLIT_T: no memory for 3D coefficients\n"); } a0UM = f3tensor(0, LmaxUM, 0, LmaxUM, 0, n_radUM+1); a0LM = f3tensor(0, LmaxLM, 0, LmaxLM, 0, n_radLM+1); if (a0UM == NULL || a0LM == NULL) { fprintf(stderr, "BSP_SPLIT_T: no memory for 3D coefficients\n"); exitLog("BSP_SPLIT_T: no memory for 3D coefficients\n"); } bUM = f3tensor(0, LmaxUM, 0, LmaxUM, 0, n_radUM+1); bLM = f3tensor(0, LmaxLM, 0, LmaxLM, 0, n_radLM+1); if (bUM == NULL || bLM == NULL) { fprintf(stderr, "BSP_SPLIT_T: no memory for 3D coefficients\n"); exitLog("BSP_SPLIT_T: no memory for 3D coefficients\n"); } b0UM = f3tensor(0, LmaxUM, 0, LmaxUM, 0, n_radUM+1); b0LM = f3tensor(0, LmaxLM, 0, LmaxLM, 0, n_radLM+1); if (b0UM == NULL || b0LM == NULL) { fprintf(stderr, "BSP_SPLIT_T: no memory for 3D coefficients\n"); exitLog("BSP_SPLIT_T: no memory for 3D coefficients\n"); } radialBSplineUM = new RadialBSpline(co_radUM, n_radUM); radialBSplineLM = new RadialBSpline(co_radLM, n_radLM); int n_radp1 = n_radUM + 1; // set unneccessary coefficients to zero int L; for (L = 0; L <= LmaxUM; L++) { for (int m = 0; m <= L; m++) { aUM[L][m][0] = 0.0; bUM[L][m][0] = 0.0; aUM[L][m][n_radp1] = 0.0; bUM[L][m][n_radp1] = 0.0; a0UM[L][m][0] = 0.0; b0UM[L][m][0] = 0.0; a0UM[L][m][n_radp1] = 0.0; b0UM[L][m][n_radp1] = 0.0; } } n_radp1 = n_radLM + 1; for (L = 0; L <= LmaxLM; L++) { for (int m = 0; m <= L; m++) { aLM[L][m][0] = 0.0; bLM[L][m][0] = 0.0; aLM[L][m][n_radp1] = 0.0; bLM[L][m][n_radp1] = 0.0; a0LM[L][m][0] = 0.0; b0LM[L][m][0] = 0.0; a0LM[L][m][n_radp1] = 0.0; b0LM[L][m][n_radp1] = 0.0; } } return; } void BSP_WHOLE_T::free_model(void) { free_vector(co_rad, -1, n_rad + 2); free_f3tensor(a, 0, Lmax, 0, Lmax, 0, n_rad + 1); free_f3tensor(b, 0, Lmax, 0, Lmax, 0, n_rad + 1); free_f3tensor(a0, 0, Lmax, 0, Lmax, 0, n_rad + 1); free_f3tensor(b0, 0, Lmax, 0, Lmax, 0, n_rad + 1); } BSP_WHOLE_T::~BSP_WHOLE_T() { free_vector(co_rad, -1, n_rad + 2); free_f3tensor(a, 0, Lmax, 0, Lmax, 0, n_rad + 1); free_f3tensor(b, 0, Lmax, 0, Lmax, 0, n_rad + 1); free_f3tensor(a0, 0, Lmax, 0, Lmax, 0, n_rad + 1); free_f3tensor(b0, 0, Lmax, 0, Lmax, 0, n_rad + 1); } BSP_SPLIT_T::~BSP_SPLIT_T() { free_model(); } CHEB_SPLIT_T::~CHEB_SPLIT_T() { free_model(); } void BSP_SPLIT_T::free_model(void) { free_vector(co_radUM, -1, n_radUM + 2); free_f3tensor(aUM, 0, LmaxUM, 0, LmaxUM, 0, n_radUM + 1); free_f3tensor(bUM, 0, LmaxUM, 0, LmaxUM, 0, n_radUM + 1); free_f3tensor(a0UM, 0, LmaxUM, 0, LmaxUM, 0, n_radUM + 1); free_f3tensor(b0UM, 0, LmaxUM, 0, LmaxUM, 0, n_radUM + 1); free_vector(co_radLM, -1, n_radLM + 2); free_f3tensor(aLM, 0, LmaxLM, 0, LmaxLM, 0, n_radLM + 1); free_f3tensor(bLM, 0, LmaxLM, 0, LmaxLM, 0, n_radLM + 1); free_f3tensor(a0LM, 0, LmaxLM, 0, LmaxLM, 0, n_radLM + 1); free_f3tensor(b0LM, 0, LmaxLM, 0, LmaxLM, 0, n_radLM + 1); return; } void CHEB_SPLIT_T::free_model(void) { free_f3tensor(aUM, 0, LmaxUM, 0, LmaxUM, 0, kmaxUM + 1); free_f3tensor(bUM, 0, LmaxUM, 0, LmaxUM, 0, kmaxUM + 1); free_f3tensor(a0UM, 0, LmaxUM, 0, LmaxUM, 0, kmaxUM + 1); free_f3tensor(b0UM, 0, LmaxUM, 0, LmaxUM, 0, kmaxUM + 1); free_f3tensor(aLM, 0, LmaxLM, 0, LmaxLM, 0, kmaxLM + 1); free_f3tensor(bLM, 0, LmaxLM, 0, LmaxLM, 0, kmaxLM + 1); free_f3tensor(a0LM, 0, LmaxLM, 0, LmaxLM, 0, kmaxLM + 1); free_f3tensor(b0LM, 0, LmaxLM, 0, LmaxLM, 0, kmaxLM + 1); return; } void CHEB_WHOLE_T::read_model3D() { FILE * fp; int kmaxp2, irad; int i, k, L, m; double **ak, **bk, **akm, **bkm; char Line[MAXSTR], RadialMask[MAXSTR], Lmask[MAXSTR]; if ((fp = fopen(name, "r")) == NULL) { fprintf(stderr, "model3D: can't open: %s\n", name); char msg[MAXSTR]; sprintf(msg, "model3D: can't open: %s\n", name); exitLog(msg); } fgets(Line, MAXSTR, fp); sscanf(Line, "%d %s %d %s %d", &Lmax, Lmask, &kmaxp2, RadialMask, &irad); kmax = kmaxp2 - 3; if (kmaxp2 <= 0 || Lmax <= 0) { fprintf(stderr, "Model type is not what expected!\n"); exitLog("Model type is not what expected!\n"); } // Model type verification for (i = 0; i <= Lmax; i++) { if (Lmask[i] != '0' && Lmask[i] != '1') { fprintf(stderr, "Either degree mask or the model type"); fprintf(stderr, " not what expected!\n"); exitLog("Either degree mask or the model type not what expected!\n"); } } if (irad == 4) { printf("B-Spline is used in radial function.\n"); radfun_t = B_SPLINE; } else if (irad == 2) { printf("ChebyChev Polynomial is used in radial function.\n"); radfun_t = CHEB_WHOLE; } else if (irad == 3) { printf("No radial function is used in this model.\n"); radfun_t = NO_RADFUN; } else { fprintf(stderr, "radial function type is unknown!\n"); fprintf(stderr, "Or model type is what expected!\n"); exitLog("radial function type is unknown! Or model type is what expected!\n"); } printf("Maximum degree: %d, Maximum Radial Order: %d\n", Lmax, kmax); double dummy; // skip crust and 670 discontinuity coefficients for (k = 0; k < 2; k++) { if (RadialMask[k] =='0') continue; if (k == 0) printf("Skip crustal coefficients\n"); else printf("Skip 670km discontinuity coefficients\n"); for (L = 0; L <= Lmax; L++) { if (Lmask[L] == '0') continue; switch (L) { case 0: fscanf(fp, "%le", &dummy); break; default: for (m = 0; m <= L; m++) { switch (m) { case 0: fscanf(fp, "%le", &dummy); break; default: fscanf(fp, "%le %le", &dummy, &dummy); break; } } } } } ak = matrix(0, Lmax, 0, Lmax); bk = matrix(0, Lmax, 0, Lmax); akm = matrix(0, Lmax, 0, Lmax); bkm = matrix(0, Lmax, 0, Lmax); CreateModel3D(); for (k = 0; k <= kmax; k++) { if (RadialMask[k+2] == '0') continue; else if (RadialMask[k+2] !='1') { fprintf(stderr, "Your input model is not whole-mantle earth model\n"); fprintf(stderr, "Model type is not what I expected!\n"); exitLog("Model type is not what I expected!\n"); } for (L = 0; L <= Lmax; L++) { if (Lmask[L] == '0') continue; switch (L) { case 0: fscanf(fp, "%le", &ak[L][0]); break; default: for (m = 0; m <= L; m++) { switch (m) { case 0: fscanf(fp, "%le", &ak[L][0]); break; default: fscanf(fp, "%le %le", &ak[L][m], &bk[L][m]); break; } } } a0[L][0][k] = ak[L][0]; if (L > 0) { for (m = 1; m <= L; m++) { a0[L][m][k] = ak[L][m]; b0[L][m][k] = bk[L][m]; } } } if (alpha != 0.0 || beta != 0.0 || gamma != 0.0) { rotate_velocity(alpha, beta, gamma, Lmax, ak, bk, akm, bkm); for (L = 0; L <= Lmax; L++) { a[L][0][k] = akm[L][0]; if (L > 0) { for (m = 1; m <= L; m++) { a[L][m][k] = akm[L][m]; b[L][m][k] = bkm[L][m]; } } } } else { for (L = 0; L <= Lmax; L++) { a[L][0][k] = ak[L][0]; if (L > 0) { for (m = 1; m <= L; m++) { a[L][m][k] = ak[L][m]; b[L][m][k] = bk[L][m]; } } } } } free_matrix(ak, 0, Lmax, 0, Lmax); free_matrix(bk, 0, Lmax, 0, Lmax); free_matrix(akm, 0, Lmax, 0, Lmax); free_matrix(bkm, 0, Lmax, 0, Lmax); fclose(fp); } void SPH_SPLINE_T::read_model3D() { FILE * fp; int i, j, n, dum[7]; double lat, lon; double th, ph; double dist, theta1, phi1, theta2, phi2; int np, nb; char Line[MAXSTR]; if ((fp = fopen(name, "r")) == NULL) { fprintf(stderr, "model3D: can't open: %s\n", name); char msg[MAXSTR]; sprintf(msg, "model3D: can't open: %s\n", name); exitLog(msg); } fgets(Line, MAXSTR, fp); sscanf(Line, "%d %d", &n_node, &n_rad); if (n_rad <= 0 || n_rad >= 100000 || n_node <= 0) { fprintf(stderr, "The model type is not what expected!\n"); exitLog("The model type is not what expected!\n"); } printf("B-spline model (both radial and horizontal directions).\n"); printf("Number of horizontal nodes: %d\n", n_node); printf("Number of radial nodes: %d\n", n_rad); /* Allocate memory for model parameters, e.g. coefficients */ CreateModel3D(); /* read coordinates of horizontal nodes */ for (i = 1; i <= n_node; i++) { fgets(Line, MAXSTR, fp); sscanf(Line, "%*d %lf %lf %d %d %d %d %d %d %d", &lat, &lon, &n, &dum[1], &dum[2], &dum[3], &dum[4], &dum[5], &dum[6]); // Model type verification if (fabs(lat) > 90.0 || fabs(lon) > 180.0) { fprintf(stderr, "model type is not what expected\n"); exitLog("model type is not what expected\n"); } aj_node[i][0] = n; for (j = 1; j <= n; j++) { if (dum[j] < 0 || dum[j] > n_node) { fprintf(stderr, "Either model type or model format "); fprintf(stderr, "is not what expected\n"); exitLog("Either model type or model format is not what expected\n"); } aj_node[i][j] = dum[j]; } if (n < 6) aj_node[i][6] = 0; /* transform the geocentral coordinate into spherical coordinate */ th = Geocen2theta(lat); co_node[i][1] = th; /* unrotated coordinates */ ph = long2phi(lon); co_node[i][2] = ph; rot(alpha, beta, gamma, &th, &ph); co_node[i][3] = th; /* rotated coordinates */ co_node[i][4] = ph; } for (j = 1; j <= n_rad; j++) { fscanf(fp, "%*d %lf", &co_rad[j]); /* radial knots (depth) */ if (type == XV_SPH_SPLINE_MODEL) // tranform into depth co_rad[j] = RADIUS - co_rad[j]; if (co_rad[j] > RCMB) { fprintf(stderr, "Depth in the inner core: model type not expected\n"); exitLog("Depth in the inner core: model type not expected\n"); } for (i = 1; i <= n_node; i++) fscanf(fp, "%lf", &coeff[j][i]); /* model coefficients */ } co_rad[-1] = co_rad[0] = co_rad[1]; co_rad[n_rad+2] = co_rad[n_rad+1] = co_rad[n_rad]; /* Calculate average distance between nodes */ for (n = 0, dist = 0.0, i = 1; i <= n_node; i++) { double sinth, costh; np = aj_node[i][0]; /* total number of neighboring points */ theta1 = co_node[i][1]; phi1 = co_node[i][2]; costh = cos(theta1); sinth = sin(theta1); for (j = 1; j <= np; j++) { nb = aj_node[i][j]; if (nb > i) { theta2 = co_node[nb][1]; phi2 = co_node[nb][2]; dist += ARC_DIST(costh, sinth, phi1, theta2, phi2); n++; } } } avg_dist = dist / (double) n; printf("Average distance: %f\n", avg_dist * REPRAD); fclose(fp); } /* function CreateModel3D-- * allocate memory for a new 3-D model; * return a pointer to the new 3-D model. */ void CHEB_WHOLE_T::CreateModel3D() { a = f3tensor(0, Lmax, 0, Lmax, 0, kmax); if (a == NULL) { fprintf(stderr, "CHEB_WHOLE_create: no memory for 3D coefficients\n"); exitLog("CHEB_WHOLE_create: no memory for 3D coefficients\n"); } a0 = f3tensor(0, Lmax, 0, Lmax, 0, kmax); if (a0 == NULL) { fprintf(stderr, "CHEB_WHOLE_create: no memory for 3D coefficients\n"); exitLog("CHEB_WHOLE_create: no memory for 3D coefficients\n"); } b = f3tensor(0, Lmax, 0, Lmax, 0, kmax); if (b == NULL) { fprintf(stderr, "CHEB_WHOLE_create: no memory for 3D coefficients\n"); exitLog("CHEB_WHOLE_create: no memory for 3D coefficients\n"); } b0 = f3tensor(0, Lmax, 0, Lmax, 0, kmax); if (b0 == NULL) { fprintf(stderr, "CHEB_WHOLE_create: no memory for 3D coefficients\n"); exitLog("CHEB_WHOLE_create: no memory for 3D coefficients\n"); } // This value (B for m=0) does not exist in spherical harmonics. // I initialize it in case of bugs in summation of spherical harmonics. for (int L = 0; L <= Lmax; L++) { for (int k = 0; k <= kmax; k++) b[L][0][k] = b0[L][0][k] = 0.0; } return; } /* function CreateModel3D-- * allocate memory for a new 3-D model; * return a pointer to the new 3-D model. */ void BSP_WHOLE_T::CreateModel3D() { co_rad = vector(-1, n_rad + 2); a = f3tensor(0, Lmax, 0, Lmax, 0, n_rad+1); if (a == NULL) { fprintf(stderr, "BSP_WHOLE_create: no memory for 3D coefficients\n"); exitLog("BSP_WHOLE_create: no memory for 3D coefficients\n"); } a0 = f3tensor(0, Lmax, 0, Lmax, 0, n_rad+1); if (a0 == NULL) { fprintf(stderr, "BSP_WHOLE_create: no memory for 3D coefficients\n"); exitLog("BSP_WHOLE_create: no memory for 3D coefficients\n"); } b = f3tensor(0, Lmax, 0, Lmax, 0, n_rad+1); if (b == NULL) { fprintf(stderr, "BSP_WHOLE_create: no memory for 3D coefficients\n"); exitLog("BSP_WHOLE_create: no memory for 3D coefficients\n"); } b0 = f3tensor(0, Lmax, 0, Lmax, 0, n_rad+1); if (b0 == NULL) { fprintf(stderr, "BSP_WHOLE_create: no memory for 3D coefficients\n"); exitLog("BSP_WHOLE_create: no memory for 3D coefficients\n"); } int n_radp1 = n_rad + 1; // set unneccessary coefficients to zero for (int L = 0; L <= Lmax; L++) { for (int m = 0; m <= L; m++) { a[L][m][0] = 0.0; b[L][m][0] = 0.0; a[L][m][n_radp1] = 0.0; b[L][m][n_radp1] = 0.0; a0[L][m][0] = 0.0; b0[L][m][0] = 0.0; a0[L][m][n_radp1] = 0.0; b0[L][m][n_radp1] = 0.0; } } return; } /* function CreateModel3D-- * allocate memory for a new 3-D model; * return a pointer to the new 3-D model. */ void BSP_SPLIT_T::CreateModel3D() { co_radUM = vector(-1, n_radUM + 2); co_radLM = vector(-1, n_radLM + 2); aUM = f3tensor(0, LmaxUM, 0, LmaxUM, 0, n_radUM+1); aLM = f3tensor(0, LmaxLM, 0, LmaxLM, 0, n_radLM+1); if (aUM == NULL || aLM == NULL) { fprintf(stderr, "BSP_SPLIT_create: no memory for 3D coefficients\n"); exitLog("BSP_SPLIT_create: no memory for 3D coefficients\n"); } a0UM = f3tensor(0, LmaxUM, 0, LmaxUM, 0, n_radUM+1); a0LM = f3tensor(0, LmaxLM, 0, LmaxLM, 0, n_radLM+1); if (a0UM == NULL || a0LM == NULL) { fprintf(stderr, "BSP_SPLIT_create: no memory for 3D coefficients\n"); exitLog("BSP_SPLIT_create: no memory for 3D coefficients\n"); } bUM = f3tensor(0, LmaxUM, 0, LmaxUM, 0, n_radUM+1); bLM = f3tensor(0, LmaxLM, 0, LmaxLM, 0, n_radLM+1); if (bUM == NULL || bLM == NULL) { fprintf(stderr, "BSP_SPLIT_create: no memory for 3D coefficients\n"); exitLog("BSP_SPLIT_create: no memory for 3D coefficients\n"); } b0UM = f3tensor(0, LmaxUM, 0, LmaxUM, 0, n_radUM+1); b0LM = f3tensor(0, LmaxLM, 0, LmaxLM, 0, n_radLM+1); if (b0UM == NULL || b0LM == NULL) { fprintf(stderr, "BSP_SPLIT_create: no memory for 3D coefficients\n"); exitLog("BSP_SPLIT_create: no memory for 3D coefficients\n"); } int n_radp1 = n_radUM + 1; // set unneccessary coefficients to zero int L; for (L = 0; L <= LmaxUM; L++) { for (int m = 0; m <= L; m++) { aUM[L][m][0] = 0.0; bUM[L][m][0] = 0.0; aUM[L][m][n_radp1] = 0.0; bUM[L][m][n_radp1] = 0.0; a0UM[L][m][0] = 0.0; b0UM[L][m][0] = 0.0; a0UM[L][m][n_radp1] = 0.0; b0UM[L][m][n_radp1] = 0.0; } } n_radp1 = n_radLM + 1; for (L = 0; L <= LmaxLM; L++) { for (int m = 0; m <= L; m++) { aLM[L][m][0] = 0.0; bLM[L][m][0] = 0.0; aLM[L][m][n_radp1] = 0.0; bLM[L][m][n_radp1] = 0.0; a0LM[L][m][0] = 0.0; b0LM[L][m][0] = 0.0; a0LM[L][m][n_radp1] = 0.0; b0LM[L][m][n_radp1] = 0.0; } } return; } /* function CreateModel3D-- * allocate memory for a new 3-D model; * return a pointer to the new 3-D model. */ void CHEB_SPLIT_T::CreateModel3D() { aUM = f3tensor(0, LmaxUM, 0, LmaxUM, 0, kmaxUM+1); aLM = f3tensor(0, LmaxLM, 0, LmaxLM, 0, kmaxLM+1); if (aUM == NULL || aLM == NULL) { fprintf(stderr, "BSP_SPLIT_create: no memory for 3D coefficients\n"); exitLog("BSP_SPLIT_create: no memory for 3D coefficients\n"); } a0UM = f3tensor(0, LmaxUM, 0, LmaxUM, 0, kmaxUM+1); a0LM = f3tensor(0, LmaxLM, 0, LmaxLM, 0, kmaxLM+1); if (a0UM == NULL || a0LM == NULL) { fprintf(stderr, "CHEB_SPLIT_create: no memory for 3D coefficients\n"); exitLog("CHEB_SPLIT_create: no memory for 3D coefficients\n"); } bUM = f3tensor(0, LmaxUM, 0, LmaxUM, 0, kmaxUM+1); bLM = f3tensor(0, LmaxLM, 0, LmaxLM, 0, kmaxLM+1); if (bUM == NULL || bLM == NULL) { fprintf(stderr, "CHEB_SPLIT_create: no memory for 3D coefficients\n"); exitLog("CHEB_SPLIT_create: no memory for 3D coefficients\n"); } b0UM = f3tensor(0, LmaxUM, 0, LmaxUM, 0, kmaxUM+1); b0LM = f3tensor(0, LmaxLM, 0, LmaxLM, 0, kmaxLM+1); if (b0UM == NULL || b0LM == NULL) { fprintf(stderr, "CHEB_SPLIT_create: no memory for 3D coefficients\n"); exitLog("CHEB_SPLIT_create: no memory for 3D coefficients\n"); } int kmaxp1 = kmaxUM + 1; // set unneccessary coefficients to zero int L; for (L = 0; L <= LmaxUM; L++) { for (int m = 0; m <= L; m++) { aUM[L][m][0] = 0.0; bUM[L][m][0] = 0.0; aUM[L][m][kmaxp1] = 0.0; bUM[L][m][kmaxp1] = 0.0; a0UM[L][m][0] = 0.0; b0UM[L][m][0] = 0.0; a0UM[L][m][kmaxp1] = 0.0; b0UM[L][m][kmaxp1] = 0.0; } } kmaxp1 = kmaxLM + 1; for (L = 0; L <= LmaxLM; L++) { for (int m = 0; m <= L; m++) { aLM[L][m][0] = 0.0; bLM[L][m][0] = 0.0; aLM[L][m][kmaxp1] = 0.0; bLM[L][m][kmaxp1] = 0.0; a0LM[L][m][0] = 0.0; b0LM[L][m][0] = 0.0; a0LM[L][m][kmaxp1] = 0.0; b0LM[L][m][kmaxp1] = 0.0; } } return; } /* function CHEB_WHOLE_free-- * free memory for an old 3-D model; */ void CHEB_WHOLE_T::free_model() { free_f3tensor(a, 0, Lmax, 0, Lmax, 0, kmax); free_f3tensor(b, 0, Lmax, 0, Lmax, 0, kmax); free_f3tensor(a0, 0, Lmax, 0, Lmax, 0, kmax); free_f3tensor(b0, 0, Lmax, 0, Lmax, 0, kmax); return; } /* function CHEB_WHOLE_free-- * free memory for an old 3-D model; */ CHEB_WHOLE_T::~CHEB_WHOLE_T() { free_f3tensor(a, 0, Lmax, 0, Lmax, 0, kmax); free_f3tensor(b, 0, Lmax, 0, Lmax, 0, kmax); free_f3tensor(a0, 0, Lmax, 0, Lmax, 0, kmax); free_f3tensor(b0, 0, Lmax, 0, Lmax, 0, kmax); } /* function CreateModel3D-- * allocate memory for a new SPH_SPLINE_T model; * return a pointer to the new model. */ void SPH_SPLINE_T::CreateModel3D(void) { radfun_t = B_SPLINE; co_rad = vector(-1, n_rad + 2); if (co_rad == NULL) { fprintf(stderr, "SPH_SPLINE_create: no memory for radial coordinates.\n"); exitLog("SPH_SPLINE_create: no memory for radial coordinates.\n"); } co_node = matrix(1, n_node, 1, 4); if (co_node == NULL) { fprintf(stderr, "SPH_SPLINE_create: no memory for horizontal coordinates.\n"); exitLog("SPH_SPLINE_create: no memory for horizontal coordinates.\n"); } /* The exact matrix should be matrix(1, n_rad, 1, n_node). * I wasted a little memory here to gain speed later on. * I don't need any logic sentences on 1 or n_node. */ coeff = matrix(0, n_rad, 1, n_node + 1); if (coeff == NULL) { fprintf(stderr, "SPH_SPLINE_create: no memory for coefficients.\n"); exitLog("SPH_SPLINE_create: no memory for coefficients.\n"); } aj_node = imatrix(1, n_node, 0, 7); if (aj_node == NULL) { fprintf(stderr, "SPH_SPLINE_create: no memory for adjacent nodes.\n"); exitLog("SPH_SPLINE_create: no memory for adjacent nodes.\n"); } return; } /* function SPH_SPLINE_free-- * free memory for an old SPH_SPLINE_T model; */ void SPH_SPLINE_T::free_model(void) { free_vector(co_rad, -1, n_rad + 2); free_matrix(co_node, 1, n_node, 1, 4); free_matrix(coeff, 0, n_rad, 1, n_node + 1); free_matrix(aj_node, 1, n_node, 0, 7); return; } /* Function dv2vonly-- * radius (km) * xlat,xlong: co-latitude, co-longitude in radians * output: dVs = velocity anomaly divided by average velocity (PREM) */ void SPH_SPLINE_T::dv2vonly(double xr, double xlat, double xlong, double *dVs) const { double rr = xr * RADIUS; double B0, B1, B2, B3; double B[4], *ar; int km1, kp1, kp2; int k, i; *dVs = 0.0; if (rr < RCMB || rr > RMOHO) return; k = radialBSpline->BSplKer(B, RADIUS - rr); if (k == 0) return; ar = vector(1, n_node); kp1 = k + 1; kp2 = k + 2; km1 = k - 1; B0 = B[0]; B1 = B[1]; B2 = B[2]; B3 = B[3]; double tmp; for (i = 1; i <= n_node; i++) { /* for all horizontal nodes */ tmp = coeff[k][i] * B1 + coeff[kp1][i] * B2; if (k != 1) tmp += coeff[km1][i] * B0; if (kp1 != n_rad) tmp += coeff[kp2][i] * B3; ar[i] = tmp; } *dVs = SumSphBspline(ar, xlat, xlong); free_vector(ar, 1, n_node); return; } // Do the radial summation for spherical harmonics model. // The result is in coeffcients Alm and Blm. void CHEB_SPLIT_T::SumRad(double xr) { double rr = xr * RADIUS; int m, L, k, Lmax, kmax; double ***a, ***b; if (rr < RCMB || rr > RMOHO) return; if (rr >= R670) { // upper mantle kmax = kmaxUM; a = aUM; b = bUM; Lmax = LmaxUM; } else { // Lower mantle kmax = kmaxLM; a = aLM; b = bLM; Lmax = LmaxLM; } double *func = new double[kmax+1]; radfun(radfun_t, func, rr, kmax); CreateAB(xr); double *arl, *brl; for (L = 0; L <= Lmax; L++) { arl = Alm[L]; brl = Blm[L]; for (m = 0; m <= L; m++) { for (k = 0; k <= kmax; k++) { arl[m] += func[k] * a[L][m][k]; brl[m] += func[k] * b[L][m][k]; } } } delete [] func; } // Do the radial summation for spherical harmonics model. // The result is in coeffcients Alm and Blm. void CHEB_WHOLE_T::SumRad(double xr) { double rr = xr * RADIUS; int m, L, k; if (rr < RCMB || rr > RMOHO) return; double *func = new double[kmax+1]; radfun(radfun_t, func, rr, kmax); CreateAB(xr); double *arl, *brl; for (L = 0; L <= Lmax; L++) { arl = Alm[L]; brl = Blm[L]; for (m = 0; m <= L; m++) { for (k = 0; k <= kmax; k++) { arl[m] += func[k] * a[L][m][k]; brl[m] += func[k] * b[L][m][k]; } } } delete [] func; } // Do the radial summation for spherical harmonics model. // The result is in coeffcients Alm and Blm. void BSP_WHOLE_T::SumRad(double xr) { double rr = xr * RADIUS; int m, L, k; CreateAB(xr); if (rr < RCMB || rr > RMOHO) return; double depth = RADIUS - rr; double B[4]; k = radialBSpline->BSplKer(B, depth); if (k <= 0 || k >= n_rad) return; int km1 = k - 1; int kp1 = k + 1; int kp2 = k + 2; double B0 = B[0]; double B1 = B[1]; double B2 = B[2]; double B3 = B[3]; double *arl, *brl; double *alm, *blm, **al, **bl; // for efficiency double Almk, Blmk, Almkp1, Blmkp1; for (L = 0; L <= Lmax; L++) { al = a[L]; bl = b[L]; // Summation of radial functions arl = Alm[L]; brl = Blm[L]; for (m = 0; m <= L; m++) { alm = al[m]; blm = bl[m]; Almk = alm[k]; Blmk = blm[k]; Almkp1 = alm[kp1]; Blmkp1 = blm[kp1]; arl[m] = B1 * Almk + B2 * Almkp1; brl[m] = B1 * Blmk + B2 * Blmkp1; if (k != 1) { arl[m] += B0 * alm[km1]; brl[m] += B0 * blm[km1]; } if (kp1 != n_rad) { arl[m] += B3 * alm[kp2]; brl[m] += B3 * blm[kp2]; } } } } /* Function dv2vonly-- * radius (km) * xlat,xlong: co-latitude, co-longitude in radians * output: dVs = velocity anomaly * divided by average velocity (PREM) */ void BSP_WHOLE_T::dv2vonly(double xr, double xlat, double xlong, double *dVs) const { double rr = xr * RADIUS; double *cosm, *sinm, **p; int m, L, k; double phim, plm, cosphi, sinphi; *dVs = 0.0; if (rr < RCMB || rr > RMOHO) return; double depth = RADIUS - rr; double B[4]; k = radialBSpline->BSplKer(B, depth); if (k == 0 || k > n_rad) return; int km1 = k - 1; int kp1 = k + 1; int kp2 = k + 2; double B0 = B[0]; double B1 = B[1]; double B2 = B[2]; double B3 = B[3]; cosm = vector(0, Lmax); sinm = vector(0, Lmax); for (m = 0; m <= Lmax; m++) { phim = m * xlong; cosm[m] = cos(phim); sinm[m] = sin(phim); } p = matrix(0, Lmax, 0, Lmax); Ylm(xlat, Lmax, p); for (L = 0; L <= Lmax; L++) { for (m = 0; m <= L; m++) { plm = p[L][m]; cosphi = cosm[m]; sinphi = sinm[m]; *dVs += B1 * plm * (a[L][m][k] * cosphi + b[L][m][k] * sinphi) + B2 * plm * (a[L][m][kp1] * cosphi + b[L][m][kp1] * sinphi); if (k != 1) *dVs += B0 * plm * (a[L][m][km1] * cosphi + b[L][m][km1] * sinphi); if (kp1 != n_rad) *dVs += B3 * plm * (a[L][m][kp2] * cosphi + b[L][m][kp2] * sinphi); } } free_matrix(p, 0, Lmax, 0, Lmax); free_vector(cosm, 0, Lmax); free_vector(sinm, 0, Lmax); return; } /* Function dv2vonly-- * radius (km) * xlat,xlong: co-latitude, co-longitude in radians * output: dVs = velocity anomaly * divided by average velocity (PREM) */ void BSP_SPLIT_T::dv2vonly(double xr, double xlat, double xlong, double *dVs) const { double rr = xr * RADIUS; double *cosm, *sinm, **p; double *co_rad; double ***a, ***b; RadialBSpline *radialBSpline; int Lmax, n_rad; int m, L, k; double phim, plm, cosphi, sinphi; *dVs = 0.0; if (rr < RCMB || rr > RMOHO) return; if (rr >= R670) { // upper mantle co_rad = co_radUM; n_rad = n_radUM; a = aUM; b = bUM; Lmax = LmaxUM; radialBSpline = radialBSplineUM; } else { // Lower mantle co_rad = co_radLM; n_rad = n_radLM; a = aLM; b = bLM; Lmax = LmaxLM; radialBSpline = radialBSplineLM; } double depth = RADIUS - rr; double B[4]; k = radialBSpline->BSplKer(B, depth); //printf("k=%d B=%f %f %f %f\n", k, B[0], B[1], B[2], B[3]); if (k == 0 || k > n_rad) return; int km1 = k - 1; int kp1 = k + 1; int kp2 = k + 2; double B0 = B[0]; double B1 = B[1]; double B2 = B[2]; double B3 = B[3]; cosm = vector(0, Lmax); sinm = vector(0, Lmax); for (m = 0; m <= Lmax; m++) { phim = m * xlong; cosm[m] = cos(phim); sinm[m] = sin(phim); } p = matrix(0, Lmax, 0, Lmax); Ylm(xlat, Lmax, p); for (L = 0; L <= Lmax; L++) { for (m = 0; m <= L; m++) { plm = p[L][m]; cosphi = cosm[m]; sinphi = sinm[m]; *dVs += B1 * plm * (a[L][m][k] * cosphi + b[L][m][k] * sinphi) + B2 * plm * (a[L][m][kp1] * cosphi + b[L][m][kp1] * sinphi); if (k != 1) *dVs += B0 * plm * (a[L][m][km1] * cosphi + b[L][m][km1] * sinphi); if (kp1 != n_rad) *dVs += B3 * plm * (a[L][m][kp2] * cosphi + b[L][m][kp2] * sinphi); } } free_matrix(p, 0, Lmax, 0, Lmax); free_vector(cosm, 0, Lmax); free_vector(sinm, 0, Lmax); return; } /* Function dv2vonly-- * radius (km) * xlat,xlong: co-latitude, co-longitude in radians * output: dVs =velocity anomaly divided by average velocity (PREM) * */ void CHEB_WHOLE_T::dv2vonly(double xr, double xlat, double xlong, double *dVs) const { double rr = xr * RADIUS; double *cosm, *sinm, **p; double *func; int m, L, k; double phim, plm, cosphi, sinphi; *dVs = 0.0; if (rr < RCMB || rr > RMOHO) return; func = vector(0, kmax); cosm = vector(0, Lmax); sinm = vector(0, Lmax); radfun(radfun_t, func, rr, kmax); for (m = 0; m <= Lmax; m++) { phim = m * xlong; cosm[m] = cos(phim); sinm[m] = sin(phim); } p = matrix(0, Lmax, 0, Lmax); Ylm(xlat, Lmax, p); for (L = 0; L <= Lmax; L++) { for (m = 0; m <= L; m++) { plm = p[L][m]; if (m != 0) { cosphi = cosm[m]; sinphi = sinm[m]; } for (k = 0; k <= kmax; k++) { switch (m) { case 0: *dVs += func[k] * plm * a[L][0][k]; break; default: /* add all spherical harmonics together */ *dVs += func[k] * plm * (a[L][m][k] * cosphi + b[L][m][k] * sinphi); break; } } } } free_matrix(p, 0, Lmax, 0, Lmax); free_vector(func, 0, kmax); free_vector(cosm, 0, Lmax); free_vector(sinm, 0, Lmax); return; } /* Function dv2vonly-- * radius (km) * xlat,xlong: co-latitude, co-longitude in radians * output: dVs =velocity anomaly divided by average velocity (PREM) * */ void CHEB_SPLIT_T::dv2vonly(double xr, double xlat, double xlong, double *dVs) const { double rr = xr * RADIUS; double *cosm, *sinm, **p; double ***a, ***b; double *func; int m, L, k, Lmax, kmax; double phim, plm, cosphi, sinphi; *dVs = 0.0; if (rr < RCMB || rr > RMOHO) return; if (rr >= R670) { // upper mantle kmax = kmaxUM; a = aUM; b = bUM; Lmax = LmaxUM; } else { // Lower mantle kmax = kmaxLM; a = aLM; b = bLM; Lmax = LmaxLM; } func = vector(0, kmax); cosm = vector(0, Lmax); sinm = vector(0, Lmax); radfun(radfun_t, func, rr, kmax); for (m = 0; m <= Lmax; m++) { phim = m * xlong; cosm[m] = cos(phim); sinm[m] = sin(phim); } p = matrix(0, Lmax, 0, Lmax); Ylm(xlat, Lmax, p); for (L = 0; L <= Lmax; L++) { for (m = 0; m <= L; m++) { plm = p[L][m]; if (m != 0) { cosphi = cosm[m]; sinphi = sinm[m]; } for (k = 0; k <= kmax; k++) { switch (m) { case 0: *dVs += func[k] * plm * a[L][0][k]; break; default: /* add all spherical harmonics together */ *dVs += func[k] * plm * (a[L][m][k] * cosphi + b[L][m][k] * sinphi); break; } } } } free_matrix(p, 0, Lmax, 0, Lmax); free_vector(func, 0, kmax); free_vector(cosm, 0, Lmax); free_vector(sinm, 0, Lmax); return; } /* Function dv2v for whole mantle spherical harmonics model (such as * S12WM13 by Su et al. (1994). * Input: radius (km) * xlat,xlong: co-latitude, co-longitude in radians * output: dv =velocity anomaly divided by average velocity (PREM) * dvdx =d(dv)/dx (x=r/6371, normalized radius) * dvdtheta =d(dv)/dtheta * dvdphi =d(dv)/dphi * Warning: hacks ahead. This function is not quite readable. * Efficiency is the number one concern here. This function * consumes most (more than 70%) of the CPU time in ray tracing. */ void CHEB_WHOLE_T::dv2v(double xr, double xlat, double xlong, double *dv, double *dvdx, double *dvdtheta, double *dvdphi) const { double rr = xr * RADIUS; double *cosm, *sinm, **p, **dp; double *func, *dfunc; int m, L, k; double phim, plm, dplm, cosphi, sinphi, sinphim; double cosphim, funck, dfunck, ac, pac, dpac, Aklm, Bklm; double acbs, pacbs, dpacbs, asbc, pasbc; *dv = 0.0; *dvdx = 0.0; *dvdtheta = 0.0; *dvdphi = 0.0; if (rr < RCMB || rr > RMOHO) return; double np_dv = 0.0; // for efficiency double np_dvdx = 0.0; double np_dvdtheta = 0.0; double np_dvdphi = 0.0; func = vector(0, kmax); dfunc = vector(0, kmax); cosm = vector(0, Lmax); sinm = vector(0, Lmax); dradfun(radfun_t, func, dfunc, rr, kmax); // We don't need to calculate m = 0 for (m = 0; m <= Lmax; m++) { phim = m * xlong; cosm[m] = cos(phim); sinm[m] = sin(phim); } p = matrix(0, Lmax, 0, Lmax); dp = matrix(0, Lmax, 0, Lmax); dYlm(xlat, Lmax, p, dp); double *alm, *blm, **al, **bl, *pl, *dpl; // for efficiency for (L = 0; L <= Lmax; L++) { al = a[L]; // for efficiency bl = b[L]; pl = p[L]; dpl = dp[L]; // m = 0 first plm = pl[0]; dplm = dpl[0]; alm = al[0]; // for efficiency for (k = 0; k <= kmax; k++) { funck = func[k]; dfunck = dfunc[k]; ac = alm[k]; pac = plm * ac; dpac = dplm * ac; np_dv += funck * pac; np_dvdx += dfunck * pac; np_dvdtheta += funck * dpac; } // m > 0 for (m = 1; m <= L; m++) { plm = pl[m]; // p[L][m] dplm = dpl[m]; // dp[L][m] cosphi = cosm[m]; sinphi = sinm[m]; sinphim = sinphi * m; cosphim = cosphi * m; alm = al[m]; blm = bl[m]; for (k = 0; k <= kmax; k++) { funck = func[k]; dfunck = dfunc[k]; Aklm = alm[k]; // a[L][m][k] Bklm = blm[k]; // b[L][m][k] acbs = Aklm * cosphi + Bklm * sinphi; pacbs = plm * acbs; dpacbs = dplm * acbs; asbc = -Aklm * sinphim + Bklm * cosphim; pasbc = asbc * plm; np_dv += funck * pacbs; np_dvdx += dfunck * pacbs; np_dvdtheta += funck * dpacbs; np_dvdphi += funck * pasbc; } } } free_matrix(p, 0, Lmax, 0, Lmax); free_matrix(dp, 0, Lmax, 0, Lmax); free_vector(func, 0, kmax); free_vector(dfunc, 0, kmax); free_vector(cosm, 0, Lmax); free_vector(sinm, 0, Lmax); *dv = np_dv; *dvdx = np_dvdx; *dvdtheta = np_dvdtheta; *dvdphi = np_dvdphi; return; } /* Function dv2v for split mantle spherical harmonics model (such as * S12WM13 by Su et al. (1994). * Input: radius (km) * xlat,xlong: co-latitude, co-longitude in radians * output: dv =velocity anomaly divided by average velocity (PREM) * dvdx =d(dv)/dx (x=r/6371, normalized radius) * dvdtheta =d(dv)/dtheta * dvdphi =d(dv)/dphi * Warning: hacks ahead. This function is not quite readable. * Efficiency is the number one concern here. This function * consumes most (more than 70%) of the CPU time in ray tracing. */ void CHEB_SPLIT_T::dv2v(double xr, double xlat, double xlong, double *dv, double *dvdx, double *dvdtheta, double *dvdphi) const { double rr = xr * RADIUS; double *cosm, *sinm, **p, **dp; double *func, *dfunc; int m, L, k, Lmax, kmax; double phim, plm, dplm, cosphi, sinphi, sinphim; double cosphim, funck, dfunck, ac, pac, dpac, Aklm, Bklm; double acbs, pacbs, dpacbs, asbc, pasbc; double ***a, ***b; *dv = 0.0; *dvdx = 0.0; *dvdtheta = 0.0; *dvdphi = 0.0; if (rr < RCMB || rr > RMOHO) return; if (rr >= R670) { // upper mantle kmax = kmaxUM; a = aUM; b = bUM; Lmax = LmaxUM; } else { // Lower mantle kmax = kmaxLM; a = aLM; b = bLM; Lmax = LmaxLM; } double np_dv = 0.0; // for efficiency double np_dvdx = 0.0; double np_dvdtheta = 0.0; double np_dvdphi = 0.0; func = vector(0, kmax); dfunc = vector(0, kmax); cosm = vector(0, Lmax); sinm = vector(0, Lmax); dradfun(radfun_t, func, dfunc, rr, kmax); // We don't need to calculate m = 0 for (m = 0; m <= Lmax; m++) { phim = m * xlong; cosm[m] = cos(phim); sinm[m] = sin(phim); } p = matrix(0, Lmax, 0, Lmax); dp = matrix(0, Lmax, 0, Lmax); dYlm(xlat, Lmax, p, dp); double *alm, *blm, **al, **bl, *pl, *dpl; // for efficiency for (L = 0; L <= Lmax; L++) { al = a[L]; // for efficiency bl = b[L]; pl = p[L]; dpl = dp[L]; // m = 0 first plm = pl[0]; dplm = dpl[0]; alm = al[0]; // for efficiency for (k = 0; k <= kmax; k++) { funck = func[k]; dfunck = dfunc[k]; ac = alm[k]; pac = plm * ac; dpac = dplm * ac; np_dv += funck * pac; np_dvdx += dfunck * pac; np_dvdtheta += funck * dpac; } // m > 0 for (m = 1; m <= L; m++) { plm = pl[m]; // p[L][m] dplm = dpl[m]; // dp[L][m] cosphi = cosm[m]; sinphi = sinm[m]; sinphim = sinphi * m; cosphim = cosphi * m; alm = al[m]; blm = bl[m]; for (k = 0; k <= kmax; k++) { funck = func[k]; dfunck = dfunc[k]; Aklm = alm[k]; // a[L][m][k] Bklm = blm[k]; // b[L][m][k] acbs = Aklm * cosphi + Bklm * sinphi; pacbs = plm * acbs; dpacbs = dplm * acbs; asbc = -Aklm * sinphim + Bklm * cosphim; pasbc = asbc * plm; np_dv += funck * pacbs; np_dvdx += dfunck * pacbs; np_dvdtheta += funck * dpacbs; np_dvdphi += funck * pasbc; } } } free_matrix(p, 0, Lmax, 0, Lmax); free_matrix(dp, 0, Lmax, 0, Lmax); free_vector(func, 0, kmax); free_vector(dfunc, 0, kmax); free_vector(cosm, 0, Lmax); free_vector(sinm, 0, Lmax); *dv = np_dv; *dvdx = np_dvdx; *dvdtheta = np_dvdtheta; *dvdphi = np_dvdphi; return; } /* Function dv2v for whole mantle spherical harmonics model (such as * S12WM13 by Su et al. (1994). * Input: radius (km) * xlat,xlong: co-latitude, co-longitude in radians * output: dv =velocity anomaly divided by average velocity (PREM) * dvdx =d(dv)/dx (x=r/6371, normalized radius) * dvdtheta =d(dv)/dtheta * dvdphi =d(dv)/dphi * Warning: hacks ahead. This function is not quite readable. * Efficiency is the number one concern here. This function * consumes most (more than 70%) of the CPU time in ray tracing. */ void BSP_WHOLE_T::dv2v(double xr, double xlat, double xlong, double *dv, double *dvdx, double *dvdtheta, double *dvdphi) const { double rr = xr * RADIUS; double *cosm, *sinm, **p, **dp; int m, L; double phim, plm, dplm, cosphi, sinphi, sinphim; double cosphim; *dv = 0.0; *dvdx = 0.0; *dvdtheta = 0.0; *dvdphi = 0.0; if (rr < RCMB || rr > RMOHO) return; double np_dv = 0.0; // for efficiency double np_dvdx = 0.0; double np_dvdtheta = 0.0; double np_dvdphi = 0.0; double B[4], DB[4]; int k = radialBSpline->DBBSplKer(B, DB, (RADIUS - rr)); if (k <= 0 || k >= n_rad) return; int km1 = k - 1; int kp1 = k + 1; int kp2 = k + 2; double B0 = B[0]; double B1 = B[1]; double B2 = B[2]; double B3 = B[3]; double DB0 = DB[0]; double DB1 = DB[1]; double DB2 = DB[2]; double DB3 = DB[3]; cosm = vector(0, Lmax); sinm = vector(0, Lmax); for (m = 0; m <= Lmax; m++) { phim = m * xlong; cosm[m] = cos(phim); sinm[m] = sin(phim); } p = matrix(0, Lmax, 0, Lmax); dp = matrix(0, Lmax, 0, Lmax); dYlm(xlat, Lmax, p, dp); double *alm, *blm, **al, **bl, *pl, *dpl; // for efficiency double Almk, Blmk, Almkp1, Blmkp1, acbs, acbsp1; double Almkm1, Blmkm1, acbsm1; double Almkp2, Blmkp2, acbsp2; for (L = 0; L <= Lmax; L++) { pl = p[L]; dpl = dp[L]; al = a[L]; bl = b[L]; for (m = 0; m <= L; m++) { alm = al[m]; blm = bl[m]; plm = pl[m]; dplm = dpl[m]; cosphi = cosm[m]; sinphi = sinm[m]; sinphim = sinphi * m; cosphim = cosphi * m; Almk = alm[k]; Blmk = blm[k]; Almkp1 = alm[kp1]; Blmkp1 = blm[kp1]; acbs = Almk * cosphi + Blmk * sinphi; acbsp1 = (Almkp1 * cosphi + Blmkp1 * sinphi); np_dv += B1 * plm * acbs + B2 * plm * acbsp1; np_dvdx += DB1 * plm * acbs + DB2 * plm * acbsp1; np_dvdtheta += B1 * dplm * acbs + B2 * dplm * acbsp1; np_dvdphi += B1 * plm * (-Almk * sinphim + Blmk * cosphim) + B2 * plm * (-Almkp1 * sinphim + Blmkp1 * cosphim); if (k != 1) { Almkm1 = alm[km1]; Blmkm1 = blm[km1]; acbsm1 = (Almkm1 * cosphi + Blmkm1 * sinphi); np_dv += B0 * plm * acbsm1; np_dvdx += DB0 * plm * acbsm1; np_dvdtheta += B0 * dplm * acbsm1; np_dvdphi += B0 * plm * (-Almkm1 * sinphim + Blmkm1 * cosphim); } if (kp1 != n_rad) { Almkp2 = alm[kp2]; Blmkp2 = blm[kp2]; acbsp2 = (Almkp2 * cosphi + Blmkp2 * sinphi); np_dv += B3 * plm * acbsp2; np_dvdx += DB3 * plm * acbsp2; np_dvdtheta += B3 * dplm * acbsp2; np_dvdphi += B3 * plm * (-Almkp2 * sinphim + Blmkp2 * cosphim); } } } free_matrix(p, 0, Lmax, 0, Lmax); free_matrix(dp, 0, Lmax, 0, Lmax); free_vector(cosm, 0, Lmax); free_vector(sinm, 0, Lmax); *dv = np_dv; *dvdx = -np_dvdx * RADIUS; *dvdtheta = np_dvdtheta; *dvdphi = np_dvdphi; return; } /* Function dv2v for split mantle (radially B-spline) * spherical harmonics model. * Input: radius (km) * xlat,xlong: co-latitude, co-longitude in radians * output: dv =velocity anomaly divided by average velocity (PREM) * dvdx =d(dv)/dx (x=r/6371, normalized radius) * dvdtheta =d(dv)/dtheta * dvdphi =d(dv)/dphi * Warning: hacks ahead. This function is not quite readable. * Efficiency is the number one concern here. This function * consumes most (more than 70%) of the CPU time in ray tracing. */ void BSP_SPLIT_T::dv2v(double xr, double xlat, double xlong, double *dv, double *dvdx, double *dvdtheta, double *dvdphi) const { double rr = xr * RADIUS; double *cosm, *sinm, **p, **dp; double *co_rad; double ***a, ***b; RadialBSpline *radialBSpline; int Lmax, n_rad; int m, L; double phim, plm, dplm, cosphi, sinphi, sinphim; double cosphim; *dv = 0.0; *dvdx = 0.0; *dvdtheta = 0.0; *dvdphi = 0.0; if (rr < RCMB || rr > RMOHO) return; if (rr >= R670) { // upper mantle co_rad = co_radUM; n_rad = n_radUM; a = aUM; b = bUM; Lmax = LmaxUM; radialBSpline = radialBSplineUM; } else { // Lower mantle co_rad = co_radLM; n_rad = n_radLM; a = aLM; b = bLM; Lmax = LmaxLM; radialBSpline = radialBSplineLM; } double np_dv = 0.0; // for efficiency double np_dvdx = 0.0; double np_dvdtheta = 0.0; double np_dvdphi = 0.0; double B[4], DB[4]; int k = radialBSpline->DBBSplKer(B, DB, (RADIUS - rr)); if (k <= 0 || k >= n_rad) return; int km1 = k - 1; int kp1 = k + 1; int kp2 = k + 2; double B0 = B[0]; double B1 = B[1]; double B2 = B[2]; double B3 = B[3]; double DB0 = DB[0]; double DB1 = DB[1]; double DB2 = DB[2]; double DB3 = DB[3]; cosm = vector(0, Lmax); sinm = vector(0, Lmax); for (m = 0; m <= Lmax; m++) { phim = m * xlong; cosm[m] = cos(phim); sinm[m] = sin(phim); } p = matrix(0, Lmax, 0, Lmax); dp = matrix(0, Lmax, 0, Lmax); dYlm(xlat, Lmax, p, dp); double *alm, *blm, **al, **bl, *pl, *dpl; // for efficiency double Almk, Blmk, Almkp1, Blmkp1, acbs, acbsp1; double Almkm1, Blmkm1, acbsm1; double Almkp2, Blmkp2, acbsp2; for (L = 0; L <= Lmax; L++) { pl = p[L]; dpl = dp[L]; al = a[L]; bl = b[L]; for (m = 0; m <= L; m++) { alm = al[m]; blm = bl[m]; plm = pl[m]; dplm = dpl[m]; cosphi = cosm[m]; sinphi = sinm[m]; sinphim = sinphi * m; cosphim = cosphi * m; Almk = alm[k]; Blmk = blm[k]; Almkp1 = alm[kp1]; Blmkp1 = blm[kp1]; acbs = Almk * cosphi + Blmk * sinphi; acbsp1 = (Almkp1 * cosphi + Blmkp1 * sinphi); np_dv += B1 * plm * acbs + B2 * plm * acbsp1; np_dvdx += DB1 * plm * acbs + DB2 * plm * acbsp1; np_dvdtheta += B1 * dplm * acbs + B2 * dplm * acbsp1; np_dvdphi += B1 * plm * (-Almk * sinphim + Blmk * cosphim) + B2 * plm * (-Almkp1 * sinphim + Blmkp1 * cosphim); if (k != 1) { Almkm1 = alm[km1]; Blmkm1 = blm[km1]; acbsm1 = (Almkm1 * cosphi + Blmkm1 * sinphi); np_dv += B0 * plm * acbsm1; np_dvdx += DB0 * plm * acbsm1; np_dvdtheta += B0 * dplm * acbsm1; np_dvdphi += B0 * plm * (-Almkm1 * sinphim + Blmkm1 * cosphim); } if (kp1 != n_rad) { Almkp2 = alm[kp2]; Blmkp2 = blm[kp2]; acbsp2 = (Almkp2 * cosphi + Blmkp2 * sinphi); np_dv += B3 * plm * acbsp2; np_dvdx += DB3 * plm * acbsp2; np_dvdtheta += B3 * dplm * acbsp2; np_dvdphi += B3 * plm * (-Almkp2 * sinphim + Blmkp2 * cosphim); } } } free_matrix(p, 0, Lmax, 0, Lmax); free_matrix(dp, 0, Lmax, 0, Lmax); free_vector(cosm, 0, Lmax); free_vector(sinm, 0, Lmax); *dv = np_dv; *dvdx = -np_dvdx * RADIUS; *dvdtheta = np_dvdtheta; *dvdphi = np_dvdphi; return; } // Convert a 1-D vector into 3-D coefficients a[][][] and b[][][] void BSP_SPLIT_T::ConvCoeff(double x[]) { int k, L, m; int index = 0; for (k = 1; k <= n_radUM; k++) { for (L = 0; L <= LmaxUM; L++) { aUM[L][0][k] = x[++index]; for (m = 1; m <= L; m++) { aUM[L][m][k] = x[++index]; bUM[L][m][k] = x[++index]; } } } for (k = 1; k <= n_radLM; k++) { for (L = 0; L <= LmaxLM; L++) { aLM[L][0][k] = x[++index]; for (m = 1; m <= L; m++) { aLM[L][m][k] = x[++index]; bLM[L][m][k] = x[++index]; } } } return; } // Convert a 1-D vector into 3-D coefficients a[][][] and b[][][] void BSP_WHOLE_T::ConvCoeff(double x[]) { int index = 0; for (int k = 1; k <= n_rad; k++) { for (int L = 0; L <= Lmax; L++) { a[L][0][k] = x[++index]; for (int m = 1; m <= L; m++) { a[L][m][k] = x[++index]; b[L][m][k] = x[++index]; } } } return; } // read the model coefficents into a file. void BSP_WHOLE_T::read_model3D(void) { char DegreeMask[MAXSTR], RadialMask[MAXSTR]; int n_rad0; int rty; FILE *fin = fopen(name, "r"); if (fin == NULL) { fprintf(stderr, "cannot open file %s\n", name); exit(1); } fscanf(fin, "%d %s %d %s %d", &Lmax, DegreeMask, &n_rad0, RadialMask, &rty); n_rad = n_rad0 - 2; // minus crust and 670 coefficients printf("%d %s %d %s %d\n", Lmax, DegreeMask, n_rad, RadialMask, rty); if (rty != 4) { fprintf(stderr, "Model type is not what expected!\n"); exitLog("Model type is not what expected!\n"); } if (n_rad <= 0 || Lmax <= 0) { fprintf(stderr, "Model type is not what expected!\n"); exitLog("Model type is not what expected!\n"); } int k, L; double dummy; // skip crust and 670 discontinuity coefficients for (k = 0; k < 2; k++) { if (RadialMask[k] =='0') continue; if (k == 0) printf("Skip crustal coefficients\n"); else printf("Skip 670km discontinuity coefficients\n"); for (L = 0; L <= Lmax; L++) { if (DegreeMask[L] == '0') continue; switch (L) { case 0: fscanf(fin, "%le", &dummy); break; default: for (int m = 0; m <= L; m++) { switch (m) { case 0: fscanf(fin, "%le", &dummy); break; default: fscanf(fin, "%le %le", &dummy, &dummy); break; } } } } } double **ak = matrix(0, Lmax, 0, Lmax); double **bk = matrix(0, Lmax, 0, Lmax); double **akm = matrix(0, Lmax, 0, Lmax); double **bkm = matrix(0, Lmax, 0, Lmax); CreateModel3D(); for (k = 1; k <= n_rad; k++) { if (RadialMask[k+2] == '0') { printf("Radial k = %d is skipped\n", k); continue; } if (RadialMask[k+2] == '2') { printf("Discontinuous model is not expected! %d\n", k); exitLog("Discontinuous model is not expected!\n"); } for (L = 0; L <= Lmax; L++) { if (DegreeMask[L] == '0') { printf("Degree L = %d is skipped\n", L); continue; } fscanf(fin, "%le", &ak[L][0]); for (int m = 1; m <= L; m++) fscanf(fin, "%le %le", &ak[L][m], &bk[L][m]); a0[L][0][k] = ak[L][0]; if (L > 0) { for (m = 1; m <= L; m++) { a0[L][m][k] = ak[L][m]; b0[L][m][k] = bk[L][m]; } } } if (alpha != 0.0 || beta != 0.0 || gamma != 0.0) { rotate_velocity(alpha, beta, gamma, Lmax, ak, bk, akm, bkm); for (L = 0; L <= Lmax; L++) { a[L][0][k] = akm[L][0]; if (L > 0) { for (int m = 1; m <= L; m++) { a[L][m][k] = akm[L][m]; b[L][m][k] = bkm[L][m]; } } } } else { for (L = 0; L <= Lmax; L++) { a[L][0][k] = ak[L][0]; if (L > 0) { for (int m = 1; m <= L; m++) { a[L][m][k] = ak[L][m]; b[L][m][k] = bk[L][m]; } } } } } for (k = 1; k <= n_rad; k++) { fscanf(fin, "%lf", &co_rad[k]); if (co_rad[k] > (RADIUS-RCMB) || co_rad[k] < (RADIUS-RMOHO)) { fprintf(stderr, "Wrong with your 3-D model coefficients file\n"); fprintf(stderr, "Radial knot value co_rad[%d] = %f\n", k, co_rad[k]); char msg[MAXSTR]; sprintf(msg, "Wrong with your 3-D model coefficients file, co_rad[%d] = %f\n", k, co_rad[k]); exitLog(msg); } } free_matrix(ak, 0, Lmax, 0, Lmax); free_matrix(bk, 0, Lmax, 0, Lmax); free_matrix(akm, 0, Lmax, 0, Lmax); free_matrix(bkm, 0, Lmax, 0, Lmax); fclose(fin); co_rad[-1] = co_rad[0] = co_rad[1]; co_rad[n_rad+2] = co_rad[n_rad+1] = co_rad[n_rad]; int n_radp1 = n_rad + 1; // set unneccessary coefficients to zero for (L = 0; L <= Lmax; L++) { for (int n = 0; n <= n_rad; n++) { b[L][0][n] = 0.0; b0[L][0][n] = 0.0; } for (int m = 0; m <= L; m++) { a[L][m][0] = 0.0; b[L][m][0] = 0.0; a[L][m][n_radp1] = 0.0; b[L][m][n_radp1] = 0.0; a0[L][m][0] = 0.0; b0[L][m][0] = 0.0; a0[L][m][n_radp1] = 0.0; b0[L][m][n_radp1] = 0.0; } } return; } // read the model coefficents into a file. void BSP_SPLIT_T::read_model3D(void) { char DegreeMask[MAXSTR], RadialMask[MAXSTR]; int n_rad0; int rty; FILE *fin = fopen(name, "r"); if (fin == NULL) { fprintf(stderr, "cannot open file %s\n", name); exit(1); } fscanf(fin, "%d %s %d %s %d", &LmaxUM, DegreeMask, &n_rad0, RadialMask, &rty); LmaxLM = LmaxUM; if (rty != 4) { fprintf(stderr, "Model type is not what expected!\n"); exitLog("Model type is not what expected!\n"); } if (n_rad0 <= 0 || LmaxUM <= 0 || LmaxLM <= 0) { fprintf(stderr, "Model type is not what expected!\n"); exitLog("Model type is not what expected!\n"); } n_radUM = 0; n_radLM = 0; for (int i = 0; i < n_rad0; i++) { if (RadialMask[i] == '1' || RadialMask[i] == '0') { n_radUM++; } else if (RadialMask[i] == '2') { // lower mantle flag n_radLM++; } else if (RadialMask[i] != '0') { // Model type verification fprintf(stderr, "Either radial mask or the model type"); fprintf(stderr, " is not what expected!\n"); exitLog("Either radial mask or the model type is not what expected!\n"); } } n_radUM -= 2; // minus the radial mask for crust and 670 // Error checking if (n_radUM < 0 || n_radLM < 0) { fprintf(stderr, " n_radUM = %d n_radLM = %d\n", n_radUM, n_radLM); exit(1); } int k, L; double dummy; // skip crust and 670 discontinuity coefficients for (k = 0; k < 2; k++) { if (RadialMask[k] =='0') continue; if (k == 0) printf("Skip crustal coefficients\n"); else printf("Skip 670km discontinuity coefficients\n"); for (L = 0; L <= LmaxUM; L++) { if (DegreeMask[L] == '0') continue; switch (L) { case 0: fscanf(fin, "%le", &dummy); break; default: for (int m = 0; m <= L; m++) { switch (m) { case 0: fscanf(fin, "%le", &dummy); break; default: fscanf(fin, "%le %le", &dummy, &dummy); break; } } } } } double **ak = matrix(0, LmaxUM, 0, LmaxUM); double **bk = matrix(0, LmaxUM, 0, LmaxUM); double **akm = matrix(0, LmaxUM, 0, LmaxUM); double **bkm = matrix(0, LmaxUM, 0, LmaxUM); CreateModel3D(); for (k = 1; k <= n_radUM; k++) { if (RadialMask[k+1] == '0') continue; for (L = 0; L <= LmaxUM; L++) { if (DegreeMask[L] == '0') continue; fscanf(fin, "%le", &ak[L][0]); for (int m = 1; m <= L; m++) fscanf(fin, "%le %le", &ak[L][m], &bk[L][m]); a0UM[L][0][k] = ak[L][0]; if (L > 0) { for (m = 1; m <= L; m++) { a0UM[L][m][k] = ak[L][m]; b0UM[L][m][k] = bk[L][m]; } } } if (alpha != 0.0 || beta != 0.0 || gamma != 0.0) { rotate_velocity(alpha, beta, gamma, LmaxUM, ak, bk, akm, bkm); for (L = 0; L <= LmaxUM; L++) { aUM[L][0][k] = akm[L][0]; if (L > 0) { for (int m = 1; m <= L; m++) { aUM[L][m][k] = akm[L][m]; bUM[L][m][k] = bkm[L][m]; } } } } else { for (L = 0; L <= LmaxUM; L++) { aUM[L][0][k] = ak[L][0]; if (L > 0) { for (int m = 1; m <= L; m++) { aUM[L][m][k] = ak[L][m]; bUM[L][m][k] = bk[L][m]; } } } } } free_matrix(ak, 0, LmaxUM, 0, LmaxUM); free_matrix(bk, 0, LmaxUM, 0, LmaxUM); free_matrix(akm, 0, LmaxUM, 0, LmaxUM); free_matrix(bkm, 0, LmaxUM, 0, LmaxUM); int n_radp1 = n_radUM + 1; // set unneccessary coefficients to zero for (L = 0; L <= LmaxUM; L++) { for (int n = 0; n <= n_radUM; n++) { bUM[L][0][n] = 0.0; b0UM[L][0][n] = 0.0; } for (int m = 0; m <= L; m++) { aUM[L][m][0] = 0.0; bUM[L][m][0] = 0.0; aUM[L][m][n_radp1] = 0.0; bUM[L][m][n_radp1] = 0.0; a0UM[L][m][0] = 0.0; b0UM[L][m][0] = 0.0; a0UM[L][m][n_radp1] = 0.0; b0UM[L][m][n_radp1] = 0.0; } } // Now lower mantle ak = matrix(0, LmaxLM, 0, LmaxLM); bk = matrix(0, LmaxLM, 0, LmaxLM); akm = matrix(0, LmaxLM, 0, LmaxLM); bkm = matrix(0, LmaxLM, 0, LmaxLM); for (k = 1; k <= n_radLM; k++) { for (L = 0; L <= LmaxLM; L++) { fscanf(fin, "%le", &ak[L][0]); for (int m = 1; m <= L; m++) fscanf(fin, "%le %le", &ak[L][m], &bk[L][m]); a0LM[L][0][k] = ak[L][0]; if (L > 0) { for (m = 1; m <= L; m++) { a0LM[L][m][k] = ak[L][m]; b0LM[L][m][k] = bk[L][m]; } } } if (alpha != 0.0 || beta != 0.0 || gamma != 0.0) { rotate_velocity(alpha, beta, gamma, LmaxLM, ak, bk, akm, bkm); for (L = 0; L <= LmaxLM; L++) { aLM[L][0][k] = akm[L][0]; if (L > 0) { for (int m = 1; m <= L; m++) { aLM[L][m][k] = akm[L][m]; bLM[L][m][k] = bkm[L][m]; } } } } else { for (L = 0; L <= LmaxLM; L++) { aLM[L][0][k] = ak[L][0]; if (L > 0) { for (int m = 1; m <= L; m++) { aLM[L][m][k] = ak[L][m]; bLM[L][m][k] = bk[L][m]; } } } } } for (k = 1; k <= n_radUM; k++) { fscanf(fin, "%lf", &co_radUM[k]); if (co_radUM[k] > (RADIUS-RCMB) || co_radUM[k] < (RADIUS-RMOHO)) { fprintf(stderr, "Wrong with your 3-D model coefficients file\n"); fprintf(stderr, "Radial knot value co_radUM[%d] = %f\n", k, co_radUM[k]); char msg[MAXSTR]; sprintf(msg, "Wrong with your 3-D model coefficients file, co_radUM[%d] = %f\n", k, co_radUM[k]); exitLog(msg); } } for (k = 1; k <= n_radLM; k++) { fscanf(fin, "%lf", &co_radLM[k]); if (co_radLM[k] > (RADIUS-RCMB) || co_radLM[k] < (RADIUS-R670)) { fprintf(stderr, "Wrong with your 3-D model coefficients file\n"); fprintf(stderr, "Radial knot value co_radLM[%d] = %f\n", k, co_radLM[k]); char msg[MAXSTR]; sprintf(msg, "Wrong with your 3-D model coefficients file, co_radLM[%d] = %f\n", k, co_radLM[k]); exitLog(msg); } } if (co_radLM[1] != (RADIUS-R670)) { fprintf(stderr, "Wrong with your 3-D model coefficients file\n"); fprintf(stderr, "Radial knot value R670 = %f\n", co_radLM[1]); char msg[MAXSTR]; sprintf(msg, "Wrong with your 3-D model coefficients file, R670 = %f\n", co_radLM[1]); exitLog(msg); } free_matrix(ak, 0, LmaxLM, 0, LmaxLM); free_matrix(bk, 0, LmaxLM, 0, LmaxLM); free_matrix(akm, 0, LmaxLM, 0, LmaxLM); free_matrix(bkm, 0, LmaxLM, 0, LmaxLM); fclose(fin); co_radUM[0] = co_radUM[1]; co_radUM[-1] = co_radUM[0]; co_radUM[n_radUM+1] = co_radUM[n_radUM]; co_radUM[n_radUM+2] = co_radUM[n_radUM+1]; co_radLM[-1] = co_radLM[0] = co_radLM[1]; co_radLM[n_radLM+2] = co_radLM[n_radLM+1] = co_radLM[n_radLM]; n_radp1 = n_radLM + 1; // set unneccessary coefficients to zero for (L = 0; L <= LmaxLM; L++) { for (int n = 0; n <= n_radLM; n++) { bLM[L][0][n] = 0.0; b0LM[L][0][n] = 0.0; } for (int m = 0; m <= L; m++) { aLM[L][m][0] = 0.0; bLM[L][m][0] = 0.0; aLM[L][m][n_radp1] = 0.0; bLM[L][m][n_radp1] = 0.0; a0LM[L][m][0] = 0.0; b0LM[L][m][0] = 0.0; a0LM[L][m][n_radp1] = 0.0; b0LM[L][m][n_radp1] = 0.0; } } return; } // read the model coefficents into a file. void CHEB_SPLIT_T::read_model3D(void) { char DegreeMask[MAXSTR], RadialMask[MAXSTR]; int kmax0; int irad; FILE *fin = fopen(name, "r"); if (fin == NULL) { fprintf(stderr, "cannot open file %s\n", name); exit(1); } fscanf(fin, "%d %s %d %s %d", &LmaxUM, DegreeMask, &kmax0, RadialMask, &irad); LmaxLM = LmaxUM; if (irad != 2) { fprintf(stderr, "Model type is not what expected!\n"); exitLog("Model type is not what expected!\n"); } if (kmax0 <= 0 || LmaxUM <= 0 || LmaxLM <= 0) { fprintf(stderr, "Model type is not what expected!\n"); exitLog("Model type is not what expected!\n"); } if (irad == 4) { printf("B-Spline is used in radial function.\n"); radfun_t = B_SPLINE; } else if (irad == 2) { printf("ChebyChev Polynomial is used in radial function.\n"); radfun_t = CHEB_WHOLE; } else { fprintf(stderr, "radial function type is unknown!\n"); fprintf(stderr, "Or model type is what expected!\n"); exitLog("radial function type is unknown! Or model type is what expected!\n"); } kmaxUM = 0; kmaxLM = 0; for (int i = 0; i < kmax0; i++) { if (RadialMask[i] == '1' || RadialMask[i] == '0') { kmaxUM++; } else if (RadialMask[i] == '2') { // lower mantle flag kmaxLM++; } else if (RadialMask[i] != '0') { // Model type verification fprintf(stderr, "Either radial mask or the model type"); fprintf(stderr, " is not what expected!\n"); exitLog("Either radial mask or the model type is not what expected!\n"); } } kmaxUM -= 2; // minus the radial mask for crust and 670 // Error checking if (kmaxUM <= 0 || kmaxLM <= 0) { fprintf(stderr, " kmaxUM = %d kmaxLM = %d\n", kmaxUM, kmaxLM); fprintf(stderr, "Model type is not what I expected\n"); fprintf(stderr, " I expected a split earth model\n"); exit(1); } printf("Highest order of Chebchev polynomial for upper mantle: %d\n", kmaxUM); printf("Highest order of Chebchev polynomial for lower mantle: %d\n", kmaxLM); int k, L; double dummy; // skip crust and 670 discontinuity coefficients for (k = 0; k < 2; k++) { if (RadialMask[k] =='0') continue; if (k == 0) printf("Skip crustal coefficients\n"); else printf("Skip 670km discontinuity coefficients\n"); for (L = 0; L <= LmaxUM; L++) { if (DegreeMask[L] == '0') continue; switch (L) { case 0: fscanf(fin, "%le", &dummy); break; default: for (int m = 0; m <= L; m++) { switch (m) { case 0: fscanf(fin, "%le", &dummy); break; default: fscanf(fin, "%le %le", &dummy, &dummy); break; } } } } } double **ak = matrix(0, LmaxUM, 0, LmaxUM); double **bk = matrix(0, LmaxUM, 0, LmaxUM); double **akm = matrix(0, LmaxUM, 0, LmaxUM); double **bkm = matrix(0, LmaxUM, 0, LmaxUM); CreateModel3D(); for (k = 1; k <= kmaxUM; k++) { if (RadialMask[k+1] == '0') continue; for (L = 0; L <= LmaxUM; L++) { if (DegreeMask[L] == '0') continue; fscanf(fin, "%le", &ak[L][0]); for (int m = 1; m <= L; m++) fscanf(fin, "%le %le", &ak[L][m], &bk[L][m]); a0UM[L][0][k] = ak[L][0]; if (L > 0) { for (m = 1; m <= L; m++) { a0UM[L][m][k] = ak[L][m]; b0UM[L][m][k] = bk[L][m]; } } } if (alpha != 0.0 || beta != 0.0 || gamma != 0.0) { rotate_velocity(alpha, beta, gamma, LmaxUM, ak, bk, akm, bkm); for (L = 0; L <= LmaxUM; L++) { aUM[L][0][k] = akm[L][0]; if (L > 0) { for (int m = 1; m <= L; m++) { aUM[L][m][k] = akm[L][m]; bUM[L][m][k] = bkm[L][m]; } } } } else { for (L = 0; L <= LmaxUM; L++) { aUM[L][0][k] = ak[L][0]; if (L > 0) { for (int m = 1; m <= L; m++) { aUM[L][m][k] = ak[L][m]; bUM[L][m][k] = bk[L][m]; } } } } } free_matrix(ak, 0, LmaxUM, 0, LmaxUM); free_matrix(bk, 0, LmaxUM, 0, LmaxUM); free_matrix(akm, 0, LmaxUM, 0, LmaxUM); free_matrix(bkm, 0, LmaxUM, 0, LmaxUM); int kmaxp1 = kmaxUM + 1; // set unneccessary coefficients to zero for (L = 0; L <= LmaxUM; L++) { for (int n = 0; n <= kmaxUM; n++) { bUM[L][0][n] = 0.0; b0UM[L][0][n] = 0.0; } for (int m = 0; m <= L; m++) { aUM[L][m][0] = 0.0; bUM[L][m][0] = 0.0; aUM[L][m][kmaxp1] = 0.0; bUM[L][m][kmaxp1] = 0.0; a0UM[L][m][0] = 0.0; b0UM[L][m][0] = 0.0; a0UM[L][m][kmaxp1] = 0.0; b0UM[L][m][kmaxp1] = 0.0; } } // Now lower mantle ak = matrix(0, LmaxLM, 0, LmaxLM); bk = matrix(0, LmaxLM, 0, LmaxLM); akm = matrix(0, LmaxLM, 0, LmaxLM); bkm = matrix(0, LmaxLM, 0, LmaxLM); for (k = 1; k <= kmaxLM; k++) { for (L = 0; L <= LmaxLM; L++) { fscanf(fin, "%le", &ak[L][0]); for (int m = 1; m <= L; m++) fscanf(fin, "%le %le", &ak[L][m], &bk[L][m]); a0LM[L][0][k] = ak[L][0]; if (L > 0) { for (m = 1; m <= L; m++) { a0LM[L][m][k] = ak[L][m]; b0LM[L][m][k] = bk[L][m]; } } } if (alpha != 0.0 || beta != 0.0 || gamma != 0.0) { rotate_velocity(alpha, beta, gamma, LmaxLM, ak, bk, akm, bkm); for (L = 0; L <= LmaxLM; L++) { aLM[L][0][k] = akm[L][0]; if (L > 0) { for (int m = 1; m <= L; m++) { aLM[L][m][k] = akm[L][m]; bLM[L][m][k] = bkm[L][m]; } } } } else { for (L = 0; L <= LmaxLM; L++) { aLM[L][0][k] = ak[L][0]; if (L > 0) { for (int m = 1; m <= L; m++) { aLM[L][m][k] = ak[L][m]; bLM[L][m][k] = bk[L][m]; } } } } } free_matrix(ak, 0, LmaxLM, 0, LmaxLM); free_matrix(bk, 0, LmaxLM, 0, LmaxLM); free_matrix(akm, 0, LmaxLM, 0, LmaxLM); free_matrix(bkm, 0, LmaxLM, 0, LmaxLM); fclose(fin); kmaxp1 = kmaxLM + 1; // set unneccessary coefficients to zero for (L = 0; L <= LmaxLM; L++) { for (int n = 0; n <= kmaxLM; n++) { bLM[L][0][n] = 0.0; b0LM[L][0][n] = 0.0; } for (int m = 0; m <= L; m++) { aLM[L][m][0] = 0.0; bLM[L][m][0] = 0.0; aLM[L][m][kmaxp1] = 0.0; bLM[L][m][kmaxp1] = 0.0; a0LM[L][m][0] = 0.0; b0LM[L][m][0] = 0.0; a0LM[L][m][kmaxp1] = 0.0; b0LM[L][m][kmaxp1] = 0.0; } } return; } // read the model coefficents into a file. // This model comes directly from decinv*.f. // Its header information is not correct, and it does not // contain any radial values for split B-spline model. void BSP_SPLIT_T::read_NO_RADIUS_model3D(void) { char DegreeMask[MAXSTR], RadialMask[MAXSTR]; char Line[MAXSTR]; int n_rad0; FILE *fin = fopen(name, "r"); if (fin == NULL) { fprintf(stderr, "cannot open file %s\n", name); exit(1); } if (fgets(Line, MAXSTR, fin) == NULL) { fprintf(stderr, "end of file\n"); exitLog("end of file\n"); } fprintf(stderr, "first line: %s", Line); sscanf(Line, "%d %s %d %s", &LmaxUM, DegreeMask, &n_rad0, RadialMask); for (int i = 0; i < MAXSTR; i++) { if (RadialMask[i] == '\0') break; if (RadialMask[i] == '\n') { RadialMask[i] = '\0'; break; } } printf("%d %s %d %s\n", LmaxUM, DegreeMask, n_rad0, RadialMask); LmaxLM = LmaxUM; if (n_rad0 <= 0 || LmaxUM <= 0 || LmaxLM <= 0) { fprintf(stderr, "%d %d %d\n", n_rad0, LmaxUM, LmaxLM); fprintf(stderr, "Model type is not what expected!\n"); exitLog("Model type is not what expected!\n"); } n_radUM = 6; n_radLM = 8; if (n_rad0 != (n_radUM+n_radLM)) { fprintf(stderr, "inconsistent model\n"); exitLog("inconsistent model\n"); } printf("n_radUM = %d, n_radLM = %d, Lmax = %d\n", n_radUM, n_radLM, LmaxUM); int k, L; double **ak = matrix(0, LmaxUM, 0, LmaxUM); double **bk = matrix(0, LmaxUM, 0, LmaxUM); double **akm = matrix(0, LmaxUM, 0, LmaxUM); double **bkm = matrix(0, LmaxUM, 0, LmaxUM); CreateModel3D(); for (k = 1; k <= n_radUM; k++) { if (RadialMask[k+1] == '0') continue; for (L = 0; L <= LmaxUM; L++) { if (DegreeMask[L] == '0') continue; fscanf(fin, "%le", &ak[L][0]); for (int m = 1; m <= L; m++) fscanf(fin, "%le %le", &ak[L][m], &bk[L][m]); a0UM[L][0][k] = ak[L][0]; if (L > 0) { for (m = 1; m <= L; m++) { a0UM[L][m][k] = ak[L][m]; b0UM[L][m][k] = bk[L][m]; } } } if (alpha != 0.0 || beta != 0.0 || gamma != 0.0) { rotate_velocity(alpha, beta, gamma, LmaxUM, ak, bk, akm, bkm); for (L = 0; L <= LmaxUM; L++) { aUM[L][0][k] = akm[L][0]; if (L > 0) { for (int m = 1; m <= L; m++) { aUM[L][m][k] = akm[L][m]; bUM[L][m][k] = bkm[L][m]; } } } } else { for (L = 0; L <= LmaxUM; L++) { aUM[L][0][k] = ak[L][0]; if (L > 0) { for (int m = 1; m <= L; m++) { aUM[L][m][k] = ak[L][m]; bUM[L][m][k] = bk[L][m]; } } } } } free_matrix(ak, 0, LmaxUM, 0, LmaxUM); free_matrix(bk, 0, LmaxUM, 0, LmaxUM); free_matrix(akm, 0, LmaxUM, 0, LmaxUM); free_matrix(bkm, 0, LmaxUM, 0, LmaxUM); co_radUM[0] = co_radUM[1]; co_radUM[-1] = co_radUM[0]; co_radUM[n_radUM+1] = co_radUM[n_radUM]; co_radUM[n_radUM+2] = co_radUM[n_radUM+1]; int n_radp1 = n_radUM + 1; // set unneccessary coefficients to zero for (L = 0; L <= LmaxUM; L++) { for (int n = 0; n <= n_radUM; n++) { bUM[L][0][n] = 0.0; b0UM[L][0][n] = 0.0; } for (int m = 0; m <= L; m++) { aUM[L][m][0] = 0.0; bUM[L][m][0] = 0.0; aUM[L][m][n_radp1] = 0.0; bUM[L][m][n_radp1] = 0.0; a0UM[L][m][0] = 0.0; b0UM[L][m][0] = 0.0; a0UM[L][m][n_radp1] = 0.0; b0UM[L][m][n_radp1] = 0.0; } } // Now lower mantle ak = matrix(0, LmaxLM, 0, LmaxLM); bk = matrix(0, LmaxLM, 0, LmaxLM); akm = matrix(0, LmaxLM, 0, LmaxLM); bkm = matrix(0, LmaxLM, 0, LmaxLM); for (k = 1; k <= n_radLM; k++) { for (L = 0; L <= LmaxLM; L++) { fscanf(fin, "%le", &ak[L][0]); for (int m = 1; m <= L; m++) fscanf(fin, "%le %le", &ak[L][m], &bk[L][m]); a0LM[L][0][k] = ak[L][0]; if (L > 0) { for (m = 1; m <= L; m++) { a0LM[L][m][k] = ak[L][m]; b0LM[L][m][k] = bk[L][m]; } } } if (alpha != 0.0 || beta != 0.0 || gamma != 0.0) { rotate_velocity(alpha, beta, gamma, LmaxLM, ak, bk, akm, bkm); for (L = 0; L <= LmaxLM; L++) { aLM[L][0][k] = akm[L][0]; if (L > 0) { for (int m = 1; m <= L; m++) { aLM[L][m][k] = akm[L][m]; bLM[L][m][k] = bkm[L][m]; } } } } else { for (L = 0; L <= LmaxLM; L++) { aLM[L][0][k] = ak[L][0]; if (L > 0) { for (int m = 1; m <= L; m++) { aLM[L][m][k] = ak[L][m]; bLM[L][m][k] = bk[L][m]; } } } } } free_matrix(ak, 0, LmaxLM, 0, LmaxLM); free_matrix(bk, 0, LmaxLM, 0, LmaxLM); free_matrix(akm, 0, LmaxLM, 0, LmaxLM); free_matrix(bkm, 0, LmaxLM, 0, LmaxLM); fclose(fin); SetSplitNode(co_radUM, co_radLM); n_radp1 = n_radLM + 1; // set unneccessary coefficients to zero for (L = 0; L <= LmaxLM; L++) { for (int n = 0; n <= n_radLM; n++) { bLM[L][0][n] = 0.0; b0LM[L][0][n] = 0.0; } for (int m = 0; m <= L; m++) { aLM[L][m][0] = 0.0; bLM[L][m][0] = 0.0; aLM[L][m][n_radp1] = 0.0; bLM[L][m][n_radp1] = 0.0; a0LM[L][m][0] = 0.0; b0LM[L][m][0] = 0.0; a0LM[L][m][n_radp1] = 0.0; b0LM[L][m][n_radp1] = 0.0; } } return; } void BSP_WHOLE_T::rotate(double alp, double bet, double gam) { alpha = alp; beta = bet; gamma = gam; double **ak = matrix(0, Lmax, 0, Lmax); double **bk = matrix(0, Lmax, 0, Lmax); double **akm = matrix(0, Lmax, 0, Lmax); double **bkm = matrix(0, Lmax, 0, Lmax); int k, L; for (k = 1; k <= n_rad; k++) { for (L = 0; L <= Lmax; L++) { ak[L][0] = a0[L][0][k]; if (L > 0) { for (int m = 1; m <= L; m++) { ak[L][m] = a0[L][m][k]; bk[L][m] = b0[L][m][k]; } } } rotate_velocity(alpha, beta, gamma, Lmax, ak, bk, akm, bkm); for (L = 0; L <= Lmax; L++) { a[L][0][k] = akm[L][0]; if (L > 0) { for (int m = 1; m <= L; m++) { a[L][m][k] = akm[L][m]; b[L][m][k] = bkm[L][m]; } } } } free_matrix(ak, 0, Lmax, 0, Lmax); free_matrix(bk, 0, Lmax, 0, Lmax); free_matrix(akm, 0, Lmax, 0, Lmax); free_matrix(bkm, 0, Lmax, 0, Lmax); } void BSP_SPLIT_T::rotate(double alp, double bet, double gam) { alpha = alp; beta = bet; gamma = gam; double **ak = matrix(0, LmaxUM, 0, LmaxUM); double **bk = matrix(0, LmaxUM, 0, LmaxUM); double **akm = matrix(0, LmaxUM, 0, LmaxUM); double **bkm = matrix(0, LmaxUM, 0, LmaxUM); int k, L; for (k = 1; k <= n_radUM; k++) { for (L = 0; L <= LmaxUM; L++) { ak[L][0] = a0UM[L][0][k]; if (L > 0) { for (int m = 1; m <= L; m++) { ak[L][m] = a0UM[L][m][k]; bk[L][m] = b0UM[L][m][k]; } } } rotate_velocity(alpha, beta, gamma, LmaxUM, ak, bk, akm, bkm); for (L = 0; L <= LmaxUM; L++) { aUM[L][0][k] = akm[L][0]; if (L > 0) { for (int m = 1; m <= L; m++) { aUM[L][m][k] = akm[L][m]; bUM[L][m][k] = bkm[L][m]; } } } } free_matrix(ak, 0, LmaxUM, 0, LmaxUM); free_matrix(bk, 0, LmaxUM, 0, LmaxUM); free_matrix(akm, 0, LmaxUM, 0, LmaxUM); free_matrix(bkm, 0, LmaxUM, 0, LmaxUM); ak = matrix(0, LmaxLM, 0, LmaxLM); bk = matrix(0, LmaxLM, 0, LmaxLM); akm = matrix(0, LmaxLM, 0, LmaxLM); bkm = matrix(0, LmaxLM, 0, LmaxLM); for (k = 1; k <= n_radLM; k++) { for (L = 0; L <= LmaxLM; L++) { ak[L][0] = a0LM[L][0][k]; if (L > 0) { for (int m = 1; m <= L; m++) { ak[L][m] = a0LM[L][m][k]; bk[L][m] = b0LM[L][m][k]; } } } rotate_velocity(alpha, beta, gamma, LmaxLM, ak, bk, akm, bkm); for (L = 0; L <= LmaxLM; L++) { aLM[L][0][k] = akm[L][0]; if (L > 0) { for (int m = 1; m <= L; m++) { aLM[L][m][k] = akm[L][m]; bLM[L][m][k] = bkm[L][m]; } } } } free_matrix(ak, 0, LmaxLM, 0, LmaxLM); free_matrix(bk, 0, LmaxLM, 0, LmaxLM); free_matrix(akm, 0, LmaxLM, 0, LmaxLM); free_matrix(bkm, 0, LmaxLM, 0, LmaxLM); } void CHEB_SPLIT_T::rotate(double alp, double bet, double gam) { alpha = alp; beta = bet; gamma = gam; double **ak = matrix(0, LmaxUM, 0, LmaxUM); double **bk = matrix(0, LmaxUM, 0, LmaxUM); double **akm = matrix(0, LmaxUM, 0, LmaxUM); double **bkm = matrix(0, LmaxUM, 0, LmaxUM); int k, L; for (k = 1; k <= kmaxUM; k++) { for (L = 0; L <= LmaxUM; L++) { ak[L][0] = a0UM[L][0][k]; if (L > 0) { for (int m = 1; m <= L; m++) { ak[L][m] = a0UM[L][m][k]; bk[L][m] = b0UM[L][m][k]; } } } rotate_velocity(alpha, beta, gamma, LmaxUM, ak, bk, akm, bkm); for (L = 0; L <= LmaxUM; L++) { aUM[L][0][k] = akm[L][0]; if (L > 0) { for (int m = 1; m <= L; m++) { aUM[L][m][k] = akm[L][m]; bUM[L][m][k] = bkm[L][m]; } } } } free_matrix(ak, 0, LmaxUM, 0, LmaxUM); free_matrix(bk, 0, LmaxUM, 0, LmaxUM); free_matrix(akm, 0, LmaxUM, 0, LmaxUM); free_matrix(bkm, 0, LmaxUM, 0, LmaxUM); ak = matrix(0, LmaxLM, 0, LmaxLM); bk = matrix(0, LmaxLM, 0, LmaxLM); akm = matrix(0, LmaxLM, 0, LmaxLM); bkm = matrix(0, LmaxLM, 0, LmaxLM); for (k = 1; k <= kmaxLM; k++) { for (L = 0; L <= LmaxLM; L++) { ak[L][0] = a0LM[L][0][k]; if (L > 0) { for (int m = 1; m <= L; m++) { ak[L][m] = a0LM[L][m][k]; bk[L][m] = b0LM[L][m][k]; } } } rotate_velocity(alpha, beta, gamma, LmaxLM, ak, bk, akm, bkm); for (L = 0; L <= LmaxLM; L++) { aLM[L][0][k] = akm[L][0]; if (L > 0) { for (int m = 1; m <= L; m++) { aLM[L][m][k] = akm[L][m]; bLM[L][m][k] = bkm[L][m]; } } } } free_matrix(ak, 0, LmaxLM, 0, LmaxLM); free_matrix(bk, 0, LmaxLM, 0, LmaxLM); free_matrix(akm, 0, LmaxLM, 0, LmaxLM); free_matrix(bkm, 0, LmaxLM, 0, LmaxLM); } void CHEB_WHOLE_T::rotate(double alp, double bet, double gam) { alpha = alp; beta = bet; gamma = gam; double **ak = matrix(0, Lmax, 0, Lmax); double **bk = matrix(0, Lmax, 0, Lmax); double **akm = matrix(0, Lmax, 0, Lmax); double **bkm = matrix(0, Lmax, 0, Lmax); int k, L; for (k = 0; k <= kmax; k++) { for (L = 0; L <= Lmax; L++) { ak[L][0] = a0[L][0][k]; if (L > 0) { for (int m = 1; m <= L; m++) { ak[L][m] = a0[L][m][k]; bk[L][m] = b0[L][m][k]; } } } rotate_velocity(alpha, beta, gamma, Lmax, ak, bk, akm, bkm); for (L = 0; L <= Lmax; L++) { a[L][0][k] = akm[L][0]; if (L > 0) { for (int m = 1; m <= L; m++) { a[L][m][k] = akm[L][m]; b[L][m][k] = bkm[L][m]; } } } } free_matrix(ak, 0, Lmax, 0, Lmax); free_matrix(bk, 0, Lmax, 0, Lmax); free_matrix(akm, 0, Lmax, 0, Lmax); free_matrix(bkm, 0, Lmax, 0, Lmax); } /* function dv2v for spherical spline model. * Input: radius (km) * xlat,xlong: co-latitude, co-longitude in radians * output: dv =velocity anomaly divided by average velocity (PREM) * dvdx =d(dv)/dx (x=r/6371, normalized radius) * dvdtheta =d(dv)/dtheta * dvdphi =d(dv)/dphi * Warning: hacks ahead. This function is not quite readable. * Efficiency is the number one concern here. This function * consumes most (more than 70%) of the CPU time in ray tracing. */ void SPH_SPLINE_T::dv2v(double xr, double xlat, double xlong, double *dv, double *dvdx, double *dvdtheta, double *dvdphi) const { double rr = xr * RADIUS; double B[4], DB[4], *ar, *ar_d; double B0, B1, B2, B3, DB0, DB1, DB2, DB3; int kp1, kp2, km1, i; int k; *dv = 0.0; *dvdx = 0.0; *dvdtheta = 0.0; *dvdphi = 0.0; if (rr < RCMB || rr > RMOHO) return; k = radialBSpline->DBBSplKer(B, DB, RADIUS - rr); if (k == 0 || k > n_rad) return; ar = vector(1, n_node); ar_d = vector(1, n_node); kp1 = k + 1; kp2 = k + 2; km1 = k - 1; B0 = B[0]; B1 = B[1]; B2 = B[2]; B3 = B[3]; DB0 = DB[0]; DB1 = DB[1]; DB2 = DB[2]; DB3 = DB[3]; double *coeff_k = coeff[k]; double *coeff_kp1 = coeff[kp1]; double *coeff_km1, *coeff_kp2; coeff_km1 = coeff[km1]; coeff_kp2 = coeff[kp2]; double coeff_k_i, coeff_kp1_i, coeff_km1_i, coeff_kp2_i; /* The following codes may looks stupid in term of number of lines. * However, it is quite efficient. */ if (k == 1) { for (i = 1; i <= n_node; i++) { /* for all horizontal nodes */ coeff_k_i = coeff_k[i]; coeff_kp1_i = coeff_kp1[i]; coeff_kp2_i = coeff_kp2[i]; ar[i] = coeff_k_i * B1 + coeff_kp1_i * B2 + coeff_kp2_i * B3; ar_d[i] = coeff_k_i * DB1 + coeff_kp1_i * DB2 + coeff_kp2_i * DB3; } } else if (kp1 != n_rad) { for (i = 1; i <= n_node; i++) { /* for all horizontal nodes */ coeff_k_i = coeff_k[i]; coeff_kp1_i = coeff_kp1[i]; coeff_km1_i = coeff_km1[i]; coeff_kp2_i = coeff_kp2[i]; ar[i] = coeff_km1_i * B0 + coeff_k_i * B1 + coeff_kp1_i * B2 + coeff_kp2_i * B3; ar_d[i] = coeff_km1_i * DB0 + coeff_k_i * DB1 + coeff_kp1_i * DB2 + coeff_kp2_i * DB3; } } else { // k == n_rad - 1 for (i = 1; i <= n_node; i++) { /* for all horizontal nodes */ coeff_k_i = coeff_k[i]; coeff_kp1_i = coeff_kp1[i]; coeff_km1_i = coeff_km1[i]; ar[i] = coeff_km1_i * B0 + coeff_k_i * B1 + coeff_kp1_i * B2; ar_d[i] = coeff_km1_i * DB0 + coeff_k_i * DB1 + coeff_kp1_i * DB2; } } *dv = d_SumSphBspline(ar, ar_d, xlat, xlong, dvdx, dvdtheta, dvdphi); free_vector(ar, 1, n_node); free_vector(ar_d, 1, n_node); return; } /* Function alldv2v-- * xr: normalized radius (0-1) * xlat,xlong: co-latitude, co-longitude in radians * output: dv =velocity anomaly divided by average velocity (PREM) * dvdx =d(dv)/dx (x=r/6371, normalized radius) * dvdtheta =d(dv)/dtheta * dvdphi =d(dv)/dphi * d2vd2x =d^2(dv)/dx^2 * d2vd2theta =d^2(dv)/dtheta^2 * d2vdxtheta =d^2(dv)/dxdtheta * d2vdxdphi =d^2(dv)/dxdphi * d2vdthetadphi=d^2(dv)/dthetadphi * Warning: hacks ahead. This function is not quite readable. * Efficiency is the number one concern here. */ void CHEB_WHOLE_T::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 rr = xr * RADIUS; double *cosm, *sinm, **p, **dp, **ddp; double *func, *dfunc, *ddfunc; int Lmaxp1; int m, L, k; double phim, plm, dplm, ddplm, cosphi, sinphi, sinphim, cosphim; double funck, dfunck, ddfunck, ac, pac, dpac, Aklm, Bklm; double acbs, pacbs, dpacbs, asbc, pasbc; *dv = 0.0; *dvdx = 0.0; *dvdtheta = 0.0; *dvdphi = 0.0; *d2vd2x = 0.0; *d2vdxdtheta = 0.0; *d2vdxdphi = 0.0; *d2vdthetadphi = 0.0; *d2vd2theta = 0.0; if (rr < RCMB || rr > RMOHO) return; double np_dv = 0.0; // for efficiency double np_dvdx = 0.0; double np_dvdtheta = 0.0; double np_dvdphi = 0.0; double np_d2vd2x = 0.0; double np_d2vdxdtheta = 0.0; double np_d2vdxdphi = 0.0; double np_d2vdthetadphi = 0.0; double np_d2vd2theta = 0.0; Lmaxp1 = Lmax + 1; func = vector(0, kmax); dfunc = vector(0, kmax); ddfunc = vector(0, kmax); cosm = vector(0, Lmax); sinm = vector(0, Lmax); ddradfun(radfun_t, func, dfunc, ddfunc, rr, kmax); // We don't need to calculate cos and sin at m = 0 for (m = 1; m <= Lmax; m++) { phim = m * xlong; cosm[m] = cos(phim); sinm[m] = sin(phim); } p = matrix(0, Lmaxp1, 0, Lmaxp1); dp = matrix(0, Lmaxp1, 0, Lmaxp1); ddp = matrix(0, Lmax, 0, Lmax); ddYlm(xlat, Lmax, p, dp, ddp); double *alm, *blm, **al, **bl, *pl, *dpl, *ddpl; // for efficiency for (L = 0; L <= Lmax; L++) { al = a[L]; // for efficiency bl = b[L]; pl = p[L]; dpl = dp[L]; ddpl = ddp[L]; // m = 0 first plm = pl[0]; dplm = dpl[0]; ddplm = ddpl[0]; alm = al[0]; // for efficiency for (k = 0; k <= kmax; k++) { funck = func[k]; dfunck = dfunc[k]; ddfunck = ddfunc[k]; ac = alm[k]; pac = plm * ac; dpac = dplm * ac; np_dv += funck * pac; np_dvdx += dfunck * pac; np_d2vd2x += ddfunck * pac; np_dvdtheta += funck * dpac; np_d2vd2theta += funck * ddplm * ac; np_d2vdxdtheta += dfunck * dpac; } // m > 0 for (m = 1; m <= L; m++) { plm = pl[m]; dplm = dpl[m]; ddplm = ddpl[m]; cosphi = cosm[m]; sinphi = sinm[m]; sinphim = sinphi * m; cosphim = cosphi * m; alm = al[m]; // for efficiency blm = bl[m]; for (k = 0; k <= kmax; k++) { funck = func[k]; dfunck = dfunc[k]; ddfunck = ddfunc[k]; /* add all spherical harmonics together */ /* Sorry, it is not easy to read this part, but !efficiency is more important here */ Aklm = alm[k]; Bklm = blm[k]; acbs = Aklm * cosphi + Bklm * sinphi; pacbs = plm * acbs; dpacbs = dplm * acbs; asbc = -Aklm * sinphim + Bklm * cosphim; pasbc = asbc * plm; np_dv += funck * pacbs; np_dvdx += dfunck * pacbs; np_d2vd2x += ddfunck * pacbs; np_dvdtheta += funck * dpacbs; np_d2vdxdtheta += dfunck * dpacbs; np_dvdphi += funck * pasbc; np_d2vdxdphi += dfunck * pasbc; np_d2vdthetadphi += funck * dplm * asbc; np_d2vd2theta += funck * ddplm * acbs; } } } free_matrix(p, 0, Lmaxp1, 0, Lmaxp1); free_matrix(dp, 0, Lmaxp1, 0, Lmaxp1); free_matrix(ddp, 0, Lmax, 0, Lmax); free_vector(func, 0, kmax); free_vector(dfunc, 0, kmax); free_vector(ddfunc, 0, kmax); free_vector(cosm, 0, Lmax); free_vector(sinm, 0, Lmax); *dv = np_dv; *dvdx = np_dvdx; *dvdtheta = np_dvdtheta; *dvdphi = np_dvdphi; *d2vd2x = np_d2vd2x; *d2vdxdtheta = np_d2vdxdtheta; *d2vdxdphi = np_d2vdxdphi; *d2vdthetadphi = np_d2vdthetadphi; *d2vd2theta = np_d2vd2theta; return; } /* Function alldv2v-- * xr: normalized radius (0-1) * xlat,xlong: co-latitude, co-longitude in radians * output: dv =velocity anomaly divided by average velocity (PREM) * dvdx =d(dv)/dx (x=r/6371, normalized radius) * dvdtheta =d(dv)/dtheta * dvdphi =d(dv)/dphi * d2vd2x =d^2(dv)/dx^2 * d2vd2theta =d^2(dv)/dtheta^2 * d2vdxtheta =d^2(dv)/dxdtheta * d2vdxdphi =d^2(dv)/dxdphi * d2vdthetadphi=d^2(dv)/dthetadphi * Warning: hacks ahead. This function is not quite readable. * Efficiency is the number one concern here. */ void CHEB_SPLIT_T::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 rr = xr * RADIUS; double *cosm, *sinm, **p, **dp, **ddp; double *func, *dfunc, *ddfunc; int Lmaxp1; int m, L, k, Lmax, kmax; double phim, plm, dplm, ddplm, cosphi, sinphi, sinphim, cosphim; double funck, dfunck, ddfunck, ac, pac, dpac, Aklm, Bklm; double acbs, pacbs, dpacbs, asbc, pasbc; double ***a, ***b; *dv = 0.0; *dvdx = 0.0; *dvdtheta = 0.0; *dvdphi = 0.0; *d2vd2x = 0.0; *d2vdxdtheta = 0.0; *d2vdxdphi = 0.0; *d2vdthetadphi = 0.0; *d2vd2theta = 0.0; if (rr < RCMB || rr > RMOHO) return; if (rr >= R670) { // upper mantle kmax = kmaxUM; a = aUM; b = bUM; Lmax = LmaxUM; } else { // Lower mantle kmax = kmaxLM; a = aLM; b = bLM; Lmax = LmaxLM; } double np_dv = 0.0; // for efficiency double np_dvdx = 0.0; double np_dvdtheta = 0.0; double np_dvdphi = 0.0; double np_d2vd2x = 0.0; double np_d2vdxdtheta = 0.0; double np_d2vdxdphi = 0.0; double np_d2vdthetadphi = 0.0; double np_d2vd2theta = 0.0; Lmaxp1 = Lmax + 1; func = vector(0, kmax); dfunc = vector(0, kmax); ddfunc = vector(0, kmax); cosm = vector(0, Lmax); sinm = vector(0, Lmax); ddradfun(radfun_t, func, dfunc, ddfunc, rr, kmax); // We don't need to calculate cos and sin at m = 0 for (m = 1; m <= Lmax; m++) { phim = m * xlong; cosm[m] = cos(phim); sinm[m] = sin(phim); } p = matrix(0, Lmaxp1, 0, Lmaxp1); dp = matrix(0, Lmaxp1, 0, Lmaxp1); ddp = matrix(0, Lmax, 0, Lmax); ddYlm(xlat, Lmax, p, dp, ddp); double *alm, *blm, **al, **bl, *pl, *dpl, *ddpl; // for efficiency for (L = 0; L <= Lmax; L++) { al = a[L]; // for efficiency bl = b[L]; pl = p[L]; dpl = dp[L]; ddpl = ddp[L]; // m = 0 first plm = pl[0]; dplm = dpl[0]; ddplm = ddpl[0]; alm = al[0]; // for efficiency for (k = 0; k <= kmax; k++) { funck = func[k]; dfunck = dfunc[k]; ddfunck = ddfunc[k]; ac = alm[k]; pac = plm * ac; dpac = dplm * ac; np_dv += funck * pac; np_dvdx += dfunck * pac; np_d2vd2x += ddfunck * pac; np_dvdtheta += funck * dpac; np_d2vd2theta += funck * ddplm * ac; np_d2vdxdtheta += dfunck * dpac; } // m > 0 for (m = 1; m <= L; m++) { plm = pl[m]; dplm = dpl[m]; ddplm = ddpl[m]; cosphi = cosm[m]; sinphi = sinm[m]; sinphim = sinphi * m; cosphim = cosphi * m; alm = al[m]; // for efficiency blm = bl[m]; for (k = 0; k <= kmax; k++) { funck = func[k]; dfunck = dfunc[k]; ddfunck = ddfunc[k]; /* add all spherical harmonics together */ /* Sorry, it is not easy to read this part, but !efficiency is more important here */ Aklm = alm[k]; Bklm = blm[k]; acbs = Aklm * cosphi + Bklm * sinphi; pacbs = plm * acbs; dpacbs = dplm * acbs; asbc = -Aklm * sinphim + Bklm * cosphim; pasbc = asbc * plm; np_dv += funck * pacbs; np_dvdx += dfunck * pacbs; np_d2vd2x += ddfunck * pacbs; np_dvdtheta += funck * dpacbs; np_d2vdxdtheta += dfunck * dpacbs; np_dvdphi += funck * pasbc; np_d2vdxdphi += dfunck * pasbc; np_d2vdthetadphi += funck * dplm * asbc; np_d2vd2theta += funck * ddplm * acbs; } } } free_matrix(p, 0, Lmaxp1, 0, Lmaxp1); free_matrix(dp, 0, Lmaxp1, 0, Lmaxp1); free_matrix(ddp, 0, Lmax, 0, Lmax); free_vector(func, 0, kmax); free_vector(dfunc, 0, kmax); free_vector(ddfunc, 0, kmax); free_vector(cosm, 0, Lmax); free_vector(sinm, 0, Lmax); *dv = np_dv; *dvdx = np_dvdx; *dvdtheta = np_dvdtheta; *dvdphi = np_dvdphi; *d2vd2x = np_d2vd2x; *d2vdxdtheta = np_d2vdxdtheta; *d2vdxdphi = np_d2vdxdphi; *d2vdthetadphi = np_d2vdthetadphi; *d2vd2theta = np_d2vd2theta; return; } /* Function alldv2v-- * xr: normalized radius (0-1) * xlat,xlong: co-latitude, co-longitude in radians * output: dv =velocity anomaly divided by average velocity (PREM) * dvdx =d(dv)/dx (x=r/6371, normalized radius) * dvdtheta =d(dv)/dtheta * dvdphi =d(dv)/dphi * d2vd2x =d^2(dv)/dx^2 * d2vd2theta =d^2(dv)/dtheta^2 * d2vdxtheta =d^2(dv)/dxdtheta * d2vdxdphi =d^2(dv)/dxdphi * d2vdthetadphi=d^2(dv)/dthetadphi * Warning: hacks ahead. This function is not quite readable. * Efficiency is the number one concern here. */ void BSP_WHOLE_T::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 rr = xr * RADIUS; double *cosm, *sinm, **p, **dp, **ddp; int Lmaxp1; int m, L; double phim, plm, dplm, ddplm, cosphi, sinphi, sinphim, cosphim; *dv = 0.0; *dvdx = 0.0; *dvdtheta = 0.0; *dvdphi = 0.0; *d2vd2x = 0.0; *d2vdxdtheta = 0.0; *d2vdxdphi = 0.0; *d2vdthetadphi = 0.0; *d2vd2theta = 0.0; if (rr < RCMB || rr > RMOHO) return; double np_dv = 0.0; // for efficiency double np_dvdx = 0.0; double np_dvdtheta = 0.0; double np_dvdphi = 0.0; double np_d2vd2x = 0.0; double np_d2vdxdtheta = 0.0; double np_d2vdxdphi = 0.0; double np_d2vdthetadphi = 0.0; double np_d2vd2theta = 0.0; Lmaxp1 = Lmax + 1; double B[4], DB[4], DDB[4]; int k = radialBSpline->DDBBSplKer(B, DB, DDB, (RADIUS - rr)); if (k == 0 || k > n_rad) return; int km1 = k - 1; int kp1 = k + 1; int kp2 = k + 2; double B0 = B[0]; double B1 = B[1]; double B2 = B[2]; double B3 = B[3]; double DB0 = DB[0]; double DB1 = DB[1]; double DB2 = DB[2]; double DB3 = DB[3]; double DDB0 = DDB[0]; double DDB1 = DDB[1]; double DDB2 = DDB[2]; double DDB3 = DDB[3]; cosm = vector(0, Lmax); sinm = vector(0, Lmax); for (m = 0; m <= Lmax; m++) { phim = m * xlong; cosm[m] = cos(phim); sinm[m] = sin(phim); } p = matrix(0, Lmaxp1, 0, Lmaxp1); dp = matrix(0, Lmaxp1, 0, Lmaxp1); ddp = matrix(0, Lmax, 0, Lmax); ddYlm(xlat, Lmax, p, dp, ddp); double *alm, *blm, **al, **bl, *pl, *dpl, *ddpl; // for efficiency double Almk, Blmk, Almkp1, Blmkp1, acbs, acbsp1, asbc, asbcp1; double Almkm1, Blmkm1, acbsm1, asbcm1; double Almkp2, Blmkp2, acbsp2, asbcp2; for (L = 0; L <= Lmax; L++) { pl = p[L]; dpl = dp[L]; ddpl = ddp[L]; al = a[L]; bl = b[L]; for (m = 0; m <= L; m++) { alm = al[m]; blm = bl[m]; plm = pl[m]; dplm = dpl[m]; ddplm = ddpl[m]; cosphi = cosm[m]; sinphi = sinm[m]; sinphim = sinphi * m; cosphim = cosphi * m; Almk = alm[k]; Blmk = blm[k]; Almkp1 = alm[kp1]; Blmkp1 = blm[kp1]; acbs = Almk * cosphi + Blmk * sinphi; acbsp1 = (Almkp1 * cosphi + Blmkp1 * sinphi); asbc = (-Almk * sinphim + Blmk * cosphim); asbcp1 = (-Almkp1 * sinphim + Blmkp1 * cosphim); np_dv += B1 * plm * acbs + B2 * plm * acbsp1; np_dvdx += DB1 * plm * acbs + DB2 * plm * acbsp1; np_d2vdxdtheta += DB1 * dplm * acbs + DB2 * dplm * acbsp1; np_d2vd2x += DDB1 * plm * acbs + DDB2 * plm * acbsp1; np_dvdtheta += B1 * dplm * acbs + B2 * dplm * acbsp1; np_d2vd2theta += B1 * ddplm * acbs + B2 * ddplm * acbsp1; np_dvdphi += B1 * plm * asbc + B2 * plm * asbcp1; np_d2vdxdphi += DB1 * plm * asbc + DB2 * plm * asbcp1; np_d2vdthetadphi += B1 * dplm * asbc + B2 * dplm * asbcp1; if (k != 1) { Almkm1 = a[L][m][km1]; Blmkm1 = b[L][m][km1]; acbsm1 = (Almkm1 * cosphi + Blmkm1 * sinphi); asbcm1 = (-Almkm1 * sinphim + Blmkm1 * cosphim); np_dv += B0 * plm * acbsm1; np_dvdx += DB0 * plm * acbsm1; np_d2vd2x += DDB0 * plm * acbsm1; np_dvdtheta += B0 * dplm * acbsm1; np_d2vdxdtheta += DB0 * dplm * acbsm1; np_d2vd2theta += B0 * ddplm * acbsm1; np_dvdphi += B0 * plm * asbcm1; np_d2vdxdphi += DB0 * plm * asbcm1; np_d2vdthetadphi += B0 * dplm * asbcm1; } if (kp1 != n_rad) { Almkp2 = a[L][m][kp2]; Blmkp2 = b[L][m][kp2]; acbsp2 = (Almkp2 * cosphi + Blmkp2 * sinphi); asbcp2 = (-Almkp2 * sinphim + Blmkp2 * cosphim); np_dv += B3 * plm * acbsp2; np_dvdx += DB3 * plm * acbsp2; np_d2vd2x += DDB3 * plm * acbsp2; np_dvdtheta += B3 * dplm * acbsp2; np_d2vdxdtheta += DB3 * dplm * acbsp2; np_d2vd2theta += B3 * ddplm * acbsp2; np_dvdphi += B3 * plm * asbcp2; np_d2vdxdphi += DB3 * plm * asbcp2; np_d2vdthetadphi += B3 * dplm * asbcp2; } } } free_matrix(p, 0, Lmaxp1, 0, Lmaxp1); free_matrix(dp, 0, Lmaxp1, 0, Lmaxp1); free_matrix(ddp, 0, Lmax, 0, Lmax); free_vector(cosm, 0, Lmax); free_vector(sinm, 0, Lmax); *dv = np_dv; *dvdx = -RADIUS * np_dvdx; *dvdtheta = np_dvdtheta; *dvdphi = np_dvdphi; *d2vd2x = RADIUS * RADIUS * np_d2vd2x; *d2vdxdtheta = -RADIUS * np_d2vdxdtheta; *d2vdxdphi = -RADIUS * np_d2vdxdphi; *d2vdthetadphi = np_d2vdthetadphi; *d2vd2theta = np_d2vd2theta; return; } /* Function alldv2v-- * xr: normalized radius (0-1) * xlat,xlong: co-latitude, co-longitude in radians * output: dv =velocity anomaly divided by average velocity (PREM) * dvdx =d(dv)/dx (x=r/6371, normalized radius) * dvdtheta =d(dv)/dtheta * dvdphi =d(dv)/dphi * d2vd2x =d^2(dv)/dx^2 * d2vd2theta =d^2(dv)/dtheta^2 * d2vdxtheta =d^2(dv)/dxdtheta * d2vdxdphi =d^2(dv)/dxdphi * d2vdthetadphi=d^2(dv)/dthetadphi * Warning: hacks ahead. This function is not quite readable. * Efficiency is the number one concern here. */ void BSP_SPLIT_T::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 rr = xr * RADIUS; double *cosm, *sinm, **p, **dp, **ddp; double *co_rad; RadialBSpline *radialBSpline; double ***a, ***b; int Lmax, n_rad; int Lmaxp1; int m, L; double phim, plm, dplm, ddplm, cosphi, sinphi, sinphim, cosphim; *dv = 0.0; *dvdx = 0.0; *dvdtheta = 0.0; *dvdphi = 0.0; *d2vd2x = 0.0; *d2vdxdtheta = 0.0; *d2vdxdphi = 0.0; *d2vdthetadphi = 0.0; *d2vd2theta = 0.0; if (rr < RCMB || rr > RMOHO) return; if (rr >= R670) { // upper mantle co_rad = co_radUM; n_rad = n_radUM; a = aUM; b = bUM; Lmax = LmaxUM; radialBSpline = radialBSplineUM; } else { // Lower mantle co_rad = co_radLM; n_rad = n_radLM; a = aLM; b = bLM; Lmax = LmaxLM; radialBSpline = radialBSplineLM; } double np_dv = 0.0; // for efficiency double np_dvdx = 0.0; double np_dvdtheta = 0.0; double np_dvdphi = 0.0; double np_d2vd2x = 0.0; double np_d2vdxdtheta = 0.0; double np_d2vdxdphi = 0.0; double np_d2vdthetadphi = 0.0; double np_d2vd2theta = 0.0; Lmaxp1 = Lmax + 1; double B[4], DB[4], DDB[4]; int k = radialBSpline->DDBBSplKer(B, DB, DDB, (RADIUS - rr)); if (k == 0 || k > n_rad) return; int km1 = k - 1; int kp1 = k + 1; int kp2 = k + 2; double B0 = B[0]; double B1 = B[1]; double B2 = B[2]; double B3 = B[3]; double DB0 = DB[0]; double DB1 = DB[1]; double DB2 = DB[2]; double DB3 = DB[3]; double DDB0 = DDB[0]; double DDB1 = DDB[1]; double DDB2 = DDB[2]; double DDB3 = DDB[3]; cosm = vector(0, Lmax); sinm = vector(0, Lmax); for (m = 0; m <= Lmax; m++) { phim = m * xlong; cosm[m] = cos(phim); sinm[m] = sin(phim); } p = matrix(0, Lmaxp1, 0, Lmaxp1); dp = matrix(0, Lmaxp1, 0, Lmaxp1); ddp = matrix(0, Lmax, 0, Lmax); ddYlm(xlat, Lmax, p, dp, ddp); double *alm, *blm, **al, **bl, *pl, *dpl, *ddpl; // for efficiency double Almk, Blmk, Almkp1, Blmkp1, acbs, acbsp1, asbc, asbcp1; double Almkm1, Blmkm1, acbsm1, asbcm1; double Almkp2, Blmkp2, acbsp2, asbcp2; for (L = 0; L <= Lmax; L++) { pl = p[L]; dpl = dp[L]; ddpl = ddp[L]; al = a[L]; bl = b[L]; for (m = 0; m <= L; m++) { alm = al[m]; blm = bl[m]; plm = pl[m]; dplm = dpl[m]; ddplm = ddpl[m]; cosphi = cosm[m]; sinphi = sinm[m]; sinphim = sinphi * m; cosphim = cosphi * m; Almk = alm[k]; Blmk = blm[k]; Almkp1 = alm[kp1]; Blmkp1 = blm[kp1]; acbs = Almk * cosphi + Blmk * sinphi; acbsp1 = (Almkp1 * cosphi + Blmkp1 * sinphi); asbc = (-Almk * sinphim + Blmk * cosphim); asbcp1 = (-Almkp1 * sinphim + Blmkp1 * cosphim); np_dv += B1 * plm * acbs + B2 * plm * acbsp1; np_dvdx += DB1 * plm * acbs + DB2 * plm * acbsp1; np_d2vdxdtheta += DB1 * dplm * acbs + DB2 * dplm * acbsp1; np_d2vd2x += DDB1 * plm * acbs + DDB2 * plm * acbsp1; np_dvdtheta += B1 * dplm * acbs + B2 * dplm * acbsp1; np_d2vd2theta += B1 * ddplm * acbs + B2 * ddplm * acbsp1; np_dvdphi += B1 * plm * asbc + B2 * plm * asbcp1; np_d2vdxdphi += DB1 * plm * asbc + DB2 * plm * asbcp1; np_d2vdthetadphi += B1 * dplm * asbc + B2 * dplm * asbcp1; if (k != 1) { Almkm1 = a[L][m][km1]; Blmkm1 = b[L][m][km1]; acbsm1 = (Almkm1 * cosphi + Blmkm1 * sinphi); asbcm1 = (-Almkm1 * sinphim + Blmkm1 * cosphim); np_dv += B0 * plm * acbsm1; np_dvdx += DB0 * plm * acbsm1; np_d2vd2x += DDB0 * plm * acbsm1; np_dvdtheta += B0 * dplm * acbsm1; np_d2vdxdtheta += DB0 * dplm * acbsm1; np_d2vd2theta += B0 * ddplm * acbsm1; np_dvdphi += B0 * plm * asbcm1; np_d2vdxdphi += DB0 * plm * asbcm1; np_d2vdthetadphi += B0 * dplm * asbcm1; } if (kp1 != n_rad) { Almkp2 = a[L][m][kp2]; Blmkp2 = b[L][m][kp2]; acbsp2 = (Almkp2 * cosphi + Blmkp2 * sinphi); asbcp2 = (-Almkp2 * sinphim + Blmkp2 * cosphim); np_dv += B3 * plm * acbsp2; np_dvdx += DB3 * plm * acbsp2; np_d2vd2x += DDB3 * plm * acbsp2; np_dvdtheta += B3 * dplm * acbsp2; np_d2vdxdtheta += DB3 * dplm * acbsp2; np_d2vd2theta += B3 * ddplm * acbsp2; np_dvdphi += B3 * plm * asbcp2; np_d2vdxdphi += DB3 * plm * asbcp2; np_d2vdthetadphi += B3 * dplm * asbcp2; } } } free_matrix(p, 0, Lmaxp1, 0, Lmaxp1); free_matrix(dp, 0, Lmaxp1, 0, Lmaxp1); free_matrix(ddp, 0, Lmax, 0, Lmax); free_vector(cosm, 0, Lmax); free_vector(sinm, 0, Lmax); *dv = np_dv; *dvdx = -RADIUS * np_dvdx; *dvdtheta = np_dvdtheta; *dvdphi = np_dvdphi; *d2vd2x = RADIUS * RADIUS * np_d2vd2x; *d2vdxdtheta = -RADIUS * np_d2vdxdtheta; *d2vdxdphi = -RADIUS * np_d2vdxdphi; *d2vdthetadphi = np_d2vdthetadphi; *d2vd2theta = np_d2vd2theta; return; } /* Function alldv2v-- * xr: normalized radius (0-1) * xlat,xlong: co-latitude, co-longitude in radians * output: dv =velocity anomaly divided by average velocity (PREM) * dvdx =d(dv)/dx (x=r/6371, normalized radius) * dvdtheta =d(dv)/dtheta * dvdphi =d(dv)/dphi * d2vd2x =d^2(dv)/dx^2 * d2vd2theta =d^2(dv)/dtheta^2 * d2vdxtheta =d^2(dv)/dxdtheta * d2vdxdphi =d^2(dv)/dxdphi * d2vdthetadphi=d^2(dv)/dthetadphi * Warning: not quite readable, hacks ahead. */ void SPH_SPLINE_T::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 rr = xr * RADIUS; int kp1, kp2, km1, i; int k; double B[4], DB[4], DDB[4], *ar, *ar_d, *ar_dd; double B0, B1, B2, B3, DB0, DB1, DB2, DB3; double DDB0, DDB1, DDB2, DDB3; *dv = 0.0; *dvdx = 0.0; *dvdtheta = 0.0; *dvdphi = 0.0; *d2vd2x = 0.0; *d2vdxdtheta = 0.0; *d2vdxdphi = 0.0; *d2vdthetadphi = 0.0; *d2vd2theta = 0.0; if (rr < RCMB || rr > RMOHO) return; k = radialBSpline->DDBBSplKer(B, DB, DDB, RADIUS - rr); if (k == 0) return; ar = vector(1, n_node); ar_d = vector(1, n_node); ar_dd = vector(1, n_node); kp1 = k + 1; kp2 = k + 2; km1 = k - 1; B0 = B[0]; B1 = B[1]; B2 = B[2]; B3 = B[3]; DB0 = DB[0]; DB1 = DB[1]; DB2 = DB[2]; DB3 = DB[3]; DDB0 = DDB[0]; DDB1 = DDB[1]; DDB2 = DDB[2]; DDB3 = DDB[3]; double *coeff_k = coeff[k]; double *coeff_kp1 = coeff[kp1]; double *coeff_km1, *coeff_kp2; coeff_km1 = coeff[km1]; coeff_kp2 = coeff[kp2]; double coeff_k_i, coeff_kp1_i, coeff_km1_i, coeff_kp2_i; /* The following codes may looks stupid in term of number of lines. * However, it is efficiency-oriented. */ if (k == 1) { for (i = 1; i <= n_node; i++) { /* for all horizontal nodes */ coeff_k_i = coeff_k[i]; coeff_kp1_i = coeff_kp1[i]; coeff_kp2_i = coeff_kp2[i]; ar[i] = coeff_k_i * B1 + coeff_kp1_i * B2 + coeff_kp2_i * B3; ar_d[i] = coeff_k_i * DB1 + coeff_kp1_i * DB2 + coeff_kp2_i * DB3; ar_dd[i] = coeff_k_i * DDB1 + coeff_kp1_i * DDB2 + coeff_kp2_i * DDB3; } } else if (kp1 != n_rad) { for (i = 1; i <= n_node; i++) { /* for all horizontal nodes */ coeff_k_i = coeff_k[i]; coeff_kp1_i = coeff_kp1[i]; coeff_km1_i = coeff_km1[i]; coeff_kp2_i = coeff_kp2[i]; ar[i] = coeff_km1_i * B0 + coeff_k_i * B1 + coeff_kp1_i * B2 + coeff_kp2_i * B3; ar_d[i] = coeff_km1_i * DB0 + coeff_k_i * DB1 + coeff_kp1_i * DB2 + coeff_kp2_i * DB3; ar_dd[i] = coeff_km1_i * DDB0 + coeff_k_i * DDB1 + coeff_kp1_i * DDB2 + coeff_kp2_i * DDB3; } } else { // k == n_rad - 1 for (i = 1; i <= n_node; i++) { /* for all horizontal nodes */ coeff_k_i = coeff_k[i]; coeff_kp1_i = coeff_kp1[i]; coeff_km1_i = coeff_km1[i]; ar[i] = coeff_km1_i * B0 + coeff_k_i * B1 + coeff_kp1_i * B2; ar_d[i] = coeff_km1_i * DB0 + coeff_k_i * DB1 + coeff_kp1_i * DB2; ar_dd[i] = coeff_km1_i * DDB0 + coeff_k_i * DDB1 + coeff_kp1_i * DDB2; } } *dv = dd_SumSphBspline(ar, ar_d, ar_dd, xlat, xlong, dvdx, dvdtheta, dvdphi, d2vd2x, d2vd2theta, d2vdxdtheta, d2vdxdphi, d2vdthetadphi); free_vector(ar, 1, n_node); free_vector(ar_d, 1, n_node); free_vector(ar_dd, 1, n_node); return; } /* Returns values of radial functions at specified radius. * For some iopt's, it computes all radial orders (e.g. 4 for M84A) * even though the calling code may only use a single radial order * at a time. * * input: * iopt = choice of radial functions * 1 = Legendre polynomials in upper mantle, as in M84A * 2 = " " " lower " , as in LO2.56 * 3 = constants, 3 layer upper mantle * 4 = constants, crustal portion of M84C * 5 = constant, single layer upper mantle * 6 = single layer crust, from Smith and Masters * 7 = Tanimoto model, 10 depth knots, lin. interp. between knots * 8 = 4 knots, linear interpolation between knots * 9 = Tanimoto whole mantle model, 11 layers * 10= Inoue et al. whole mantle model * 11= 4 layer (5 knot) upper mantle (lin. interp. between knots) * 12= Legendre polynomial l. m. with D" layer * 13= 4 knots, linear interpolation between knots (PREM crust) * 14= Chebyshev upper mantle k=0 to 4 * 15= Chebyshev lower mantle k=0 to 8 * 16= Chebyshev whole mantle k=0 to 8 (increased to 13) * 17= 14 layer whole mantle * 18= 14 layer upper mantle (Zhang's model - reduced from 26 layers) * 19= experimental * r = radius (km) * output: * func[0:kmax] = radial function values * dfunc[0:kmax] = radial function derivatives * ddfunc[0:kmax] = second derivatives */ void ddradfun(RADFUN_T iopt, double *func, double *dfunc, double *ddfunc, double r, int kmax) { double u, u2, u2r, u2x, u2xsq; double *xnode, B[4], DB[4], DDB[4]; const int maxknot = 20; int i, k0, k; int kmin0, kmax0; static double conc[maxknot] = { 0.70710678118655, 1.2247448713916, 1.0350983390135, 1.0145993123918, 1.0080322575484, 1.0050890913907, 1.0035149493262, 1.0025740068320, 1.0019665702378, 1.0015515913133, 1.0012554932754, 1.0010368069141, 1.0008707010792, 1.0007415648034, 1.0, 1.0, // not exactly right 1.0, 1.0, 1.0, 1.0 }; for (k = 0; k <= kmax; k++) ddfunc[k] = dfunc[k] = func[k] = 0.0; if (kmax < 0) { fprintf(stderr, "radial function order underflow! kmax=%d\n", kmax); exitLog("radial function order underflow!"); } switch (iopt) { case CHEB_UPPER: /* Chebyshev polynomial in upper mantle */ if (r > RMOHO || r < R670) return; u = (r + r - RMOHO - R670) / (RMOHO - R670); u2r = 2.0 / (RMOHO - R670); break; case CHEB_LOWWER: /* Chebyshev polynomial in lower mantle */ if (r > R670 || r < RCMB) return; u = (r + r - R670 - RCMB) / (R670 - RCMB); u2r = 2.0 / (R670 - RCMB); break; case CHEB_WHOLE: /* Chebyshev polynomial for whole mantle */ if (r > RMOHO || r < RCMB) return; u = (r + r - RMOHO - RCMB) / (RMOHO - RCMB); u2r = 2.0 / (RMOHO - RCMB); break; case B_SPLINE: if (r > RMOHO || r < RCMB) return; xnode = vector(-1, kmax+2); SetNode(xnode); k0 = DDBBSplKer(B, DB, DDB, RADIUS - r, xnode, kmax); free_vector(xnode, -1, kmax+2); if (k0 <= 0 || k0 >= kmax) return; if (k0 >= 2) kmin0 = k0; else kmin0 = k0 + 1; if (k0 != (kmax-1)) kmax0 = k0 + 3; else kmax0 = k0 + 2; kmin0--; // Difference between Chebchev and B-spline kmax0--; k0--; for (k = kmin0; k <= kmax0; k++) { func[k-1] = B[k-k0]; dfunc[k-1] = -DB[k-k0] * RADIUS; ddfunc[k-1] = DDB[k-k0] * RADIUS * RADIUS; } return; default: fprintf(stderr, "unknown radial function of dradfun! iopt=%d\n", iopt); exitLog("unknown radial function of dradfun!\n"); } /* compute Chebyshev polynomials (iopt=14 or 15 or 16) */ u2 = 2.0 *u; func[0] = 1.0; dfunc[0] = 0.0; ddfunc[0] = 0.0; func[1] = u; dfunc[1] = 1.0; ddfunc[1] = 0.0; for (i = 2; i <= kmax; i++) { double tmpm1, dtmpm1; int im1, im2; im1 = i - 1; im2 = i - 2; tmpm1 = func[im1]; dtmpm1 = dfunc[im1]; func[i] = u2 * tmpm1 - func[im2]; dfunc[i] = u2 * dtmpm1 + 2.0 * tmpm1 - dfunc[im2]; ddfunc[i] = u2 * ddfunc[im1] + 4.0 * dtmpm1 - ddfunc[im2]; } u2x = u2r*RADIUS; u2xsq = u2x*u2x; for (i = 0; i <= kmax; i++) { if (i <= maxknot) { double tmpc = conc[i]; func[i] *= tmpc; dfunc[i] *= tmpc * u2x; ddfunc[i] *= tmpc * u2xsq; } else { dfunc[i] *= u2x; ddfunc[i] *= u2xsq; } } } /* file radfun.c * Returns values of radial functions at specified radius. * For some iopt's, it computes all radial orders (e.g. 4 for M84A) * even though the calling code may only use a single radial order * at a time. * * * input: * iopt = choice of radial functions * 1 = Legendre polynomials in upper mantle, as in M84A * 2 = " " " lower " , as in LO2.56 * 3 = constants, 3 layer upper mantle * 4 = constants, crustal portion of M84C * 5 = constant, single layer upper mantle * 6 = single layer crust, from Smith and Masters * 7 = Tanimoto model, 10 depth knots, lin. interp. between knots * 8 = 4 knots, linear interpolation between knots * 9 = Tanimoto whole mantle model, 11 layers * 10= Inoue et al. whole mantle model * 11= 4 layer (5 knot) upper mantle (lin. interp. between knots) * 12= Legendre polynomial l. m. with D" layer * 13= 4 knots, linear interpolation between knots (PREM crust) * 14= Chebyshev upper mantle k=0 to 4 * 15= Chebyshev lower mantle k=0 to 8 * 16= Chebyshev whole mantle k=0 to 8 (increased to 13) * 17= 14 layer whole mantle * 18= 14 layer upper mantle (Zhang's model - reduced from 26 layers) * 19= experimental * r = radius (km) * output: * func[0:kmax] = radial function values */ void radfun(RADFUN_T iopt, double *func, double r, int kmax) { double u, u2; double *xnode, B[4]; const int maxknot = 20; int i, k0, k; int kmin0, kmax0; // We assume that any conc[] beyond maxknot is 1 static double conc[maxknot] = { 0.70710678118655, 1.2247448713916, 1.0350983390135, 1.0145993123918, 1.0080322575484, 1.0050890913907, 1.0035149493262, 1.0025740068320, 1.0019665702378, 1.0015515913133, 1.0012554932754, 1.0010368069141, 1.0008707010792, 1.0007415648034, 1.0, 1.0, // not exactly right 1.0, 1.0, 1.0, 1.0 }; for (k = 0; k <= kmax; k++) func[k] = 0.0; if (kmax < 0) { fprintf(stderr, "radial function order underflow! kmax=%d\n", kmax); exitLog("radial function order underflow!\n"); } switch (iopt) { case CHEB_UPPER: /* Chebyshev polynomial in upper mantle */ if (r > RMOHO || r < R670) return; u = (r + r - RMOHO - R670) / (RMOHO - R670); break; case CHEB_LOWWER: /* Chebyshev polynomial in lower mantle */ if (r > R670 || r < RCMB) return; u = (r + r - R670 - RCMB) / (R670 - RCMB); break; case CHEB_WHOLE: /* Chebyshev polynomial for whole mantle */ if (r > RMOHO || r < RCMB) return; u = (r + r - RMOHO - RCMB) / (RMOHO - RCMB); break; case B_SPLINE: if (r > RMOHO || r < RCMB) return; xnode = vector(-1, kmax+2); SetNode(xnode); k0 = BSplKer(B, RADIUS - r, xnode, kmax); free_vector(xnode, -1, kmax+2); if (k0 <= 0 || k0 >= kmax) return; if (k0 >= 2) kmin0 = k0; else kmin0 = k0 + 1; if (k0 != (kmax-1)) kmax0 = k0 + 3; else kmax0 = k0 + 2; kmin0--; // Difference between Chebchev and B-spline kmax0--; k0--; for (k = kmin0; k <= kmax0; k++) func[k-1] = B[k-k0]; return; default: fprintf(stderr, "unknown radial function! iopt= %d\n", iopt); exitLog("unknown radial function!\n"); } /* compute Chebyshev polynomials (iopt=14 or 15 or 16) */ u2 = 2.0 *u; func[0] = 1.0; func[1] = u; for (i = 2; i <= kmax; i++) func[i] = u2 * func[i-1] - func[i-2]; for (i = 0; i <= kmax; i++) if (i <= maxknot) func[i] *= conc[i]; return; } void SetNode(double xnode[]) { xnode[1] = 24.40; xnode[2] = 75.00; xnode[3] = 125.00; xnode[4] = 225.00; xnode[5] = 350.00; xnode[6] = 500.00; xnode[7] = 670.00; xnode[8] = 820.00; xnode[9] = 1320.00; xnode[10] = 1820.00; xnode[11] = 2320.00; xnode[12] = 2550.00; xnode[13] = 2791.00; xnode[14] = 2891.00; xnode[-1] = xnode[0] = xnode[1]; xnode[16] = xnode[15] = xnode[14]; return; } void SetSplitNode(double xnodeUM[], double xnodeLM[]) { xnodeUM[1] = 24.40; xnodeUM[2] = 100.0; xnodeUM[3] = 225.00; xnodeUM[4] = 350.00; xnodeUM[5] = 500.00; xnodeUM[6] = 670.00; xnodeUM[-1] = xnodeUM[0] = xnodeUM[1]; xnodeUM[8] = xnodeUM[7] = xnodeUM[6]; xnodeLM[1] = 670.00; xnodeLM[2] = 820.00; xnodeLM[3] = 1320.00; xnodeLM[4] = 1820.00; xnodeLM[5] = 2320.00; xnodeLM[6] = 2550.00; xnodeLM[7] = 2791.00; xnodeLM[8] = 2891.00; xnodeLM[-1] = xnodeLM[0] = xnodeLM[1]; xnodeLM[10] = xnodeLM[9] = xnodeLM[8]; return; } /* Returns values of radial functions at specified radius. * For some iopt's, it computes all radial orders (e.g. 4 for M84A) * even though the calling code may only use a single radial order * at a time. */ /* * input: * iopt = choice of radial functions * 1 = Legendre polynomials in upper mantle, as in M84A * 2 = " " " lower " , as in LO2.56 * 3 = constants, 3 layer upper mantle * 4 = constants, crustal portion of M84C * 5 = constant, single layer upper mantle * 6 = single layer crust, from Smith and Masters * 7 = Tanimoto model, 10 depth knots, lin. interp. between knots * 8 = 4 knots, linear interpolation between knots * 9 = Tanimoto whole mantle model, 11 layers * 10= Inoue et al. whole mantle model * 11= 4 layer (5 knot) upper mantle (lin. interp. between knots) * 12= Legendre polynomial l. m. with D" layer * 13= 4 knots, linear interpolation between knots (PREM crust) * 14= Chebyshev upper mantle k=0 to 4 * 15= Chebyshev lower mantle k=0 to 8 * 16= Chebyshev whole mantle k=0 to 8 (increased to 13) * 17= 14 layer whole mantle * 18= 14 layer upper mantle (Zhang's model - reduced from 26 layers) * 19= experimental * r = radius (km) * output: * func[0:kmax] = radial function values * dfunc[0:kmax] = radial function derivatives */ void dradfun(RADFUN_T iopt, double *func, double *dfunc, double r, int kmax) { double u, u2r, u2, u2x; double *xnode, B[4], DB[4]; const int maxknot = 20; int i, k0, k; int kmin0, kmax0; static double conc[maxknot] = { 0.70710678118655, 1.2247448713916, 1.0350983390135, 1.0145993123918, 1.0080322575484, 1.0050890913907, 1.0035149493262, 1.0025740068320, 1.0019665702378, 1.0015515913133, 1.0012554932754, 1.0010368069141, 1.0008707010792, 1.0007415648034, 1.0, 1.0, // not exactly right 1.0, 1.0, 1.0, 1.0 }; for (k = 0; k <= kmax; k++) dfunc[k] = func[k] = 0.0; if (kmax > maxknot) { fprintf(stderr, "radial function order overflow! kmax=%d\n", kmax); exitLog("radial function order overflow!\n"); } switch (iopt) { case CHEB_UPPER: /* Chebyshev polynomial in upper mantle */ if (r > RMOHO || r < R670) return; u = (r + r - RMOHO - R670) / (RMOHO - R670); u2r = 2.0 / (RMOHO - R670); break; case CHEB_LOWWER: /* Chebyshev polynomial in lower mantle */ if (r > R670 || r < RCMB) return; u = (r + r - R670 - RCMB) / (R670 - RCMB); u2r = 2.0 / (R670 - RCMB); break; case CHEB_WHOLE: /* Chebyshev polynomial for whole mantle */ if (r > RMOHO || r < RCMB) return; u = (r + r - RMOHO - RCMB) / (RMOHO - RCMB); u2r = 2.0 / (RMOHO - RCMB); break; case B_SPLINE: if (r > RMOHO || r < RCMB) return; xnode = vector(-1, kmax+2); SetNode(xnode); k0 = DBBSplKer(B, DB, RADIUS - r, xnode, kmax); free_vector(xnode, -1, kmax+2); if (k0 <= 0 || k0 >= kmax) return; if (k0 >= 2) kmin0 = k0; else kmin0 = k0 + 1; if (k0 != (kmax-1)) kmax0 = k0 + 3; else kmax0 = k0 + 2; kmin0--; // Difference between Chebchev and B-spline kmax0--; k0--; for (k = kmin0; k <= kmax0; k++) { func[k-1] = B[k-k0]; dfunc[k-1] = -DB[k-k0] * RADIUS; } return; default: fprintf(stderr, "unknown radial function of dradfun! iopt=%d\n", iopt); exitLog("unknown radial function of dradfun!\n"); } /* compute Chebyshev polynomials (iopt=14 or 15 or 16) */ u2 = 2.0 *u; func[0] = 1.0; dfunc[0] = 0.0; func[1] = u; dfunc[1] = 1.0; for (i = 2; i <= kmax; i++) { int im1 = i - 1; int im2 = i - 2; double tmpm1 = func[im1]; func[i] = u2 * tmpm1 - func[im2]; dfunc[i] = u2 * dfunc[im1] + 2.0 * tmpm1 - dfunc[im2]; } u2x = u2r*RADIUS; for (i = 0; i <= kmax; i++) { if (i <= maxknot) { double tmpc = conc[i]; func[i] *= tmpc; dfunc[i] *= tmpc * u2x; } else { dfunc[i] *= u2x; } } } /* The subroutine rotate_velocity calculates the phase velocity expansion coefficients ap and bp in the equatorial reference frame from the expansion coefficients a and b in the geographical reference frame for a given triplet of Euler angles alpha, beta and gamma. */ void rotate_velocity(double alpha, double beta, double gamma, int nl, double **a, double **b, double **ap, double **bp) { int l, m, mp; double **c_r, **c_i; double **cp_r, **cp_i; double Dmmp[2]; double alp = alpha, bet = beta, gam = gamma; c_r = matrix(0, nl, 0, 2 * nl); c_i = matrix(0, nl, 0, 2 * nl); cp_r = matrix(0, nl, 0, 2 * nl); cp_i = matrix(0, nl, 0, 2 * nl); convert_coef(nl, a, b, c_r, c_i); /* Calculate complex coefficients c'_l^m in the equatorial reference frame from the complex coefficients c_l^m in the geographic frame: c'_l^m = Sum_m' D_mm' c_l^m' This is basically the fourth equation on page 25 of Edmonds (1960). */ for (l = 0; l <= nl; l++) { for (m = -l; m <= l; m++) { cp_r[l][l+m] = cp_i[l][l+m] = 0.0; for (mp = -l; mp <= l; mp++) { rotylm(alp, bet, gam, l, m, mp, Dmmp); cp_r[l][l+m] += c_r[l][l+mp] * Dmmp[0] - c_i[l][l+mp] * Dmmp[1]; cp_i[l][l+m] += c_r[l][l+mp] * Dmmp[1] + c_i[l][l+mp] * Dmmp[0]; } } } /* Switch back to real coefficients a'_l^m and b'_l^m: a'_l^0 = Re c'_l^0 b'_l^0 = 0 a'_l^m = 2 Re c'_l^0, m>0, b'_l^m = -2 Im c'_l^0, m>0. */ for (l = 0; l <= nl; l++) { m = 0; ap[l][m] = cp_r[l][l+m]; bp[l][m] = 0.0; for (m = 1; m <= l; m++) { ap[l][m] = 2.0 * cp_r[l][l+m]; bp[l][m] = -2.0 * cp_i[l][l+m]; } } free_matrix(c_r, 0, nl, 0, 2 * nl); free_matrix(c_i, 0, nl, 0, 2 * nl); free_matrix(cp_r, 0, nl, 0, 2 * nl); free_matrix(cp_i, 0, nl, 0, 2 * nl); } /* Switch from real coefficients a_l^m and b_l^m to complex coefficients c_l^m. We have c_l^0 = a_l^0 c_l^m = (a_l^m - i b_l^m)/2, m>0, c_l^-m = (-1)^m (a_l^m + i b_l^m)/2, m<0. */ void convert_coef(int nl, double **a, double **b, double **c_r, double **c_i) { int l, m; for (l = 0; l <= nl; l++) { for (m = -l; m <= l; m++) { if (m != 0) { if (m > 0) { c_r[l][l+m] = 0.5 * a[l][m]; c_i[l][l+m] = -0.5 * b[l][m]; } else { if (EVEN(-m)) { c_r[l][l+m] = 0.5 * a[l][-m]; c_i[l][l+m] = 0.5 * b[l][-m]; } else { c_r[l][l+m] = -0.5 * a[l][-m]; c_i[l][l+m] = -0.5 * b[l][-m]; } } } else { c_r[l][l+m] = a[l][m]; c_i[l][l+m] = 0.0; } } } } /* Function BSplKer-- * to calculate cubic B-spline. * written by Wei-jia Su (in wslib.c) * modified by Xian-feng Liu * Input: x (given depth or radius); * xnode[-1:N+2]: depth (or radius) on each knot; * N: number of knots; * Output: B[0],B[1], B[2], B[3] spline coefficients. * Also return the knot value (beginning from 1). * Warning: This function uses static values which * is bad for two different kinds of spline coefficients. * */ int BSplKer(double *B, double x, double *xnode, int N) { double *xn, temp, temp1, X1, X2, X3, X4; int i, ib, ie; static double A1, A2, A3, A4; static int i0 = 0; Boolean increase; if (xnode[N] < xnode[1]) increase = FALSE; else increase = TRUE; if (((x < xnode[1] || x > xnode[N]) && increase) || ((x < xnode[N] || x > xnode[1]) && !increase)) { B[0] = B[1] = B[2] = B[3] = 0.0; return 0; } ib = 1; ie = N; while (ie - ib > 1) { i = ib + ((ie - ib) >> 1); if ((x < xnode[i] && increase) || (x > xnode[i] && !increase)) ie = i; else ib = i; } i = ib; xn = &xnode[i-2]; if (i != i0) { temp = 1.0 / (xn[3] - xn[2]); A1 = temp / (xn[3] - xn[1]); A3 = temp / (xn[4] - xn[2]); A4 = A3 / (xn[5] - xn[2]); temp = xn[4] - xn[1]; A2 = A1 / temp; A3 /= temp; A1 /= xn[3] - xn[0]; i0 = i; } X1 = x - xn[1]; X2 = x - xn[2]; X3 = x - xn[3]; X4 = x - xn[4]; temp = X3 * X3 * A1; temp1 = X3 * X1 * A2 + X2 * X4 * A3; B[0] = -temp * X3; B[1] = (x - xn[0]) * temp + X4 * temp1; temp = X2 * X2 * A4; B[3] = X2 * temp; B[2] = -X1 * temp1 - (x - xn[5]) * temp; /* * Imposing boundary conditions: * for 2nd derivative = 0 */ if (i == 1) { B[2] += B[1] / 3.0; B[1] = B[1] / 1.5 + B[0]; B[0] = 0.0; } else if (i == 2) { B[1] += B[0] / 3.0; B[0] /= 1.5; } else if (i == N - 1) { B[1] += B[2] / 3.0; B[2] = B[2] / 1.5 + B[3]; B[3] = 0.0; } else if (i == N - 2) { B[2] += B[3] / 3.0; B[3] /= 1.5; } return i; } /* * to calculate cubic B-spline Kernel and its derivative */ int DBBSplKer(double *B, double *DB, double x, double *xnode, int N) { static double A1, A2, A3, A4; double *xn, temp, X0, X1, X2, X3, X4, X5; double XX1, XX2, XX3, XX4; static int i0 = 0; int i, ib, ie; Boolean increase; if (xnode[N] < xnode[1]) increase = FALSE; else increase = TRUE; if (((x < xnode[1] || x > xnode[N]) && increase) || ((x < xnode[N] || x > xnode[1]) && !increase)) { B[0] = B[1] = B[2] = B[3] = 0.0; DB[0] = DB[1] = DB[2] = DB[3] = 0.0; return 0; } ib = 1; ie = N; while (ie - ib > 1) { i = ib + ((ie - ib) >> 1); if ((x < xnode[i] && increase) || (x > xnode[i] && !increase)) ie = i; else ib = i; } i = ib; xn = &xnode[i-2]; if (i != i0) { temp = 1.0 / (xn[3] - xn[2]); A1 = temp / (xn[3] - xn[1]); A3 = temp / (xn[4] - xn[2]); A4 = A3 / (xn[5] - xn[2]); temp = xn[4] - xn[1]; A2 = A1 / temp; A3 /= temp; A1 /= xn[3] - xn[0]; i0 = i; } X0 = x - xn[0]; X1 = x - xn[1]; X2 = x - xn[2]; X3 = x - xn[3]; X4 = x - xn[4]; X5 = x - xn[5]; XX1 = X1 * X1; XX2 = X2 * X2; XX3 = X3 * X3; XX4 = X4 * X4; temp = XX3 * A1; B[0] = -temp * X3; B[1] = X0 * temp; DB[0] = -3.0 * temp; DB[1] = temp; temp = XX2 * A4; B[3] = X2 * temp; DB[3] = 3.0 * temp; B[2] = -X5 * temp; DB[2] = -temp; temp = XX1 * A2; B[2] -= X3 * temp; DB[2] -= temp; temp = XX4 * A3; B[1] += X2 * temp; DB[1] += temp; temp = X1 * X3 * A2; B[1] += X4 * temp; DB[1] += temp; DB[2] -= temp + temp; temp = X1 * X4; DB[1] += temp * A2; temp *= A3; B[2] -= X2 * temp; DB[2] -= temp; temp = X2 * X4 * A3; DB[1] += 2.0 * (temp + X0 * X3 * A1) + X3 * X4 * A2; DB[2] -= temp + X2 * (X1 * A3 + 2.0 * X5 * A4); /* Imposing boundary conditions: * for 2nd derivative = 0 */ if (i == 1) { B[2] += B[1] / 3.0; B[1] = B[1] / 1.5 + B[0]; DB[2] += DB[1] / 3.0; DB[1] = DB[1] / 1.5 + DB[0]; B[0] = DB[0] = 0.0; } else if (i == 2) { B[1] += B[0] / 3.0; B[0] /= 1.5; DB[1] += DB[0] / 3.0; DB[0] /= 1.5; } else if ( i == N - 1) { B[1] += B[2] / 3.0; B[2] = B[2] / 1.5 + B[3]; DB[1] += DB[2] / 3.0; DB[2] = DB[2] / 1.5 + DB[3]; B[3] = DB[3] = 0.0; } else if (i == N - 2) { B[2] += B[3] / 3.0; B[3] /= 1.5; DB[2] += DB[3] / 3.0; DB[3] /= 1.5; } return i; } /* * to calculate cubic B-spline Kernel and its derivative */ int DDBBSplKer(double *B, double *DB, double *DDB, double x, double *xnode, int N) { static double A1, A2, A3, A4; double *xn, temp, X0, X1, X2, X3, X4, X5; double XX1, XX2, XX3, XX4; static int i0 = 0; int i, ib, ie; Boolean increase; if (xnode[N] < xnode[1]) increase = FALSE; else increase = TRUE; if (((x < xnode[1] || x > xnode[N]) && increase) || ((x < xnode[N] || x > xnode[1]) && !increase)) { B[0] = B[1] = B[2] = B[3] = 0.0; DB[0] = DB[1] = DB[2] = DB[3] = 0.0; DDB[0] = DDB[1] = DDB[2] = DDB[3] = 0.0; return 0; } ib = 1; ie = N; while (ie - ib > 1) { i = ib + ((ie - ib) >> 1); if ((x < xnode[i] && increase) || (x > xnode[i] && !increase)) ie = i; else ib = i; } i = ib; xn = &xnode[i-2]; if (i != i0) { temp = 1.0 / (xn[3] - xn[2]); A1 = temp / (xn[3] - xn[1]); A3 = temp / (xn[4] - xn[2]); A4 = A3 / (xn[5] - xn[2]); temp = xn[4] - xn[1]; A2 = A1 / temp; A3 /= temp; A1 /= xn[3] - xn[0]; i0 = i; } X0 = x - xn[0]; X1 = x - xn[1]; X2 = x - xn[2]; X3 = x - xn[3]; X4 = x - xn[4]; X5 = x - xn[5]; XX1 = X1 * X1; XX2 = X2 * X2; XX3 = X3 * X3; XX4 = X4 * X4; temp = XX3 * A1; B[0] = -temp * X3; B[1] = X0 * temp; DB[0] = -3.0 * temp; DB[1] = temp; temp = XX2 * A4; B[3] = X2 * temp; DB[3] = 3.0 * temp; B[2] = -X5 * temp; DB[2] = -temp; temp = XX1 * A2; B[2] -= X3 * temp; DB[2] -= temp; temp = XX4 * A3; B[1] += X2 * temp; DB[1] += temp; temp = X1 * X3 * A2; B[1] += X4 * temp; DB[1] += temp; DB[2] -= 2.0 * temp; temp = X1 * X4; DB[1] += temp * A2; temp *= A3; B[2] -= X2 * temp; DB[2] -= temp; temp = X2 * X4 * A3; DB[1] += 2.0 * (temp + X0 * X3 * A1) + X3 * X4 * A2; DB[2] -= temp + X2 * (X1 * A3 + 2.0 * X5 * A4); /* Second derivative */ DDB[0] = -6.0 * X3 * A1; DDB[1] = 2.0 * (X3 + X3 + X0) * A1 + 2.0 * (X1 + X3 + X4) * A2 + 2.0 * (X2 + X4 + X4) * A3; DDB[2] = -2.0 * (X1 + X1 + X3) * A2 - 2.0 * (X1 + X2 + X4) * A3 -2.0 * (X2 + X2 + X5) * A4; DDB[3] = 6.0 * X2 * A4; /* Imposing boundary conditions: * for 2nd derivative = 0 */ if (i == 1) { B[2] += B[1] / 3.0; B[1] = B[1] / 1.5 + B[0]; DB[2] += DB[1] / 3.0; DB[1] = DB[1] / 1.5 + DB[0]; DDB[2] += DDB[1] / 3.0; DDB[1] = DDB[1] / 1.5 + DDB[0]; B[0] = DB[0] = DDB[0] = 0.0; } else if (i == 2) { B[1] += B[0] / 3.0; B[0] /= 1.5; DB[1] += DB[0] / 3.0; DB[0] /= 1.5; DDB[1] += DDB[0] / 3.0; DDB[0] /= 1.5; } else if (i == N - 1) { B[1] += B[2] / 3.0; B[2] = B[2] / 1.5 + B[3]; DB[1] += DB[2] / 3.0; DB[2] = DB[2] / 1.5 + DB[3]; DDB[1] += DDB[2] / 3.0; DDB[2] = DDB[2] / 1.5 + DDB[3]; B[3] = DB[3] = DDB[3] = 0.0; } else if (i == N - 2) { B[2] += B[3] / 3.0; B[3] /= 1.5; DB[2] += DB[3] / 3.0; DB[3] /= 1.5; DDB[2] += DDB[3] / 3.0; DDB[3] /= 1.5; } return i; } /* Function SumSphBspline-- * written by Wei-jia Su * modified by Xian-feng Liu * Input: spherical spline model: model; * summation of radial coefficients at given depth: ar[1:n_node]; * theta, phi: spherical coordinates. * Output: summation (return value). */ double SPH_SPLINE_T::SumSphBspline(double *ar, double theta, double phi) const { double val, delta; int i, node, n, jj, k, j, m; double theta2, phi2; double costh, sinth; char *pbits; int max_char; BSpline bSpline(avg_dist); /* printf("%f %d %f %f\n",h,n_node,theta,phi); */ max_char = n_node + 1; pbits = (char *) malloc(max_char * sizeof(char)); memset(pbits, 0, max_char); /* clear set */ val = 0.0; node = search_node(theta, phi); sinth = sin(theta); costh = cos(theta); n = aj_node[node][0]; /* number of neighboring nodes */ for (jj = 1; jj <= n; jj++) { j = aj_node[node][jj]; /* index of neighboring nodes */ m = aj_node[j][0]; /* number of neighboring nodes again */ for (k = 1; k <= m; k++) { i = aj_node[j][k]; /* index of neighboring nodes again */ if (pbits[i]) /* already calculated */ continue; else /* uncalculated */ pbits[i] = 1; /* set the flag (calculated) */ theta2 = co_node[i][3]; phi2 = co_node[i][4]; delta = ARC_DIST(costh, sinth, phi, theta2, phi2); val += bSpline.fdelta(delta) * ar[i]; } } free(pbits); return val; } /* Function d_SumSphBspline-- * modified by Xian-feng Liu from SumSphBspline; * Input: spherical spline model: model; * summation of radial coefficients at given depth: ar[1:n_node]; * theta, phi: spherical coordinates. * Output: summation (return value); * dv_dx: dv/dx; * dv_delta: dv/ddelta. */ double SPH_SPLINE_T::d_SumSphBspline(double *ar, double *ar_d, double theta, double phi, double *dv_dx, double *dv_dth, double *dv_dphi) const { double val, delta; int i, node, n, jj, k, j, m; double theta2, phi2, ar_i, dv_delta, sin_delta; double costh2, sinth2, cosph_diff, cos_arc; double ddelta_dth, ddelta_dphi, dvdr, dvdth, dvdphi; double costh, sinth, fdel, fdel_d; char *pbits; int max_char; BSpline bSpline(avg_dist); max_char = n_node + 1; pbits = (char *) malloc(max_char * sizeof(char)); memset(pbits, 0, max_char); /* clear set */ val = 0.0; dvdr = 0.0; dvdth = 0.0; dvdphi = 0.0; node = search_node(theta, phi); sinth = sin(theta); costh = cos(theta); int *aj_node_node = aj_node[node]; n = aj_node_node[0]; /* number of neighboring nodes */ int *aj_node_j; double *co_node_i; double phi_diff; for (jj = 1; jj <= n; jj++) { j = aj_node_node[jj]; /* index of neighboring nodes */ aj_node_j = aj_node[j]; m = aj_node_j[0]; /* number of neighboring nodes again */ for (k = 1; k <= m; k++) { i = aj_node_j[k]; /* index of neighboring nodes again */ if (pbits[i]) continue; else pbits[i] = 1; co_node_i = co_node[i]; theta2 = co_node_i[3]; phi2 = co_node_i[4]; costh2 = cos(theta2); sinth2 = sin(theta2); phi_diff = phi - phi2; cosph_diff = cos(phi_diff); cos_arc = costh * costh2 + sinth * sinth2 * cosph_diff; delta = acos(cos_arc); // fdel = fdelta(delta, h); fdel = bSpline.f_df_delta(delta, &fdel_d); ar_i = ar[i]; val += fdel * ar_i; dvdr += fdel * ar_d[i]; dv_delta = fdel_d * ar_i; sin_delta = sin(delta); ddelta_dth = (sinth * costh2 - costh * sinth2 * cosph_diff) / sin_delta; ddelta_dphi = sinth * sinth2 * sin(phi_diff) / sin_delta; dvdth += dv_delta * ddelta_dth; dvdphi += dv_delta * ddelta_dphi; } } free(pbits); *dv_dx = -dvdr * RADIUS; *dv_dth = dvdth; *dv_dphi = dvdphi; return val; } /* Function dd_SumSphBspline-- * modified by Xian-feng Liu from SumSphBspline; * Input: spherical spline model: model; * summation of radial coefficients at given depth: ar[1:n_node]; * theta, phi: spherical coordinates. * Output: summation (return value); * dv_r: dv/dr; * dv_delta: dv/ddelta. * dv_dx =d(dv)/dx (x=r/6371, normalized radius) * dv_dth =d(dv)/dtheta * dv_dphi =d(dv)/dphi * ddv_d2x =d^2(dv)/dx^2 * ddv_dthdth =d^2(dv)/dtheta^2 * ddv_dxth =d^2(dv)/dxdtheta * ddv_dxdph =d^2(dv)/dxdphi * ddv_dthdph =d^2(dv)/dthetadphi */ double SPH_SPLINE_T::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 { double val, delta; int i, node, n, jj, k, j, m; double theta2, phi2, ar_i, dv_delta, sinx; double costh2, sinth2, cosph_diff, cosx; double ddelta_dth, ddelta_dph, dvdr, dvdth, dvdphi; double cotx, sinph_diff, fdel, fdel_dd; double d2v_d2r, d2v_d2th, d2v_dthdph; double d2v_drdph, d2v_drdth, fdel_d, ar_d_i; double d2v_delta, d2x_d2th, d2x_dthdph; double d2delta_d2th, d2delta_dthdph; double costh, sinth; char *pbits; int max_char; BSpline bSpline(avg_dist); max_char = n_node + 1; pbits = (char *) malloc(max_char * sizeof(char)); memset(pbits, 0, max_char); /* clear set */ val = 0.0; dvdr = 0.0; d2v_d2r = 0.0; dvdth = 0.0; dvdphi = 0.0; d2v_d2th = 0.0; d2v_dthdph = 0.0; d2v_drdph = 0.0; d2v_drdth = 0.0; node = search_node(theta, phi); sinth = sin(theta); costh = cos(theta); n = aj_node[node][0]; /* number of neighboring nodes */ int *aj_node_node = aj_node[node]; int *aj_node_j; double *co_node_i, fdel_ar_d, tmp; double costh_sinth2, sinth_sinth2, d_th_cotx; for (jj = 1; jj <= n; jj++) { j = aj_node_node[jj]; /* index of neighboring nodes */ aj_node_j = aj_node[j]; m = aj_node_j[0]; /* number of neighboring nodes again */ for (k = 1; k <= m; k++) { i = aj_node_j[k]; /* index of neighboring nodes again */ if (pbits[i]) continue; else pbits[i] = 1; co_node_i = co_node[i]; theta2 = co_node_i[3]; phi2 = co_node_i[4]; costh2 = cos(theta2); sinth2 = sin(theta2); cosph_diff = cos(phi - phi2); cosx = costh * costh2 + sinth * sinth2 * cosph_diff; delta = acos(cosx); sinph_diff = sin(phi - phi2); fdel = bSpline.f_df_ddf_delta(delta, &fdel_d, &fdel_dd); ar_i = ar[i]; ar_d_i = ar_d[i]; val += fdel * ar_i; dvdr += fdel * ar_d_i; d2v_d2r += fdel * ar_dd[i]; dv_delta = fdel_d * ar_i; d2v_delta = fdel_dd * ar_i; sinx = sin(delta); cotx = cosx / sinx; costh_sinth2 = costh * sinth2; sinth_sinth2 = sinth * sinth2; ddelta_dth = (sinth * costh2 - costh_sinth2 * cosph_diff) / sinx; d2x_d2th = -costh * costh2 - sinth_sinth2 * cosph_diff; ddelta_dph = sinth_sinth2 * sinph_diff / sinx; d2x_dthdph = -costh_sinth2 * sinph_diff; d_th_cotx = ddelta_dth * cotx; d2delta_d2th = -d2x_d2th / sinx - ddelta_dth * d_th_cotx; d2delta_dthdph = -d2x_dthdph / sinx - ddelta_dph * d_th_cotx; dvdth += dv_delta * ddelta_dth; dvdphi += dv_delta * ddelta_dph; tmp = d2v_delta * ddelta_dth; d2v_d2th += tmp * ddelta_dth + dv_delta * d2delta_d2th; d2v_dthdph += tmp * ddelta_dph + dv_delta * d2delta_dthdph; fdel_ar_d = fdel_d * ar_d_i; d2v_drdth += fdel_ar_d * ddelta_dth; d2v_drdph += fdel_ar_d * ddelta_dph; } } free(pbits); *dv_dx = -dvdr * RADIUS; *dv_dth = dvdth; *dv_dphi = dvdphi; *ddv_d2x = d2v_d2r * RADIUS * RADIUS; *ddv_dthth = d2v_d2th; *ddv_dthdph = d2v_dthdph; *ddv_dxdth = -d2v_drdth * RADIUS; *ddv_dxdph = -d2v_drdph * RADIUS; return val; } /* search_node -- * an efficient way to search the node for spherical spline model */ int SPH_SPLINE_T::search_node(double theta, double phi) const { int n_strip = (int) sqrt((n_node - 2) / 10); double th_step = PI / (double) (3 * n_strip); int th_index; int ind_min, total_phi; double h = avg_dist; double start_phi, step_phi; rot(-gamma, -beta, -alpha, &theta, &phi); th_index = (int) (theta / th_step + 0.5); if (theta < 0.7 * h) return 1; if (PI - theta < 0.7 * h) return n_node; ind_min = th_to_index(n_strip, th_index, &total_phi); start_phi = co_node[ind_min][2]; step_phi = TWO_PI / total_phi; if ((phi - start_phi) > (0.5 * step_phi)) start_phi += TWO_PI; ind_min += (int) ((phi - start_phi) / step_phi + 0.5); return ind_min; } /* function th_to_index-- * utility function for search_node. */ int th_to_index(int n_strip, int th_index, int *total_phi) { int n2_strip = n_strip + n_strip; int n3_strip = n2_strip + n_strip; int rest; int ind = 1; if (th_index > n2_strip) { rest = th_index - n2_strip; ind += (n_strip * rest - rest * (rest + 1) / 2) * 5; *total_phi = 5 * (n3_strip - th_index); } if (th_index > n_strip) { if (th_index > n2_strip) rest = n_strip; else { rest = th_index - n_strip; *total_phi = 5 * n_strip; } ind += n_strip * rest * 5; } if (th_index > n_strip) rest = n_strip; else { rest = th_index; *total_phi = 5 * th_index; } ind += 5 * (rest + 1) * rest / 2; return ind; } void BSP_WHOLE_T::convert_to_complex() { cr = f3tensor(1, n_rad, 0, Lmax, 0, 2*Lmax); ci = f3tensor(1, n_rad, 0, Lmax, 0, 2*Lmax); if (cr == NULL || ci == NULL) { cerr << "BSP_WHOLE_T: no memory for 3D coefficients" << endl; exit(1); } double **ak = matrix(0, Lmax, 0, Lmax); double **bk = matrix(0, Lmax, 0, Lmax); if (ak == NULL || bk == NULL) { cerr << "BSP_WHOLE_T: no memory for 2D coefficients" << endl; exit(1); } int k, L, m; for (k = 1; k <= n_rad; k++) { ak[0][0] = a[0][0][k]; for (L = 1; L <= Lmax; L++) { for (m = 0; m <= L; m++) { ak[L][m] = a[L][m][k]; bk[L][m] = b[L][m][k]; } } convert_coef(Lmax, ak, bk, cr[k], ci[k]); } free_matrix(ak, 0, Lmax, 0, Lmax); free_matrix(bk, 0, Lmax, 0, Lmax); return; } // Function MakeModel3D serves as a virtual constructor for // abstract base class MODEL_3D. Remember that there is no virtual // constructor in the C++. Whenever you add a new kind of // derived class (i.e. a new kind of three-dimensional model) // from MODEL_3D, you have to change this function. I know // this is unfortunate, but there is probably no better way in C++. MODEL_3D *MakeModel3D(char *model3D_name, MODEL3D_T type, double alpha, double beta, double gamma, char *descrip) { switch (type) { case SPH_SPLINE_MODEL: return (new SPH_SPLINE_T(model3D_name, type, alpha, beta, gamma, descrip)); break; case XV_SPH_SPLINE_MODEL: return (new SPH_SPLINE_T(model3D_name, type, alpha, beta, gamma, descrip)); break; case CHEB_WHOLE_MODEL: return (new CHEB_WHOLE_T(model3D_name, type, alpha, beta, gamma, descrip)); break; case BSP_WHOLE_MODEL: return (new BSP_WHOLE_T(model3D_name, type, alpha, beta, gamma, descrip)); break; case BSP_SPLIT_MODEL: return (new BSP_SPLIT_T(model3D_name, type, alpha, beta, gamma, descrip)); break; case CHEB_SPLIT_MODEL: return (new CHEB_SPLIT_T(model3D_name, type, alpha, beta, gamma, descrip)); break; default: fprintf(stderr, "unknown model type.\n"); return NULL; } } // Function MakeSphericModel3D serves as a virtual constructor for // abstract base class SPHERIC_MODEL_3D. Remember that there is no // virtual constructor in the C++. Whenever you add a new kind of // derived class (i.e. a new kind of three-dimensional model) // from SPHERIC_MODEL_3D, you have to change this function. I know // this is unfortunate, but there is probably no better way in C++. // If the given model is not spherical harmonics, return NULL. SPHERIC_MODEL_3D *MakeSphericModel3D(char *model3D_name, MODEL3D_T type, double alpha, double beta, double gamma) { switch (type) { case CHEB_WHOLE_MODEL: return (new CHEB_WHOLE_T(model3D_name, type, alpha, beta, gamma)); break; case BSP_WHOLE_MODEL: return (new BSP_WHOLE_T(model3D_name, type, alpha, beta, gamma)); break; case BSP_SPLIT_MODEL: return (new BSP_SPLIT_T(model3D_name, type, alpha, beta, gamma)); break; case CHEB_SPLIT_MODEL: return (new CHEB_SPLIT_T(model3D_name, type, alpha, beta, gamma)); break; case SPH_SPLINE_MODEL: // Not a spherical harmonics model fprintf(stderr, "not spherical harmonics model.\n"); return NULL; break; case XV_SPH_SPLINE_MODEL: // Not a spherical harmonics model fprintf(stderr, "not spherical harmonics model.\n"); return NULL; break; default: fprintf(stderr, "unknown model type.\n"); return NULL; } } void SPHERIC_MODEL_3D::CreateAB(double xr) { int Lmax = GetLmax(xr); Alm = matrix(0, Lmax, 0, Lmax); Blm = matrix(0, Lmax, 0, Lmax); for (int L = 0; L <= Lmax; L++) for (int m = 0; m <= L; m++) Alm[L][m] = Blm[L][m] = 0.0; } // Free memory allocated by CreateAB. It is only used after CreateAB. void SPHERIC_MODEL_3D::FreeAB(double xr) { int Lmax = GetLmax(xr); free_matrix(Alm, 0, Lmax, 0, Lmax); free_matrix(Blm, 0, Lmax, 0, Lmax); } // Do the radial summation for spherical harmonics model. // The result is in coeffcients Alm and Blm. void BSP_SPLIT_T::SumRad(double xr) { double rr = xr * RADIUS; double *co_rad; RadialBSpline *radialBSpline; int Lmax, n_rad; double ***a, ***b; int m, L, k; Lmax = GetLmax(xr); if (rr >= R670) { // upper mantle co_rad = co_radUM; n_rad = n_radUM; a = aUM; b = bUM; radialBSpline = radialBSplineUM; } else { // Lower mantle co_rad = co_radLM; n_rad = n_radLM; a = aLM; b = bLM; radialBSpline = radialBSplineLM; } if (rr < RCMB || rr > RMOHO) return; double depth = RADIUS - rr; double B[4]; k = radialBSpline->BSplKer(B, depth); //printf("k=%d B=%f %f %f %f\n", k, B[0], B[1], B[2], B[3]); if (k <= 0 || k >= n_rad) return; CreateAB(xr); double *arl, *brl; int km1 = k - 1; int kp1 = k + 1; int kp2 = k + 2; double B0 = B[0]; double B1 = B[1]; double B2 = B[2]; double B3 = B[3]; double *alm, *blm, **al, **bl; // for efficiency double Almk, Blmk, Almkp1, Blmkp1; // FILE *fp = fopen("ar.dat", "w"); for (L = 0; L <= Lmax; L++) { al = a[L]; bl = b[L]; arl = Alm[L]; brl = Blm[L]; for (m = 0; m <= L; m++) { alm = al[m]; blm = bl[m]; Almk = alm[k]; Blmk = blm[k]; Almkp1 = alm[kp1]; Blmkp1 = blm[kp1]; arl[m] = B1 * Almk + B2 * Almkp1; brl[m] = B1 * Blmk + B2 * Blmkp1; if (k != 1) { arl[m] += B0 * alm[km1]; brl[m] += B0 * blm[km1]; } if (kp1 != n_rad) { arl[m] += B3 * alm[kp2]; brl[m] += B3 * blm[kp2]; } // fprintf(fp, "%d %d %f %f\n", L, m, arl[m], brl[m]); } } // fclose(fp); } // Given two spherical harmonics model at radius xr (normalized), // this function returns its cross-correlation coefficient at // this radius and different degrees. The overall correlation coefficient // is put at the end of the returned vector. corr[L] represents the // cross-correlation coefficient for each degree. // A friend of abstract data class SPHERIC_MODEL_3D. double *Correlation(SPHERIC_MODEL_3D *m1, SPHERIC_MODEL_3D *m2, double xr) { if (xr < 0.0 || xr > 1.0) { fprintf(stderr, "Corrleation: radius = %f\n", xr * RADIUS); exit(1); } m1->SumRad(xr); m2->SumRad(xr); int Lmax = m1->GetLmax(xr); int Lmax2 = m2->GetLmax(xr); if (Lmax2 < Lmax) // Only comparison between common degrees Lmax = Lmax2; double *corr = new double[Lmax+2]; double totalpower1 = 0.0; double totalpower2 = 0.0; double totalpower12 = 0.0; double a1, b1, a2, b2; for (int L = 0; L <= Lmax; L++) { double power1 = 0.0; double power2 = 0.0; double power12 = 0.0; for (int m = 0; m <= L; m++) { a1 = m1->Alm[L][m]; b1 = m1->Blm[L][m]; a2 = m2->Alm[L][m]; b2 = m2->Blm[L][m]; power1 += a1 * a1 + b1 * b1; power2 += a2 * a2 + b2 * b2; power12 += a1 * a2 + b1 * b2; totalpower1 += a1 * a1 + b1 * b1; totalpower2 += a2 * a2 + b2 * b2; totalpower12 += a1 * a2 + b1 * b2; } corr[L] = (power12 / sqrt(power1 * power2)); } corr[Lmax+1] = totalpower12 / sqrt(totalpower1 * totalpower2); m1->FreeAB(xr); m2->FreeAB(xr); return corr; } // Given two spherical harmonics model at different // radius xr (normalized), // this function returns its cross-correlation coefficient at // this radius. // A friend of abstract data class SPHERIC_MODEL_3D. double *Correlation(SPHERIC_MODEL_3D *m1, SPHERIC_MODEL_3D *m2, double xr1, double xr2) { if (xr1 < 0.0 || xr1 > 1.0) { fprintf(stderr, "Corrleation: radius = %f\n", xr1 * RADIUS); exit(1); } if (xr2 < 0.0 || xr2 > 1.0) { fprintf(stderr, "Corrleation: radius = %f\n", xr2 * RADIUS); exit(1); } m1->SumRad(xr1); m2->SumRad(xr2); int Lmax = m1->GetLmax(xr1); int Lmax2 = m2->GetLmax(xr2); if (Lmax2 < Lmax) // Only comparison between common degrees Lmax = Lmax2; double *corr = new double[Lmax+2]; double totalpower1 = 0.0; double totalpower2 = 0.0; double totalpower12 = 0.0; double a1, b1, a2, b2; for (int L = 0; L <= Lmax; L++) { double power1 = 0.0; double power2 = 0.0; double power12 = 0.0; for (int m = 0; m <= L; m++) { a1 = m1->Alm[L][m]; b1 = m1->Blm[L][m]; a2 = m2->Alm[L][m]; b2 = m2->Blm[L][m]; power1 += a1 * a1 + b1 * b1; power2 += a2 * a2 + b2 * b2; power12 += a1 * a2 + b1 * b2; totalpower1 += a1 * a1 + b1 * b1; totalpower2 += a2 * a2 + b2 * b2; totalpower12 += a1 * a2 + b1 * b2; } corr[L] = (power12 / sqrt(power1 * power2)); } corr[Lmax+1] = totalpower12 / sqrt(totalpower1 * totalpower2); m1->FreeAB(xr1); m2->FreeAB(xr2); return corr; } // Implementation of classes BSpline, RadialBSpline, RadialBSplineModel // RadialBSplineSplitModel const double BSpline::sixth = 0.16666666666666667; // 1/6 BSpline::BSpline (double h0) : h(h0) {} // Calculate f(delta) double BSpline::fdelta(double delta) const { double deltil, value; deltil = delta / h; if (deltil <= 1.0) { value = (4.0 + 3.0 * deltil * deltil * (deltil - 2.0)); } else if (deltil > 2.0) { return 0.0; } else { value = 2.0 - deltil; value *= value * value; } return (value * sixth); } /* combines fdelta and dfddelta together. * Hard to read. * Efficiency is NO. 1 concern here. */ double BSpline::f_df_delta(double delta, double *df) const { double deltil, value; double value_sq, deltil_sq; deltil = delta / h; if (deltil <= 1.0) { deltil_sq = deltil * deltil; value = 4.0 + 3.0 * deltil_sq * (deltil - 2.0); *df = (1.5 * deltil_sq - deltil - deltil) / h; } else if (deltil > 2.0) { *df = 0.0; return 0.0; } else { value = 2.0 - deltil; value_sq = value * value; value *= value_sq; *df = -0.5 * value_sq / h; } return (value * sixth); } // Calculate derivative df/ddelta double BSpline::dfddelta(double delta) const { double deltil, value; deltil = delta / h; if (deltil <= 1.0) { value = deltil * (-2.0 + 1.5 * deltil) / h; } else if (deltil > 2.0) { return 0.0; } else { value = 2.0 - deltil; value *= -0.5 * value / h; } return value; } /* combines fdelta, dfddelta and dfdddelta together. * Hard to read. * Efficiency is NO. 1 concern here. */ double BSpline::f_df_ddf_delta(double delta, double *df, double *ddf) const { double deltil, value; double value_sq, deltil_sq; deltil = delta / h; if (deltil <= 1.0) { deltil_sq = deltil * deltil; value = 4.0 + 3.0 * deltil_sq * (deltil - 2.0); *df = (1.5 * deltil_sq - deltil - deltil) / h; *ddf = (3.0 * deltil - 2.0) / (h * h); } else if (deltil > 2.0) { *df = 0.0; *ddf = 0.0; return 0.0; } else { value = 2.0 - deltil; *ddf = value / (h * h); value_sq = value * value; value *= value_sq; *df = -0.5 * value_sq / h; } return (value * sixth); } // Calculate second derivative ddf/dddelta double BSpline::dfdddelta(double delta) const { double deltil, value; deltil = delta / h; if (deltil <= 1.0) { value = (3.0 * deltil - 2.0) / (h * h); } else if (deltil > 2.0) { return 0.0; } else { value = (2.0 - deltil) / (h * h); } return value; } /* Function BSplKer-- * to calculate cubic B-spline. * written by Wei-jia Su (in wslib.c) * modified by Xian-feng Liu * Input: x (given depth or radius); * Output: B[0],B[1], B[2], B[3] spline coefficients. * Also return the knot value (beginning from 1). * Warning: This function uses static values which * is bad for two different kinds of spline coefficients. * */ int RadialBSpline::BSplKer(double *B, double x) const { double *xn, temp, temp1, X1, X2, X3, X4; int i, ib, ie; static double A1, A2, A3, A4; static int i0 = 0; if (((x < xnode[1] || x > xnode[N]) && increase) || ((x < xnode[N] || x > xnode[1]) && !increase)) { B[0] = B[1] = B[2] = B[3] = 0.0; return 0; } ib = 1; ie = N; while (ie - ib > 1) { i = ib + ((ie - ib) >> 1); if ((x < xnode[i] && increase) || (x > xnode[i] && !increase)) ie = i; else ib = i; } i = ib; xn = &xnode[i-2]; if (i != i0) { temp = 1.0 / (xn[3] - xn[2]); A1 = temp / (xn[3] - xn[1]); A3 = temp / (xn[4] - xn[2]); A4 = A3 / (xn[5] - xn[2]); temp = xn[4] - xn[1]; A2 = A1 / temp; A3 /= temp; A1 /= xn[3] - xn[0]; i0 = i; } X1 = x - xn[1]; X2 = x - xn[2]; X3 = x - xn[3]; X4 = x - xn[4]; temp = X3 * X3 * A1; temp1 = X3 * X1 * A2 + X2 * X4 * A3; B[0] = -temp * X3; B[1] = (x - xn[0]) * temp + X4 * temp1; temp = X2 * X2 * A4; B[3] = X2 * temp; B[2] = -X1 * temp1 - (x - xn[5]) * temp; /* * Imposing boundary conditions: * for 2nd derivative = 0 */ if (i == 1) { B[2] += B[1] / 3.0; B[1] = B[1] / 1.5 + B[0]; B[0] = 0.0; } else if (i == 2) { B[1] += B[0] / 3.0; B[0] /= 1.5; } else if (i == N - 1) { B[1] += B[2] / 3.0; B[2] = B[2] / 1.5 + B[3]; B[3] = 0.0; } else if (i == N - 2) { B[2] += B[3] / 3.0; B[3] /= 1.5; } return i; } /* Function BSplKer1-- * for testing purpose: current does not work * to calculate cubic B-spline using recursive relation. * Input: x (given depth or radius); * Output: B[0],B[1], B[2], B[3] spline coefficients. * Also return the knot value (beginning from 1). * Warning: This function uses static values which * is bad for two different kinds of spline coefficients. * */ int RadialBSpline::BSplKer1(double **coeff, double x) const { double *xn, temp, temp1, X1, X2, X3, X4; int i, ib, ie; static double A1, A2, A3, A4; static int i0 = 0; if (((x < xnode[1] || x > xnode[N]) && increase) || ((x < xnode[N] || x > xnode[1]) && !increase)) { (*coeff)[0] = (*coeff)[1] = (*coeff)[2] = (*coeff)[3] = 0.0; return 0; } ib = 1; ie = N; while (ie - ib > 1) { i = ib + ((ie - ib) >> 1); if ((x < xnode[i] && increase) || (x > xnode[i] && !increase)) ie = i; else ib = i; } i = ib; int loc = i; int order = 3; int node = N; double **B = matrix(1, order+1, -1, node+2); for (i = -1; i <= node+1; i++) { if ((increase && (x < xnode[i+1] && x >= xnode[i])) || (!increase && (x >= xnode[i+1] && x < xnode[i]))){ B[1][i] = 1.0; } else B[1][i] = 0.0; } double ci, cip1; for (int k = 2; k <= order+1; k++) { for (i = 0; i <= node; i++) { if ((i+k) <= node+2 && xnode[i+k] != xnode[i]) ci = (x - xnode[i]) / (xnode[i+k] - xnode[i]); else ci = 0.0; if ((i+k+1) <= node+2 && xnode[i+k+1] != xnode[i+1]) cip1 = (xnode[i+k+1] - x) / (xnode[i+k+1] - xnode[i+1]); else cip1 = 0.0; B[k][i] = ci * B[k-1][i] + cip1 * B[k-1][i+1]; } B[k][node+1] = 0.0; B[k][-1] = B[k][0] = 0.0; } *coeff = &(B[order+1][loc-3]); return loc; } /* * to calculate cubic B-spline Kernel and its derivative */ int RadialBSpline::DBBSplKer(double *B, double *DB, double x) const { static double A1, A2, A3, A4; double *xn, temp, X0, X1, X2, X3, X4, X5; double XX1, XX2, XX3, XX4; static int i0 = 0; int i, ib, ie; if (((x < xnode[1] || x > xnode[N]) && increase) || ((x < xnode[N] || x > xnode[1]) && !increase)) { B[0] = B[1] = B[2] = B[3] = 0.0; DB[0] = DB[1] = DB[2] = DB[3] = 0.0; return 0; } ib = 1; ie = N; while (ie - ib > 1) { i = ib + ((ie - ib) >> 1); if ((x < xnode[i] && increase) || (x > xnode[i] && !increase)) ie = i; else ib = i; } i = ib; xn = &xnode[i-2]; if (i != i0) { temp = 1.0 / (xn[3] - xn[2]); A1 = temp / (xn[3] - xn[1]); A3 = temp / (xn[4] - xn[2]); A4 = A3 / (xn[5] - xn[2]); temp = xn[4] - xn[1]; A2 = A1 / temp; A3 /= temp; A1 /= xn[3] - xn[0]; i0 = i; } X0 = x - xn[0]; X1 = x - xn[1]; X2 = x - xn[2]; X3 = x - xn[3]; X4 = x - xn[4]; X5 = x - xn[5]; XX1 = X1 * X1; XX2 = X2 * X2; XX3 = X3 * X3; XX4 = X4 * X4; temp = XX3 * A1; B[0] = -temp * X3; B[1] = X0 * temp; DB[0] = -3.0 * temp; DB[1] = temp + 2.0 * X0 * X3 * A1; temp = XX2 * A4; B[3] = X2 * temp; DB[3] = 3.0 * temp; B[2] = -X5 * temp; DB[2] = -temp - 2.0 * X5 * X2 * A4; temp = XX1 * A2; B[2] -= X3 * temp; DB[2] -= temp; temp = XX4 * A3; B[1] += X2 * temp; DB[1] += temp + 2.0 * X2 * X4 * A3; temp = X1 * X3 * A2; B[1] += X4 * temp; DB[1] += temp + X4 * (X1 + X3) * A2; temp = X1 * X4 * A3; B[2] -= X2 * temp; DB[2] -= temp + X2 * (X1 + X4) * A3; /* Imposing boundary conditions: * for 2nd derivative = 0 */ if (i == 1) { B[2] += B[1] / 3.0; B[1] = B[1] / 1.5 + B[0]; DB[2] += DB[1] / 3.0; DB[1] = DB[1] / 1.5 + DB[0]; B[0] = DB[0] = 0.0; } else if (i == 2) { B[1] += B[0] / 3.0; B[0] /= 1.5; DB[1] += DB[0] / 3.0; DB[0] /= 1.5; } else if ( i == N - 1) { B[1] += B[2] / 3.0; B[2] = B[2] / 1.5 + B[3]; DB[1] += DB[2] / 3.0; DB[2] = DB[2] / 1.5 + DB[3]; B[3] = DB[3] = 0.0; } else if (i == N - 2) { B[2] += B[3] / 3.0; B[3] /= 1.5; DB[2] += DB[3] / 3.0; DB[3] /= 1.5; } return i; } /* * to calculate cubic B-spline Kernel and its derivative */ int RadialBSpline::DDBBSplKer(double *B, double *DB, double *DDB, double x) const { static double A1, A2, A3, A4; double *xn, temp, X0, X1, X2, X3, X4, X5; double XX1, XX2, XX3, XX4; static int i0 = 0; int i, ib, ie; if (((x < xnode[1] || x > xnode[N]) && increase) || ((x < xnode[N] || x > xnode[1]) && !increase)) { B[0] = B[1] = B[2] = B[3] = 0.0; DB[0] = DB[1] = DB[2] = DB[3] = 0.0; DDB[0] = DDB[1] = DDB[2] = DDB[3] = 0.0; return 0; } ib = 1; ie = N; while (ie - ib > 1) { i = ib + ((ie - ib) >> 1); if ((x < xnode[i] && increase) || (x > xnode[i] && !increase)) ie = i; else ib = i; } i = ib; xn = &xnode[i-2]; if (i != i0) { temp = 1.0 / (xn[3] - xn[2]); A1 = temp / (xn[3] - xn[1]); A3 = temp / (xn[4] - xn[2]); A4 = A3 / (xn[5] - xn[2]); temp = xn[4] - xn[1]; A2 = A1 / temp; A3 /= temp; A1 /= xn[3] - xn[0]; i0 = i; } X0 = x - xn[0]; X1 = x - xn[1]; X2 = x - xn[2]; X3 = x - xn[3]; X4 = x - xn[4]; X5 = x - xn[5]; XX1 = X1 * X1; XX2 = X2 * X2; XX3 = X3 * X3; XX4 = X4 * X4; temp = XX3 * A1; B[0] = -temp * X3; B[1] = X0 * temp; DB[0] = -3.0 * temp; DB[1] = temp; temp = XX2 * A4; B[3] = X2 * temp; DB[3] = 3.0 * temp; B[2] = -X5 * temp; DB[2] = -temp; temp = XX1 * A2; B[2] -= X3 * temp; DB[2] -= temp; temp = XX4 * A3; B[1] += X2 * temp; DB[1] += temp; temp = X1 * X3 * A2; B[1] += X4 * temp; DB[1] += temp; DB[2] -= 2.0 * temp; temp = X1 * X4; DB[1] += temp * A2; temp *= A3; B[2] -= X2 * temp; DB[2] -= temp; temp = X2 * X4 * A3; DB[1] += 2.0 * (temp + X0 * X3 * A1) + X3 * X4 * A2; DB[2] -= temp + X2 * (X1 * A3 + 2.0 * X5 * A4); /* Second derivative */ DDB[0] = -6.0 * X3 * A1; DDB[1] = 2.0 * (X3 + X3 + X0) * A1 + 2.0 * (X1 + X3 + X4) * A2 + 2.0 * (X2 + X4 + X4) * A3; DDB[2] = -2.0 * (X1 + X1 + X3) * A2 - 2.0 * (X1 + X2 + X4) * A3 -2.0 * (X2 + X2 + X5) * A4; DDB[3] = 6.0 * X2 * A4; /* Imposing boundary conditions: * for 2nd derivative = 0 */ if (i == 1) { B[2] += B[1] / 3.0; B[1] = B[1] / 1.5 + B[0]; DB[2] += DB[1] / 3.0; DB[1] = DB[1] / 1.5 + DB[0]; DDB[2] += DDB[1] / 3.0; DDB[1] = DDB[1] / 1.5 + DDB[0]; B[0] = DB[0] = DDB[0] = 0.0; } else if (i == 2) { B[1] += B[0] / 3.0; B[0] /= 1.5; DB[1] += DB[0] / 3.0; DB[0] /= 1.5; DDB[1] += DDB[0] / 3.0; DDB[0] /= 1.5; } else if (i == N - 1) { B[1] += B[2] / 3.0; B[2] = B[2] / 1.5 + B[3]; DB[1] += DB[2] / 3.0; DB[2] = DB[2] / 1.5 + DB[3]; DDB[1] += DDB[2] / 3.0; DDB[2] = DDB[2] / 1.5 + DDB[3]; B[3] = DB[3] = DDB[3] = 0.0; } else if (i == N - 2) { B[2] += B[3] / 3.0; B[3] /= 1.5; DB[2] += DB[3] / 3.0; DB[3] /= 1.5; DDB[2] += DDB[3] / 3.0; DDB[3] /= 1.5; } return i; } // radialKnot: an array of radial knot values; // numKnots: number of knots in radialKnot RadialBSplineModel::RadialBSplineModel (double radialKnot[], int numKnots) { n_rad = numKnots; co_rad = radialKnot; radialBSpline = new RadialBSpline(radialKnot, n_rad); } // radialKnot: an array of radial knot values; // numKnots: number of knots in radialKnot RadialBSplineSplitModel::RadialBSplineSplitModel (double radialKnotUM[], int numKnotsUM, double radialKnotLM[], int numKnotsLM) { n_radUM = numKnotsUM; co_radUM = radialKnotUM; radialBSplineUM = new RadialBSpline(radialKnotUM, n_radUM); n_radLM = numKnotsLM; co_radLM = radialKnotLM; radialBSplineLM = new RadialBSpline(radialKnotLM, n_radLM); } /* convert a latitude into a colatitude */ double lat2co(double latitude) { double theta; theta = HALF_PI - atan(GEOCO * tan(latitude * RADIAN)); return (theta); } /* convert a geocentral latitude into a spherical coordinate theta */ double Geocen2theta(double geocen_latitude) { double theta; theta = HALF_PI - geocen_latitude * RADIAN; return (theta); } /* convert a longitude into a spherical phi */ double long2phi(double longitude) { double phi; phi = longitude * RADIAN; if (phi < 0.0) phi += TWO_PI; return (phi); } /* Given the source colatitude scol, the source longitude slon, the receiver colatitude rcol, the receiver longitude rlon, the subroutine euler calculates the angular distance delta between the source and the receiver, in addition to the Euler angles alpha, beta and gamma. */ void euler(double scol, double slon, double rcol, double rlon, double *delta, double *alpha, double *beta, double *gamma) { double az, sina, cosa, eta, mu, ep; *delta = acos(cos(scol) * cos(rcol) + sin(scol) * sin(rcol) * cos(rlon - slon)); sina = sin(rlon - slon) * sin(rcol) / sin(*delta); cosa = (cos(rcol) - cos(scol) * cos(*delta)) / (sin(scol) * sin(*delta)); az = atan2(sina, cosa); sina = sin(az); cosa = cos(az); if (fabs(cosa) < 1.0e-5) cosa = 1.0e-5; eta = atan(cos(scol) / (cosa * sin(scol))); if (fabs(eta) < 1.0e-5) eta = 1.0e-5; else if (eta < 0.0) eta += PI; *gamma = HALF_PI + eta; mu = atan2(sina * sin(scol), cos(scol) / sin(eta)); *beta = HALF_PI - mu; ep = atan2(sina * sin(eta), cos(eta) / sin(scol)); *alpha = - HALF_PI + slon - ep; } /* The subroutine rot rotates coordinate frames based on the Euler angles alpha, beta and gamma (see Edmonds page 7 and page 53). The input theta and phi are the coordinates in the geographic coordinate system. The angle zeta is measured counterclockwise from due east. The returned values of theta, phi, and zeta are in the equatorial coordinate system. The returned value of theta is between 0 and pi, and the returned value of phi is between 0 and 2*pi. The inverse operation, i.e. from the rotated (equatorial) frame back to the original (geographic) frame is performed by replacing alpha by -gamma, beta by -beta, and gamma by -alpha in the call to the subroutine rot. */ void rot(double alpha, double beta, double gamma, double *theta, double *phi, double *zeta) { int i, j; double r_hat[3], rp_hat[3], k_hat[3], kp_hat[3], rotation[3][3]; double sint, cost, sinp, cosp, sina, cosa, sinb; double cosb, sing, cosg, sinz, cosz; sint = sin(*theta); cost = cos(*theta); sinp = sin(*phi); cosp = cos(*phi); sina = sin(alpha); cosa = cos(alpha); sinb = sin(beta); cosb = cos(beta); sing = sin(gamma); cosg = cos(gamma); sinz = sin(*zeta); cosz = cos(*zeta); r_hat[0] = sint * cosp; r_hat[1] = sint * sinp; r_hat[2] = cost; k_hat[0] = -(sinz * cost * cosp + cosz * sinp); k_hat[1] = cosz * cosp - sinz * cost * sinp; k_hat[2] = sinz * sint; rotation[0][0] = cosg * cosb * cosa - sing * sina; rotation[0][1] = cosg * cosb * sina + sing * cosa; rotation[0][2] = -cosg * sinb; rotation[1][0] = -sing * cosb * cosa - cosg * sina; rotation[1][1] = -sing * cosb * sina + cosg * cosa; rotation[1][2] = sing * sinb; rotation[2][0] = sinb * cosa; rotation[2][1] = sinb * sina; rotation[2][2] = cosb; for (i = 0; i < 3; i++) { rp_hat[i] = 0.0; kp_hat[i] = 0.0; for (j = 0; j < 3; j++) { rp_hat[i] += rotation[i][j] * r_hat[j]; kp_hat[i] += rotation[i][j] * k_hat[j]; } } *theta = atan2(sqrt(rp_hat[0] * rp_hat[0] + rp_hat[1] * rp_hat[1]), rp_hat[2]); *phi = atan2(rp_hat[1], rp_hat[0]); reduce(theta, phi); sint = sin(*theta); cost = cos(*theta); sinp = sin(*phi); cosp = cos(*phi); sinz = -(kp_hat[0] * cost * cosp + kp_hat[1] * cost * sinp - kp_hat[2] * sint); cosz = -kp_hat[0] * sinp + kp_hat[1] * cosp; *zeta = atan2(sinz, cosz); } /* The subroutine rot rotates coordinate frames based on the Euler angles alpha, beta and gamma (see Edmonds page 7 and page 53). The input theta and phi are the coordinates in the geographic coordinate system. The angle zeta is measured counterclockwise from due east. The returned values of theta, phi, and zeta are in the equatorial coordinate system. The returned value of theta is between 0 and pi, and the returned value of phi is between 0 and 2*pi. The inverse operation, i.e. from the rotated (equatorial) frame back to the original (geographic) frame is performed by replacing alpha by -gamma, beta by -beta, and gamma by -alpha in the call to the subroutine rot. */ void rot(double alpha, double beta, double gamma, double *theta, double *phi) { int i, j; double r_hat[3], rp_hat[3]; double sint, cost, sinp, cosp; static Boolean initialize = FALSE; static double rotation[3][3]; static double alp, bet, gam; sint = sin(*theta); cost = cos(*theta); sinp = sin(*phi); cosp = cos(*phi); if (!initialize || alpha != alp || beta != bet || gamma != gam) { double sina = sin(alpha); double cosa = cos(alpha); double sinb = sin(beta); double cosb = cos(beta); double sing = sin(gamma); double cosg = cos(gamma); rotation[0][0] = cosg * cosb * cosa - sing * sina; rotation[0][1] = cosg * cosb * sina + sing * cosa; rotation[0][2] = -cosg * sinb; rotation[1][0] = -sing * cosb * cosa - cosg * sina; rotation[1][1] = -sing * cosb * sina + cosg * cosa; rotation[1][2] = sing * sinb; rotation[2][0] = sinb * cosa; rotation[2][1] = sinb * sina; rotation[2][2] = cosb; alp = alpha; bet = beta; gam = gamma; initialize = TRUE; } r_hat[0] = sint * cosp; r_hat[1] = sint * sinp; r_hat[2] = cost; double rp_hat_i, *rotation_i; for (i = 0; i < 3; i++) { rp_hat_i = 0.0; rotation_i = rotation[i]; for (j = 0; j < 3; j++) { rp_hat_i += rotation_i[j] * r_hat[j]; } rp_hat[i] = rp_hat_i; } double rp_hat0 = rp_hat[0]; double rp_hat1 = rp_hat[1]; *theta = atan2(sqrt(rp_hat0 * rp_hat0 + rp_hat1 * rp_hat1), rp_hat[2]); *phi = atan2(rp_hat1, rp_hat0); reduce(theta, phi); } void reduce(double *pth, double *pph) { /* force theta between 0 and pi, and phi between 0 and 2*pi */ int i; double th, ph; th = *pth; ph = *pph; i = abs((int)(ph / TWO_PI)); if (ph < 0.0) ph += (i + 1) * TWO_PI; else if (ph > TWO_PI) ph -= i * TWO_PI; *pph = ph; if (th < 0.0 || th > PI) { i = (int)(th / PI); if (th > 0.0) { if (ODD(i)) { th = (i + 1) * PI - th; ph = (ph < PI) ? ph + PI : ph - PI; } else th -= i * PI; } else { if (EVEN(i)) { th = -th + i * PI; ph = (ph < PI) ? ph + PI : ph - PI; } else th -= i * PI; } *pth = th; *pph = ph; } } /* * Given a positive-definite symmetric matrix a[1..n][1..n], this routine * constructs its Cholesky decomposition, A=L*LT. On input, only * the upper triangle of a need be given; it is not modified. The * Cholesky factor L is returned in the lower triangle of a, * except for its diagonal elements which are returned in p[1..n]. */ void choldc(double **a, int n, double p[]) { void nrerror(char error_text[]); int i, j, k; double sum; for (i = 1; i <= n; i++) { for (j = i; j <= n; j++) { for (sum = a[i][j], k = i - 1; k >= 1; k--) sum -= a[i][k] * a[j][k]; if (i == j) { if (sum <= 0.0) nrerror("choldc failed"); p[i] = sqrt(sum); } else a[j][i] = sum / p[i]; } } } /* * Solves the set of n linear equations A*x=b, where a is a positive-definte * symmetric matrix. a[1..n][1..n] and p[1..n] are input as the output of * the routine choldc. Only the lower triangle of a is accessed. b[1..n] is * input as the right-hand side vector. The solution vector is returned * in x[1..n]. a, n and p are not modified and can be left in place for * successive calls with different right-hand sides b. b is not modifed * unless you identify b and x in the calling sequence, which is allowed. */ void cholsl(double **a, int n, double p[], double b[], double x[]) { int i, k; double sum; for (i = 1; i <= n; i++) { for (sum = b[i], k = i - 1; k >= 1; k--) sum -= a[i][k] * x[k]; x[i] = sum / p[i]; } for (i = n; i >= 1; i--) { for (sum = x[i], k = i + 1; k <= n; k++) sum -= a[k][i] * x[k]; x[i] = sum / p[i]; } } /* (C) Copr. 1986-92 Numerical Recipes Software Q26>$. */ /* function binary_search-- * find x in v[0] <= v[1] <= ... <= v[n-1] * or in v[0] >= v[1] >= ... >= v[n-1] */ int binary_search(double x, double v[], int n) { int low, high, mid; Boolean increase; if (v[n-1] > v[0]) increase = TRUE; else increase = FALSE; low = 0; high = n - 1; if (increase) { //find x in v[0] <= v[1] <= ... <= v[n-1] while (low <= high) { mid = (low + high) / 2; if (x < v[mid]) high = mid - 1; else if ( x > v[mid+1]) low = mid + 1; else return mid; } } else { //find x in v[0] >= v[1] >= ... >= v[n-1] while (low <= high) { mid = (low + high) / 2; if (x > v[mid]) high = mid - 1; else if ( x < v[mid+1]) low = mid + 1; else return mid; } } return (-1); /* no match */ } // Given angular order L, calculate Xln(PI/2) from n= -L to L // in equation 15 of Dziewonski (1984, JGR, Page 5934) // Xln is a vector from -L to L // Convention used by Robert Woodward // Translated from Fortran which was written by Bob Woodward // It is equivalent to a call Ylm(HALF_PI, L, y), the output // y[L][] is half of Xln[]. void XlnHpi(double Xln[], int L) { double L2p1 = L + L + 1; double p = 1.0; int i, m; for (i = 1; i <= L; i++) p *= (i + i - 1.0) / (i + i); Xln[L] = sqrt(L2p1 * p / FOUR_PI) * (1 - 2 * ( L % 2)); double x1 = 0.0; for (i = 1; i <= L; i++) { m = L - i; p = sqrt(i * (L2p1 - i)); Xln[m] = -x1 / p; x1 = Xln[m+1] * p; } double sgn = -1.0; for (m = 1; m <= L; m++) { Xln[-m] = sgn * Xln[m]; sgn = -sgn; } return; } /* function rotylm-- * This subroutine calculates bmp = Dmpm as defined by * equation (4.1.12) of Edmonds (1960). * alpha, beta and gamma are the Euler angles (radians). * s = angular degree l, t = azimuthal order m, and mp = m'. * The relation between Y_l^m in the geographical (i.e. not * rotated) reference frame and Y'_l^m' in the equatorial (i.e. * rotated) reference frame is * Y_l^m (geographic) = Sum_m' D_m'm Y'_l^m' (equatorial) * After the rotation, the source-receiver path is along the * equator (theta=pi/2) with the source at phi=0 and receiver * at phi=DELTA. The real and imaginary parts of Dmpm are * returned in bmp(2). */ static double pjac0, pjac1; void rotylm(double alpha, double beta, double gamma, int s, int mp, int t, double *bmp) { int facsym, tsave, i1, i2, mpsave, mpp, iferror; int i3, i4, i5, i6, n0, n1, i1s, i2s; double cb2, sb2; double cosb, a1, a2, anorm, xs, pjac; double hbeta, dmag1, dmag2, dmpm, arg; i1 = mp - t; i2 = mp + t; /* For negative values of i1 or i2, must use symmetry properties of Edmonds, 1960, eqn's (4.2.5) and (4.2.6).*/ facsym = 1; mpsave = mp; tsave = t; iferror = 0; if (i1 >= 0 && i2 < 0) { mpp = mpsave; mpsave = -tsave; tsave = -mpp; i1 = mpsave - tsave; i2 = mpsave + tsave; facsym = 1; } else if (i1 < 0 && i2 >= 0) { mpp = mpsave; mpsave = tsave; tsave = mpp; i1 = mpsave - tsave; i2 = mpsave + tsave; if ((i1 - 2 * (i1 / 2)) == 0) facsym = 1; else facsym = -1; } else if (i1 < 0 && i2 < 0) { mpsave = - mpsave; tsave = -tsave; i1 = mpsave - tsave; i2 = mpsave + tsave; if ((i1 - 2 * (i1 / 2)) == 0) facsym = 1; else facsym = -1; } i3 = s - mpsave; i4 = s + mpsave; i5 = s + tsave; i6 = s - tsave; cosb = cos(beta); fact(i4, i5, &a1); fact(i3, i6, &a2); anorm = a1 + a2; /* First determine n=0 and n=1 values of pjac. Then recurrence relations will get actual value. */ n0 = 0; i1s = i1; i2s = i2; xs = cosb; pnab(n0, i1s, i2s, anorm, &xs, &pjac); if (i3 != 0) { n1 = 1; i1s = i1; i2s = i2; xs = cosb; pnab(n1, i1s, i2s, anorm, &xs, &pjac); } if (i3 != 1) { pnab(i3, i1, i2, anorm, &cosb, &pjac); /*Guard against underflow.*/ } hbeta = beta * 0.5; cb2 = cos(hbeta); sb2 = sin(hbeta); if (cb2 != 0.0 && sb2 != 0.0) { dmag1 = fabs((mpsave + tsave)) * log(fabs(cb2)); dmag2 = fabs((mpsave - tsave)) * log(fabs(sb2)); if (dmag1 < (-30.0) || dmag2 < (-30.0)) iferror = 1; } if (iferror != 1) { double pow_cb2 = 1.0, pow_sb2 = 1.0; int tmp_p = mpsave + tsave; int tmp_m = mpsave - tsave; if (tmp_p != 0) pow_cb2 = pow(cb2, tmp_p); if (tmp_m != 0) pow_sb2 = pow(sb2, tmp_m); dmpm = pow_cb2 * pow_sb2 * pjac; arg = (float) (mp) * gamma + (float) (t) * alpha; bmp[0] = cos(arg) * dmpm * (facsym); bmp[1] = sin(arg) * dmpm * (facsym); } else { bmp[0] = 0.0; bmp[1] = 0.0; } } /* Returns Jacobi polynomial P(n,[alpha,beta])(x) -- Edmonds, 1960, p. 58. * multiplied by [(s+mp)!(s-mp)!/(s+t)!(s-t)!]**0.5 */ void pnab(int n, int nalpha, int nbeta, double anorm, double *x, double *pjac) { double flog2, a1, a2, a3, a4, opx, omx, dmag1, dmag2, xlg; double apb, a2b2, fn0, fn1, a1n, a2n, a3n, a4n, fn2; double sum; int nbeta1, i; int ix, np1, nu, npa, npb, nmnu, npamnu, nupb, n1, n2; if (n <= 1) { flog2 = (n) * log(2.); ix = 1; if (*x < 0) { /* Use symmetry relation (4.1.20) */ if ((n - 2 * (n / 2)) == 0) ix = 1; else ix = -1; nbeta1 = nbeta; nbeta = nalpha; nalpha = nbeta1; *x = -*x; } sum = 0.0; np1 = n + 1; for (i = 1; i <= np1; i++) { nu = i - 1; npa = n + nalpha; npb = n + nbeta; nmnu = n - nu; npamnu = n + nalpha - nu; nupb = nu + nbeta; fact(npa, npamnu, &a1); n1 = 1; fact(n1, nu, &a2); fact(npb, nupb, &a3); fact(n1, nmnu, &a4); opx = 1. + *x; omx = 1. - *x; if (opx == 0.0) opx = 1.0E-8; if (omx == 0.0) omx = 1.0E-8; dmag1 = (nu) * (log(opx)); dmag2 = (n - nu) * (log(omx)); xlg = 0.5 * anorm + a1 + a2 + a3 + a4 + dmag1 + dmag2 - flog2; if (xlg >= (-20.0)) { if (EVEN(n - nu)) sum += exp(xlg); else sum += -exp(xlg); } } *pjac = sum * (ix); if (n == 0) pjac0 = *pjac; if (n == 1) pjac1 = *pjac; } else { /* Use recurrance relations. */ apb = (nalpha + nbeta); a2b2 = (nalpha * nalpha - nbeta * nbeta); fn0 = 0.0; fn1 = 1.0; for (n2 = 2; ; n2++) { fn2 = n2; a1n = 2. * fn2 * (fn1 + apb + 1.) * (2. * fn1 + apb); a2n = (2. * fn1 + apb + 1.) * (a2b2); a3n = (2. * fn1 + apb) * (2. * fn1 + apb + 1.) * (2. * fn1 + apb + 2.); a4n = 2. * (fn1 + (nalpha)) * (fn1 + (nbeta)) * (2. * fn1 + apb + 2.); *pjac = ((a2n + a3n * (*x)) / a1n) * (pjac1) - (a4n / a1n) * (pjac0); if (n2 == n) break; fn0 = fn1; fn1 = fn2; pjac0 = pjac1; pjac1 = *pjac; } } } /* Finds log[n1(factorial)/n2(factorial)], returned in a. */ void fact(int n1, int n2, double *a) { int m1, m2, m, i; double b; m1 = n1; m2 = n2; if (m1 == 0) m1 = 1; if (m2 == 0) m2 = 1; if (m2 <= m1) { *a = -log((m2)); m = m1 - m2 + 1; for (i = 1; i <= m; i++) (*a) += log((m2 + i - 1)); } else { b = log((m1)); m = m2 - m1 + 1; for (i = 1; i <= m; i++) b -= log((m1 + i - 1)); *a = b; } } // Exit, leave a message in file "dead.log" void exitLog(char *logMessage) { if (logMessage) { FILE *fp = fopen("dead.log", "w"); fprintf(fp, logMessage); fclose(fp); } exit(1); } /* Given arrays x[0..n] and y[0..n] containing a tabulated function * y_i = f(x_i), this routine returns an array of coefficients * cof[0..n] such that y_i = sum_j (cof_j * (x_i)^j) * From Numerical Recipes, page 120-121, second edition */ void polcoe(double x[], double y[], int n, double cof[]) { int k, j, i; double phi, ff, b, *s; s = new double[n+1]; for (i = 0; i <= n; i++) s[i] = cof[i] = 0.0; s[n] = -x[0]; for (i = 1; i <= n; i++) { for (j = n - i; j <= n - 1; j++) s[j] -= x[i] * s[j+1]; s[n] -= x[i]; } for (j = 0; j <= n; j++) { phi = n + 1; for (k = n; k >= 1; k--) phi = k * s[k] + x[j] * phi; ff = y[j] / phi; b = 1.0; for (k = n; k >= 0; k--) { cof[k] += b * ff; b = s[k] + x[j] * b; } } delete[] s; } MODEL3D_T GetModelType(void) { MODEL3D_T type; int tp; printf("Enter your model type->\n"); printf("8-Harvard Whole-Mantle earth model (e.g. S12WM13 or SKS12WM13);\n"); printf("9-Harvard Spherical Spline Model (starting at Moho): \n"); printf("10-Harvard Spline-Harmonics Mixed Whole Mantle Model (e.g. S16U6L8): \n"); printf("11-Harvard Spline-Harmonics Mixed Split Model (e.g. S16VB14): \n"); printf("12-Harvard Chebchev Split mantle model (e.g. S20U7L5):\n"); printf("13-Harvard Spherical Spline Model (starting at CMB, compatible with XGeomap format):\n"); printf("others-exit;\n"); scanf("%d", &tp); switch (tp) { case 8: type = CHEB_WHOLE_MODEL; break; case 9: type = SPH_SPLINE_MODEL; break; case 10: type = BSP_WHOLE_MODEL; break; case 11: type = BSP_SPLIT_MODEL; break; case 12: type = CHEB_SPLIT_MODEL; break; case 13: type = XV_SPH_SPLINE_MODEL; break; default: exit(1); } return type; } /* utility function from Numerical Recipes (Version 2) */ #define FREE_ARG char* static const int NR_END = 1; void nrerror(char error_text[]) /* Numerical Recipes standard error handler */ { fprintf(stderr,"run-time error...\n"); fprintf(stderr,"%s\n",error_text); fprintf(stderr,"...now exiting to system...\n"); exit(1); } void error(char error_text[]) { fprintf(stderr, "%s\n", error_text); exit(1); } void error(char error_text[], int error_value) { fprintf(stderr, "%s %d\n", error_text, error_value); exit(1); } void error(char error_text[], double error_value) { fprintf(stderr, "%s %f\n", error_text, error_value); exit(1); } double *vector(long nl, long nh) /* allocate a double vector with subscript range v[nl..nh] */ { double *v; v=(double *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(double))); if (!v) nrerror("allocation failure in vector()"); return v-nl+NR_END; } int *ivector(long nl, long nh) /* allocate an int vector with subscript range v[nl..nh] */ { int *v; v=(int *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(int))); if (!v) nrerror("allocation failure in ivector()"); return v-nl+NR_END; } unsigned char *cvector(long nl, long nh) /* allocate an unsigned char vector with subscript range v[nl..nh] */ { unsigned char *v; v=(unsigned char *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(unsigned char))); if (!v) nrerror("allocation failure in cvector()"); return v-nl+NR_END; } unsigned long *lvector(long nl, long nh) /* allocate an unsigned long vector with subscript range v[nl..nh] */ { unsigned long *v; v=(unsigned long *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(long))); if (!v) nrerror("allocation failure in lvector()"); return v-nl+NR_END; } float *fvector(long nl, long nh) /* allocate a double vector with subscript range v[nl..nh] */ { float *v; v=(float *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(float))); if (!v) nrerror("allocation failure in fvector()"); return v-nl+NR_END; } double **matrix(long nrl, long nrh, long ncl, long nch) /* allocate a double matrix with subscript range m[nrl..nrh][ncl..nch] */ { long i, nrow=nrh-nrl+1,ncol=nch-ncl+1; double **m; /* allocate pointers to rows */ m=(double **) malloc((size_t)((nrow+NR_END)*sizeof(double*))); if (!m) nrerror("allocation failure 1 in matrix()"); m += NR_END; m -= nrl; /* allocate rows and set pointers to them */ m[nrl]=(double *) malloc((size_t)((nrow*ncol+NR_END)*sizeof(double))); if (!m[nrl]) nrerror("allocation failure 2 in matrix()"); m[nrl] += NR_END; m[nrl] -= ncl; for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol; /* return pointer to array of pointers to rows */ return m; } float **fmatrix(long nrl, long nrh, long ncl, long nch) /* allocate a float matrix with subscript range m[nrl..nrh][ncl..nch] */ { long i, nrow=nrh-nrl+1,ncol=nch-ncl+1; float **m; /* allocate pointers to rows */ m=(float **) malloc((size_t)((nrow+NR_END)*sizeof(float*))); if (!m) nrerror("allocation failure 1 in matrix()"); m += NR_END; m -= nrl; /* allocate rows and set pointers to them */ m[nrl]=(float *) malloc((size_t)((nrow*ncol+NR_END)*sizeof(float))); if (!m[nrl]) nrerror("allocation failure 2 in matrix()"); m[nrl] += NR_END; m[nrl] -= ncl; for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol; /* return pointer to array of pointers to rows */ return m; } int **imatrix(long nrl, long nrh, long ncl, long nch) /* allocate a int matrix with subscript range m[nrl..nrh][ncl..nch] */ { long i, nrow=nrh-nrl+1,ncol=nch-ncl+1; int **m; /* allocate pointers to rows */ m=(int **) malloc((size_t)((nrow+NR_END)*sizeof(int*))); if (!m) nrerror("allocation failure 1 in matrix()"); m += NR_END; m -= nrl; /* allocate rows and set pointers to them */ m[nrl]=(int *) malloc((size_t)((nrow*ncol+NR_END)*sizeof(int))); if (!m[nrl]) nrerror("allocation failure 2 in matrix()"); m[nrl] += NR_END; m[nrl] -= ncl; for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol; /* return pointer to array of pointers to rows */ return m; } double **submatrix(double **a, long oldrl, long oldrh, long oldcl, long oldch, long newrl, long newcl) /* point a submatrix [newrl..][newcl..] to a[oldrl..oldrh][oldcl..oldch] */ { long i,j,nrow=oldrh-oldrl+1,ncol=oldcl-newcl; double **m; /* allocate array of pointers to rows */ m=(double **) malloc((size_t) ((nrow+NR_END)*sizeof(double*))); if (!m) nrerror("allocation failure in submatrix()"); m += NR_END; m -= newrl; /* set pointers to rows */ for(i=oldrl,j=newrl;i<=oldrh;i++,j++) m[j]=a[i]+ncol; /* return pointer to array of pointers to rows */ return m; } double **convert_matrix(double *a, long nrl, long nrh, long ncl, long nch) /* allocate a double matrix m[nrl..nrh][ncl..nch] that points to the matrix declared in the standard C manner as a[nrow][ncol], where nrow=nrh-nrl+1 and ncol=nch-ncl+1. The routine should be called with the address &a[0][0] as the first argument. */ { long i,j,nrow=nrh-nrl+1,ncol=nch-ncl+1; double **m; /* allocate pointers to rows */ m=(double **) malloc((size_t) ((nrow+NR_END)*sizeof(double*))); if (!m) nrerror("allocation failure in convert_matrix()"); m += NR_END; m -= nrl; /* set pointers to rows */ m[nrl]=a-ncl; for(i=1,j=nrl+1;i