/***************************************************************

TIMSLIP

This program was begun in July 1988 and used initially for
determining the calibration factors of the WEPOCS II Moana
Wave ADCP.  Since the DAS then was based on the HP 9816 and
recorded time only to the nearest minute, there were unknown
discrepancies between the ADCP time and true time.  This
time offset causes glitches in the calculation of absolute
velocity, primarily during ship accelerations (including
turns).  Therefore, before trying to calibrate the ADCP, we
attempt to determine the time offset at each major
acceleration (ship going on or off station) when there was
GPS data.  These offsets are then applied before the
calibration factors are calculated.

The time offset is determined by trying to minimize some
measure of the acceleration glitch.  These measures can
include the variance of the absolute reference layer
velocity on the ensemble time grid (every 5 minutes) or on
the fix grid (about every minute for raw GPS), or it can be
something like the difference between the absolute reference
layer velocity during the acceleration and the average of
its values before and after the acceleration.

As of Sun  07-31-1988, this program is still in an
experimental state.  Parameters and options are set in the
source code rather than via a control file.

It is expected that components of this program may have more
general usefulness.   This may require rewriting some of the
routines for greater generality or modularity.

Wed  08-10-1988  changed to control file input.
Thu  01-05-1989  switched to the read_ref from refabs, which
                 filters out BADFLOAT values.
89/01/10  Added the ability to work with course changes as
          well as speed changes.  This requires a change in
          the control file structure.

Fri  12-15-1989  Added ability to specify a fixed time shift in
         seconds in place of "yes" or "no" in response to
         question "use_shifted_times";
         Added output columns for delta-u and delta-v, to be
         used in determining the effect, if any, of Schuler oscillations,
         presumed to be proportional to the North acceleration.
         A corresponding change was made in ADCPCAL.M

CONTROL FILE STRUCTURE:


      fix_file_type:          { simple or HIG }
      fix_file:

      reference_file:
      output_file:

      year_base=            { year base for decimal day calculations }
      min_n_fixes=          { minimum number of fixes for calculation
                              of time shift to proceed. }

                             { The jump in velocity will be
                             between indices (n_refs-1)/2
                             and (n_refs-1)/2 +1.  For
                             example, it will be between 6
                             and 7 if n_refs is 13.  Note
                             that all these indices start
                             from 0, C-style, not from 1.
                             However, the lowest index
                             specified must be >= 1, since
                             the time of the start of
                             ensemble 1 is (approximately)
                             the time recorded with ensemble
                             0.}

      n_refs=                { number of ensembles used in calculation
                               of time shift }
      i_ref_l0=              { first and last indices of ensembles used }
      i_ref_l1=              {  in calibration calculation }
      i_ref_r0=              { "l" is before jump, "r" is after }
      i_ref_r1=              { minimum is 1, not zero }


      up_thresh=             { m/s, for jump detection }
      down_thresh=

      turn_speed=            { m/s, for turn detection }
      turn_thresh=           { degrees }

      dtmax=     360         { max seconds between ensembles, to screen gaps }

      tolerance= 5.e-5       { for BRENT, in days }
      grid:                  { fix or ensemble; calculate
                                 reference layer velocities
                                 averaged between fix times,
                                 or over ensemble intervals}
      use_shifted_times?     { yes or no or fixed time shift in seconds;
                               use the best-fit time shift in the calibration }

*/

#include "common.h"
#include "dbext.h"
#include "use_db.h"
#include "ioserv.h"
#include "misc.h"
#include "vstat.h"
#include "cal.h"    /* FIX_TYPE, VELOCITY_TYPE, fix_to_velocity(), ... */

#define TIMSLIP
#define TSHIFTMAX  10.0/MINUTES_PER_DAY

#define NRMAX 120
#define NFMAX 120


