#include <stdio.h>
#include <stdlib.h>		/* calloc(), free() */
#include <math.h>		/* fabs() */
#include <ctype.h>		/* isspace() */

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

void
write_kraken_setup (FILE *fp, SVP *s, double f, double c0,
		    double rdepth_start, double rdepth_delta, double rdepth_stop, double bottom_depth,
		    double max_range)
{
  /*
   * Writes a KRAKEN input file (used in the mode projection code)
   */
  int i;
  double min_c;

  /* Line 1: Title */

  fprintf (fp, "'KRAKEN run for mode projection at FFP array'\n");

  /* Line 2: Frequency */

  fprintf (fp, "%f\n",
	   f);

  /* Line 3: Number of Media.  For the SBCX environment, this is limited to 3 */

  fprintf (fp, "3\n");

  /* Line 4: Options: 

     Type of interpolation for the SVP: C for C-linear.
     Type of top boundary condition:    V for vacuum.
     Attenuation units:                 W for dB/wavelength */

  fprintf (fp, "CVW\n");

  /* WATER LAYER */

  /* Medium info.  We'll use 0 for the number of mesh points; this allows KRAKEN
     to calculate the number of mesh points automatically.  Sigma (the RMS roughness)
     is also zero.  The depth of the water column is specified by the user. */

  fprintf (fp, "0 0 %f\n",
	   bottom_depth);

  /* Now write the sound speed profile */

  min_c = s->c[0];
  for (i = 0; i < s->nz; i++) {
    if (s->z[i] < bottom_depth) {
      fprintf (fp, "%f %f /\n",
	       s->z[i],
	       s->c[i]);
    }
    min_c = (s->c[i] < min_c) ? s->c[i] : min_c;
  }

  /* Since the bottom may be variable, assume the SVP is more or less a 
     straight line to the bottom */

  fprintf (fp, "%f %f /\n",
	   bottom_depth, s->c[s->nz - 1]);

  /* SEDIMENT LAYER 1 */

  fprintf (fp, "0 0 %f\n",
	   bottom_depth + 89.0);

  fprintf (fp, "%f 1560.0 0.0 1.85 0.18 0.0\n",
	   bottom_depth);
  fprintf (fp, "%f 1821.0 /\n",
	   bottom_depth + 89.0);

  /* Sediment layer 2 */

  fprintf (fp, "0 0 %f\n",
	   bottom_depth + 389.0);
  fprintf (fp, "%f 1862.0 0.0 1.88 0.03 0.0\n",
	   bottom_depth + 89.0);
  fprintf (fp, "%f 2374.0 /\n",
	   bottom_depth + 389.0);

  /* Bottom halfspace: Acoustic w/RMS roughness of 0.0 */

  fprintf (fp, "A 0.0\n");
  fprintf (fp, "%f 2374.0 0.0 2.03 0.04 0.0\n",
	   bottom_depth + 389.0);

  /* Phase speed limits: We'll set these to be the 
     the minimum sound speed on the low end,
     and the maximum (bottom halfspace) speed on the high end */

  fprintf (fp, "%f %f \n",
	   min_c,
	   2374.0);

  /* Maximum range: Even though we're evaluating this only at the
     FFP array, it is a good idea to account for processing away from the FFP
     array */

  fprintf (fp, "%f\n",
	   max_range);

  /* Source/Receiver depth information.  We'll put a dummy source in,
     and put in the receiver depths */

  fprintf (fp, "1 30\n");
  fprintf (fp, "%d %f %f /\n",
	   (int) ((rdepth_stop - rdepth_start) / rdepth_delta + 1),
	   rdepth_start,
	   rdepth_stop);

  /* All done! */

  return;
}

