Package jsky.util

Source Code of jsky.util.SkyCalc

/**
* $Id: SkyCalc.java,v 1.1.1.1 2009/02/17 22:49:36 abrighto Exp $
*/

package jsky.util;

import jsky.coords.WorldCoords;
import jsky.coords.SiteDesc;
import jsky.plot.ElevationPlotUtil;

import java.util.Date;
import java.util.Calendar;
import java.util.TimeZone;


/**
* Uses a collection of utility routines ported from skycalc.v5.c
* (original C code by John Thorstensen, Dartmouth College) to calculate
* the altitude, azimuth, airmass, and parallactic angle for a given object,
* date, and site.
*/
public class SkyCalc {

    // defined quantities for apparent place transforms ..
    private static final int XFORM_FROMSTD = 1;
    private static final int XFORM_TOSTDEP = -1;
    private static final int XFORM_JUSTPRE = 1;
    private static final int XFORM_DOAPPAR = 0;

    // some (not all) physical, mathematical, and astronomical constants used are defined here.
    private static final double TWOPI = 6.28318530717959;
    private static final double PI_OVER_2 = 1.57079632679490; // From Abramowitz & Stegun
    private static final double ARCSEC_IN_RADIAN = 206264.8062471;
    private static final double DEG_IN_RADIAN = 57.2957795130823;
    private static final double HRS_IN_RADIAN = 3.819718634205;
    private static final double KMS_AUDAY = 1731.45683633// km per sec in 1 AU/day
    private static final double SPEED_OF_LIGHT = 299792.458// in km per sec ... exact.
    private static final double J2000 = 2451545.; // Julian date at standard epoch
    private static final double SEC_IN_DAY = 86400.;
    private static final double FLATTEN = 0.003352813; // flattening of earth, 1/298.257
    private static final double EQUAT_RAD = 6378137.// equatorial radius of earth, meters
    private static final double ASTRO_UNIT = 1.4959787066e11; // 1 AU in meters

    // used in numerical differentiation to find earth velocity -- this value gives
    // about 8 digits of numerical accuracy on the VAX, but is
    // about 3 orders of magnitude larger than the value where roundoff
    // errors become apparent.
    private static final double EARTH_DIFF = 0.05;

    // Constants needed in etcorr method
    static final double[] DELTS = new double[]{
        -2.72,
        3.86,
        10.46,
        17.20,
        21.16,
        23.62,
        24.02,
        23.93,
        24.33,
        26.77,
        29.15,
        31.07,
        33.15,
        35.73,
        40.18,
        45.48,
        50.54,
        54.34,
        56.86,
        60.78,
        62.97,
    };


    /**
     * Used to hold reference params ported from C code
     */
    private static class DoubleRef {
        private double d;

        private DoubleRef() {
            d = 0;
        }
    };

    private static class DateTime {
        short y;
        short mo;
        short d;
        short h;
        short mn;
        double s;

        private DateTime(Date date) {
            Calendar cal = Calendar.getInstance(UT);
            cal.setTime(date);
            y = (short)cal.get(Calendar.YEAR);
            mo = (short)(cal.get(Calendar.MONTH) + 1);
            d = (short)cal.get(Calendar.DAY_OF_MONTH);
            h = (short)cal.get(Calendar.HOUR_OF_DAY);
            mn = (short)cal.get(Calendar.MINUTE);
            s = cal.get(Calendar.SECOND) + cal.get(Calendar.MILLISECOND)/1000.;
        }
    };

    /** UT TimeZone */
    public static final TimeZone UT = TimeZone.getTimeZone("UT");

    // Site parameters
    private double _hoursLongitude;
    private double _degreesLatitude;

//    // Saved date parameter to last call to calculate()
//    private Date _date;

    // calculated results
    private double _altitude;
    private double _azimuth;
    private double _parallacticAngle;
    private double _airmass;

    /**
     * Initialize for the given telescope site.
     * @param site describes the location of the telescope
     */
    public SkyCalc(SiteDesc site) {
        _hoursLongitude = -site.getLongitude()/15.;
        _degreesLatitude = site.getLatitude();
    }

    /**
     * Given the WCS position of an object, the time/date, and a description of the telescope
     * location, calculate the altitude of the object above the horizon, the azimuth,
     * parallactic angle and airmass for the given time. The results are available via
     * properties of this class.
     *
     * @param obj a world coordinates position of an object
     * @param date the date for the calculation
     */
    public void calculate(WorldCoords obj, Date date) {
//        _date = date;
        DateTime dateTime = new DateTime(date);
        DoubleRef jdut = new DoubleRef();
        DoubleRef sid = new DoubleRef();
        DoubleRef curepoch = new DoubleRef();

        setup_time_place(dateTime, _hoursLongitude, jdut, sid, curepoch);

        double objra = obj.getRaDeg()/15;
        double objdec = obj.getDecDeg();
        double objepoch = 2000.;
        getCircumstances(objra, objdec, objepoch, curepoch.d, sid.d, _degreesLatitude);
    }


    /**
     * Return the LST time for the given UT time at the given site.
     */
    public Date getLst(Date date) {
        DateTime dateTime = new DateTime(date);
        double jd = date_to_jd(dateTime);
        double lstHours = lst(jd, _hoursLongitude);
        return _getLst(lstHours, date);
    }

