//=== File Prolog =============================================================
//	This code was developed by NASA, Goddard Space Flight Center, Code 588
//	for the Scientist's Expert Assistant (SEA) project.
//
//--- Contents ----------------------------------------------------------------
//	class platepos
//
//--- Description -------------------------------------------------------------
//
//--- Notes -------------------------------------------------------------------
//
//--- Development History -----------------------------------------------------
//
//	07/16/98	J. Jones / 588
//
//		Original implementation.
//
//--- DISCLAIMER---------------------------------------------------------------
//
//	This software is provided "as is" without any warranty of any kind, either
//	express, implied, or statutory, including, but not limited to, any
//	warranty that the software will conform to specification, any implied
//	warranties of merchantability, fitness for a particular purpose, and
//	freedom from infringement, and any warranty that the documentation will
//	conform to the program, or any warranty that the software will be error
//	free.
//
//	In no event shall NASA be liable for any damages, including, but not
//	limited to direct, indirect, special or consequential damages, arising out
//	of, resulting from, or in any way connected with this software, whether or
//	not based upon warranty, contract, tort or otherwise, whether or not
//	injury was sustained by persons or property or otherwise, and whether or
//	not loss was sustained from or arose out of the results of, or use of,
//	their software or services provided hereunder.
//
//=== End File Prolog =========================================================

/* File saoimage/wcslib/platepos.c
 * February 25, 1996
 * By Doug Mink, Harvard-Smithsonian Center for Astrophysics

 * Module:	platepos.c (Plate solution WCS conversion
 * Purpose:	Compute WCS from Digital Sky Survey plate fit
 * Subroutine:	platepos() converts from pixel location to RA,Dec
 * Subroutine:	platepix() converts from RA,Dec to pixel location

    These functions are based on the astrmcal.c portion of GETIMAGE by
    J. Doggett and the documentation distributed with the Digital Sky Survey.

*/

package jsky.coords;

import java.awt.geom.*;

public class platepos {

    // Terms used by all methods
    private static final double COND2R = 1.745329252e-2;

    // Terms used by getPosition()
    private static final double CONS2R = 206264.8062470964;
    private static final double TWOPI = 6.28318530717959;

    // Terms used by getPixels()
    private static final int MAX_ITERATIONS = 50;
    private static final double TOLERANCE = 0.0000005;
    private static final double CONR2S = 206264.8062470964;

