/***************************************************************
   ADCPSECT.C   (derived from PROFSTAT.C)

   This program does lots of things with ADCP data in
   a CODAS database.  It is still evolving, and I have not
   yet written a manual for it.

   ADCPSECT must be compiled with the LARGE MEMORY MODEL.


                     Eric Firing

Thu  11-03-1988  fixed bug in navigation option, so missing reference
                       layer data is set to BADFLOAT.
Sat  12-03-1988  Changed reference layer calculation to
                        exclude data that are not within the
                        good range as specified by the
                        access variable structure (which
                        counts bins starting from 1).
                        Changed bin_range[] so that it
                        counts bins from 1; subtract 1 when
                        using as an index.  Added time as an
                        option for contouring x-coordinate.
                        Made year_base a required parameter
                        in the adcpsect parameter list,
                        and removed it from the sequential_correlation
                        and navigation option parameter lists.
89/01/18 Changed select_good() to mark velocities up to
            first_good_bin as BADFLOAT.

Tue  02-07-1989  Broke the transport calculation out into a separate
            module.
Sat  02-11-1989  Moved file name generation into a function, new_extension.
            Broke the sequential correlation function into a
            new module, seqcor.c; same for navigation and reference options,
            now in module shipref.c.
Sun  02-12-1989  Moved contour file output to new module: confile.c.
            It still needs to be refined and reorganized.
Sat  02-25-1989  More reorganization, including ctd.c;
            Added the option "nominal_svel".
            Added "integrate_top" to the regridding options,
            to cause the top bin value to be extrapolated to
            the surface, and all integration to start from
            the surface.
Sat  03-11-1989  Added MATfile output to flux.c and to adcpsect.c,
            so that the ensemble mean times and positions
            will go in a file with extension ".mxy", the
            velocity will be in ".muv", and the transports
            in ".mtr".
            ".mxy": contains array "xyt" with a column for
               each ensemble, and three rows: longitude,
               latitude, decimal day; and vectors "z" and
               "zc" with the depth grid and the grid
               interval centers, respectively.
            ".muv": contains array "uv" with two columns for
               each ensemble, first u and then v.  The row
               corresponds to the depth grid.
            ".mtr": contains array "transport" with two
               columns for each ensemble, first the integral
               of u dy, then -v dx.  The vertical integral
               (cumulative summation) has not been done; do
               it in MATLAB with cumsum(transport).  Note
               that the sign of the horizontal integral
               depends on the direction of ship motion.  For
               example, with an eastward current the
               integral will be positive if the ship is
               going north, negative if southbound.

Fri  03-17-1989  Fixed errors that had been made in adding the
               integrate_top option, which affected the gridding
               and editing regardless of which regrid option
               was used.
               Changed the handling of "output:" so that any
               extension is ignored, and the statistics
               output go to a file with extension ".sta".
Wed  05-03-1989  Changed "write_navigation()" in SHIPREF.C
               so that all types of reference calculations
               can be used for .NAV output; previously
               "final_ship_ref" was omitted.  In making the
               change I changed the algorithm for this
               function, making it much more compact.

Tue  01-23-90  Changed time structure variable name from 'end' to
     (JR)      'end_time' because of problems on the Sun.

Thu  08-16-90  Added PROFILE_FLAGS editing to select_good_velocity;
     (JR)      Added flag_mask option to control file.
               Put call to select_good_velocity before call to
               get_ship.velocity

Tue  09-24-91  For the rotate option, set ship heading = mean heading
     (JR)      if available (nonzero for databases loaded after
               9/13/91 correction to loadping), else use last heading

Tue  07-28-92  Added SIGMA_T option to CTD module of adcpsect.
     (FB)

Wed  09-22-93  Changed matlab output file extensions to _uv.mat, _xy.mat
     (JR)      and _tr.mat because Matlab 4 assumes ASCII-formatted file
               for filenames with extension other than .mat
Tue  08-09-94  Added "ascii" output option to derive ASCII subset with
     (PC)      columns time, lon, lat, U & V at each depth layer
Wed  12-07-94  Modified "ascii" output option to also include
     (PC)      mean transducer temperature and ship velocity

Fri  98/07/17  1) default for .sta ouput is none; "mean"
               must be specified explicitly.
               2) New "time_ranges:" value "single" causes
               each individual profile to be written out.
Tue  99/08/24  Added year_base to matlab output file.

2000/02/22 (EF) Added detide: option.  This may be changed later for
               improved flexibility; at present it requires
               a specific ascii file format.

Needed:
Add "dir_position" code.
Add Richardson number calculation, in module richno.c
Add checking of data processing mask!
----------------------------------------------------------------
{control file structure:}
   {Comments and explanations are in curly brackets, and
      examples of control variables are in square brackets.
      Do not include the square brackets (or curly ones) in
      the non-commented (active) part of an actual control
      file.  C-style comments can be used anywhere, so long
      as they are not nested.  Anything that is not in
      either kind of brackets is a key word.}

     dbname: [WEP]
     output: [wepsect]  {no extension;
                        all extensions will be added automatically}
     step_size: [1]  {number of profiles to advance}
     ndepth: [60] {number of depth bins to read and process}
     time_ranges: [combined] { or separate or single )
     year_base= [1988]
     option_list:
        verbose   { prints lots of debugging and auxiliary output }
        underway? [ yes ] limit= [ 3 ]  { limit in m/s }
        pg_min= [ 80 ]
        AGC_margin= [ 7 ]
        skew_limit:
               agc_limit=  [100]
               max_velocity_difference= [ 0.6 ] {m/s}
        reference: [reference_bins] [ 5 ] to  [ 20 ]
                                    { bin range; first bin is 1 }
                   { or one of the following:
                    none
                    final_ship_ref
                    nav_position
                    nav_velocity
                    bottom_track
                    dir_position }    {dir_position not available yet}

        navigation:

                   [reference_bins] [ 5 ] to  [ 20 ]
                                    { bin range; first bin is 1 }
                   { and/or any of the following:
                    final_ship_ref
                    nav_position
                    nav_velocity
                    bottom_track
                    dir_position }    {dir_position not available yet}
                    end      {required to end navigation option list}

        detide: filename: [ tide.asc ]
                    { ascii file with columns:  }
                    {     decimal day           }
                    {     decimal longitude     }
                    {     decimal latitude      }
                    {     predicted tide U, mm/s}
                    {     predicted tide V, mm/s}
                    {     predicted tide h, mm  }
        rotate:  [ 45 ]       { angle in degrees, old x-axis to new x-axis }
                              { directions are counterclockwise }
                              { or substitute "ship_coordinates"}
        ctd:  dbname: [ \wepocs\ctd\c2us ]
        nominal_svel= [1470.0] {default, for RDI; use 1500 for Ametek}
        regrid:  [integrate ] { or integrate_top, interpolate,
                              centered_average, average, per_grid_interval }
          [ depth ]  { or svan, pot_svan, or sigma_theta
                           (if ctd option selected) }
          {must have one and only one of the following: }
          grid      number= [ 50 ] origin= [ 16 ]  increment= [ 8 ]
          grid_list number= [ 5 ] boundaries: [20] [40] [60] [80] [100]
        contour: [ latitude ] { or longitude or time}
                 [ middle ]    { or mean }
            minimum_npts= [ 10 ]
            units= [0.01]
        vector:
            minimum_npts= [ 10 ]
        ascii:                { time, position, u v repeated }
            minimum_percent= [ 0.50 ] { %ensembles in range with good data }
                                      { less than this, program prints bad = 99999 }
               end        {required to end ascii option list}
        transport:
            minimum_npts= [ 10 ]
            units= [1.e6] {for Sverdrups}
            direction= [all] {or left or right}
        sequential_cor:
            bin_range: [3] to [40]  {starting from 1, inclusive}

        statistics:   { followed by list of options, terminated by end }
               [mean]   { and/or any or all of the following:
                          variance, extrema, norm_extrema, indices, all}
               units= [0.01]
               end        {required to end statistics option list}
        flag_mask:  [ bit-mask ] ...        {see adcp/include/adcpflag.h}
               end        {required to end flag mask list}
     end   {required to end adcpsect option list}

     {list of YMDHMS time pairs:}
        {(yy/mm/dd hh:mm:ss) to (yy/mm/dd hh:mm:ss)}
        {etc.}
****************************************************************/

