PlanetaryEphemeris.java

// Calculate the positions of the
// major planets using the current "osculating elements"
// from the Astronomical Almanac.
//
// A simple elliptical orbit
// is assumed for both the planet and the Earth - no
// corrections are made from within the program as the
// osculating elements will already take account of
// perturbations.
//
// The method used here is based on finding the rectangular
// coordinates of the planet and of the Earth, and then applying
// successive coordinate transformations to find the rectangular
// gocentric equatorial coordinates of the planet.
//
// From "Practical Astronomy with your calculator or spreadsheet
// 4th edition"; Duffet-Smith & Zwart.

import java.util.*;

class PlanetaryEphemeris {

  RaDec getPlanet(int planet) {

    // Orbital periods of the planets, in tropical years.
    double[] op = new double[8];
    op[0] = 0.24085;   // Mercury.
    op[1] = 0.615207;  // Venus.
    op[2] = 0.999996;  // Earth.
    op[3] = 1.880765;  // Mars.
    op[4] = 11.857911; // Jupiter.
    op[5] = 29.310579; // Saturn.
    op[6] = 84.039492; // Uranus.
    op[7] = 165.84539; // Neptune.

    // Longitude of the planet's J2000 epoch, in degrees.
    double[] le = new double[8];
    le[0] = 75.5671;    // Mercury.
    le[1] = 272.30044;  // Venus.
    le[2] = 99.556772;  // Earth.
    le[3] = 109.09646;  // Mars.
    le[4] = 337.917132; // Jupiter.
    le[5] = 172.398316; // Saturn.
    // le[6] = 271.063148; // Uranus.
    // le[6] = 313.238105; // Uranus.
    le[6] = 356.0;       // Uranus.
    le[7] = 326.895127; // Neptune.

    // Longitude of the planet's perihelion. in degrees.
    double[] lp = new double[8];
    lp[0] = 77.5671;     // Mercury.
    lp[1] = 131.54;      // Venus.
    lp[2] = 103.2055;    // Earth.
    lp[3] = 336.217;     // Mars.
    lp[4] = 14.6633;     // Jupiter.
    lp[5] = 89.567;      // Saturn.
    lp[6] = 172.8844833; // Uranus.
    lp[7] = 23.07;       // Neptune.

    // Eccentricity of the planet's orbit.
    double[] ep = new double[8];
    ep[0] = 0.205627; // Mercury.
    ep[1] = 0.006812; // Venus.
    ep[2] = 0.016671; // Earth.
    ep[3] = 0.093384; // Mars.
    ep[4] = 0.048907; // Jupiter.
    ep[5] = 0.053853; // Saturn.
    ep[6] = 0.046321; // Uranus.
    ep[7] = 0.010483; // Neptune.

    // Semi major axis of the planet's orbit in AU.
    double[] sm = new double[8];
    sm[0] = 0.387098; // Mercury.
    sm[1] = 0.723329; // Venus.
    sm[2] = 0.999985; // Earth.
    sm[3] = 1.523689; // Mars.
    sm[4] = 5.20278;  // Jupiter.
    sm[5] = 9.51134;  // Saturn.
    sm[6] = 19.21814; // Uranus.
    sm[7] = 30.1985;  // Neptune.

    // Inclination of the planet's orbit. in degrees.
    double[] ip = new double[8];
    ip[0] = 7.0051;   // Mercury.
    ip[1] = 3.3947;   // Venus.
    ip[2] = 0.0;      // Earth.
    ip[3] = 1.8497;   // Mars.
    ip[4] = 1.3035;   // Jupiter.
    ip[5] = 2.4873;   // Saturn.
    ip[6] = 0.773059; // Uranus.
    ip[7] = 1.7673;   // Neptune.

    // Longitude of the planet's ascending node, in degrees.
    double[] la = new double[8];
    la[0] = 48.449;    // Mercury.
    la[1] = 76.769;    // Venus.
    la[2] = 0.0;       // Earth.
    la[3] = 49.632;    // Mars.
    la[4] = 100.595;   // Jupiter.
    la[5] = 113.752;   // Saturn.
    la[6] = 73.926961; // Uranus.
    la[7] = 131.879;   // Neptune.

    JavaIO q = new JavaIO();        // Debug.
    AstroUtil au = new AstroUtil();
    AU_II a2 = new AU_II();

    double r2d = 180 / Math.PI,
           d2r = Math.PI / 180;

    // Get the number of days from 2010.
    GregorianCalendar clndr = new GregorianCalendar();
    long unixTime = clndr.getTimeInMillis() / 1000L;
    double d = (double) (unixTime - 1262304000L) / 86400.0;

    // Calculate the mean anomaly of the planet.
    double ma = ((360 / 365.242191) * (d / op[planet]));
    ma = fixRange(ma);

    double maa = ma + le[planet] - lp[planet];

    double ta = maa + (360 / Math.PI) * ep[planet] * Math.sin(maa * d2r);		
    ta = fixRange(ta);
    
    // Heliocentric longitude.
    double hl = ta + lp[planet]; 
    hl = fixRange(hl);

    // The distance of the planet from the sun.
    double rd = sm[planet] * (1 - (ep[planet] * ep[planet]));
    double rn = 1 + (ep[planet] * Math.cos(ta * d2r));
    double r = rd / rn;

    // q.p("The distance of the planet from the sun (r): " + r);

    // Find the position of the earth. Begin with the mean anomaly.
    double meE = (360 / 365.242191) * (d / op[2]);
    meE = fixRange(meE);

    // q.p("meE: " + meE);

    double maE = meE + le[2] - lp[2];
    maE = fixRange(maE);

    // q.p("Mean anomaly of earth (maE/Me): " + maE);

    // The Earth's true anomaly.
    // double taE = a2.solveKepler(ep[2], (maE * d2r));
    // taE = fixRange(taE * r2d);
    // q.p("Iterated True Earth anomaly (Vp): " + t2);  //TEST

    double taE = maE + ((360 / Math.PI) * ep[2] * Math.sin(maE * d2r));
    taE = fixRange(taE);
    // q.p("True anomaly of earth (taE/Ve): " + taE);

    // Earth's heliocentric longitude.
    double hlE = taE + lp[2]; 
    hlE = fixRange(hlE);

    // q.p("Earth's heliocentric longitude (hlE/L): " + hlE);

    // The distance of the Earth from the sun.
    double rEn = sm[2] * (1 - (ep[2] * ep[2]));
    double rEd = 1 + (ep[2] * Math.cos(taE * d2r));
    double rE = rEn / rEd;

    // q.p("The distance of the Earth from the sun (rE/R): " + rE);

    // Planet's heliocentric latitude.
    double hlp = Math.asin((Math.sin((hl - la[planet]) * d2r)) *
                Math.sin(ip[planet] * d2r));

    // q.p("Planet's heliocentric latitude (hlp/phi): " + (hlp * r2d));

    // Planet's coordinates along the ecliptic.
    double y = Math.sin((hl - la[planet]) * d2r) * Math.cos(ip[planet] * d2r);
    double x = Math.cos((hl - la[planet]) * d2r);

    // q.p("X: " + x);
    // q.p("Y: " + y);

    // Find the arc tangent of y/x.
    double angle = Math.atan(y / x);
    if (x < 0) { angle += Math.PI; }
    else if (y < 0) { angle += Math.PI * 2; }

    // q.p("Planet's coordinates wrt ecliptic (tan-1(x/y)): " + (angle * r2d));

    // The projected heliocentric longitude.
    double phl = (angle * r2d) + la[planet];

    // q.p("The projected heliocentric longitude (phl/l'): " + phl);
 
    // The projected radius vector.
    double prv = r * Math.cos(hlp);   

    // q.p("The projected radius vector (prv/r'): " + prv);

    // The geocentric ecliptic longitude.
    double geln = rE * Math.sin((phl - hlE) * d2r);
    double geld = prv - (rE * Math.cos((phl - hlE) * d2r)); 
    double gel = Math.atan(geln / geld) + (phl * d2r);
    gel = fixRange(gel * r2d);

    // q.p("The geocentric ecliptic latitude (gel/lamda): " + gel);

    // The geocentric ecliptic latitude.
    double gelan = prv * Math.tan(hlp) * Math.sin((gel - phl) * d2r); 
    double gelad = rE * Math.sin((phl - hlE) * d2r); 
    double gela = Math.atan(gelan / gelad);
    gela = fixRange(gela * r2d);

    // q.p("The geocentric ecliptic latitude (gela/beta): " + gela);

    // From here on out, all angles will be in radians.
    gel *= d2r;
    gela *= d2r;

    // Convert gel and gela to ra and dec. Calculate the obliquity of the
    // ecliptic first.
    double te = (d - 3655) / 36525; // Centuries since 2000 Jan 1.
    double de = (46.815 * te) + (0.0006 * te * te) + (0.00181 * te * te * te);
    double oe = (23.439292 - (de / 3600)) * d2r;

    // q.p("The obliquity of the ecliptic (oe/epsilon): " + (oe * r2d));

    // gel = 139.686111 * d2r;
    // gela = 4.875278 * d2r;

    // Right ascension:
    y = (Math.sin(gel) * Math.cos(oe)) - (Math.tan(gela) * Math.sin(oe));
    x = Math.cos(gel);

    // q.p("y: " + y + ", x: " + x);

    double ra = Math.atan(y / x);
    if (x < 0) { ra += Math.PI; }
    else if (y < 0) { ra += Math.PI * 2; }

    // Declination:
    double dec1 = Math.sin(gela) * Math.cos(oe);
    double dec2 = Math.cos(gela) * Math.sin(oe) * Math.sin(gel);

    // q.p("sin dec: " + (dec1 + dec2));

    double dec = Math.asin(dec1 + dec2);

    RaDec coord = new RaDec();
    coord.ra =  ra;
    coord.dec = dec;
    coord.sra =  au.r2ra(ra);
    coord.sdec = au.r2dec(dec);
    coord.name = getName(planet);

    return coord;
  }

  // Ensure that an angle is between 0 and 360 degrees.
  double fixRange(double ma) {
    if (ma > 360) {
      while (ma > 360) { ma -= 360; }
    } else if (ma < 0) {
      while (ma < 0) { ma += 360; }
    }
    return ma;
  }

  String getName(int planet) {
    String name = null;
    if      (planet == 0) { name = "Mercury"; }
    else if (planet == 1) { name = "Venus"; }
    else if (planet == 3) { name = "Mars"; }
    else if (planet == 4) { name = "Jupiter"; }
    else if (planet == 5) { name = "Saturn"; }
    else if (planet == 6) { name = "Uranus"; }
    else if (planet == 7) { name = "Neptune"; }
    return name;
  }

  // Unit test.
  public static void main(String[] s) {
    PlanetaryEphemeris pe = new PlanetaryEphemeris();
    for (int i = 0; i < 8; i++) {
      if (i != 2) { pe.getPlanet(i); }
    }
    // pe.getPlanet(4);
  }
}