    /**
     * Return the LST time as a Date object for the given LST hours and date.
     */
    private Date _getLst(double lstHours, Date date) {
        Calendar cal = Calendar.getInstance(UT);
        cal.setTime(date);
        int h = cal.get(Calendar.HOUR_OF_DAY);
        boolean nextDay = (lstHours < h);
        CalendarUtil.setHours(cal, lstHours, nextDay);
        return cal.getTime();
    }

    /**
     * Return the airmass for the given altitude in degrees.
     */
    public static double getAirmass(double alt) {
        double secz = secant_z(alt);
        if (secz >= 0.) {
            if (secz < 12.) {
               return true_airmass(secz);
            } else if (secz <= 99.) {
                return secz;
            }
        }
        return 0.;
    }


    /**
     * This takes the date (which contains the time), and the site parameters,
     * and prints out a banner giving the various dates and times; also
     * computes and returns various jd's, the sidereal time, and the epoch.
     * Returns negative number to signal error if date is out of range of
     * validity of algorithms, or if you specify a bad time during daylight-time
     * change; returns zero if successful.
     */
    private static short setup_time_place(DateTime date, double longit,
                                   DoubleRef jdut,
                                   DoubleRef sid,
                                   DoubleRef curepoch) {
        double jd = date_to_jd(date);
        sid.d = lst(jd, longit);
        jdut.d = jd;
        curepoch.d = 2000. + (jd - J2000) / 365.25;
        return (0);
    }


    /*
     * Given object, site, and time information, prints the circumstances
     * of an observation (based on print_circumstances() in skycalc.v5.c).
     */
    private void getCircumstances(double objra, double objdec, double objepoch,
                                  double curep, double sid, double lat) {

        double ha, alt;
        DoubleRef az = new DoubleRef();
        DoubleRef par = new DoubleRef();
        DoubleRef curra = new DoubleRef();
        DoubleRef curdec = new DoubleRef();

        cooxform(objra, objdec, objepoch, curep, curra, curdec, XFORM_JUSTPRE,
                XFORM_FROMSTD);

        ha = adj_time(sid - curra.d);
        alt = altit(curdec.d, ha, lat, az, par);

        _airmass = getAirmass(alt);
//        _lstHours = sid;
//        _lstDate = _getLst(_lstHours, _date);
        _altitude = alt;
        _azimuth = az.d;
        _parallacticAngle = par.d;
    }

    /*
     * General routine for precession and apparent place. Either
     * transforms from current epoch (given by jd) to a standard
     * epoch or back again, depending on value of the switch
     * "from_std"; 1 transforms from standard to current, -1 goes
     * the other way.  Optionally does apparent place including
     * nutation and annual aberration
     * (but neglecting diurnal aberration, parallax, proper motion,
     * and GR deflection of light); switch for this is "just_precess",
     * 1 does only precession, 0 includes other aberration & nutation.
     * <p>
     * Precession uses a matrix procedures
     * as outlined in Taff's Computational Spherical Astronomy book.
     * This is the so-called 'rigorous' method which should give very
     * accurate answers all over the sky over an interval of several
     * centuries.  Naked eye accuracy holds to ancient times, too.
     * Precession constants used are the new IAU1976 -- the 'J2000'
     * system.
     * <p>
     * Nutation is incorporated into matrix formalism by constructing an
     * approximate nutation matrix and taking a matrix product with
     * precession matrix.
     * <p>
     * Aberration is done by adding the vector velocity of the earth to
     * the velocity of the light ray .... not kosher relativistically,
     * but empirically correct to a high order for the angle.
     */
    private static void cooxform(double rin, double din, double std_epoch,
                          double date_epoch, DoubleRef rout, DoubleRef dout,
                          int just_precess, int from_std) {

        /* all the 3-d stuff is declared as [4] 'cause I'm not using the
          zeroth element. */

        double ti, tf, zeta, z, theta;  /* all as per  Taff */
        double cosz, coszeta, costheta, sinz, sinzeta, sintheta;  /* ftns */
        double[][] p = new double[4][4];
        /* elements of the rotation matrix */
        double[][] n = new double[4][4];
        /* elements of the nutation matrix */
        double[][] r = new double[4][4];
        /* their product */
        double[][] t = new double[4][4]/* temporary matrix for inversion .... */
        double radian_ra, radian_dec;

        /* nutation angles in radians */
        DoubleRef del_psi = new DoubleRef();
        DoubleRef del_eps = new DoubleRef();
        double eps;

        double[] orig = new double[4];   /* original unit vector */
        double[] fin = new double[4];   /* final unit vector */
        int i, j, k;


        ti = (std_epoch - 2000.) / 100.;
        tf = (date_epoch - 2000. - 100. * ti) / 100.;

        zeta = (2306.2181 + 1.39656 * ti + 0.000139 * ti * ti) * tf +
                (0.30188 - 0.000344 * ti) * tf * tf + 0.017998 * tf * tf * tf;
        z = zeta + (0.79280 + 0.000410 * ti) * tf * tf + 0.000205 * tf * tf * tf;
        theta = (2004.3109 - 0.8533 * ti - 0.000217 * ti * ti) * tf
                - (0.42665 + 0.000217 * ti) * tf * tf - 0.041833 * tf * tf * tf;

        /* convert to radians */

        zeta = zeta / ARCSEC_IN_RADIAN;
        z = z / ARCSEC_IN_RADIAN;
        theta = theta / ARCSEC_IN_RADIAN;

        /* compute the necessary trig functions for speed and simplicity */

        cosz = Math.cos(z);
        coszeta = Math.cos(zeta);
        costheta = Math.cos(theta);
        sinz = Math.sin(z);
        sinzeta = Math.sin(zeta);
        sintheta = Math.sin(theta);

        /* compute the elements of the precession matrix -- set up
           here as *from* standard epoch *to* input jd. */

        p[1][1] = coszeta * cosz * costheta - sinzeta * sinz;
        p[1][2] = -1. * sinzeta * cosz * costheta - coszeta * sinz;
        p[1][3] = -1. * cosz * sintheta;

        p[2][1] = coszeta * sinz * costheta + sinzeta * cosz;
        p[2][2] = -1. * sinzeta * sinz * costheta + coszeta * cosz;
        p[2][3] = -1. * sinz * sintheta;

        p[3][1] = coszeta * sintheta;
        p[3][2] = -1. * sinzeta * sintheta;
        p[3][3] = costheta;

        if (just_precess == XFORM_DOAPPAR) {  /* if apparent place called for */

            /* do the same for the nutation matrix. */

            nutation_params(date_epoch, del_psi, del_eps);
            eps = 0.409105/* rough obliquity of ecliptic in radians */

            n[1][1] = 1.;
            n[2][2] = 1.;
            n[3][3] = 1.;
            n[1][2] = -1. * del_psi.d * Math.cos(eps);
            n[1][3] = -1. * del_psi.d * Math.sin(eps);
            n[2][1] = -1. * n[1][2];
            n[2][3] = -1. * del_eps.d;
            n[3][1] = -1. * n[1][3];
            n[3][2] = -1. * n[2][3];

            /* form product of precession and nutation matrices ... */
            for (i = 1; i <= 3; i++) {
                for (j = 1; j <= 3; j++) {
                    r[i][j] = 0.;
                    for (k = 1; k <= 3; k++) {
                        r[i][j] += p[i][k] * n[k][j];
                    }
                }
            }
        } else /* if you're just precessing .... */
            for (i = 1; i <= 3; i++) {
                for (j = 1; j <= 3; j++) {
                    r[i][j] = p[i][j]/* simply copy precession matrix */
                }
            }
        }

