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

   FLUX.C module for ADCPSECT

   This module integrates in the horizontal, from ensemble
   to ensemble, using trapezoidal integration.  If the input
   has already been integrated in the vertical by selecting
   the "integrate" regridding option, then the result is the
   transport.

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

   Tue  03-05-91 Added a transport direction option "along"
                which results in the line integral along the ship track
                being calculated: udx+vdy.

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


#include <stdio.h>
#include <string.h>
#include <math.h>
#include "dbext.h"
#include "ioserv.h"
#include "vector.h"
#include "use_db.h"
#include "vstat.h"
#include "pos_to_m.h"
#include "adcp.h"
#include "adcpsect.h"

#define DIR_ALONG 2
#define DIR_ALL   0

extern  YMDHMS_TIME_TYPE start, end_time;    /* for time range */


int transport = FALSE;             /* global flag */
extern char out_filename[];
extern char time_stamp[];
static char tran_filename[NPATHMAX];
static FILE *fp_tran;
static FILE *fp_tran_mat;
static MAT_ARRAY_FP_TYPE *mat_ptr;

/*--> for transport calculation using trapezoidal rule */
static double u_prev[NDMAX], v_prev[NDMAX], tran_inc[NDMAX];
static double u_this[NDMAX], v_this[NDMAX];
static double u_nan[NDMAX]; /* temporary array with BADFLOAT converted to NaN */
static DOUBLE_LL_TYPE pos_prev;
static UNISTAT_TYPE    u_tran, v_tran;
static double tran_units = 1.0;
static int tran_first_ensemble = 1;  /* flag to suppress integration */
static int tran_min_npts;         /* minimum number for transport output */
                                  /* This threshold will determine whether
                                     integration proceeds at a given bin.
                                     If it stops, the sum is not zeroed or
                                     changed to bad; rather, the number of
                                     points in the integral is reduced.
                                     Hence, the user must look at N in the
                                     printout to determine whether all profiles
                                     were included in the integration at a
                                     given depth.
                                   */
static int tran_dir;                 /* direction of transport rel to track */

NAME_LIST_ENTRY_TYPE tran_directions[] =
{
   {"right"   ,    1},
   {"left"    ,   -1},
   {"all"     ,    DIR_ALL},
   {"along"   ,    DIR_ALONG},
   {NULL      ,    0}
};

NAME_LIST_TYPE transport_name_lists[] =
{
   tran_directions
};

PARAMETER_LIST_ENTRY_TYPE transport_params[] =
{
   {" minimum_npts=" , " %d"     , &tran_min_npts  , TYPE_INT},
   {" units="        , " %lf"    , &tran_units     , TYPE_DOUBLE},
   {" direction="    , " %s"     , &tran_dir       , TYPE_TRANS},
   {NULL             , NULL      , NULL            , 0}
};

#if PROTOTYPE_ALLOWED
int init_transport(FILE *fp_cnt)
#else
int init_transport(fp_cnt)
FILE *fp_cnt;
#endif
{
   check_error(get_parameters(fp_cnt, transport_params, transport_name_lists),
                  "transport parameters");
   transport = TRUE;
   /* Make a transport output file name, and open the file. */
   new_extension(tran_filename, out_filename, ".trn");
   fp_tran = check_fopen(tran_filename, "w");
   new_extension(tran_filename, out_filename, "_tr.mat");
   fp_tran_mat = check_fopen(tran_filename, "wb");
   savemat(fp_tran_mat, 't', "direction", 1, 1, 0,
                        get_name(tran_directions, tran_dir), (char *)NULL);
   savemat(fp_tran_mat, 't', "adcpsect_flux_time", 1, 1, 0, time_stamp,
                           (char *)NULL);
   return(0);
}                       /* init_transport */

/* The allocation of arrays must be done in a separate function
   because it must wait until adj_ndepth is known.
*/
#if PROTOTYPE_ALLOWED
void allocate_transport(int adj_nu)
#else
void allocate_transport(adj_nu)
int adj_nu;
#endif
{
   allocate_unistat(&u_tran, adj_nu, "U transport", U_MEAN);
   allocate_unistat(&v_tran, adj_nu, "V transport", U_MEAN);
   /* Transport calculation will be cumulative for the
      whole run of adcpsect, so this is the only place
      that zeroing is needed.  If malloc zeros the new
      memory, it would not really be needed even here.
   */
   zero_unistat(&u_tran);
   zero_unistat(&v_tran);
   mat_ptr = start_mat_array(fp_tran_mat, 'd', "transport", 'c', adj_nu);
}                       /* allocate_transport */