#if PROTOTYPE_ALLOWED
static int interp_fix(FIX_TYPE *f);
int calc_ship(VELOCITY_TYPE *refs,int i_ref_0,int i_ref_1,VELOCITY_TYPE *ship);
double abs_ref_grid(VELOCITY_TYPE *refs);
int average_velocity(VELOCITY_TYPE *varray, VELOCITY_TYPE *vavg,
                     int n, double t0, double dtmax);
double abs_ref_grid(VELOCITY_TYPE *refs);
double shifted_abs_ref_grid(double t);
double shifted_vrefabs(double t);
double spike_size(int i0, int i1, int i2, int i3, VELOCITY_TYPE *refs);
int is_jump(VELOCITY_TYPE refs[], int n_refs, double up_thresh, double down_thresh);
int is_turn(VELOCITY_TYPE refs[], int n_refs, double turn_speed, double turn_thresh);
int get_ref(FILE *fp_ref, VELOCITY_TYPE refs[], int n_refs);
int get_fixes(double t0, double t1);
int calc_ddop(VELOCITY_TYPE *refs);
int calc_dship(VELOCITY_TYPE *refs);
int adcp_calib(VELOCITY_TYPE *dop, VELOCITY_TYPE *ship,
					double *amp, double *phase);
double vrefabs(VELOCITY_TYPE *arefs, FIX_TYPE *afixes, VELOCITY_TYPE *arefabs,
               int an_ref, int an_fix, int *an_refabs);

#else
static int interp_fix();
int calc_ship();
double abs_ref_grid();
int average_velocity();
double abs_ref_grid();
double shifted_abs_ref_grid();
double shifted_vrefabs();
double spike_size();
int is_jump();
int is_turn();
int get_ref();
int get_fixes();
int calc_ddop();
int calc_dship();
int adcp_calib();
double vrefabs();
#endif


/*--------------------------------------------------------------
   Global variables for this module:
*/

static int (*read_fix_ptr)();
FILE *fp_ref, *fp_fix, *fp_out;
static int startup = 1;
VELOCITY_TYPE refs[NRMAX];
static VELOCITY_TYPE absrefs[NRMAX];  /* on ref grid */
static VELOCITY_TYPE shifted_refs[NRMAX];
int n_refs;
static FIX_TYPE fixes[NFMAX];
static int n_fixes;
static VELOCITY_TYPE refabs[NFMAX];   /* on fix grid */
static int n_refabs;
int year_base;

double up_thresh, down_thresh;
double turn_speed, turn_thresh;

/* ref indices of left, right ranges for calib: */
static int i_ref_l0, i_ref_l1, i_ref_r0, i_ref_r1;

static VELOCITY_TYPE  dship, ddop;
static double dtmax;                    /* reject longer intervals */

static UNISTAT_TYPE refstat;         /* for calculating variance */


/*--------------------------------------------------------------
   Functions:
*/

/*
   interp_fix: input with f->t set to desired time;
               at return, f->x and f->y are interpolated.

      returns: 1 on error (interpolation time out of array range)
               0 otherwise

      This function is very specific to this application,
      and quite dumb.  It always does a simple linear search
      from the beginning of the array.
*/
#if PROTOTYPE_ALLOWED
static int interp_fix(FIX_TYPE *f)
#else
static int interp_fix(f)
FIX_TYPE *f;
#endif
{
   extern FIX_TYPE fixes[];
   extern int n_fixes;

   int i1 = 1;
   int i0;
   double df, dt;

   while ( i1 < n_fixes && fixes[i1].t < f->t )
      i1++;
   i0 = i1 - 1;
   if (i1 == n_fixes || fixes[i0].t > f->t)
      return 1;  /* error */

   dt = fixes[i1].t - fixes[i0].t;
   df = fixes[i1].x - fixes[i0].x;
   f->x = fixes[i0].x + df/dt * (f->t - fixes[i0].t);

   df = fixes[i1].y - fixes[i0].y;
   f->y = fixes[i0].y + df/dt * (f->t - fixes[i0].t);
   return 0;
}

