// $Id: Station.cc,v 1.8 2003/02/20 15:19:56 flaterco Exp $
/*  Station  Abstract superclass of ReferenceStation and SubordinateStation.

    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"

Station::Station (TideContext *in_context) {
  assert (in_context);
  context = in_context;
  markLevel = NULL;
  aspect = context->settings->ga;
  step = Interval (HOURSECONDS);
}

PredictionValue Station::minLevel() const {
  return constants->datum() - constants->maxAmplitude;
}

PredictionValue Station::maxLevel() const {
  return constants->datum() + constants->maxAmplitude;
}

#if 0
ostream &operator<< (ostream &out, const Station &s) {
  out << "Station name: " << s.name << endl;
  out << "Station time zone: " << s.timeZone << endl;
  out << "Station coordinates: " << s.coordinates << endl;
  out << "Is current: " << s.isCurrent << endl;
  out << "Max level: " << s.maxLevel() << endl;
  out << "Min level: " << s.minLevel() << endl;
  out << "Constants:" << endl;
  out << *(s.constants);
  return out;
}
#endif

// 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
 *
 ************************************************************************/

/* TIDE_TIME_PREC
 *   Precision (in seconds) to which we will find roots
 *   See config.hh for def_TIDE_TIME_PREC.
 */
// This initialization fails on AIX 4.2.1, IBM RS6000, 43P-140,
// yielding a zero interval (Alan J. Wylie).  See _predictExactTideEvent
// for workaround.
static Interval TIDE_TIME_PREC (def_TIDE_TIME_PREC);

/* TIDE_TIME_STEP
 *   We are guaranteed to find all high and low tides as long as their
 * spacing is greater than this value (in seconds).
 */
#define TIDE_TIME_STEP (TIDE_TIME_PREC)

// Functions to zero out

// Option #1 -- find maxima and minima
PredictionValue
Station::f_hiorlo (Timestamp t, unsigned deriv) {
  return constants->predictHeight (t, deriv+1);
}

// Option #2 -- find mark crossings or slack water
// Marklev is a class-level variable that is set by predictExactTideEvent.
// ** Marklev must be made compatible with the tide as returned by
// predictHeight, i.e., no datum, no conversion from KnotsSquared.
PredictionValue
Station::f_mark (Timestamp t, unsigned deriv)
{
  PredictionValue pv_out = constants->predictHeight (t, deriv);
  if (deriv == 0)
    pv_out -= marklev;
  return pv_out;
}

/* find_zero (time_t t1, time_t t2, double (*f)(time_t t, int deriv))
 *   Find a zero of the function f, which is bracketed by t1 and t2.
 *   Returns a value which is either an exact zero of f, or slightly
 *   past the zero of f.
 */

/*
 * Here's a root finder based upon a modified Newton-Raphson method.
 */

// The direction parameter just controls which side of the bracket
// we return.

