// Static SLALIB like routines using Eclipse
package net.mar;

import java.util.Arrays;

public class SLALIB {
    static final double d2s=86400.0;
    
    
    public static double DSEP (double A1, double B1, double A2, double B2) {
        /*
         * *+
*     - - - - -
*      D S E P
*     - - - - -
*
*  Angle between two points on a sphere.
*
*  (double precision)
*
*  Given:
*     A1,B1    d     spherical coordinates of one point
*     A2,B2    d     spherical coordinates of the other point
*
*  (The spherical coordinates are [RA,Dec], [Long,Lat] etc, in radians.)
*
*  The result is the angle, in radians, between the two points.  It
*  is always positive.
*
*  Called:  sla_DCS2C, sla_DSEPV
*
*  Last revision:   7 May 2000
*
*  Copyright P.T.Wallace.  All rights reserved.
*
*-
         */
      
        double[] V1 = new double[3];
        double[] V2 = new double[3];
 //  Convert coordinates from spherical to Cartesian.
        V1=SLALIB.DCS2C(A1,B1);
        V2=SLALIB.DCS2C(A2,B2);
        
        return SLALIB.DSEPV(V1,V2);
        
     
    }
    
    public static double DSEPV (double [] V1,double [] V2) {
        /*
         * *+
*     - - - - - -
*      D S E P V
*     - - - - - -
*
*  Angle between two vectors.
*
*  (double precision)
*
*  Given:
*     V1      d(3)    first vector
*     V2      d(3)    second vector
*
*  The result is the angle, in radians, between the two vectors.  It
*  is always positive.
*
*  Notes:
*
*  1  There is no requirement for the vectors to be unit length.
*
*  2  If either vector is null, zero is returned.
*
*  3  The simplest formulation would use dot product alone.  However,
*     this would reduce the accuracy for angles near zero and pi.  The
*     algorithm uses both cross product and dot product, which maintains
*     accuracy for all sizes of angle.
*
*  Called:  sla_DVXV, sla_DVN, sla_DVDV
         */
        double[] V1XV2= new double[3];
        double[] WV= new double[4];
        double S,C ; //sla_DVDV;
        double sla_DSEPV;
        V1XV2=SLALIB.DVXV(V1,V2);
        WV=SLALIB.DVN(V1XV2); // Wv[3] is S
        S=WV[3];
        //*  Dot product = cosine multiplied by the two moduli.
        C = SLALIB.DVDV(V1,V2);
        //*  Angle between the vectors.
        if (S != 0.0) {
           sla_DSEPV = Math.atan2(S,C);
        }
        else {
           sla_DSEPV = 0.0;
        }
        return sla_DSEPV;
        
        
    }
    
    public static double DVDV (double [] VA, double [] VB) {
        /*
         * *+
*     - - - - -
*      D V D V
*     - - - - -
*
*  Scalar product of two 3-vectors  (double precision)
*
*  Given:
*      VA      dp(3)     first vector
*      VB      dp(3)     second vector
*
*  The result is the scalar product VA.VB (double precision)
*
*  P.T.Wallace   Starlink   November 1984
*
*  Copyright (C) 1995 Rutherford Appleton Laboratory
*
         */
        return VA[0]*VB[0]+VA[1]*VB[1]+VA[2]*VB[2];
        
    }
    public static double[] DVN(double [] V) {
        /*
         * *+
*     - - - -
*      D V N
*     - - - -
*
*  Normalizes a 3-vector also giving the modulus (double precision)
*
*  Given:
*     V       dp(3)      vector
*
*  Returned:
*     UV      dp(3)      unit vector in direction of V
*     VM      dp         modulus of V
*
*  If the modulus of V is zero, UV is set to zero as well
*
*  P.T.Wallace   Starlink   23 November 1995
         */
        double [] UV = new double [4]; // 4 not 3 
        double VM;
        int I;
        double W1,W2;
   //*  Modulus
        W1=0.0;
        for (I=0;I<3;I++) {
          W2=V[I];
          W1=W1+W2*W2;
        }     
        W1=Math.sqrt(W1);
        VM=W1;
 // *  Normalize the vector
        if (W1 <= 0.0) {
            W1=1.0;
        }
        for (I=0;I<3;I++) {
           UV[I]=V[I]/W1;
        }
        UV[3]=VM;
        return UV;
 
        
    }
    public static double [] DVXV(double[] VA, double [] VB) {
        /*
         * *+
*     - - - - -
*      D V X V
*     - - - - -
*
*  Vector product of two 3-vectors  (double precision)
*
*  Given:
*      VA      dp(3)     first vector
*      VB      dp(3)     second vector
*
*  Returned:
*      VC      dp(3)     vector result
*
*  P.T.Wallace   Starlink   March 1986
*
*  Copyright (C) 1995 Rutherford Appleton Laboratory
*
         */
       
        double [] VC=new double [3];
        double [] VW=new double [3];
        int I;


  //*  Form the vector product VA cross VB
        VW[0]=VA[1]*VB[2]-VA[2]*VB[1];
        VW[1]=VA[2]*VB[0]-VA[0]*VB[2];
        VW[2]=VA[0]*VB[1]-VA[1]*VB[0];

  //*  Return the result
        for (I=0; I<3;I++) {     
           VC[I]=VW[I];
        }
        return VC;
    }
    
    
    
    
	public static double[] DCS2C(double A, double B) {
/*
 * * Spherical coordinates to direction cosines (double precision)
 * 
 * Given: A,B dp spherical coordinates in radians (RA,Dec), (Long,Lat) etc
 * 
 * Returned: V dp(3) x,y,z unit vector
 * 
 * The spherical coordinates are longitude (+ve anticlockwise looking from the
 * +ve latitude pole) and latitude. The Cartesian coordinates are right handed,
 * with the x axis at zero longitude and latitude, and the z axis at the +ve
 * latitude pole.
 *  
 */
	    double[] V = new double[3];
		double COSB = Math.cos(B);
		V[0] = Math.cos(A) * COSB;
		V[1] = Math.sin(A) * COSB;
		V[2] = Math.sin(B);
		return V;
	}
	