/*
   calc_ship:
   calculates the average ship velocity over the
   reference velocity ensembles i_ref_0 to i_ref_1.  The
   former must be greater than zero, since the time of the
   first ensemble is the recorded time of the previous
   ensemble.

   RETURNS: 0 on success
            1 otherwise

*/
#if PROTOTYPE_ALLOWED
int calc_ship(VELOCITY_TYPE *refs, int i_ref_0, int i_ref_1, VELOCITY_TYPE *ship)
#else
int calc_ship(refs, i_ref_0, i_ref_1, ship)
VELOCITY_TYPE *refs;
int i_ref_0, i_ref_1;
VELOCITY_TYPE *ship;
#endif
{
   FIX_TYPE pos0, pos1;  /* interpolated positions */

   pos0.t = refs[i_ref_0-1].t;
   if(interp_fix(&pos0))
      return 1;
   pos1.t = refs[i_ref_1].t;
   if(interp_fix(&pos1))
      return 1;
   if(fix_to_velocity(&pos0, &pos1, ship))
      return 1;
   else
      return 0;
}                       /* calc_ship() */


/*--------------------------------------------------------------
   FUNCTION: average_velocity

   PURPOSE: given an array of ADCP ensemble velocities, calculate
            the average over n ensembles.
            t0 is the start time of the first ensemble, usually
            taken as the recorded time of the previous ensemble.
            dtmax is usually set to a little more than the ensemble
            duration; if any interval between ensembles exceeds it,
            the function terminates with an error.

   RETURNS: 0 if successful
            1 otherwise

----------------------------------------------------------------*/

#if PROTOTYPE_ALLOWED
int average_velocity(VELOCITY_TYPE *varray, VELOCITY_TYPE *vavg, int n,
  double t0, double dtmax)
#else
int average_velocity(varray, vavg, n, t0, dtmax)
VELOCITY_TYPE *varray, *vavg;
int n;
double t0, dtmax;
#endif
{
   int i;
   double dt;

   vavg->u = 0;
   vavg->v = 0;
   dt = varray[0].t - t0;

   for (i=0; i<n; i++)
   {
      if (i != 0)
         dt = varray[i].t - varray[i-1].t;
      if (dt > dtmax)
         return(1);  /* error */
      vavg->u += varray[i].u * dt;
      vavg->v += varray[i].v * dt;
   }
   dt = varray[n-1].t - t0;
   vavg->t = 0.5 * (varray[n-1].t + t0);
   vavg->u = vavg->u / dt;
   vavg->v = vavg->v / dt;
   return 0;
}                       /* average_velocity() */


/*--------------------------------------------------------------

   FUNCTION:   abs_ref_grid

      Calculates the absolute reference layer velocity for
      each ADCP ensemble, assuming the ship velocity is
      constant between fixes.  Note that this does not make
      much sense unless the interval between fixes is
      smaller than or equal to the ensemble interval.

   PARAMETERS:

      refs = the array of VELOCITY_TYPE containing the
               reference layer velocities

   RETURNS:

      variance of the reference layer velocities.

--------------------------------------------------------------*/

#if PROTOTYPE_ALLOWED
double abs_ref_grid(VELOCITY_TYPE *refs)
#else
double abs_ref_grid(refs)
VELOCITY_TYPE *refs;
#endif
{
   VELOCITY_TYPE ship;
   double variance;
   int i;

   /* initialize the statistics routines */
   allocate_unistat(&refstat, 2, "abs ref layer", U_VARIANCE);
   zero_unistat(&refstat);

   for (i=1; i<n_refs; i++)
   {
      if (!calc_ship(refs, i, i, &ship))
      {
         absrefs[i].u = ship.u - refs[i].u;
         absrefs[i].v = ship.v - refs[i].v;
         update_unistat_d(&refstat, (double *) &(absrefs[i].u), 2);
      }
      else
         return BADFLOAT;
   }
   calculate_unistat(&refstat);
   variance =  (refstat.sumsq[0] + refstat.sumsq[1]);  /* Var(u) + Var(v) */
   free_unistat(&refstat);
   return variance;
}