#include "geninc.h"
#include <time.h>
#include "ioserv.h"
#include "vector.h"
#include "use_db.h"
#include "vstat.h"
#include "ocean.h"
#include "pos_to_m.h"
#include "adcp.h"
#include "adcpflag.h"
#include "adcpsect.h"
#include "data_id.h"

#ifndef M_PI
#define M_PI            3.14159265358979323846
#endif
/*--> ADCP general parameters */

YMDHMS_TIME_TYPE start, end_time;       /* for time range */
int verbose = FALSE;            /* flag for debugging output to stdout */
int first_profile;              /* flag: true if first profile in ensemble */

char time_stamp[30];            /* global variable set to start time of
                                 * calculation. */
/* It will be set in do_it, using ctime(time). */

static FILE_NAME_TYPE DBNAME;
static int STEPS;               /* used by write_ensemble() */
static int ndepth;              /* max number of depth bins to process */
static int combined = 0;    /* flag: 1 if all time ranges in one ensemble */
static int separate = 0;    /* flag: 1 if one time range per ensemble */
static int single_ens = 0;      /* flag: 1 if one profile per ensemble */
                                /* (Note: ident "single" conflicts with a
                                   Solaris header file declaration! */
static int rotate = FALSE;      /* flag */
static int on_station = FALSE;  /* flag: select on-station profiles */
extern int use_ctd;             /* flag; ctd data in use */
extern int detide;              /* flag; detiding is active */
static int statistics = FALSE;  /* former default was U_MEAN;  multi-bit flag */
static int vector = FALSE;      /* flag */

extern int transport;           /* flag in flux.c */
extern int seq_cor;             /* flag in seqcor.c */
extern int contour;             /* flag */
static int ascii;               /* flag */

static ADCP_PROFILE_TYPE adcp;

static CTD_PROFILE_TYPE ctd;

/*--> ADCP statistics */
static UNISTAT_TYPE u_stat, v_stat;
static double stat_units = 1.0;
UNISTAT_TYPE adcp_pos_stat;     /* longitude, latitude, decimal day */
UNISTAT_TYPE adcp_ship_stat;    /* ship's velocity u and v */
UNISTAT_TYPE tr_temp_stat;      /* transducer temperatute */
double mid_adcp_pos[3];         /* avg of min and max */
static double *temp_d_array;    /* for conversion to NaN flagging */

/*--> ASCII standard subset file */
static float asc_min_pct = 0.50;/* minimum number for ascii output */

NAME_LIST_ENTRY_TYPE lat_lon_names[] =
{
   {"longitude", LONGITUDE},
   {"latitude", LATITUDE},
   {"time", DDTIME},
   {NULL, 0}
};
static OPTION_TYPE unistat_options[] =
{
   {"end", 0, NULL, NULL},
   {"mean", U_MEAN, logical_or, (char *) &statistics},
   {"variance", U_VARIANCE, logical_or, (char *) &statistics},
   {"extrema", U_EXTREMA, logical_or, (char *) &statistics},
   {"norm_extrema", U_NORM_EXTREMA, logical_or, (char *) &statistics},
   {"indices", U_INDICES, logical_or, (char *) &statistics},
   {"all", U_ALL, logical_or, (char *) &statistics},
   {"units=", TYPE_DOUBLE, op_get_param, (char *) &stat_units},
   {NULL, 0, NULL, NULL}
};

static OPTION_TYPE ascii_options[] =
{
   {"end", 0, NULL, NULL},
   {"minimum_percent=", TYPE_FLOAT, op_get_param, (char *) &asc_min_pct},
   {NULL, 0, NULL, NULL}
};


/*--> rotation parameters */
static double angle_old_to_new; /* angle between old and new coords */
#define  ROTATE_ABSOLUTE  1
#define  ROTATE_SHIP      2

/*--> parameters for profile selection and editing */
static double speed_limit = 0.0;/* criterion for on or off station */
static int pg_min = 0;          /* threshold for data acceptance */
static int amp_margin = 0;      /* threshold for data acceptance */
static int agc_limit = 0;
static UBYTE flag_mask = ALL_BITS;      /* PROFILE_FLAGS bits to consider jr+ */
static double bad_dif = 0.0;