	public static double[] DCC2S(double[] V) {
/*
 * * Direction cosines to spherical coordinates (double precision)
 * 
 * Given: V d(3) x,y,z vector
 * 
 * Returned: A,B d spherical coordinates in radians
 * 
 * The spherical coordinates are longitude (+ve anticlockwise looking from the
 * +ve latitude pole) and latitude. The Cartesian coordinates are right handed,
 * with the x axis at zero longitude and latitude, and the z axis at the +ve
 * latitude pole.
 * 
 * If V is null, zero A and B are returned. At either pole, zero A is returned.
 *  
 */
	    double X = V[0];
		double Y = V[1];
		double Z = V[2];
		double A;
		double B;
		double R = Math.sqrt(X * X + Y * Y);
		if (R == 0.0) {
			A = 0.0;
		} else {
			A = Math.atan2(Y, X);
		}
		if (Z == 0.0) {
			B = 0.0;
		} else {
			B = Math.atan2(Z, R);
		}

		double[] coords = new double[2];
		coords[0] = A;
		coords[1] = B;
		return coords;
	}

	public static double[] FK45Z(double R1950, double D1950, double BEPOCH) {
/*
 * * Convert B1950.0 FK4 star data to J2000.0 FK5 assuming zero proper motion in
 * the FK5 frame (double precision)
 * 
 * This routine converts stars from the old, Bessel-Newcomb, FK4 system to the
 * new, IAU 1976, FK5, Fricke system, in such a way that the FK5 proper motion
 * is zero. Because such a star has, in general, a non-zero proper motion in the
 * FK4 system, the routine requires the epoch at which the position in the FK4
 * system was determined.
 * 
 * The method is from Appendix 2 of Ref 1, but using the constants of Ref 4.
 * 
 * Given: R1950,D1950 dp B1950.0 FK4 RA,Dec at epoch (rad) BEPOCH dp Besselian
 * epoch (e.g. 1979.3D0)
 * 
 * Returned: R2000,D2000 dp J2000.0 FK5 RA,Dec (rad)
 * 
 * Notes:
 * 
 * 1) The epoch BEPOCH is strictly speaking Besselian, but if a Julian epoch is
 * supplied the result will be affected only to a negligible extent.
 * 
 * 2) Conversion from Besselian epoch 1950.0 to Julian epoch 2000.0 only is
 * provided for. Conversions involving other epochs will require use of the
 * appropriate precession, proper motion, and E-terms routines before and/or
 * after FK45Z is called.
 * 
 * 3) In the FK4 catalogue the proper motions of stars within 10 degrees of the
 * poles do not embody the differential E-term effect and should, strictly
 * speaking, be handled in a different manner from stars outside these regions.
 * However, given the general lack of homogeneity of the star data available for
 * routine astrometry, the difficulties of handling positions that may have been
 * determined from astrometric fields spanning the polar and non-polar regions,
 * the likelihood that the differential E-terms effect was not taken into
 * account when allowing for proper motion in past astrometry, and the
 * undesirability of a discontinuity in the algorithm, the decision has been
 * made in this routine to include the effect of differential E-terms on the
 * proper motions for all stars, whether polar or not. At epoch 2000, and
 * measuring on the sky rather than in terms of dRA, the errors resulting from
 * this simplification are less than 1 milliarcsecond in position and 1
 * milliarcsecond per century in proper motion.
 * 
 * References:
 * 
 * 1 Aoki,S., et al, 1983. Astron.Astrophys., 128, 263.
 * 
 * 2 Smith, C.A. et al, 1989. "The transformation of astrometric catalog systems
 * to the equinox J2000.0". Astron.J. 97, 265.
 * 
 * 3 Yallop, B.D. et al, 1989. "Transformation of mean star places from FK4
 * B1950.0 to FK5 J2000.0 using matrices in 6-space". Astron.J. 97, 274.
 * 
 * 4 Seidelmann, P.K. (ed), 1992. "Explanatory Supplement to the Astronomical
 * Almanac", ISBN 0-935702-68-7.
 * 
 * Called: sla_DCS2C, sla_EPJ, sla_EPB2D, sla_DCC2S, sla_DRANRM
 */
	    
	    
	    
	    double R2000;
		double D2000;
		double D2PI = 2.0 * Math.PI;
		double PMF = 100.0 * 60.0 * 60.0 * 360.0 / D2PI;
		double W;
		int I, J;
		double[] R0 = new double[3];
		double[] A1 = new double[3];
		double[] V1 = new double[3];
		double[] V2 = new double[6];
		double[] VECT = new double[3];
		double[] A = { -1.62557E-6, -0.31919E-6, -0.13843E-6 };
		double[] AD = { +1.245E-3, -1.580E-3, -0.659E-3 };
		double[][] EM = new double[6][3];
		EM[0][0] = +0.9999256782;
		EM[1][0] = +0.0111820610;
		EM[2][0] = +0.0048579479;
		EM[3][0] = -0.000551;
		EM[4][0] = +0.238514;
		EM[5][0] = -0.435623;

		EM[0][1] = -0.0111820611;
		EM[1][1] = +0.9999374784;
		EM[2][1] = -0.0000271474;
		EM[3][1] = -0.238565;
		EM[4][1] = -0.002667;
		EM[5][1] = +0.012254;

		EM[0][2] = -0.0048579477;
		EM[1][2] = -0.0000271765;
		EM[2][2] = +0.9999881997;
		EM[3][2] = +0.435739;
		EM[4][2] = -0.008541;
		EM[5][2] = +0.002117;

		R0 = DCS2C(R1950, D1950);

		W = (BEPOCH - 1950.0) / PMF;

		for (I = 0; I < 3; I++) {
			A1[I] = A[I] + W * AD[I];

		}

		W = R0[0] * A1[0] + R0[1] * A1[1] + R0[2] * A1[2];

		for (I = 0; I < 3; I++) {
			V1[I] = R0[I] - A1[I] + W * R0[I];

		}
		for (I = 0; I < 6; I++) {
			W = 0.0;
			for (J = 0; J < 3; J++) {
				W = W + EM[I][J] * V1[J];

			}
			V2[I] = W;
		}

		W = (EPJ(EPB2D(BEPOCH)) - 2000.0) / PMF;
		for (I = 0; I < 3; I++) {
			V2[I] = V2[I] + W * V2[I + 3];

		}

		double[] temp = new double[2];
		temp = DCC2S(V2);
		D2000 = temp[1];
		R2000 = DRANRM(temp[0]);

		double[] coords = new double[2];
		coords[0] = R2000;
		coords[1] = D2000;
		return coords;
	}