Timestamp
Station::find_zero (Timestamp tl, Timestamp tr,
PredictionValue (Station::*f)(Timestamp, unsigned deriv), Direction d) {
  PredictionValue fl = (this->*f)(tl, 0);
  PredictionValue fr = (this->*f)(tr, 0);
  double scale = 1.0;
  Interval dt;
  Timestamp t;
  PredictionValue fp, ft, f_thresh;
#ifdef SUPER_ULTRA_VERBOSE_DEBUGGING
  long itctr = 0;
#endif

  assert(fl.val() != 0.0 && fr.val() != 0.0);
  assert(tl < tr);
  if (fl.val() > 0) {
    scale = -1.0;
    fl = -fl;
    fr = -fr;
  }
  assert(fl.val() < 0.0 && fr.val() > 0.0);

  while (tr - tl > TIDE_TIME_PREC) {
#ifdef SUPER_ULTRA_VERBOSE_DEBUGGING
    itctr++;
    cout << tr << " - " << tl << " is > " << TIDE_TIME_PREC << endl;
    cout << tr.timet() << " " << tl.timet() << endl;
#endif

    if (t.isNull())
      dt = Interval(0); // Force bisection on first step
    else if (abs(ft) > f_thresh   /* not decreasing fast enough */
	|| (ft.val() > 0.0 ?    /* newton step would go outside bracket */
	    (fp <= ft / (long)((t - tl).in_seconds())) :
	    (fp <= -ft / (long)((tr - t).in_seconds()))))
      dt = Interval(0); /* Force bisection */
    else {
      /* Attempt a newton step */
      assert (fp.val() != 0.0);
      dt = Interval((long)floor(-ft/fp + 0.5));

      /* Since our goal specifically is to reduce our bracket size
	 as quickly as possible (rather than getting as close to
	 the zero as possible) we should ensure that we don't take
	 steps which are too small.  (We'd much rather step over
	 the root than take a series of steps which approach the
	 root rapidly but from only one side.) */
      if (abs(dt) < TIDE_TIME_PREC)
	  dt = ft.val() < 0.0 ? TIDE_TIME_PREC : -TIDE_TIME_PREC;

      if ((t += dt) >= tr || t <= tl)
	  dt = Interval(0);   /* Force bisection if outside bracket */
      f_thresh = abs(ft) / 2.0;
    }

    if (dt.in_seconds() == 0) {
      /* Newton step failed, do bisection */
      t = tl + (tr - tl) / 2;
      f_thresh = fr > -fl ? fr : -fl;
    }

    if ((ft = scale * (this->*f)(t,0)).val() == 0.0) {
#ifdef SUPER_ULTRA_VERBOSE_DEBUGGING
        cerr << "find_zero took " << itctr << " iterations" << endl;
#endif
	return t;             /* Exact zero */
    } else if (ft.val() > 0.0)
	tr = t, fr = ft;
    else
	tl = t, fl = ft;

    fp = scale * (this->*f)(t,1);
  }

#ifdef SUPER_ULTRA_VERBOSE_DEBUGGING
  cerr << "find_zero took " << itctr << " iterations" << endl;
#endif
  if (d == forward)
    return tr;
  return tl;
}

/* next_zero(time_t t, double (*f)(), double max_fp, double max_fpp)
 *   Find the next zero of the function f which occurs after time t.
 *   The arguments max_fp and max_fpp give the maximum possible magnitudes
 *   that the first and second derivative of f can achieve.
 *
 *   Algorithm:  Our goal here is to bracket the next zero of f ---
 *     then we can use find_zero() to quickly refine the root.
 *     So, we will step forward in time until the sign of f changes,
 *     at which point we know we have bracketed a root.
 *     The trick is to use large steps in our search, making
 *     sure the steps are not so large that we inadvertently
 *     step over more than one root.
 *
 *     The big trick, is that since the tides (and derivatives of
 *     the tides) are all just harmonic series', it is easy to place
 *     absolute bounds on their values.
 */

// I've munged this, via clone and hack, to search backwards and forwards.

// This function is ONLY used for finding maxima and minima.  Since by
// definition the tide cannot change direction between maxima and
// minima, there is at most one crossing of a given mark level between
// min/max points.  Therefore, we already have a bracket of the mark
// level to give to find_zero.

