/*
 * The program solves the ODE describing underwater sound ray propagation
 *
 *   z'' = - (dv/dz)/(a^2 v^3)
 *
 * where z is the vertical coordinate, x is the horizontal coordinate,
 * "'" denotes the derivative with respect to x, v(z) is the (depth-dependent)
 * speed of the sound waves, v(z) is a given function of z; A is a constant
 * determined by the initial conditions.
 * 
 *   a = ...
 *
 * This can be converted into a first order system 
 * by introducing a separate variable for the velocities, 
 * vx = x'(t) and vy = y'(t)
 *
 *   z'  = w
 *   w'  = - (dv(z)/dz)/(a^2 v(z)^3)
 *
 * The program evolves the solution for from (z, w) = (1, 0.01) at x=0 to x=100. 
 * The step-size h is automatically adjusted by the controller 
 * to maintain an absolute accuracy of 10^{-6} in the function 
 * values (z, w). 
 */

#include <stdio.h>
#include <math.h>

#include <gsl/gsl_errno.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_odeiv2.h>

int
func (double x, const double yy[], double f[], void *params)
{
  double asq = ((double *) params)[0];
  double delta = ((double *) params)[1];
  double z = yy[0], w = yy[1];
  double t, v;

  t = delta*(z - 1.);
  v = 1. + t*t/delta;

  f[0] = w;
  f[1] = -2*t/(asq*v*v*v);

  return GSL_SUCCESS;
}

int
main (void)
{
  size_t neqs = 2;
  double eps_abs = 1.e-6, eps_rel = 0.;

  gsl_odeiv2_step    *s = gsl_odeiv2_step_alloc (gsl_odeiv2_step_rkf45, neqs);
  gsl_odeiv2_control *c = gsl_odeiv2_control_y_new (eps_abs, eps_rel);
  gsl_odeiv2_evolve  *e = gsl_odeiv2_evolve_alloc (neqs);
    
  double params[2];
  gsl_odeiv2_system sys = {func, NULL, neqs, &params};
    
  double x = 0., x1 = 100.;
  double stepsize = 1e-6;

  double theta0 = -0.19; /* initial angle in radians */
  double z0 = 1.2;      /* initial depth in dimensionless units */
  double delta = 0.05;
  double v0 = 1. + delta*pow((z0-1.),2); /* sound speed at the initial depth */
  params[0] = pow(cos(theta0)/v0, 2);
  params[1] = delta;
  double yy[2] = { z0, tan(theta0) };
    
  while (x < x1) {
      int status = gsl_odeiv2_evolve_apply (e, c, s, &sys, &x, x1, &stepsize, yy);

      if (status != GSL_SUCCESS)
          break;

      printf ("% .5e  % .5e\n", x, yy[0]);
  }

  gsl_odeiv2_evolve_free (e);
  gsl_odeiv2_control_free (c);
  gsl_odeiv2_step_free (s);

  return 0;
}