/*--> regridding variables and constants: */
extern float sigma_sva_grid[NDMAX_NG];
extern float adj_depth[NDMAX_NG];
extern float adj_u[NDMAX_NG];
extern float adj_v[NDMAX_NG];
extern int adj_ndepth;          /* new grid; = rg_number if regridding is
                                 * done */
extern int adj_nu;              /* adj_ndepth or one less; npts in array of
                                 * centers */
extern VERTICAL_GRID_TYPE vgrid;
extern int regrid_var;
extern int regrid_method;
extern NAME_LIST_ENTRY_TYPE regrid_methods[];
extern NAME_LIST_ENTRY_TYPE regrid_variables[];

/*
 * parameter in ctd.c used for depth correction.  It is global so that it can
 * be set easily by an option in the main option list. It would more
 * logically be local to ctd.c, and part of a local option list there, but
 * that would require a change to the control file format.
 */
extern double nominal_svel;


static int year_base;
extern PROFILE_REF_TYPE ref;
extern PROFILE_REF_TYPE nav;


/*--> output parameters and variables */
FILE_NAME_TYPE out_filename;    /* global, so other modules can use it to
                                 * construct file names.  In the other
                                 * modules, it must be declared as extern
                                 * char out_filename[]; declaring it as a
                                 * pointer does not work. */
static FILE_NAME_TYPE vec_filename, stat_filename, mat_filename, asc_filename;
static FILE *fp_out, *fp_zxyt_mat, *fp_uv_mat, *fp_vec, *fp_asc;
static int vec_min_npts;        /* minimum number for vector output */
static MAT_ARRAY_FP_TYPE *mat_zxyt_ptr, *mat_uv_ptr;


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

   Parameter input structures for  do_it in ADCPSECT:

*/

NAME_LIST_ENTRY_TYPE time_ranges[] =
{
   {"separate", FALSE},
   {"combined", TRUE},
   {"single", -1},
   {NULL, 0}
};
NAME_LIST_TYPE adcpsect_name_lists[] =
{
   time_ranges
};
static PARAMETER_LIST_ENTRY_TYPE adcpsect_params[] =
{
   {" dbname:", " %79s", DBNAME, TYPE_STRING},
   {" output:", " %79s", out_filename, TYPE_STRING},
   {" step_size:", " %d", &STEPS, TYPE_INT},
   {" ndepth:", " %d", &ndepth, TYPE_INT},
   {" time_ranges:", " %19s", &combined, TYPE_TRANS},
   {" year_base=", " %d", &year_base, TYPE_INT},
   {" option_list:", NULL, NULL, 0},
   {NULL, NULL, NULL, 0}
};

#if PROTOTYPE_ALLOWED
int speed_in_range(ADCP_PROFILE_TYPE * p, double ship_u, double ship_v)
#else
int speed_in_range(p, ship_u, ship_v)
   ADCP_PROFILE_TYPE *p;
   double ship_u, ship_v;
#endif
/* returns TRUE if ship speed is in range */
{
   double speed;

   if (speed_limit == 0.0)
      return (TRUE);            /* no speed checking */

   if (good_float(ship_u))      /* a ship velocity is available; use it */
      speed = sqrt(sqr(ship_u) + sqr(ship_v));
   else if (good_float(p->u[2]))/* otherwise use relative velocity in 3rd bin */
      speed = sqrt(sqr(p->u[2]) + sqr(p->v[2]));
   else
      return (FALSE);           /* profile is not much good anyway */

   if (on_station)              /* select only on-station profiles */
   {
      if (speed < speed_limit)
         return (TRUE);
      else
         return (FALSE);
   }
   else                         /* select underway profiles */
   {
      if (speed >= speed_limit)
         return (TRUE);
      else
         return (FALSE);
   }
}                               /* speed_in_range */

#if PROTOTYPE_ALLOWED
int index_bad_tracking(ADCP_PROFILE_TYPE * p, int agc_limit, double bad_dif)
#else
int index_bad_tracking(p, agc_limit, bad_dif)
   ADCP_PROFILE_TYPE *p;
   int agc_limit;
   double bad_dif;
#endif
{
   int i;
   int itop, n0;
   struct
   {
      double u, v;
   }      pmin, pmax, pdif;
   double magdif = 0.0;         /* magnitude of difference */

   pmin.u = pmin.v = BADFLOAT;
   pmax.u = pmax.v = -BADFLOAT;
   n0 = min_val(p->namp, p->nu);
   for (itop = 1; itop < n0 - 1 && p->amp[itop + 1] > agc_limit; itop++);       /* ef */

   for (i = itop; i < n0 && magdif <= bad_dif; i++)
   {
      if (p->u[i] < ADJ_BADFLOAT)
      {
         pmin.u = min_val(pmin.u, p->u[i]);
         pmin.v = min_val(pmin.v, p->v[i]);
         pmax.u = max_val(pmax.u, p->u[i]);
         pmax.v = max_val(pmax.v, p->v[i]);
         pdif.u = pmax.u - pmin.u;
         pdif.v = pmax.v - pmin.v;
         magdif = sqrt((pdif.u * pdif.u) + (pdif.v * pdif.v));
      }
   }
   if (i == n0)
      return (0);               /* threshold was not reached */
   else
      return (i);               /* first index where threshold was met or
                                 * exceeded */
}                               /* index_bad_tracking() */


#if PROTOTYPE_ALLOWED
void select_good_velocity(ADCP_PROFILE_TYPE * p)
#else
void select_good_velocity(p)
   ADCP_PROFILE_TYPE *p;