    /* Routine to determine accurate position for pixel coordinates */
    /* returns 0 if successful otherwise 1 = angle too large for projection; */
    /* based on amdpos() from getimage */
    /* Input: */
    /* x pixel number  (RA or long without rotation) */
    /* y pixel number  (dec or lat without rotation) */
    /* Output: */
    /* Right ascension or longitude in degrees */
    /* Declination or latitude in degrees */
    public static Point2D.Double getPosition(double xpix, double ypix, WCSTransform wcs) {
        double x, y, xmm, ymm, xmm2, ymm2, xmm3, ymm3, x2y2;
        double xi, xir, eta, etar, raoff, ra, dec;
        double ctan, ccos;

        /*  Ignore magnitude and color terms
        double mag = 0.0;
        double color = 0.0; */

        /* Convert from image pixels to plate pixels */
        x = xpix + wcs.x_pixel_offset - 1.0 + 0.5;
        y = ypix + wcs.y_pixel_offset - 1.0 + 0.5;

        /* Convert from pixels to millimeters */
        xmm = (wcs.ppo_coeff[2] - x * wcs.x_pixel_size) / 1000.0;
        ymm = (y * wcs.y_pixel_size - wcs.ppo_coeff[5]) / 1000.0;
        xmm2 = xmm * xmm;
        ymm2 = ymm * ymm;
        xmm3 = xmm * xmm2;
        ymm3 = ymm * ymm2;
        x2y2 = xmm2 + ymm2;

        /*  Compute coordinates from x,y and plate model */

        xi = wcs.amd_x_coeff[0] * xmm + wcs.amd_x_coeff[1] * ymm +
                wcs.amd_x_coeff[2] + wcs.amd_x_coeff[3] * xmm2 +
                wcs.amd_x_coeff[4] * xmm * ymm + wcs.amd_x_coeff[5] * ymm2 +
                wcs.amd_x_coeff[6] * (x2y2) + wcs.amd_x_coeff[7] * xmm3 +
                wcs.amd_x_coeff[8] * xmm2 * ymm + wcs.amd_x_coeff[9] * xmm * ymm2 +
                wcs.amd_x_coeff[10] * ymm3 + wcs.amd_x_coeff[11] * xmm * (x2y2) +
                wcs.amd_x_coeff[12] * xmm * x2y2 * x2y2;

        /*  Ignore magnitude and color terms
        + wcs.amd_x_coeff[13]*mag	+ wcs.amd_x_coeff[14]*mag*mag +
        wcs.amd_x_coeff[15]*mag*mag*mag + wcs.amd_x_coeff[16]*mag*xmm +
        wcs.amd_x_coeff[17]*mag*x2y2	+ wcs.amd_x_coeff[18]*mag*xmm*x2y2 +
        wcs.amd_x_coeff[19]*color; */

        eta = wcs.amd_y_coeff[0] * ymm + wcs.amd_y_coeff[1] * xmm +
                wcs.amd_y_coeff[2] + wcs.amd_y_coeff[3] * ymm2 +
                wcs.amd_y_coeff[4] * xmm * ymm + wcs.amd_y_coeff[5] * xmm2 +
                wcs.amd_y_coeff[6] * (x2y2) + wcs.amd_y_coeff[7] * ymm3 +
                wcs.amd_y_coeff[8] * ymm2 * xmm + wcs.amd_y_coeff[9] * ymm * xmm2 +
                wcs.amd_y_coeff[10] * xmm3 + wcs.amd_y_coeff[11] * ymm * (x2y2) +
                wcs.amd_y_coeff[12] * ymm * x2y2 * x2y2;

        /*  Ignore magnitude and color terms
        + wcs.amd_y_coeff[13]*mag	+ wcs.amd_y_coeff[14]*mag*mag +
        wcs.amd_y_coeff[15]*mag*mag*mag + wcs.amd_y_coeff[16]*mag*ymm +
        wcs.amd_y_coeff[17]*mag*x2y2)	+ wcs.amd_y_coeff[18]*mag*ymm*x2y2 +
        wcs.amd_y_coeff[19]*color; */

        /* Convert to radians */

        xir = xi / CONS2R;
        etar = eta / CONS2R;

        /* Convert to RA and Dec */

        ctan = Math.tan(wcs.plate_dec);
        ccos = Math.cos(wcs.plate_dec);
        raoff = Math.atan2(xir / ccos, 1.0 - etar * ctan);
        ra = raoff + wcs.plate_ra;
        if (ra < 0.0) ra = ra + TWOPI;

        dec = Math.atan(Math.cos(raoff) / ((1.0 - (etar * ctan)) / (etar + ctan)));

        return new Point2D.Double(ra / COND2R, dec / COND2R);
    }

