/* 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 <stdio.h>
#include <string.h>
#include <math.h>
#include <stdlib.h>
#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");
   }
}