#endif
/* this can be expanded to use various criteria */
{
   int i;
   unsigned int npf;            /* jr+ */

   p->namp = -1;                /* flag to show it has not been read *//* ef */

   if (p->access.last_good_bin < MAXSHORT)  /* check added  97/04/25 */
      p->nu = min_val(p->nu, p->access.last_good_bin + 1); /* ef */
      /* This check is needed only because for a time,
      updscan and dbupdate could set last_good_bin to
      MAXSHORT; they have now been changed to set it to its
      default, conf1.num_bins, instead.
      */
   for (i = 1; i < (p->access.first_good_bin); i++)     /* ef */
   {
      p->u[i] = BADFLOAT;
      p->v[i] = BADFLOAT;
   }

   if (pg_min)
   {
      p->npg = 1 + check_dbget(PERCENT_GOOD, p->percent_good + 1, ndepth,       /* ef */
                               "DBGET, percent good");
      p->nu = min_val(p->nu, p->npg);

      /*
       * Reduce nu until all bad points at the end of the profile are
       * eliminated.  The minimum number of points is 1, u[0], which is not
       * read in from the data but may be used to extrapolate to the surface
       * (by setting it equal to the top bin).
       */
      while ((p->percent_good[p->nu - 1] < pg_min) && p->nu > 1)
         p->nu--;

      /* Now flag any bad points in the middle of the profile. */
      for (i = 1; i < p->nu; i++)       /* ef */
      {
         if (p->percent_good[i] < pg_min)
         {
            p->u[i] = BADFLOAT;
            p->v[i] = BADFLOAT;
         }
      }
   }                            /* end pg_min screening */
   if (amp_margin)
   {
      int noise_floor = 0;
      int count = 0;
      int min_amp;

      p->namp = 1 + check_dbget(AMP_SOUND_SCAT, (char *) (p->amp + 1), ndepth,
                                "DBGET, amplitude");
      /* Define the noise floor as average AGC over last 10 bins. */
      for (i = p->namp - 10; i < p->namp; i++)  /* ef */
      {
         if (p->amp[i] != BADUBYTE)
         {
            noise_floor += p->amp[i];
            count++;
         }
      }
      if (count)
         noise_floor /= count;
      min_amp = noise_floor + amp_margin;
      p->nu = min_val(p->nu, p->namp);
      while ((p->amp[p->nu - 1] < min_amp) && p->nu > 1)
         p->nu--;
      for (i = 1; i < p->nu; i++)       /* ef */
      {
         if (p->amp[i] < min_amp)
         {
            p->u[i] = BADFLOAT;
            p->v[i] = BADFLOAT;
         }
      }
   }                            /* end amp_margin screening */
   if (agc_limit)
   {
      int ibad;

      if (p->namp < 0)          /* ef */
         p->namp = 1 + check_dbget(AMP_SOUND_SCAT, (char *) (p->amp + 1), ndepth,
                                   "DBGET, amplitude");
      ibad = index_bad_tracking(p, agc_limit, bad_dif);
      if (ibad)
      {
         p->nu = ibad;
         printf("\nskew error starting at %d  ", ibad);
      }
   }                            /* end skew error screening */
   if (flag_mask && p->nu > 0 &&/* jr+ */
     (0 < (npf = check_dbget(PROFILE_FLAGS, (char *) (p->profile_flags + 1),
                             (p->nu), "DBGET, profile flags"))))        /* use p->nu to
                                                                         * immediately limit npf */
   {
      flag_data(p->u + 1, p->profile_flags + 1, npf, flag_mask);
      flag_data(p->v + 1, p->profile_flags + 1, npf, flag_mask);
   }                            /* end profile_flags screening */
}                               /* select_good_velocity */

#if PROTOTYPE_ALLOWED
void initialize_ensemble(void)
#else
void initialize_ensemble()
#endif
{
   zero_unistat(&u_stat);
   zero_unistat(&v_stat);
   zero_unistat(&adcp_pos_stat);
   zero_unistat(&adcp_ship_stat);
   zero_unistat(&tr_temp_stat);
   adcp.depth[0] = 0.0;
   first_profile = TRUE;
}                               /* initialize_ensemble */

/*
 * 91/09/24 JR- : use mean_heading if available (loaded after 91/09/13
 * correction to loadping) else use last_heading #define SHIP_HEADING
 * adcp.ancillary_2.last_heading
 */
#if PROTOTYPE_ALLOWED
void increment_ensemble(void)
#else
void increment_ensemble()
#endif
{
   int i;
   int nu, nv;
   float *u_ptr, *v_ptr;
   static FIX_TYPE last_fix;
   struct
   {
      double u, v;
   }      ship;
   double ship_heading;

   check_dbget(ANCILLARY_1, (char *) &adcp.ancillary_1, sizeof(ANCILLARY_1_TYPE),
               "DBGET, ancillary_1");
   check_dbget(ANCILLARY_2, (char *) &adcp.ancillary_2, sizeof(ANCILLARY_2_TYPE),
               "DBGET, ancillary_2");
   check_dbget(ACCESS_VARIABLES, (char *) &adcp.access, sizeof(ACCESS_VARIABLES_TYPE),
               "DBGET, access");
   /*
    * We have to get u and v arrays in order to get the ship velocity and
    * check for speed range (in case navigation is not available, or the bin
    * reference option has been chosen).
    */
   nu = check_dbget_f(U, adcp.u + 1, ndepth, "DBGET_F, u");
   nv = check_dbget_f(V, adcp.v + 1, ndepth, "DBGET_F, v");
   adcp.nu = 1 + min_val(nu, nv);
   adcp.u[0] = BADFLOAT;
   adcp.v[0] = BADFLOAT;
   if (detide)
      subtract_tide(&adcp);


   get_ship_velocity(&adcp, &ref, &last_fix, &ship.u, &ship.v);
   if (verbose)
      printf(" ship U= %f V= %f; ", ship.u, ship.v);

   if (speed_in_range(&adcp, ship.u, ship.v))   /* otherwise,skip everything */
   {
      adcp.nd = 1 + check_dbget_f(DEPTH, adcp.depth + 1,        /* ef */
                                  ndepth, "DBGET_F, depth");
      if (use_ctd)              /* We should also check the data processing
                                 * mask. */
         correct_depths(&ctd, adcp.depth, &adcp.nd);
      if (adcp.nd < adcp.nu)
         adcp.nu = adcp.nd;     /* 94/02/10 jr+ */

      select_good_velocity(&adcp);
      if (nav.method != REF_NONE)
         write_navigation(&adcp, &nav);

      get_latlon(&adcp.position);
      update_unistat_d(&adcp_pos_stat, (double *) &adcp.position, 3);
      update_unistat_d(&adcp_ship_stat, (double *) &ship, 2);
      update_unistat(&tr_temp_stat, &adcp.ancillary_1.tr_temp, 1);

      if (ref.method != REF_NONE)
         /* we should check data proc mask also */
      {
         if (good_float(ship.u))
         {
            for (i = 0, u_ptr = adcp.u; i < adcp.nu; i++, u_ptr++)
               if (*u_ptr < ADJ_BADFLOAT)
                  *u_ptr += ship.u;
            for (i = 0, v_ptr = adcp.v; i < adcp.nu; i++, v_ptr++)
               if (*v_ptr < ADJ_BADFLOAT)
                  *v_ptr += ship.v;
         }
         else
         {
            for (i = 0, u_ptr = adcp.u; i < adcp.nu; i++, u_ptr++)
               *u_ptr = BADFLOAT;
            for (i = 0, v_ptr = adcp.v; i < adcp.nu; i++, v_ptr++)
               *v_ptr = BADFLOAT;
         }
      }
      if (rotate)
      {
         if (rotate == ROTATE_SHIP)     /* v is forward, u is starboard */
         {
            if (adcp.ancillary_1.mn_heading == 0.0)
               ship_heading = (double) adcp.ancillary_2.last_heading;
            else
               ship_heading = (double) adcp.ancillary_1.mn_heading;
            angle_old_to_new = -(M_PI / 180) * ship_heading;
         }
         if (verbose)
            printf("  rotation angle: %7.2f radians ", angle_old_to_new);
         rotate_vector(angle_old_to_new, adcp.u, adcp.v, adcp.nu);
      }

      regrid_velocity(&adcp, &ctd);

      if (seq_cor)
         calc_seqcor(adj_u, adj_v, adj_nu, adcp.ddtime);
      update_unistat(&u_stat, adj_u, adj_nu);
      update_unistat(&v_stat, adj_v, adj_nu);
      first_profile = FALSE;
   }

}                               /* increment_ensemble */