#if PROTOTYPE_ALLOWED
void integrate_ds(double *last_ptr, double *now_ptr, int adj_nu, 
		  double distance, double *inc_ptr, int sign)
#else
void integrate_ds(last_ptr, now_ptr, adj_nu, distance, inc_ptr, sign)
double *last_ptr, *now_ptr;
int adj_nu, sign;
double distance, *inc_ptr;
#endif
{
   double d, sum;
   int i;
   int last_pos, now_pos;

   /* On entry, sign is the desired sign of the product, u*ds.
      By reversing sign if ds < 0, we make sign into the
      desired sign of u alone.
   */
   if (distance < 0.0)
      sign *= -1;

   distance *= 0.5;   /* for trapezoidal rule */
   for (i=0; i<adj_nu; i++, now_ptr++, last_ptr++, inc_ptr++)
   {
      if (good_float(*now_ptr) && good_float(*last_ptr))
      {
         if (sign != 0)
         {

            last_pos = (*last_ptr > 0.0);
            now_pos  = (*now_ptr  > 0.0);
            sum = fabs(*last_ptr) + fabs(*now_ptr);

            if (   (sign > 0 && !last_pos && !now_pos)
                || (sign < 0 && last_pos && now_pos) )
               *inc_ptr = 0.0;  /* not in an integration region */


            else if ( last_pos && !now_pos)
            {
               /* transition plus to minus */
               if (sign > 0)  /* ending an integration region */
               {
                  d = distance * fabs(*last_ptr) / sum;
                  *inc_ptr = d * (*last_ptr);
               }
               else           /* starting an integration region */
               {
                  d = distance * fabs(*now_ptr) / sum;
                  *inc_ptr = d * (*now_ptr);
               }
            }
            else if ( !last_pos && now_pos)
            {
               /* transition minus to plus */
               if (sign < 0)  /* ending an integration region */
               {
                  d = distance * fabs(*last_ptr) / sum;
                  *inc_ptr = d * (*last_ptr);
               }
               else           /* starting an integration region */
               {
                  d = distance * fabs(*now_ptr) / sum;
                  *inc_ptr = d * (*now_ptr);
               }
            }
            else              /* entirely in an integration region */
            {
               d = distance;
               *inc_ptr = d * (*last_ptr + *now_ptr);
            }
         }
         else  /* sign == 0 */
         {
            d = distance;
            *inc_ptr = d * (*last_ptr + *now_ptr);
         }

      }
      else
         *inc_ptr = BADFLOAT;
   }
}                       /* integrate_ds() */

#if PROTOTYPE_ALLOWED
void calculate_transport(UNISTAT_TYPE *u_stat, UNISTAT_TYPE *v_stat, UNISTAT_TYPE *adcp_pos_stat)
#else
void calculate_transport(u_stat, v_stat, adcp_pos_stat)
UNISTAT_TYPE *u_stat, *v_stat, *adcp_pos_stat;
#endif
{
   double dx, dy;
   int i;

   memcpy(u_this, u_stat->sum, u_stat->n * sizeof(double));
   memcpy(v_this, v_stat->sum, v_stat->n * sizeof(double));

   if (!tran_first_ensemble)
   {
      dx = dlon_to_meters(adcp_pos_stat->sum[LONGITUDE] - pos_prev.lon,
                              pos_prev.lat);
      dy = dlat_to_meters(adcp_pos_stat->sum[LATITUDE]  - pos_prev.lat,
                              pos_prev.lat);

      /* To block the integration if there are not enough points in
         the present ensemble, mark u and v as bad.
      */
      for (i=0; i<u_stat->n; i++)
      {
         if (u_stat->npts[i] < tran_min_npts)
         {
            u_this[i] = BADFLOAT;
            v_this[i] = BADFLOAT;
         }
      }

      if (tran_dir == DIR_ALONG)
         integrate_ds(u_prev, u_this, u_stat->n, dx, tran_inc, DIR_ALL);
      else
         integrate_ds(u_prev, u_this, u_stat->n, dy, tran_inc, tran_dir);
      update_unistat_d(&u_tran, tran_inc, u_stat->n);

      /* Since we are passing -dx, not dx, into integrate_ds
         for the v component, we pass tran_dir, not -tran_dir.
      */
      if (tran_dir == DIR_ALONG)
         integrate_ds(v_prev, v_this, u_stat->n, dy, tran_inc, DIR_ALL);
      else
         integrate_ds(v_prev, v_this, u_stat->n, -dx, tran_inc, tran_dir);
      update_unistat_d(&v_tran, tran_inc, u_stat->n);
   }
   pos_prev.lon = adcp_pos_stat->sum[LONGITUDE];
   pos_prev.lat = adcp_pos_stat->sum[LATITUDE];
   memcpy(u_prev, u_this, u_stat->n * sizeof(double));
   memcpy(v_prev, v_this, u_stat->n * sizeof(double));

   tran_first_ensemble = 0;
}                       /* calculate_transport() */