	public static double[] EQGAL(double DR, double DD) {
/*
 * * Transformation from J2000.0 equatorial coordinates to IAU 1958 galactic
 * coordinates (double precision)
 * 
 * Given: DR,DD dp J2000.0 RA,Dec
 * 
 * Returned: DL,DB dp galactic longitude and latitude L2,B2
 * 
 * (all arguments are radians)
 * 
 * Called: sla_DCS2C, sla_DMXV, sla_DCC2S, sla_DRANRM, sla_DRANGE
 * 
 * Note: The equatorial coordinates are J2000.0. Use the routine sla_EG50 if
 * conversion from B1950.0 'FK4' coordinates is required.
 */
	    double[][] RMAT = new double[3][3];
		RMAT[0][0] = -0.054875539726;
		RMAT[0][1] = -0.873437108010;
		RMAT[0][2] = -0.483834985808;
		RMAT[1][0] = +0.494109453312;
		RMAT[1][1] = -0.444829589425;
		RMAT[1][2] = +0.746982251810;
		RMAT[2][0] = -0.867666135858;
		RMAT[2][1] = -0.198076386122;
		RMAT[2][2] = +0.455983795705;
		double[] V1 = new double[3];
		double[] V2 = new double[3];
		V1 = SLALIB.DCS2C(DR, DD);
		V2 = SLALIB.DMXV(RMAT, V1);
		double[] temp = new double[2];
		temp = DCC2S(V2);
		double[] coords = new double[2];
		coords[0] = SLALIB.DRANRM(temp[0]);
		coords[1] = SLALIB.DRANGE(temp[1]);
		return coords;

	}

	public static double[] DMXV(double[][] DM, double[] VA) {
/*
 * * Performs the 3-D forward unitary transformation:
 * 
 * vector VB = matrix DM * vector VA
 * 
 * (double precision)
 * 
 * Given: DM dp(3,3) matrix VA dp(3) vector
 * 
 * Returned: VB dp(3) result vector
 *  
 */
	    double W;
		double[] VW = new double[3];
		double[] VB = new double[3];
		for (int J = 0; J < 3; J++) {
			W = 0.0;
			for (int I = 0; I < 3; I++) {
				W = W + DM[J][I] * VA[I];
			}
			VW[J] = W;
		}
		for (int J = 0; J < 3; J++) {
			VB[J] = VW[J];
		}
		return VB;
	}

