/*
 * ESO Archive
 *
 * $Id: WorldCoords.java,v 1.2 2009/02/20 23:10:11 abrighto Exp $
 *
 * who             when        what
 * --------------  ----------  ----------------------------------------
 * Allan Brighton  1999/05/03  Created, based on C++ version
 */

package jsky.coords;

import java.awt.geom.Point2D;
import java.io.Serializable;

/**
 * Class representing world coordinates (right-ascension, declination,
 * stored as J2000 internally).
 */
public class WorldCoords implements WorldCoordinates, Serializable {

    /** The RA coordinate */
    protected HMS ra;

    /** The DEC coordinate */
    protected DMS dec;

    /** Default constructor: initialize null coordinates. */
    public WorldCoords() {
        ra = new HMS();
        dec = new DMS();
    }

    /**
     * Initialize from ra (hours) and dec (degrees).
     *
     * @param ra The RA value in hours
     * @param dec The DEC value in degrees
     * @param equinox The equinox of the input coordinates.
     */
    public WorldCoords(HMS ra, DMS dec, double equinox) {
        this.ra = ra;
        this.dec = dec;
        init(equinox);
    }

    /**
     * Initialize from ra (hours) and dec (degrees) in J2000.
     *
     * @param ra The RA value in hours
     * @param dec The DEC value in degrees
     */
    public WorldCoords(HMS ra, DMS dec) {
        this(ra, dec, 2000.);
    }

    /**
     * Initialize from ra and dec (in degrees).
     *
     * @param ra The RA value in degrees
     * @param dec The DEC value in degrees
     * @param equinox The equinox of the input coordinates.
     */
    public WorldCoords(double ra, double dec, double equinox) {
        this.ra = new HMS(ra / 15);
        this.dec = new DMS(dec);
        init(equinox);
    }

    /**
     * Initialize from ra and dec (in degrees).
     *
     * @param ra The RA value in degrees
     * @param dec The DEC value in degrees
     * @param equinox The equinox of the input coordinates.
     */
    public WorldCoords(Double ra, Double dec, double equinox) {
        this.ra = new HMS(ra / 15);
        this.dec = new DMS(dec);
        init(equinox);
    }

    /**
     * Initialize from ra and dec (in degrees J2000).
     *
     * @param ra The RA value in degrees
     * @param dec The DEC value in degrees
     */
    public WorldCoords(double ra, double dec) {
        this(ra, dec, 2000.);
    }

    /**
     * Initialize from an ra, dec point in the given equinox.
     *
     * @param p The RA,Dec point in degrees
     */
    public WorldCoords(Point2D.Double p, double equinox) {
        this(p.getX(), p.getY(), equinox);
    }

    /**
     * Initialize from an ra, dec point (in degrees J2000).
     *
     * @param p The RA,Dec point in degrees
     */
    public WorldCoords(Point2D.Double p) {
        this(p.getX(), p.getY(), 2000.);
    }

    /**
     * Initialize from ra and dec (in degrees J2000).
     *
     * @param ra The RA value in degrees
     * @param dec The DEC value in degrees
     */
    public WorldCoords(Double ra, Double dec) {
        this(ra.doubleValue(), dec.doubleValue(), 2000.);
    }

    /**
     * Initialize from RA hours, minutes, seconds and and DEC degress, minutes, seconds.
     *
     * @param rh RA hours
     * @param rm RA minutes
     * @param rs RA seconds
     * @param dd DEC degrees
     * @param dm DEC minutes
     * @param ds DEC seconds
     * @param equinox The equinox of the input coordinates.
     */
    public WorldCoords(double rh, int rm, double rs, double dd, int dm, double ds, double equinox) {
        ra = new HMS(rh, rm, rs);
        dec = new DMS(dd, dm, ds);
        init(equinox);
    }

    /**
     * Initialize from RA hours, minutes, seconds and and DEC degress, minutes, seconds
     * in J2000.
     *
     * @param rh RA hours
     * @param rm RA minutes
     * @param rs RA seconds
     * @param dd DEC degrees
     * @param dm DEC minutes
     * @param ds DEC seconds
     */
    public WorldCoords(double rh, int rm, double rs, double dd, int dm, double ds) {
        this(rh, rm, rs, dd, dm, ds, 2000.);
    }

    /**
     * Parse RA and DEC in string format.
     *
     * @param ra The RA value in the form "hh mm ss.s", "hh:mm:ss.s",
     *           or just "hh.hhh".
     *
     * @param dec The DEC value in the form "[+/-]dd mm ss.s", "[+/-]dd:mm:ss.s"
     *            or just "[+/-]dd.ddd".
     *
     * @param equinox The equinox of the input coordinates.
     */
    public WorldCoords(String ra, String dec, double equinox) {
        this.ra = new HMS(ra);
        this.dec = new DMS(dec);
        init(equinox);
    }