#if PROTOTYPE_ALLOWED
void calculate_ensemble(void)
#else
void calculate_ensemble()
#endif
{
   int i;

   calculate_unistat(&u_stat);
   calculate_unistat(&v_stat);
   calculate_unistat(&adcp_pos_stat);
   for (i = 0; i < 3; i++)
   {
      mid_adcp_pos[i] = adcp_pos_stat.npts[i] > 0 ?
         0.5 * (adcp_pos_stat.min[i] + adcp_pos_stat.max[i]) :
         BADFLOAT;
   }
   if (transport)
      calculate_transport(&u_stat, &v_stat, &adcp_pos_stat);

   if (ascii)
   {
      calculate_unistat(&adcp_ship_stat);
      calculate_unistat(&tr_temp_stat);
      if (tr_temp_stat.sum[0] <= ADJ_BADFLOAT)
         if (tr_temp_stat.sum[0] < -10. || tr_temp_stat.sum[0] > 37.)
         {
            printf("\n\n EXTREME TEMPERATURE = %8.2f time=%9.5f\n",
                   tr_temp_stat.sum[0], adcp_pos_stat.sum[DDTIME]);
            /* tr_temp_stat.sum[0] = ADJ_BADFLOAT; */
         }
      if (verbose)
      {
         printf(" \nship_stat=%f sd=%f\n", adcp_ship_stat.sum[0],
                adcp_ship_stat.sumsq[0]);
         printf("\n tr=%8.2f sd=%8.2f\n", tr_temp_stat.sum[0],
                tr_temp_stat.sumsq[0]);
      }
   }

}                               /* calculate_ensemble */

void end_ensemble(void)
{
   if (verbose)
      printf("\nready to calculate\n");
   calculate_ensemble();
   if (verbose)
      printf("\nready to print\n");
   print_ensemble();
}


#if PROTOTYPE_ALLOWED
void print_header(FILE * fp)
#else
void print_header(fp)
   FILE *fp;
#endif
{
   fprintf(fp, "\nADCPSECT  %s\n", time_stamp);

   print_parameters(fp, adcpsect_params, adcpsect_name_lists, "\n");

   if (use_ctd)
      write_ctd_params(fp);

   write_vgrid_params(fp);

   if (ref.method == REF_BINS)
      fprintf(fp, "\nReference layer: bins %d through %d",
              ref.bin_range[0], ref.bin_range[1]);
   else
      fprintf(fp, "\nNo reference layer removed.");
   if (speed_limit != 0.0)
      fprintf(fp, "\nProfiles from %s, speed limit = %f m/s.",
              (on_station ? "on station" : "underway"), speed_limit);
   if (!pg_min)
      fprintf(fp, "\n minimum percent good = %d", pg_min);

   if (!amp_margin)
      fprintf(fp, "\n minimum AGC counts above noise floor = %d", amp_margin);

   if (rotate == ROTATE_ABSOLUTE)
      fprintf(fp, "\n x-axis rotated %f degrees CCW from east.",
              angle_old_to_new * 180.0 / M_PI);
   else if (rotate == ROTATE_SHIP)
      fprintf(fp, "\n Rotation: x-axis is starboard, y is forward.");
   if (detide)
      write_detide_params(fp);
}                               /* print_header() */



