// $Id: ConstantSetWrapper.cc,v 1.3 2002/12/19 18:40:45 flaterco Exp $
/*  ConstantSetWrapper

    Copyright (C) 1998  David Flater.

    This program is free software; you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation; either version 2 of the License, or
    (at your option) any later version.

    This program is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with this program; if not, write to the Free Software
    Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
*/

#include "common.hh"

Amplitude
ConstantSetWrapper::dt_tide_max (unsigned deriv) {
  /* We need to be able to calculate max tide derivatives for one
   * derivative higher than we actually need to know the tides.
   */
  assert(deriv <= TIDE_MAX_DERIV+1);
  // This is initialized in the constructor.
  return maxdt[deriv];
}

ConstantSetWrapper::ConstantSetWrapper (ConstituentSet *in_constituents,
                      ConstantSet *in_constants) {
  unsigned i, tempyear;

  assert (in_constituents->length == in_constants->length);
  length = in_constants->length;
  constituents = in_constituents;
  origConstants = in_constants;

  {
    // Squeeze out null constituents
    unsigned newlength = 0;
    for (i=0; i<length; i++) {
      if (origConstants->amplitudes[i].val() > 0.0) {
        if (i > newlength) {
          origConstants->amplitudes[newlength] = origConstants->amplitudes[i];
          origConstants->phases[newlength] = origConstants->phases[i];
          constituents->constituents[newlength] = constituents->constituents[i];
        }
        newlength++;
      }
    }
    assert (newlength > 0);
    length = origConstants->length = constituents->length = newlength;
  }

  speeds = new double[length];
  amplitudes = new double[length];
  phases = new double[length];

  // Prefetch speeds
  for (i=0; i<length; i++)
    speeds[i] = (*in_constituents)[i].speed().rad(Speed::SECOND);

  // Nasty loop to figure orig_maxdt and origMaxAmplitude
  unsigned deriv;
  for (deriv=0; deriv<=TIDE_MAX_DERIV+1; deriv++) {
    for (tempyear=(*constituents)[0].firstvalidyear().val();
    tempyear<=(*constituents)[0].lastvalidyear().val(); tempyear++) {
      Year ty (tempyear);
      Amplitude max;
      for (i=0;i<length;i++)
        max += origConstants->amplitudes[i] * (*constituents)[i].nod(ty)
               * pow(speeds[i], (double)deriv);
      if (orig_maxdt[deriv].val() == 0.0)
        orig_maxdt[deriv] = max; // Get units
      else if (max > orig_maxdt[deriv])
        orig_maxdt[deriv] = max;
    }
    if (deriv == 0)
      origMaxAmplitude = orig_maxdt[deriv];
    orig_maxdt[deriv] *= 1.1;      /* Add a little safety margin... */
  }
  assert (origMaxAmplitude.val() > 0.0);
  if (origMaxAmplitude.Units().mytype == PredictionValue::Unit::KnotsSquared)
    origMaxAmplitude.Units (PredictionValue::Unit::Knots);

#ifdef SUPER_ULTRA_VERBOSE_DEBUGGING
  cerr << "Calculated origMaxAmplitude: " << origMaxAmplitude << endl;
  cerr << "Calculated orig_maxdt: [";
  for (deriv=0; deriv<=TIDE_MAX_DERIV+1; deriv++) {
    if (deriv)
      cerr << ",";
    cerr << orig_maxdt[deriv];
  }
  cerr << "]" << endl;
#endif

  // Initialize to valid values.
  myUnits = in_constants->datum.Units();
  mySimpleOffsets.levelAdd.Units (myUnits);
  maxAmplitude = origMaxAmplitude;

  // Harmonics file range of years may exceed that of this platform.
  // Try valiantly to find a safe initial value.
  {
    unsigned b = (*in_constituents)[0].firstvalidyear().val();
    unsigned e = (*in_constituents)[0].lastvalidyear().val();
    if (b <= 2000 && e >= 2000)
      currentYear = Year(2000);
    else if (b <= 1970 && e >= 1970)
      currentYear = Year(1970);
    else if (b <= 2037 && e >= 2037)
      currentYear = Year(2037);
    else
      currentYear = Year((b+e)/2);
  }

  refresh_adjConstants_yearConstants ();
}

