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

    FILE:  REFABS.C

    This file is the source for both the refabs and refabsbt
    programs.

    REFABS:  The purpose of this program is to calculate the raw absolute
             reference layer velocity averaged between navigation fixes.

             For each fix interval, the output line has the time of
             the FIRST fix, lon and lat of that fix, absolute u and
             v for that fix interval, length of the interval in minutes,
             and three accumulated gap numbers: u, v, and t, also in
             minutes.  If there is no velocity info in the fix interval,
             then the gap length equals the length of the fix interval.

                     Eric Firing

    REFABSBT:  The purpose of this program is to produce the information
               needed for bottom track calibration.

               For each fix interval, the output line has the time of the
               SECOND fix, the x and y displacements in meters between
               fixes, the x and y displacements in meters based on the
               bottom track, and the same 3 gap info items as in REFABS.
               In Matlab, use the cumsum operation down the displacement
               columns before doing the regression.

               Mon  07-03-1989  EF

    UPDATES:

    Fri  12-02-1988  E.F.
    Modified to work with different fix file formats,
    initially including only the HIG format and the simple format.

    Thu  10-05-1989  E.F.
    Fixed bug which led to extrapolation of ADCP reference layer
    velocity before the beginning of ADCP data.

    Sun  95/04/30  E.F.
    Fixed initialization bug found by Erik Gottlieb in
    refabsbt, which caused the first refabsbt (and probably
    first refabs, though no example was found) velocity
    values (second line in refabsbt output) to be incorrect.


--------------------------------------------------
control file structure:

fix_file_type:          { simple or HIG }
reference_file:  [name of nav option file from ADCPSECT]
fix_file:        [name of fix file]
output:          [output file name]
year_base=            { year base for decimal day calculations }

ensemble_length= [in seconds; usually 300]
gap_tolerance=   [in seconds; allow a little, perhaps 60]

*/
#include "common.h"
#include "dbext.h"    /* ADJ_BADFLOAT, YMDHMS_TIME_TYPE */
#include "cal.h"      /* FIX_TYPE, VELOCITY_TYPE */
#include "time_.h"    /* TIMDIF() */
#include "use_db.h"
#include "ioserv.h"

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

static FIX_TYPE fix1, fix0;
static VELOCITY_TYPE r1, r0, rg, integral, gap, gap_accum;
static double interval;      /* ensemble interval, days */
static double min_gap;       /* ignore gaps less than this, days */
static FILE *fp_ref, *fp_fix, *fp_out;
static FILE_NAME_TYPE fn_ref, fn_fix, fn_out;
static int first_ref = 1;
static int first_fix = 1;
static int year_base;
static int i_fix_file;

#if PROTOTYPE_ALLOWED
int add_ref(VELOCITY_TYPE *r, double t0);
#else
int add_ref();
#endif

#include "add_ref.c"

static NAME_LIST_TYPE name_list[]=
{
   {fix_file_types}
};

static PARAMETER_LIST_ENTRY_TYPE refabs_params[] =
{
   {" fix_file_type:"      , " %s"    , &i_fix_file      , TYPE_TRANS},
   {" reference_file:"     , " %79s"  , fn_ref           , TYPE_STRING},
   {" fix_file:"           , " %79s"  , fn_fix           , TYPE_STRING},
   {" output:"             , " %79s"  , fn_out           , TYPE_STRING},
   {" year_base="          , " %d"    , &year_base       , TYPE_INT},
   /* The following are input as seconds, and scaled to days. */
   {" ensemble_length="    , " %lf"   , &interval        , TYPE_DOUBLE},
   {" gap_tolerance="      , " %lf"   , &min_gap         , TYPE_DOUBLE},
   {NULL                   , NULL     , NULL             , 0}
};

#if PROTOTYPE_ALLOWED
void do_it(FILE *fp_cnt)
#else
void do_it(fp_cnt)
FILE *fp_cnt;
#endif
{
   double gap_length;               /* total length of a gap */
   int finished = 0;

   check_error(get_parameters(fp_cnt, refabs_params, name_list),
                        "getting refabs parameters");
   read_fix_ptr = read_fix_ptr_array[i_fix_file];

   fp_ref = check_fopen(fn_ref, "r");
   fp_fix = check_fopen(fn_fix, "r");
   fp_out = check_fopen(fn_out, "w");

   interval /= SECONDS_PER_DAY;
   min_gap  /= SECONDS_PER_DAY;

   finished = read_ref(fp_ref, &r1);
   do
   {
      finished += (*read_fix_ptr)(fp_fix, &fix1);
   } while (!finished && (fix1.t < r1.t - interval));
   /*
      fix1 is now the first fix that is after the START of
      the first ensemble.  If it is also after the END of
      the first ensemble, advance to the first ensemble for
      which this is not true.  Then we will have a fix that
      falls within an ensemble.
   */
   while (!finished && (fix1.t > r1.t))
   {
      finished = read_ref(fp_ref, &r1);
   }

   while (!finished)
   {
      if (!first_ref)
      {
         gap_length = r1.t - r0.t - interval;
         if (gap_length > min_gap)
         {
            /* fill the gap with a fake ensemble */
            rg.u = 0.5 * (r1.u + r0.u);
            rg.v = 0.5 * (r1.v + r0.v);
            rg.t = r1.t - interval;

            /* set the gap velocity jump */
            gap.u = fabs(r1.u - r0.u);
            gap.v = fabs(r1.v - r0.v);

            finished = add_ref(&rg, r0.t);

            /* reset the gap back to zero */
            gap.u = 0;
            gap.v = 0;

            if (!finished)
               finished = add_ref(&r1, rg.t);
         }
         else
         {
            finished = add_ref(&r1, r0.t);
         }
      } /* end "if (!first_ref)" */
      else  /* first_ref: start integrating from the first fix */
      {
         finished = add_ref(&r1, fix1.t);
      }
      r0 = r1;
      if (!finished)
         finished = read_ref(fp_ref, &r1);
      first_ref = 0;
   }

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