#define print_fixed(x,w,p) x>=ADJ_BADFLOAT ? fprintf(fp_vec, " %*s", w, "1E38") : fprintf(fp_vec, " %*.*f", w, p, x)
#if PROTOTYPE_ALLOWED
void print_ensemble(void)
#else
void print_ensemble()
#endif
{
   int i;
   int min_matuv = 1;           /* This could be made a parameter. */

   if (verbose)
      printf("\nentering print_ensemble\n");

   if (contour)
      write_confile(&u_stat, &v_stat, vgrid.center_array,
                    vgrid.increase_down);

   if (vector)
      /*
       * This is a minimal section to put out a file that can be input to a
       * program to plot velocity (or transport) arrows.  The main ASCII file
       * will have to be used as documentation of the parameters.  If
       * adj_ndepth is not small (like a half-dozen), the lines will be very
       * long!  The format for each line is:
       *
       * longitude latitude u1 v1 u2 v2 .......
       *
       * MKS units are used.
       */
   {
      if (adcp_pos_stat.sum[LONGITUDE] <= ADJ_BADFLOAT)
         fprintf(fp_vec, "%9.4f %9.4f",
                 adcp_pos_stat.sum[LONGITUDE], adcp_pos_stat.sum[LATITUDE]);
      else
         fprintf(fp_vec, "     1E38      1E38");
      for (i = 0; i < adj_nu; i++)
      {
         if (u_stat.npts[i] >= vec_min_npts)
         {
            fprintf(fp_vec, "  %6.2f", u_stat.sum[i]);
            fprintf(fp_vec, " %6.2f", v_stat.sum[i]);
         }
         else
         {
            fprintf(fp_vec, "    1E38   1E38");
         }
      }
      fprintf(fp_vec, "\n");
   }

   if (ascii)
      /*
       * MKS units are used.
       */
   {
      /* time as decimal day */
      if (adcp_pos_stat.sum[DDTIME] <= ADJ_BADFLOAT)
         fprintf(fp_asc, "%9.5f", adcp_pos_stat.sum[DDTIME]);
      else                      /* use middle of input time range as
                     * placeholder for missing data */
         fprintf(fp_asc, "%9.5f",
            ((year_day(&start, year_base) + year_day(&end_time, year_base)) / 2.));

      /* write out position as longitude and latitude; missing = 1E38 */
      if (adcp_pos_stat.sum[LONGITUDE] <= ADJ_BADFLOAT)
      {
         fprintf(fp_asc, " %9.4f %8.4f",
            (adcp_pos_stat.sum[LONGITUDE] <= 180 ? adcp_pos_stat.sum[LONGITUDE] :
             adcp_pos_stat.sum[LONGITUDE] - 360.0), adcp_pos_stat.sum[LATITUDE]);
      }
      else
         fprintf(fp_asc, "     1E38      1E38");

      /*
       * write out transducer temperature (deg C) mean and variance missing =
       * 1E38
       */
      if (tr_temp_stat.sum[0] <= ADJ_BADFLOAT)
         fprintf(fp_asc, " %4.1f %5.2f", tr_temp_stat.sum[0],
            sqrt(tr_temp_stat.sumsq[0]));
      else
         fprintf(fp_asc, " 1E38  1E38");

      /* write out ship mean speed and variance; missing = 1E38 (m/sec) */
      if (adcp_ship_stat.sum[0] <= ADJ_BADFLOAT)
         fprintf(fp_asc, " %5.1f %5.2f %5.1f %5.2f",
            adcp_ship_stat.sum[0], sqrt(adcp_ship_stat.sumsq[0]),
            adcp_ship_stat.sum[1], sqrt(adcp_ship_stat.sumsq[1]));
      else
         fprintf(fp_asc, "  1E38  1E38  1E38  1E38");


      /* write out each depth pair of u and v components */
      for (i = 0; i < adj_nu; i++)
      {
         if ( ((float) u_stat.npts[i]) / u_stat.n_increments >= asc_min_pct &&
            u_stat.sum[i] < ADJ_BADFLOAT)
         {
            fprintf(fp_asc, " %5.0f %5.0f",
               u_stat.sum[i] / 0.001, v_stat.sum[i] / 0.001);
         }
         else
         {
            fprintf(fp_asc, " 99999 99999");
         }
      }
      fprintf(fp_asc, "\n");
   }  /* if (ascii) */

   if (statistics)
   {
      /* always print the time and position information */
      if (separate)
         print_time_range(fp_out, &start, &end_time);
      fprintf(fp_out, "\n ADCP mean position: %g lat, %g long",
              adcp_pos_stat.sum[LATITUDE],
              adcp_pos_stat.sum[LONGITUDE]);
      fprintf(fp_out, "\n ADCP middle position: %g lat, %g long, %g days",
              mid_adcp_pos[LATITUDE],
              mid_adcp_pos[LONGITUDE],
              mid_adcp_pos[DDTIME]);
      fprintf(fp_out, "\n ADCP latitude range: %g to %g",
              adcp_pos_stat.min[LATITUDE],
              adcp_pos_stat.max[LATITUDE]);
      fprintf(fp_out, "\n ADCP longitude range: %g to %g",
              adcp_pos_stat.min[LONGITUDE],
              adcp_pos_stat.max[LONGITUDE]);
      if (use_ctd)
      {
         fprintf(fp_out, "\n CTD position: %f  %f    time: ",
                 ctd.position.lat, ctd.position.lon);
         write_ymdhms_time(fp_out, &(ctd.time));
      }
      fprintf(fp_out, "\n No. of ensembles in time range:  %d\n", u_stat.n_increments);
      write_unistat(fp_out, &u_stat, stat_units, vgrid.array, vgrid.label, vgrid.precision);
      write_unistat(fp_out, &v_stat, stat_units, vgrid.array, vgrid.label, vgrid.precision);
      fprintf(fp_out, "\n\n");
   } /* if (statistics) */

   if (transport)
      print_transport(&adcp_pos_stat, vgrid.array, vgrid.label, vgrid.precision);

   tonan_d(temp_d_array, adcp_pos_stat.sum, 3);
   add_to_mat_array(mat_zxyt_ptr, (char *) temp_d_array);
   tonan_mask_d(temp_d_array, u_stat.sum, u_stat.npts, min_matuv, adj_nu);
   add_to_mat_array(mat_uv_ptr, (char *) temp_d_array);
   tonan_mask_d(temp_d_array, v_stat.sum, v_stat.npts, min_matuv, adj_nu);
   add_to_mat_array(mat_uv_ptr, (char *) temp_d_array);

}                               /* print_ensemble */

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

   Option initialization functions:

*/
#if PROTOTYPE_ALLOWED
int init_stats(FILE * fp_cnt)
#else
int init_stats(fp_cnt)
   FILE *fp_cnt;
#endif
{
   check_error((execute_options(fp_cnt, unistat_options, ECHO) <= 0),
               "reading statistics options");
   return (0);
}

#if PROTOTYPE_ALLOWED
int init_uw(FILE * fp_cnt)
#else
int init_uw(fp_cnt)
   FILE *fp_cnt;
#endif
{
   char sub_option[5];
   check_error((fscanf(fp_cnt, " %4s limit= %lf",
                    sub_option, &speed_limit) != 2), "reading speed limit");
   if (strcmp(sub_option, "yes") == 0)
      on_station = FALSE;
   else if (strcmp(sub_option, "no") == 0)
      on_station = TRUE;
   else
      check_error(1, "underway? yes or no");
   return (0);
}

#if PROTOTYPE_ALLOWED
int init_rotate(FILE * fp_cnt)
#else
int init_rotate(fp_cnt)
   FILE *fp_cnt;