    /**
     * Parse RA and DEC in string format.
     *
     * @param ra The RA value in the form "hh mm ss.s", "hh:mm:ss.s",
     *           or just "hh.hhh".
     *
     * @param dec The DEC value in the form "[+/-]dd mm ss.s", "[+/-]dd:mm:ss.s"
     *            or just "[+/-]dd.ddd".
     *
     * @param equinox The equinox of the input coordinates.
     * @param hflag if true, always assume ra is in hours, otherwise, assume deg if
     *        ra is a decimal value.
     */
    public WorldCoords(String ra, String dec, double equinox, boolean hflag) {
        this.ra = new HMS(ra, hflag);
        this.dec = new DMS(dec);
        init(equinox);
    }

    /**
     * Parse RA and DEC in string format (assume J2000).
     *
     * @param ra The RA value in the form "hh mm ss.s", "hh:mm:ss.s",
     *           or just "hh.hhh".
     *
     * @param dec The DEC value in the form "[+/-]dd mm ss.s", "[+/-]dd:mm:ss.s"
     *            or just "[+/-]dd.ddd".
     */
    public WorldCoords(String ra, String dec) {
        this(ra, dec, 2000.);
    }

    /**
     * Parse RA and DEC in string format (assume J2000).
     *
     * @param ra The RA value in the form "hh mm ss.s", "hh:mm:ss.s",
     *           or just "hh.hhh".
     *
     * @param dec The DEC value in the form "[+/-]dd mm ss.s", "[+/-]dd:mm:ss.s"
     *            or just "[+/-]dd.ddd".
     *
     * @param hflag if true, always assume ra is in hours, otherwise, assume deg if
     *        ra is a decimal value.
     */
    public WorldCoords(String ra, String dec, boolean hflag) {
        this(ra, dec, 2000., hflag);
    }

    /** return the RA value */
    public HMS getRA() {
        return ra;
    }

    /** return the DEC value */
    public DMS getDec() {
        return dec;
    }

    /** return the RA value in deg */
    public double getRaDeg() {
        return ra.getVal() * 15.;
    }

    /** return the DEC value in deg */
    public double getDecDeg() {
        return dec.getVal();
    }

    /** return the X (ra) coordinate in deg */
    public double getX() {
        return ra.getVal() * 15.;
    }

    /** return the Y *dec) coordinate in deg */
    public double getY() {
        return dec.getVal();
    }

    /**
     * Return the coordinates as a string in h:m:s [+-]d:m:s format
     */
    public String toString() {
        return ra.toString() + ", " + dec.toString() + " J2000";
    }

    /**
     * Format RA and DEC in h:m:s [+-]d:m:s format in the given equinox
     * and return them as a 2 element String array.
     */
    public String[] format(double equinox) {
        String[] ar = new String[2];
        if (equinox == 2000.0) {
            ar[0] = ra.toString();
            ar[1] = dec.toString();
        } else {
            // make tmp copy and convert equinox before printing
            WorldCoords tmp = new WorldCoords(ra, dec);
            tmp.convertEquinox(2000.0, equinox);
            ar[0] = tmp.ra.toString();
            ar[1] = tmp.dec.toString();
        }
        return ar;
    }

    /**
     * Format RA and DEC in h:m:s [+-]d:m:s format in J2000
     * and return them as a 2 element String array.
     */
    public String[] format() {
        return format(2000.);
    }

    /**
     * Return RA and DEC in degrees as an array of 2 doubles {ra, dec} in the given equinox.
     */
    public double[] getRaDec(double equinox) {
        double[] ar = new double[2];
        if (equinox == 2000.0) {
            ar[0] = getRaDeg();
            ar[1] = getDecDeg();
        } else {
            // make tmp copy and convert equinox before printing
            WorldCoords tmp = new WorldCoords(ra, dec);
            tmp.convertEquinox(2000.0, equinox);
            ar[0] = tmp.getRaDeg();
            ar[1] = tmp.getDecDeg();
        }
        return ar;
    }