        /* The inverse of a rotation matrix is its transpose ... */

        if (from_std == XFORM_TOSTDEP) {    /* if you're transforming back to std
                         epoch, rather than forward from std */
            for (i = 1; i <= 3; i++) {
                for (j = 1; j <= 3; j++) {
                    t[i][j] = r[j][i]/* store transpose ... */
                }
            }
            for (i = 1; i <= 3; i++) {
                for (j = 1; j <= 3; j++) {
                    r[i][j] = t[i][j]/* replace original w/ transpose.*/
                }
            }
        }

        /* finally, transform original coordinates */

        radian_ra = rin / HRS_IN_RADIAN;
        radian_dec = din / DEG_IN_RADIAN;

        orig[1] = Math.cos(radian_dec) * Math.cos(radian_ra);
        orig[2] = Math.cos(radian_dec) * Math.sin(radian_ra);
        orig[3] = Math.sin(radian_dec);


        if (from_std == XFORM_TOSTDEP && just_precess == XFORM_DOAPPAR)
        /* if you're transforming from jd to std epoch, and doing apparent place,
       first step is to de-aberrate while still in epoch of date ... */ {
            aberrate(date_epoch, orig, from_std);
        }


        for (i = 1; i <= 3; i++) {
            fin[i] = 0.;
            for (j = 1; j <= 3; j++) {
                fin[i] += r[i][j] * orig[j];
            }
        }

        if (from_std == XFORM_FROMSTD && just_precess == XFORM_DOAPPAR)
        /* if you're transforming from std epoch to jd,
             last step is to apply aberration correction once you're in
             equinox of that jd. */ {
            aberrate(date_epoch, fin, from_std);
        }

        /* convert back to spherical polar coords */

        xyz_cel(fin[1], fin[2], fin[3], rout, dout);

