// $Id: SubordinateStation.cc,v 1.2 2003/02/04 15:53:46 flaterco Exp $
/*  SubordinateStation  Implementation for offset station.

    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"

SubordinateStation::SubordinateStation (ReferenceStation *rs,
SubordinateStationRef *ssr): Station (rs->context) {
  stationRef = ssr;
  name = ssr->name;
  timeZone = ssr->timezone;
  note = ssr->note;
  if (timeZone.isNull())
    timeZone = rs->timeZone;
  coordinates = ssr->coordinates;
  constants = rs->constants;
  isCurrent = rs->isCurrent;
  isHydraulicCurrent = rs->isHydraulicCurrent;
  myUnits = rs->myUnits;
  origOffsets = ssr->offsets;
  residualOffsets = NULL;
  reduce_offsets();

  rs->constants = NULL;
  delete rs;
}

PredictionValue SubordinateStation::minLevel() const {
  PredictionValue pv = Station::minLevel();
  if (residualOffsets) {
    // Yet again...
    if (residualOffsets->min.levelAdd.Units() != pv.Units())
      residualOffsets->min.levelAdd.Units (pv.Units());
    pv = pv * residualOffsets->min.levelMultiply() +
              residualOffsets->min.levelAdd;
  }
  return pv;
}

PredictionValue SubordinateStation::maxLevel() const {
  PredictionValue pv = Station::maxLevel();
  if (residualOffsets) {
    // And again.
    if (residualOffsets->max.levelAdd.Units() != pv.Units())
      residualOffsets->max.levelAdd.Units (pv.Units());
    pv = pv * residualOffsets->max.levelMultiply() +
              residualOffsets->max.levelAdd;
  }
  return pv;
}

void
SubordinateStation::reduce_offsets () {
  if (residualOffsets)
    delete residualOffsets;
  residualOffsets = NULL;
  assert (origOffsets);

  // Step 1:  reduce hairy offsets to simple wherever possible (catch
  // hydraulics in step 2).  While we're at it, extract current angles.

  Offsets *mid = NULL;
  if (origOffsets->is_simple())
    mid = new SimpleOffsets (*((SimpleOffsets *)origOffsets));
  else {
    HairyOffsets *h = (HairyOffsets *)origOffsets;
    if (isCurrent) { // Extract current angles.
      minCurrentAngle = h->min.currentAngle;
      maxCurrentAngle = h->max.currentAngle;
    }
    if (h->max == h->min &&
         ((!isCurrent) || (h->max.timeAdd == h->floodbegins.timeAdd &&
                           h->max.timeAdd == h->ebbbegins.timeAdd))) {
      mid = new SimpleOffsets (h->max);
#ifdef SUPER_ULTRA_VERBOSE_DEBUGGING
      cerr << "Reduced HairyOffsets to simple" << endl;
#endif
    } else
      mid = new HairyOffsets (*h);
  }

  // Step 2:  extract whatever SimpleOffsets are relevant, and determine
  // if any residual offsets remain.

  SimpleOffsets s;
  if (mid->is_simple()) {
    SimpleOffsets *mp = (SimpleOffsets*)mid;
    if ((!isHydraulicCurrent)||(mp->levelMultiply() == 1.0))
      s = *mp;
    else {
      // We actually have to "un-simplify" this special case.  In the
      // case of hydraulic currents, level multipliers cannot be
      // treated as simple and applied to the constants because of the
      // square root operation.
      residualOffsets = new HairyOffsets(*mp);
      residualOffsets->max = residualOffsets->min = *mp;
      residualOffsets->floodbegins.timeAdd =
        residualOffsets->ebbbegins.timeAdd =  mp->timeAdd;
#ifdef SUPER_ULTRA_VERBOSE_DEBUGGING
      cerr << "Un-simplified offsets for hydraulic station" << endl;
#endif
    }
  } else {
    HairyOffsets *mp = (HairyOffsets*)mid;
    residualOffsets = new HairyOffsets(*mp);

    // We already know that this can't be reduced to SimpleOffsets.
    // There isn't any point in factoring out individual fields since
    // hairy offsets need to be processed in the prediction routine in
    // both cases.  However, it's important to center the time to
    // minimize the inaccuracy in predictApproximate for slack
    // offsets.  (If you were to factor out levelAdd and
    // levelMultiply, you would still need to avoid the pitfall of
    // messing up the order of operations.  For example, factoring out
    // a levelAdd but not a levelMultiply would yield incorrect
    // results.)
  }

  // Step 3:  if we have residual offsets, center the timeAdd.
  if (residualOffsets) {
    // If current, then center time for slack.  Otherwise, center
    // time over max and min.
    if (isCurrent) {
      s.timeAdd = (residualOffsets->floodbegins.timeAdd + residualOffsets->ebbbegins.timeAdd) / 2;
      residualOffsets->floodbegins.timeAdd -= s.timeAdd;
      residualOffsets->ebbbegins.timeAdd -= s.timeAdd;
    } else
      s.timeAdd = (residualOffsets->max.timeAdd + residualOffsets->min.timeAdd) / 2;
    residualOffsets->max.timeAdd -= s.timeAdd;
    residualOffsets->min.timeAdd -= s.timeAdd;
  }

  // Step 4:  apply SimpleOffsets.
  constants->setSimpleOffsets (s);

#ifdef SUPER_ULTRA_VERBOSE_DEBUGGING
  cerr << "Orig offsets: " << origOffsets << endl;
  cerr << "Mid-offsets: " << mid << endl;
  cerr << "Applied offsets: " << &s << endl;
  cerr << "Residual offsets: " << residualOffsets << endl;
#endif

  // Just "delete mid" yields a segmentation violation.
  if (mid->is_simple()) {
    SimpleOffsets *mp = (SimpleOffsets*)mid;
    delete mp;
  } else {
    HairyOffsets *mp = (HairyOffsets*)mid;
    delete mp;
  }
}

SubordinateStation::~SubordinateStation() {
  if (constants)
    delete constants;
  if (markLevel)
    delete markLevel;
}

// tm is "uncorrected" or "internal" timestamp
// t_out is "corrected" timestamp (with offsets applied)
// For reference station or simple offsets, they are identical.
// pv_out is undefined if the event is a sun or moon event.
void SubordinateStation::predictExactTideEvent (Timestamp &tm, Direction d,
Timestamp &t_out, EventType &etype_out, Dstr &etype_desc_long,
Dstr &etype_desc_short, PredictionValue &pv_out) {
  _predictExactTideEvent (tm, d, etype_out);

  t_out = tm;
  if (!isSunMoonEvent (etype_out)) {
    // Explicit Station:: gets the original, un-corrected value
    pv_out = Station::predictApproximate (tm);

    // It may be inefficient to do this here, but it's a lot more robust.
    if (residualOffsets) {
      if (residualOffsets->max.levelAdd.Units() != pv_out.Units())
        residualOffsets->max.levelAdd.Units (pv_out.Units());
      if (residualOffsets->min.levelAdd.Units() != pv_out.Units())
        residualOffsets->min.levelAdd.Units (pv_out.Units());
    }
  }

  // Apply residual offsets.
  if (residualOffsets) {
    switch (etype_out) {
    case max:
      t_out += residualOffsets->max.timeAdd;
      pv_out = pv_out * residualOffsets->max.levelMultiply()
		      + residualOffsets->max.levelAdd;
      break;
    case min:
      t_out += residualOffsets->min.timeAdd;
      pv_out = pv_out * residualOffsets->min.levelMultiply()
		      + residualOffsets->min.levelAdd;
      break;
    case slackrise:
      t_out += residualOffsets->floodbegins.timeAdd;
      break;
    case slackfall:
      t_out += residualOffsets->ebbbegins.timeAdd;
      break;
    // Mark should never be allowed if there are residual offsets.
    case markrise:
    case markfall:
      {
        Dstr details ("The offending station is ");
        details += name;
        details += ".";
        barf (ILLEGAL_MARK, details);
      }
    default:
      ;
    }
  }

  // The etype_desc figuring needs to be delayed until after everything
  // so that "Max Flood" versus "Min Ebb" will be correct when offsets
  // bump something over the line.
  etypedesc (etype_out, pv_out, etype_desc_long, etype_desc_short);
}


// Interpolate approximate tide levels for subordinate stations.  This is
// the most difficult algorithm in the whole program.  When NOAA needed to
// do this, they just used a simple sine curve to interpolate.  Nobody
// liked that answer.
//
// The original algorithm submitted by Jean-Pierre Lapointe
// (scipur@collegenotre-dame.qc.ca), which is in block comments far
// below, got me past the question of how to do the interpolation, but
// assumed the availability of usable high and low tides for the day
// of prediction.  Getting these tides turned out to be highly
// non-trivial.  It is not as simple as just using the two nearest
// bracketing tide events for the following reasons:
//
//    1.  Some locations like to produce lots of maxima that are
//    actually below the mean, and minima that are above the mean.
//    These data points cause the algorithm to blow up, and must be
//    thrown out.
//
//    2.  The lower high and higher low tide events produced by mostly
//    diurnal tides, as well as the various other "meaningless" maxima
//    and minima that happen with such tides, may be closer together
//    than the offset intervals.  In these cases it is better to just
//    apply offsets to the higher high and lower low, and leave the
//    others to be interpolated as they will.  Ed Wallner makes the
//    following relevant quote in his diurnal.txt file:
//
//       The NOS tide table notes on using tidal differences include
//       the statement that " ... at subordinate stations where the
//       tide is chiefly diurnal, the tide corrections are intended to
//       be used to approximate the times and heights of the higher
//       high and lower low waters. When the lower high water and
//       higher low water at the reference station are nearly the same
//       height, great reliance should not be placed on the calculated
//       corresponding tides at the subordinate station."
//
//    In these rare cases, it would actually be more accurate to
//    derive tide events from predictApproximate than to faithfully
//    apply the offsets with predictExact, but it is more important to
//    produce consistent results than to try to make sense of a
//    degenerate case.
//
//    3.  The choices of highs and lows to use must be deterministic
//    and consistent based on a given Timestamp input.  This
//    complicates the caching of values, but failure to cache here
//    would lead to totally unacceptable performance.
//
// The implementation of this algorithm in XTide 1 was not too smooth
// in its approach to solving these problems.  It used a fixed search
// interval (originally 26 hours, later expanded to 30) and used the
// highest high and lowest low tides within that interval, with
// kludges to stretch the interval as needed to solve problem #1.  To
// make a long story short, it didn't really work; it unfairly
// penalized the accuracy of regular old semi-diurnal stations, and it
// was not consistent.  The interval had to be expanded because Dublon
// managed to produce a nearly 29-hour period with no high tide events
// and crash the algorithm (from 1997-09-10 19:18 to 1997-09-12 00:13
// NCST).
//
// A large part of the reason for using this crummy approach in XTide
// 1 was the lack of a facility for searching backwards for tide
// events.  Finding previous tide events without a fixed interval to
// search in would have been prohibitive.  This problem has been
// corrected in XTide 2.
//
// XTide 2 doesn't use a fixed interval.  It just keeps looking for
// tide events until it finds some that are a "reasonable distance"
// away from the mean.  "Reasonable" is defined as a fraction of the
// total theoretical amplitude in config.hh.
//
// Since I use a localized mean instead of a global one, we shouldn't
// have to look too far to find a reasonable answer (provided that
// reasonable_fraction is reasonably defined).  The caching problem is
// solved by using the actual two nearest tide events as a bracket,
// and forcing a refresh when the bracket is exited.
//
// NOTE 1:  It is still theoretically possible to obtain inconsistent
// results if a tide is on the very border of "reasonable" and gets
// thrown over by inaccuracy in the third decimal of a tide
// prediction.  If you can demonstrate this bug in practice, you
// deserve a medal or something.
//
// NOTE 2:  This algorithm does nothing about slack water offsets.

// Following is Jean-Pierre's explanation of his original algorithm.
// Things that I changed are:
//
//    (1) Use a moving mean instead of the datum.  Using the datum
//    gives unusable results in summer and winter for locations with
//    large seasonal variations.
//
//    (2) Interpolate multiplicative offsets as well as additive ones.
//    Jean-Pierre figured that the changeover from low-tide
//    multipliers to high-tide multipliers always happened at zero,
//    which is not the case here.  We need to interpolate to avoid
//    causing a discontinuity.

/************** COMMENT BLOCK ************************

1-  Store the following data
     Z  (datum at reference station)
     t  (time of prediction at the secondary station)
     TH  (correction for high tide time at secondary station)
     TL  (correction for low tide time at secondary station)
     RH  (correction ratio for high tide height at secondary station)
     RL  (correction ratio for low tide height at secondary station)
     HH  (height correction for high tide height at secondary station)
     HL  (height correction for low tide height at secondary station)

2-  Run XTIDE for the reference station for the day of prediction
     to find the height of the higher tide and store it as MAX and
     to find the height of the lower tide and store it as MIN.

3-  Run XTIDE for the reference station at a time T defined as
     T =  t - (TH + TL) / 2
    Store the height calculated by XTIDE as HI (intermediate height)

4-  Store S defined as S =HI - Z

5-  Run XTIDE for the reference station at a time T defined as:
     if S > 0 then
     T = t - (TH + TL) / 2 - (TH-TL) / 2 x S / absolute value(MAX - Z)
     else
     T = t - (TH + TL) / 2 - (TH-TL) / 2 x S / absolute value(MIN - Z)

    Store the height calculated by XTIDE as HI (intermediate height)
    and calculate HS (height at secondary station at time t) defined as:
    if S > 0 then
    HS =  HI x RH + (HH + HL) / 2 + (HH-HL) / 2 x S/absolute value(MAX - Z)
    else
    HS =  HI x RL + (HH + HL) / 2 + (HH-HL) / 2 x S/absolute value(MIN - Z)

    You now have HS the height of the tide at the secondary station at a
    time t for this station.

********  END COMMENT BLOCK  **********************************************/

