risesettime.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 RiseSetTime {
    
    final static public int SUNRISESET = 0;
    final static public int MOONRISESET = 1;
    final static public int CIVILIANTWILIGHT = 2;
    final static public int NAUTICALTWILIGHT = 3;
    final static public int ASTRONOMICALTWILIGHT = 4;

    final static private Float sinh0[] = {
        Float.sin(Float.toRadians(new Float(-50).Div(60))), // Sunrise at h=-50'
        Float.sin(Float.toRadians(new Float(8).Div(60))), // Moonrise at h= +8' 
        Float.sin(Float.toRadians(new Float(-6))), // Civil twilight at h=-6 deg
        Float.sin(Float.toRadians(new Float(-12))), // Nautical twilight at h=-12deg
        Float.sin(Float.toRadians(new Float(-18))) // Astronomical twilight at h=-18deg
      };

    private Float mRiseTime;
    private Float mSetTime;
    private boolean mRises;
    private boolean mSets;
    private boolean mAbove;

    public RiseSetTime(int Event, Float MJD, Float zone, EarthPosition earth) {
        Float phi = Float.toRadians(earth.getLatitude());
        Float lambda = Float.toRadians(earth.getLongitude());
        Float Cphi = Float.cos(phi);
        Float Sphi = Float.sin(phi);

        int hour = 1;
        Float root1 = new Float(0);
        Float root2 = new Float(0);
        Float xe, ye;
        int nRoot;
        // MJD0h = floor(MJD + zone / 24.0) - zone / 24.0;
        Float MJD0h = Float.floor(MJD.Add(zone.Div(24))).Sub(zone.Div(24));

        // Initialize for search
        Float y_minus = SinAlt(Event, MJD0h, hour - 1, lambda, Cphi, Sphi).Sub(sinh0[Event]);

        mAbove = (y_minus.Great(new Float(0)));
        mRises = false;
        mSets  = false;

        // loop over search intervals from [0h-2h] to [22h-24h]
        do {
            Float y_0 = SinAlt(Event, MJD0h, hour, lambda, Cphi, Sphi).Sub(sinh0[Event]);
            Float y_plus = SinAlt(Event, MJD0h, hour + 1, lambda, Cphi, Sphi).Sub(sinh0[Event]);

            // find parabola through three values y_minus,y_0,y_plus
            nRoot = 0;

            // Coefficients of interpolating parabola y=a*x^2+b*x+c
            // a  = 0.5 * (y_plus + y_minus) - y_0;
            Float a  = y_plus.Add(y_minus).Div(2).Sub(y_0);
            // b  = 0.5 * (y_plus - y_minus);
            Float b  = y_plus.Sub(y_minus).Div(2);
            Float c  = y_0;

            // Find extreme value
            // xe = -b / (2.0 * a);
            xe = b.Div(a.Mul(-2)); 
            // ye = (a * xe + b) * xe + c;
            ye = a.Mul(xe).Add(b).Mul(xe).Add(c);
            // dis = b * b - 4.0 * a * c; // Discriminant of y=a*x^2+b*x+c
            Float dis = b.Mul(b).Sub(a.Mul(c).Mul(4));

            if (!dis.Less(new Float(0))) // Parabola has roots 
            {
                // dx = 0.5 * sqrt (dis) / fabs (a);
                Float dx = Float.sqrt(dis).Div(Float.abs(a)).Div(2);
                // root1 = xe - dx; 
                root1 = xe.Sub(dx);
                //root2 = xe + dx;
                root2 = xe.Add(dx);

                if (!Float.abs(root1).Great(new Float(1)))
                    ++nRoot;  
                if (!Float.abs(root2).Great(new Float(1)))
                    ++nRoot;
                if (root1.Less(new Float(-1)))
                    root1 = root2;
            }

            if (nRoot == 1)
            {
                if (y_minus.Less(new Float(0)))
                {
                    mRiseTime = root1.Add(new Float(hour));
                    mRises = true;
                } else {
                    mSetTime  = root1.Add(new Float(hour));
                    mSets  = true;
                }
            }

            if (nRoot == 2)
            {
                if (ye.Less(new Float(0)))
                {
                    mRiseTime = root2.Add(new Float(hour));
                    mSetTime = root1.Add(new Float(hour));
                } else {
                    mRiseTime = root1.Add(new Float(hour));
                    mSetTime = root2.Add(new Float(hour));
                }
                mRises = true;
                mSets = true;
            }          

            y_minus = y_plus; // prepare for next interval
            hour += 2;         
        } while (!((hour == 25) || (mRises && mSets)));
    }

    public String getRiseString()
    {
        return mRises ? timeString(mRiseTime) : "none";
    }
    
    public String getSetString()
    {
        return mSets ? timeString(mSetTime) : "none";
    }

    public String getTransitString()
    {
        return (mRises && mSets) ? timeString(mRiseTime.Add(mSetTime).Div(2)) : "none";
    }

    // SinAlt: Sine of the altitude of Sun or Moon
    // Input:
    //   Event     Indicates event to find
    //   MJD0      0h at date to investigate (as Modified Julian Date)
    //   Hour      Hour
    //   lambda    Geographic east longitude in [rad]
    //   Cphi      Cosine of geographic latitude
    //   Sphi      Sine of geographic latitude
    // <return>:   Sine of the altitude of Sun or Moon at instant of Event
    private static Float SinAlt(int Event, Float MJD0, int Hour, Float lambda, Float Cphi, Float Sphi)
    {
        // MJD = MJD0 + Hour / 24.0;
        Float MJD = MJD0.Add(new Float(Hour).Div(24));

        EclipticPosition ecliptic;
        if (Event == MOONRISESET)
            ecliptic = Astrometric.moonPosition(MJD);
        else
            ecliptic = Astrometric.sunPosition(MJD);
        EquitorialPosition equitorial = ecliptic.toEquitorialPosition();

        // tau = GMST(MJD) + lambda - RA;
        Float tau = UtcDate.GMST(MJD).Add(lambda).Sub(equitorial.getRightAscension());

        // return (Sphi * sin(Dec) + Cphi * cos(Dec) * cos(tau));
        return Sphi.Mul(Float.sin(equitorial.getDeclination()))
            .Add(Cphi.Mul(Float.cos(equitorial.getDeclination())).Mul(Float.cos(tau)));
    }

    private static String timeString(Float time)
    {
        Float x = Float.abs(Float.floor(time.Mul(60).Add(new Float(5, -1))).Div(60).Add(new Float(1, -5)));
        int H = (int)x.toLong();
        x = x.Sub(new Float(H)).Mul(60);
        int M = (int)x.toLong();
        return (H < 10 ? "0" + H : "" + H) + ":" + (M < 10 ? "0" + M : "" + M);
    }

}