    /**
     * Computes distance and position angle solving a spherical
     * triangle (no approximations).
     * The arguments are the coordinates in decimal degrees and
     * the result is an array containing the 2 values.
     *
     * @param dra0 center RA
     * @param decd0 center DEC
     * @param dra point RA
     * @param decd point DEC
     *
     * @return An array of 2 doubles ar[2], where ar[0] is the distance in arcmin and
     *         ar[1] is the position angle phi in degrees (East of North).
     *
     * (Based on the C version from A. P. Martinez.)
     */
    public static double[] dispos(double dra0, double decd0, double dra, double decd) {
        double[] ar = new double[2];		// return array
        double alf, alf0, del, del0, phi;
        double sd, sd0, cd, cd0, cosda, cosd, sind, sinpa, cospa;
        double radian = 180. / Math.PI;

        // coo transformed in radiants
        alf = dra / radian;
        alf0 = dra0 / radian;
        del = decd / radian;
        del0 = decd0 / radian;

        sd0 = Math.sin(del0);
        sd = Math.sin(del);
        cd0 = Math.cos(del0);
        cd = Math.cos(del);
        cosda = Math.cos(alf - alf0);
        cosd = sd0 * sd + cd0 * cd * cosda;
        double dist = Math.acos(cosd);
        phi = 0.0;
        if (dist > 0.0000004) {
            sind = Math.sin(dist);
            cospa = (sd * cd0 - cd * sd0 * cosda) / sind;
//            if (cospa > 1.0)
//                cospa = 1.0;
            if (Math.abs(cospa) > 1.0)
                cospa=cospa/Math.abs(cospa); // 2005-06-02: fix from awicenec@eso.org
            sinpa = cd * Math.sin(alf - alf0) / sind;
            phi = Math.acos(cospa) * radian;
            if (sinpa < 0.0)
                phi = 360.0 - phi;
        }
        dist *= radian;
        dist *= 60.0;
        if (decd0 == 90.)
            phi = 180.0;
        if (decd0 == -90.)
            phi = 0.0;

        ar[0] = dist;
        ar[1] = phi;
        return ar;
    }

    /**
     * return the distance between this position and the given one in arcmin
     * and also set the position angle
     *
     * @param pos The other point.
     *
     * @return An array of 2 doubles ar[2], where ar[0] is the distance in arcmin and
     *         ar[1] is the position angle in degrees (East of North).
     */
    public double[] dispos(WorldCoordinates pos) {
        return dispos(getRaDeg(), getDecDeg(), pos.getRaDeg(), pos.getDecDeg());
    }

    /**
     * Compute the distance in degrees between the two given coordinates.
     * (Based on the C version in D. Mink's wcssubs package.)
     */
    public static double wcsdist(double x1, double y1, double x2, double y2) {
        double xr1, xr2, yr1, yr2, w, diff, cosb;
        double[] pos1 = new double[3], pos2 = new double[3];
        int i;

        // Convert two vectors to direction cosines
        double deg = Math.PI / 180.;
        xr1 = x1 * deg;
        yr1 = y1 * deg;
        cosb = Math.cos(yr1);
        pos1[0] = Math.cos(xr1) * cosb;
        pos1[1] = Math.sin(xr1) * cosb;
        pos1[2] = Math.sin(yr1);

        xr2 = x2 * deg;
        yr2 = y2 * deg;
        cosb = Math.cos(yr2);
        pos2[0] = Math.cos(xr2) * cosb;
        pos2[1] = Math.sin(xr2) * cosb;
        pos2[2] = Math.sin(yr2);

        // Modulus squared of half the difference vector
        w = 0.0;
        for (i = 0; i < 3; i++)
            w = w + (pos1[i] - pos2[i]) * (pos1[i] - pos2[i]);

        w = w / 4.0;
        if (w > 1.0)
            w = 1.0;

        // Angle beween the vectors
        diff = 2.0 * Math.atan2(Math.sqrt(w), Math.sqrt(1.0 - w));

        // convert to deg and return
        return diff * 180. / Math.PI;
    }


    /**
     * Return the distance between this position and the given one in arcmin.
     * @param pos The other point.
     * @return The distance to the given point in arcmin.
     */
    public double dist(WorldCoords pos) {
        double[] ar = dispos(getRaDeg(), getDecDeg(), pos.getRaDeg(), pos.getDecDeg());
        return ar[0];
    }

    /**
     * Return the distance between this position and the given one in arcmin.
     * @param pos The other point.
     * @return The distance to the given point in arcmin.
     */
    public double dist(Coordinates pos) {
        return dist((WorldCoords) pos);
    }

    /**
     * static member to get the distance between 2 points in arcmin
     *
     * @param ra0 The first RA value
     * @param dec0 The first DEC value
     * @param ra1 The second RA value
     * @param dec1 The second DEC value
     *
     * @return The distance in arcmin.
     */
    static public double dist(double ra0, double dec0, double ra1, double dec1) {
        double[] ar = dispos(ra0, dec0, ra1, dec1);
        return ar[0];
    }