Timestamp
Station::next_zero (Timestamp t,
PredictionValue (Station::*f)(Timestamp, unsigned deriv),
int &risingflag, Amplitude max_fp, Amplitude max_fpp, Direction d) {
  Timestamp t_left, t_right;
  Interval step, step1, step2;
  PredictionValue f_left, df_left, f_right, df_right;
  double scale = 1.0;
#ifdef SUPER_ULTRA_VERBOSE_DEBUGGING
  long itctr = 0;
#endif

  if (d == forward) {
    // This is the original forwards version.

    t_left = t;

    /* If we start at a zero, step forward until we're past it. */
    while ((f_left = (this->*f)(t_left,0)).val() == 0.0)
      t_left += TIDE_TIME_PREC;

    if (!(risingflag = f_left.val() < 0.0)) {
      scale = -1.0;
      f_left = -f_left;
    }

    while (1) {
#ifdef SUPER_ULTRA_VERBOSE_DEBUGGING
      itctr++;
#endif

      /* Minimum time to next zero: */
      step1 = Interval((long)(abs(f_left) / max_fp));

      /* Minimum time to next turning point: */
      df_left = scale * (this->*f)(t_left,1);
      step2 = Interval((long)(abs(df_left) / max_fpp));

      if (df_left.val() < 0.0)
	/* Derivative is in the wrong direction. */
	step = step1 + step2;
      else
	step = step1 > step2 ? step1 : step2;

      if (step < TIDE_TIME_STEP)
	  step = TIDE_TIME_STEP; /* No ridiculously small steps */

      t_right = t_left + step;
      /*
       * If we hit upon an exact zero, step right until we're off
       * the zero.  If the sign has changed, we are bracketing a desired
       * root, if the sign hasn't changed, then the zero was at
       * an inflection point (i.e. a double-zero to within TIDE_TIME_PREC)
       * and we want to ignore it.
       */
      while ((f_right = scale * (this->*f)(t_right, 0)).val() == 0.0)
	t_right += TIDE_TIME_PREC;
      if (f_right.val() > 0.0) {
#ifdef SUPER_ULTRA_VERBOSE_DEBUGGING
        cerr << "next_zero took " << itctr << " iterations" << endl;
#endif
	return find_zero(t_left, t_right, f, d); /* Found a bracket */
      }

      t_left = t_right, f_left = f_right;
    }

  } else {
    // This is the mirror image of the above.

    t_right = t;

    /* If we start at a zero, step backward until we're past it. */
    while ((f_right = (this->*f)(t_right,0)).val() == 0.0)
      t_right -= TIDE_TIME_PREC;

    if ((risingflag = f_right.val() > 0.0)) {
      scale = -1.0;
      f_right = -f_right;
    }

    while (1) {
#ifdef SUPER_ULTRA_VERBOSE_DEBUGGING
      itctr++;
#endif

      /* Minimum time to next zero: */
      step1 = Interval((long)(abs(f_right) / max_fp));

      /* Minimum time to next turning point: */
      df_right = scale * (this->*f)(t_right,1);
      step2 = Interval((long)(abs(df_right) / max_fpp));

      if (df_right.val() > 0.0)
	/* Derivative is in the wrong direction. */
	step = step1 + step2;
      else
	step = step1 > step2 ? step1 : step2;

      if (step < TIDE_TIME_STEP)
	  step = TIDE_TIME_STEP; /* No ridiculously small steps */

      t_left = t_right - step;
      /*
       * If we hit upon an exact zero, step left until we're off
       * the zero.  If the sign has changed, we are bracketing a desired
       * root, if the sign hasn't changed, then the zero was at
       * an inflection point (i.e. a double-zero to within TIDE_TIME_PREC)
       * and we want to ignore it.
       */
      while ((f_left = scale * (this->*f)(t_left, 0)).val() == 0.0)
	t_left -= TIDE_TIME_PREC;
      if (f_left.val() > 0.0) {
#ifdef SUPER_ULTRA_VERBOSE_DEBUGGING
        cerr << "next_zero took " << itctr << " iterations" << endl;
#endif
	return find_zero(t_left, t_right, f, d); /* Found a bracket */
      }

      t_right = t_left, f_right = f_left;
    }
  }
}

Timestamp
Station::next_high_or_low_tide (Timestamp t, int &hiflag, Direction d) {
  int rising;
  Timestamp thilo = next_zero (t, &Station::f_hiorlo, rising,
    constants->dt_tide_max(2), constants->dt_tide_max(3), d);
  hiflag = !rising;
  return thilo;
}