void ConstantSetWrapper::refresh_adjConstants_yearConstants () {
  adjConstants = *origConstants;

  // FIXME if there are ever units of velocity other than knots.
  // (This was just a kludge to avoid ruining knots squared.)
  // Two more fixmes below...
  if (adjConstants.datum.Units().mytype != PredictionValue::Unit::Knots &&
    adjConstants.datum.Units() != myUnits)
    adjConstants.setUnits (myUnits);

  adjConstants.adjust (mySimpleOffsets, *constituents);
  maxAmplitude = origMaxAmplitude * mySimpleOffsets.levelMultiply();
  // FIXME if there are ever units of velocity other than knots.
  // maxAmplitude will never be KnotsSquared
  if (maxAmplitude.Units().mytype != PredictionValue::Unit::Knots &&
      maxAmplitude.Units() != myUnits)
    maxAmplitude.Units (myUnits);
  // Update maxdt the same way
  for (unsigned deriv=0; deriv<=TIDE_MAX_DERIV+1; deriv++) {
    maxdt[deriv] = orig_maxdt[deriv] * mySimpleOffsets.levelMultiply();
    // FIXME if there are ever units of velocity other than knots.
    if (maxdt[deriv].Units().mytype != PredictionValue::Unit::Knots &&
        maxdt[deriv].Units().mytype != PredictionValue::Unit::KnotsSquared &&
        maxdt[deriv].Units() != myUnits)
      maxdt[deriv].Units (myUnits);
  }
  refresh_yearConstants();
}

void ConstantSetWrapper::refresh_yearConstants () {
  yearConstants = adjConstants;
  for (unsigned i=0; i<length; i++) {
    // Apply node factor
    yearConstants.amplitudes[i] *= (*constituents)[i].nod(currentYear);
    amplitudes[i] = yearConstants.amplitudes[i].val();
    // Apply equilibrium argument.  Recall that phases have been pre-negated
    // per -k'.
    yearConstants.phases[i] += (*constituents)[i].arg(currentYear);
    phases[i] = yearConstants.phases[i].rad();
  }
  epoch = Timestamp (currentYear);
  next_epoch = Timestamp (currentYear + 1);
}

ConstantSetWrapper::~ConstantSetWrapper () {
  delete constituents;
  delete origConstants;
  delete [] speeds;
  delete [] amplitudes;
  delete [] phases;
}

void ConstantSetWrapper::setSimpleOffsets (SimpleOffsets in_offsets) {
  mySimpleOffsets = in_offsets;
  refresh_adjConstants_yearConstants ();
}

void ConstantSetWrapper::setUnits (PredictionValue::Unit in_units) {
  // PredictionValue::setUnits enforces possible conversions.
  myUnits = in_units;
  refresh_adjConstants_yearConstants ();
}

PredictionValue::Unit ConstantSetWrapper::predictUnits () {
  // Kludgy?
  return yearConstants.amplitudes[0].Units();
}

PredictionValue ConstantSetWrapper::predictHeight (Timestamp in_timestamp,
unsigned deriv) {
  // We need to take this double and make it civilized
  PredictionValue a (predictUnits(), time2dt_tide (in_timestamp, deriv));

  // Don't do this here.
  // if (a.Units().mytype == PredictionValue::Unit::KnotsSquared)
  //   a.Units (PredictionValue::Unit::Knots);

  return a;
}

PredictionValue ConstantSetWrapper::movingMean (Timestamp in_timestamp) {
  PredictionValue a (predictUnits(), time2mean (in_timestamp));
  return a;
}

PredictionValue ConstantSetWrapper::datum() const {
  return yearConstants.datum;
}

// The following block of functions is slightly revised from the code
// delivered by Geoffrey T. Dairiki for XTide 1.  I have refrained
// from renaming the functions (much) so that Jeff's comments might
// still make (some) sense.