#endif
{
   char sub_option[20];
   check_error((fscanf(fp_cnt, " %19s", sub_option) != 1),
               "reading rotation option");
   if (!strcmp(sub_option, "ship_coordinates"))
   {
      rotate = ROTATE_SHIP;
   }
   else                         /* should be absolute rotation angle */
   {
      check_error( (sscanf(sub_option, " %lf", &angle_old_to_new) != 1),
         "reading rotation angle");
      angle_old_to_new *= M_PI/180.0;
      rotate = ROTATE_ABSOLUTE;
   }
   return(0);
}

#if PROTOTYPE_ALLOWED
int init_vector(FILE *fp_cnt)
#else
int init_vector(fp_cnt)
   FILE *fp_cnt;
#endif
{
   check_error( (fscanf(fp_cnt, " minimum_npts= %d", &vec_min_npts) != 1),
      "reading vector parameter");
   vector = TRUE;
   /* Make a vector output file name, and open the file. */
   new_extension(vec_filename, out_filename, ".vec");
   fp_vec = check_fopen(vec_filename, "w");
   return (0);
}

#if PROTOTYPE_ALLOWED
int init_ascii(FILE * fp_cnt)
#else
int init_ascii(fp_cnt)
   FILE *fp_cnt;
#endif
{
   int i;

   check_error((execute_options(fp_cnt, ascii_options, ECHO) <= 0),
               "reading ascii options");
   ascii = TRUE;

   /* Make a ascii output file name, and open the file. */
   new_extension(asc_filename, out_filename, ".sub");
   fp_asc = check_fopen(asc_filename, "w");

   /*
    * Scan each block file to find maximum number of levels and minimum start
    * level
    */

   /* Make header */
   for (i = 1; i < NDMAX_NG; i++)
      if (adj_depth[i] == 0.0)
         break;

   /* 2002/01/17 Pierre Jaccard
    *            Prepend a '%' to the header
    */
   if (ref.method == REF_FINAL_SHIP)
      fprintf(fp_asc, "%%sac_id=%s yr_base=%d start_lev=%3.0fm num_lev=%3d absolute\n",
              out_filename, year_base, adj_depth[0], i);
   else
      fprintf(fp_asc, "%%sac_id=%s yr_base=%d start_lev=%3.0fm num_lev=%3d relative\n",
              out_filename, year_base, adj_depth[0], i);

   return (0);
}

#if PROTOTYPE_ALLOWED
int init_skew(FILE * fp_cnt)
#else
int init_skew(fp_cnt)
   FILE *fp_cnt;
#endif
{
   check_error((fscanf(fp_cnt, " agc_limit= %d max_velocity_difference= %lf",
                       &agc_limit, &bad_dif) != 2),
               "reading skew parameters");
   return (0);
}

#if PROTOTYPE_ALLOWED
int init_mask(FILE * fp_cnt)    /* jr+ */
#else
int init_mask(fp_cnt)           /* jr+ */
   FILE *fp_cnt;
#endif
{
   char flag_name[40];
   int bit_flag;

   flag_mask = 0;
   do
   {
      check_error(!fscanf(fp_cnt, " %s", flag_name), "reading flag_mask");
      if (!strcmp(flag_name, "end"))
         break;
      check_error((BADINT == (bit_flag = get_code(adcp_profflags, flag_name))),
                  "unrecognized flag_mask");
      flag_mask |= (UBYTE) bit_flag;    /* cumulate flag bit-by-bit */
   } while (1);
   return (0);
}

static OPTION_TYPE adcpsect_options[] =
{
   {"end", 0, NULL, (char *) NULL},
   {"reference:", 1, init_ref, (char *) NULL},
   {"navigation:", 2, init_nav, (char *) NULL},
   {"pg_min=", TYPE_INT, op_get_param, (char *) &pg_min},
   {"verbose", 4, set_code, (char *) &verbose},
   {"underway?", 5, init_uw, (char *) NULL},
   {"rotate:", 6, init_rotate, (char *) NULL},
   {"ctd:", 7, init_ctd, (char *) NULL},
   {"nominal_svel=", TYPE_DOUBLE, op_get_param, (char *) &nominal_svel},
   {"regrid:", 8, init_regrid, (char *) NULL},
   {"statistics:", 9, init_stats, (char *) NULL},
   {"contour:", 10, init_contour, (char *) NULL},
   {"vector:", 11, init_vector, (char *) NULL},
   {"ascii:", 16, init_ascii, (char *) NULL},
   {"transport:", 13, init_transport, (char *) NULL},
   {"sequential_cor:", 14, init_seqcor, (char *) NULL},
   {"AGC_margin=", TYPE_INT, op_get_param, (char *) &amp_margin},
   {"skew_limit:", 12, init_skew, (char *) NULL},
   {"flag_mask:", 15, init_mask, (char *) NULL},
   {"detide:", 17, init_detide, (char *) NULL},
   {NULL, 0, NULL, NULL}
};

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

Main routine: read the control file, process the data
accordingly.
*/

#if PROTOTYPE_ALLOWED
void do_it(FILE * fp_cnt)
#else
void do_it(fp_cnt)
   FILE *fp_cnt;