Mode *
read_all_kraken_modes (char *filename, float *rd, int nrd, char *comp,
		       char swapbytes)
{
  Mode *mode;
  Modeheader *mhead;
  int err;
  int m;
  int *z_index, i, j, depth_flag;

  /* Allocate space for mode and modeheader */

  mode = (Mode *) calloc (1, sizeof(Mode));
  if (mode == NULL) {
    perror ("calloc mode");
    return(NULL);
  }

  mode->nrd = nrd;
  mode->rd = (float *) calloc (mode->nrd, sizeof(float));
  if (mode->rd == NULL) {
    perror ("calloc mode->rd");
    return(NULL);
  }
  memcpy (mode->rd, rd, mode->nrd * sizeof(float));

  mode->comp = (char *) calloc (strlen(comp) + 1, sizeof(char));
  if (mode->comp == NULL) {
    perror ("calloc mode->comp");
    return(NULL);
  }
  memcpy (mode->comp, comp, strlen(comp));
  
  mode->filename = (char *) calloc (strlen(filename) + 1, sizeof(char));
  if (mode->filename == NULL) {
    perror ("calloc mode->filename");
    return(NULL);
  }
  memcpy(mode->filename, filename, strlen(filename));

  /* Read the mode header from the mode file */

  mhead = read_kraken_mode_header (mode->filename, swapbytes);
  if (mhead == NULL) {
    fprintf (stderr, "error reading mode header from file: %s\n",
             filename);
    return(NULL);
  }

  /* Next, compare the depths which are reported in the mode file to 
     the depths which have been requested */

  z_index = (int *) calloc (mode->nrd, sizeof(int));
  if (z_index == NULL) {
    perror ("calloc z_index");
    return(NULL);
  }

  for (i = 0; i < mode->nrd; i++) {
    depth_flag = j = 0;

    do {
      depth_flag = (fabs(mode->rd[i] - mhead->z[j]) <
		    DEPTH_JITTER);
    } while (!depth_flag && ++j < mhead->ntot);

    if (!depth_flag) {
      fprintf (stderr, "Requested depth %f is not located ",
	       mode->rd[i]);
      fprintf (stderr, "in this modefile\n");
      return(NULL);
    } else 
      z_index[i] = j;
  }

  mode->mhead = mhead;
  mode->freq = mhead->freq;

  /* Copy over the title of the plot */

  mode->title = (char *) calloc (strlen (mhead->title) + 1, sizeof(char));
  if (mode->title == NULL) {
    perror ("calloc mode->title");
    return(NULL);
  }
  memcpy (mode->title, mhead->title, strlen(mhead->title));

  /* Copy over the eigenvalues */

  mode->ck = (complex *) calloc (mhead->m, sizeof(complex));
  if (mode->ck == NULL) {
    perror ("calloc mode->ck");
    return(NULL);
  }
  memcpy (mode->ck, mhead->ck, mhead->m * sizeof(complex));

  for (i = 0; i < mhead->m; i++) {
    if (complex_magnitude(mode->ck[i]) > 1e10 || complex_magnitude(mode->ck[i]) < 1e-10)
      fprintf (stderr, "(after memcopy) HEY: mode: %d freq: %f mag: %e\n",
	       i, mode->freq,
	       complex_magnitude(mode->ck[i]));
  }

  /* Allocate space for the eigenfunctions */

  mode->phir = (complex **) calloc (mhead->m, sizeof(complex *));
  if (mode->phir == NULL) {
    perror ("calloc mode->phir");
    return(NULL);
  }


  /* Allocate space for the indices of receiver points */

  mode->w = (float *) calloc (mode->nrd, sizeof(float));
  if (mode->w == NULL) {
    perror ("calloc mode->w");
    return(NULL);
  }

  mode->ird = (int *) calloc (mode->nrd, sizeof(int));
  if (mode->ird == NULL) {
    perror ("calloc ird");
    return(NULL);
  }

  /* Read in the modes */

  for (m = 0; m < mhead->m; m++) {

    err = get_one_kraken_mode (m, mode, z_index);
    if (err < 0) {
      fprintf (stderr, "error getting mode %d, bailing\n", m);
      return (NULL);
    }
  }

  free (z_index);
  return (mode);
}