/* Calculate the variance of the reference layer velocities
   when the reference layer times have been shifted.
*/
#if PROTOTYPE_ALLOWED
double shifted_abs_ref_grid(double t)
#else
double shifted_abs_ref_grid(t)
double t;
#endif
{
   int i;
   double ref_var;

   for (i=0; i<n_refs; i++)
   {
      shifted_refs[i] = refs[i];
      shifted_refs[i].t += t;
   }
   ref_var = abs_ref_grid(shifted_refs);
   return(ref_var);
}
/*
   This calculates the variance of the reference layer
   velocity calculated between fixes (in the external array
   "fixes"), using the reference layer velocities in the
   array "refs" after copying them to "shifted_refs" where
   their times are increased by the offset "t".  For
   example, if "t" is positive, we increase all the
   reference times; this would correct for an ADCP DAS clock
   that was behind the navigation clock.
*/

#if PROTOTYPE_ALLOWED
double shifted_vrefabs(double t)
#else
double shifted_vrefabs(t)
double t;
#endif
{
   int i;
   extern int n_refabs;
   double ref_var;

   for (i=0; i<n_refs; i++)
   {
      shifted_refs[i] = refs[i];               /* structure assignment */
      shifted_refs[i].t += t;
   }
   ref_var = vrefabs(shifted_refs, fixes, refabs, n_refs, n_fixes, &n_refabs);
   return(ref_var);
}


#if 0
/*xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx

   spike_size returns the square of the difference between
   the velocity in the middle of an interval and the average
   of the velocities on either side.  This can be minimized to
   zero, and is probably not a very good measure of the presence
   of spikes due to ADCP time errors.  It is left here as an
   alternative, but probably should not be used routinely.
*/
#if PROTOTYPE_ALLOWED
double spike_size(int i0, int i1, int i2, int i3, VELOCITY_TYPE *refs)
#else
double spike_size(i0, i1, i2, i3, refs)
int i0, i1, i2, i3;
VELOCITY_TYPE *refs;
#endif
{
   VELOCITY_TYPE v_left, v_right, v_middle;
   VELOCITY_TYPE ship, dop;
   double u, v;

   if (!calc_ship(refs, i0, i1, &ship) &&
       !average_velocity(&refs[i0], &dop, (i1-i0+1), refs[i0-1].t, dtmax))
   {
      v_left.u = ship.u - dop.u;
      v_left.v = ship.v - dop.v;
   }
   else
      return BADFLOAT;
   if (!calc_ship(refs, i1, i2, &ship) &&
       !average_velocity(&refs[i1], &dop, (i2-i1+1), refs[i1-1].t, dtmax))
   {
      v_middle.u = ship.u - dop.u;
      v_middle.v = ship.v - dop.v;
   }
   else
      return BADFLOAT;
   if (!calc_ship(refs, i2, i3, &ship) &&
       !average_velocity(&refs[i2], &dop, (i3-i2+1), refs[i2-1].t, dtmax))
   {
      v_right.u = ship.u - dop.u;
      v_right.v = ship.v - dop.v;
      u = 0.5 * (v_right.u + v_left.u) - v_middle.u;
      v = 0.5 * (v_right.v + v_left.v) - v_middle.v;
      return (u * u  +  v * v);
   }
   else
      return BADFLOAT;
}                       /* spike_size */

#if PROTOTYPE_ALLOWED
double shifted_spike_size(double t)
#else
double shifted_spike_size(t)
double t;
#endif
{
   int i;
   double ref_var;

   for (i=0; i<n_refs; i++)
   {
      shifted_refs[i] = refs[i];
      shifted_refs[i].t += t;
   }
   ref_var = spike_size(i_ref_l0, i_ref_l1, i_ref_r0, i_ref_r1, shifted_refs);
   return(ref_var);
}
/*xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx*/
#endif

