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

   FIX_TEMP.C

   Within a time range, this program corrects the transducer
   temperature (ancil1.tr_temp), which recorded incorrectly
   earler and those variables calculated from the transducer
   temperature including sound speed (ancil1.snd_spd_used),
   U, V, W and Error velocity.

   Set offset_temperature= 0 if input true_temperature will be
   used.  Otherwise

     true_temperature = original_temp + offset_temperature

   will be applied in velocity correction. Where original_temp
   is the temperature original recorded in the ping data file.

                       Willa Zhu
                       03-07-1990

Modified Mon  07-09-1990  by EF; new parameter original_soundspeed
is the soundspeed that was input to the DAS if the DAS was not set
to calculate soundspeed from temperature, or 0 otherwise.  In the
former case one can compensate for errors in manual input of
soundspeed to the DAS; in the latter one compensates for errors in the
soundspeed due to temperature errors.
------
Modified to add offset_temperature input parameter (WZ, 06/10/91)

92/01/15 (JR) Modified to allow user to provide temperature &
              salinity values from an external file with 2 columns:
              ADCP time (in decimal days) and the true temperature,
              a temperature offset, or salinity.
93/07/27 (JR) Moved n_vel initialization inside the while loop so
              all bins get grabbed and rotated even in the case
              where the conf.num_bins may change within the
              database
93/10/29 (JR) Added rescaling of bottom track velocities also.
              Previously, only U, V, W, and E were rescaled.
2000/06/23 (EF) Added ability to specify true_soundspeed= as a file.
2000/07/27 (EF) Changed MIN_SAL to 0, MAX_SAL to 40.

------------------------------------------------------------------

control file structure:

     dbname:    (e.g. WEP)
     original_soundspeed=
     offset_temperature=
     true_temperature=
     true_salinity=
     year_base=
     end
     time_ranges:
     (list of YMDHMS time pairs:)
        (yy/mm/dd hh:mm:ss) to (yy/mm/dd hh:mm:ss)
        .
        .
****************************************************************/

#include "geninc.h"
#include "ioserv.h"
#include "use_db.h"
#include "adcp.h"
#include "data_id.h"
#include "dpmask.h"

#define  sqr(x)  (x)*(x)
#define  equal(x,y)  (fabs((x)-(y)) < 0.0001)
#define  good_float(x)  ( (x) < ADJ_BADFLOAT )
#define  MAX_NBINS 128

#define OFFSET 1
#define TRUE_TEMP 2
#define TRUE_SALINITY 3
#define TRUE_SOUNDSPEED 4

#define MIN_TEMP -5.0 /* deg C */
#define MAX_TEMP 45.0
#define MIN_SAL  0.0 /* PSU */
#define MAX_SAL  40.0
#define MIN_SOUNDSPEED 1400 /* m/s */
#define MAX_SOUNDSPEED 1600

#define FV_UNUSED 0
#define FV_CONSTANT 1
#define FV_FILE 2

typedef struct
{
   int source;
   double value;
   FILE *fp;
} FILE_VALUE_TYPE;

#if PROTOTYPE_ALLOWED
int check_file(FILE *fp, RANGE_TYPE range, int year_base, int type);
int read_param(FILE *fp_cnt, char *arg, int code);
#endif

