#include <stdio.h>
#include <math.h>		/* atan2(), sqrt(), cos(), sin() */
#include <stdlib.h>		/* calloc() */

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

#ifndef BATH_UTIL_DEF
#include "bath-util.h"
#endif

#ifdef CRAY
#include <mpp/pvm3.h>
#endif

int
main (int argc, char **argv)
{
  /* This program is designed to generate a time series which can be
     used to simulate MFP.  Here are the steps:

     Generate input file:
     - need frequency
     - need bathymetry from source to receiver
     - need source depth

     Process output file:
     - need AEL positions (XYZ)
     - need source position to calculate ranges
     - obtain complex pressure at receivers
     - invert complex pressure to form time series
     - write time series out to file

  */

  int arg_counter = 0;
  char *bottom_receiver_file, *baseline_svp_file, *output_file_prefix, *output_file;
  char *ael_file, *reverse_bathymetry_file, *bathymetry_file;
  int acoid, cn;
  double tide_offset, sx, sy, sdepth, freq;
  double bx, by, sbearing, srange;
  double ri_depth, range_step = -1.0, depth_step = -1.0;
  unsigned int range_independent_flag = 0;
  unsigned int recip_flag = 0;
  unsigned int reverse_bathymetry_flag = 0, bathymetry_file_flag = 0;
  unsigned int aco_flag = 0, grid_flag = 0, scale_flag = 0, project_flag = 0;
  double c0 = 1500.0;

  double grid_min_r, grid_delta_r, grid_max_r, grid_r;
  double grid_min_z, grid_delta_z, grid_max_z, grid_z;
  double grid_bearing;
  int grid_rindex, grid_zindex;
  int grid_nr, grid_nz;

  double recip_bearing, center_range, recip_r, recip_base;
  double cx, cy;
  int num_ael;
  int i, j;

  FILE *sfp, *bfp, *ofp, *afp;

  SVP *s;
  double **p;
  Bathymetric *b;
  AEL *ael;

#ifdef CRAY
  int tid, mype;
#endif

  if (argc == 1) {
    fprintf (stderr, "usage: %s [-ri depth] baseline-svp-file",
	     argv[0]);
    fprintf (stderr, " tide-offset");
    fprintf (stderr, " [-ns] [-rs range-step] [-ds depth-step] [-c0 speed] \n");
    fprintf (stderr, " [-aco acoid bottom-phone-file] [-point x y] \n");
    fprintf (stderr, " [-cart x y depth] \n");
    fprintf (stderr, " [-grid min_r delta_r max_r r_index min_z delta_z max_z z_index bearing] \n");
    fprintf (stderr, " [-recip ael-file cn bearing center_range] \n");
    fprintf (stderr, " [-project bottom-phone-file min_z delta_z z_index bearing range] \n");
    fprintf (stderr, " [-bath bathymetry-file]");
    fprintf (stderr, " [-rbath reverse-bathymetry-file]");
    fprintf (stderr, " freq output-file-prefix\n");

    fprintf (stderr, "\n");

    fprintf (stderr, "This program generates a RAM input file for simulating propagation\n");
    fprintf (stderr, "between a sound source and a vertical line array.  The program can\n");
    fprintf (stderr, "operate in a variety of modes.  First, the [-ri depth] parameter\n");
    fprintf (stderr, "specifies range-independent propagation, with specified bottom depth.\n");
    fprintf (stderr, "If this parameter is not specified, then a range depdendent\n");
    fprintf (stderr, "environment is assumed. \n");
    fprintf (stderr, "\n");
    fprintf (stderr, "The target for the MFP (where the receivers are) can be specified in two ways;\n");
    fprintf (stderr, "either through the base of a single FFP array, using the  [-aco acoid]\n");
    fprintf (stderr, "parameter, or through an arbitrary x,y point [-point x y]. \n");
    fprintf (stderr, "\n");
    fprintf (stderr, "The start of the MFP processing (where the source is) can be specified\n");
    fprintf (stderr, "in two ways.  If this is a 'one shot' propagation, a the [-cart x y depth] can\n");
    fprintf (stderr, "be used to specify the x, y position of the source, as well as the depth\n");
    fprintf (stderr, "of the source.  If this is part of an MFP grid search, the limits of the\n");
    fprintf (stderr, "grid can be specified in the [-grid] line (in range and depth) as well\n");
    fprintf (stderr, "as the current indices into the range and depth grids.  The relative bearing\n");
    fprintf (stderr, "is needed (in radians from the receiver x, y position) as well.\n");

    return(-1);
  }

  if (++arg_counter < argc) {
    if (strcmp(argv[arg_counter], "-ri") == 0) {
      range_independent_flag = 1;
      if (++arg_counter < argc && sscanf(argv[arg_counter], "%lf", &ri_depth) != 1) {
	fprintf (stderr, "depth must be specified for range-independent propagation\n");
	fprintf (stderr, "depth should be a floating point number, in meters\n");
	return(-1);
      }
    } else
      arg_counter--;
  }

  if (++arg_counter < argc)
    baseline_svp_file = argv[arg_counter];
  else {
    fprintf (stderr, "sorry, need baseline SVP file to generate an input file.\n");
    return(-1);
  }

  if (++arg_counter < argc) {
    if (sscanf(argv[arg_counter], "%lf", &tide_offset) != 1) {
      fprintf (stderr, "tide offset should be a floating point number (in meters), not [%s]\n",
	       argv[arg_counter]);
      return(-1);
    }
  } else {
    fprintf (stderr, "sorry, need tide offset to write propagation code input file\n");
    return(-1);
  }

  scale_flag = 1;
  if (++arg_counter < argc && strcmp(argv[arg_counter], "-ns") == 0) 
    scale_flag = 2;
  else 
    arg_counter--;


  if (++arg_counter < argc && strcmp(argv[arg_counter], "-rs") == 0) {
    if (++arg_counter < argc) {
      if (sscanf(argv[arg_counter], "%lf", &range_step) != 1) {
	fprintf (stderr, "range step should be a floating point number, not [%s]\n",
		 argv[arg_counter]);
	return(-1);
      }
    } 
  } else
    arg_counter--;


  if (++arg_counter < argc && strcmp(argv[arg_counter], "-ds") == 0) {
    if (++arg_counter < argc) {
      if (sscanf(argv[arg_counter], "%lf", &depth_step) != 1) {
	fprintf (stderr, "depth step should be a floating point number, not [%s]\n",
		 argv[arg_counter]);
	return(-1);
      }
    } 
  } else
    arg_counter--;


  if (++arg_counter < argc && strcmp(argv[arg_counter], "-c0") == 0) {
    if (++arg_counter < argc) {
      if (sscanf(argv[arg_counter], "%lf", &c0) != 1) {
	fprintf (stderr, "c0 should be a floating point number, not [%s]\n",
		 argv[arg_counter]);
	return(-1);
      }
    } 
  } else
    arg_counter--;

  aco_flag = 0;
  if (++arg_counter < argc && strcmp(argv[arg_counter], "-aco") == 0) {
    aco_flag = 1;
    if (++arg_counter < argc) {
      if (sscanf(argv[arg_counter], "%d", &acoid) != 1) {
	fprintf (stderr, "ACOID should be a decimal integer, not [%s]\n",
		 argv[arg_counter]);
	return(-1);
      }
      if (acoid < 1 || acoid > 5) {
	fprintf (stderr, "ACOID must be 1, 2, 3, 4 or 5, not [%s]\n",
		 argv[arg_counter]);
	return(-1);
      }
    } else {
      fprintf (stderr, "sorry, need ACOID to figure out which array to generate input file for\n");
      return(-1);
    }
    if (++arg_counter < argc)
      bottom_receiver_file = argv[arg_counter];
    else {
      fprintf (stderr, "sorry, need bottom phone file to generate an input file.\n");
      return(-1);
    }
  } else 
    arg_counter--;

  if (++arg_counter < argc && strcmp(argv[arg_counter], "-point") == 0) {
    if (++arg_counter < argc) {
      if (sscanf(argv[arg_counter], "%lf", &bx) != 1) {
	fprintf (stderr, "receiver base x must be a floating point number, not [%s]\n",
		 argv[arg_counter]);
	return(-1);
      }
    } else {
      fprintf (stderr, "sorry, need receiver base x to figure out where to run simulation to\n");
      return(-1);
    }
    if (++arg_counter < argc) {
      if (sscanf(argv[arg_counter], "%lf", &by) != 1) {
	fprintf (stderr, "receiver base y must be a floating point number. not [%s]\n",
		 argv[arg_counter]);
	return(-1);
      }
    } else {
      fprintf (stderr, "sorry, need receiver base y to figure out where to run simulation to\n");
      return(-1);
    }
  } else 
    arg_counter--;

  if (++arg_counter < argc && strcmp(argv[arg_counter], "-cart") == 0) {
    if (++arg_counter < argc) {
      if (sscanf(argv[arg_counter], "%lf", &sx) != 1) {
	fprintf (stderr, "source x must be a floating point number, not [%s]\n",
		 argv[arg_counter]);
	return(-1);
      }
    } else {
      fprintf (stderr, "sorry, need source x to figure out where to run simulation from\n");
      return(-1);
    }
    if (++arg_counter < argc) {
      if (sscanf(argv[arg_counter], "%lf", &sy) != 1) {
	fprintf (stderr, "source y must be a floating point number. not [%s]\n",
		 argv[arg_counter]);
	return(-1);
      }
    } else {
      fprintf (stderr, "sorry, need source y to figure out where to run simulation from\n");
      return(-1);
    }
    if (++arg_counter < argc) {
      if (sscanf(argv[arg_counter], "%lf", &sdepth) != 1) {
	fprintf (stderr, "source depth must be a floating point number. not [%s]\n",
		 argv[arg_counter]);
	return(-1);
      }
    } else {
      fprintf (stderr, "sorry, need source depth to figure out where to run simulation from\n");
      return(-1);
    }
  } else
    arg_counter--;

  grid_flag = 0;
  if (++arg_counter < argc && strcmp(argv[arg_counter], "-grid") == 0) {
    grid_flag = 1;
    if (++arg_counter < argc && sscanf(argv[arg_counter], "%lf", &grid_min_r) != 1) {
      fprintf (stderr, "minimum radius must be a floating point number, not [%s]\n",
	       argv[arg_counter]);
      return(-1);
    } else if (arg_counter >= argc) {
      fprintf (stderr, "need minimum radius to determine trial source position\n");
      return(-1);
    }
    if (++arg_counter < argc && sscanf(argv[arg_counter], "%lf", &grid_delta_r) != 1) {
      fprintf (stderr, "delta radius must be a floating point number, not [%s]\n",
	       argv[arg_counter]);
      return(-1);
    } else if (arg_counter >= argc) {
      fprintf (stderr, "need delta radius to determine trial source position\n");
      return(-1);
    }
    if (++arg_counter < argc && sscanf(argv[arg_counter], "%lf", &grid_max_r) != 1) {
      fprintf (stderr, "maximum radius must be a floating point number, not [%s]\n",
	       argv[arg_counter]);
      return(-1);
    } else if (arg_counter >= argc) {
      fprintf (stderr, "need maximum radius to determine trial source position\n");
      return(-1);
    }
    if (++arg_counter < argc && sscanf(argv[arg_counter], "%d", &grid_rindex) != 1) {
      fprintf (stderr, "radius index must be an integer, not [%s]\n",
	       argv[arg_counter]);
      return(-1);
    } else if (arg_counter >= argc) {
      fprintf (stderr, "need radius index to determine trial source position\n");
      return(-1);
    }

    if (++arg_counter < argc && sscanf(argv[arg_counter], "%lf", &grid_min_z) != 1) {
      fprintf (stderr, "minimum depth must be a floating point number, not [%s]\n",
	       argv[arg_counter]);
      return(-1);
    } else if (arg_counter >= argc) {
      fprintf (stderr, "need minimum depth to determine trial source position\n");
      return(-1);
    }
    if (++arg_counter < argc && sscanf(argv[arg_counter], "%lf", &grid_delta_z) != 1) {
      fprintf (stderr, "delta depth must be a floating point number, not [%s]\n",
	       argv[arg_counter]);
      return(-1);
    } else if (arg_counter >= argc) {
      fprintf (stderr, "need delta depth to determine trial source position\n");
      return(-1);
    }
    if (++arg_counter < argc && sscanf(argv[arg_counter], "%lf", &grid_max_z) != 1) {
      fprintf (stderr, "maximum depth must be a floating point number, not [%s]\n",
	       argv[arg_counter]);
      return(-1);
    } else if (arg_counter >= argc) {
      fprintf (stderr, "need maximum depth to determine trial source position\n");
      return(-1);
    }
    if (++arg_counter < argc && sscanf(argv[arg_counter], "%d", &grid_zindex) != 1) {
      fprintf (stderr, "depth index must be an integer, not [%s]\n",
	       argv[arg_counter]);
      return(-1);
    } else if (arg_counter >= argc){
      fprintf (stderr, "need depth index to determine trial source position\n");
      return(-1);
    }

    if (++arg_counter < argc && sscanf(argv[arg_counter], "%lf", &grid_bearing) != 1) {
      fprintf (stderr, "bearing index must be a floating point number, not [%s]\n",
	       argv[arg_counter]);
      return(-1);
    } else if (arg_counter >= argc) {
      fprintf (stderr, "need bearing to determine trial source position\n");
      return(-1);
    }

  } else
    arg_counter--;

  recip_flag = 0;
  if (++arg_counter < argc && strcmp(argv[arg_counter], "-recip") == 0) {
    recip_flag = 1;

    if (++arg_counter < argc)
      ael_file = argv[arg_counter];
    else {
      fprintf (stderr, "need AEL file to determine where source/receiver pairs are!\n");
      return(-1);
    }

    if (++arg_counter < argc) {
      if (sscanf(argv[arg_counter], "%d", &cn) != 1) {
	fprintf (stderr, "channel must be an integer, not [%s]\n",
		 argv[arg_counter]);
	return(-1);
      }
    } else {
      fprintf (stderr, "sorry, need channel to figure out which source/receiver pair to use.\n");
      return(-1);
    }

    if (++arg_counter < argc) {
      if (sscanf(argv[arg_counter], "%lf", &recip_bearing) != 1) {
	fprintf (stderr, "recip bearing must be a floating point number, not [%s]\n",
		 argv[arg_counter]);
	return(-1);
      }
    } else {
      fprintf (stderr, "sorry, need recip bearing to figure out where to run simulation out to\n");
      return(-1);
    }

    if (++arg_counter < argc) {
      if (sscanf(argv[arg_counter], "%lf", &center_range) != 1) {
	fprintf (stderr, "max range must be a floating point number. not [%s]\n",
		 argv[arg_counter]);
	return(-1);
      }
    } else {
      fprintf (stderr, "sorry, need max range figure out where to run simulation out to\n");
      return(-1);
    }
  } else
    arg_counter--;

  project_flag = 0;
  if (++arg_counter < argc && strcmp(argv[arg_counter], "-project") == 0) {
    project_flag = acoid = 1;

    if (++arg_counter < argc)
      bottom_receiver_file = argv[arg_counter];
    else {
      fprintf (stderr, "need bottom receiver file to determine where to project virtual array from.\n");
      return(-1);
    }

    if (++arg_counter < argc) {
      if (sscanf(argv[arg_counter], "%lf", &grid_min_z) != 1) {
	fprintf (stderr, "min_z should be a floating point number, not [%s]\n",
		 argv[arg_counter]);
	return(-1);
      }
    } else {
      fprintf (stderr, "need min z to figure out where to start virtual array.\n");
      return(-1);
    }
	
    if (++arg_counter < argc) {
      if (sscanf(argv[arg_counter], "%lf", &grid_delta_z) != 1) {
	fprintf (stderr, "delta_z should be a floating point number, not [%s]\n",
		 argv[arg_counter]);
	return(-1);
      }
    } else {
      fprintf (stderr, "need delta z to figure out where to start virtual array.\n");
      return(-1);
    }
	
    if (++arg_counter < argc) {
      if (sscanf(argv[arg_counter], "%d", &grid_zindex) != 1) {
	fprintf (stderr, "z_index should be an integer, not [%s]\n",
		 argv[arg_counter]);
	return(-1);
      }
    } else {
      fprintf (stderr, "need z_index to figure out where to start virtual array.\n");
      return(-1);
    }
	
    if (++arg_counter < argc) {
      if (sscanf(argv[arg_counter], "%lf", &recip_bearing) != 1) {
	fprintf (stderr, "bearing must be a floating point number, not [%s]\n",
		 argv[arg_counter]);
	return(-1);
      }
    } else {
      fprintf (stderr, "sorry, need recip bearing to figure out where to run simulation out to\n");
      return(-1);
    }

    if (++arg_counter < argc) {
      if (sscanf(argv[arg_counter], "%lf", &center_range) != 1) {
	fprintf (stderr, "max range must be a floating point number. not [%s]\n",
		 argv[arg_counter]);
	return(-1);
      }
    } else {
      fprintf (stderr, "sorry, need max range figure out where to run simulation out to\n");
      return(-1);
    }
  } else
    arg_counter--;

  reverse_bathymetry_flag = 0;
  if (++arg_counter < argc && strcmp(argv[arg_counter], "-rbath") == 0) {
    reverse_bathymetry_flag = 1;

    if (++arg_counter < argc)
      reverse_bathymetry_file = argv[arg_counter];
    else {
      fprintf (stderr, "need reverse bathymetry file to determine bathymetry...\n");
      return(-1);
    }
  } else
    arg_counter--;

  bathymetry_file_flag = 0;
  if (++arg_counter < argc && strcmp(argv[arg_counter], "-bath") == 0) {
    bathymetry_file_flag = 1;

    if (++arg_counter < argc)
      bathymetry_file = argv[arg_counter];
    else {
      fprintf (stderr, "need bathymetry file to determine bathymetry...\n");
      return(-1);
    }
  } else
    arg_counter--;

  if (++arg_counter < argc) {
    if (sscanf(argv[arg_counter], "%lf", &freq) != 1) {
      fprintf (stderr, "frequency should be a floating point number, not [%s]\n",
	       argv[arg_counter]);
      return(-1);
    }
  } else {
    fprintf (stderr, "sorry, need frequency for input file\n");
    return(-1);
  }

  if (++arg_counter < argc)
    output_file_prefix = argv[arg_counter];
  else {
    fprintf (stderr, "sorry, need output file name to generate an input file.\n");
    return(-1);
  }

  /* Set range step if not specified */

  if (range_step < 0.0) {
    range_step = c0 / (freq * 5.0);
    fprintf (stderr, "FYI: range step set to %f meters\n",
	     range_step);
  }

  /* Set depth step if not specified */

  if (depth_step < 0.0) {
    depth_step = c0 / (freq * 5.0);
    fprintf (stderr, "FYI: depth step set to %f meters\n",
	     depth_step);
  }

  /* End of command line argument processing */

#ifdef CRAY
  /* Check which PE */

  tid = pvm_mytid();
  mype = pvm_get_PE (tid);
  if (mype != 1)
    exit(0);
  fprintf (stderr, "running on PE: %d\n",
	   mype);
#endif

  /* Begin getting data to generate input file */

  /*
   * Read in the baseline water SVP
   */
  sfp = fopen (baseline_svp_file, "r");
  if (sfp == NULL) {
    perror ("fopen baseline_svp_filen");
    fprintf (stderr, "error opening baseline SVP file [%s] for reading\n",
	     baseline_svp_file);
    return(-1);
  }
					    
  s = read_baseline_water_svp (sfp);
  if (s == NULL) {
    fprintf (stderr, "error reading baseline water svp from [%s]\n",
	     baseline_svp_file);
    return(-1);
  }
  fclose(sfp);

  /*
   * Read in the FFP array location
   */
  if (aco_flag || project_flag) {
    bfp = fopen (bottom_receiver_file, "r");
    if (bfp == NULL) {
      perror ("fopen bottom receiver file");
      fprintf (stderr, "error opening bottom receiver file [%s] for reading\n",
	       bottom_receiver_file);
      return(-1);
    }
    p = read_bottom_receiver_positions (bfp);
    if (p == NULL) {
      fprintf (stderr, "error reading bottom receiver positions from file [%s]\n",
	       bottom_receiver_file);
      return(-1);
    }
    fclose(bfp);

    /* Save the array location */

    bx = p[acoid - 1][0];
    by = p[acoid - 1][1];
  }

  /* If the grid flag were selected, convert the requested grid position to an x, y
     coordinate */

  if (grid_flag) {

    /* Check grid indices */

    grid_r = grid_min_r + grid_delta_r * grid_rindex;
    grid_nr = (grid_max_r - grid_min_r) / grid_delta_r;
    if (grid_r > grid_max_r) {
      fprintf (stderr, "WARNING: overshot maximum range in grid\n");
      fprintf (stderr, "index %d yields range of %f; max range is %f\n",
	       grid_rindex, grid_r, grid_max_r);
      return(-1);
    }
    grid_z = grid_min_z + grid_delta_z * grid_zindex;
    grid_nz = (grid_max_z - grid_min_z) / grid_delta_z;
    if (grid_z > grid_max_z) {
      fprintf (stderr, "WARNING: overshot maximum depth in grid\n");
      fprintf (stderr, "index %d yields depth of %f; max depth is %f\n",
	       grid_zindex, grid_z, grid_max_z);
      return(-1);
    }

    /* Compute source x, y, depth */

    sx = bx + grid_r * cos(grid_bearing);
    sy = by + grid_r * sin(grid_bearing);
    sdepth = grid_z;
  }

  /* If the reciprocating flag were set, 
     - project the selected ACOID/CN to a range, depth position based on AEL
     - set this as sx, sy
  */

  if (recip_flag) {

    /* Read in AEL positions */

    afp = fopen (ael_file, "r");
    if (afp == NULL) {
      perror ("fopen AEL file");
      fprintf (stderr, "error opening AEL file [%s]\n",
	       ael_file);
      return(-1);
    }
    num_ael = 0;
    ael = read_ael_output (afp, &num_ael);
    if (ael == NULL) {
      fprintf (stderr, "error reading in AEL information from [%s]\n",
	       ael_file);
      return(-1);
    }
  
    /* correct the AEL positions for tide/bottom locations */

    (void)correct_ael (ael, num_ael, tide_offset, p);

    /* Solve for the geometric center of the array (x, y) */

    cx = cy = 0;
    for (i = 0; i < NUM_ARRAYS; i++) {
      cx += (p[i][0] / (double) NUM_ARRAYS);
      cy += (p[i][1] / (double) NUM_ARRAYS);
    }

    /* Now project the element locations onto the bearing line */

    (void)project_ael (cx, cy, p, recip_bearing, ael, num_ael);

    /* Set the source x, y based on the projection - first find the actual channel */

    i = 0;
    while (i < num_ael && !(ael[i].acoid == acoid && ael[i].cn == cn)) 
      i++;

    if (i == num_ael) {
      fprintf (stderr, "could not find ACO/CN combination %d/%d in AEL file %s\n",
	       acoid, cn, ael_file);
      return(-1);
    }

    /* ... then find the bottom of the current array (channel 1) */

    j = 0;
    while (j < num_ael && !(ael[j].acoid == acoid && ael[j].cn == 1)) 
      j++;

    if (j == num_ael) {
      fprintf (stderr, "could not find ACO/CN combination %d/%d in AEL file %s\n",
	       acoid, 1, ael_file);
      return(-1);
    }

    /* for bathymetric purposes, sx and sy can be taken as the base of the array */

    sx = p[acoid - 1][0];
    sy = p[acoid - 1][1];
    sdepth = ael[i].z;

    /* the range can be taken as the range from the center (max range)
       to the bottom of the array */

    sbearing = recip_bearing;
    srange = center_range - ael[j].r;

    recip_r = ael[i].r;
    recip_base = ael[j].r;

    /* bx and by are taken as the max range from the center of
       the array down the bearing line */

    bx = cx + center_range * cos(recip_bearing);
    by = cy + center_range * sin(recip_bearing);

    /* we need to fake the grid positions.  The minimum range will
       be 500 meters from the center of the FFP array.  The 
       maximum range will be specified by the user.
       Delta-range specified by the range-step flag. */

    grid_flag = 1;
    grid_min_r = 500.0 - recip_r;
    grid_delta_r = range_step;
    grid_max_r = center_range - recip_r;
    grid_nr = (grid_max_r - grid_min_r) / grid_delta_r;
    grid_rindex = 0;

    /* the z/depth steps are far more mundane... */

    grid_min_z = 1.0;
    grid_delta_z = depth_step;
    grid_max_z = 300.0;
    grid_bearing = recip_bearing;
    grid_nz = (grid_max_z - grid_min_z) / grid_delta_z;
    grid_zindex = 0;

  } else if (project_flag) {
    /*
     * If the projection flag is set,
     * (1) solve for the virtual center of the array
     * (2) construct the virtual array
     * (3) set source to be the zindex element of the virtual array
     */

    /* Solve for the geometric center of the array (x, y) */

    cx = cy = 0;
    for (i = 0; i < NUM_ARRAYS; i++) {
      cx += (p[i][0] / (double) NUM_ARRAYS);
      cy += (p[i][1] / (double) NUM_ARRAYS);
    }

    /* Construct the "virtual" array */

    sx = cx;
    sy = cy;
    sdepth = grid_min_z + grid_zindex *  grid_delta_z;

    /* the range can be taken as the range from the center (max range)
       to the bottom of the array */

    sbearing = recip_bearing;
    srange = center_range;
    recip_r = recip_base =  0;

    /* bx and by are taken as the max range from the center of
       the array down the bearing line */

    bx = cx + center_range * cos(recip_bearing);
    by = cy + center_range * sin(recip_bearing);

    /* we need to fake the grid positions.  The minimum range will
       be 500 meters from the center of the FFP array.  The 
       maximum range will be specified by the user.
       Delta-range specified by the range-step flag. */

    grid_flag = 1;
    grid_min_r = 500.0 - recip_r;
    grid_delta_r = range_step;
    grid_max_r = center_range - recip_r;
    grid_nr = (grid_max_r - grid_min_r) / grid_delta_r;
    grid_rindex = 0;

    /* the z/depth steps are far more mundane... */

    grid_min_z = 1.0;
    grid_delta_z = depth_step;
    grid_max_z = 300.0;
    grid_bearing = recip_bearing;
    grid_nz = (grid_max_z - grid_min_z) / grid_delta_z;
    grid_zindex = 0;

  } else {

    /* Calculate the bearing and range from the source to the receiver array */
    
    sbearing = atan2(by - sy, bx - sx);

    fprintf (stderr, "bearing between receiver array and source is %f radians.\n",
	     sbearing);

    srange = sqrt((by - sy) * (by - sy) + (bx - sx) * (bx - sx));
  }

  /*
   * Obtain bathymetric data from the bottom of the array out to the maximum range
   */
  if (!range_independent_flag) {

    if (reverse_bathymetry_flag) {
      b = reverse_bathymetric_data (reverse_bathymetry_file,
				    sx, sy, bx, by);
    } else if (bathymetry_file_flag) {
      b = read_bathymetric_data (bathymetry_file,
				 sx, sy);
    } else {
      b = obtain_bathymetric_data (sx, sy, bx, by,
				   range_independent_flag); 
    }

    if (b == NULL) {
      fprintf (stderr, "error obtaining bathymetric data between %f, %f and %f, %f\n",
	       sx, sy,
	       bx, by);
      return(-1);
    }
  } else {
    fprintf (stderr, "setting range independent depth to be %f meters\n",
	     ri_depth);
    b = (Bathymetric *) calloc (1, sizeof(Bathymetric));
    if (b == NULL) {
      perror ("calloc b");
      return(-1);
    }
    b->nr = 1;
    b->bp = (BathymetricPoint *) calloc (1, sizeof(BathymetricPoint));
    if (b->bp == NULL) {
      perror ("calloc b->bp");
      return(-1);
    }
    b->bp[0].r = 0;
    b->bp[0].x = sx;
    b->bp[0].y = sy;
    b->bp[0].z = ri_depth;
  }

  /*
   * Compensate for tides
   */
  (void)tide_compensation (b, tide_offset);

  /* Allocate space to hold the output file name */

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

  /* Open the output RAM file */

  ofp = fopen (output_file, "w");
  if (ofp == NULL) {
    perror ("fopen output file");
    fprintf (stderr, "error opening RAM output file [%s] for writing\n",
	     output_file);
    return(-1);
  }

  /* Write the RAM input file */

  (void)write_ram_setup (ofp, b, s, freq, c0,
			 sx, sy, sdepth,
			 bx, by,
			 srange, range_step, sbearing,
			 scale_flag, recip_flag, project_flag,
			 recip_r, recip_base, grid_min_r);

  fclose(ofp);

  /* Write the companion file */

  sprintf (output_file, "%s.cmp",
	   output_file_prefix);
  ofp = fopen (output_file, "w");
  if (ofp == NULL) {
    perror ("fopen output companion file");
    fprintf (stderr, "error opening output companion file [%s] for writing\n",
	     output_file);
    return(-1);
  }
  if (write_companion_file (ofp,
			    acoid,
			    sx, sy, sdepth,
			    bx, by,
			    c0, range_step,
			    grid_flag, scale_flag, recip_flag, project_flag,
			    grid_min_r, grid_delta_r, grid_rindex, grid_nr,
			    grid_min_z, grid_delta_z, grid_zindex, grid_nz,
			    grid_bearing,
			    recip_base, recip_r, cx, cy,
			    center_range, cn) < 0) {
    fprintf (stderr, "error writing companion file [%s]\n",
	     output_file);
    return(-1);
  }

  fclose(ofp);

  free(output_file);

  /* All done */

  return(0);
}