        return;
    }


    /**
     * computes the nutation parameters delta psi and
     * delta epsilon at julian epoch (in years) using approximate
     * formulae given by Jean Meeus, Astronomical Formulae for
     * Calculators, Willman-Bell, 1985, pp. 69-70. Accuracy
     * appears to be a few hundredths of an arcsec or better
     * and numerics have been checked against his example.
     * Nutation parameters are returned in radians.
     */
    private static void nutation_params(double date_epoch, DoubleRef del_psi,
                                 DoubleRef del_ep) {

        double T, jd, L, Lprime, M, Mprime, Omega;

        jd = (date_epoch - 2000.) * 365.25 + J2000;
        T = (jd - 2415020.0) / 36525.;

        L = 279.6967 + (36000.7689 + 0.000303 * T) * T;
        Lprime = 270.4342 + (481267.8831 - 0.001133 * T) * T;
        M = 358.4758 + (35999.0498 - 0.000150 * T) * T;
        Mprime = 296.1046 + (477198.8491 + 0.009192 * T) * T;
        Omega = 259.1833 - (1934.1420 - 0.002078 * T) * T;

        L = L / DEG_IN_RADIAN;
        Lprime = Lprime / DEG_IN_RADIAN;
        M = M / DEG_IN_RADIAN;
        Mprime = Mprime / DEG_IN_RADIAN;
        Omega = Omega / DEG_IN_RADIAN;


        del_psi.d = -1. * (17.2327 + 0.01737 * T) * Math.sin(Omega)
                - (1.2729 + 0.00013 * T) * Math.sin(2. * L)
                + 0.2088 * Math.sin(2 * Omega)
                - 0.2037 * Math.sin(2 * Lprime)
                + (0.1261 - 0.00031 * T) * Math.sin(M)
                + 0.0675 * Math.sin(Mprime)
                - (0.0497 - 0.00012 * T) * Math.sin(2 * L + M)
                - 0.0342 * Math.sin(2 * Lprime - Omega)
                - 0.0261 * Math.sin(2 * Lprime + Mprime)
                + 0.0214 * Math.sin(2 * L - M)
                - 0.0149 * Math.sin(2 * L - 2 * Lprime + Mprime)
                + 0.0124 * Math.sin(2 * L - Omega)
                + 0.0114 * Math.sin(2 * Lprime - Mprime);

        del_ep.d = (9.2100 + 0.00091 * T) * Math.cos(Omega)
                + (0.5522 - 0.00029 * T) * Math.cos(2 * L)
                - 0.0904 * Math.cos(2 * Omega)
                + 0.0884 * Math.cos(2. * Lprime)
                + 0.0216 * Math.cos(2 * L + M)
                + 0.0183 * Math.cos(2 * Lprime - Omega)
                + 0.0113 * Math.cos(2 * Lprime + Mprime)
                - 0.0093 * Math.cos(2 * L - M)
                - 0.0066 * Math.cos(2 * L - Omega);

        del_psi.d = del_psi.d / ARCSEC_IN_RADIAN;
        del_ep.d = del_ep.d / ARCSEC_IN_RADIAN;
    }


    /**
     * A much cleaner rewrite of the original skycalc code for this,
     * which was transcribed from a PL/I routine ....
     */
    private static void xyz_cel(double x, double y, double z, /* cartesian coordinate triplet */
                         DoubleRef ra, DoubleRef dec) { /* corresponding right ascension and declination,
                                          returned in decimal hours and decimal degrees. */
        double mod;    /* modulus */
        double xy;     /* component in xy plane */

        /* normalize explicitly and check for bad input */

        mod = Math.sqrt(x * x + y * y + z * z);
        if (mod > 0.) {
            x = x / mod;
            y = y / mod;
            z = z / mod;
        } else {   /* this has never happened ... */
            System.out.println(
                    "Bad data in xyz_cel .... zero modulus position vector.\n");
            ra.d = 0.;
            dec.d = 0.;
            return;
        }

        xy = Math.sqrt(x * x + y * y);

        if (xy < 1.0e-11) {   /* practically on a pole -- limit is arbitrary ...  */
            ra.d = 0.;   /* degenerate anyway */
            dec.d = PI_OVER_2;
            if (z < 0.) {
                dec.d *= -1.;
            }
        } else { /* in a normal part of the sky ... */
            dec.d = Math.asin(z);
            ra.d = atan_circ(x, y);
        }

        ra.d *= HRS_IN_RADIAN;
        dec.d *= DEG_IN_RADIAN;
    }

    /**
     * corrects celestial unit vector for aberration due to earth's motion.
     * Uses accurate sun position ... replace with crude one for more speed if
     * needed.
     */
    private static void aberrate(double epoch, /* decimal year ...  */
                          double[] vec, /* celestial unit vector ...  */
                          int from_std) {  /* 1 = apply aberration, -1 = take aberration out. */

        double jd, jd1, jd2, Xdot, Ydot, Zdot;   /* page C24 */

        /* throwaways */
        DoubleRef ras = new DoubleRef();
        DoubleRef decs = new DoubleRef();
        DoubleRef dists = new DoubleRef();
        DoubleRef topora = new DoubleRef();
        DoubleRef topodec = new DoubleRef();

        DoubleRef x = new DoubleRef();
        DoubleRef y = new DoubleRef();
        DoubleRef z = new DoubleRef();
        DoubleRef x1 = new DoubleRef();
        DoubleRef y1 = new DoubleRef();
        DoubleRef z1 = new DoubleRef();
        DoubleRef x2 = new DoubleRef();
        DoubleRef y2 = new DoubleRef();
        DoubleRef z2 = new DoubleRef();

        double norm;

        /* find heliocentric velocity of earth as a fraction of the speed of light ... */

        jd = J2000 + (epoch - 2000.) * 365.25;
        jd1 = jd - EARTH_DIFF;
        jd2 = jd + EARTH_DIFF;

        accusun(jd1, 0., 0., ras, decs, dists, topora, topodec, x1, y1, z1);
        accusun(jd2, 0., 0., ras, decs, dists, topora, topodec, x2, y2, z2);
        accusun(jd, 0., 0., ras, decs, dists, topora, topodec, x, y, z);

        Xdot = KMS_AUDAY * (x2.d - x1.d) / (2. * EARTH_DIFF * SPEED_OF_LIGHT)/* numerical differentiation */
        Ydot = KMS_AUDAY * (y2.d - y1.d) / (2. * EARTH_DIFF * SPEED_OF_LIGHT)/* crude but accurate */
        Zdot = KMS_AUDAY * (z2.d - z1.d) / (2. * EARTH_DIFF * SPEED_OF_LIGHT);

        /* approximate correction ... non-relativistic but very close.  */

        vec[1] += from_std * Xdot;
        vec[2] += from_std * Ydot;
        vec[3] += from_std * Zdot;

        norm = Math.pow((vec[1] * vec[1] + vec[2] * vec[2] + vec[3] * vec[3]), 0.5);

        vec[1] = vec[1] / norm;
        vec[2] = vec[2] / norm;
        vec[3] = vec[3] / norm;
    }


    /**
     * returns radian angle 0 to 2pi for coords x, y --
     * get that quadrant right !!
     */
    private static double atan_circ(double x, double y) {
        double theta;

        if ((x == 0.) && (y == 0.)) {
            return (0.)/* guard ... */
        }

        theta = Math.atan2(y, x);
        while (theta < 0.) {
            theta += TWOPI;
        }
        return (theta);
    }


    /**
     * implemenataion of Jean Meeus' more accurate solar
     * ephemeris.  For ultimate use in helio correction! From
     * Astronomical Formulae for Calculators, pp. 79 ff.  This
     * gives sun's position wrt *mean* equinox of date, not
     * apparent*.  Accuracy is << 1 arcmin.  Positions given are
     * geocentric ... parallax due to observer's position on earth is
     * ignored. This is up to 8 arcsec; routine is usually a little
     * better than that.
     * // -- topocentric correction *is* included now. -- //
     * Light travel time is apparently taken into
     * account for the ra and dec, but I don't know if aberration is
     * and I don't know if distance is simlarly antedated.
     * <p/>
     * x, y, and z are heliocentric equatorial coordinates of the
     * EARTH, referred to mean equator and equinox of date.
     */
    private static void accusun(double jd, double lst, double geolat,
                         DoubleRef ra, DoubleRef dec, DoubleRef dist, DoubleRef topora, DoubleRef topodec,
                         DoubleRef x, DoubleRef y, DoubleRef z) {

        double L, T, Tsq, Tcb;
        double M, e, Cent, nu, sunlong;
        double Mrad, nurad, R;
        double A, B, C, D, E, H;
        double xtop, ytop, ztop, topodist, l, m, n;
        DoubleRef xgeo = new DoubleRef();
        DoubleRef ygeo = new DoubleRef();
        DoubleRef zgeo = new DoubleRef();

        jd = jd + etcorr(jd) / SEC_IN_DAY;  /* might as well do it right .... */
        T = (jd - 2415020.) / 36525./* 1900 --- this is an oldish theory*/
        Tsq = T * T;
        Tcb = T * Tsq;
        L = 279.69668 + 36000.76892 * T + 0.0003025 * Tsq;
        M = 358.47583 + 35999.04975 * T - 0.000150 * Tsq - 0.0000033 * Tcb;
        e = 0.01675104 - 0.0000418 * T - 0.000000126 * Tsq;

        L = circulo(L);
        M = circulo(M);
/*      printf("raw L, M: %15.8f, %15.8f\n",L,M); */

        A = 153.23 + 22518.7541 * T;  /* A, B due to Venus */
        B = 216.57 + 45037.5082 * T;
        C = 312.69 + 32964.3577 * T;  /* C due to Jupiter */
        /* D -- rough correction from earth-moon
            barycenter to center of earth. */
        D = 350.74 + 445267.1142 * T - 0.00144 * Tsq;
        E = 231.19 + 20.20 * T;    /* "inequality of long period .. */
        H = 353.40 + 65928.7155 * T;  /* Jupiter. */

        A = circulo(A) / DEG_IN_RADIAN;
        B = circulo(B) / DEG_IN_RADIAN;
        C = circulo(C) / DEG_IN_RADIAN;
        D = circulo(D) / DEG_IN_RADIAN;
        E = circulo(E) / DEG_IN_RADIAN;
        H = circulo(H) / DEG_IN_RADIAN;

        L = L + 0.00134 * Math.cos(A)
                + 0.00154 * Math.cos(B)
                + 0.00200 * Math.cos(C)
                + 0.00179 * Math.sin(D)
                + 0.00178 * Math.sin(E);

        Mrad = M / DEG_IN_RADIAN;

        Cent = (1.919460 - 0.004789 * T - 0.000014 * Tsq) * Math.sin(Mrad)
                + (0.020094 - 0.000100 * T) * Math.sin(2.0 * Mrad)
                + 0.000293 * Math.sin(3.0 * Mrad);
        sunlong = L + Cent;


        nu = M + Cent;
        nurad = nu / DEG_IN_RADIAN;

        R = (1.0000002 * (1 - e * e)) / (1. + e * Math.cos(nurad));
        R = R + 0.00000543 * Math.sin(A)
                + 0.00001575 * Math.sin(B)
                + 0.00001627 * Math.sin(C)
                + 0.00003076 * Math.cos(D)
                + 0.00000927 * Math.sin(H);

        sunlong = sunlong / DEG_IN_RADIAN;

        dist.d = R;
        x.d = Math.cos(sunlong)/* geocentric */
        y.d = Math.sin(sunlong);
        z.d = 0.;
        eclrot(jd, y, z);

        /*      --- code to include topocentric correction for sun .... */

        geocent(lst, geolat, 0., xgeo, ygeo, zgeo);

        xtop = x.d - xgeo.d * EQUAT_RAD / ASTRO_UNIT;
        ytop = y.d - ygeo.d * EQUAT_RAD / ASTRO_UNIT;
        ztop = z.d - zgeo.d * EQUAT_RAD / ASTRO_UNIT;

        topodist = Math.sqrt(xtop * xtop + ytop * ytop + ztop * ztop);

        l = xtop / (topodist);
        m = ytop / (topodist);
        n = ztop / (topodist);

        topora.d = atan_circ(l, m) * HRS_IN_RADIAN;
        topodec.d = Math.asin(n) * DEG_IN_RADIAN;

        ra.d = atan_circ(x.d, y.d) * HRS_IN_RADIAN;
        dec.d = Math.asin(z.d) * DEG_IN_RADIAN;

        x.d = x.d * R * -1/* heliocentric */
        y.d = y.d * R * -1;
        z.d = z.d * R * -1;
    }

    /**
     * Given a julian date in 1900-2100, returns the correction
     * delta t which is:
     * TDT - UT (after 1983 and before 1998)
     * ET - UT (before 1983)
     * an extrapolated guess  (after 1998).
     * <p/>
     * For dates in the past (<= 1998 and after 1900) the value is linearly
     * interpolated on 5-year intervals; for dates after the present,
     * an extrapolation is used, because the true value of delta t
     * cannot be predicted precisely.  Note that TDT is essentially the
     * modern version of ephemeris time with a slightly cleaner
     * definition.
     * <p/>
     * Where the algorithm shifts there is an approximately 0.1 second
     * discontinuity.  Also, the 5-year linear interpolation scheme can
     * lead to errors as large as 0.5 seconds in some cases, though
     * usually rather smaller.
     */
    private static double etcorr(double jd) {

        double[] dates = new double[22];
        double year, delt = 0.;
        int i;

        for (i = 0; i <= 19; i++) {
            dates[i] = 1900 + (double) i * 5.;
        }
        dates[20] = 1998./* the last accurately tabulated one in the
                                2000 Almanac ... */

        year = 1900. + (jd - 2415019.5) / 365.25;

        if (year < 1998. && year >= 1900.) {
            i = (int) (year - 1900) / 5;
            delt = DELTS[i] +
                    ((DELTS[i + 1] - DELTS[i]) / (dates[i + 1] - dates[i])) *
                    (year - dates[i]);
        } else if (year >= 1998. && year < 2100.) {
            delt = 33.15 + (2.164e-3) * (jd - 2436935.4)/* rough extrapolation */
        } else if (year < 1900) {
            System.out.println("etcorr ... no ephemeris time data for < 1900.\n");
            delt = 0.;
        } else if (year >= 2100.) {
            System.out.println(
                    "etcorr .. very long extrapolation in delta T - inaccurate.\n");
            delt = 180.; /* who knows? */
        }

        return (delt);
    }

    /**
     * assuming x is an angle in degrees, returns
     * modulo 360 degrees.
     */
    private static double circulo(double x) {
        int n = (int) (x / 360.);
        return (x - 360. * n);
    }

    /**
     * rotates ecliptic rectangular coords x, y, z to
     * equatorial (all assumed of date.)
     */
    private static void eclrot(double jd, DoubleRef y, DoubleRef z) {
        double incl;
//        double xpr;
        double ypr;
        double zpr;
        double T;

        T = (jd - J2000) / 36525/* centuries since J2000 */

        incl = (23.439291 + T * (-0.0130042 - 0.00000016 * T)) / DEG_IN_RADIAN;
        /* 1992 Astron Almanac, p. B18, dropping the
           cubic term, which is 2 milli-arcsec! */
        ypr = Math.cos(incl) * y.d - Math.sin(incl) * z.d;
        zpr = Math.sin(incl) * y.d + Math.cos(incl) * z.d;
        y.d = ypr;
        z.d = zpr;
        /* x remains the same. */
    }


    /**
     * computes the geocentric coordinates from the geodetic
     * (standard map-type) longitude, latitude, and height.
     * These are assumed to be in decimal hours, decimal degrees, and
     * meters respectively.  Notation generally follows 1992 Astr Almanac,
     * p. K11
     */
    private static void geocent(double geolong, double geolat, double height,
                         DoubleRef x_geo, DoubleRef y_geo, DoubleRef z_geo) {

        double denom, C_geo, S_geo;

        geolat = geolat / DEG_IN_RADIAN;
        geolong = geolong / HRS_IN_RADIAN;
        denom = (1. - FLATTEN) * Math.sin(geolat);
        denom = Math.cos(geolat) * Math.cos(geolat) + denom * denom;
        C_geo = 1. / Math.sqrt(denom);
        S_geo = (1. - FLATTEN) * (1. - FLATTEN) * C_geo;
        C_geo = C_geo + height / EQUAT_RAD;  /* deviation from almanac
                   notation -- include height here. */
        S_geo = S_geo + height / EQUAT_RAD;
        x_geo.d = C_geo * Math.cos(geolat) * Math.cos(geolong);
        y_geo.d = C_geo * Math.cos(geolat) * Math.sin(geolong);
        z_geo.d = S_geo * Math.sin(geolat);
    }


    /**
     * adjusts a time (decimal hours) to be between -12 and 12,
     * generally used for hour angles.
     */
    private static double adj_time(double x) {
        if (Math.abs(x) < 100000.) {  /* too inefficient for this! */
            while (x > 12.) {
                x = x - 24.;
            }
            while (x < -12.) {
                x = x + 24.;
            }
        } else {
            System.out.println("Out of bounds in adj_time!\n");
        }
        return (x);
    }

    /**
     * returns altitude(degr) for dec, ha, lat (decimal degr, hr, degr);
     * also computes and returns azimuth through pointer argument,
     * and as an extra added bonus returns parallactic angle (decimal degr)
     * through another pointer argument.
     *
     * @param dec target declination in degrees
     * @param ha  the hour angle in hours
     * @param lat the observer's latitude in radians
     * @return the parallactic angle in degrees
     */
    private static double altit(double dec, double ha, double lat,
                         DoubleRef az, DoubleRef parang) {

        double x, y, z;
        double sinp, cosp;  /* sin and cos of parallactic angle */
        double cosdec, sindec, cosha, sinha, coslat, sinlat;
        /* time-savers ... */

        dec = dec / DEG_IN_RADIAN;
        ha = ha / HRS_IN_RADIAN;
        lat = lat / DEG_IN_RADIAN;  /* thank heavens for pass-by-value */
        cosdec = Math.cos(dec);
        sindec = Math.sin(dec);
        cosha = Math.cos(ha);
        sinha = Math.sin(ha);
        coslat = Math.cos(lat);
        sinlat = Math.sin(lat);
        x = DEG_IN_RADIAN * Math.asin(cosdec * cosha * coslat + sindec * sinlat);
        y = sindec * coslat - cosdec * cosha * sinlat; /* due N comp. */
        z = -1. * cosdec * sinha; /* due east comp. */
        az.d = Math.atan2(z, y);

        /* as it turns out, having knowledge of the altitude and
               azimuth makes the spherical trig of the parallactic angle
               less ambiguous ... so do it here!  Method uses the
           "astronomical triangle" connecting celestial pole, object,
               and zenith ... now know all the other sides and angles,
               so we can crush it ... */

        if (cosdec != 0.) { /* protect divide by zero ... */
            sinp = -1. * Math.sin(az.d) * coslat / cosdec;
            /* spherical law of sines .. note cosdec = sin of codec,
                coslat = sin of colat .... */
            cosp = -1. * Math.cos(az.d) * cosha - Math.sin(az.d) * sinha * sinlat;
            /* spherical law of cosines ... also transformed to local
                          available variables. */
            parang.d = Math.atan2(sinp, cosp) * DEG_IN_RADIAN;
            /* let the library function find the quadrant ... */
        } else { /* you're on the pole */
            if (lat >= 0.) {
                parang.d = 180.;
            } else {
                parang.d = 0.;
            }
        }

        az.d *= DEG_IN_RADIAN;  /* done with taking trig functions of it ... */
        while (az.d < 0.) {
            az.d += 360./* force 0 -> 360 */
        }
        while (az.d >= 360.) {
            az.d -= 360.;
        }

        return (x);
    }

    /**
     * Computes the secant of z, assuming the object is not
     * too low to the horizon; returns 100. if the object is
     * low but above the horizon, -100. if the object is just
     * below the horizon.
     */
    private static double secant_z(double alt) {

        double secz;
        if (alt != 0) {
            secz = 1. / Math.sin(alt / DEG_IN_RADIAN);
        } else {
            secz = 100.;
        }
        if (secz > 100.) {
            secz = 100.;
        }
        if (secz < -100.) {
            secz = -100.;
        }
        return (secz);
    }


    /**
     * returns the true airmass for a given secant z.
     * The expression used is based on a tabulation of the mean KPNO
     * atmosphere given by C. M. Snell & A. M. Heiser, 1968,
     * PASP, 80, 336.  They tabulated the airmass at 5 degr
     * intervals from z = 60 to 85 degrees; I fit the data with
     * a fourth order poly for (secz - airmass) as a function of
     * (secz - 1) using the IRAF curfit routine, then adjusted the
     * zeroth order term to force (secz - airmass) to zero at
     * z = 0.  The poly fit is very close to the tabulated points
     * (largest difference is 3.2e-4) and appears smooth.
     * This 85-degree point is at secz = 11.47, so for secz > 12
     * I just return secz - 1.5 ... about the largest offset
     * properly determined.
     */
    private static double true_airmass(double secz) {

        double seczmin1;
        int i, ord = 4;
        double[] coef = new double[5];
        double result = 0;

        coef[1] = 2.879465E-3/* sun compilers do not allow automatic
                initializations of arrays. */
        coef[2] = 3.033104E-3;
        coef[3] = 1.351167E-3;
        coef[4] = -4.716679E-5;
        if (secz < 0.) {
            return (-1.)/* out of range. */
        }
        if (secz > 12) {
            return (secz - 1.5)/* shouldn't happen .... */
        }
        seczmin1 = secz - 1.;
        /* evaluate polynomial ... */
        for (i = ord; i > 0; i--) {
            result = (result + coef[i]) * seczmin1;
        }
        /* no zeroth order term. */
        result = secz - result;
        return (result);

    }

    /*
    From Meeus' Astronomical Formulae for Calculators.  The two JD
    conversion routines routines were replaced 1998 November 29 to
    avoid inclusion of copyrighted "Numerical Recipes" code.  A test
    of 1 million random JDs between 1585 and 3200 AD gave the same
    conversions as the NR routines.
    */
    private static double date_to_jd(DateTime date) {
        double jd;
        int y, m;
        long A, B;

        if (date.mo <= 2) {
            y = date.y - 1;
            m = date.mo + 12;
        } else {
            y = date.y;
            m = date.mo;
        }

        A = (long) (y / 100.);
        B = 2 - A + (long) (A / 4.);

        jd = (long) (365.25 * y) + (long) (30.6001 * (m + 1)) + date.d +
                1720994.5;

        jd += date.h / 24. + date.mn / 1440. + date.s / 86400.;

        if (date.y > 1583) {
            return (jd + B);
        } else {
            return (jd);
        }
        /* Not quite right, since Gregorian calendar first
        adopted around Oct 1582.  But fine for modern. */
    }


    /**
     * returns the local MEAN sidereal time (dec hrs) at julian date jd
     * at west longitude long (decimal hours).  Follows
     * definitions in 1992 Astronomical Almanac, pp. B7 and L2.
     * Expression for GMST at 0h ut referenced to Aoki et al, A&A 105,
     * p.359, 1982.  On workstations, accuracy (numerical only!)
     * is about a millisecond in the 1990s.
     */
    private static double lst(double jd, double longit) {

        double t, ut, jdmid, jdint, jdfrac, sid_g;
        long jdin, sid_int;

        jdin = (long) jd;         /* fossil code from earlier package which
                                   split jd into integer and fractional parts ... */
        jdint = jdin;
        jdfrac = jd - jdint;
        if (jdfrac < 0.5) {
            jdmid = jdint - 0.5;
            ut = jdfrac + 0.5;
        } else {
            jdmid = jdint + 0.5;
            ut = jdfrac - 0.5;
        }
        t = (jdmid - J2000) / 36525;
        sid_g = (24110.54841 + 8640184.812866 * t + 0.093104 * t * t -
                6.2e-6 * t * t * t) / SEC_IN_DAY;
        sid_int = (long) sid_g;
        sid_g = sid_g - (double) sid_int;
        sid_g = sid_g + 1.0027379093 * ut - longit / 24.;
        sid_int = (long) sid_g;
        sid_g = (sid_g - (double) sid_int) * 24.;
//        if (sid_g < 0.) {
//            sid_g = sid_g + 24.;
//        }
        return (sid_g);
    }

    /**
     * Returns the altitude in degrees calculated in the last call
     * to {@link #calculate(WorldCoords, Date)}
     */
    public double getAltitude() {
        return _altitude;
    }

    /**
     * Returns the azimuth calculated in the last call
     * to {@link #calculate(WorldCoords, Date)}
     */
    public double getAzimuth() {
        return _azimuth;
    }

    /**
     * Returns the parallactic angle in degrees calculated in the last call
     * to {@link #calculate(WorldCoords, Date)}
     */
    public double getParallacticAngle() {
        return _parallacticAngle;
    }

    /**
     * Returns the airmass calculated in the last call
     * to {@link #calculate(WorldCoords, Date)}
     */
    public double getAirmass() {
        return _airmass;
    }


    /** Test main. */
    public static void main(String[] args) {
        SiteDesc site = ElevationPlotUtil.MAUNA_KEA;
        Calendar cal = Calendar.getInstance(UT);

        cal.set(2005, Calendar.MARCH, 9, 11, 0, 0);
        Date date = cal.getTime();
        WorldCoords obj = new WorldCoords("10 00 00", "+19 49 26");
        SkyCalc skyCalc = new SkyCalc(site);
        skyCalc.calculate(obj, date);
        System.out.println("For Site = Mauna Kea, date = 2005-03-09, obj = " + obj + ":");
        System.out.println("   Altitude = " + skyCalc.getAltitude());
        System.out.println("   Azimuth = " + skyCalc.getAzimuth());
        System.out.println("   Airmass = " + skyCalc.getAirmass());
        System.out.println("   Parallactic Angle = " + skyCalc.getParallacticAngle());
    }
}
TOP

Related Classes of jsky.util.SkyCalc

TOP
Copyright © 2018 www.massapi.com. All rights reserved.
All source code are property of their respective owners. Java is a trademark of Sun Microsystems, Inc and owned by ORACLE Inc. Contact coftware#gmail.com.