// Find the mark crossing in this bracket, or return NULL if there is
// none.
Timestamp
Station::find_mark_crossing (Timestamp t1, Timestamp t2, int &risingflag,
Direction d)
{
  if (t1 > t2) {
    Timestamp t = t1;
    t1 = t2;
    t2 = t;
  }

  PredictionValue f1 = f_mark(t1,0);
  PredictionValue f2 = f_mark(t2,0);

  assert (f1 != f2);

  if (!(risingflag = f1.val() < 0.0 || f2.val() > 0.0)) {
     f1 = -f1;
     f2 = -f2;
  }

  // Since f1 != f2, we can't get two zeros, so it doesn't matter which
  // one we check first.
  if (f1.val() == 0.0)
    return t1;
  else if (f2.val() == 0.0)
    return t2;

  if (f1.val() < 0.0 && f2.val() > 0.0)
    return find_zero (t1, t2, &Station::f_mark, d);

  return Timestamp(); // Don't have a bracket, return null timestamp
}

// This has changed radically and gotten really complicated.
// (Sun and moon didn't help.)
//
// 2003-02-04
// Considered replacing the whole thing with a generic sort but
// decided just to make it more complicated again.  Extended for
// moonrise and moonset.
//
// tm is "uncorrected" or "internal" timestamp.
void Station::_predictExactTideEvent (Timestamp &tm, Direction d,
EventType &etype_out) {

  int ns = (context->settings->ns != 'n');

  // Workaround for failed initialization of TIDE_TIME_PREC on
  // AIX 4.2.1, IBM RS6000, 43P-140 (Alan J. Wylie)
  TIDE_TIME_PREC = Interval (def_TIDE_TIME_PREC);

  // Invalidate cache if necessary.
  if (d != cache_d || tm != last_tm) {
    t_hilo.make_null();
    t_slack.make_null();
    t_mark.make_null();
    t_moonrs.make_null();
    t_moonphase.make_null();
    t_sunrs.make_null();
  }

  if (!ns) {
    // Update sun
    if (t_sunrs.isNull() && (!(coordinates.isNull()))) {
      t_sunrs = tm;
      find_next_rise_or_set (t_sunrs, coordinates, sunrs_t, d, 0);
    }
    // Update moon
    if (t_moonrs.isNull() && (!(coordinates.isNull()))) {
      // Danger, danger.  The moonrise and moonset logic blows up if you
      // go before 1900 or after 2099.  Leave a safety margin.
      // This logic is straight out of skycal.cc.
      double danger_y = 1900. + (tm.jd() - 2415019.5) / 365.25;
      if (danger_y > 1900.1 && danger_y < 2099.9) {
        t_moonrs = tm;
        find_next_rise_or_set (t_moonrs, coordinates, moonrs_t, d, 1);
      }
    }
    if (t_moonphase.isNull()) {
      t_moonphase = tm;
      find_next_moon_phase (t_moonphase, moonphase_t, d);
    }
  }

  // If t_hilo is a valid cached value, then we needn't check again
  // for t_slack or t_mark -- they will also be up-to-date.  See the
  // next big block of comments below for more explanation.

  if (t_hilo.isNull()) {
    int risingflag;

    /* Find next high/low tide */
    t_hilo = next_high_or_low_tide(tm, risingflag, d);
    assert ((d == forward && t_hilo > tm) ||
            (d == backward && t_hilo < tm));
    if (risingflag)
      hilo_t = max;
    else
      hilo_t = min;

    // Check for slack water, if applicable.
    assert (t_slack.isNull());
    if (isCurrent) {
      // marklev must compensate for datum and KnotsSquared.  See f_mark.
      marklev = -(constants->datum());
      if (constants->predictUnits() != marklev.Units())
	marklev.Units (constants->predictUnits());
      t_slack = find_mark_crossing (tm, t_hilo, risingflag, d);
      if (risingflag)
	slack_t = slackrise;
      else
	slack_t = slackfall;
    }

    // Check for mark, if applicable.
    assert (t_mark.isNull());
    if (markLevel) {
      // marklev must compensate for datum and KnotsSquared.  See f_mark.
      marklev = *markLevel;
      // Correct meters / feet
      if (marklev.Units() != constants->datum().Units())
        marklev.Units (constants->datum().Units());
      marklev -= constants->datum();
      // Correct knots / knots squared
      if (constants->predictUnits() != marklev.Units())
	marklev.Units (constants->predictUnits());
      t_mark = find_mark_crossing (tm, t_hilo, risingflag, d);
      if (risingflag)
	mark_t = markrise;
      else
	mark_t = markfall;
    }
  }

  // Sun and moon -- let's try to keep this simple.  (Ha ha ha!)
  if (!ns) {
    Timestamp *t_sunmoon = &t_moonphase;
    EventType sunmoon_t = moonphase_t;

    // Make t_sunmoon point to whichever of sunrs, moonrs, or moonphase
    // would win.  N.B., the sun and moon might never rise (particularly
    // if we have no coordinates), but the moon always has a phase.
    
    if (!(t_sunrs.isNull())) {
      if ((d == forward && t_sunrs < *t_sunmoon) ||
	  (d == backward && t_sunrs > *t_sunmoon)) {
	t_sunmoon = &t_sunrs;
	sunmoon_t = sunrs_t;
      }
    }
    if (!(t_moonrs.isNull())) {
      if ((d == forward && t_moonrs < *t_sunmoon) ||
	  (d == backward && t_moonrs > *t_sunmoon)) {
	t_sunmoon = &t_moonrs;
	sunmoon_t = moonrs_t;
      }
    }

    // While we're waiting for goto to come back into style...
    while (1) {
      if (d == forward) {
	if (!(t_hilo.isNull()))
	  if (t_hilo < *t_sunmoon)
	    break;
	if (!(t_mark.isNull()))
	  if (t_mark < *t_sunmoon)
	    break;
	if (!(t_slack.isNull()))
	  if (t_slack < *t_sunmoon)
	    break;
      } else {
	if (!(t_hilo.isNull()))
	  if (t_hilo > *t_sunmoon)
	    break;
	if (!(t_mark.isNull()))
	  if (t_mark > *t_sunmoon)
	    break;
	if (!(t_slack.isNull()))
	  if (t_slack > *t_sunmoon)
	    break;
      }

      // Sun/moon wins.
      last_tm = tm = *t_sunmoon;
      cache_d = d;
      etype_out = sunmoon_t;
      t_sunmoon->make_null(); // This is why I needed a pointer.
      return;
    }
  }

  // Okay, sun and moon are out of the picture.
  // If hilo is the next event, then t_slack and t_mark must both be
  // null, since find_zero returns null if they are not bracketed by
  // tm and t_hilo, and we nullify the cached values after they are
  // used.  Furthermore, if t_slack or t_mark are set, then they are
  // necessarily ahead of t_hilo in line.  So the only real comparison
  // that needs to be made is whether t_slack or t_mark comes first.

  // Did that make sense?

  int do_mark = 0;
  if (!(t_mark.isNull())) {
    if (!(t_slack.isNull())) {
      if ((t_mark < t_slack && d == forward) ||
          (t_mark > t_slack && d == backward))
        do_mark = 1;
    } else do_mark = 1;
  }
  if (do_mark) {
    tm = t_mark;
    etype_out = mark_t;
    t_mark.make_null();
  } else {
    if (!(t_slack.isNull())) {
      tm = t_slack;
      etype_out = slack_t;
      t_slack.make_null();
    } else {
      tm = t_hilo;
      etype_out = hilo_t;
      t_hilo.make_null();
    }
  }

  // Validate cache.
  last_tm = tm;
  cache_d = d;
}