    /**
     * Given a radius in arcmin, return an array {pos1, pos2} with the 2 endpoints
     * that form a box with center at "this" position.
     *
     * @param radius The radius in arcmin.
     *
     * @return Array of 2 WorldCoords objects that are the endpoints of a box
     *         with the given radius and centered at "this" position.
     */
    public WorldCoordinates[] box(double radius) {
        WorldCoordinates[] ar = new WorldCoordinates[2]; // return array

        // get units in degrees
        double ra = this.ra.getVal(), dec = this.dec.getVal();
        radius /= 60.0;

        // get width of square
        double width = Math.sqrt(2.0 * radius * radius);
        double r1, r2, d1, d2, cosdec;

        d1 = dec - width / 2.0;
        if (d1 <= -90.0) {
            d1 = -90.0;
            d2 = dec + width / 2.0;
            r1 = 0.0;
            r2 = 24.0;
        } else {
            d2 = dec + width / 2.0;
            if (d2 >= 90.0) {
                d1 = dec - width / 2.0;
                d2 = 90.0;
                r1 = 0.0;
                r2 = 24.0;
            } else {
                if (dec > 0.0)
                    cosdec = Math.abs(Math.cos(d1 * Math.PI / 180.));
                else
                    cosdec = Math.abs(Math.cos(d2 * Math.PI / 180.));

                r1 = ra - width / 15. / 2 / cosdec;
                r2 = ra + width / 15. / 2 / cosdec;

                if (r1 < 0.0)
                    r1 += 24;
                if (r2 > 24.0)
                    r2 -= 24;
            }
        }

        // set the result array
        ar[0] = new WorldCoords(r1 * 15., d1);
        ar[1] = new WorldCoords(r2 * 15., d2);

        return ar;
    }


    /**
     * Given the endpoints of a box (pos1, pos2), return an array containing the
     * center ra, dec in degrees, as well as the width, height, and radius of the
     * box in arcmin.
     *
     * @param pos1 The first endpoint of the box.
     * @param pos2 The second endpoint of the box.
     *
     * @return An array of 5 doubles: {ra, dec, width, height, radius}
     *         where (ra, dec) gives the center position in deg,
     *         width and height are the size of the box in arcmin,
     *         and radius is the distance from the center to a corner in arcmin.
     */
    public double[] center(WorldCoordinates pos1, WorldCoordinates pos2) {
        double ar[] = new double[5]; // result

        // get center pos
        double ra1 = pos1.getRaDeg(), dec1 = pos1.getDecDeg();
        double ra2 = pos2.getRaDeg(), dec2 = pos2.getDecDeg();
        ar[0] = (ra1 + ra2) / 2.0;
        ar[1] = (dec1 + dec2) / 2.0;

        // get width and height of box
        ar[2] = (wcsdist(ra1, dec1, ra2, dec1) * 60.);
        ar[3] = (wcsdist(ra1, dec1, ra1, dec2) * 60.);

        // radius is half the distance from pos1 to pos2
        ar[4] = (wcsdist(ra1, dec1, ra2, dec2) * 60.) / 2.;

        return ar;
    }


    /** Returns the name of the coordinate system as a string. */
    public String getCoordinateSystemName() {
        return "wcs";
    }

    /**
     * Called by constructors to initialize the new object.
     * Set dec to include the sign and check the ra and dec ranges.
     * The equinox argument indicates the equinox of ra and dec.
     * Internally, ra and dec are converted if needed and stored in J2000.
     */
    private void init(double equinox) {
        checkRange();
        convertEquinox(equinox, 2000.);
    }

    /**
     * Convert the coordinates from fromEquinox to toEquinox.
     */
    private void convertEquinox(double fromEquinox, double toEquinox) {
        if (fromEquinox == toEquinox)
            return;

        double[] q0 = new double[2], q1 = new double[2];
        q0[0] = ra.getVal() * 15;	// hours to degrees
        q0[1] = dec.getVal();
        JPrec.prej_q(q0, q1, fromEquinox, toEquinox);
        ra = new HMS(q1[0] / 15);	// degrees to hours
        dec = new DMS(q1[1]);
    }

    /** check range of ra,dec values */
    private void checkRange() {
        double ra = this.ra.getVal(), dec = this.dec.getVal();

        if (ra < -0.001 || ra >= 25.0) {
            // System.out.println("XXX RA value " + ra + " out of range (0..24 hours)");
            // throw new IllegalArgumentException("RA value " + ra + " out of range (0..24 hours)");
        }

        if (dec < -90. || dec > 90.) {
            // System.out.println("XXX DEC value " + dec + " out of range (-90..+90 deg)");
            // throw new IllegalArgumentException("DEC value " + dec + " out of range (-90..+90 deg)");
        }
    }