#if PROTOTYPE_ALLOWED
int read_param(FILE *fp_cnt, char *arg, int code)
#else
int read_param(fp_cnt, arg, code)
FILE *fp_cnt;
char *arg;
int code;
#endif
{
   FILE_VALUE_TYPE *fv;
   FILE_NAME_TYPE filename;
   double value;

   fv = (FILE_VALUE_TYPE *) arg;

   if (fscanf(fp_cnt, "%s", filename) != 1)
   {
      fprintf(stderr, " ERROR: Unable to read input parameter from control file\n\n");
      exit(-1);
   }

   if (sscanf(filename, " %lf", &value) != 1) /* user gave filename, not double value */
   {
      if ( (fv->fp = fopen(filename, "rt")) == NULL )
      {
         fprintf(stderr, " ERROR: Unable to open input file %s\n\n", filename);
         exit(-1);
      }
      if (ECHO) fprintf(stdout, " file: %s", filename);
      fv->source = FV_FILE;
   }
   else
   {
      fv->value = value;
      fv->source = FV_CONSTANT;
      if (ECHO) fprintf(stdout, " %f", fv->value);
      switch (code) /* check for reasonable values; in the case where values are
                       supplied in a file, the check is done later together with
                       the check for 1-1 match with ADCP times--see check_file() */
      {
         case TRUE_TEMP:
            if (fv->value != 0.0)  /* 0.0 allowed for backward compatibility */
            {
               if (fv->value < MIN_TEMP || fv->value > MAX_TEMP)
               {
                  fprintf(stderr, " ERROR: Unreasonable value for temperature %f\n\n",
                     fv->value);
                  exit(-1);
               }
            }
            break;
         case TRUE_SALINITY:
            if (fv->value < MIN_SAL || fv->value > MAX_SAL)
            {
                fprintf(stderr, " ERROR: Unreasonable value for salinity %f\n\n",
                   fv->value);
                exit(-1);
            }
            break;
         case OFFSET:
            /* too much hassle checking for reasonable offset:  would need to
               actually step thru every tr_temp in the database and see if result
               is reasonable */
            break;
         case TRUE_SOUNDSPEED:
            break;
         default:
            break;
      }
   }
   return(0);
}

/*
Main routine
*/