	public static double[] DS2TP(double RA, double DEC, double RAZ, double DECZ) {
/*
 * * Projection of spherical coordinates onto tangent plane: "gnomonic"
 * projection - "standard coordinates" (double precision)
 * 
 * Given: RA,DEC dp spherical coordinates of point to be projected RAZ,DECZ dp
 * spherical coordinates of tangent point
 * 
 * Returned: XI,ETA dp rectangular coordinates on tangent plane J int status: 0 =
 * OK, star on tangent plane 1 = error, star too far from axis 2 = error,
 * antistar on tangent plane 3 = error, antistar too far from axis
 *  
 */
	    double[] results = new double[3];
		double TINY = 1E-6;
		double SDECZ = Math.sin(DECZ);
		double SDEC = Math.sin(DEC);
		double CDECZ = Math.cos(DECZ);
		double CDEC = Math.cos(DEC);
		double RADIF = RA - RAZ;
		double SRADIF = Math.sin(RADIF);
		double CRADIF = Math.cos(RADIF);
		double DENOM = SDEC * SDECZ + CDEC * CDECZ * CRADIF;
		if (DENOM > TINY) {
			results[2] = 0.0;
		} else if (DENOM >= 0.0) {
			results[2] = 1.0;
			DENOM = TINY;
			
		} else if (DENOM > -TINY) {
			results[2] = 2.0;
			DENOM = -TINY;
			
		} else {
			results[2] = 3.0;
		}
		
		results[0] = CDEC * SRADIF / DENOM; //XI
		results[1] = (SDEC * CDECZ - CDEC * SDECZ * CRADIF) / DENOM; //ETA

		return results;
	}

	public static double[] DTP2S(double XI, double ETA, double RAZ, double DECZ) {
    
    /*
     * + - - - - - - D T P 2 S - - - - - -
     * 
     * Transform tangent plane coordinates into spherical (double precision)
     * 
     * Given: XI,ETA dp tangent plane rectangular coordinates RAZ,DECZ dp
     * spherical coordinates of tangent point
     * 
     * Returned: RA,DEC dp spherical coordinates (0-2pi,+/-pi/2)
     * 
     * Called: sla_DRANRM
     * 
     * P.T.Wallace Starlink 24 July 1995
     *  
     */

	    double[] coords = new double[2];
        double SDECZ=Math.sin(DECZ);
        double CDECZ=Math.cos(DECZ);
        double DENOM=CDECZ-ETA*SDECZ;
          
        coords[0]=SLALIB.DRANRM(Math.atan2(XI,DENOM)+RAZ);
        coords[1]=Math.atan2(SDECZ+ETA*CDECZ,Math.sqrt(XI*XI+DENOM*DENOM));
        return coords;
	}
	
	
	
	
	public static double EPB2D(double EPB) {
/*
 * * Conversion of Besselian Epoch to Modified Julian Date (double precision)
 * 
 * Given: EPB dp Besselian Epoch
 * 
 * The result is the Modified Julian Date (JD - 2400000.5).
 */
	    return 15019.81352 + (EPB - 1900.0) * 365.242198781;
	}

	public static double EPJ(double DATE) {
/*
 * * Conversion of Modified Julian Date to Julian Epoch (double precision)
 * 
 * Given: DATE dp Modified Julian Date (JD - 2400000.5)
 * 
 * The result is the Julian Epoch.
 */
	    return 2000.0 + (DATE - 51544.5) / 365.25;
	}

	public static double DRANRM(double ANGLE) {
/*
 * * Normalize angle into range 0-2 pi (double precision)
 * 
 * Given: ANGLE dp the angle in radians
 * 
 * The result is ANGLE expressed in the range 0-2 pi (double precision).
 */
	    double D2PI = Math.PI * 2;
		double dranrm = ANGLE % D2PI;
		if (dranrm < 0) {
			dranrm = dranrm + D2PI;
		}
		return dranrm;
	}

	public static double DRANGE(double ANGLE) {
/*
 * * Normalize angle into range +/- pi (double precision)
 * 
 * Given: ANGLE dp the angle in radians
 * 
 * The result (double precision) is ANGLE expressed in the range +/- pi.
 */
	    double D2PI = Math.PI * 2;
		double DPI = Math.PI;
		double drange = ANGLE % D2PI;
		if (Math.abs(drange) >= DPI) {
			if (ANGLE >= 0.0) {
				drange = drange - D2PI;
			} else {
				drange = drange + D2PI;
			}
		}
		return drange;
	}
	