    /* Routine to determine pixel coordinates for sky position */
    /* returns 0 if successful otherwise 1 = angle too large for projection; */
    /* based on amdinv() from getimage */
    /* Input: */
    /* Right ascension or longitude in degrees */
    /* Declination or latitude in degrees */
    /* Output: */
    /* x pixel number  (RA or long without rotation) */
    /* y pixel number  (dec or lat without rotation) */
    public static Point2D.Double getPixels(double xpos, double ypos, WCSTransform wcs) {
        double div,xi,eta,x,y,xy,x2,y2,x2y,y2x,x3,y3,x4,y4,x2y2,cjunk,dx,dy;
        double sypos,cypos,syplate,cyplate,sxdiff,cxdiff;
        double f,fx,fy,g,gx,gy, xmm, ymm;
        int i;
        double xr, yr; 	/* position in radians */

        /* Convert RA and Dec in radians to standard coordinates on a plate */
        xr = xpos * COND2R;
        yr = ypos * COND2R;
        sypos = Math.sin(yr);
        cypos = Math.cos(yr);
        syplate = Math.sin(wcs.plate_dec);
        cyplate = Math.cos(wcs.plate_dec);
        sxdiff = Math.sin(xr - wcs.plate_ra);
        cxdiff = Math.cos(xr - wcs.plate_ra);
        div = (sypos * syplate) + (cypos * cyplate * cxdiff);
        xi = cypos * sxdiff * CONR2S / div;
        eta = ((sypos * cyplate) - (cypos * syplate * cxdiff)) * CONR2S / div;

        /* Set initial value for x,y */
        xmm = xi / wcs.plate_scale;
        ymm = eta / wcs.plate_scale;

        /* Iterate by Newton's method */
        for (i = 0; i < MAX_ITERATIONS; i++) {
            /* X plate model */
            xy = xmm * ymm;
            x2 = xmm * xmm;
            y2 = ymm * ymm;
            x2y = x2 * ymm;
            y2x = y2 * xmm;
            x2y2 = x2 + y2;
            cjunk = x2y2 * x2y2;
            x3 = x2 * xmm;
            y3 = y2 * ymm;
            x4 = x2 * x2;
            y4 = y2 * y2;
            f = wcs.amd_x_coeff[0] * xmm + wcs.amd_x_coeff[1] * ymm +
                    wcs.amd_x_coeff[2] + wcs.amd_x_coeff[3] * x2 +
                    wcs.amd_x_coeff[4] * xy + wcs.amd_x_coeff[5] * y2 +
                    wcs.amd_x_coeff[6] * x2y2 + wcs.amd_x_coeff[7] * x3 +
                    wcs.amd_x_coeff[8] * x2y + wcs.amd_x_coeff[9] * y2x +
                    wcs.amd_x_coeff[10] * y3 + wcs.amd_x_coeff[11] * xmm * x2y2 +
                    wcs.amd_x_coeff[12] * xmm * cjunk;
            /* magnitude and color terms ignored
            + wcs.amd_x_coeff[13]*mag +
            wcs.amd_x_coeff[14]*mag*mag   + wcs.amd_x_coeff[15]*mag*mag*mag +
            wcs.amd_x_coeff[16]*mag*xmm   + wcs.amd_x_coeff[17]*mag*(x2+y2) +
            wcs.amd_x_coeff[18]*mag*xmm*(x2+y2)  + wcs.amd_x_coeff[19]*color;
            */

            /*  Derivative of X model wrt x */
            fx = wcs.amd_x_coeff[0] + wcs.amd_x_coeff[3] * 2.0 * xmm +
                    wcs.amd_x_coeff[4] * ymm + wcs.amd_x_coeff[6] * 2.0 * xmm +
                    wcs.amd_x_coeff[7] * 3.0 * x2 + wcs.amd_x_coeff[8] * 2.0 * xy +
                    wcs.amd_x_coeff[9] * y2 + wcs.amd_x_coeff[11] * (3.0 * x2 + y2) +
                    wcs.amd_x_coeff[12] * (5.0 * x4 + 6.0 * x2 * y2 + y4);
            /* magnitude and color terms ignored
            wcs.amd_x_coeff[16]*mag      + wcs.amd_x_coeff[17]*mag*2.0*xmm +
            wcs.amd_x_coeff[18]*mag*(3.0*x2+y2);
            */

            /* Derivative of X model wrt y */
            fy = wcs.amd_x_coeff[1] + wcs.amd_x_coeff[4] * xmm +
                    wcs.amd_x_coeff[5] * 2.0 * ymm + wcs.amd_x_coeff[6] * 2.0 * ymm +
                    wcs.amd_x_coeff[8] * x2 + wcs.amd_x_coeff[9] * 2.0 * xy +
                    wcs.amd_x_coeff[10] * 3.0 * y2 + wcs.amd_x_coeff[11] * 2.0 * xy +
                    wcs.amd_x_coeff[12] * 4.0 * xy * x2y2;
            /* magnitude and color terms ignored
            wcs.amd_x_coeff[17]*mag*2.0*ymm +
            wcs.amd_x_coeff[18]*mag*2.0*xy;
            */

            /* Y plate model */
            g = wcs.amd_y_coeff[0] * ymm + wcs.amd_y_coeff[1] * xmm +
                    wcs.amd_y_coeff[2] + wcs.amd_y_coeff[3] * y2 +
                    wcs.amd_y_coeff[4] * xy + wcs.amd_y_coeff[5] * x2 +
                    wcs.amd_y_coeff[6] * x2y2 + wcs.amd_y_coeff[7] * y3 +
                    wcs.amd_y_coeff[8] * y2x + wcs.amd_y_coeff[9] * x2y +
                    wcs.amd_y_coeff[10] * x3 + wcs.amd_y_coeff[11] * ymm * x2y2 +
                    wcs.amd_y_coeff[12] * ymm * cjunk;
            /* magnitude and color terms ignored
            wcs.amd_y_coeff[13]*mag        + wcs.amd_y_coeff[14]*mag*mag +
            wcs.amd_y_coeff[15]*mag*mag*mag + wcs.amd_y_coeff[16]*mag*ymm +
            wcs.amd_y_coeff[17]*mag*x2y2 +
            wcs.amd_y_coeff[18]*mag*ymm*x2y2 + wcs.amd_y_coeff[19]*color;
            */

            /* Derivative of Y model wrt x */
            gx = wcs.amd_y_coeff[1] + wcs.amd_y_coeff[4] * ymm +
                    wcs.amd_y_coeff[5] * 2.0 * xmm + wcs.amd_y_coeff[6] * 2.0 * xmm +
                    wcs.amd_y_coeff[8] * y2 + wcs.amd_y_coeff[9] * 2.0 * xy +
                    wcs.amd_y_coeff[10] * 3.0 * x2 + wcs.amd_y_coeff[11] * 2.0 * xy +
                    wcs.amd_y_coeff[12] * 4.0 * xy * x2y2;
            /* magnitude and color terms ignored
            wcs.amd_y_coeff[17]*mag*2.0*xmm +
            wcs.amd_y_coeff[18]*mag*ymm*2.0*xmm;
            */

            /* Derivative of Y model wrt y */
            gy = wcs.amd_y_coeff[0] + wcs.amd_y_coeff[3] * 2.0 * ymm +
                    wcs.amd_y_coeff[4] * xmm + wcs.amd_y_coeff[6] * 2.0 * ymm +
                    wcs.amd_y_coeff[7] * 3.0 * y2 + wcs.amd_y_coeff[8] * 2.0 * xy +
                    wcs.amd_y_coeff[9] * x2 + wcs.amd_y_coeff[11] * (x2 + 3.0 * y2) +
                    wcs.amd_y_coeff[12] * (5.0 * y4 + 6.0 * x2 * y2 + x4);
            /* magnitude and color terms ignored
            wcs.amd_y_coeff[16]*mag       + wcs.amd_y_coeff[17]*mag*2.0*ymm +
            wcs.amd_y_coeff[18]*mag*(x2+3.0*y2);
            */

            f = f - xi;
            g = g - eta;
            dx = ((-f * gy) + (g * fy)) / ((fx * gy) - (fy * gx));
            dy = ((-g * fx) + (f * gx)) / ((fx * gy) - (fy * gx));
            xmm = xmm + dx;
            ymm = ymm + dy;

            if ((Math.abs(dx) < TOLERANCE) && (Math.abs(dy) < TOLERANCE)) break;
        }

        /* Convert mm from plate center to plate pixels */
        x = (wcs.ppo_coeff[2] - xmm * 1000.0) / wcs.x_pixel_size;
        y = (wcs.ppo_coeff[5] + ymm * 1000.0) / wcs.y_pixel_size;

        /* Convert from plate pixels to image pixels */
        double xpix = x - wcs.x_pixel_offset + 1.0 - 0.5;
        double ypix = y - wcs.y_pixel_offset + 1.0 - 0.5;

        /* If position is off of the image, return offscale code */
        // XXX if (xpix < 0.5 || xpix > wcs.nxpix+0.5)
        //	return null;
        // XXXif (ypix < 0.5 || ypix > wcs.nypix+0.5)
        //	return null;

        return new Point2D.Double(xpix, ypix);
    }

    /* Mar  6 1995	Original version of this code
       May  4 1995	Fix eta cross terms which were all in y
       Jun 21 1995	Add inverse routine
       Oct 17 1995	Fix inverse routine (degrees -> radians)
       Nov  7 1995	Add half pixel to image coordinates to get astrometric
                      plate coordinates
       Feb 26 1996	Fix plate to image pixel conversion error
     */
}