/***************************************************************
      write_transport()

      (derived from write_unistat)

      The input variable d is a pointer to a float array
      used for labelling the rows of output.  If the
      statistics are calculated from velocity profiles, then
      d is the array of depths.  If fourier transformed
      profiles were used, then d is the array of vertical
      wavenumbers.  If d is a NULL pointer, then the array
      index, starting from 1, will be printed.

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

#define depth  (d==NULL ? (float) i+2 : d[i+1])

#if PROTOTYPE_ALLOWED
void write_transport(FILE *fpt, UNISTAT_TYPE *u, UNISTAT_TYPE *v, 
		     double scale, float *d, char *d_name, int nd_prec)
#else
void write_transport(fpt, u, v, scale, d, d_name, nd_prec)
   FILE *fpt;
   UNISTAT_TYPE *u, *v;
   double scale;  /* e.g., set scale to 0.01 to convert m/s to cm/s */
   float *d;       /* if it is a null pointer, then the index will be used */
   char *d_name;   /* for labelling the "depth" column */
   int  nd_prec;   /* precision (decimal places) for "depth" column */
#endif
{
   int i, np;
   double uu, vv;
   double uus = 0.0, vvs = 0.0;

   check_error( (scale == 0.0), "write_transport; scale is zero");
   fprintf(fpt, "\nADCPSECT %s", time_stamp);
   fprintf(fpt, "\nTransports scaled by %e", scale);
   fprintf(fpt, "\nDepth, bin, or density variable gives lower boundary of integration");
   fprintf(fpt, "\n                     --------horizontal--------");
   fprintf(fpt,                    "   --horizontal and vertical--");
   if (tran_dir == DIR_ALONG)
      fprintf(fpt, "\n%12s    N       udx       vdy       sum       udx       vdy       sum ", d_name);
   else
      fprintf(fpt, "\n%12s    N       udy      -vdx       sum       udy      -vdx       sum ", d_name);
   for (i=0; i<u->n; i++)  if (depth <= ADJ_BADFLOAT)
   {
      if ( (np = (u->npts)[i]) >0 )
      {
         uu = (u->sum)[i]/scale;
         vv = (v->sum)[i]/scale;
         uus += uu;
         vvs += vv;
         fprintf(fpt, "\n%10.*f %6d %9.3f %9.3f %9.3f %9.3f %9.3f %9.3f", nd_prec,
                       depth, np, uu, vv, uu+vv, uus, vvs, uus+vvs);
      }
   }
   fprintf(fpt, "\n");
}                       /* write_transport() */
#undef depth

#if PROTOTYPE_ALLOWED
void write_tran_mat(MAT_ARRAY_FP_TYPE *mat_ptr, double *u, double *v, int n)
#else
void write_tran_mat(mat_ptr, u, v, n)
   MAT_ARRAY_FP_TYPE *mat_ptr;
   double *u, *v;
   int n;
#endif
{
   tonan_d(u_nan, u, n);
   add_to_mat_array(mat_ptr, (char *)u_nan);
   tonan_d(u_nan, v, n);
   add_to_mat_array(mat_ptr, (char *)u_nan);
}

#if PROTOTYPE_ALLOWED
void print_transport(UNISTAT_TYPE *adcp_pos_stat, float *y, char *label, int n_zdigits)
#else
void print_transport(adcp_pos_stat, y, label, n_zdigits)
   UNISTAT_TYPE *adcp_pos_stat;
   float *y;
   char *label;
   int n_zdigits;
#endif
{
   extern YMDHMS_TIME_TYPE start, end_time;
   write_vgrid_params(fp_tran);
   print_time_range(fp_tran, &start, &end_time);
   fprintf(fp_tran, "\n ADCP mean position: %5.2f lat, %6.2f long\n",
                    adcp_pos_stat->sum[LATITUDE],
                    adcp_pos_stat->sum[LONGITUDE]);

   write_transport(fp_tran, &u_tran, &v_tran, tran_units, y, label, n_zdigits);
   write_tran_mat(mat_ptr, u_tran.sum, v_tran.sum, u_tran.n);

   fprintf(fp_tran, "\n\n");
}

#if PROTOTYPE_ALLOWED
void close_transport(void)
#else
void close_transport()
#endif
{
   fclose(fp_tran);
   end_mat_array(mat_ptr);
   fclose(fp_tran_mat);
   free_unistat(&u_tran);
   free_unistat(&v_tran);
}

