/* Given a model type and its model name (produced by Harvard),
* this program calculate its velocity perturbation
* (include 1st, 2nd partial derivatives) for
* given location (latitude, longitude, and depth).
*
* Most people ask Harvard seismology group
* a program to calculate the velocity perturbation at
* given locations. Although this program meets this
* purpose, it "over-kills" everything.
* In fact, it can also calculate 1st, 2nd partial derivativs
* with respect to r, theta and phi which you need to calculate
* synthetic seismograms based on classic ray theory or
* Maslov theory. It is much more complicated to calculate
* 1st and 2nd order partial derivatives than velocity perturbations.
* It can also rotate the 3-D model to arbitrary Euler angles.
* This explains why the program is over 7,000 lines.
*
* If you only want to calculate velocity perturbations, you only
* need to call: dv2vonly(...).
* If you only want to calculate velocity perturbations and their
* 1st order derivatives, you only need to call: dv2v(...).
* If you only want to calculate velocity pertubations, their
* 1st order derivatives and their 2nd derivatives, you only
* need to call alldv2v(...).
*
* The 3-D model must be produced by Harvard seismology group.
* It supports only the following models:
* 8-Harvard S12WM13 or SKS12WM13 (file s12wm13.dat or sks12wm13.dat)
* 9-Harvard Spherical Spline Model (file s12wm13.XBS or
* sks12wm13.XBS, starting at Moho discontinuity)
* 10-Harvard Spline-Harmonics Mixed Whole Mantle Model (S16VB14)
* (s16vb14.dat)
* 11-Harvard Spline-Harmonics Mixed Split Model (S16U6L8)
* (s16u6l8.dat)
* 12-Harvard Chebchev Split mantle model (S20U7L5)
* (s20u7l5.dat).
* 13-Harvard Spherical Spline Model (reversed order in radial direction
* of type-9 model, starting at CMB):
* At this time (Feb., 1998), type-9 and type-13 models (they are equivalent)
* are not born yet. I expand the degree 12 models in terms of these
* kinds of models. I am sure that Harvard is going to publish this kind
* of model in the near future.
*
* You can distribute and modify this program as long as
* you keep the original comments made by Harvard seismology group.
*
* Very Technical: the program is written using C++ which is still
* the industry standard at this time.
* It has been tested by Sun C++ compiler 4.0
* and GNU C++ compiler 2.6.3. Due to the incompatibility of
* GNU C++ 2.6.x and 2.7.x, you may have trouble compiling
* the code in GNU C++ 2.7.x. The orginal author has no
* access to the 2.7.x version.
*
* You need 3 files: this file, model3Dlib.cc and model3Dlib.h.
* To compile this command, type "CC testmodel3D.cc model3Dlib.cc -lm"
* in a UNIX platform. If you use GNU C++ compiler, try
* "g++ testmodel3D.cc model3Dlib.cc -lm" instead.
*
* IMPORTANT. As an industry convention, I disclaim that:
* No Liability. IN NO EVENT WILL HARVARD SEISMOLOGY
* GROUP OR THE ORIGINAL AUTHOR BE LIABLE FOR ANY LOST REVENUE,
* PROFIT OR DATA, OR FOR SPECIAL, INDIRECT, CONSEQUENTIAL,
* INCIDENTAL OR PUNITIVE DAMAGES HOWEVER CAUSED AND REGARDLESS
* OF THEORY OF LIABILITY ARISING OUT OF THE USE OF OR INABILITY
* TO USE SOFTWARE, EVEN IF HARVARD SEISMOLOGY OR THE ORGINAL
* AUTHOR HAS BEEN ADVISED OF THE POSSIBLITY OF SUCH DAMAGES.
* --- X.-F. Liu, 1998.
* Any comments or feedback should be sent to Harvard seismology group,
* instead to the original author who has left Academia
* to make a living in the street. ^_^ I recommend that you should
* contact either Yu Gu (gu@seismology.harvard.edu) or
* Professor Dziewonski (dziewons@seismology.harvard.edu).
*/
#include
#include
#include
#include
#include "model3Dlib.h"
#ifndef EXIT_FAILURE // should be defined in ANSI C
const int EXIT_FAILURE = 1;
#endif
int main(void)
{
char modelname[MAXSTR];
MODEL_3D *Vmodel;
double locLat, locLong;
double scol, slon;
double dvdx, dvdth, dvdphi, d2vdthdph0;
double xr, dv;
double d2vd2x0, d2vd2th0, d2vdxdth0, d2vdxdph0;
double depth, dV;
double alpha, beta, gam; // Euler angles
int tp;
tp = GetModelType();
printf("Enter your model name: ");
scanf("%s", modelname);
// Euler angles: alpha, beta, gamma.
// If you want rotate the coordinates, you should
// should change these angles.
// In any case (no matter what alpha, beta, gamma are),
// it won't affect input or output at all.
alpha = beta = gam = 0.0;
// Make a generic 3-D model.
// Technical: serves a virtual constructor
Vmodel = MakeModel3D(modelname, (MODEL3D_T) tp);
if (Vmodel == NULL) {
fprintf(stderr, "Model %s does not exist ", modelname);
fprintf(stderr, " or you have no permission to read it\n");
return EXIT_FAILURE;
}
while (TRUE) {
// Asking given locations.
printf("Latitude, longitude, depth of your requested location: \n");
if ((scanf("%lf %lf %lf", &locLat, &locLong, &depth)) == EOF)
break;
scol = lat2co(locLat);
slon = long2phi(locLong);
xr = 1.0 - depth / RADIUS; // normalized radius
printf("locLat = %f locLong = %f\n", locLat, locLong);
// If you only want to calculate velocity perturbations, you only
// need to call: dv2vonly(...).
Vmodel->dv2vonly(xr, scol, slon, &dV);
printf("dV/V = %f at depth %d km for %s\n", dV,
(int) depth, modelname);
// If you only want to calculate velocity perturbations and their
// 1st order derivatives, you only need to call: dv2v(...).
Vmodel->dv2v(xr, scol, slon, &dv, &dvdx, &dvdth, &dvdphi);
// All these partial derivatives are derivatives of dV/V.
// e.g. dv/dx = d(dV/V) / dx; where x is the normalized radius.
// e.g. dv/dth = d(dV/V) / dTheta
// e.g. d2vd2x = dd(dV/V) / d(dx); second order derivatives
// etc.
printf("dV/V = %f, dv/dx = %f, dv/dth = %f dv/dphi = %f\n",
dv, dvdx, dvdth, dvdphi);
// If you only want to calculate velocity pertubations, their
// 1st order derivatives and their 2nd derivatives, you only
// need to call alldv2v(...).
Vmodel->alldv2v(xr, scol, slon, &dv, &dvdx, &dvdth, &dvdphi,
&d2vd2x0, &d2vd2th0, &d2vdxdth0, &d2vdxdph0,
&d2vdthdph0);
printf("dV/V = %f dv/dx = %f dv/dth = %f dv/dphi = %f\n",
dv, dvdx, dvdth, dvdphi);
printf("d2vd2x = %f, d2vd2th = %f, d2vdxdth = %f\n",
d2vd2x0, d2vd2th0, d2vdxdth0);
printf("d2vdxdph = %f, d2vdthdph = %f\n", d2vdxdph0, d2vdthdph0);
printf("\n");
}
}