    /**
     * Test cases
     */
    public static void main(String[] args) {
        WorldCoords c1 = new WorldCoords(49.95096, 41.51173);
        WorldCoords c2 = new WorldCoords(3, 19, 48.2304, 41, 30, 42.228);
        WorldCoords c3 = new WorldCoords(new HMS(3, 19, 48.2304), new DMS(41, 30, 42.228));
        WorldCoords c4 = new WorldCoords(new HMS(c1.getRA()), new DMS(c1.getDec()));
        WorldCoords c5 = new WorldCoords("3 19 48.2304", "+41 30 42.228", 2000.0);
        WorldCoords c6 = new WorldCoords("3:19:48.2304", "+41:30:42.228", 2000.0);
        WorldCoords c7 = new WorldCoords(Double.toString(49.95096 / 15.), "41.51173", 2000.0);

        System.out.println("these coords should all be the same (or very close):");
        System.out.println(c1);
        System.out.println(c2);
        System.out.println(c3);
        System.out.println(c4);
        System.out.println(c5);
        System.out.println(c6);
        System.out.println(c7);

        c1 = new WorldCoords(49.95096, -41.51173);
        c2 = new WorldCoords(3, 19, 48.2304, -41, 30, 42.228);
        c3 = new WorldCoords(new HMS(3, 19, 48.2304), new DMS(-41, 30, 42.228));
        c4 = new WorldCoords(new HMS(c1.getRA()), new DMS(c1.getDec()));
        c5 = new WorldCoords("3 19 48.2304", "-41 30 42.228", 2000.0);
        c6 = new WorldCoords("3:19:48.2304", "-41:30:42.228", 2000.0);
        c7 = new WorldCoords(Double.toString(49.95096 / 15.), "-41.51173", 2000.0);

        System.out.println("Here is the same with negative dec:");
        System.out.println(c1);
        System.out.println(c2);
        System.out.println(c3);
        System.out.println(c4);
        System.out.println(c5);
        System.out.println(c6);
        System.out.println(c7);

        WorldCoords c8 = new WorldCoords("3:19", "+41:30", 2000.0);
        WorldCoords c9 = new WorldCoords("3", "+41", 2000.0);
        System.out.println("And with missing minutes, ... seconds, ...:");
        System.out.println(c8);
        System.out.println(c9);

        // test the "box" method (get 2 points given a radius)
        WorldCoordinates c10 = new WorldCoords("03:19:48.243", "+41:30:40.31"), c11, c12;
        WorldCoordinates[] ar1 = c10.box(7.05);
        c11 = ar1[0];
        c12 = ar1[1];
        System.out.println("box of radius 7.05 with center at (03:19:48.243, +41:30:40.31):");
        System.out.println(c11);
        System.out.println(c12);

        // test values at or near 0,0
        WorldCoords c13 = new WorldCoords("0", "+41:30:40.31");
        System.out.println("With ra = 0.0: " + c13 + " ("
                           + c13.getRA().getVal() + ", " + c13.getDec().getVal() + ")");

        WorldCoords c14 = new WorldCoords("0.0", "-0.0");
        System.out.println("With ra = 0.0, dec = -0.0: " + c14 + " ("
                           + c14.getRA().getVal() + ", " + c14.getDec().getVal() + ")");

        WorldCoords c15 = new WorldCoords("0:0:1", "-0:1:1");
        System.out.println("With ra = 0:0:1, dec = -0:1:1: " + c15 + " ("
                           + c15.getRA().getVal() + ", " + c15.getDec().getVal() + ")");

        // test conversion between h:m:s and deg and back
        WorldCoords c16 = new WorldCoords("22:45:22.74", "-39:34:14.63");
        System.out.println("test conversion between h:m:s and deg and back");
        System.out.println("22:45:22.74 -39:34:14.63 == " + c16 + " == " + c16.toString());

        String[] ar2 = c16.format();
        System.out.println(" == " + ar2[0] + " " + ar2[1]);
        WorldCoords c17 = new WorldCoords(ar2[0], ar2[1]);
        System.out.println(" == " + c17);

        // test equinox conversion
        WorldCoords c18 = new WorldCoords(0.0, 0.0, 1950.);
        String[] ar3 = c18.format(1950.);
        System.out.println("00:00:00 B1950 == " + c18 + " J2000 == " + ar3[0] + " " + ar3[1] + " B1950");
    }
}