/* bug_in_interp2d.c
 *
 * 29.6.2018 Brijesh Upadhaya 
 *
 * Test for bug 
 * Typo mistake at line 97 in file gsl-2.5/interpolation/gsl_interp2d.h)
 * 
 * Missing function "int gsl_interp2d_eval_extrap_e ()" 
 * The function "int gsl_interp2d_eval_extrap_e ()" is mentioned in gsl-ref.pdf
 *
 * However the actual function name is "int gsl_interp2d_eval_e_extrap ()"
 */

#include <stdio.h>
#include <stdlib.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_interp2d.h>


int 
main (void)
{
  const gsl_interp2d_type *T = gsl_interp2d_bilinear;
  const size_t N = 100;
  const double xa[] = { 0.0, 1.0 };
  const double ya[] = { 0.0, 1.0 };
  const size_t nx = sizeof (xa) / sizeof (double);
  const size_t ny = sizeof (ya) / sizeof (double);
    
  double *za = (double *) malloc (nx * ny * sizeof (double));
  
  gsl_interp_accel *xacc = gsl_interp_accel_alloc ();
  gsl_interp_accel *yacc = gsl_interp_accel_alloc ();

  gsl_interp2d *sp = gsl_interp2d_alloc (T, nx, ny);
  gsl_interp2d_set (sp, za, 0, 0, 0.0);
  gsl_interp2d_set (sp, za, 1, 0, 1.0);
  gsl_interp2d_set (sp, za, 0, 1, 1.0);
  gsl_interp2d_set (sp, za, 1, 1, 0.5);

  gsl_interp2d_init (sp, xa, ya, za, nx, ny);

  size_t i, j;
  for (i = 0; i < N; ++i)
    {
      double xi = (double) i / ((double) N - 1.0);
      for (j = 0; j < N; ++j)
        {
          double yi = (double) j / ((double) N - 1.0);
          double zi = gsl_interp2d_eval (sp, xa, ya, za, xi, yi, xacc, yacc);
          fprintf (stdout, "%f %f %f\n", xi, yi, zi);
        }
      fprintf (stdout, "\n");
    }

  double ze = gsl_interp2d_eval_extrap (sp, xa, ya, za, 1.05, 1.05, xacc, yacc);

  /* One of the following function is defined in gsl-ref.pdf */
  double ze0 = 0.0;
  int c = gsl_interp2d_eval_extrap_e (sp, xa, ya, za, 1.05, 1.05, xacc, yacc, &ze0);
  //int c = gsl_interp2d_eval_e_extrap (sp, xa, ya, za, 1.05, 1.05, xacc, yacc, &ze0);
  (void) c;

  fprintf (stderr, "%f %f\n", ze, ze0);

  gsl_interp2d_free (sp);
  gsl_interp_accel_free (yacc);
  gsl_interp_accel_free (xacc);
  free (za);
  return 0; 
}