#include <stdio.h>
#include <math.h>
#include <stdlib.h>		/* calloc() */
#include <string.h>		/* memset() */

#ifndef KRAKEN_DEF
#include "kraken-util.h"
#endif

int
main (int argc, char **argv) 
{
  char *kraken_modefile;
  double source_depth, source_range, receiver_start, receiver_delta, receiver_stop;
  double *z_orig, *z_new;
  complex *phis, *g, **phi_new;
  int arg_counter = 0, i;
  int nz_new;

  Modeheader *mhead;
  Mode *mode;

  /*
   * Parse command line arguments
   */
  if (argc == 1) {
    fprintf (stderr, "usage: %s kraken-mode-file source-depth source-range receiver-start receiver-delta recever-stop\n",
	     argv[0]);
    return(-1);
  }

  /* get KRAKEN mode file name */

  if (++arg_counter < argc) 
    kraken_modefile = argv[arg_counter];
  else {
    fprintf (stderr, "need KRAKEN mode file to calculate Green's function\n");
    return(-1);
  }

  if (++arg_counter < argc) {
    if (sscanf (argv[arg_counter], "%lf", &source_depth) != 1) {
      fprintf (stderr, "source depth should be a floating point number, not [%s]\n",
	       argv[arg_counter]);
      return(-1);
    }
  } else {
    fprintf (stderr, "need source depth to calculate Green's function\n");
    return(-1);
  }

  if (++arg_counter < argc) {
    if (sscanf (argv[arg_counter], "%lf", &source_range) != 1) {
      fprintf (stderr, "source range should be a floating point number, not [%s]\n",
	       argv[arg_counter]);
      return(-1);
    }
  } else {
    fprintf (stderr, "need source range to calculate Green's function\n");
    return(-1);
  }

  if (++arg_counter < argc) {
    if (sscanf (argv[arg_counter], "%lf", &receiver_start) != 1) {
      fprintf (stderr, "receiver start should be a floating point number, not [%s]\n",
	       argv[arg_counter]);
      return(-1);
    }
  } else {
    fprintf (stderr, "need receiver start to calculate Green's function\n");
    return(-1);
  }

  if (++arg_counter < argc) {
    if (sscanf (argv[arg_counter], "%lf", &receiver_delta) != 1) {
      fprintf (stderr, "receiver delta should be a floating point number, not [%s]\n",
	       argv[arg_counter]);
      return(-1);
    }
  } else {
    fprintf (stderr, "need receiver delta to calculate Green's function\n");
    return(-1);
  }

  if (++arg_counter < argc) {
    if (sscanf (argv[arg_counter], "%lf", &receiver_stop) != 1) {
      fprintf (stderr, "receiver stop should be a floating point number, not [%s]\n",
	       argv[arg_counter]);
      return(-1);
    }
  } else {
    fprintf (stderr, "need receiver stop to calculate Green's function\n");
    return(-1);
  }

  /*
   * End of command line argument parsing
   */

  /* 
   * Read in the KRAKEN header
   */
  mhead = read_kraken_mode_header (kraken_modefile, 0);
  if (mhead == NULL) {
    fprintf (stderr, "error reading in KRAKEN mode file header: [%s]\n",
	     kraken_modefile);
    return(-1);
  }

  /*
   * Read in the KRAKEN modeshapes
   */
  mode = read_all_kraken_modes (kraken_modefile, mhead->z, mhead->ntot,
				"N", 0);
  if (mode == NULL) {
    fprintf (stderr, "error reading in KRAKEN mode information from: [%s]\n",
	     kraken_modefile);
    return(-1);
  }

  /*
   * Get the modefunction evaluated at the source depth
   */
  z_new = (double *) calloc (1, sizeof(double));
  if (z_new == NULL) {
    perror ("calloc z_new");
    return(-1);
  }
  z_new[0] = source_depth;
  nz_new = 1;

  phi_new = (complex **) calloc (mhead->m, sizeof(complex *));
  if (phi_new == NULL) {
    perror ("calloc phi_new");
    return(-1);
  }
  for (i = 0; i < mhead->m; i++) {
    phi_new[i] = (complex *) calloc (nz_new, sizeof(complex));
    if (phi_new[i] == NULL) {
      perror ("calloc phi_new[i]");
      return(-1);
    }
  }

  z_orig = (double *) calloc (mhead->ntot, sizeof(double));
  if (z_orig == NULL) {
    perror ("caloc z_orig");
    return(-1);
  }
  for (i = 0; i < mhead->ntot; i++)
    z_orig[i] = mhead->z[i];

  if (spline_modeshapes (mhead->m,
			 mhead->ntot, z_orig, mode->phir,
			 nz_new, z_new, phi_new) < 0) {
    fprintf (stderr, "error splining source modeshapes\n");
    return(-1);
  }

  phis = (complex *) calloc (mhead->m, sizeof(complex));
  if (phis == NULL) {
    perror ("calloc phis");
    return(-1);
  }
  for (i = 0; i < mhead->m; i++) 
    phis[i] = phi_new[i][0];

  /*
   * Get the modefunctions evaluated at the desired receiver depths
   */
  nz_new = (int) ((receiver_stop - receiver_start) / receiver_delta + 1.0);
  z_new = (double *) realloc (z_new, nz_new * sizeof(double));
  if (z_new == NULL) {
    perror ("realloc z_new");
    return(-1);
  }
  for (i = 0; i < nz_new; i++)
    z_new[i] = receiver_start + i * receiver_delta;

  for (i = 0; i < mhead->m; i++) {
    phi_new[i] = (complex *) realloc(phi_new[i], nz_new * sizeof(complex));
    if (phi_new[i] == NULL) {
      perror ("realloc phi_new[i]");
      return(-1);
    }
  }

  if (spline_modeshapes (mhead->m,
			 mhead->ntot, z_orig, mode->phir,
			 nz_new, z_new, phi_new) < 0) {
    fprintf (stderr, "error splining receiver modeshapes\n");
    return(-1);
  }

  /*
   * Calculate Green's function at the selected range
   */
  g = (complex *) calloc (nz_new, sizeof(complex));
  if (g == NULL) {
    perror ("calloc g");
    return(-1);
  }

  (void)calc_greens_function (mhead->m, mhead->ck, nz_new, phi_new, phis,
			      source_range, g);

  /* Print the result out */

  for (i = 0; i < nz_new; i++)
    printf ("%f %e %e\n",
	    z_new[i],
	    complex_real(g[i]),
	    complex_imag(g[i]));

  /* All done! */

  return(0);
}