/*************************************************************************
 *
 * Geoffrey T. Dairiki Fri Jul 19 15:44:21 PDT 1996
 *
 ************************************************************************/

/*
 * We will need a function for tidal height as a function of time
 * which is continuous (and has continuous first and second derivatives)
 * for all times.
 *
 * Since the epochs & multipliers for the tidal constituents change
 * with the year, the regular time2tide(t) function has small
 * discontinuities at new years.  These discontinuities really
 * fry the fast root-finders.
 *
 * We will eliminate the new-years discontinuities by smoothly
 * interpolating (or "blending") between the tides calculated with one
 * year's coefficients, and the tides calculated with the next year's
 * coefficients.
 *
 * i.e. for times near a new years, we will "blend" a tide
 * as follows:
 *
 * tide(t) = tide(year-1, t)
 *                  + w((t - t0) / Tblend) * (tide(year,t) - tide(year-1,t))
 *
 * Here:  t0 is the time of the nearest new-year.
 *        tide(year-1, t) is the tide calculated using the coefficients
 *           for the year just preceding t0.
 *        tide(year, t) is the tide calculated using the coefficients
 *           for the year which starts at t0.
 *        Tblend is the "blending" time scale.  This is set by
 *           the macro TIDE_BLEND_TIME, currently one hour.
 *        w(x) is the "blending function", whice varies smoothly
 *           from 0, for x < -1 to 1 for x > 1.
 *
 * Derivatives of the blended tide can be evaluated in terms of derivatives
 * of w(x), tide(year-1, t), and tide(year, t).  The blended tide is
 * guaranteed to have as many continuous derivatives as w(x).  */

/* time2dt_tide(time_t t, unsigned n)
 *
 *   Calculate nth time derivative the normalized tide.
 *
 * Notes: This function does not check for changes in year.
 *  This is important to our algorithm, since for times near
 *  new years, we interpolate between the tides calculated
 *  using one years coefficients, and the next years coefficients.
 *
 *  Except for this detail, time2dt_tide(t,0) should return a value
 *  identical to time2tide(t).
 */

// t has been changed to be seconds since the epoch.
// (This is especially good since time_t might be a long long.)

double
ConstantSetWrapper::_time2dt_tide (long t, unsigned deriv) {
  double dt_tide = 0.0;
  unsigned a;
  int b;
  double term;

  double tempd = M_PI / 2.0 * deriv;
  for (a=0; a<length; a++) {
    term = amplitudes[a] * cos (tempd + speeds[a] * t + phases[a]);
    for (b = deriv; b > 0; b--)
      term *= speeds[a];
    dt_tide += term;
  }
  return dt_tide;
}

/* blend_weight (double x, unsigned deriv)
 *
 * Returns the value nth derivative of the "blending function" w(x):
 *
 *   w(x) =  0,     for x <= -1
 *
 *   w(x) =  1/2 + (15/16) x - (5/8) x^3 + (3/16) x^5,
 *                  for  -1 < x < 1
 *
 *   w(x) =  1,     for x >= 1
 *
 * This function has the following desirable properties:
 *
 *    w(x) is exactly either 0 or 1 for |x| > 1
 *
 *    w(x), as well as its first two derivatives are continuous for all x.
 */
double
ConstantSetWrapper::blend_weight (double x, unsigned deriv)
{
  double x2 = x * x;

  if (x2 >= 1.0)
      return deriv == 0 && x > 0.0 ? 1.0 : 0.0;

  switch (deriv) {
  case 0:
      return ((3.0 * x2 -10.0) * x2 + 15.0) * x / 16.0 + 0.5;
  case 1:
      return ((x2 - 2.0) * x2 + 1.0) * (15.0/16.0);
  case 2:
      return (x2 - 1.0) * x * (15.0/4.0);
  }
  assert(0);
  // Silence bogus SGI compiler warning.
  return 0.0;
}

/*
 * This function does the actual "blending" of the tide
 * and its derivatives.
 */