	public static double[] GALEQ(double DL,double DB) {
	    
    /*
    *+
    *     - - - - - -
    *      G A L E Q
    *     - - - - - -
    *
    *  Transformation from IAU 1958 galactic coordinates to
    *  J2000.0 equatorial coordinates (double precision)
    *
    *  Given:
    *     DL,DB       dp       galactic longitude and latitude L2,B2
    *
    *  Returned:
    *     DR,DD       dp       J2000.0 RA,Dec
    *
    *  (all arguments are radians)
    *
    *  Called:
    *     sla_DCS2C, sla_DIMXV, sla_DCC2S, sla_DRANRM, sla_DRANGE
    *
    *  Note:
    *     The equatorial coordinates are J2000.0.  Use the routine
    *     sla_GE50 if conversion to B1950.0 'FK4' coordinates is
    *     required.
    *
    *  Reference:
    *     Blaauw et al, Mon.Not.R.Astron.Soc.,121,123 (1960)
    *
    *  P.T.Wallace   Starlink   21 September 1998
    *
    *  Copyright (C) 1998 Rutherford Appleton Laboratory
    *
    *
    *  L2,B2 system of galactic coordinates
    *
    *  P = 192.25       RA of galactic north pole (mean B1950.0)
    *  Q =  62.6        inclination of galactic to mean B1950.0 equator
    *  R =  33          longitude of ascending node
    *
    *  P,Q,R are degrees
    *
    *  Equatorial to galactic rotation matrix (J2000.0), obtained by
    *  applying the standard FK4 to FK5 transformation, for zero proper
    *  motion in FK5, to the columns of the B1950 equatorial to
    *  galactic rotation matrix:
    *
    **/
	    double[][] RMAT = new double[3][3];
		RMAT[0][0] = -0.054875539726;
		RMAT[0][1] = -0.873437108010;
		RMAT[0][2] = -0.483834985808;
		RMAT[1][0] = +0.494109453312;
		RMAT[1][1] = -0.444829589425;
		RMAT[1][2] = +0.746982251810;
		RMAT[2][0] = -0.867666135858;
		RMAT[2][1] = -0.198076386122;
		RMAT[2][2] = +0.455983795705;
		double[] V1 = new double[3];
		double[] V2 = new double[3];

    //  Spherical to Cartesian
		V1=SLALIB.DCS2C(DL,DB);
    //  Galactic to equatorial		
        V2=SLALIB.DIMXV(RMAT,V1);
		double[] temp = new double[2];
    //  Cartesian to spherical
        temp = DCC2S(V2); 
		double[] coords = new double[2];
    //  Express in conventional ranges
		coords[0] = SLALIB.DRANRM(temp[0]);
		coords[1] = SLALIB.DRANGE(temp[1]);
		return coords;         
	}
	public static double[] DIMXV(double[][] DM, double[] VA) {
	  /*
	   * *+
*     - - - - - -
*      D I M X V
*     - - - - - -
*
*  Performs the 3-D backward unitary transformation:
*
*     vector VB = (inverse of matrix DM) * vector VA
*
*  (double precision)
*
*  (n.b.  the matrix must be unitary, as this routine assumes that
*   the inverse and transpose are identical)
*
*  Given:
*     DM       dp(3,3)    matrix
*     VA       dp(3)      vector
*
*  Returned:
*     VB       dp(3)      result vector
* 
*/
	    double W;
	    double[] VW = new double[3];
		double[] VB = new double[3];
	    for (int J = 0; J < 3; J++) {
	        W=0.0;
	        for (int I = 0; I < 3; I++) {
	            W=W+DM[I][J]*VA[I];
	        }
	        VW[J]=W;  
	    }
		for (int J = 0; J < 3; J++) {
			VB[J] = VW[J];
		}
		return VB;
	}

