// Copyright (C) 2007 Michael Kosowsky  All rights reserved.



// 1/2010  Rewrote to closer match the C++ classes in coords.h and sight.h,
//         and be more explciit about geocentric and vernal equinox based cartesian coordinates
//         Planet tranges and radii now in meters, not km
//
//    We have the notion of earth centered cartesian coordinates (geo_cartesian*)
//    where x-axis runs from center of earth through the equator at longitude 0, and
//    celestial cartesian coordinates where the x-axis runs to the vernal equinox.
//    You convert between them by rotating by GAST (greenwich apparent sidereal time);
//    at a given moment the celestial meridan
//    whose RA is GAST is lined up with terrestrial longitude 0.
//  
//    The POV matrix, which we copied from the C++ Sight class, converts back and
//    forth between a viewer's local coordinate system and geocentric cartesian.  (And
//    I think it handles a particularly tricky issue, that to look at the meridian
//    corresponding to my own longitude, I look south).
//  
//    The only gotcha is that azimuth and the spherical coordinate theta run in
//    opposite directions (from above, azimuth runs clockwise and theta ccw),
//    and in the C code all cartesian is earht based so we don't use the geo_ prefix
//


Sidereal = function(lat, lon, date) {
  this.lat  = lat;
  this.lon  = lon;
  this.date = date || new Date();
					// Greenwich apparent sidereal time
  this.gast = gast(this.date);
					// local apparent sidereal time
  this.last = normalize_hours(this.gast + this.lon/15);

  this.sinlat = sin_d(this.lat);
  this.coslat = cos_d(this.lat);
  this.sinlon = sin_d(this.lon);
  this.coslon = cos_d(this.lon);

  this.set_observer();
};


function ra_dec_to_gmaps_latlon(ra, dec) {
  return { lat: dec, lon: -(15 * ra - 180) };
}

function ra_dec_to_ge_latlon(ra, dec) {
  return { lat: dec, lon: 15 * ra - 180 };
}


Sidereal.prototype.DEGREES_PER_RADIAN = 180 / Math.PI;
Sidereal.prototype.HOURS_PER_RADIAN   =  12 / Math.PI;


function sin_d(x) {
    return Math.sin(x / Sidereal.prototype.DEGREES_PER_RADIAN);
}

function cos_d(x) {
    return Math.cos(x / Sidereal.prototype.DEGREES_PER_RADIAN);
}

function normalize_degrees(x) {
    y = x % 360 + (x < 0? 360 : 0);
    if (y >= 180) y -= 360;
    return y;
}

function sin_h(x) {
    return Math.sin(x / Sidereal.prototype.HOURS_PER_RADIAN);
}

function cos_h(x) {
    return Math.cos(x / Sidereal.prototype.HOURS_PER_RADIAN);
}

function normalize_hours(x) {
	// % returns value with same sign as left operand
    return x % 24 + (x < 0? 24 : 0);
}



Sidereal.prototype.ra_dec = function(az, alt) {
    var v = this.az_alt_to_geo_cartesian(az, alt);
    var o = cartesian_to_spherical(v[0], v[1], v[2]);
    return { ra:  o.theta * this.HOURS_PER_RADIAN + this.gast,
	     dec: o.phi   * this.DEGREES_PER_RADIAN            };
}


Sidereal.prototype.az_alt = function(ra, dec) {
   var c = spherical_to_cartesian(1, 15 * (ra - this.gast), dec);
   return this.geo_cartesian_to_az_alt(c[0], c[1], c[2]);
}




	// determine Greenwich apparent sidereal time
	// cf. aa.usno.navy.mil/faq/docs/GAST.html
	// gast seems to be within a second (of time, which is 1/4 minute of arc) of gast_precisely
	// see Sidereal.pm for mor precision

Sidereal.prototype.TIME_AT_JD_2451545_0 = 946728000;  // J2000.0, 1/1/2000 12:00 UTC in seconds

function gast(d) {
    return _gast(d.getTime()/1000);
}

function _gast(t) {
    return normalize_hours(18.697374558 + 24.06570982441908 * (t - Sidereal.prototype.TIME_AT_JD_2451545_0)/(24 * 60 * 60));
}


	// this is 0 to 180
function angular_distance_ra_dec(pra, pdec, qra, qdec) {
    var sinpdec = sin_d(pdec);
    var cospdec = cos_d(pdec);
    var sinqdec = sin_d(qdec);
    var cosqdec = cos_d(qdec);
    var cos_diff = Math.cos((pra - qra) / Sidereal.prototype.HOURS_PER_RADIAN) * cospdec * cosqdec + sinpdec * sinqdec;
    return Math.atan2(Math.sqrt(1 - cos_diff * cos_diff), cos_diff) * Sidereal.prototype.DEGREES_PER_RADIAN;
}

	// this is -180 to 180, where the sign is simply whether first object is ahead or behind in RA
function angular_difference_ra_dec(pra, pdec, qra, qdec) {
    return (pra >= qra? 1 : -1) * angular_distance_ra_dec(pra, pdec, qra, qdec);
}



    /* TOPOCENTRIC PLACE
       determine vector in earth-centered cartesian space from our position to the object,
       then use that vector to determine adjusted ra/dec

       BUG: doesn't take refraction into account, which is bigger than this geocentric parallax,
            but since it moves everything together the relative placements are right, and we can't change
            the base sky map
       BUG: doesn't use observer's elevation
    */