/*--------------------------------------------------------------

   FUNCTIONS:
                     calc_ddop
                     calc_dship

      Calculate the difference between doppler velocities
      (ship relative to reference layer) and ship velocities
      (absolute), respectively, in two time ranges specified
      by index ranges into the array refs.  The index ranges
      are inclusive.

   PARAMETERS:

      refs = array of reference layer velocities
             (calc_dship uses only the times from this
             array, and gets the fixes from a global array.)

   RETURN:

      1     error
      0     success

----------------------------------------------------------------*/
#if PROTOTYPE_ALLOWED
int calc_ddop(VELOCITY_TYPE *refs)
#else
int calc_ddop(refs)
VELOCITY_TYPE *refs;
#endif
{
   extern VELOCITY_TYPE  ddop;
   extern double dtmax;
   extern int i_ref_l1, i_ref_l0, i_ref_r1, i_ref_r0;

   int error;
   VELOCITY_TYPE dop0, dop1;

   error = average_velocity(&refs[i_ref_l0], &dop0, (i_ref_l1 - i_ref_l0 + 1),
                            refs[i_ref_l0 - 1].t, dtmax);
   error += average_velocity(&refs[i_ref_r0], &dop1, (i_ref_r1 - i_ref_r0 + 1),
                            refs[i_ref_r0 - 1].t, dtmax);

   if (error)
      return (1);
   else
   {
      ddop.u = dop1.u - dop0.u;
      ddop.v = dop1.v - dop0.v;
      return (0);
   }
}                       /* calc_ddop() */


#if PROTOTYPE_ALLOWED
int calc_dship(VELOCITY_TYPE *refs)
#else
int calc_dship(refs)
VELOCITY_TYPE *refs;
#endif
{
   VELOCITY_TYPE ship_l, ship_r;

   /* Find ship velocity over earlier ("left") and later ("right")
      parts of time range.
   */
   if (!calc_ship(refs, i_ref_l0, i_ref_l1, &ship_l) &&
       !calc_ship(refs, i_ref_r0, i_ref_r1, &ship_r) )
   {
      dship.u = ship_r.u - ship_l.u;
      dship.v = ship_r.v - ship_l.v;
      return 0;
   }
   else
      return 1;
}                       /* calc_dship() */

/*--------------------------------------------------------------

   FUNCTION:   adcp_calib

      This uses the method of Pollard and Read (1988) to
      calculate the misalignment angle of the transducer.
      This is measured counterclockwise from the gyro to to
      the transducer.

   PARAMETERS:

      dop = pointer to difference in velocity (ship relative
            to reference layer), where the difference has
            been taken between two close time intervals
            separated by a significant change in ship's
            velocity.

      ship = pointer to difference in absolute ship velocity
            over the same intervals.

      amp (output) = pointer to amplitude; if >1, the ADCP
            is underestimating the magnitude of the flow;
            correct by multiplying ADCP velocity by
            amplitude.

      phase (output) = pointer to misalignment angle in
            degrees; correct by rotating the measured
            velocity by this angle.  The full correction to
            a measured velocity expressed as a complex
            number, u+iv, is

            U+iV =  (u+iv) * amp * exp(i*phase)

   RETURN:

      1     error
      0     success

----------------------------------------------------------------*/


