// -*- c++ -*-
// $Id: Congen 2632 2007-08-18 21:30:29Z flaterco $

/*  congen:  constituent generator.
    Copyright (C) 1997  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 3 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
    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, see <>.


#include <stdint.h>
#include <valarray>
#include <string>
#include <vector>
#include <istream>

// This file is (or was) automatically mangled to get rid of UTF-8
// characters like Greek letters and subscripts.  As of g++ 4.2.0,
// Greek letters in identifiers can be made to work (technically, if
// not readably) by translating them to \uNNNN syntax and using the
// -fextended-identifiers switch, but subscripts, primes, and double
// primes are still rejected by the compiler.
// "error: universal character \u2081 is not valid in an identifier"

/* The reader is assumed to be familiar with the following
   publication, hereafter referred to as "SP 98:"

   Manual of Harmonic Analysis and Prediction of Tides.  Special
   Publication No. 98, Revised (1940) Edition (reprinted 1958 with
   corrections; reprinted again 1994).  United States Government
   Printing Office, 1994.

   SP 98 does not cover Doodson-style constituents.  If using those,
   the following publication may be helpful:

   Foreman, M.G.G., 1977.  Manual for Tidal Heights Analysis and
   Prediction.  Pacific Marine Science Report 77-10, Institute of
   Ocean Sciences, Patricia Bay, Sidney, B.C. (2004 revision).

// Apart from the names of the constituents, the following
// astronomical terms from SP 98 are used in this header file:

// (Vₒ+u)   Constituent argument at beginning of a tidal series
// f        Node factor

// V        Principal portion of constituent argument
// T        Hour angle of mean sun
// s        Mean longitude of moon
// h        Mean longitude of sun
// p        Mean longitude of lunar perigee
// p₁       Mean longitude of solar perigee

// u        Part of constituent argument depending upon variations in
//          obliquity of lunar orbit
// ξ        Longitude in moon's orbit of lunar intersection
// ν        Right ascension of lunar intersection
// ν′       Term in argument of lunisolar constituent K₁
// 2ν″      Term in argument of lunisolar constituent K₂
// Q        Term in argument of constituent M₁
// R        Term in argument of constituent L₂
// Qᵤ       Term in argument of constituent M₁
// N        Longitude of moon's node

namespace Congen {

  // The interface to libcongen makes use of valarrays whose lengths
  // and contents are prescribed, but not checkable at compile time.
  // Thus there is the risk that future changes to the interface could
  // result in older programs compiling OK but running incorrectly.
  // To prevent this, client programs should contain a preprocessor
  // check like the following:
  //   #if Congen_interfaceRevision != 0
  //   #error trying to compile with an incompatible version of libcongen
  //   #endif
  // The value will increment only if and when one of those risky
  // changes is made to the interface.

  #define Congen_interfaceRevision 0

  // Similarly, the next function can be used in autoconf
  // ( scripts as follows to test that the library that
  // you are about to link with is the right version:
  //   AC_MSG_CHECKING([for libcongen (interfaceRevision_0)])
  //   LIBS="$LIBS -lcongen"
  //    [AC_LANG_PROGRAM([[#include <Congen>]],
  //                     [[Congen::interfaceRevision_0();]])],
  //     AC_MSG_RESULT([yes]),
  //    [AC_MSG_RESULT([NO])
  //     AC_MSG_ERROR([either cannot find libcongen or found wrong version; try setting LDFLAGS.])
  //    ])
  // If invoked it does nothing.

  void interfaceRevision_0 ();

  // Where valarrays are used to contain coefficients of the following
  // terms of V and u, these are the indices to the respective
  // coefficients.

  // V portion of argument, evaluated at beginning of year.
  const uint_fast8_t T_index   = 0;
  const uint_fast8_t s_index   = 1;
  const uint_fast8_t h_index   = 2;
  const uint_fast8_t p_index   = 3;
  const uint_fast8_t p₁_index  = 4;
  const uint_fast8_t c_index   = 5;    // Constant angle (degrees).
  const uint_fast8_t numVterms = 6;

  // u portion of argument, evaluated at mid-year.
  const uint_fast8_t ξ_index   = 0;
  const uint_fast8_t ν_index   = 1;
  const uint_fast8_t ν′_index  = 2;
  const uint_fast8_t 2ν″_index = 3;
  const uint_fast8_t Q_index   = 4;
  const uint_fast8_t R_index   = 5;
  const uint_fast8_t Qᵤ_index  = 6;
  const uint_fast8_t numuterms = 7;

  // With the exception of "unity," for which we use 1, Congen and SP
  // 98 Table 2 identify node factor formulas by their formula numbers
  // in SP 98.  The following constants are optional syntactic sugar
  // that may at least clarify matters for those without access to SP
  // 98.

  const uint_fast8_t f_1   = 1;
  const uint_fast8_t f_Mm  = 73;
  const uint_fast8_t f_Mf  = 74;
  const uint_fast8_t f_O₁  = 75;
  const uint_fast8_t f_J₁  = 76;
  const uint_fast8_t f_OO₁ = 77;
  const uint_fast8_t f_M₂  = 78;
  const uint_fast8_t f_KJ₂ = 79;
  const uint_fast8_t f_M₁C = 144;  // A71 in SP 98.
  const uint_fast8_t f_M₃  = 149;
  const uint_fast8_t f_M₁  = 206;  // Unique U.S. definition.
  const uint_fast8_t f_L₂  = 215;
  const uint_fast8_t f_K₁  = 227;
  const uint_fast8_t f_K₂  = 235;

  // Compound constituents may be specified using coefficients of the
  // following popular constituents.

  const uint_fast8_t O₁_index = 0;
  const uint_fast8_t K₁_index = 1;
  const uint_fast8_t P₁_index = 2;
  const uint_fast8_t M₂_index = 3;
  const uint_fast8_t S₂_index = 4;
  const uint_fast8_t N₂_index = 5;
  const uint_fast8_t L₂_index = 6;
  const uint_fast8_t K₂_index = 7;
  const uint_fast8_t Q₁_index = 8;
  const uint_fast8_t ν₂_index = 9;
  const uint_fast8_t S₁_index = 10;
  const uint_fast8_t M₁_DUTCH_index = 11;
  const uint_fast8_t λ₂_index = 12;
  const uint_fast8_t numCompoundBases = 13;

  typedef uint16_t year_t;  // Range 1..4000

  // Satellites for Doodson-style constituents get specified with a
  // vector of Satellite structs.  Congen does not perform adjustments
  // for latitude-dependent satellites, so the client should either
  // apply the latitude adjustment to r before passing it to Congen or
  // omit those satellites entirely.

  // It appears to be traditional Doodson usage to work with the
  // negative of N instead of N itself.  Those deltas need to be
  // negated before passing to Congen.

  struct Satellite {
    double r;    // amplitude ratio vs. the main constituent
    double Δp;   // delta of coefficient of p vs. the main constituent
    double ΔN;   // delta of coefficient of N vs. the main constituent
    double Δp₁;  // delta of coefficient of p₁ vs. the main constituent
    double α;    // phase correction vs. the main constituent (degrees)

  struct Constituent {

    // A Constituent is maintained using valarrays of equilibrium
    // arguments and node factors for the range of years that was
    // specified at the time of its construction.  Operations such as
    // adding two Constituents together are performed as vector
    // operations on those valarrays.  The client program must ensure
    // that every Constituent participating in such operations was
    // created using the same range of years; otherwise, garbage in,
    // garbage out.

    // libcongen will never normalize the (Vₒ+u) attribute of any
    // Constituent to [0, 360) or (-180, 180] degrees.  The helper
    // functions normalize, snormalize, and makeArrays supply
    // normalized results as needed for export.

    std::string name;
    double speed;                     // degrees per hour
    std::valarray<double> (Vₒ+u);     // degrees; [0] is firstYear
    std::valarray<double> f;          // [0] is firstYear

    // Default constructed constituents are only placeholders.  Speed
    // is initalized to zero; (Vₒ+u) and f are left in their default
    // initial states (zero length).  Name is set to "default".
    Constituent ();

    // This constructor is semantically equivalent to default
    // construction followed by resize(numYears).
    explicit Constituent (uint16_t numYears);

    // Copy constructor is semantically equivalent to default
    // construction followed by assignment.
    Constituent (const Constituent &x);

    // Constructor for Darwin-style constituents.

    // Coefficients for terms of V and u are passed in valarrays of
    // size 6 and 7, respectively.  See above for the ordering.

    // f_formula is the node factor formula number from SP 98; see above.

    // Speed is calculated as of the beginning of the year specified
    // by epochForSpeed.  For NOS compatibility, use 1900.

    Constituent (const std::string &name_in,
		 const std::valarray<double> &V_coefficients,
		 const std::valarray<double> &u_coefficients,
		 uint_fast8_t f_formula,
		 year_t firstYear,
		 year_t lastYear,
                 year_t epochForSpeed);

    // Constructor for Doodson-style constituents.  A vector of
    // Satellites substitutes for the u coefficients and node factor
    // formula, but otherwise the operation is exactly the same as the
    // previous constructor.

    Constituent (const std::string &name_in,
		 const std::valarray<double> &V_coefficients,
		 const std::vector<Satellite> &satellites,
		 year_t firstYear,
		 year_t lastYear,
                 year_t epochForSpeed);

    // The other constructors and the operations on Constituent are
    // sufficient to enable clients to define base constituents and
    // compute compound constituents from them.  However, constructing
    // the same base constituents in every client is painful and
    // redundant.  This constructor allows compound constituents to be
    // specified in terms of the 13 "popular" constituents listed
    // above.  Their definitions are compiled in, so each of those 13
    // may be specified in lazy fashion using a coefficient of 1 for
    // itself.

    Constituent (const std::string &name_in,
		 const std::valarray<double> &coefficients,
		 year_t firstYear,
		 year_t lastYear,
                 year_t epochForSpeed);

    // Assignment duplicates every field regardless of current state.
    Constituent &operator=  (const Constituent &x);

    // Add constituents together to produce a compound constituent.
    // The name is reset to "nameless".
    Constituent &operator+= (const Constituent &x);

    // Multiply by a scalar.
    // The name is reset to "nameless".
    Constituent &operator*= (double x);

    // Resize destroys any data currently stored in the Constituent.
    //   - (Vₒ+u) is set to the specified size and filled with zeroes.
    //   - f is set to the specified size and filled with ones.
    //   - speed is set to zero.
    //   - name is set to "zero".
    void resize (uint16_t numYears);


  // Reminder to self:  The C++ argument-dependent name lookup rules
  // conveniently cause these operations to work without a Congen::
  // qualifier by virtue of the fact that a Constituent is involved.

  Constituent operator+ (const Constituent &x, const Constituent &y);
  Constituent operator* (double x, const Constituent &y);
  Constituent operator* (const Constituent &y, double x);

  // As a sanity check, generate tables from SP 98 and send to stdout.
  // Make sure your terminal is using a monospace font, UTF-8, and at
  // least 87 characters wide.  LANG=en_US.utf8 xterm -fn
  //   -misc-fixed-medium-r-normal--20-200-75-75-c-100-iso10646-1 &
  void tables();

  // ----- Helper functions -----

  // Normalize an angle to [0, 360) degrees, round to the specified
  // number of decimals (> 0, ≤ 20), and return as a string of length
  // decimals+4.  Halfway cases are rounded according to the
  // implementation-dependent behavior of the %f printf format, which
  // for GNU/Linux is the weird "nearest even number" rule.
  std::string normalize (double degrees, uint_fast8_t decimals);

  // Normalize an angle to (-180, 180] degrees, round to the specified
  // number of decimals (> 0, ≤ 20), and return as a string of length
  // decimals+5.  Halfway cases are rounded according to the
  // implementation-dependent behavior of the %f printf format, which
  // for GNU/Linux is the weird "nearest even number" rule.
  std::string snormalize (double degrees, uint_fast8_t decimals);

  // Parse the text format used for congen_input.txt in Congen 1.5 and
  // earlier and append any Constituents specified therein to
  // constituents.

  // A Darwin-style constituent is specified as a line of this form:
  //   name  Basic   T s h p p₁ c   ξ ν ν′ 2ν″ Q R   f_formula

  // A Doodson-style constituent begins with a line of this form:
  //   name  Doodson   T s h p p₁ c   #-satellites
  // followed by additional lines that specify the indicated number of
  // satellites, possibly more than one per line.  Each satellite has
  // the following form, which is identical to the format used in the
  // IOS tidal package's tide3.dat file:
  //       Δp -ΔN Δp₁ α r[RN]
  // (A capital letter 'R' and a digit are optionally suffixed to r.)
  // Please note:
  //   Consistent with IOS, ΔN is NEGATED.
  //   Consistent with IOS, α is specified in ROTATIONS, not degrees.
  //   Latitude-dependent satellites (those with RN) are IGNORED.

  // A compound constituent is specified as a line of this form:
  //   name  Compound   O₁ K₁ P₁ M₂ S₂ N₂ L₂ K₂ Q₁ ν₂ S₁ M₁-DUTCH λ₂
  // For compound constituents, trailing zeroes may be omitted.

  // Comment lines contain # in the first column.

  // Returns 0 if input parsed OK, line number containing error if not.

  unsigned parseLegacyInput (std::istream &is,
			     year_t firstYear,
			     year_t lastYear,
			     year_t epochForSpeed,
			     std::vector<Constituent> &constituents);

  // Convert a vector of Constituents into the plain C arrays that are
  // expected by libtcd's create_tide_db function.  It is the client's
  // responsibility to deallocate the memory using delete [] (or
  // choose not to).  Equilibrium arguments are normalized to [0.00,
  // 359.99] as required by libtcd.

  void makeArrays (const std::vector<Constituent> &constituents,
		   char **&names,
		   double *&speeds,
		   float **&equilibriumArgs,
		   float **&nodeFactors);