void
calc_greens_function (int nmodes, complex *k, int nz, complex **phi, complex *phis,
		      double range, complex *g)
     /* nmodes: number of modes
	k: complex mode numbers
	nz: number of receivers
	phi[i][j]: modeshape which goes with mode i evaluated at depth j
	phis[i]: modeshape which goes with mode i evaluated at source depth */
	
{
  /*
   * Calculate Green's function and place the answer in g
   * (previously allocated)
   */

  int i, j;
  complex a, e;

  memset (g, 0, nz * sizeof(complex));
  for (i = 0; i < nmodes; i++) {
    a = complex_mult_component (k[i], range);
    e = complex_div (complex_exp (complex_mult (make_complex(0.0, -1.0), a)), complex_sqrt(a));

    for (j = 0; j < nz; j++) {
      g[j] = complex_add (g[j],
			  complex_mult (complex_mult(phi[i][j], phis[i]), e));
    } /* z loop */
  } /* mode loop */

  for (j = 0; j < nz; j++) 
    g[j] = complex_mult_component (g[j], sqrt(2.0 * M_PI));

  /* All done... */
}

int
spline_modeshapes (int nmodes,
		   int nz_orig, double *z_orig, complex **phi_orig,
		   int nz_new, double *z_new, complex **phi_new)
     /* nmodes: number of modes,
	nz_orig: number of receivers (original),
	z_orig: receiver depths (original),
	phi_orig: modeshapes (original),
	nz_new: number of receivers (new),
	z_new: receiver depths (new),
	phi_new: modeshapes (new) */
{
  /*
   * Splines modeshapes.
   */
  double *spline_component_in, *spline_component_out;
  int i, j;
  double *y2;

  spline_component_in = (double *) calloc (nz_orig, sizeof(double));
  if (spline_component_in == NULL) {
    perror ("calloc spline_component_in in spline_modeshape()");
    return(-1);
  }

  spline_component_out = (double *) calloc (nz_new, sizeof(double));
  if (spline_component_out == NULL) {
    perror ("calloc spline_component_out in spline_modeshape()");
    return(-1);
  }

  for (i = 0; i < nmodes; i++) {
    for (j = 0; j < nz_orig; j++) 
      spline_component_in[j] = complex_real (phi_orig[i][j]);

    y2 = real_spline (z_orig, spline_component_in, nz_orig,
		      9.0e99, 9.0e99);
    if (y2 == NULL) {
      fprintf (stderr, "error performing real spline on real part\n");
      return(-1);
    }

    for (j = 0; j < nz_new; j++) 
      spline_component_out[j] = real_splint(z_orig, spline_component_in,
					    y2, nz_orig, z_new[j]);

    free(y2);

    for (j = 0; j < nz_orig; j++) 
      spline_component_in[j] = complex_imag (phi_orig[i][j]);

    y2 = real_spline (z_orig, spline_component_in, nz_orig,
		      9.0e99, 9.0e99);
    if (y2 == NULL) {
      fprintf (stderr, "error performing real spline on imag part\n");
      return(-1);
    }

    for (j = 0; j < nz_new; j++) 
      phi_new[i][j] = make_complex (spline_component_out[j],
				    real_splint(z_orig, spline_component_in,
						y2, nz_orig, z_new[j]));
    free(y2);
  } /* mode loop */

  free (spline_component_in);
  free (spline_component_out);

  return(0);
}

/* todo: (hazy)

   1. Range independent propagation run in RAM
   2. Same for KRAKEN

   3. Calculate Green's functions at 1, 5, 10, 15, 20 km
      Compare results to make sure they agree.

   4. Create N sensor array, run in reciprocity
   5. Derive weights for each case (100 modes at 250 Hz?!?)

   Store:
   weights
   mode numbers

   6. Write mode-shooting routine.

   7. Write Green's function derivation routine
   (given weights, mode numbers, source and receiver depths, generate Green's function)
*/