#ifndef cal_included
#include "cal.h"  /* VELOCITY_TYPE, *_PER_DAY, fix_to_velocity(), read_fix_... */
#endif
/*--------------------------------------------------------------

   FUNCTIONS: add_ref(), add_ref_piece()

   Together these calculate the reference layer velocity
   between fixes and write out the results to the output
   file.  They were originally written as a single recursive
   function, which was more elegant but raised the risk of
   stack overflow.

   PARAMETERS:    r     pointer to VELOCITY_TYPE containing
                        the ship velocity relative to a
                        reference layer, and end
                        time, for an ensemble

                  t0    time in decimal days of the presumed
                        start of the ensemble

   RETURN:        0     if successful
                  1     if EOF or error has been encountered
                 -1     add_ref_piece only: repeat the function

   CALLED BY:  REFABS, vrefabs()
   NOTE:   This file uses preprocessors to distinguish between code
           specific to REFABS, REFABSBT, and VREFABS.  It must therefore
           be #included and re-compiled depending on which of these
           programs it is used with.

   UPDATES:  from REFABSBT.C

   Wed  05-03-1989  Added warning of out-of-sequence fixes.

*/

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

#include "read_fix.c"

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

   finished = add_ref_piece(r, t0);
   while (finished == -1)
      finished = add_ref_piece(r, fix0.t);
   return(finished);
}

#if PROTOTYPE_ALLOWED
void accumulate(VELOCITY_TYPE *r, double dt);
void reset_accum(void);
#ifdef VREFABS
static int read_fix(FIX_TYPE *fix);
#endif
#else
void accumulate();
void reset_accum();
#ifdef VREFABS
static int read_fix();
#endif
#endif

#if PROTOTYPE_ALLOWED
int add_ref_piece(VELOCITY_TYPE *r, double t0)
#else
int add_ref_piece(r, t0)
VELOCITY_TYPE *r;
double t0;
#endif
{
   double dt, fix_dt;
#ifndef REFABSBT
   VELOCITY_TYPE absolute;
#endif
   VELOCITY_TYPE ship;
   int finished = 0;

   if (r->t >= fix1.t)  /* end this fix interval */
   {
      if (!first_fix)
      {
         dt = (fix1.t - t0) * SECONDS_PER_DAY;
         accumulate(r, dt);

#ifndef VREFABS
         /* find average gap velocities */
         if (gap_accum.t != 0)
         {
            gap_accum.u /= gap_accum.t;
            gap_accum.v /= gap_accum.t;
         }
#endif
#ifndef REFABSBT
         /*  find average reference velocity relative to ship */
         if (integral.t != 0)
         {
            integral.u /= integral.t;
            integral.v /= integral.t;
         }
#endif
/* NOTE:  VREFABS didn't have ff. 2 statements */
         fix_dt = fix1.t-fix0.t;
         if (fix_dt <= 0.0)
            printf("WARNING: fixes out of sequence at %12.7f, %12.7f\n",
                        fix0.t, fix1.t);

         /* find ship velocity between fixes */
         fix_to_velocity(&fix0, &fix1, &ship);
#ifndef REFABSBT
         /* find absolute reference velocity */
         absolute.u = -integral.u + ship.u;
         absolute.v = -integral.v + ship.v;
#endif

#ifdef REFABSBT
         /**** output to file, refabsbt ****/
         fprintf(fp_out, "%12.7f  %9.0f %9.0f  %9.0f %9.0f ",
                        fix1.t,
                        ship.u * fix_dt * SECONDS_PER_DAY,
                        ship.v * fix_dt * SECONDS_PER_DAY,
                        integral.u,
                        integral.v);
         fprintf(fp_out, "  %5.2f %5.2f %5.1f\n",
                        gap_accum.u, gap_accum.v, gap_accum.t / 60.0);
      }
      else  /* first fix; refabsbt */
      {
         fprintf(fp_out, "%12.7f  %9.0f %9.0f  %9.0f %9.0f ",
                        fix1.t, 0.0, 0.0, 0.0, 0.0);
         fprintf(fp_out, "  %5.2f %5.2f %5.1f\n",  0.0, 0.0, 0.0);
#else
#ifdef REFABS
         /**** output to file ****/
         fprintf(fp_out, "%12.7f %11.6f %10.6f %7.3f %7.3f %7.2f",
                        fix0.t, fix0.x, fix0.y, absolute.u, absolute.v,
                        fix_dt * MINUTES_PER_DAY);
         fprintf(fp_out, "  %5.2f %5.2f %5.1f\n",
                        gap_accum.u, gap_accum.v, gap_accum.t / 60.0);
#else
#ifdef VREFABS
         /* output to array */
         refabs[n_refabs].t = fix0.t;
         refabs[n_refabs].u = absolute.u;
         refabs[n_refabs].v = absolute.v;
         n_refabs++;
         /* gather statistics */
         update_unistat_d(&refstat, (double *) &(absolute.u), 2);
#endif
#endif
#endif
      }
      /* start again, a new fix interval, or the first fix interval */
      reset_accum();
      fix0 = fix1;
#ifdef VREFABS
      finished = read_fix(&fix1);
#else
      while ( ((finished = (*read_fix_ptr)(fp_fix, &fix1)) == 0) && fix1.t == fix0.t)
      {
         printf("\nWARNING: duplicate fix at time %f\n", fix1.t);
      }
#endif
      first_fix = 0;
      if (!finished)
         return(-1);    /* need to repeat this function */
   }
   else     /* still within a fix interval; just accumulate sums */
   {
      dt = (r->t - t0) * SECONDS_PER_DAY;
      accumulate(r, dt);
   }
   return(finished);
}                       /* add_ref_piece() */

/*

    FUNCTIONS:  accumulate() and reset_accum()

    CALLED BY:  add_ref_piece()

*/
#if PROTOTYPE_ALLOWED
void accumulate(VELOCITY_TYPE *r, double dt)
#else
void accumulate(r, dt)
VELOCITY_TYPE *r;
double dt;
#endif
{
   integral.u  += r->u * dt;
   integral.v  += r->v * dt;
   integral.t  += dt;
#ifndef VREFABS
   if (gap.u != 0.0)
   {
      gap_accum.u += gap.u * dt;
      gap_accum.v += gap.v * dt;
      gap_accum.t += dt;
   }
#endif
}

#if PROTOTYPE_ALLOWED
void reset_accum(void)
#else
void reset_accum()
#endif
{
   integral.u  = 0;
   integral.v  = 0;
   integral.t  = 0;
#ifndef VREFABS
   gap_accum.u = 0;
   gap_accum.v = 0;
   gap_accum.t = 0;
#endif
}