#include <iostream>
#include <iomanip>
#include <cmath>

#include "JAcoustics/JSoundVelocity.hh"

#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"


int main(int argc, char **argv)
{
  using namespace std;
  using namespace JPP;

  double     depth;
  double     precision;
  int        debug;

  try {

    JParser<> zap;

    zap['D'] = make_field(depth, "depth relative to seabed [m] (< 0)")  = -2500.0;
    zap['e'] = make_field(precision)      = 1e-6;
    zap['d'] = make_field(debug)          =  3;

    zap(argc, argv);
  }
  catch(const exception &error) {
    FATAL(error.what() << endl);
  }


  JSoundVelocity V = getSoundVelocity[depth];

  for (double z = 0.0; z < 1000.0; z += 100.0) {
    DEBUG("Sound velocity at " << FIXED(5,1) << z << " " << FIXED(10,3) << V(z) << " [m/s] " << endl);
  }

  const double D = 500.0;
  const int    N = 1000;

  for (double ct = +1.0; ct >= 0.0; ct -= 0.1) {

    DEBUG("distance " << FIXED(5,2) << ct << ' ' << FIXED(5,0) << D << endl);

    const double z0 = 0.0;
    const double z1 = z0 + D * ct;
    const double ds = D / (double) N;

    double t1 = 0.0;

    for (double s = 0.0 + 0.5 * ds; s < D; s += ds) {

      const double z = z0 + ct * s;

      t1 += ds / V(z);
    }

    const double eps = 1.0e-6;

    const double Vi  = (V.getTime(D + 0.5*eps, z0, z1) - V.getTime(D - 0.5*eps, z0, z1)) / eps;

    DEBUG("Time [s]               " 
	  << FIXED(12,7) << t1                                   << ' ' 
	  << FIXED(12,7) << V.getTime(D,  z0, z1)                << endl);

    DEBUG("Distance [m]           " 
          << FIXED(12,5) << D                                    << ' '
          << FIXED(12,5) << V.getDistance(t1, z0, z1)            << endl);

    DEBUG("Inverse velocity [m/s] " 
	  << SCIENTIFIC(12,5) << Vi                              << ' ' 
	  << SCIENTIFIC(12,5) << V.getInverseVelocity(D, z0, z1) << endl);

    ASSERT(fabs(t1 - V.getTime           (D,  z0, z1)) < precision * t1, "Test of propagation time.");
    ASSERT(fabs(D  - V.getDistance       (t1, z0, z1)) < precision * D,  "Test of propagation distance.");
    ASSERT(fabs(Vi - V.getInverseVelocity(D,  z0, z1)) < precision * Vi, "Test of inverse velocity.");
  }

  return 0;
}