#if PROTOTYPE_ALLOWED
int adcp_calib(VELOCITY_TYPE *dop, VELOCITY_TYPE *ship, double *amp, double *phase)
#else
int adcp_calib(dop, ship, amp, phase)
VELOCITY_TYPE *dop, *ship;
double *amp, *phase;
#endif
{
   double uuvv, duv, u2v2d, a_denom;

   duv   = (dop->v * ship->u) - (dop->u * ship->v);
   uuvv  = (dop->u * ship->u) + (dop->v * ship->v);
   u2v2d = sqr(dop->u) + sqr(dop->v);

   if (uuvv == 0.0)
      return 1;
   *phase = -atan( duv / uuvv );
   if ((a_denom = (cos(*phase) * u2v2d)) == 0.0)
      return 1;
   *amp   = uuvv / a_denom;
   /* The sign on *amp is opposite to that given by Pollard and
      Read, because here we are using velocity of ship relative
      to water, instead of the reverse.  The only effect is on
      the sign of the amplitude.
      The sign of the phase is reversed so that it agrees with
      the convention of Joyce (1988) and with the way Ramon's
      original rotate program worked: The angle is measured
      counterclockwise from the gyro axis to the transducer axis.
   */


   *phase *= (180.0 / M_PI);
   return 0;
}

#if PROTOTYPE_ALLOWED
int get_fixes(double t0, double t1)
#else
int get_fixes(t0, t1)
double t0, t1;
#endif
{
   int finished;
   static int n;

   /* The pair of fixes bracketing the last time of the
      previous call become the first pair of fixes to
      test for bracketing of the first time of the present
      call.
   */
   if (startup)
   {
      finished  = (*read_fix_ptr)(fp_fix, &(fixes[0]));
      finished |= (*read_fix_ptr)(fp_fix, &(fixes[1]));
      n = 2;
      startup = 0;
   }
   /* Make sure the first fix is before t0.  If not, t0 cannot
      be bracketed, so there is no point in looking through
      more fixes.  The program should look for a later velocity
      jump instead.
   */
   if (fixes[0].t > t0)
      return (0);

   /* Look first in the existing array of fixes for the pair
      bracketing the first time, t0.
   */
   while ( (n>2) && !( (fixes[0].t < t0) && (fixes[1].t >= t0) ) )
   {
      memcpy(&(fixes[0]), &(fixes[1]), (n-1)*sizeof(FIX_TYPE));
      n--;
   }
   /* If the bracketing pair is not in memory, get it from the fix file. */
   while ( !( (fixes[0].t <= t0) && (fixes[1].t > t0) ))
   {
      fixes[0] = fixes[1];
      finished = (*read_fix_ptr)(fp_fix, &(fixes[1]));
      if (finished)
         return (EOF);   /* end of file exit */
   }
   /* Now the first time is bracketed by fixes */

   while ( !(fixes[n-1].t > t1) && (n < NFMAX))
   {
      finished = (*read_fix_ptr)(fp_fix, &(fixes[n]));
      if (finished)
         return (EOF);       /* end of file exit */
      n++;
   }
   return (n);
}                       /* get_fixes() */


/*============================================================*/