PredictionValue SubordinateStation::predictApproximate (Timestamp tm) {
  // Get rid of the trivial cases.
  if (!residualOffsets)
    return Station::predictApproximate (tm);

  // Invalidate cache?
  if (!bracketleft.isNull())
    if (tm < bracketleft || tm > bracketright || usemax.Units() != myUnits)
      bracketleft.make_null();

  // Refresh cache?

  // This is tricky.  A criterion of "nearest reasonable high/low
  // tide" would violate the requirement for consistency -- you could
  // get different results based on where in the interval you are when
  // you refresh the cache.  One might instead constrain the direction
  // of the search so that we always get one tide on each side, of the
  // same type as the bracketing tides.  But this causes
  // discontinuities since the value being used in lieu of an
  // unreasonable max or min will change as you cross that max or min.
  // So, we have to be extra-cynical, and just always search in the
  // same direction.
  if (bracketleft.isNull()) {
#ifdef SUPER_ULTRA_VERBOSE_DEBUGGING
    cerr << "Refreshing predictApproximate cache" << endl;
#endif
    int bhiflag=0, fhiflag=0;
    Timestamp b, f;
    PredictionValue pv;

    // It's plausible that we might miss a tide if we are sitting right
    // on top of it -- detect that and correct for it.
    {
      Timestamp ttm = tm;
      int itcnt = 0;
      do {
        b = bracketleft = next_high_or_low_tide (ttm, bhiflag, backward);
        f = bracketright = next_high_or_low_tide (ttm, fhiflag, forward);
        ttm += Interval(60);
        assert (itcnt++ < 2);
      } while (bhiflag == fhiflag);
#ifdef SUPER_ULTRA_VERBOSE_DEBUGGING
      if (itcnt > 1)
        cerr << "  ** Corrected for bhiflag == bloflag!" << endl;
#endif
    }

    // Make the left side reasonable.
    pv = Station::predictApproximate(b);
    while (!(isReasonable (b, pv, bhiflag))) {
#ifdef SUPER_ULTRA_VERBOSE_DEBUGGING
      cerr << "  * Left side unreasonable" << endl;
#endif
      int tempflag;
      do
        b = next_high_or_low_tide (b, tempflag, forward); // See comments.
      while (tempflag != bhiflag);
      pv = Station::predictApproximate(b);
    }
    if (bhiflag)
      usemax = pv;
    else
      usemin = pv;

    // Make the right side reasonable.
    pv = Station::predictApproximate(f);
    while (!(isReasonable (f, pv, fhiflag))) {
#ifdef SUPER_ULTRA_VERBOSE_DEBUGGING
      cerr << "  * Right side unreasonable" << endl;
#endif
      int tempflag;
      do
        f = next_high_or_low_tide (f, tempflag, forward);
      while (tempflag != fhiflag);
      pv = Station::predictApproximate(f);
    }
    if (fhiflag)
      usemax = pv;
    else
      usemin = pv;
    assert (bhiflag != fhiflag);

#ifdef SUPER_ULTRA_VERBOSE_DEBUGGING
    cerr << "predictApproximate cache refreshed" << endl;
    cerr << "   timestamp: " << tm << endl;
    cerr << "   l bracket: " << bracketleft << endl;
    cerr << "   r bracket: " << bracketright << endl;
    cerr << "   l using: " << b << endl;
    cerr << "   r using: " << f << endl;
    cerr << "   usemax: " << usemax << endl;
    cerr << "   usemin: " << usemin << endl;
#endif
  }

  // Do the interpolation.  The capital letter variable names here
  // match up with those used in the algorithm in the comment block
  // above.

  /* This is the initial guess (average of time offsets) */
  Timestamp T = tm - (residualOffsets->max.timeAdd +
                      residualOffsets->min.timeAdd) / 2;
  PredictionValue Z = movingMean(T);
  PredictionValue S = Station::predictApproximate(T) - Z;

  /* Improve the initial guess. */
  double magicnum;
  if (S.val() > 0.0)
    magicnum = 0.5 * S.val() / fabs(usemax.val() - Z.val());
  else
    magicnum = 0.5 * S.val() / fabs(usemin.val() - Z.val());
  T -= (residualOffsets->max.timeAdd - residualOffsets->min.timeAdd) * magicnum;
  PredictionValue HI = Station::predictApproximate(T);

  /* Apply the height offsets. */
  // It may be inefficient to do this here, but it's a lot more robust.
  if (residualOffsets->max.levelAdd.Units() != HI.Units())
    residualOffsets->max.levelAdd.Units (HI.Units());
  if (residualOffsets->min.levelAdd.Units() != HI.Units())
    residualOffsets->min.levelAdd.Units (HI.Units());
  return HI * ((residualOffsets->max.levelMultiply() +
                residualOffsets->min.levelMultiply()) / 2.0 +
               (residualOffsets->max.levelMultiply() -
                residualOffsets->min.levelMultiply()) * magicnum)
             + (residualOffsets->max.levelAdd +
                residualOffsets->min.levelAdd) / 2.0 +
               (residualOffsets->max.levelAdd -
                residualOffsets->min.levelAdd) * magicnum;
}

int SubordinateStation::isReasonable (Timestamp tm, PredictionValue
pv, int ishigh) {
  // Calculate a fresh mean when testing for reasonableness so that
  // the decision will be most consistent.
  PredictionValue mean = movingMean (tm);
  if (ishigh)
    return ((pv - mean) / (constants->maxAmplitude) > reasonable_fraction);
  else
    return ((mean - pv) / (constants->maxAmplitude) > reasonable_fraction);
}

int SubordinateStation::is_reference_station() {
  return 0;
}