/*
 *  coord() is a function that generates equatorial (J2000) and
 *  galactic coordinates (II) from AMANDA coordinates.
 *
 *  by Patrick Mock, July 24, 1998
 *
 *  VARIABLE: UNITS   :  DEFINITION             :  I/O
 *  gpssec  : seconds :  seconds into the day   :  in
 *  gpsday  : days    :  days into the year     :  in
 *  gpsyear : years   :  year                   :  in
 *  th      : degrees :  theta in AMANDA        :  in
 *  ph      : degrees :  phi in AMANDA          :  out
 *  ra      : hours   :  right ascension J2000  :  out
 *  dec     : degrees :  declination J2000      :  out
 *  ll      : degrees :  galactic longitude II  :  out
 *  bb      : degrees :  galactic latitude II   :  out
 *
 *  AMANDA coordinate system is (x,y,z) <---> (r,theta,phi)
 *  where positive y = grid north = Greenwich Meridian
 *           "     x = grid east
 *           "     z = up (towards the equatorial south pole)
 *   theta = 0 is up
 *   phi   = 0 is grid east
 *
 *  DEFINITIONS:
 *  jdyr = (Julian date of 1200 UT 1/1/year)
 *         derived from Zombeck p107 equation for julian date
 *
 *  GMST from derived NOVAS package sidereal_time() function.
 *  NOVAS is online at http://aa.usno.navy.mil/aa/
 *
 *  galactic coordinates from J2000
 *  Reference:  Murray, C.A., Astronomy and Astrophysics (1989) 218:325.
 *
 *  Comment on precision:
 *  This function is accurate to ~1 arcminute.  You can improve the
 *  accuracy to ~0.01 arcseconds by switching to double precision
 *  floating point, incorporating the equation of the equinoxes, precession,
 *  and correcting for AMANDA's ~0.5 arcminute offset from the south pole.
 *  In other words the formulae and numerical constants are sufficent
 *  0.01 arcsecond accuracy.
 */
#include <math.h>

#define JDJ2000   (2451545)
#define JCENTURY  (36525.0)
#define RADTODEG  (57.29577951308232300l)
#define DEGTORAD  ( 0.01745329251994330l)
#define HOURTORAD ( 0.26179938779914941l)

void coord(float gpssec, int gpsday, int gpsyear, float th, float ph,
	   float *ra, float *dec, float *ll, float *bb)
{
  int jdyr, dayj2000, rai;
  float hour_angle, tdt, gmst, ra0;
  float rar,dcr,cosdec,eq1,eq2,eq3,gl1,gl2,gl3;
/*
 * Declination and Hour Angle
 */
  *dec       = th - 90.0;
  hour_angle = (ph - 90.0)/15.0;
/*
 * Julian Day
 * Note: GPS has no day zero 
 *       jdyr = (Julian date of 1200 UT 1/1/year)
 *       jd = gpssec/86400.0 + (jdyr + gpsday) - 1.5;  
 */
  jdyr = (1461*(gpsyear+4799))/4 - (3*((gpsyear+4899)/100))/4 - 31738;
  dayj2000 = jdyr + gpsday - JDJ2000;
/*
 * Sidereal Time in seconds
 * Note: tdt = (jd - JDJ2000)/JCENTURY;
 * This formula is accurate to 0.0001 s, if the calculation uses
 * double precison.
 */
  tdt  = ((float)dayj2000 + gpssec/86400.0 - 1.5)/JCENTURY;
  gmst = ((((-6.2e-6 * tdt + 0.093104) * tdt + 184.812866) * tdt
	   + 67310.54841) + 3.1644e9 * tdt);
/*
 * Right Ascension
 * Should use GAST, but (GAST-GMST) is small (<1.7 sec for 1990 - 2020)
 */
  ra0 = gmst/3600.0 - hour_angle;
  rai = (int) ra0;
  if      (rai < 0)  ra0 += (1-rai/24)*24;
  else if (rai > 24) ra0 -= (  rai/24)*24;
  *ra = ra0;
/*
 *  galactic coordinates from J2000
 *  Reference:  Murray, C.A., Astronomy and Astrophysics (1989) 218:325.
 */
  rar = ra0  * HOURTORAD;
  dcr = *dec * DEGTORAD;
  cosdec = cosf(dcr);
  eq1 = cosdec * cosf(rar);
  eq2 = cosdec * sinf(rar);
  eq3 = sinf(dcr);
  gl1 = -0.054875539 * eq1 + -0.873437105 * eq2 + -0.483834992 * eq3;
  gl2 =  0.494109454 * eq1 + -0.444829594 * eq2 +  0.746982249 * eq3;
  gl3 = -0.867666136 * eq1 + -0.198076390 * eq2 +  0.455983795 * eq3;
  *bb = asinf(gl3) * RADTODEG;
  *ll = atan2f(gl2,gl1) * RADTODEG;
}