#if PROTOTYPE_ALLOWED
void do_it(FILE *fp_cnt)
#else
void do_it(fp_cnt)
FILE *fp_cnt;
#endif
{
   int finished = 0;
   int i;
   int jump_dir;     /* 1 is up, -1 is down in speed */
   double ref_var;   /* reference layer variance */
   double ref_var0, ref_var1;   /* with - and + TSHIFTMAX */
   double ref_varc;   /* with const_tshift_d */
   double min_ref_var;
   double brent();   /* function that minimizes a function of one variable */
   static double tol = 5.e-5,   /* about 5 seconds */
          xmin;    /* parameters for brent */

   double amp, phase;   /* for calibration */

   static FILE_NAME_TYPE ref_filename, fix_filename, out_filename;

   static VELOCITY_TYPE *selected_refs;
   static int use_shift;
   static char shift_spec[20];
   static double const_tshift_d;
   static int const_tshift = 0;

   static int min_n_fixes;

   static int i_fix_file;        /* index into array of function pointers */
   static int i_grid;            /* index into array of function pointers */

   /*---------------------------------------------------------
      Order of elements must be coordinated between the following
      two arrays, since the second names the index into the first.
   */
   static double (*shift_var_func_ptr[])()=
   {
      shifted_vrefabs,
      shifted_abs_ref_grid
   };
   static NAME_LIST_ENTRY_TYPE grid_names[]=
   {
      {"fix",        0},
      {"ensemble",   1},
      {NULL,         0}
   };
   static NAME_LIST_ENTRY_TYPE yes_no[]=
   {
      {"no",        0},
      {"yes",       1},
      {NULL,         0}
   };

   /*------------------------------------------------------*/

   extern NAME_LIST_ENTRY_TYPE fix_file_types[];
   extern int (*read_fix_ptr_array[])();
   static NAME_LIST_TYPE name_list[]=
   {
      fix_file_types,
      grid_names,
      yes_no
   };

   static PARAMETER_LIST_ENTRY_TYPE parameters[] =
   {
      {" fix_file_type:",   " %s",    &i_fix_file,   TYPE_TRANS},
      {" fix_file:",        " %79s",  fix_filename,  TYPE_STRING},
      {" reference_file:",  " %79s",  ref_filename,  TYPE_STRING},
      {" output_file:",     " %79s",  out_filename,  TYPE_STRING},
      {" year_base=",       " %d",    &year_base,    TYPE_INT},
      {" min_n_fixes=",     " %d",    &min_n_fixes,  TYPE_INT},
      {" n_refs=",          " %d",    &n_refs,       TYPE_INT},
      {" i_ref_l0=",        " %d",    &i_ref_l0,     TYPE_INT},
      {" i_ref_l1=",        " %d",    &i_ref_l1,     TYPE_INT},
      {" i_ref_r0=",        " %d",    &i_ref_r0,     TYPE_INT},
      {" i_ref_r1=",        " %d",    &i_ref_r1,     TYPE_INT},
      {" up_thresh=",       " %lf",   &up_thresh,    TYPE_DOUBLE},
      {" down_thresh=",     " %lf",   &down_thresh,  TYPE_DOUBLE},
      {" turn_speed=",      " %lf",   &turn_speed,   TYPE_DOUBLE},
      {" turn_thresh=",     " %lf",   &turn_thresh,  TYPE_DOUBLE},
      {" dtmax=",           " %lf",   &dtmax,        TYPE_DOUBLE},
      {" tolerance=",       " %lf",   &tol,          TYPE_DOUBLE},
      {" grid:",            " %s",    &i_grid,       TYPE_TRANS + 256},
      {" use_shifted_times?"," %s",   shift_spec,    TYPE_STRING},
      { NULL,              NULL,    NULL,          0           }
   };                   /* parameters[] */



   if (get_parameters(fp_cnt, parameters, name_list))  exit(-1);

   if (i_ref_l0 < 1 || i_ref_l1 < i_ref_l0 ||
                       i_ref_r0 < i_ref_l1 || i_ref_r1 < i_ref_r0)
   {
      printf("\nError in i_ref indices; minimum value is 1, they must increase.\n");
      exit(-1);
   }

   if ( (use_shift = get_code(yes_no, shift_spec)) == BADINT)  /* not yes or no */
   {
      check_error( (const_tshift = atoi(shift_spec)) == 0,
            "\nuse_shifted_times? <yes>, <no>, or <t>; <t> is shift in sec.\n");
      use_shift = 1;
      const_tshift_d = (double)const_tshift / (double)SECONDS_PER_DAY;
   }

   selected_refs = (use_shift ? shifted_refs : refs);

   /* If the ensemble grid is used, n_refabs will always be
      the same as n_refs; otherwise, n_refabs will be set
      for each calibration point by shifted_vrefabs().
   */
   n_refabs = n_refs;

   read_fix_ptr = read_fix_ptr_array[i_fix_file];

   fp_ref = check_fopen(ref_filename, "r");
   fp_fix = check_fopen(fix_filename, "r");
   fp_out = check_fopen(out_filename, "w");

   /* initialize by filling the refs array and reading two fixes */
   for (i=0; i<n_refs && !finished; i++)
      finished = get_ref(fp_ref, refs, n_refs);

   fprintf(fp_out, "%%");
   print_parameters(fp_out, parameters, name_list, "\n%");

   fprintf(fp_out,"\n\n%%");
   fprintf(fp_out,"u/d     time #abs #fix  deltaU deltaV variance    vmin shift     amp   phase\n");
   while (!finished)
   {
      finished |= get_ref(fp_ref, refs, n_refs);
      if (!finished &&
         ( (jump_dir = is_jump(refs, n_refs, up_thresh, down_thresh)) != 0
         ||(jump_dir = is_turn(refs, n_refs, turn_speed, turn_thresh)) != 0 ) )
      {
/*
         fprintf(fp_out, "jump: %s\n", ( (jump_dir == 1) ? "up" : "down") );
*/
#if 0
         for (i=0; i<n_refs; i++)
         {
            fprintf(fp_out, "%10.5f %10.3f %10.3f\n",
                     refs[i].t, refs[i].u, refs[i].v);
         }
#endif

         n_fixes = get_fixes(refs[0].t        - TSHIFTMAX,
                             refs[n_refs-1].t + TSHIFTMAX);
         if (n_fixes == EOF)
            finished = 1;

#if 0
         fprintf(fp_out, "%d fixes:\n", n_fixes);
         if (!finished) for (i=0; i<n_fixes; i++)
         {
            fprintf(fp_out, "%10.5f %10.3f %10.3f\n",
                      fixes[i].t, fixes[i].x, fixes[i].y);
         }
#endif
         if (n_fixes > min_n_fixes)
         {
            ref_var  = (*shift_var_func_ptr[i_grid])(0.0);
            ref_var1 = (*shift_var_func_ptr[i_grid])(TSHIFTMAX);
            ref_var0 = (*shift_var_func_ptr[i_grid])(-TSHIFTMAX);

#if 0
            fprintf(fp_out, "variance of %d absolutes: %f\n",
                     n_refabs);
            fprintf(fp_out, "  + TSHIFTMAX: %f,  - TSHIFTMAX: %f\n",
                     ref_var1 = (*shift_var_func_ptr[i_grid])(TSHIFTMAX),
                     ref_var0 = (*shift_var_func_ptr[i_grid])(-TSHIFTMAX));
            fprintf(fp_out, "  + 1 min: %f,  - 1 min: %f\n",
                     (*shift_var_func_ptr[i_grid])( 1.0/MINUTES_PER_DAY),
                     (*shift_var_func_ptr[i_grid])(-1.0/MINUTES_PER_DAY));
#endif
            if (ref_var < ref_var1 && ref_var < ref_var0)
            {
               min_ref_var = brent(-TSHIFTMAX, 0.0, TSHIFTMAX,
                                    shift_var_func_ptr[i_grid], tol, &xmin);
               if (const_tshift != 0)
               {
                  ref_varc  = (*shift_var_func_ptr[i_grid])(const_tshift_d);
               }

               if (calc_ddop(selected_refs))
                  printf("\nerror in calc_ddop");
               else if (calc_dship(selected_refs))
                  printf("\nerror in calc_dship");
               else if (adcp_calib(&ddop, &dship, &amp, &phase))
                  printf("\nerror in adcp_calib");
               else /* everything is ok, so print results on a single line */
               {
                   fprintf(fp_out,
                      " %2d   %7.3f   %2d  %3d  %6.2f %6.2f   %7.3f %7.3f",
                      jump_dir, (refs[(n_refs-1)/2].t), n_refabs, n_fixes,
                      dship.u, dship.v, ref_var, min_ref_var);
                   fprintf(fp_out, "  %4.0f %7.3f %7.3f\n",
                      xmin*SECONDS_PER_DAY, amp, phase);
               }
            }
         }
      }
   }

   fclose(fp_ref);
   fclose(fp_fix);
   fclose(fp_out);
}                       /* do_it() */