Modeheader *
read_kraken_mode_header (char *filename, char swapbytes) 
{
  /* This program attempts to parse a mode header */

  Modeheader *mhead;
  FILE *mfp;
  int i, fortran_record_number = 0;
  int err = 1;
  complex ctmp[2];

  mhead = (Modeheader *) calloc (1, sizeof(Modeheader));
  if (mhead == NULL) {
    perror ("calloc mhead");
    return(NULL);
  }
  
  mhead->filename = (char *) calloc (strlen(filename) + 1, sizeof(char));
  if (mhead->filename == NULL) {
    perror ("calloc mhead->filename");
    return(NULL);
  }
  memcpy (mhead->filename, filename, strlen(filename));

  /* Open file */

  mfp = fopen (mhead->filename, "r");
  if (mfp == NULL) {
    perror ("fopen mode file");
    fprintf (stderr, "error opening %s, terminating\n", mhead->filename);
    exit(1);
  }

  /* Get the record length */

  mhead->lrecl = read_fortran_integer (mfp, swapbytes);

  /* PE says that the record length should be multiplied by 4 */
 
  (mhead->lrecl)*=4;

  /* Get the title of the plot */
  
  fread (mhead->title, sizeof(char), 80, mfp);
  i = 79;
  while (i > 0 && isspace(mhead->title[i]))
    mhead->title[i--] = 0;

  /* read the frequency of the plot */

  mhead->freq = read_fortran_float (mfp);

  /* read the number of media */

  mhead->nmedia = read_fortran_integer (mfp, swapbytes);

  /* read the "ntot" variable */

  mhead->ntot = read_fortran_integer (mfp, swapbytes);

  /* Read the "nmat" variable */

  mhead->nmat = read_fortran_integer (mfp, swapbytes);

  /* Allocate space for additional, media dependent variables */

  mhead->n = (int *) calloc (mhead->nmedia, sizeof(int));
  if (mhead->n == NULL) {
    perror ("calloc mhead->n");
    return(NULL);
  }

  mhead->mater = (char **) calloc (mhead->nmedia, sizeof(char *));
  if (mhead->mater == NULL) {
    perror ("calloc mhead->mater");
    return(NULL);
  }

  /* PE suggests 8 characters per media name */

  for (i = 0; i < mhead->nmedia; i++) {
    mhead->mater[i] = (char *) calloc (9, sizeof (char));
    if (mhead->mater[i] == NULL) {
      perror ("calloc mhead->mater[i]");
      return(NULL);
    }
  }

  /* Seek to the next record */

  fseek (mfp, (++fortran_record_number) * (mhead->lrecl), 0);

  /* Start record 2 */

  for (i = 0; i < mhead->nmedia; i++) {
    mhead->n[i] = read_fortran_integer (mfp, swapbytes);
    err = fread (mhead->mater[i], sizeof(char), 8, mfp);
  }
  if (err == 0) {
    perror ("fread");
    return(NULL);
  }

  /* Seek to the next record */

  fseek (mfp, (++fortran_record_number) * (mhead->lrecl), 0);
  
  /* Start record 3 */

  err = fread (&(mhead->bctop), sizeof(char), 1, mfp);
  if (err == 0) {
    perror ("fread");
    return(NULL);
  }
  
  mhead->cpt = read_fortran_scomplex(mfp);
  mhead->cst = read_fortran_scomplex(mfp);
  mhead->rhot = read_fortran_float(mfp);
  mhead->deptht = read_fortran_float(mfp);

  err = fread (&(mhead->bcbot), sizeof(char), 1, mfp);
  if (err == 0) {
    perror ("fread");
    return(NULL);
  }
  
  mhead->cpb = read_fortran_scomplex(mfp);
  mhead->csb = read_fortran_scomplex(mfp);
  mhead->rhob = read_fortran_float(mfp);
  mhead->depthb = read_fortran_float(mfp);

  /* Seek to the next record */

  fseek (mfp, (++fortran_record_number) * (mhead->lrecl), 0);

  /* Start record 4 */

  /* Allocate space for depth and density */

  mhead->depth = (float *) calloc (mhead->nmedia, sizeof(float));
  if (mhead->depth == NULL) {
    perror ("calloc mhead->depth");
    return(NULL);
  }

  mhead->rho = (float *) calloc (mhead->nmedia, sizeof(float));
  if (mhead->rho == NULL) {
    perror ("calloc mhead->rho");
    return(NULL);
  }

  /* Read the data */

  for (i = 0; i < mhead->nmedia; i++) {
    mhead->depth[i] = read_fortran_float (mfp);
    mhead->rho[i] = read_fortran_float (mfp);
  }

  /* Seek to the next record */

  fseek (mfp, (++fortran_record_number) * (mhead->lrecl), 0);

  /* Start record 5 */

  mhead->m = read_fortran_integer (mfp, swapbytes);

  /* Sanity check: are there any modes to read? */

  if (mhead->m == 0) {
    fprintf (stderr, "This file has no modes, bailing\n");
    return(mhead);
  }

  /* Seek to the next record */

  fseek (mfp, (++fortran_record_number) * (mhead->lrecl), 0);

  /* Start record 6 */

  /* Allocate space for z */

  mhead->z = (float *) calloc (mhead->ntot, sizeof(float));
  if (mhead->z == NULL) {
    perror ("calloc mhead->z");
    return(NULL);
  }

  /* Read in z */

  err = fread (mhead->z, sizeof(float), mhead->ntot, mfp);
  if (err == 0) {
    perror ("fread mhead->z");
    return(NULL);
  }

  /* Allocate space for eigenvalues and eigenfunctions */

  mhead->ck = (complex *) calloc (mhead->m, sizeof(complex));
  if (mhead->ck == NULL) {
    perror ("calloc ck");
    return(NULL);
  }

  /* Read the eigenvalues */

  fseek (mfp, (6 + mhead->m) * (mhead->lrecl), 0);
  for (i = 0; i < mhead->m; i++) {
    mhead->ck[i] = read_fortran_scomplex(mfp);

    if (complex_magnitude(mhead->ck[i]) > 1e10 || complex_magnitude(mhead->ck[i]) < 1e-10)
      fprintf (stderr, "HEY: mode: %d freq: %f mag: %e\n",
	       i, mhead->freq,
	       complex_magnitude(mhead->ck[i]));
  }
  /* Next, solve for the top and bottom "k" if acoustic top/bottom */

  if (mhead->bctop == 'A' || mhead->bcbot == 'A') {
    ctmp[0].re = (2 * M_PI * mhead->freq);
    ctmp[0].im = 0.0;
  }

  if (mhead->bctop == 'A') { 
    ctmp[1] = complex_div (ctmp[0], mhead->cpt);
    mhead->ktop2 = complex_mult (ctmp[1], ctmp[1]);
  }
  if (mhead->bcbot == 'A') { 
    ctmp[1] = complex_div (ctmp[0], mhead->cpb);
    mhead->kbot2 = complex_mult (ctmp[1], ctmp[1]);
  }

  /* Close the file */

  fclose(mfp);

  return(mhead);
}

