JEphem Informatic Trail JEphem source code VSOP87.java
//*********************************************************************************
// class jephem.astro.planets.vsop87.VSOP87
// Software released under the General Public License (version 2 or later), available at
// http://www.gnu.org/copyleft/gpl.html
//*********************************************************************************

package jephem.astro.solarsystem.vsop87;

import jephem.astro.solarsystem.PlanetaryTheory;
import jephem.astro.solarsystem.SolarSystemConstants;
import jephem.astro.solarsystem.ComputationException;
import jephem.astro.Body;
import jephem.astro.AstroException;
import jephem.astro.spacetime.SpaceConstants;
import jephem.astro.spacetime.Time;
import jephem.astro.spacetime.TimeConstants;
import jephem.astro.spacetime.UnitsConstants;

import tig.maths.Maths;
import tig.maths.Vector3;
import tig.GeneralConstants;

import java.lang.reflect.Method;
import java.io.*;
/******************************************************************************
Low-level class performing calculations of planet coordinates from VSOP87 theory,

@author Thierry Graff.
@history nov 26 2000 : Creation
@history dec 17 2000 : Integration to the rest of JEphem.
@history aug 15 2001 : In calcCoord(), replaced 'filter' parameter by 'precision' ; implemented Exception mechanism
@history nov 14 2001 : Adapted the code to the evolution of other parts of the API.
@history feb 10 2002 : Adapted to load the data from binary files.
                       Introduced _dataFull and _dataJEphem (to load data only when necessary);

@todo internationalize error messages.
@todo handle correctly precision.
@todo see why XEphem coefs are different from BDL;
@todo in calcCoord : check date and precision validity
*********************************************************************************/
public abstract class VSOP87
    implements PlanetaryTheory, SolarSystemConstants, GeneralConstants, TimeConstants{

  //=================================================================================
  //                   STATIC VARIABLES
  //=================================================================================

  private static final int NB_PLANETS = 8;
  // _data[iBody][iTerm][iABC]
  /** Contains the terms of full precision version. */
  private static double[][][] _dataFull = new double[NB_PLANETS][][];
  /** Contains the terms of full version truncated using JEphem truncation. */
  private static double[][][] _dataJEphem = new double[NB_PLANETS][][];

  /** Path to the directory where VSOP87 files are stored (for full precision version). */
  private static String _dataPath;

  //=================================================================================
  //                     CONSTANTS
  //=================================================================================
  // For constants nbTerms and totalTerms, see at the end of class

  /** VSOP Version used to do the computation */
  private static final String VSOP87_VERSION = "A";

  /** Strings used to build the filenames and class to retrieve the data */
  private static final String FILENAME_PREFIX = "DataVSOP87" + VSOP87_VERSION + "_Full_";
  private static final String CLASSNAME_PREFIX = "jephem.astro.planets.vsop87.DataVSOP87"
                                               + VSOP87_VERSION + "_JEphem_";
  private static final String STR_DATA = "data";

  /** Max degree of time */
  private static final int ALPHA_MAX	= 5;
  /** Nb of coordinates */
  private static final int NB_COORDS	= 3;

  /** Semimajor axes; only for VSOP87 planets; for precision control.
      Positions in the array must correspond to SolarSystemConstants constants. */
  private static final double[] a0 = {0.0, 0.0, 0.3871, 0.7233, 1.0000, 1.5237,
                                5.2026, 9.5547, 19.2181, 30.1096};

  /** Planet names used for naming data files or classes ; these names are local to VSOP87.
  <BR>Their range in the array correspond to constants of
  <CODE>jephem.astro.SolarSystemConstants</CODE>.*/
  private static final String[] planetNames	= {"", "", "Mercury", "Venus", "Earth", "Mars",
                                               "Jupiter", "Saturn", "Uranus", "Neptune"};

  //*************************************************************
  // ******* Variables for precision and validity dates. ********
  //*************************************************************

  /** Represents the best precision the truncated version can reach. */
  private static final double LIMIT_TRUNCATED_PRECISION = 4.0;

  /** Initial time of validity */
  private static final double[] ti;
  /** Final time of validity */
  private static final double[] tf;
  /** Precision at final time of validity */
  private static final double[] pNow;
  /** Precision at initial and final times of validity */
  private static final double pLim = 1; // arc second
  static{
    // took J2000 and 365250 days per millenium
    double M6000 = 260045; // J2000 - 6000 years
    double M4000 = 990545; // J2000 - 4000 years
    double M2000 = 1721045; // J2000 - 2000 years
    double P2000 = 3182045; // J2000 + 2000 years
    double P4000 = 3912545; // J2000 + 4000 years
    double P6000 = 4643045; // J2000 + 6000 years
    // for all the arrays, the two first values are meaningless
    // this permits to use them with SolarSystem constants.
    ti = new double[]{0, 0, M4000, M4000, M4000, M4000, M2000, M2000, M6000, M6000};
    tf = new double[]{0, 0, P4000, P4000, P4000, P4000, P2000, P2000, P6000, P6000};
    pNow = new double[]{0, 0, 0.001, 0.006, 0.005, 0.023, 0.020, 0.100, 0.016, 0.030};
  };

  //=================================================================================
  //                                 PUBLIC METHODS
  //=================================================================================

  //******************* getPrecision *************
  // (implementation of PlanetaryTheory)
  /** Returns the precision of VSOP87 for a given body and a given julian day. */
  public static double getPrecision(int bodyIndex, double jd){
    if (bodyIndex < MERCURY || bodyIndex > NEPTUNE)
      throw new IllegalArgumentException("'bodyIndex' not valid - doesn't represent a body computed by VSOP87");

    if(jd >= JD1900 || jd <= JD2100)
      return pNow[bodyIndex];
    else if(jd < JD1900){
      double a = (pLim - pNow[bodyIndex]) / (ti[bodyIndex] - JD1900);
      return a * jd + pLim - a*ti[bodyIndex];
    }
    else{ // jd > 2100
      double a = (pLim - pNow[bodyIndex]) / (tf[bodyIndex] - JD2100);
      return a * jd + pLim - a*tf[bodyIndex];
    }
  }// end getPrecision

  //******************* setDataPath *************
  /** Sets the path where the VSOP87 files are located ; MUST be called before using VSOP87. */
  public static void setDataPath(String dataPath){
    _dataPath = dataPath;
  }// end setDataPath

  //******************* calcCoord(jd, body, precision, velocities) *************
  // This method only know how to perform the summations.
  // Choice of data depending on the precision, loading of data are handled by getData() end getNbTerms()
  /** Calculation of planetary positions, from Mercury to Neptune, using VSOP87 theory.
  <BR>If the date asked for the computation is not handled by VSOP87 theory or if parameter 'precision' can't be
  handled by the theory, a {@link jephem.astro.Body#setComputationException(ComputationException)} is called.
  @param jd julian date ; time scale : dynamical time TDB.
  @param body Index of the planet to calculate (using constants of
  {@link jephem.astro.solarsystem.SolarSystemConstants}).
  @param prec Precision required for calculations (in arc seconds).
  @param velocities Flag indicating if velocities must be also calculated.

  @throws AstroException if an IO or reflection error occurs.
  */
  public static void calcCoord(double jd, Body body, double precision, boolean velocities) throws AstroException{
    try{
      // Check validity of parameters
      int iBody = body.getIndex();
      if (iBody < MERCURY || iBody > NEPTUNE)
        throw new IllegalArgumentException("'body' not valid : '" + iBody + "' doesn't represent a body computed by VSOP87");
      if(precision < 0)
        throw new IllegalArgumentException("'precision' must be positive");

      // ComputationException checking
      if(precision < getPrecision(iBody, jd))
        body.setComputationException(new ComputationException(ComputationException.PRECISION_ERROR,
                                                              iBody, jd, TimeConstants.TT_TDB));
      if(jd > JD2000 + getValidityInterval(iBody) || jd < JD2000 - getValidityInterval(iBody))
        body.setComputationException(new ComputationException(ComputationException.DATE_LIMIT_ERROR,
                                                              iBody, jd, TimeConstants.TT_TDB));

      int[][] nbTerms = null;
      double[][] data = null;
      data = getData(iBody, precision);
      nbTerms = getNbTerms(iBody, precision);

      int coord, alpha, n;                            // index variables
      int nbegin, nend;                               // limit of index n
      double p, q;                                    // for precision control
      double prec = 0.0;
      double a=0, b=0, c=0, arg, term, termdot;       // for term calculation
      double res[] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0};  // to stock results
      int i;                                          // used sevral times.

      // Time and its powers ; t[i] = t^i
      double[] t = new double[ALPHA_MAX + 1];
      double[] t_abs = new double[ALPHA_MAX + 1];
      t[0] = 1;
      t[1] = (jd - JD2000)/DAYS_PER_MILLENIUM;
      for (i = 2; i<= ALPHA_MAX; i++) t[i] = t[i-1] * t[1];
      for (i = 0; i <= ALPHA_MAX; ++i) t_abs[i] = Math.abs(t[i]);

      // Precision control
      q = Math.max(3.0, -Maths.log10(prec + 1e-50));
      //System.out.println("q = " + q);

      // *** Term summation ***
      nbegin = 0;
      // for each coordinate (order : X, Y, Z)
      for (coord = 0; coord < NB_COORDS; coord++){

        // sum on the powers of time
        for (alpha = 0 ; alpha <= ALPHA_MAX ; alpha ++){

            // reduce threshold progressively for higher precision
            p = a0[iBody] * prec / 10.0 / (q-2) /
                (t_abs[alpha] + alpha * (alpha>0 ? t_abs[alpha-1] : 0.0) * 1e-4 + 1e-50);
            //System.out.println("p = " + p);

          nend = nbTerms[coord][alpha];
          if (nend == 0)
            continue; // no term for this couple (coord, alpha)
          nend += nbegin;

          term = termdot = 0.0;

          // summation of terms for this couple (coord, alpha)
          for (n = nbegin ; n < nend ; n++){
            a = data[n][0];
            if (a < p) continue;	// ignore small terms
            b = data[n][1];
            c = data[n][2];
            arg = b + c * t[1];
            term += a * Math.cos(arg);
            if (velocities) termdot += -c * a * Math.sin(arg);
          }// end for n

          res[coord] += t[alpha] * term;
          if (velocities) res[coord + NB_COORDS] += t[alpha] * termdot
                          + ((alpha > 0) ? alpha * t[alpha - 1] * term : 0.0);
          nbegin = nend;

        }// end for alpha

      }// end for coord

      // convert millenium rate to day rate
      if (velocities)
        for (i = 3; i < 6; ++i) res[i] /= DAYS_PER_MILLENIUM;

      // Set the fields of Body
      body.setFrame(SpaceConstants.FRAME_THEORY);
      body.setCoordinateExpression(SpaceConstants.CARTESIAN);
      body.setPositionUnits(UnitsConstants.UNITGROUP_AU_AU_AU);
      if (velocities) body.setVelocityUnits(UnitsConstants.UNITGROUP_AUD_AUD_AUD);
      body.setPositionCoords(res[0], res[1], res[2]);
      if (velocities) body.setVelocityCoords(res[3], res[4], res[5]);
    }
    catch(Exception e){
      throw new AstroException(e);
    }
  }// end calcCoord(jd, body, precision, velocities)

  //=================================================================================
  //                                 PRIVATE METHODS
  //=================================================================================

  //******************* getValidityInterval *************
  /** Returns the validity interval of VSOP87 for a given body, in days.
  The theory is valid from J2000 - interval to J2000 + interval
  */
  private static double getValidityInterval(int bodyIndex){
    switch(bodyIndex){
      case MERCURY:
      case VENUS:
      case EARTH:
      case MARS:
        return 1461000;// 4000 years
      case JUPITER:
      case SATURN:
        return 730500; // 2000 years
      case URANUS:
      case NEPTUNE:
        return 2191500; // 6000 years
      default:
        throw new IllegalArgumentException("'bodyIndex' not valid - doesn't represent a body computed by VSOP87");
    }
  }// end getValidityInterval

  //******************* getData(iBody, precision) *************
  // Choice of data depending on the precision, loading of data are handled by getData() end getNbTerms()
  /** Returns the terms for a given body ;
  param iBody Index of a body, using <CODE>jephem.astro.SolarSystemConstants</CODE> constants.
  */
  private static double[][] getData(int iBody, double precision) throws AstroException{
    if (_dataPath == null){
      throw new AstroException("Before Using VSOP87, you must indicate where VSOP data are located with setDataPath()");
    }

    double[][] data; // returned value

    try{
      // 1 - Not in applet context - data retrieved from binary files
      if(precision < LIMIT_TRUNCATED_PRECISION){
        if (_dataFull[iBody - MERCURY] == null){
          // need to fill data from binary files
          ObjectInputStream  ois = new ObjectInputStream(new FileInputStream(_dataPath + FILENAME_PREFIX + planetNames[iBody]));
          int totalTerms = getTotalTerms(iBody, precision);
          data = new double[totalTerms][3]; // 3 for A, B and C
          for (int i = 0; i < totalTerms; i++){
            data[i][0] = ois.readDouble();
            data[i][1] = ois.readDouble();
            data[i][2] = ois.readDouble();
          }// end for
          ois.close();
          // put data in static variable
          _dataFull[iBody - MERCURY] = data;
        }
        else{ // Already loaded, nothing to do
          data = _dataFull[iBody - MERCURY];
        }
      }// end if(precision < LIMIT_TRUNCATED_PRECISION)
      else{
        if (_dataJEphem[iBody - MERCURY] == null){
          // data[][] must be retrieved from a data class
          // Use reflection, as we don't know which data class to use until execution time.
          Class dataClass = Class.forName(CLASSNAME_PREFIX + planetNames[iBody]);
          data = (double[][])dataClass.getDeclaredField(STR_DATA).get(null);
        // put data in static variable
        _dataJEphem[iBody - MERCURY] = data;
        }
        else{ // Already loaded, nothing to do
          data = _dataJEphem[iBody - MERCURY];
        }
      }
      // return results
      return data;
    }
    catch (Exception e) {
      throw new AstroException(e);
    }

  }// end getData

  //**********************************************************
  /** Returns the nb of terms for a planet - use : nbTerms[iCoord][iAlpha]*/
  private static int[][] getNbTerms(int iBody, double precision) throws Exception{
    // reflection not used as slower (17%)
    if(precision < LIMIT_TRUNCATED_PRECISION){
      switch(iBody){
        case MERCURY : return nbTerms_A_Full_Mercury;
        case VENUS : return nbTerms_A_Full_Venus;
        case EARTH : return nbTerms_A_Full_Earth;
        case MARS : return nbTerms_A_Full_Mars;
        case JUPITER : return nbTerms_A_Full_Jupiter;
        case SATURN : return nbTerms_A_Full_Saturn;
        case URANUS : return nbTerms_A_Full_Uranus;
        case NEPTUNE : return nbTerms_A_Full_Neptune;
      }
    }
    else{
      switch(iBody){
        case MERCURY : return nbTerms_A_JEphem_Mercury;
        case VENUS : return nbTerms_A_JEphem_Venus;
        case EARTH : return nbTerms_A_JEphem_Earth;
        case MARS : return nbTerms_A_JEphem_Mars;
        case JUPITER : return nbTerms_A_JEphem_Jupiter;
        case SATURN : return nbTerms_A_JEphem_Saturn;
        case URANUS : return nbTerms_A_JEphem_Uranus;
        case NEPTUNE : return nbTerms_A_JEphem_Neptune;
      }
    }
    throw new IllegalArgumentException("'body' parameter incorrect");
  }// end getNbTerms

  //**********************************************************
  /** Returns the total nb of terms for a planet. */
  private static int getTotalTerms(int iBody, double precision) throws Exception{
    // reflection not used as slower (17%)
    if(precision < LIMIT_TRUNCATED_PRECISION){
      switch(iBody){
        case MERCURY : return totalTerms_A_Full_Mercury;
        case VENUS : return totalTerms_A_Full_Venus;
        case EARTH : return totalTerms_A_Full_Earth;
        case MARS : return totalTerms_A_Full_Mars;
        case JUPITER : return totalTerms_A_Full_Jupiter;
        case SATURN : return totalTerms_A_Full_Saturn;
        case URANUS : return totalTerms_A_Full_Uranus;
        case NEPTUNE : return totalTerms_A_Full_Neptune;
      }
    }
    else{
      switch(iBody){
        case MERCURY : return totalTerms_A_JEphem_Mercury;
        case VENUS : return totalTerms_A_JEphem_Venus;
        case EARTH : return totalTerms_A_JEphem_Earth;
        case MARS : return totalTerms_A_JEphem_Mars;
        case JUPITER : return totalTerms_A_JEphem_Jupiter;
        case SATURN : return totalTerms_A_JEphem_Saturn;
        case URANUS : return totalTerms_A_JEphem_Uranus;
        case NEPTUNE : return totalTerms_A_JEphem_Neptune;
      }
    }
    throw new IllegalArgumentException("'body' parameter incorrect");
  }// end getTotalTerms

  //=================================================================================
  //                     CONSTANTS
  //=================================================================================


  //***********************************************************************
  // Number of terms for version A, when generated with "JEPHEM" filter
  //***********************************************************************
  private static final int[][] nbTerms_A_JEphem_Mercury = {
    {25, 11, 8, 7, 1, 0},
    {26, 11, 8, 7, 2, 0},
    {10, 7, 6, 4, 2, 0}
  }; // end nbTerms_A_JEphem_Mercury[][]
  private static final int totalTerms_A_JEphem_Mercury = 135;

  private static final int[][] nbTerms_A_JEphem_Venus = {
    {85, 40, 22, 4, 4, 3},
    {84, 42, 22, 4, 4, 3},
    {30, 17, 6, 3, 3, 2}
  }; // end nbTerms_A_JEphem_Venus[][]
  private static final int totalTerms_A_JEphem_Venus = 378;

  private static final int[][] nbTerms_A_JEphem_Earth = {
    {57, 28, 13, 6, 4, 1},
    {57, 28, 13, 6, 4, 1},
    {3, 4, 3, 3, 2, 1}
  }; // end nbTerms_A_JEphem_Earth[][]
  private static final int totalTerms_A_JEphem_Earth = 234;

  private static final int[][] nbTerms_A_JEphem_Mars = {
    {128, 106, 65, 32, 16, 11},
    {128, 106, 65, 32, 16, 11},
    {34, 20, 10, 7, 4, 2}
  }; // end nbTerms_A_JEphem_Mars[][]
  private static final int totalTerms_A_JEphem_Mars = 793;

  private static final int[][] nbTerms_A_JEphem_Jupiter = {
    {57, 48, 32, 13, 8, 4},
    {56, 47, 32, 13, 8, 4},
    {18, 11, 6, 3, 0, 0}
  }; // end nbTerms_A_JEphem_Jupiter[][]
  private static final int totalTerms_A_JEphem_Jupiter = 360;

  private static final int[][] nbTerms_A_JEphem_Saturn = {
    {77, 48, 30, 21, 13, 5},
    {76, 49, 31, 21, 13, 5},
    {23, 15, 11, 8, 3, 0}
  }; // end nbTerms_A_JEphem_Saturn[][]
  private static final int totalTerms_A_JEphem_Saturn = 449;

  private static final int[][] nbTerms_A_JEphem_Uranus = {
    {75, 58, 41, 28, 8, 0},
    {73, 59, 40, 31, 8, 0},
    {19, 11, 5, 2, 0, 0}
  }; // end nbTerms_A_JEphem_Uranus[][]
  private static final int totalTerms_A_JEphem_Uranus = 458;

  private static final int[][] nbTerms_A_JEphem_Neptune = {
    {25, 15, 8, 5, 0, 0},
    {23, 15, 7, 5, 0, 0},
    {13, 5, 1, 1, 0, 0}
  }; // end nbTerms_A_JEphem_Neptune[][]
  private static final int totalTerms_A_JEphem_Neptune = 123;

  //***********************************************************************
  // Number of terms for version A, when generated with "FULL" filter
  //***********************************************************************

  private static final int[][] nbTerms_A_Full_Mercury = {
    {1449, 792, 299, 54, 15, 10},
    {1438, 782, 299, 59, 15, 10},
    {598, 351, 143, 28, 10, 7}
  }; // end nbTerms_A_Full_Mercury[][]
  private static final int totalTerms_A_Full_Mercury = 6359;

  private static final int[][] nbTerms_A_Full_Venus = {
    {548, 338, 99, 5, 4, 3},
    {565, 325, 99, 5, 4, 3},
    {190, 108, 45, 10, 3, 3}
  }; // end nbTerms_A_Full_Venus[][]
  private static final int totalTerms_A_Full_Venus = 2357;

  private static final int[][] nbTerms_A_Full_Earth = {
    {843, 491, 204, 18, 15, 6},
    {854, 496, 202, 17, 15, 6},
    {178, 120, 53, 12, 6, 2}
  }; // end nbTerms_A_Full_Earth[][]
  private static final int totalTerms_A_Full_Earth = 3538;

  private static final int[][] nbTerms_A_Full_Mars = {
    {1584, 956, 387, 135, 41, 21},
    {1612, 969, 384, 136, 44, 21},
    {355, 232, 122, 51, 16, 7}
  }; // end nbTerms_A_Full_Mars[][]
  private static final int totalTerms_A_Full_Mars = 7073;

  private static final int[][] nbTerms_A_Full_Jupiter = {
    {1055, 488, 255, 140, 58, 11},
    {1037, 499, 259, 136, 60, 11},
    {216, 104, 65, 27, 10, 3}
  }; // end nbTerms_A_Full_Jupiter[][]
  private static final int totalTerms_A_Full_Jupiter = 4434;

  private static final int[][] nbTerms_A_Full_Saturn = {
    {1652, 892, 481, 215, 87, 31},
    {1658, 917, 465, 201, 88, 32},
    {420, 217, 87, 44, 19, 6}
  }; // end nbTerms_A_Full_Saturn[][]
 private static final int totalTerms_A_Full_Saturn = 7512;

  private static final int[][] nbTerms_A_Full_Uranus = {
    {1464, 649, 249, 84, 12, 0},
    {1447, 659, 255, 80, 12, 0},
    {235, 98, 33, 12, 0, 0}
  }; // end nbTerms_A_Full_Uranus[][]
  private static final int totalTerms_A_Full_Uranus = 5289;

  private static final int[][] nbTerms_A_Full_Neptune = {
    {772, 330, 102, 33, 7, 0},
    {746, 325, 97, 34, 7, 0},
    {133, 37, 11, 2, 0, 0}
  }; // end nbTerms_A_Full_Neptune[][]
  private static final int totalTerms_A_Full_Neptune = 2636;

  //=================================================================================
  //                     BUILD (to comment after execution)
  //=================================================================================
/*
  public static void main(String[] args){
    buildTotalNbTerms();
  }// end main

  private static void buildTotalNbTerms(){
    try{
      int[][] nbTerms;
      int total, i, j, k;

      String[] strFilter = {"Full", "JEphem"};
      for(k = 0; k < 2; k++){
        System.out.println(strFilter[k]);
        Class dataClass = Class.forName("jephem.astro.planets.vsop87.VSOP87");
        for (int iBody = 0; iBody < NB_PLANETS; iBody++){
          nbTerms = (int[][])dataClass.getDeclaredField("nbTerms_A_" + strFilter[k] + "_" + planetNames[iBody + MERCURY]).get(null);
          total = 0;
          for(i = 0; i < NB_COORDS; i++){
            for(j = 0; j <= ALPHA_MAX; j++){
              total += nbTerms[i][j];
            }
          }
          System.out.println("  private static final int totalTerms_A_" + strFilter[k] + "_" + planetNames[iBody + MERCURY] + " = " + total + ";");
        } // iBody
      } // k
    }
    catch(Exception e){
      e.printStackTrace();
    }
  } // end buildTotalNbTerms
*/
} //end class VSOP87