JEphem Informatic Trail build classes source code TestWithSwiss.java
//*********************************************************************************
// class TestWithSwiss // default package
// Software released under the General Public License (version 2 or later), available at
// http://www.gnu.org/copyleft/gpl.html
//*********************************************************************************

import jephem.astro.*;
import jephem.astro.spacetime.*;
import jephem.astro.planets.vsop87.VSOP87;
import jephem.util.JEphemException;
import jephem.astro.AstroException;

import tig.*;

import java.io.*;
import java.util.Date;
import java.text.NumberFormat;

/******************************************************************
Contains methods used to compare the results of JEphem and SwissEphemeris.
<BR>JEphem program doesn't use it ; it was used to test JEphem.
<BR>This class is dependant from <CODE>jephem</CODE>, <CODE>tig</CODE> and <CODE>swisseph</CODE> packages, which must be in the classpath.
<BR>In this class, when arrays characterizing planetary data were needed (ex : for "planetNames"), the indexes of <CODE>jephem.astro.SolarSystemConstants</CODE> were used.
@author Thierry Graff
@history feb 17 2002 creation from TestVSOP87.java
@todo
*****************************************************************/
public abstract class TestWithSwiss 
  implements GeneralConstants, SolarSystemConstants, SpaceConstants, TimeConstants, UnitsConstants{

  private static final int NB_COORD = 6;

  /** Differences between the values of JEphem and BDL [planet][jd][coord]. */
  private static double[][][] delta;

  /** To know computation time. */
  private static Date startCalcTime, endCalcTime;

  // planet names ; indexes correspond to jephem.astro.ISolarSystem.java constants.
  private static String[] planetNames	= {"", "", "Mercury", "Venus", "Earth", "Mars", "Jupiter",
                                  "Saturn", "Uranus", "Neptune"};

  // Absolute path to the files containing terms for summation.
  private final static String JEPHEM_DATA_PATH = "d:\\b_dvpt\\jephem\\data\\astro\\planets\\vsop87\\VSOP87A\\";
  private final static String SWISS_DATA_PATH = "d:\\b_dvpt\\jephem\\data\\astro\\swissEphem\\";
  private final static String DELTA = "<FONT FACE=\"SYMBOL\">D</FONT>";
  private final static String XXX = "xxx";

  private static final NumberFormat NF = NumberFormat.getNumberInstance();
  private static final NumberFormat NF2 = NumberFormat.getNumberInstance();
  static{
    NF.setMaximumFractionDigits(10);
    NF.setMinimumFractionDigits(10);
    NF2.setMaximumFractionDigits(3);
    NF2.setMinimumFractionDigits(3);
  };


  //********************************************************************************
  //******************************** main() ****************************************
  //********************************************************************************
  public static void main(String[] args){
//    for (int i = 0; i < args.length; i++) System.out.println(args[i]);
    if(args[0].equalsIgnoreCase("compareWithSwiss")){
      System.out.println("calling compareWithSwiss");
      compareWithSwiss(args[1], args[2], args[3], args[4], args[5], args[6], args[7]);
    }
  }// end main

  //**************** compareWithSwiss() ****************************************
  /**
  WARNING : 'strUnits' must have a length of 6.
  */
  public static void compareWithSwiss(String strJd, 
                                      String strBodies,
                                      String strFrame, 
                                      String strTimeFrame, 
                                      String strSphereCart, 
                                      String strUnits, 
                                      String outputFile){

      int i, j;
      String strTmp;
      
    try{
      // *** 1 - Parameter parsing - no checking
      // all the 6 coords are imposed
      int[] coords = new int[]{COORD_X0, COORD_X1, COORD_X2, COORD_V0, COORD_V1, COORD_V2};
      int nbCoords = 6;

      double jd = Double.parseDouble(strJd);
      System.out.println("JD = " + jd);

      int[] bodyIndices = Strings.stringToIntArray(strBodies);
      int nbBodies = bodyIndices.length; 
      strTmp="Bodies : ";
      for(i=0; i<nbBodies; i++){
        strTmp += SolarSystem.getBodyLabel(bodyIndices[i]) + " ";
      }
      System.out.println(strTmp);

      int frame = Integer.parseInt(strFrame);
      System.out.println("frame : " + Space.getFrameLabel(frame));

      int timeFrame = Integer.parseInt(strTimeFrame);
      System.out.println("timeFrame : " + jephem.astro.spacetime.Time.getTimeFrameLabel(timeFrame));

      int sphereCart = Integer.parseInt(strSphereCart);
      System.out.println("sphereCart : " + Space.getCoordinateExpressionLabel(sphereCart));

      int[] units = Strings.stringToIntArray(strUnits);
      int nbUnits = units.length; 
      strTmp="Units : ";
      for(i=0; i<nbUnits; i++){
        strTmp += Units.getUnitLabel(units[i]) + " ";
      }
      System.out.println(strTmp);

      // *** 2 - Compute the values
      VSOP87.setDataPath(JEPHEM_DATA_PATH);
      SwissEphemeris.setDataPath(SWISS_DATA_PATH);
      double precision = 0.0;
      boolean velocities = Space.containsVelocityCoord(coords);
      double[][] jephemData = new double[nbBodies][nbCoords];
      double[][] swissData = new double[nbBodies][nbCoords];
      AstroContext ac;
      Body b;
      // compute JEphem
      ac = new AstroContext(jd, timeFrame, bodyIndices);
      ac.setAstroEngine(AstroPrefs.JEPHEM);
      try{
        ac.calcBodyCoords(frame, sphereCart, precision, velocities, units);
      }
      catch(AstroException ae){System.out.println(ae.toString());}
   
      for(i = 0; i < nbBodies; i++){
        b = ac.getBody(bodyIndices[i]);
        for(j = 0; j < nbCoords; j++){
          jephemData[i][j] = b.getCoord(j);
        }
      }
      // Compute SwissEphemeris
      ac = new AstroContext(jd, timeFrame, bodyIndices);
      ac.setAstroEngine(AstroPrefs.SWISS_EPHEMERIS);
      try{
        ac.calcBodyCoords(frame, sphereCart, precision, velocities, units);
      }
      catch(AstroException ae){System.out.println(ae.toString());}
      for(i = 0; i < nbBodies; i++){
        b = ac.getBody(bodyIndices[i]);
        for(j = 0; j < nbCoords; j++){
          swissData[i][j] = b.getCoord(j);
        }
      }

      // *** 3 - Generate the output
      FileOutputStream fos = new FileOutputStream(new File(outputFile));
      String strRes = BLANK;

      // 3.1 - Title
      strRes += buildOutputBeginning();
      strRes += "<BR><B>JD</B> = " + jd;
      strTmp = BLANK;
      for(i=0; i<nbBodies; i++){
        strTmp += SolarSystem.getBodyLabel(bodyIndices[i]) + " ";
      }
      strRes += "<BR><B>Bodies</B> : " + strTmp;
      strRes += "<BR><B>frame</B> : " + Space.getFrameLabel(frame);
      strRes += "<BR><B>timeFrame</B> : " + jephem.astro.spacetime.Time.getTimeFrameLabel(timeFrame);
      strRes += "<BR><B>sphereCart</B> : " + Space.getCoordinateExpressionLabel(sphereCart);
      strTmp="<B>Units</B> : ";
      for(i=0; i<nbUnits; i++){
        strTmp += Units.getUnitLabel(units[i]) + "   ";
      }
      strRes += "<BR>" + strTmp;

      // 3.2 - Result table
      // 3.2.1 - Headers
      strRes += "<BR><BR><TABLE CLASS=\"Small\" BORDER = \"1\">" + LS;
      strRes += "<TR><TD COLSPAN = \"2\"></TD>" + LS;
      strRes += "<TH>Swiss</TH><TH>JEphem</TH><TH>" + DELTA + "</TH>" 
             + "<TH>100 * " + DELTA + " / JEphem</TH>" + "<TH>100 * " + DELTA + " / Swiss</TH>" + LS;
      strRes += "</TR>";
      double delta, deltaJ, deltaS;
      String tdClass, tdClassSuffix = "2", tdTag;
      String strDouble;
      boolean testDeg; // do we display degrees ?
      String[] tmp = Space.getCoordGroupLabels(Space.getCoordGroup(frame));
      String[] coordLabels = new String[6];
      for(i=0; i<3; i++) coordLabels[i] = tmp[i];
      for(i=3; i<6; i++) coordLabels[i] = tmp[i-3] + "'";
      for(i=0; i<6; i++) coordLabels[i] += "&nbsp;&nbsp;(" + Units.getUnitLabel(units[i]) + ")";      for(i=0; i<6; i++) coordLabels[i] += "  (" + Units.getUnitLabel(units[i]) + ")";

      // Main loop
      for (int iBody = 0; iBody < nbBodies; iBody ++){
        //tdClassSuffix = ((iBody%2 == 0) ? "1" : "2");

        // write planet name
        strRes += "<TR><TD ROWSPAN = \"6\">" + SolarSystem.getBodyLabel(bodyIndices[iBody]) + "</TD>";
        for (int iCoord = 0; iCoord < 6; iCoord++){
          testDeg = (units[iCoord] == ANGULAR_UNIT_DEG || units[iCoord] == ANGULAR_SPEED_UNIT_DEG_PER_DAY || units[iCoord] == ANGULAR_SPEED_UNIT_DEG_PER_S);
          tdClass = ((iCoord < 3) ? "pos" : "vel") + tdClassSuffix;
          tdTag = "<TD CLASS = \"" + tdClass + "\">";

          // display label and unit
          strRes += tdTag + coordLabels[iCoord] + "</TD>";

          // display swiss eph data
          if (testDeg)
            strDouble = Formats.doubleToDMS(swissData[iBody][iCoord]);
          else
            strDouble = NF.format(swissData[iBody][iCoord]);
          strRes += tdTag + strDouble + "</TD>";

          // display jephem data
          if (testDeg)
            strDouble = Formats.doubleToDMS(jephemData[iBody][iCoord]);
          else
            strDouble = NF.format(jephemData[iBody][iCoord]);
          strRes += tdTag + strDouble + "</TD>";

          //display delta
          delta = Math.abs(swissData[iBody][iCoord] - jephemData[iBody][iCoord]);
          if (testDeg)
            strDouble = Formats.doubleToDMS(delta);
          else
            strDouble = NF.format(delta);
          strRes += tdTag + strDouble + "</TD>";

          //display 100 * delta / JEphem
          if (jephemData[iBody][iCoord] != 0){
            deltaJ = 100.0 * delta / jephemData[iBody][iCoord];
            strDouble = NF2.format(deltaJ);
          }
          else
            strDouble = XXX;
          strRes += tdTag + strDouble + "</TD>";

          //display 100 * delta / Swiss
          if (swissData[iBody][iCoord] != 0){
            deltaJ = 100.0 * delta / swissData[iBody][iCoord];
            strDouble = NF2.format(deltaJ);
          }
          else
            strDouble = XXX;
          strRes += tdTag + strDouble + "</TD>";

          strRes += "</TR>";

        }// end for iCoord
      }// end for iBody
      strRes += "<TABLE/>" + LS;

      // 3.3 - End
      strRes += "</BODY></HTML>";
      fos.write(strRes.getBytes());
      fos.close();
    }
    catch(Exception e){
      e.printStackTrace();
      System.exit(0);
    }

  }// end testJEphemPrecision

  //******************************************************************************
  //******************************** buildOutputBeginning() **************************
  //******************************************************************************
  private static String buildOutputBeginning(){

    String strRes = new String();

    strRes += "<HTML>" + LS;
    strRes += "<HEAD>" + LS;
    strRes += "<STYLE>" + LS;
    strRes += "TABLE {text-align:center}" + LS;
    strRes += "TABLE.Small {font-size:smaller}" + LS;
    strRes += "TD.pos1 {background-color : #DEEBF4;}" + LS;
    strRes += "TD.vel1 {background-color : #99CCFF;}" + LS;
    strRes += "TD.pos2 {background-color : #FFCC99;}" + LS;
    strRes += "TD.vel2 {background-color : #F9CF3F;}" + LS;
    strRes += "</STYLE>" + LS;
    strRes += "</HEAD>" + LS;
    strRes += "<BODY>" + LS;

    strRes += "<CENTER><H1>Tests JEphem - SwissEphemeris</H1>" +LS;
    strRes += "Generated by <CODE>TestWithSwiss.java</CODE></CENTER>" + LS;
    return strRes;
  }// end buildOutputBeginning()

  //******************************************************************************
  //******************************** buildOutputDetailByJD() *********************
  // Tables showing all the results, grouped by julian days
  //******************************************************************************
/*  private static void buildOutputDetailByJD(FileOutputStream fos){
    try{
      int iJD, iPlanet, iCoord; // main indexes
      String strRes = new String();

      strRes += LS + LS + "<!-- ************************************************** -->" + LS;
      strRes += "<H2><A NAME=\"DetailByJD\"></A>Detailed tables by JD</H2>" + LS;
      // for each julian day
      for (iJD = 0; iJD < _nbJD; iJD ++){
        // title
        strRes += "<TABLE CLASS=\"Small\" BORDER = \"1\">" + LS;
        strRes += "<TR><TD COLSPAN = \"2\"></TD>" + LS;
        strRes += "<TD COLSPAN = \"6\"><B>" + LS;
        strRes += "JD = " + Double.toString(bdlT[0][iJD].jd); // 0 because JDs are the same
        strRes += "    " + bdlT[0][iJD].strDate; //for different planets
        strRes += "</B></TD></TR>" + LS;
        // legends
        strRes += "<TR><TD COLSPAN = \"2\">Variables</TD>" + LS;
        for (int i=0; i < NB_COORD; i++){
          strRes += "<TD>" + coordNames[i] + " (" + strUnitPos[i] + ")"
                 + "</TD>" + LS;
        }
        strRes += "</TR>" + LS;

        //for each planet
        for (iPlanet = 0; iPlanet < _nbPlan; iPlanet ++){
//System.out.println(" - iJD = " + iJD + "iPlanet = " + iPlanet);
          // write planet name
          strRes += "<TR BGCOLOR =\"#FAD6AD\">";
          strRes += "<TD ROWSPAN = \"3\">";
          for (int i = 0; i < 3; i++)
            strRes += planetNames[bdlT[iPlanet][0].body.getIndex()].substring(i, i+1) + "<BR>";
          strRes += "</TD>" + LS;

          // *** BDL line
          strRes += "<TD>BDL</TD>" + LS;
          for (iCoord = 0; iCoord < NB_COORD; iCoord ++){
            strRes += "<TD>"
                   + NF.format(bdlT[iPlanet][iJD].body.getCoord(iCoord))
                   + "</TD>" + LS;
          } // end for iCoord
          strRes += "</TR>" + LS;

          // *** JEphem line
          strRes += "<TR>" + LS;
          strRes += "<TD>JEphem full</TD>" + LS;
          for (iCoord = 0; iCoord < NB_COORD; iCoord ++){
            strRes += "<TD>"
                   +  NF.format(jephemT[iPlanet][iJD].body.getCoord(iCoord))
                   + "</TD>" + LS;
          } // end for iCoord
          strRes += "</TR>" + LS;

          // *** Difference JEphem  - BDL line
          strRes += "<TR BGCOLOR =\"yellow\">" + LS;
          strRes += "<TD>|J. full - BDL|</TD>" + LS;
          for (iCoord = 0; iCoord < NB_COORD; iCoord ++){
            strRes += "<TD>"
                   + NF.format(delta[iPlanet][iJD][iCoord])
                   + "</TD>" + LS;
          } // end for iCoord
          strRes += "</TR>" + LS;

          fos.write(strRes.getBytes());
          strRes ="";
        } // end for each planet
        strRes += "<TABLE/><BR><BR>" + LS;
        strRes += "<BR>" + LS;

      } // end for each JD

      fos.write(strRes.getBytes());
    }
    catch(Exception e){
      System.out.println("Exception" + e.toString());
      e.printStackTrace();
      System.exit(0);
    }
  }// end buildOutputDetailByJD()
*/
  //******************************************************************************
  //******************************** buildOutputDiffByPlanet() *********************
  //******************************************************************************
/*
  private static void buildOutputDiffByPlanet(FileOutputStream fos){
    try{
      int iJD, iPlanet, iCoord; // main indexes
      String strRes = new String("");

      double[][] deltaMax = new double[_nbPlan][NB_COORD];

      strRes += LS + LS + "<!-- ************************************************** -->" + LS;
      strRes += "<H2><A NAME=\"DiffByPlan\"></A>Differences tables by planet</H2>" + LS;

      // for each planet
      for (iPlanet = 0; iPlanet < _nbPlan; iPlanet ++){
        // title
        strRes += "<TABLE CLASS=\"2\" BORDER = \"1\">" + LS;
        strRes += "<TR><TD></TD>" + LS;
        strRes += "<TD COLSPAN = \"8\"><B>" + LS;
        strRes += planetNames[bdlT[iPlanet][0].body.getIndex()];
        strRes += "</B></TD></TR>" + LS;

        // coordinate legends
        strRes += "<TR><TD>Dates</TD>" + LS;
        for (int i = 0; i < NB_COORD; i++){
          strRes += "<TD>Delta " + coordNames[i] + " (" + strUnitPos[i]
                 + ")</TD>" + LS;
          if (i == 2)
          strRes += "<TD>Norm (" + strUnitPos[0]
                 + ")</TD>" + LS;
        }
        strRes += "<TD>Norm (" + strUnitPos[NB_COORD - 1]
               + ")</TD>" + LS;

        strRes += "</TR>" + LS;

        //for each julian day, write a line of data
        for (iJD = 0; iJD < _nbJD; iJD ++){
//System.out.println(" - iJD = " + iJD + "iPlanet = " + iPlanet);
          strRes += "<TR>" + LS;
          // write date
          int idx = bdlT[iPlanet][iJD].strDate.indexOf(" ");
          strRes += "<TD>" + bdlT[iPlanet][iJD].strDate.substring(0, idx).trim()
          + "</TD>" + LS;

          // write data for delta pos and vel
          for (iCoord = 0; iCoord < NB_COORD; iCoord ++){
            strRes += "<TD>" + NF.format(delta[iPlanet][iJD][iCoord])
                   + "</TD>" + LS;
            if (delta[iPlanet][iJD][iCoord] > deltaMax[iPlanet][iCoord])
              deltaMax[iPlanet][iCoord] = delta[iPlanet][iJD][iCoord];
            if (iCoord==2)
              // write data for delta norm pos
              strRes += "<TD><B>" + NF.format(deltaNormPos[iPlanet][iJD])
                     + "</B></TD>" + LS;

          } // end for iCoord
          // write data for delta norm vel
          strRes += "<TD><B>" + NF.format(deltaNormVel[iPlanet][iJD])
                 + "</B></TD>" + LS;

          strRes += "</TR>" + LS;

          fos.write(strRes.getBytes());
          strRes ="";
        } // end for each JD

        // *** Delta max
        strRes += "<TR BGCOLOR = \"yellow\"><TD>Delta max</TD>" + LS;
        for (iCoord = 0; iCoord < NB_COORD; iCoord ++){
          strRes += "<TD>" + NF.format(deltaMax[iPlanet][iCoord])
                 + "</TD>" + LS;
          if (iCoord==2)
            // write data for delta max norm pos
            strRes += "<TD><B>" + NF.format(deltaMaxNormPos[iPlanet])
                   + "</B></TD>" + LS;
        } // end for iCoord
        // write data for delta max norm vel
        strRes += "<TD><B>" + NF.format(deltaMaxNormVel[iPlanet])
               + "</B></TD>" + LS;

        strRes += "</TR>" + LS;
        strRes += "<TABLE/>" + LS;

        // *** Delta max for the norms
        int indPla = bdlT[iPlanet][0].body.getIndex();
        strRes += "Delta max of vector norms (built only with delta max) for <B>"
               + planetNames[indPla] + "</B> :";
        strRes += "<BR><B>Position</B> : " + NF.format(deltaMaxNormPos[iPlanet]) + " au = "
               + NF2.format(deltaMaxNormPos[iPlanet] * KM_PER_AU) + " km.";

        strRes += "<BR>Seen from Earth, this represents an error of <B>"
               + NF2.format(calcGeocentricError(deltaMaxNormPos[iPlanet], indPla))
               + "</B> arc sec.";
        strRes += "<BR><B>Velocity</B> : " + NF.format(deltaMaxNormVel[iPlanet]) + " au/d = "
               + NF2.format(deltaMaxNormVel[iPlanet] * KM_PER_AU * 1000 / 84600) + " m/s.";

        strRes += "<BR><BR>" + LS;

      } // end for each planet

      fos.write(strRes.getBytes());
    }
    catch(Exception e){
      System.out.println("Exception" + e.toString());
      e.printStackTrace();
      System.exit(0);
    }

  }// end buildOutputDiffByPlanet()
*/
  //******************************************************************************
  //******************************** buildOutputSummary() ************************
  //******************************************************************************
/*
  private static void buildOutputSummary(FileOutputStream fos){
    int iPlanet;
    String strRes = new String("");
    try{
      strRes += LS + LS + "<!-- ************************************************** -->" + LS;
      strRes += "<H2><A NAME=\"Summary\"></A>Summary</H2>" + LS;
      strRes += "<TABLE BORDER = \"1\" BGCOLOR=\"#FFFFCC\">" + LS;
      strRes += "<TR><TD></TD>" + LS;
      strRes += "<TD COLSPAN = \"2\"><B>Positions</B></TD>" + LS;
      strRes += "<TD><B>Velocities</B></TD></TR>" + LS;
      strRes += "<TR><TD></TD><TD>Norm of<BR>delta max</TD>" + LS;
      strRes += "<TD>Geocentric<BR>angular error</TD>" + LS;
      strRes += "<TD>Norm of<BR>delta max</TD></TR>" + LS;

      int indPla;
      for(iPlanet = 0; iPlanet < _nbPlan; iPlanet ++){
        indPla = bdlT[iPlanet][0].body.getIndex();
        strRes += "<TR><TD>" + planetNames[bdlT[iPlanet][0].body.getIndex()] + "</TD>" + LS;
        // norm of delta max pos
        strRes+= "<TD>" + NF.format(deltaMaxNormPos[iPlanet]) + " au<BR>" + LS;
        strRes+= NF2.format(deltaMaxNormPos[iPlanet] * KM_PER_AU) + " km</TD>" + LS;
        // geocentric error
        strRes += "<TD><B>"
               + NF2.format(calcGeocentricError(deltaMaxNormPos[iPlanet], indPla))
               + " \"</B></TD>" + LS;
        // norm of delta max vel
        strRes+= "<TD>" + NF.format(deltaMaxNormVel[iPlanet]) + " au/s<BR>" + LS;
        strRes+= NF2.format(deltaMaxNormVel[iPlanet] * KM_PER_AU * 1000/84600) + " m/s - " + LS;
        strRes+= NF2.format(deltaMaxNormVel[iPlanet] * KM_PER_AU * 1000/84600 * 3.6)
              + " km/h</TD>" + LS;
        strRes+= "</TR>" + LS;
      }

      strRes += "</TABLE><BR><BR>";
      fos.write(strRes.getBytes());
    }
    catch(Exception e){
      System.out.println("Exception" + e.toString());
      e.printStackTrace();
      System.exit(0);
    }

  }// end buildOutputSummary()
*/
}//end class TestWithSwiss