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

#ifndef RAM_UTIL
#include "ram-util.h"
#endif

int spline_complex_pressure_range_depth (RAM *, RAM_COMPANION *,int, double *, complex *, double);

int
main (int argc, char **argv)
{
  /*
     1. Reads in RAM output file and companion file
     2. Extracts Green's function at a specific range (splining along the way)
     3. Writes the result to stdout.
  */

  int arg_counter = 0;
  char *ram_prefix, *ram_file;
  double receiver_start, receiver_delta, receiver_stop, range;
  double *z_new;
  int nz_new, i;

  FILE *rfp, *cfp;
  RAM *ram;
  RAM_COMPANION *rc;
  complex *spline_depth;

  if (argc == 1) {
    fprintf (stderr, "usage: %s ram-prefix receiver-start receiver-delta receiver-stop range\n",
	     argv[0]);
    return(-1);
  }

  if (++arg_counter < argc) 
    ram_prefix = argv[arg_counter];
  else {
    fprintf (stderr, "need ram prefix to read in Green's function\n");
    return(-1);
  }


  /* Receiver start */

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

  /* Receiver delta */

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

  /* Receiver stop */

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

  /* Range */

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

  /*
   * End of command line argument processing
   */

  /* Read the RAM output file */

  ram_file = (char *) calloc (strlen(ram_prefix) + BUFFER_LENGTH, sizeof(char));
  if (ram_file == NULL) {
    perror ("calloc ram_file");
    return(-1);
  }
  sprintf (ram_file, "%s.out", ram_prefix);

  rfp = fopen (ram_file, "r");
  if (rfp == NULL) {
    perror ("fopen ram output file");
    fprintf (stderr, "error reading RAM output file [%s]\n",
	     ram_file);
    return(-1);
  }
  ram = read_ram_output(rfp);
  if (ram == NULL) {
    fprintf (stderr, "error reading RAM from output file [%s]\n",
	     ram_file);
    fclose(rfp);
    return(-1);
  }
  fclose(rfp);

  /* Read the companion file */

  sprintf (ram_file, "%s.cmp", ram_prefix);

  cfp = fopen (ram_file, "r");
  if (cfp == NULL) {
    perror ("fopen ram companion file");
    fprintf (stderr, "error opening RAM companion file [%s] for reading\n",
	     ram_file);
    return(-1);
  }
  
  rc = read_ram_companion_file (cfp);
  if (rc == NULL) {
    fprintf (stderr, "error reading RAM companion file [%s]\n",
	     ram_file);
    free(ram_file);
    return(-1);
  }

  fclose(cfp);

  /*
   * Spline through pressure field to get the correct range/depth field
   */
  nz_new = (int) ((receiver_stop - receiver_start) / receiver_delta + 1.000);
  z_new = (double *) calloc (nz_new, sizeof(double));
  if (z_new == NULL) {
    perror ("calloc z_new");
    return(-1);
  }
  for (i = 0; i < nz_new; i++) 
    z_new[i] = receiver_start + i * receiver_delta;

  spline_depth = (complex *) calloc(nz_new, sizeof(complex));
  if (spline_depth == NULL) {
    perror ("calloc spline_depth");
    return(-1);
  }

  if (spline_complex_pressure_range_depth (ram, rc, nz_new, z_new, spline_depth, range) < 0) {
    fprintf (stderr, "error splining complex pressure values\n");
    return(-1);
  }

  /* Print the results */

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

  return(0);
}

int
spline_complex_pressure_range_depth (RAM *ram, RAM_COMPANION *rc, int nz_new, double *z_new, complex *spline_depth, double range)
{
  double *spline_component, *spline_range_real, *spline_depth_real;
  complex *spline_range, scale;
  int i, j;
  double *y2;
  
  /* Perform a complex spline of Green's function */

  spline_component = (double *) calloc (ram->lr, sizeof(double));
  if (spline_component == NULL) {
    perror ("calloc spline_component in spline_complex_pressure_range()");
    return(-1);
  }

  spline_range_real = (double *) calloc (ram->lz, sizeof(double));
  if (spline_range_real == NULL) {
    perror ("caloc spline_range_real in spline_complex_pressure_range()");
    return(-1);
  }

  spline_range = (complex *) calloc (ram->lz, sizeof(complex));
  if (spline_range == NULL) {
    perror ("calloc spline_range in spline_complex_pressure_range()");
    return(-1);
  }

  for (i = 0; i < ram->lz; i++) {

    /* Real part */

    for (j = 0; j < ram->lr; j++) 
      spline_component[j] = complex_real(ram->cf[j][i]);

    y2 = real_spline (ram->r, spline_component, ram->lr, 9.0e99, 9.0e99);
    if (y2 == NULL) {
      fprintf (stderr, "error performing real spline on real part, bailing\n");
      return(-1);
    }

    spline_range_real[i] = real_splint (ram->r, spline_component, y2, ram->lr, range);
    free(y2);

    /* Imaginary part */

    for (j = 0; j < ram->lr; j++) 
      spline_component[j] = complex_imag(ram->cf[j][i]);

    y2 = real_spline (ram->r, spline_component, ram->lr, 9.0e99, 9.0e99);
    if (y2 == NULL) {
      fprintf (stderr, "error performing real spline on real part, bailing\n");
      return(-1);
    }

    spline_range[i] = make_complex (spline_range_real[i], real_splint (ram->r, spline_component, y2, ram->lr, range));
    free(y2);
  }

  free(spline_range_real);

  /* We now have a splined pressure vector at the correct range.  The next step is 
     to spline in the depth direction to get the desired receivers. */

  /* Establish the receiver array */

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

  spline_component = (double *) realloc (spline_component, ram->lz * sizeof(double));
  if (spline_component == NULL) {
    perror ("realloc spline_component");
    return(-1);
  }

  /* Spline through the new depths */

  /* Real part */

  for (i = 0; i < ram->lz; i++) 
    spline_component[i] = complex_real(spline_range[i]);

  y2 = real_spline (ram->z, spline_component, ram->lz, 9.0e99, 9.0e99);
  if (y2 == NULL) {
    fprintf (stderr, "error performing real spline on real part, bailing\n");
    return(-1);
  }

  for (i = 0; i < nz_new; i++) 
    spline_depth_real[i] = real_splint (ram->z, spline_component, y2, ram->lz, z_new[i]);

  free(y2);

  /* Imaginary part */

  for (i = 0; i < ram->lz; i++) 
    spline_component[i] = complex_imag(spline_range[i]);

  y2 = real_spline (ram->z, spline_component, ram->lz, 9.0e99, 9.0e99);
  if (y2 == NULL) {
    fprintf (stderr, "error performing real spline on real part, bailing\n");
    return(-1);
  }

  for (i = 0; i < nz_new; i++) {
    spline_depth[i] = make_complex (spline_depth_real[i], real_splint (ram->z, spline_component, y2, ram->lz, z_new[i]));

    /* Put in phase and carrier information */

    if (rc->scale_flag == 2) {
      scale = complex_exp (make_complex(0.0, 2.0 * M_PI * (ram->f / rc->c0) * range));

      /* Take into account range spread */

      scale = complex_div_component (scale, sqrt(range));
      spline_depth[i] = complex_mult (scale, spline_depth[i]);
    }
  }

  free(y2);

  free (spline_depth_real);
  free (spline_component);
  free (spline_range);

  return(0);
}