int
get_one_kraken_mode (int m, Mode *mode, int *z_index)
{
  /* This procedure reads in a single eigenfunction and extracts the */
  /* receiver values */

  int i, tough_luck = 0;
  complex *phi;
  complex gamt, gamb;
  FILE *mfp;

  /*
   * Read the mode file to extract the eigenfunction
   */

  mfp = fopen (mode->filename, "r");
  if (mfp == NULL) {
    perror ("fopen mode file");
    fprintf (stderr, "error opening %s, terminating\n",
             mode->filename);
    return(-1);
  }

  fseek (mfp, (6 + m) * mode->mhead->lrecl, 0);

  /* Allocate enough space to hold the eigenfunction */

  phi = (complex *) calloc (mode->mhead->ntot, sizeof(complex));
  if (phi == NULL) {
    perror ("calloc phi");
    return(-1);
  }

  /* Read the eigenfunction in */

  for (i = 0; i < mode->mhead->nmat; i++) {
    phi[i] = read_fortran_scomplex (mfp);
  }

  /* All done with file (for now...) */

  fclose (mfp);

  /* First, check to see if there is an elastic medium in the problem */

  i = 0;
  do {
    tough_luck = (strcmp (mode->mhead->mater[i], "ELASTIC") == 0);
  } while (!tough_luck && ++i < mode->mhead->nmedia);

  /* If there is an elastic medium, extract the component specified */
  /* by the COMP parameter */

  if (tough_luck)
    extract_component (phi, mode->mhead->n, mode->mhead->mater,
                       mode->mhead->nmedia, mode->comp);

  /*
   * Now extract values at receiver depths
   */

  gamt = make_complex (0.0, 0.0);
  gamb = make_complex (0.0, 0.0);

  /* Take into account the boundary conditions, along with Pekeris */
  /* branch cuts. */

  if (mode->mhead->bctop == 'A')
    gamt = pekrt (complex_subtract (complex_real_pow(mode->mhead->ck[m], 2.0),
                                    mode->mhead->ktop2));  

  if (mode->mhead->bcbot == 'A') 
    gamb = pekrt (complex_subtract (complex_real_pow(mode->mhead->ck[m], 2.0),
                                    mode->mhead->kbot2));

  /* Allocate the appropriate space to hold the eigenfunctions */

  mode->phir[m] = (complex *) calloc (mode->nrd, sizeof(complex));
  if (mode->phir[m] == NULL) {
    perror ("calloc mode->phir[m]");
    return(-1);
  }

  /* Now loop through the receiver depths and calculate the */
  /* eigenfunction for that depth */

  for (i = 0; i < mode->nrd; i++) {
    
    if (mode->rd[i] < mode->mhead->deptht)

      /* Receiver is in upper halfspace */

      mode->phir[m][i] = complex_mult (phi[0],
	     complex_exp(complex_mult(complex_mult(make_complex(-1.0, 0.0), gamt),
              make_complex((mode->mhead->deptht - mode->rd[i]), 0.0)) ) );

    else if (mode->rd[i] > mode->mhead->depthb)

      /* Receiver is in lower halfspace */

      mode->phir[m][i] = complex_mult (phi[mode->mhead->ntot - 1],
                 complex_exp (complex_mult(complex_mult
                  (make_complex(-1.0, 0.0), gamb),
                   (make_complex(mode->rd[i] - mode->mhead->depthb,
				 0.0) ))));

    else {

      /* Receiver is not in halfspace */

      /*
      iz = mode->ird[i];
      mode->phir[m][i] = complex_add (phi[iz],
                 complex_mult(make_complex(mode->w[i], 0.0),
                  complex_subtract (phi[iz+1], phi[iz])));
		  */

      mode->phir[m][i] = phi[z_index[i]];
      /*
fprintf (stderr, "i: %d z_index[i]: %d, phi[]: (%e %ei)\n",
	 i, z_index[i],
	 phi[z_index[i]].re,
	 phi[z_index[i]].im);
	 */
    }
  }

#ifdef DEBUG
  fprintf (stderr, "%d (%f): ",
	   m,
	   mode->mhead->freq);

  for (i = 0; i < mode->nrd; i++) {
        fprintf (stderr, "(%e %ei) ",
	     mode->phir[m][i].re,
	     mode->phir[m][i].im);
  }
  fprintf (stderr, "\n");
#endif

  free(phi);
  return(0);
}