#if PROTOTYPE_ALLOWED
void do_it(FILE *fp_cnt)
#else
void do_it(fp_cnt)
FILE *fp_cnt;
#endif
{
   int              i, IERR, STEP=1;
   unsigned int     nvalues, n_vel, n_bt, ancil_size, n_time = sizeof(YMDHMS_TIME_TYPE);
   YMDHMS_TIME_TYPE prftime;
   YMDHMS_TIME_TYPE start, end,            /* for time range */
                    this_time;             /*  current profile time */
   ANCILLARY_1_TYPE ancil1;
   SHORT            u[MAX_NBINS], v[MAX_NBINS], w[MAX_NBINS], e[MAX_NBINS];
   BOTTOM_TRACK_TYPE bt;

   int const_snd_spd;  /* flag: 1 if sound speed was specified to the DAS,
                            0 if it was calculated by the DAS. */
   double old_snd_spd;
   double true_temp, corr_snd_spd, scale_factor;
   long int range_ofs;
   RANGE_TYPE range;
   char buff[80];
   double filetime, dbtime;
   int is_offset;

   static   FILE_NAME_TYPE  dbname;
   static   FILE_VALUE_TYPE temp, offset, sal, sndspd;
   static   double          orig_snd_spd;
   static   int             year_base;

   static   OPTION_TYPE options[] =
   {
      { "end"                 , 0            , NULL           , (char *) NULL},
      { "dbname:"             , TYPE_STRING  , op_get_param   , dbname       },
      { "year_base="          , TYPE_INT     , op_get_param   , (char *) &year_base},
      { "original_soundspeed=", TYPE_DOUBLE  , op_get_param   , (char *) &orig_snd_spd},
      { "offset_temperature=" , OFFSET       , read_param     , (char *) &offset},
      { "true_temperature="   , TRUE_TEMP    , read_param     , (char *) &temp},
      { "true_salinity="      , TRUE_SALINITY, read_param     , (char *) &sal},
      { "true_soundspeed="    , TRUE_SOUNDSPEED, read_param     , (char *) &sndspd},
      { NULL                  , 0            , NULL           , (char *) NULL}
   };                   /* options[] */


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

   sin(0.0);       /* to force linking of float routines */

   /* defaults: */
   orig_snd_spd = 0;
   year_base = -99;
   dbname[0] = '\0';
   is_offset = 0;  /* true temperature or offset? */
   temp.fp = sal.fp = offset.fp = sndspd.fp = NULL;
   temp.value = sal.value = offset.value = sndspd.value = 0.0;
   temp.source = sal.source = offset.source = sndspd.source = FV_UNUSED;

   fprintf(stdout, "\n Pass 1: Checking control file...\n");
   if (execute_options(fp_cnt, options, ECHO) <= 0)  exit(-1);
#if 0
   if (!( (temp.fp && temp.value == 0.0 && offset.fp == NULL && offset.value == 0.0) ||
          (offset.fp && offset.value == 0.0 && temp.fp == NULL && temp.value == 0.0) ||
          (temp.fp == NULL && temp.value != 0.0 && offset.fp == NULL && offset.value == 0.0) ||
          (offset.fp == NULL && offset.value != 0.0 && temp.fp == NULL && temp.value == 0.0) )
       || (sal.fp == NULL && sal.value == 0.0) )
#else
   if ( (sndspd.source && (temp.source || sal.source || offset.source))  ||
        (sndspd.source == UNUSED &&
            (sal.source == UNUSED ||
               (temp.source && offset.source) ||
                  (!temp.source && !offset.source))))

#endif
   {
      fprintf(stderr, " ERROR: You must supply either true_temperature OR offset_temperature\n");
      fprintf(stderr, "        AND true_salinity;  OR true_soundspeed alone\n\n");
      exit(-1);
   }
   is_offset = (offset.fp || offset.value != 0.0);

   if (orig_snd_spd < 1400.0)  const_snd_spd = 0;
   else                        const_snd_spd = 1;

   /*
   check_dbopen(1, dbname, READ_WRITE, DIR_IN_MEMORY);
   */
   check_dbopen(1, dbname, READ_WRITE, DIR_ON_DISK);
   if (year_base == -99) /* use first profile time as default */
   {
      if (dbget_cnf(TIME, (char *) &prftime, &n_time, "getting time")) exit(-1);
      year_base = prftime.year;
   }
   check_error( (fscanf(fp_cnt, " time_ranges:") != 0),
         "reading time range start");
   range.type = TIME_RANGE;
   range_ofs = ftell(fp_cnt); /* note down for processing pass later */
   while (get_time_range(fp_cnt, &(range.ru.time.start),
             &(range.ru.time.end)) == 1)
   {
      switch(check_range(&range))
      {
         case GAP:
            print_range(stdout, &range, "\n WARNING: Gap before range ");
            break;
         case OVERLAP:
            print_range(stdout, &range, "\n ERROR: Range overlaps preceding => ");
            exit(-1);
         case BAD_RANGE:
            print_range(stdout, &range, "\n ERROR: Bad range => ");
            exit(-1);
      }
      if (sndspd.fp && check_file(sndspd.fp, range, year_base, TRUE_SOUNDSPEED)) exit(-1);
      if (sal.fp && check_file(sal.fp, range, year_base, TRUE_SALINITY)) exit(-1);
      if (temp.fp && check_file(temp.fp, range, year_base, TRUE_TEMP)) exit(-1);
      else if (offset.fp && check_file(offset.fp, range, year_base, OFFSET)) exit(-1);
   }
   if (check_for_gap_at_end(&range) == GAP)
   {
      print_range(stdout, &range, "\n WARNING: Gap after range ");
   }

   /* rewind for pass 2: */
   fseek(fp_cnt, range_ofs, SEEK_SET);
   if (sal.fp) fseek(sal.fp, 0L, SEEK_SET);
   if (temp.fp) fseek(temp.fp, 0L, SEEK_SET);
   else if (offset.fp) fseek(offset.fp, 0L, SEEK_SET);
   if (sndspd.fp) fseek(sndspd.fp, 0L, SEEK_SET);
   fprintf(stdout, "\n Pass 2: Updating database...\n");
/*
---> Start the loop through time ranges listed in the control file.
*/
   while (get_time_range(fp_cnt, &start, &end) == 1)
   {
      check_dbsrch(TIME_SEARCH, (char *)&start);
      print_time_range(stdout, &start, &end);
      IERR = 0;
/*
---> Start the loop through profiles within the time range.
*/
      ancil_size = sizeof(ancil1);
      while (!IERR && check_time(&this_time, &start, &end))
      {
         dbtime = year_day(&this_time, year_base);
         n_vel = MAX_NBINS * sizeof(SHORT);
         n_bt = sizeof(BOTTOM_TRACK_TYPE);

         if (sndspd.source == FV_FILE)
            while (getline_nc(sndspd.fp, buff, 80) != EOF &&
               sscanf(buff, " %lf %lf", &filetime, &sndspd.value) == 2 &&
               !equal(filetime, dbtime));

         if (sal.source == FV_FILE)
            while (getline_nc(sal.fp, buff, 80) != EOF &&
               sscanf(buff, " %lf %lf", &filetime, &sal.value) == 2 &&
               !equal(filetime, dbtime));
         if (temp.source == FV_FILE)
            while (getline_nc(temp.fp, buff, 80) != EOF &&
               sscanf(buff, " %lf %lf", &filetime, &temp.value) == 2 &&
               !equal(filetime, dbtime));
         else if (offset.source == FV_FILE)
            while (getline_nc(offset.fp, buff, 80) != EOF &&
               sscanf(buff, " %lf %lf", &filetime, &offset.value) == 2 &&
               !equal(filetime, dbtime));

         if ((IERR = dbget_cnf(ANCILLARY_1, (char *)&ancil1, &ancil_size,
            "getting ANCILARY_1")) == 0 &&
            (IERR = dbget_cnf(U, (char *)u, &n_vel,"getting U")) == 0 &&
            (IERR = dbget_cnf(V, (char *)v, &n_vel,"getting V")) == 0 &&
            (IERR = dbget_cnf(W, (char *)w, &n_vel,"getting W")) == 0 &&
            (IERR = dbget_cnf(ERROR_VEL,(char *)e, &n_vel, "getting ERROR_VEL")) == 0 &&
            (IERR = dbget_cnf(BOTTOM_TRACK, (char *)&bt, &n_bt,"getting BOTTOM TRACK")) == 0)
         {
           
           /* 2000/12/12 PJ
            * Added a message in case n_vel becomes zero after one of the
            * calls above. This is unfortunately the case sometimes.
            */
           if(n_vel == 0){
             char str[299];
             
             sprintf(str, 
                     "\nBad number of elements in profile at time %12.6f\n",
                     dbtime);
             report_msg(str);
           }

            if (offset.source)
               true_temp = ancil1.tr_temp + offset.value;
            else
               true_temp = temp.value;

            if (sndspd.source)
            {
               corr_snd_spd = sndspd.value;
            }
            else
            {
               corr_snd_spd =   1449.2 + 4.6 * true_temp
                        - 0.055 * true_temp*true_temp
                        + 0.00029 * (true_temp*true_temp*true_temp)
                        + (1.34 - 0.01*true_temp) * (sal.value-35.0);
            }
            /* one more check for reasonable values: */
            if (corr_snd_spd < MIN_SOUNDSPEED || corr_snd_spd > MAX_SOUNDSPEED)
            {
               fprintf(stderr, " ERROR: Unreasonable value (%f) for soundspeed for dbtime %f\n\n",
                  corr_snd_spd, dbtime);
               goto close;
            }
            if (const_snd_spd)
               old_snd_spd = orig_snd_spd;
            else
               old_snd_spd = ancil1.snd_spd_used;

            nvalues = n_vel / sizeof(SHORT);
            scale_factor = corr_snd_spd/old_snd_spd;

            for (i = 0; i < nvalues; i++)
            {
               if (u[i] != BADSHORT)
                  u[i] = (double)u[i] * scale_factor;
               if (v[i] != BADSHORT)
                  v[i] = (double)v[i] * scale_factor;
               if (w[i] != BADSHORT)
                  w[i] = (double)w[i] * scale_factor;
               if (e[i] != BADSHORT)
                  e[i] = (double)e[i] * scale_factor;
            }
            if (bt.u < ADJ_BADFLOAT)
            {
               bt.u = (double) bt.u * scale_factor;
               bt.v = (double) bt.v * scale_factor;
            }
            ancil1.snd_spd_used = corr_snd_spd;
            ancil1.tr_temp = true_temp;

            ancil_size = sizeof(ancil1);
            if((IERR = dbput_cnf(ANCILLARY_1, (char *) &ancil1,
               &ancil_size, "putting ANCILLARY_1")) == 0)
               if ((IERR = dbput_cnf(U, (char *)u, &n_vel,
                  "putting U velocity")) == 0)
                  if ((IERR = dbput_cnf(V, (char *)v, &n_vel,
                       "putting V velocity")) == 0)
                     if ((IERR = dbput_cnf(W, (char *)w, &n_vel,
                          "putting W velocity")) == 0)
                        if((IERR = dbput_cnf(ERROR_VEL, (char *) e, &n_vel,
                                         "putting ERROR_VEL")) == 0)
                          if((IERR = dbput_cnf(BOTTOM_TRACK, (char *) &bt, &n_bt,
                                         "putting BOTTOM_TRACK")) == 0)
                           if ( (IERR = set_block_dpmask(SOUND_SPEED_VEL_COR)) == 0)
                              DBMOVE(&STEP, &IERR);
         }
      }  /* end of loop through profiles within a time range */
   }  /* end of loop through time ranges read from control file */

close:
   DBCLOSE(&IERR);
   check_error(IERR, "DBCLOSE");
   printf("\n fix_temp completed.\n");
   return;
}