int isSunMoonEvent (Station::EventType a) {
  switch (a) {
  case Station::sunrise:
  case Station::sunset:
  case Station::moonrise:
  case Station::moonset:
  case Station::newmoon:
  case Station::firstquarter:
  case Station::fullmoon:
  case Station::lastquarter:
    return 1;
  default:
    return 0;
  }
}

int isMaxMinEvent (Station::EventType a) {
  switch (a) {
  case Station::max:
  case Station::min:
    return 1;
  default:
    return 0;
  }
}

PredictionValue Station::predictApproximate (Timestamp tm) {
  PredictionValue pv_out;
  pv_out = constants->predictHeight (tm, 0);
  if (pv_out.Units().mytype == PredictionValue::Unit::KnotsSquared)
    pv_out.Units (PredictionValue::Unit::Knots);
  pv_out += constants->datum();
  return pv_out;
}

PredictionValue Station::movingMean (Timestamp tm) {
  PredictionValue pv_out;
  pv_out = constants->movingMean (tm);
  if (pv_out.Units().mytype == PredictionValue::Unit::KnotsSquared)
    pv_out.Units (PredictionValue::Unit::Knots);
  pv_out += constants->datum();
  return pv_out;
}

void Station::etypedesc (EventType etype, PredictionValue pv,
Dstr &etype_desc_long, Dstr &etype_desc_short) {
  switch (etype) {
  case max:
    etype_desc_long = (isCurrent ? (pv.val() >= 0.0 ? "Max Flood"
			       : "Min Ebb") : "High Tide");
    etype_desc_short = (isCurrent ? (pv.val() >= 0.0 ? "MaxFl"
			       : "MinEb") : "High");
    break;
  case min:
    etype_desc_long = (isCurrent ? (pv.val() <= 0.0 ? "Max Ebb"
                               : "Min Flood") : "Low Tide");
    etype_desc_short = (isCurrent ? (pv.val() <= 0.0 ? "MaxEb"
                               : "MinFl") : "Low");
    break;
  case slackrise:
    etype_desc_long = "Slack, Flood Begins";
    etype_desc_short = "Slack";
    break;
  case slackfall:
    etype_desc_long = "Slack, Ebb Begins";
    etype_desc_short = "Slack";
    break;
  case markrise:
    etype_desc_long = (isCurrent ? (pv.val() <= 0.0 ? "Mark, Ebb Decreasing"
			       : "Mark, Flood Increasing")
		  : "Mark Rising");
    etype_desc_short = "Mark";
    break;
  case markfall:
    etype_desc_long = (isCurrent ? (pv.val() >= 0.0 ? "Mark, Flood Decreasing"
			       : "Mark, Ebb Increasing")
		  : "Mark Falling");
    etype_desc_short = "Mark";
    break;
  case sunrise:
    etype_desc_long = "Sunrise";
    etype_desc_short = "Srise";
    break;
  case sunset:
    etype_desc_long = "Sunset";
    etype_desc_short = "Sset";
    break;
  case moonrise:
    etype_desc_long = "Moonrise";
    etype_desc_short = "Mrise";
    break;
  case moonset:
    etype_desc_long = "Moonset";
    etype_desc_short = "Mset";
    break;
  case newmoon:
    etype_desc_long = "New Moon";
    etype_desc_short = "NewM";
    break;
  case firstquarter:
    etype_desc_long = "First Quarter";
    etype_desc_short = "1stQ";
    break;
  case fullmoon:
    etype_desc_long = "Full Moon";
    etype_desc_short = "FullM";
    break;
  case lastquarter:
    etype_desc_long = "Last Quarter";
    etype_desc_short = "LastQ";
    break;
  default:
    assert (0);
  }
}

// N.B., see the .cc of xxDrawable::global_redraw() regarding an
// important dependency on how this is implemented.
Station *Station::clone() {
  Station *s = stationRef->load(context);
  if (markLevel)
    s->markLevel = new PredictionValue(*markLevel);
  else
    s->markLevel = NULL;
  s->aspect = aspect;
  s->step = step;
  if (s->myUnits != myUnits)
    s->setUnits (myUnits);
  return s;
}

void Station::setUnits (PredictionValue::Unit in_units) {
  // FIXME if there are ever units of velocity other than knots
  assert (!isCurrent);
  if (myUnits != in_units) {
    myUnits = in_units;
    constants->setUnits (in_units);
    if (markLevel)
      markLevel->Units (in_units);
  }
}

Station::~Station() {
  // I put everything in the subclasses because of some compiler
  // problem...
}