void
extract_component (complex *phi, int *n, char **mater, int nmedia, char *comp)
{
  /*
   * For elastic problems where a stress-displacement vector is
   * output, this routine extracs the desired component
   */

  int med, i, j = 0, k = 0;

  for (med = 0; med < nmedia; med++) {

    for (i = 0; i <= n[med]; i++) {
      if (strcmp (mater[med], "ACOUSTIC") == 0) 
        phi[j] = phi[k++];
      else {
        switch (comp[0]) {
        case 'H':
          phi[j] = phi[k];
          break;
        case 'V':
          phi[j] = phi[k+1];
          break;
        case 'T':
          phi[j] = phi[k+2];
          break;
        case 'N':
          phi[j] = phi[k+3];
          break;
        }
 
        k+=4;
      }
      j++;
    }
  }
  return;
}

complex
pekrt (complex z)
{
  complex p;

  if (complex_real(z) > 0.0)
    p = complex_sqrt(z);
  else
    p = complex_mult (make_complex (0.0, 1.0),
		      complex_sqrt (complex_mult (make_complex (-1.0,
								0.0),
						  z))); 

  return(p);
}

void
weight (float *x, int nx, float *xtab, int nxtab, float *w, int *ix)
{
  /* Given:
     
     x:    absicssas
     xtab: points for tabulation
     nx:   number of x points

     Compute:

     w:    weights for linear interpolation
     ix:   indices for linear interpolation

   */

  int ixtab, l, setflag = 0;

  /* Check to see if there is just one X value for interpolation */

  if (nx == 1) {
    w[0] = 0.0;
    ix[0] = 1;
    return;
  }

  /* Locate the indices and compute weights */

  l = 0;

  for (ixtab = 0; ixtab < nxtab; ixtab++) {

    do {

      setflag = 0;

      /* Does [ x[l], x[l+1] ] bracket receiver depth? */

      if (l >= (nx-1) || x[l+1] >= xtab[ixtab]) {

        /* Yes; make note of L and weight for interpolation */

        ix[ixtab] = l;
        w[ixtab] = (xtab[ixtab] - x[l]) / (x[l+1] - x[l]);
        setflag = 1;
      } 

      /* No; advance to the next depth */

    } while (!setflag && ++l < (nx-1));
  }

  return;

}

