astrometric.java

GPS Track connects to a GPS and records the path that you travel. Tracks can be uploaded to a web s

// J2ME Compass
// Copyright (C) 2006 Dana Peters
// http://www.qcontinuum.org/compass

package org.qcontinuum.astro;

import henson.midp.Float;

public class Astrometric {
    
    final static private Float c51544_5 = new Float(515445, -1);
    final static private Float c36525 = new Float(36525);
    final static private Float c99_997361 = new Float(99997361, -6);
    final static private Float c0_993133 = new Float(993133,-6);
    final static private Float c0_7859453 = new Float(7859453, -7);
    final static private Float c6893 = new Float(6893);
    final static private Float c72 = new Float(72);
    final static private Float c2 = new Float(2);
    final static private Float c6191_2 = new Float(61912, -1);
    final static private Float c1296000 = new Float(1296000);
    final static private Float c0 = new Float(0);
    final static private Float c0_606433 = new Float(606433, -6);
    final static private Float c1336_855225 = new Float(1336855225L, -6);
    final static private Float c0_374897 = new Float(374897, -6);
    final static private Float c1325_55241 = new Float(132555241, -5);
    final static private Float c0_827361 = new Float(827361, -6);
    final static private Float c1236_853086 = new Float(1236853086L, -6);
    final static private Float c0_259086 = new Float(259086, -6);
    final static private Float c1342_227825 = new Float(1342227825L, -6);
    final static private Float c412 = new Float(412);
    final static private Float c148 = new Float(148);
    final static private Float c4586 = new Float(4586);
    final static private Float c22640 = new Float(22640);
    final static private Float c55 = new Float(55);
    final static private Float c110 = new Float(110);
    final static private Float c668 = new Float(668);
    final static private Float c165 = new Float(165);
    final static private Float c192 = new Float(192);
    final static private Float c212 = new Float(212);
    final static private Float c769 = new Float(769);
    final static private Float c206 = new Float(206);
    final static private Float c2370 = new Float(2370);
    final static private Float c125 = new Float(125);
    final static private Float c541 = new Float(541);
    final static private Float c206264_8062 = new Float(2062648062L, -4);
    final static private Float c25 = new Float(25);
    final static private Float c11 = new Float(11);
    final static private Float c23 = new Float(23);
    final static private Float c21 = new Float(21);
    final static private Float c44 = new Float(44);
    final static private Float c31 = new Float(31);
    final static private Float c526 = new Float(526);
    final static private Float c18520 = new Float(18520);

    public Astrometric() {
    }
    
    public static EclipticPosition sunPosition(Float MJD)
    {
        // T = (MJD - MJD_J2000) / 36525.0;
        Float T = MJD.Sub(c51544_5).Div(c36525);
        // Mean anomaly and ecliptic longitude
        // M  = pi2 * Float.Frac(0.993133 + 99.997361 * T); 
        Float M = Float.PImul2.Mul(Float.Frac(c0_993133.Add(c99_997361.Mul(T))));
        //  lonSun  = pi2 * Float.Frac(0.7859453 + M / pi2 + 
        //            (6893.0 * sin(M) + 72.0 * sin(2.0 * M) + 6191.2 * T) / 1296.0e3);
        Float x = c0_7859453.Add(M.Div(Float.PImul2)).Add(
                c6893.Mul(Float.sin(M)).Add(c72.Mul(Float.sin(c2.Mul(M)))).Add(c6191_2.Mul(T)).Div(c1296000));
        Float lonSun = Float.PImul2.Mul(Float.Frac(x));
        return new EclipticPosition(c0, lonSun);
    }

