//*********************************************************************************
// 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] += " (" + 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