	public static String ToSexagesimal (double angle, int ndp){
	    return ToSexagesimal (angle,ndp,true);
	}
	public static String ToSexagesimal (double angle, int ndp, boolean isRA) {
	    
	/*
    SUBROUTINE sla_DD2TF (NDP, DAYS, SIGN, IHMSF)
    *+
    *     - - - - - -
    *      D D 2 T F
    *     - - - - - -
    *
    *  Convert an interval in days into hours, minutes, seconds
    *  (double precision)
    *
    *  Given:
    *     NDP      i      number of decimal places of seconds
    *     DAYS     d      interval in days
    *
    *  Returned:
    *     SIGN     c      '+' or '-'
    *     IHMSF    i(4)   hours, minutes, seconds, fraction
    *
    *  Notes:
    *
    *     1)  NDP less than zero is interpreted as zero.
    *
    *     2)  The largest useful value for NDP is determined by the size
    *         of DAYS, the format of DOUBLE PRECISION floating-point numbers
    *         on the target machine, and the risk of overflowing IHMSF(4).
    *         For example, on the VAX, for DAYS up to 1D0, the available
    *         floating-point precision corresponds roughly to NDP=12.
    *         However, the practical limit is NDP=9, set by the capacity of
    *         the 32-bit integer IHMSF(4).
    *
    *     3)  The absolute value of DAYS may exceed 1D0.  In cases where it
    *         does not, it is up to the caller to test for and handle the
    *         case where DAYS is very nearly 1D0 and rounds up to 24 hours,
    *         by testing for IHMSF(1)=24 and setting IHMSF(1-4) to zero.
    *
    *  P.T.Wallace   Starlink   19 March 1999
    *
    *  Copyright (C) 1999 Rutherford Appleton Laboratory
    *
    *  License:
    *    This program is free software; you can redistribute it and/or modify
    *    it under the terms of the GNU General Public License as published by
    *    the Free Software Foundation; either version 2 of the License, or
    *    (at your option) any later version.
    *
    *    This program is distributed in the hope that it will be useful,
    *    but WITHOUT ANY WARRANTY; without even the implied warranty of
    *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    *    GNU General Public License for more details.
    *
    *    You should have received a copy of the GNU General Public License
    *    along with this program (see SLA_CONDITIONS); if not, write to the 
    *    Free Software Foundation, Inc., 59 Temple Place, Suite 330, 
    *    Boston, MA  02111-1307  USA
    *
    *-
        */
	    
	    
	    double days = angle/(2.0*Math.PI);   
	    if (!isRA) {
	        days=days*15.0;
	    }
        char sign;
        int nrs;
        double rs,rm,rh,a,ah,am,as,af;
        String [] strArray = new String [4];
        
       /*
          INTEGER NDP
          DOUBLE PRECISION DAYS
          CHARACTER SIGN*(*)
          INTEGER IHMSF(4)

    *  Days to seconds
          DOUBLE PRECISION D2S
          PARAMETER (D2S=86400D0)

          INTEGER NRS,N
          DOUBLE PRECISION RS,RM,RH,A,AH,AM,AS,AF
*/


    //  Handle sign
          if (days >= 0.0) {
             sign='+';
          }else {
             sign='-';
          }
          

    //  Field units in terms of least significant figure
          nrs=1;
          //for (I = 0; I < 3; I++) {
          for (int n=1; n <= ndp; n++) {
              nrs=nrs*10;
          }
          
          rs=(double)nrs;
          rm=rs*60.0;
          rh=rm*60.0;
          
   //  Round interval and express in smallest units required
          a=(double)Math.round(rs*d2s*Math.abs(days));
   //       A=ANINT(RS*D2S*ABS(DAYS))

   //  Separate into fields
          ah=Math.floor(a/rh); // AH=AINT(A/RH)
          a=a-ah*rh; //A=A-AH*RH
          am=Math.floor(a/rm); // AM=AINT(A/RM)
          a=a-am*rm; // A=A-AM*RM
          as=Math.floor(a/rs); // AS=AINT(A/RS)
          af=a-as*rs; // AF=A-AS*RS

    //  Return results
          StringBuffer sb =new StringBuffer();
          sb.append(sign);
          String tempString=String.valueOf(Math.max(Math.round(ah),0)); //IHMSF(1)=MAX(NINT(AH),0)
          int l=tempString.length();
          for (int i=l;i<2;i++) {
              sb.append('0');
          }
          sb.append(tempString);  
          sb.append(':');
          tempString=String.valueOf(Math.max(Math.min(Math.round(am),59),0)); //IHMSF(2)=MAX(MIN(NINT(AM),59),0)
          l=tempString.length();
          for (int i=l;i<2;i++) {
              sb.append('0');
          }
          sb.append(tempString);
          sb.append(':');
          tempString=String.valueOf(Math.max(Math.min(Math.round(as),59),0)); //IHMSF(3)=MAX(MIN(NINT(AS),59),0)
          l=tempString.length();
          for (int i=l;i<2;i++) {
              sb.append('0');
          }
          sb.append(tempString); 
          sb.append('.');
          
          tempString=String.valueOf(Math.max(Math.round(Math.min(af,rs-1.0)),0)); //IHMSF(4)=MAX(NINT(MIN(AF,RS-1D0)),0)
          l=tempString.length();
          for (int i=l;i<ndp;i++) {
              sb.append('0');
          }
          sb.append(tempString);
          return sb.toString();
         
	}
	public static double[] PCD (double disco,double x, double y) {
/*	      SUBROUTINE sla_PCD (DISCO,X,Y)
	      *+
	      *     - - - -
	      *      P C D
	      *     - - - -
	      *
	      *  Apply pincushion/barrel distortion to a tangent-plane [x,y].
	      *
	      *  Given:
	      *     DISCO    d      pincushion/barrel distortion coefficient
	      *     X,Y      d      tangent-plane coordinates
	      *
	      *  Returned:
	      *     X,Y      d      distorted coordinates
	      *
	      *  Notes:
	      *
	      *  1)  The distortion is of the form RP = R*(1 + C*R**2), where R is
	      *      the radial distance from the tangent point, C is the DISCO
	      *      argument, and RP is the radial distance in the presence of
	      *      the distortion.
	      *
	      *  2)  For pincushion distortion, C is +ve;  for barrel distortion,
	      *      C is -ve.
	      *
	      *  3)  For X,Y in units of one projection radius (in the case of
	      *      a photographic plate, the focal length), the following
	      *      DISCO values apply:
	      *
	      *          Geometry          DISCO
	      *
	      *          astrograph         0.0
	      *          Schmidt           -0.3333
	      *          AAT PF doublet  +147.069
	      *          AAT PF triplet  +178.585
	      *          AAT f/8          +21.20
	      *          JKT f/8          +13.32
	      *
	      *  4)  There is a companion routine, sla_UNPCD, which performs the
	      *      inverse operation.
	      *
	      *  P.T.Wallace   Starlink   3 September 2000
	      *
	      *  Copyright (C) 2000 Rutherford Appleton Laboratory
	      *
	      *  License:
	      *    This program is free software; you can redistribute it and/or modify
	      *    it under the terms of the GNU General Public License as published by
	      *    the Free Software Foundation; either version 2 of the License, or
	      *    (at your option) any later version.
	      *
	      *    This program is distributed in the hope that it will be useful,
	      *    but WITHOUT ANY WARRANTY; without even the implied warranty of
	      *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
	      *    GNU General Public License for more details.
	      *
	      *    You should have received a copy of the GNU General Public License
	      *    along with this program (see SLA_CONDITIONS); if not, write to the 
	      *    Free Software Foundation, Inc., 59 Temple Place, Suite 330, 
	      *    Boston, MA  02111-1307  USA
	      *
	      *-
*/
	            double f;
	            double[] xy = new double[2];
	            f=1.0+disco*(x*x+y*y);
	            xy[0]=x*f;
	            xy[1]=y*f; 
//	            int i = Arrays.binarySearch(xy, 1.1);
	            return xy;
	    
	}
	