double
ConstantSetWrapper::blend_tide (Timestamp in_timestamp, unsigned deriv,
Year first_year, double blend) {
  double        fl[TIDE_MAX_DERIV + 1];
  double        fr[TIDE_MAX_DERIV + 1];
  double *      fp      = fl;
  double        w[TIDE_MAX_DERIV + 1];
  double        fact = 1.0;
  double        f;
  unsigned      n;
  long t;

  assert (deriv <= TIDE_MAX_DERIV);

  /*
   * If we are already set up for one of the two years
   * of interest, compute that year's tide values first.
   */
  if (currentYear == first_year + 1)
    fp = fr;
  else if (currentYear != first_year) {
    currentYear = first_year;
    refresh_yearConstants ();
  }

  // This amount is less than a year, so a long int suffices.
  t = (in_timestamp - epoch).in_seconds();
  for (n = 0; n <= deriv; n++)
    fp[n] = _time2dt_tide(t, n);

  /*
   * Compute tide values for the other year of interest,
   *  and the needed values of w(x) and its derivatives.
   */
  if (fp == fl) {
    currentYear = first_year + 1;
    fp = fr;
  } else {
    currentYear = first_year;
    fp = fl;
  }
  refresh_yearConstants ();
  t = (in_timestamp - epoch).in_seconds();
  for (n = 0; n <= deriv; n++) {
    fp[n] = _time2dt_tide(t, n);
    w[n] = blend_weight(blend, n);
  }

  /*
   * Do the blending.
   */
  f = fl[deriv];
  for (n = 0; n <= deriv; n++) {
    f += fact * w[n] * (fr[deriv-n] - fl[deriv-n]);
    fact *= (double)(deriv - n)/(n+1) * (1.0/TIDE_BLEND_SECONDS);
  }
  return f;
}

double
ConstantSetWrapper::time2dt_tide (Timestamp in_timestamp, unsigned deriv)
{
  // For starters, get us in the right year.
  Year in_year = in_timestamp.year();
  if (in_year != currentYear) {
    assert (!(currentYear.isNull()));
    currentYear = in_year;
    refresh_yearConstants ();
  }

  // This amount is less than a year, so a long int suffices.
  long t = (in_timestamp - epoch).in_seconds();

  /*
   * If we're close to either the previous or the next
   * new years we must blend the two years' tides.
   */
  if (t <= TIDE_BLEND_SECONDS) {
    return blend_tide(in_timestamp, deriv, currentYear - 1,
                     ((double)t)/TIDE_BLEND_SECONDS);
  } else {
    if (!(next_epoch.isNull())) {
      long u = (next_epoch - in_timestamp).in_seconds();
      if (u <= TIDE_BLEND_SECONDS)
        return blend_tide(in_timestamp, deriv, currentYear,
                          -((double)u)/TIDE_BLEND_SECONDS);
    }
  }

  /*
   * Else, we're far enough from newyears to ignore the blending.
   */
  return _time2dt_tide(t, deriv);
}

// This doesn't do blending; it's only used in
// SubordinateStation::predictApproximate.  **** Code duplicated from
// time2dt_tide and _time2dt_tide to avoid slowing down the inner
// loop.
double
ConstantSetWrapper::time2mean (Timestamp in_timestamp)
{
  // For starters, get us in the right year.
  Year in_year = in_timestamp.year();
  if (in_year != currentYear) {
    assert (!(currentYear.isNull()));
    currentYear = in_year;
    refresh_yearConstants ();
  }

  // This amount is less than a year, so a long int suffices.
  long t = (in_timestamp - epoch).in_seconds();

  double dt_tide = 0.0;
  unsigned a;

  for (a=0; a<length; a++)
    if (speeds[a] <= longtermspeed)
      dt_tide += amplitudes[a] * cos (speeds[a] * t + phases[a]);
  return dt_tide;
}

#if 0
ostream &operator<< (ostream &out, const ConstantSetWrapper &s) {
  out << "  Datum: " << s.datum() << endl;
  out << "  Constants:" << endl;
  for (unsigned i=0; i<s.length; i++)
    out << "    " << (s.origConstants)->amplitudes[i] << "  " <<
      (s.origConstants)->phases[i] << endl;
  return out;
}
#endif