//================================================================
// $Id: InterpolationExample.cpp,v 1.5 2005/10/29 06:29:48 harald Exp $
//================================================================

/// This file demonstrates how to use the interpolation routines
/// supplied with the initial data sets. 
///
/// COMPILATION
///
/// This executable must be linked against libSpEC_ID.a and against a
/// BLAS-library.  On my Linux-box, the command
///    g++ InterpolationExample.cpp -L. -lSpEC_ID -lblas 
/// produces an executable 'a.out'. 
///
/// libSpEC_ID.a was compiled with gcc 3.4.3 on a Linux system running
/// Red Hat Enterprise 4.  Thus, this library should link successfully
/// on most recent Linux workstations.  If you encounter problems or
/// have suggestions, please contact Harald Pfeiffer,
/// harald at tapir.caltech.edu

#include <iostream>
#include <vector>
#include <cmath>

#include "PublicID.hpp"


// see any C++ book for explanations of 'using' directive.
using std::vector;  
using std::cout;
using std::endl;


int main() {

  // Import data from the current directory on disk, store in RAM.
  // The 'double' parameter represents the orbital angular frequency for
  // the current dataset. (only applicable to BBH initial data)
  ReadData(0); 


  // assemble the coordinates of the points to interpolate to. 

  const unsigned int N=4; // four points for today

  // (see any C++ book for usage of vector).
  vector<double> x(N),y(N),z(N);
  x[0]=4;    y[0]=0;   z[0]= 0;
  x[1]=5;    y[1]=0.6; z[1]=-0.6;
  x[2]=1e4;  y[2]=3e4; z[2]= 4e4; 
  x[3]=0.5;  y[3]=1.6; z[3]= 0.7;

  // variables returning the interpolated data 
  vector<double> gxx, gxy, gxz, gyy, gyz, gzz;
  vector<double> Kxx, Kxy, Kxz, Kyy, Kyz, Kzz;
  vector<double> Shiftx, Shifty, Shiftz, Lapse;

  // Interpolate!
  InterpolateData(x,y,z,
		  gxx, gxy, gxz, gyy, gyz, gzz,
		  Kxx, Kxy, Kxz, Kyy, Kyz, Kzz,
		  Shiftx, Shifty, Shiftz, Lapse);
  

  for(unsigned int i=0; i<N; ++i) {
    // (see any C++ book for usage of cout and endl).
    cout << "At (" << x[i] << ", " << y[i] << ", " << z[i] << "):" << endl;
    cout << "gxx=" << gxx[i] << ", Kyz=" << Kyz[i]
	 << ", beta^y=" << Shifty[i] 
	 << ", N=" << Lapse[i] << endl << endl;
  };

  // now print some errors
  for(unsigned int i=0; i<x.size(); ++i) {
    const double X=x[i];
    const double Y=y[i];
    const double Z=z[i];
    const double R=sqrt(X*X+Y*Y+Z*Z);
    const double gxx_analytic=1 + 2*X*X/(R*R*R);
    const double gxy_analytic=    2*X*Y/(R*R*R);
    const double N_analytic = 1./sqrt(1+2/R);

    cout << "Difference to the analytical values for Kerr-Schild at (" 
	 << X << ", " << Y << ", " << Z << "):"  << endl;
    cout << "gxx-gxx_analytic = " << gxx[i] - gxx_analytic << endl;
    cout << "gxy-gxy_analytic = " << gxy[i] - gxy_analytic << endl;
    cout << "N  -  N_analytic = " << Lapse[i] - N_analytic << endl;
    cout << endl;
  }
};