	public static double[] XY2XY (double x1,double y1, double [] coeffs) {
	    
	/*
    SUBROUTINE sla_XY2XY (X1,Y1,COEFFS,X2,Y2)
    *+
    *     - - - - - -
    *      X Y 2 X Y
    *     - - - - - -
    *
    *  Transform one [X,Y] into another using a linear model of the type
    *  produced by the sla_FITXY routine.
    *
    *  Given:
    *     X1       d        x-coordinate
    *     Y1       d        y-coordinate
    *     COEFFS  d(6)      transformation coefficients (see note)
    *
    *  Returned:
    *     X2       d        x-coordinate
    *     Y2       d        y-coordinate
    *
    *  The model relates two sets of [X,Y] coordinates as follows.
    *  Naming the elements of COEFFS:
    *
    *     COEFFS(1) = A
    *     COEFFS(2) = B
    *     COEFFS(3) = C
    *     COEFFS(4) = D
    *     COEFFS(5) = E
    *     COEFFS(6) = F
    *
    *  the present routine performs the transformation:
    *
    *     X2 = A + B*X1 + C*Y1
    *     Y2 = D + E*X1 + F*Y1
    *
    *  See also sla_FITXY, sla_PXY, sla_INVF, sla_DCMPF
    *
    *  P.T.Wallace   Starlink   5 December 1994
    *
    *  Copyright (C) 1995 Rutherford Appleton Laboratory
    *
    *  License:
    *    This program is free software; you can redistribute it and/or modify
    *    it under the terms of the GNU General Public License as published by
    *    the Free Software Foundation; either version 2 of the License, or
    *    (at your option) any later version.
    *
    *    This program is distributed in the hope that it will be useful,
    *    but WITHOUT ANY WARRANTY; without even the implied warranty of
    *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    *    GNU General Public License for more details.
    *
    *    You should have received a copy of the GNU General Public License
    *    along with this program (see SLA_CONDITIONS); if not, write to the 
    *    Free Software Foundation, Inc., 59 Temple Place, Suite 330, 
    *    Boston, MA  02111-1307  USA
    *
    *-
    *
    */

	    double[] xy2 = new double[2];
         xy2[0]=coeffs[0]+coeffs[1]*x1+coeffs[2]*y1;
         xy2[1]=coeffs[3]+coeffs[4]*x1+coeffs[5]*y1;
          return xy2;

	}

	public static double[] INVF (double [] FWDS){
	   
	/*
	      SUBROUTINE sla_INVF (FWDS,BKWDS,J)
*+
*     - - - - -
*      I N V F
*     - - - - -
*
*  Invert a linear model of the type produced by the
*  sla_FITXY routine.
*
*  Given:
*     FWDS    d(6)      model coefficients
*
*  Returned:
*     BKWDS   d(6)      inverse model
*     J        i        status:  0 = OK, -1 = no inverse
*
*  The models relate two sets of [X,Y] coordinates as follows.
*  Naming the elements of FWDS:
*
*     FWDS(1) = A
*     FWDS(2) = B
*     FWDS(3) = C
*     FWDS(4) = D
*     FWDS(5) = E
*     FWDS(6) = F
*
*  where two sets of coordinates [X1,Y1] and [X2,Y1] are related
*  thus:
*
*     X2 = A + B*X1 + C*Y1
*     Y2 = D + E*X1 + F*Y1
*
*  the present routine generates a new set of coefficients:
*
*     BKWDS(1) = P
*     BKWDS(2) = Q
*     BKWDS(3) = R
*     BKWDS(4) = S
*     BKWDS(5) = T
*     BKWDS(6) = U
*
*  such that:
*
*     X1 = P + Q*X2 + R*Y2
*     Y1 = S + T*X2 + U*Y2
*
*  Two successive calls to sla_INVF will thus deliver a set
*  of coefficients equal to the starting values.
*
*  To comply with the ANSI Fortran standard, FWDS and BKWDS must
*  not be the same array, even though the routine is coded to
*  work on the VAX and most other computers even if this rule
*  is violated.
*
*  See also sla_FITXY, sla_PXY, sla_XY2XY, sla_DCMPF
*
*  P.T.Wallace   Starlink   11 April 1990
*
*  Copyright (C) 1995 Rutherford Appleton Laboratory
*
*  License:
*    This program is free software; you can redistribute it and/or modify
*    it under the terms of the GNU General Public License as published by
*    the Free Software Foundation; either version 2 of the License, or
*    (at your option) any later version.
*
*    This program is distributed in the hope that it will be useful,
*    but WITHOUT ANY WARRANTY; without even the implied warranty of
*    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
*    GNU General Public License for more details.
*
*    You should have received a copy of the GNU General Public License
*    along with this program (see SLA_CONDITIONS); if not, write to the 
*    Free Software Foundation, Inc., 59 Temple Place, Suite 330, 
*    Boston, MA  02111-1307  USA
*
*-

      IMPLICIT NONE
*/
      double [] BKWDS = new double[6];
      int J;

      double A,B,C,D,E,F,DET;



      A=FWDS[0];
      B=FWDS[1];
      C=FWDS[2];
      D=FWDS[3];
      E=FWDS[4];
      F=FWDS[5];
      DET=B*F-C*E;
      if (DET != 0.0) {
         BKWDS[0]=(C*D-A*F)/DET;
         BKWDS[1]=F/DET;
         BKWDS[2]=-C/DET;
         BKWDS[3]=(A*E-B*D)/DET;
         BKWDS[4]=-E/DET;
         BKWDS[5]=B/DET;
         J=0;
      }
      else {
         J=-1;
      }

      return BKWDS;
	
	}
	