void
free_mode_header (Modeheader *mhead)
{
  /*
   * This routine properly frees a mode header which is 
   * presented to it
   */

  int i;

  if (mhead == NULL) 
    return;

  free (mhead->filename);
  free (mhead->n);

  for (i = 0; i < mhead->nmedia; i++) 
    free (mhead->mater[i]);
  free (mhead->mater);

  free (mhead->depth);
  free (mhead->rho);
  free (mhead->z);
  free (mhead->ck);
  
  free (mhead);
}

void
free_mode (Mode *mode)
{
  int i;

  if (mode == NULL)
    return;

  free (mode->rd); 
  free (mode->ird); 
  free (mode->w); 
  free (mode->ck); 

  for (i = 0; i < mode->mhead->m; i++)
    free (mode->phir[i]); 
  free (mode->phir); 

  free (mode->title);
  free (mode->comp);
  free (mode->filename);

  free_mode_header (mode->mhead);

  free (mode);
}


int
read_fortran_integer (FILE *fp, char swapbytes)
{
  int i, tempi;

  if (fread (&tempi, sizeof(int), 1, fp) == 0) 
    perror ("reading fortran integer");

  i = (swapbytes) ? swap_long (tempi) : tempi;

  return(i);
}


float
read_fortran_float (FILE *fp)
{
  float tempf;

  if (fread (&tempf, sizeof(float), 1, fp) == 0)
    perror ("reading fortran float");

  return(tempf);
}

complex
read_fortran_scomplex (FILE *fp)
{
  scomplex tempc;
  complex tempc1;

  if (fread (&tempc, sizeof(scomplex), 1, fp) == 0) 
    perror ("reading fortran scomplex");

  tempc1 = scomplex_to_complex (tempc);
  return(tempc1);
}