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

This module includes a function which calculates the
variance of the reference layer absolute velocity, averaged
between navigation fixes.  It is intended to be used as part
of a scheme for finding the ADCP clock offset by looking for
glitches in absolute reference layer velocity when the ship
goes on or off station.  Fixes and ship-ref layer velocities
are passed in arrays, and an array of reference layer
absolute velocities is calculated.  Storage for these arrays
is allocated in the calling routine.

This is adapted from REFABS with minimal changes, so some
things are done less directly than they would be if this
program were written from scratch, or optimized.

The gap calculations are left in place, even though they are
not needed here.



Fri  07-22-1988

*/

#define VREFABS
#include "geninc.h"
#include "use_db.h"
#include "ioserv.h"
#include "vstat.h"   /* UNISTAT_TYPE */
#include "cal.h"     /* VELOCITY_TYPE, FIX_TYPE */

#define STAT 1

#if PROTOTYPE_ALLOWED
int add_ref(VELOCITY_TYPE *r, double t0);
static int read_ref(VELOCITY_TYPE *ref);
static int read_fix(FIX_TYPE *fix);
#else
int add_ref();
static int read_ref();
static int read_fix();
#endif

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

static FIX_TYPE fix1, fix0;
static VELOCITY_TYPE r1, r0, integral;
static int first_ref;
static int first_fix;

static FIX_TYPE *fixes;              /* input array of fixes */
static int i_fix;                    /* index */
static int n_fix;                    /* number of elements in array */
static VELOCITY_TYPE *refs;          /* input array of ship-ref layer */
static int i_ref;
static int n_ref;
static int n_refabs;

static VELOCITY_TYPE *refabs;        /* output array of absolute velocities */

static UNISTAT_TYPE refstat;         /* for calculating variance */


/*--------------------------------------------------------------
   Functions:
*/
#include "add_ref.c"

#if PROTOTYPE_ALLOWED
static int read_ref(VELOCITY_TYPE *ref)
#else
static int read_ref(ref)
VELOCITY_TYPE *ref;
#endif
{
   if (i_ref < n_ref)
   {
      *ref = refs[i_ref];   /* structure assignment */
      i_ref++;
      return (0);
   }
   else return (1);   /* end of array -> return 1 */
}

#if PROTOTYPE_ALLOWED
static int read_fix(FIX_TYPE *fix)
#else
static int read_fix(fix)
FIX_TYPE *fix;
#endif
{
   if (i_fix < n_fix)
   {
      *fix = fixes[i_fix];   /* structure assignment */
      i_fix++;
      return (0);
   }
   else return (1);   /* end of array -> return 1 */
}

#if PROTOTYPE_ALLOWED
double vrefabs(VELOCITY_TYPE *arefs, FIX_TYPE *afixes, VELOCITY_TYPE *arefabs,
  int an_ref, int an_fix, int *an_refabs)
#else
double vrefabs(arefs, afixes, arefabs, an_ref, an_fix, an_refabs)
VELOCITY_TYPE *arefs;
FIX_TYPE *afixes;
VELOCITY_TYPE *arefabs;
int an_ref, an_fix, *an_refabs;
#endif
{
   int finished = 0;
   double vel_var;     /* velocity variance; returned by this function */

   /* set global pointers for this module to point to the
      arrays passed in on the argument list, and similarly
      for the integer counts.
   */
   refs = arefs;
   fixes = afixes;
   refabs = arefabs;
   n_ref = an_ref;
   n_fix = an_fix;

   /* initialize flags and counters */
   first_ref = 1;
   first_fix = 1;
   i_fix = 0;
   i_ref = 0;
   n_refabs = 0;

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

   /* read first fix and first ref layer */
   finished = read_fix(&fix1);
   finished += read_ref(&r1);

   /* keep going until we run out of fixes or ref layer data */
   while (!finished)
   {
      if (!first_ref)
         finished = add_ref(&r1, r0.t);
      /* do the following regardless of the first_ref flag */
      r0 = r1;  /* structure assignment; OK for TC */
      /* after initialization, this is the only place a new ref is read: */
      if (!finished)
         finished = read_ref(&r1);
      first_ref = 0;
   }
   calculate_unistat(&refstat);

   vel_var =  (refstat.sumsq[0] + refstat.sumsq[1]);  /* Var(u) + Var(v) */
   free_unistat(&refstat);

   *an_refabs = n_refabs;    /* number of absolute ref layer velocities */
   return(vel_var);
}                       /* vrefabs() */