	public static double [] UNPCD (double DISCO, double X, double Y) {
	    
	/*
	      SUBROUTINE sla_UNPCD ( DISCO, X, Y )
*+
*     - - - - - -
*      U N P C D
*     - - - - - -
*
*  Remove pincushion/barrel distortion from a distorted [x,y] to give
*  tangent-plane [x,y].
*
*  Given:
*     DISCO    d      pincushion/barrel distortion coefficient
*     X,Y      d      distorted coordinates
*
*  Returned:
*     X,Y      d      tangent-plane coordinates
*
*  Notes:
*
*  1)  The distortion is of the form RP = R*(1+C*R^2), where R is
*      the radial distance from the tangent point, C is the DISCO
*      argument, and RP is the radial distance in the presence of
*      the distortion.
*
*  2)  For pincushion distortion, C is +ve;  for barrel distortion,
*      C is -ve.
*
*  3)  For X,Y in "radians" - units of one projection radius,
*      which in the case of a photograph is the focal length of
*      the camera - the following DISCO values apply:
*
*          Geometry          DISCO
*
*          astrograph         0.0
*          Schmidt           -0.3333
*          AAT PF doublet  +147.069
*          AAT PF triplet  +178.585
*          AAT f/8          +21.20
*          JKT f/8          +13.32
*
*  4)  The present routine is a rigorous inverse of the companion
*      routine sla_PCD.  The expression for RP in Note 1 is rewritten
*      in the form x^3+a*x+b=0 and solved by standard techniques.
*
*  5)  Cases where the cubic has multiple real roots can sometimes
*      occur, corresponding to extreme instances of barrel distortion
*      where up to three different undistorted [X,Y]s all produce the
*      same distorted [X,Y].  However, only one solution is returned,
*      the one that produces the smallest change in [X,Y].
*
*  P.T.Wallace   Starlink   3 September 2000
*
*  Copyright (C) 2000 Rutherford Appleton Laboratory
*
*  License:
*    This program is free software; you can redistribute it and/or modify
*    it under the terms of the GNU General Public License as published by
*    the Free Software Foundation; either version 2 of the License, or
*    (at your option) any later version.
*
*    This program is distributed in the hope that it will be useful,
*    but WITHOUT ANY WARRANTY; without even the implied warranty of
*    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
*    GNU General Public License for more details.
*
*    You should have received a copy of the GNU General Public License
*    along with this program (see SLA_CONDITIONS); if not, write to the 
*    Free Software Foundation, Inc., 59 Temple Place, Suite 330, 
*    Boston, MA  02111-1307  USA
*
*-

      IMPLICIT NONE
*/
      

      double THIRD=1.0/3.0;
      
      double D2PI=2.0*Math.PI; 
      double RP,Q,R,D,W,S,T,F,C,T3,F1,F2,F3,W1,W2,W3;


      
//  Distance of the point from the origin.
      RP = Math.sqrt(X*X+Y*Y);

//  If zero, or if no distortion, no action is necessary.
      if (RP != 0.0 && DISCO != 0.0) {

//     Begin algebraic solution.
         Q = 1.0/(3.0*DISCO);
         R = RP/(2.0*DISCO);
         W = Q*Q*Q+R*R;

//     Continue if one real root, or three of which only one is positive.
         if (W >= 0.0) {
            D = Math.sqrt(W);
            W = R+D;
            S = SIGN(Math.pow(Math.abs(W),THIRD),W);
            W = R-D;
            T = SIGN(Math.pow((Math.abs(W)),THIRD),W);
            F = S+T;
         }
         else {

//        Three different real roots:  use geometrical method instead.
            W = 2.0/Math.sqrt(-3.00*DISCO);
            C = 4.0*RP/(DISCO*W*W*W);
            S = Math.sqrt(1.0-Math.min(C*C,1.0));
            T3 = Math.atan2(S,C);

//       The three solutions.
            F1 = W*Math.cos((D2PI-T3)/3.0);
            F2 = W*Math.cos((T3)/3.0);
            F3 = W*Math.cos((D2PI+T3)/3.0);

//       Pick the one that moves [X,Y] least.
            W1 = Math.abs(F1-RP);
            W2 = Math.abs(F2-RP);
            W3 = Math.abs(F3-RP);
            if (W1 < W2) {
               if (W1 < W3) {
                  F = F1;
               }
               else {
                  F = F3;
               }
            }
            else {
               if (W2 < W3) {
                  F = F2;
               }
               else {
                  F = F3;
               }
            }
 
         }

//     Remove the distortion.
         F = F/RP;
         X = F*X;
         Y = F*Y;

      }
      double [] xy={X,Y};
      return xy;
	}
         
         public static double SIGN(double a1, double a2) {
             if (a2 >= 0.0) {
                 return Math.abs(a1);
             }
             else {
                 return -Math.abs(a1);
             }
             
         }
}