Sidereal.prototype.topocentric_ra_dec = function(ra, dec, range) {
  var c = this.ra_dec_to_geo_cartesian(ra, dec, range);
  var o = cartesian_to_spherical(c[0] - this.observer[0], c[1] - this.observer[1], c[2] - this.observer[2]);
  return { ra:    normalize_hours(o.theta * this.HOURS_PER_RADIAN + this.gast),
	   dec:   o.phi   * this.DEGREES_PER_RADIAN,
           range: o.r                                                           };
}


		// WGS 84: equatorial 6 378 137  polar 6 356 752.3142 m 
		//   for single radius we'll use ellipsoidal quadratic mean, cf. http://en.wikipedia.org/wiki/Spherical_earth
Sidereal.prototype.EARTH_RADIUS = 6372798.;
Sidereal.prototype.EARTH_SEMIMAJOR_AXIS = 6378137.;
Sidereal.prototype.EARTH_SEMIMINOR_AXIS = 6356752.3148;

Sidereal.prototype.set_observer = function() {
		// cf. coords.cpp: void EllipsoidalEarth::terrestrial_to_cartesian(const Terrestrial& t, Cartesian& c)
		// lat and lon are geodetic coordinates

  var den     = Math.sqrt(  this.EARTH_SEMIMAJOR_AXIS * this.EARTH_SEMIMAJOR_AXIS * this.coslat * this.coslat
			  + this.EARTH_SEMIMINOR_AXIS * this.EARTH_SEMIMINOR_AXIS * this.sinlat * this.sinlat);
  var rcosphi = this.coslat * (this.EARTH_SEMIMAJOR_AXIS * this.EARTH_SEMIMAJOR_AXIS / den  /* + this.elev */);

  var x = rcosphi * this.coslon;
  var y = rcosphi * this.sinlon;
  var z = this.sinlat * (this.EARTH_SEMIMINOR_AXIS * this.EARTH_SEMIMINOR_AXIS / den        /* + this.elev */ );

  this.observer = [ x, y, z];

  this.pov_matrix = [ -this.coslon * this.sinlat, -this.sinlon * this.sinlat, this.coslat,
                       this.sinlon,               -this.coslon,               0,
		       this.coslon * this.coslat,  this.sinlon * this.coslat, this.sinlat ]; 
}


	// equivalent to multiplying by transpose of pov matrix in sight.h
	// remember, azimuth is negative theta
Sidereal.prototype.az_alt_to_geo_cartesian = function(az, alt) {
    var v = spherical_to_cartesian(1, -az, alt);
    return [ this.pov_matrix[0] * v[0] + this.pov_matrix[3] * v[1] + this.pov_matrix[6] * v[2],
	     this.pov_matrix[1] * v[0] + this.pov_matrix[4] * v[1] + this.pov_matrix[7] * v[2],
	     this.pov_matrix[2] * v[0] + this.pov_matrix[5] * v[1] + this.pov_matrix[8] * v[2] ];
}


	// coords.cpp void spherical_to_cartesian(const Spherical& s, Cartesian& c) {
Sidereal.prototype.ra_dec_to_geo_cartesian = function(ra, dec, r) {
  var rcosphi  = r * cos_d(dec);

  return [ rcosphi * cos_h(ra - this.gast),
	   rcosphi * sin_h(ra - this.gast),
	   r       * sin_d(dec)
	 ];
}


function ra_dec_to_celestial_cartesian(ra, dec, r) {
  var rcosphi  = r * cos_d(dec);

  return [ rcosphi * cos_h(ra),
	   rcosphi * sin_h(ra),
	   r       * sin_d(dec)
	 ];
}

function ra_dec_to_geo_cartesian(gast, ra, dec, r) {
  var rcosphi  = r * cos_d(dec);

  return [ rcosphi * cos_h(ra - gast),
	   rcosphi * sin_h(ra - gast),
	   r       * sin_d(dec)
	 ];
}


Sidereal.prototype.geo_cartesian_to_az_alt = function(x, y, z) {
  var s = cartesian_to_spherical(this.pov_matrix[0] * x + this.pov_matrix[1] * y + this.pov_matrix[2] * z,
                                 this.pov_matrix[3] * x + this.pov_matrix[4] * y + this.pov_matrix[5] * z,
	                         this.pov_matrix[6] * x + this.pov_matrix[7] * y + this.pov_matrix[8] * z);

		// azimuth goes cw, theta ccw
  return { az:  -s.theta * this.DEGREES_PER_RADIAN,
           alt:  s.phi   * this.DEGREES_PER_RADIAN };
}


	// coords.cpp void cartesian_to_spherical(const Cartesian& c, Spherical& s) 
function cartesian_to_spherical(x, y, z) {
  var xxyy  = x * x + y * y;
  var rxy   = Math.sqrt(xxyy);
  var r     = Math.sqrt(xxyy + z * z);
  var theta = Math.atan2(y, x);
  var phi   = Math.atan2(z, rxy);

	// HACK: guarantee -180 never occurs (i.e. prefer 180)
	//  so that azimuth (== -theta) is always -180 rather than 180
  // if (s._theta <= -180.)
  //   s._theta = 180.;

  return { theta: theta,
	   phi:   phi,
	   r:     r };
}

function spherical_to_cartesian(r, theta, phi) {
  var sintheta = sin_d(theta);
  var costheta = cos_d(theta);
  var sinphi   = sin_d(phi);
  var cosphi   = cos_d(phi);
  var rcosphi  = r * cosphi;
  return [ rcosphi * costheta, rcosphi * sintheta, r * sinphi];
}