#if PROTOTYPE_ALLOWED
int check_file(FILE *fp, RANGE_TYPE range, int year_base, int type)
#else
int check_file(fp, range, year_base, type)
FILE *fp;
RANGE_TYPE range;
int year_base;
int type;   /* TRUE_SOUNDSPEED, TRUE_TEMP, OFFSET or TRUE_SALINITY */
#endif
{
   YMDHMS_TIME_TYPE current;
   char             buff[80];
   double           dbtime, filetime, value;
   ANCILLARY_1_TYPE anc1;

   if (goto_start_of_range(&range)) return(1);
   while (in_range((char *) &current, &range))
   {
      dbtime = year_day(&current, year_base);
      while (getline_nc(fp, buff, 80) != EOF &&
         sscanf(buff, " %lf %lf", &filetime, &value) == 2 &&
         !equal(filetime, dbtime)); /* match times */
      if (feof(fp))
      {
         fprintf(stderr, " ERROR: No file time to match database time %f\n\n", dbtime);
         return(-1);
      }
      /* check for reasonable value */
      switch (type)
      {
         case OFFSET:
            check_dbget(ANCILLARY_1, (char *) &anc1, sizeof(ANCILLARY_1_TYPE), "getting ANCILLARY_1");
            value += anc1.tr_temp;
         case TRUE_TEMP:
            if (value < MIN_TEMP || value > MAX_TEMP)
            {
               fprintf(stderr, " ERROR: Unreasonable temperature for time %f in file\n\n", filetime);
               exit(-1);
            }
            break;
         case TRUE_SALINITY:
            if (value < MIN_SAL || value > MAX_SAL)
            {
               fprintf(stderr, " ERROR: Unreasonable salinity for time %f in file\n\n", filetime);
               exit(-1);
            }
            break;
         case TRUE_SOUNDSPEED:
            if (value < MIN_SOUNDSPEED || value > MAX_SOUNDSPEED)
            {
               fprintf(stderr, " ERROR: Unreasonable value (%f) for soundspeed at time  %f in file\n\n",
                  value, filetime);
               exit(-1);
            }
            break;
         default:
            break;
      }
      switch(dbmove_cnf(1))
      {
         case 0:
            break;
         case SEARCH_BEYOND_END:
            return(0);
         default:
            return(-1);
      }
   }
   return(0);
}