    public static EclipticPosition moonPosition(Float MJD)
    {
        // T = (MJD - MJD_J2000) / 36525.0;
        Float T = MJD.Sub(c51544_5).Div(c36525);
        // Mean elements of lunar orbit
        // Float L_0 = Float.Frac(0.606433 + 1336.855225*T);       // mean longitude [rev]
        Float L_0 = Float.Frac(c0_606433.Add(c1336_855225.Mul(T)));
        // l  = pi2 * Float.Frac(0.374897 + 1325.552410 * T);  // Moon's mean anomaly 
        Float l  = Float.PImul2.Mul(Float.Frac(c0_374897.Add(c1325_55241.Mul(T))));
        // ls = pi2 * Float.Frac(0.993133 +   99.997361 * T);  // Sun's mean anomaly 
        Float ls = Float.PImul2.Mul(Float.Frac(c0_993133.Add(c99_997361.Mul(T))));
        // D  = pi2 * Float.Frac(0.827361 + 1236.853086 * T);  // Diff. long. Moon-Sun 
        Float D  = Float.PImul2.Mul(Float.Frac(c0_827361.Add(c1236_853086.Mul(T))));
        // F  = pi2 * Float.Frac(0.259086 + 1342.227825 * T);  // Dist. from ascending node 
        Float F  = Float.PImul2.Mul(Float.Frac(c0_259086.Add(c1342_227825.Mul(T))));

        // Perturbations in longitude and latitude
        // dL = +22640*sin(l) - 4586*sin(l-2*D) + 2370*sin(2*D) +  769*sin(2*l) 
        //      -668*sin(ls) - 412*sin(2*F) - 212*sin(2*l-2*D) - 206*sin(l+ls-2*D)
        //      +192*sin(l+2*D) - 165*sin(ls-2*D) - 125*sin(D) - 110*sin(l+ls)
        //      +148*sin(l-ls) - 55*sin(2*F-2*D);
        Float dL = c22640.Mul(Float.sin(l)).Sub(c4586.Mul(Float.sin(l.Sub(c2.Mul(D))))).Add(c2370.Mul(Float.sin(c2.Mul(D)))).Add(c769.Mul(Float.sin(c2.Mul(l))))
            .Sub(c668.Mul(Float.sin(ls))).Sub(c412.Mul(Float.sin(c2.Mul(F)))).Sub(c212.Mul(Float.sin(c2.Mul(l).Sub(c2.Mul(D))))).Sub(c206.Mul(Float.sin(l.Add(ls).Sub(c2.Mul(D)))))
            .Add(c192.Mul(Float.sin(l.Add(c2.Mul(D))))).Sub(c165.Mul(Float.sin(ls.Sub(c2.Mul(D))))).Sub(c125.Mul(Float.sin(D))).Sub(c110.Mul(Float.sin(l.Add(ls))))
            .Add(c148.Mul(Float.sin(l.Sub(ls)))).Sub(c55.Mul(Float.sin(c2.Mul(F).Sub(c2.Mul(D)))));
        // S  = F + (dL+412*sin(2*F)+541*sin(ls)) / Arcs; 
        Float S  = F.Add(dL.Add(c412.Mul(Float.sin(c2.Mul(F)))).Add(c541.Mul(Float.sin(ls))).Div(c206264_8062)); 
        // h  = F - 2 * D;
        Float h = F.Sub(c2.Mul(D));
        // N  = -526*sin(h) + 44*sin(l+h) - 31*sin(-l+h) - 23*sin(ls+h) 
        //      + 11*sin(-ls+h) - 25*sin(-2*l+F) + 21*sin(-l+F);
        Float N  = c526.Neg().Mul(Float.sin(h)).Add(c44.Mul(Float.sin(l.Add(h)))).Sub(c31.Mul(Float.sin(l.Neg().Add(h)))).Sub(c23.Mul(Float.sin(ls.Add(h))))
            .Add(c11.Mul(Float.sin(ls.Neg().Add(h)))).Sub(c25.Mul(Float.sin(c2.Neg().Mul(l).Add(F)))).Add(c21.Mul(Float.sin(l.Neg().Add(F))));

        // Ecliptic longitude and latitude
        // lonMoon = pi2 * Float.Frac(L_0 + dL / 1296.0e3); // [rad]
        Float lonMoon = Float.PImul2.Mul(Float.Frac(L_0.Add(dL.Div(c1296000))));
        // latMoon = (18520.0 * sin(S) + N) / Arcs;   // [rad]
        Float latMoon = c18520.Mul(Float.sin(S)).Add(N).Div(c206264_8062);
        return new EclipticPosition(latMoon, lonMoon);
    }

}