#endif
{
   int IERR = 0, nbin, i;
   time_t clock_start, clock_end;
   int run_time;
   /* extern long timezone;   *//* C library global variables */
   /* extern int daylight;   -- not needed & machine-dependent */

   /* timezone = 10 * 60 * 60; *//* for Hawaii */
   /* daylight = 0;  */
   time(&clock_start);
   strncpy(time_stamp, ctime(&clock_start), 24);        /* stop at 24 to strip
                                                         * \n */
   check_error(get_parameters(fp_cnt, adcpsect_params, adcpsect_name_lists),
               "getting adcpsect parameters");
   if (combined == -1) { single_ens = 1; combined = 0;}
   else if (combined == FALSE) {separate = 1; }

   check_dbopen(ADCP_DB, DBNAME, READ_ONLY, DIR_IN_MEMORY);
   printf("\n\n DATABASE:  %s SUCCESSFULLY OPENED", DBNAME);

   new_extension(stat_filename, out_filename, ".sta");
   fp_out = check_fopen(stat_filename, "w");
   new_extension(mat_filename, out_filename, "_xy.mat");
   fp_zxyt_mat = check_fopen(mat_filename, "wb");
   new_extension(mat_filename, out_filename, "_uv.mat");
   fp_uv_mat = check_fopen(mat_filename, "wb");


   /*
    * ---> Read the list of options.
    */
   if (execute_options(fp_cnt, adcpsect_options, ECHO) <= 0)
   {
      printf("\n ERROR in reading adcpsect options.");
      exit(-1);
   }

   check_error((regrid_var != RV_DEPTH && regrid_var != RV_NONE)
               && !use_ctd,
               "regridding by svan or sigma theta requires ctd option");

   /*
    * ---> Calculate array sizes and allocate the arrays.
    */
   ndepth = min_val(ndepth, NDMAX - 1); /* ef */
   set_vgrid(ndepth);           /* sets adj_ndepth and adj_nu */
   set_ctd_grid(&ctd, regrid_var);
   allocate_unistat(&u_stat, adj_nu, "U", statistics);
   allocate_unistat(&v_stat, adj_nu, "V", statistics);
   allocate_unistat(&adcp_pos_stat, 3, "ADCP positions", U_MEAN | U_EXTREMA);
   allocate_unistat(&tr_temp_stat, 1, "TR-TEMP", U_MEAN | U_VARIANCE);
   allocate_unistat(&adcp_ship_stat, 2, "SHIP SPEED", U_MEAN | U_VARIANCE);
   if (transport)
      allocate_transport(adj_nu);

   nbin = max_val(3, adj_ndepth);
   temp_d_array = (double *) calloc(nbin, sizeof(double));
   check_error((temp_d_array == NULL), "insufficient memory for temp_d_array");

   /*
    * 91/11/18 jr++:  copy all non-double to double before saving to
    * mat-file:
    */
   for (i = 0; i < adj_ndepth; i++)
      temp_d_array[i] = (double) vgrid.array[i];
   savemat(fp_zxyt_mat, 'd', "z", adj_ndepth, 1, 0, (char *) temp_d_array, (char *) NULL);
   for (i = 0; i < adj_nu; i++)
      temp_d_array[i] = (double) vgrid.center_array[i];
   savemat(fp_zxyt_mat, 'd', "zc", adj_nu, 1, 0, (char *) temp_d_array, (char *) NULL);

   savemat(fp_zxyt_mat, 't', "z_label", 1, 1, 0,
           vgrid.label, (char *) NULL);
   savemat(fp_zxyt_mat, 't', "regrid_var", 1, 1, 0,
           get_name(regrid_variables, regrid_var), (char *) NULL);
   savemat(fp_zxyt_mat, 't', "regrid_method", 1, 1, 0,
           get_name(regrid_methods, regrid_method), (char *) NULL);
   savemat(fp_zxyt_mat, 't', "dbname", 1, 1, 0, DBNAME, (char *) NULL);
   savemat(fp_zxyt_mat, 't', "adcpsect_time", 1, 1, 0, time_stamp, (char *) NULL);
   temp_d_array[0] = (double) year_base;
   savemat(fp_zxyt_mat, 'd', "year_base", 1, 1, 0,
            (char *) temp_d_array, (char *) NULL);
   mat_zxyt_ptr = start_mat_array(fp_zxyt_mat, 'd', "xyt", 'c', 3);
   mat_uv_ptr = start_mat_array(fp_uv_mat, 'd', "uv", 'c', adj_nu);

   /*
    * ---> Start the loop through time ranges listed in the control file.
    */
   print_header(fp_out);
   if (combined) initialize_ensemble();
   while (get_time_range(fp_cnt, &start, &end_time) == 1)
   {
      check_dbsrch(TIME_SEARCH, (char *) &start);

      if (use_ctd)
      {
         check_dbset(CTD_DB, "DBSET, beginning of ctd block");

         check_dbsrch(TIME_SEARCH, (char *) &start);

         get_ctd_profile(&ctd);
         printf("\n ctd position: %f  %f",
                ctd.position.lat, ctd.position.lon);

         check_dbset(ADCP_DB, "DBSET, after get_ctd_profile");
      }

      if (combined && !contour)
         print_time_range(fp_out, &start, &end_time);
      print_time_range(stdout, &start, &end_time);
      if (separate)
         initialize_ensemble();
      /*
       * ---> Start the loop through profiles within the time range.
       */
      while (!IERR && check_time(&adcp.time, &start, &end_time))
      {
         if (verbose)
         {
            printf("\n");
            write_ymdhms_time(stdout, &adcp.time);
         }
         adcp.ddtime = year_day(&adcp.time, year_base);
         if (single_ens) initialize_ensemble();
         increment_ensemble();
         DBMOVE(&STEPS, &IERR); /* IERR is checked at the start of the loop */
         if (single_ens) end_ensemble();
      }                         /* end of loop through profiles within a time
                                 * range */
      if (separate)  end_ensemble();
   }                            /* end of loop through time ranges read from
                                 * control file */
   puts("\n End of control file\n");

   if (combined)  end_ensemble();
   time(&clock_end);
   run_time = (int) (clock_end - clock_start);
   fprintf(fp_out, "\n Elapsed time: %d:%02d (min:sec)\n",
           run_time / 60, run_time % 60);
   fclose(fp_out);
   if (contour)
      close_contour();
   if (vector)
      fclose(fp_vec);
   if (transport)
      close_transport();
   if (seq_cor)
      close_seqcor();
   if (nav.method)
      close_navigation();
   DBCLOSE(&IERR);
   check_error(IERR, "DBCLOSE");
   if (use_ctd)
   {
      check_dbset(CTD_DB, "DBSET, before closing ctd");
      DBCLOSE(&IERR);
      check_error(IERR, "DBCLOSE ctd");
   }
   free_unistat(&u_stat);
   free_unistat(&v_stat);
   free_unistat(&adcp_pos_stat);
   free_unistat(&adcp_ship_stat);
   free_unistat(&tr_temp_stat);
   close_vgrid();

   free(temp_d_array);
   end_mat_array(mat_zxyt_ptr);
   end_mat_array(mat_uv_ptr);
   fclose(fp_zxyt_mat);
   fclose(fp_uv_mat);

   printf("\nADCPSECT completed.\n");
   return;
}
