// $Revision: 1.4 $
// Some constants
var TIME_UTC = 0;
var TIME_UT1 = 1;

var PI    = Math.PI;
var TWOPI = 2 * PI;
var DEGTORAD = TWOPI/360.0;
var RADTODEG = 360.0/TWOPI;
var MINSPERDAY = 1440.0;
var DAYSPERMIN = 1/MINSPERDAY;
var SECSPERDAY = 86400.0;
var DAYSPERSEC = 1/SECSPERDAY;
var MILLISECSPERDAY = 86400000.0;

var J1900 = 15019.5; // MJD of 12h January 0th 1900
var J2000 = 51544.5; // MJD of 12h January 1st 2000

// Some fixed offsets
var TTOFFSET = 32.184;   // TT = TAI + TTOFFSET
var GPSOFFSET = -19.0;   // GPS = TAI + GPSOFFSET
var LORANOFFSET = -10.0; // LORAN = TAI + LORANOFFSET;

var NTPINTERVAL = 4294967296.0; // 2^32

// TAI offsets from UTC after (and including) the given MJDs
// Now updated for the Dec 2005 leap second!
var taiOffsets = 
[ 41317, 10.0,
  41499, 11.0,
  41683, 12.0,
  42048, 13.0,
  42413, 14.0,
  42778, 15.0,
  43144, 16.0,
  43509, 17.0,
  43874, 18.0,
  44239, 19.0,
  44786, 20.0,
  45151, 21.0,
  45516, 22.0,
  46247, 23.0,
  47161, 24.0,
  47892, 25.0,
  48257, 26.0,
  48804, 27.0,
  49169, 28.0,
  49534, 29.0,
  50083, 30.0,
  50630, 31.0,
  51179, 32.0,
  53736, 33.0
];

//Utilities

// Round towards zero
function rtz(x) 
{
  if (x >= 0) return Math.floor(x);
  else return Math.ceil(x);
}

function fmod(x,y)
{
  assert(x && x < Infinity && x > -Infinity,"Bad fmod parameter: " + x);
  while (x >= y) x -= y;
  while (x < 0) x += y;
  return x;
}

function putinrange(x,lower,upper,inc)
{
  assert(x && x < Infinity && x > -Infinity,"Bad putinrange parameter: " + x);
  while (x < lower) x += inc;
  while (x >= upper) x -= inc;
  return x;
}

// Check the browser capabilities
var noregexp = false;
var noutc = false;
var notimezone = false;
var nomilliseconds = false;

function check_capabilities()
{
  var sex = "123:45:67.89";
  var dec = "3.14159";
  var a;
  // Check the regular expression support is up to it.
  // njs-js has s.match == undefined, but s.match() is OK
  a = sex.match(sexregexp);
  if (!a || a.length != 7) {
    noregexp = true;
  }
  a = dec.match(decregexp);
  if (!a || a.length != 3) {
    noregexp = true;
  }
  var d = new Date();
  if (!d.getUTCHours || !d.getUTCHours()) {
    noutc = true;
  }
  if (!d.getTimezoneOffset || !d.getTimezoneOffset()) {
    notimezone = true;
  }
  if (!d.getMilliseconds || !d.getMilliseconds()) {
    nomilliseconds = true;
  }
}
      
// Astro stuff starts here
function getTaiOffset(mjd) {
  var i = taiOffsets.length-2;
  while (i > 0) {
    if (mjd >= taiOffsets[i])
      return taiOffsets[i+1];
    i -= 2;
  }
  return null; // Return null for invalid MJDs
}

// Equation of the equinoxes
function get_eqeq(mjd)
{
  var t = (mjd - J2000)/36525.0;
  var Om = DEGTORAD * (125.04452 - 1934.136261 * t);
  var L =  DEGTORAD * (280.4665 + 36000.7698 * t);
  var L1 = DEGTORAD * (218.3165 + 481267.8813 * t);
  var eps = 23.439281 - 0.013 * t;
  var dp = -17.2*Math.sin(Om) -1.32*Math.sin(2*L) -0.23*Math.sin(2*L1) +0.21*Math.sin(2*Om);
  var de =   9.2*Math.cos(Om) +0.57*Math.cos(2*L) + 0.1*Math.cos(2*L1) -0.09*Math.cos(2*Om);
  var e =  DEGTORAD * (eps + de/3600.0);
  return dp * Math.cos(e) / 15; // Convert arcsec to sec
}

// Sidereal time conversion
function get_gmst(mjd, ut1)
{
  var t;
  var gmst;
  //  Julian centuries since J2000.
  t = (ut1 + (mjd - J2000)) / 36525.0;

  // gmst at this ut1.
  gmst = DAYSPERSEC * (24110.54841 + (8640184.812866 + (0.093104 - 6.2E-6 * t) * t) *t
		       + 86400.0 * (fmod (ut1, 1.0)+fmod(mjd, 1.0)));
  gmst = fmod(gmst,1.0);
  return gmst;
}

// Calendar to Modified Julian date
function get_mjd(y,m,d)
{
  mjd = 
    rtz((1461 * (y + 4800 + rtz((m - 14) / 12 ))) / 4) +
    rtz((367 * (m - 2 - 12 * rtz((m - 14) / 12 ))) / 12) -
    rtz((3 * rtz((y + 4900 + rtz(( m - 14) / 12)) / 100)) / 4) +
    d - 2432076;
  return mjd;
}

function earth(t, pv)
{
  //  Mean Earth:EMB distance and speed, AU and AU/s
  var remb = 3.12e-5;

  var elm,gamma,em,elt,eps0,
    e,esq,v,r,elmm,coselt,sineps,coseps,w1,selmm,celmm;

  //  Geometric mean longitude of Sun
  //  (cf 4.881627938+6.283319509911*T MOD 2PI)
  elm = fmod(4.881628 + TWOPI * (t-Math.floor(t)) + 0.00013420 * t, TWOPI);

  //  Mean longitude of perihelion
  gamma = 4.908230 + 3.0005e-4 * t;

  //  Mean anomaly
  em = elm - gamma;

  //  Mean obliquity
  eps0 = 0.40931975 - 2.27e-6 * t;

  //  Eccentricity
  e = 0.016751 - 4.2e-7 * t;
  esq = e * e;

  //  True anomaly
  v = em + 2.0 * e * Math.sin(em) + 1.25 * esq * Math.sin(2.0 * em);

  //  True ecliptic longitude
  elt = v + gamma;

  //  True distance
  r = (1.0 - esq) / (1.0 + e * Math.cos(v));

  //  Moon's mean longitude
  elmm = fmod(4.72 + 83.9971 * t, TWOPI);

  //  Useful functions
  coselt=Math.cos(elt);
  sineps=Math.sin(eps0);
  coseps=Math.cos(eps0);
  w1=-r * Math.sin(elt);
  selmm=Math.sin(elmm);
  celmm=Math.cos(elmm);

  //  Earth position and velocity
  pv[0] = -r * coselt - remb * celmm;
  pv[1] = (w1 - remb * selmm) * coseps;
  pv[2] = w1 * sineps;
}

function earth_mjd(mjd, pv)
{
  earth((mjd - J1900)/365.25, pv);
}

function get_ra(pv)
{
    return Math.atan2(-pv[1], -pv[0]);
}

function get_dec(pv)
{
    return Math.atan2(-pv[2],Math.sqrt(pv[0]*pv[0]+pv[1]*pv[1]));
}

function get_earthdistance(pv)
{
    return Math.sqrt(pv[0]*pv[0] + pv[1]*pv[1] + pv[2]*pv[2]);
}

var expArray = [1,
		10,
		100,
		1000,
		10000,
		100000,
		1000000,
		10000000,
		100000000,
		1000000000,
		10000000000,
		100000000000,
		1000000000000
];

var padArray = ["",
		"0",
		"00",
		"000",
		"0000",
		"00000",
		"000000",
		"0000000",
		"00000000",
		"000000000",
		"0000000000",
		"00000000000",
		"000000000000"
];

function float_string(x,iprec,prec)
{
  var a = new Array();
  float_string_aux(x,iprec,prec,a);
  reverse(a);
  return a.join("");
}

function float_string_signed(x,iprec,prec)
{
  var neg;
  var a = new Array();
  if (x < 0) {
    neg = true;
    x = -x;
  } else {
    neg = false;
  }
  float_string_aux(x,iprec,prec,a);
  if (neg) {
    a.push("-");
  } else {
    a.push(" ");
  }
  reverse(a);
  return a.join("");
}

function float_string_aux(r,iprec,prec,a)
{
  var precexp10 = expArray[prec];
  // Round to the appropriate precision and convert to an integer
  r = Math.floor(r * precexp10 + 0.5);
  t = (r%precexp10).toString(); // Get the fraction part
  a.push(t);
  r = Math.floor(r/precexp10);
  if (t.length < prec) {
    a.push(padArray[prec-t.length]);
  }
  a.push(".");

  t = r.toString();
  a.push(t);
  if (t.length < iprec) {
    a.push(padArray[iprec-t.length]);
  }
  return a;
}

function sexigesimal_string(r,iprec,prec)
{
  var a = new Array();
  sexigesimal_string_aux(r,iprec,prec,a);
  reverse(a);
  return a.join("");
}
function sexigesimal_string_signed(r,iprec,prec)
{
  var a = new Array();
  var neg;
  if (r < 0) {
    r = -r;
    neg = true;
  } else {
    neg = false;
  }
  sexigesimal_string_aux(r,iprec,prec,a);
  if (neg) {
    a.push("-");
  } else {
    a.push(" ");
  }
  reverse(a);
  return a.join("");
}

function sexigesimal_string_aux(r,iprec,prec,a)
{
  var precexp10 = expArray[prec];
  // Round to the appropriate precision and convert to an integer
  r = Math.floor(r * 3600.0 * precexp10 + 0.5);
  t = (r%precexp10).toString(); // Get the fraction part
  a.push(t);
  r = Math.floor(r/precexp10);
  if (t.length < prec) {
    a.push(padArray[prec-t.length]);
  }
  a.push(".");

  t = (r%60).toString();
  a.push(t);
  r = Math.floor(r/60);
  if (t.length < 2) {
    a.push(padArray[2-t.length]);
  }
  a.push(":");

  t = (r%60).toString();
  a.push(t);
  r = Math.floor(r/60);
  if (t.length < 2) {
    a.push(padArray[2-t.length]);
  }
  a.push(":");

  t = r.toString();
  a.push(t);
  if (t.length < iprec) {
    a.push(padArray[iprec-t.length]);
  }
  return a;
}

function sexigesimal_string2(r,iprec,prec)
{
  var a = new Array();
  sexigesimal_string2_aux(r,iprec,prec,a);
  reverse(a);
  return a.join("");
}
function sexigesimal_string2_signed(r,iprec,prec)
{
  var a = new Array();
  var neg;
  if (r < 0) {
    r = -r;
    neg = true;
  } else {
    neg = false;
  }
  sexigesimal_string2_aux(r,iprec,prec,a);
  if (neg) {
    a.push("-");
  } else {
    a.push(" ");
  }
  reverse(a);
  return a.join("");
}

function sexigesimal_string2_aux(r,iprec,prec,a)
{
  var precexp10 = expArray[prec];
  // Round to the appropriate precision and convert to an integer
  r = Math.floor(r * 60.0 * precexp10 + 0.5);
  t = (r%precexp10).toString(); // Get the fraction part
  a.push(t);
  r = Math.floor(r/precexp10);
  if (t.length < prec) {
    a.push(padArray[prec-t.length]);
  }
  a.push(".");

  t = (r%60).toString();
  a.push(t);
  r = Math.floor(r/60);
  if (t.length < 2) {
    a.push(padArray[2-t.length]);
  }
  a.push(":");

  t = r.toString();
  a.push(t);
  if (t.length < iprec) {
    a.push(padArray[iprec-t.length]);
  }
  return a;
}

function day_to_hms_string(t)
{
  t = fmod(t,1.0);
  return sexigesimal_string(24*t,2,3);
}

function day_to_hm_string(t)
{
  t = fmod(t,1.0);
  return sexigesimal_string2(24*t,2,3);
}

function day_to_dms_string(t)
{
  t = fmod(t,1.0);
  return rad_to_dms_string(TWOPI * t);
}


function day_to_dm_string(t)
{
  t = fmod(t,1.0);
  return rad_to_dm_string(TWOPI * t);
}

function rad_to_hms_string(r)
{
  return day_to_hms_string(r/TWOPI);
}

function rad_to_dms_string(r)
{
  r /= TWOPI;
  r = putinrange(r,-1.0,1.0,1.0);
  return sexigesimal_string_signed(360*r,3,3);
}

function rad_to_dm_string(r)
{
  r /= TWOPI;
  r = putinrange(r,-1.0,1.0,1.0);
  return sexigesimal_string2_signed(360*r,3,3);
}


var hexDigits = ["0","1","2","3","4","5","6","7","8","9","A","B","C","D","E","F"];

function hex_string_aux(n,digits,a)
{
  assert (n == Math.floor(n), "Hex string internal 1");
  assert (n >= 0, "Hex string internal error 2");
  while (digits > 0) {
    var i1 = Math.floor((n % 16));
    /*
    // This cack a workaround for njs-js bug
    var i = 0;
    while (i < i1) i++;
    */
    var c = hexDigits[i1]; 
    a.push(c);
    digits--;
    n /= 16.0;
  }
}

function reverse(a) 
{
  // Reverse a
  var i = 0;
  var j = a.length-1;
  while (i < j) {
    var t = a[i];
    a[i] = a[j];
    a[j] = t;
    i++; j--;
  }
}

function hex_string(n)
{
  var a = new Array();
  hex_string_aux(n,8,a);
  reverse(a);
  return a.join("");
}
  
function make_ntp_string(mjd,utc)
{
  // Number of ntp seconds
  var s1 = 86400.0 * (mjd - J1900 - 0.5); // Seconds for whole days
  var s2 = 86400.0 * utc; // Seconds in this day
  var t = Math.floor(s2);
  s1 += t; // Whole seconds
  s2 -= t; // fraction of a second
  // Reduce s1 to appropriate interval
  s1 = fmod(s1,NTPINTERVAL);
  // diag (s1.toString(16));
  // Now build the string, backwards
  var a = new Array();
  hex_string_aux(Math.floor(NTPINTERVAL * s2 + 0.5), 8, a);
  a.push(".");
  hex_string_aux(s1, 8, a);
  // Reverse a
  var i = 0;
  var j = a.length-1;
  while (i < j) {
    var t = a[i];
    a[i] = a[j];
    a[j] = t;
    i++; j--;
  }
  return a.join("");
}
    

var sexregexp = /([\+\-]?)([0-9]+):([0-9]+)(:([0-9]+(\.[0-9]+)?))?/;
var decregexp = /([\+\-]?[0-9]+(\.[0-9]+)?)/;

function decode(s)
{
  var a;
  if (noregexp) return null; // Can't do regexp matching here
    //Decimal
  a = s.match(decregexp);
  if (a && a[0] == s) {
    return Number(a[1]);
  } else { 
    // Sexigesimal
    a = s.match(sexregexp);
    if (a && a[0] == s) {
      var n = 0;
      if (a.length > 4) {
	n += Number(a[5]);
	n /= 60;
      }
      if (a.length > 2) {
	n += Number(a[3]);
	n /= 60;
      }
      n += Number(a[2]);
      if (a[1] == "-") {
	n = -n;
      }
      return n;
    }
  }
}

function initTimes(times,date)
{
  var mjd,y,m,d,h,min,s,ms,tz;
  if (noutc) {
    y = date.getYear();
    m = date.getMonth() + 1;
    d = date.getDate();
    h = date.getHours();
    min = date.getMinutes();
    s = date.getSeconds();
    if (y < 1900) y += 1900; // Hack for njs-js
  } else {
    y = date.getUTCFullYear();
    m = date.getUTCMonth() + 1;
    d = date.getUTCDate();
    h = date.getUTCHours();
    min = date.getUTCMinutes();
    s = date.getUTCSeconds();
  }
  if (nomilliseconds) {
    ms = 0;
  } else {
    ms = date.getMilliseconds();
  }
  if (notimezone) {
    tz = 0;
  } else {
    tz = -date.getTimezoneOffset()/MINSPERDAY;
  }
  mjd = get_mjd(y,m,d);
  time = (1000 * (60 * (60 * h + min) + s) + ms) / MILLISECSPERDAY;
  assert (mjd != null && mjd >= 0 && mjd < 100000,"Invalid mjd ");
  assert (time && time >= 0 && time < 1,"Invalid time ");
  times.mjd = mjd;
  times.time = time;
  times.tz = tz;
  times.date = date;
}

function computeTimes(times,date,lat,lon,timetype)
{
  var mjd,time,mjdt;
  var utc,ut0,ut1,ut2,lt,st,lst;
  var tai,tt,tdb,gps,loran;
  var gmst,gast,lmst,last;
  var sun_ra,sun_dec,sun_gha,esd,eqeq,eqt;
  var t,dut1ut2,dut1ut0,x,y,dut1,longtime,taioffset;
  var pv = new Array(3);
  initTimes(times,date); // Set mjd, time & timezone
  mjd = times.mjd;
  time = times.time;
  mjdt = mjd + time;

  t = 2000.0 + (mjdt - 51544.03) / 365.2422; // Besselian years
  // dut1ut2 = UT2 - UT1
  dut1ut2 = 0.022 * Math.sin(2*PI*t) - 0.012 * Math.cos(2*PI*t) - 
            0.006 * Math.sin(4*PI*t) + 0.007 * Math.cos(4*PI*t);

  // IERS Bulletin A Prediction formula for DUT1 = UT1 - UT
  // Only valid for the era of the prediction.
  {
    var pi = PI;
    var MJD = mjdt;
    var A,C;
    var cos = Math.cos;
    var sin = Math.sin;
    // Really ought have a perl script or whatever to update this automatically
    // Bulletin A: 10 August 2006 Vol. XIX No. 032    

    A = 2*pi*(MJD-53957)/365.25;
    C = 2*pi*(MJD-53957)/435;
    x = 0.0536 +  0.1652 * cos(A) +  0.0681 * sin(A) -  0.1062 * cos(C) -  0.1663 * sin(C);
    y = 0.3538 +  0.0676 * cos(A) -  0.1457 * sin(A) -  0.1663 * cos(C) +  0.1062 * sin(C);
    // dut1 = UT1-UTC
    dut1 =   0.1299 - 0.00058 * (MJD - 53962) - dut1ut2;
  }

  var tanlat = Math.tan(lat);
  if (tanlat == null || tanlat > 1000 || tanlat < -1000) {
    warn("Latitude out of range, setting UT0-UT1 = 0");
    dut1ut0 = 0;
  } else {
    dut1ut0 = (tanlat * (x * Math.sin(lon) + y * Math.cos(lon)) / 15); // asec to sec
  }
  taioffset = getTaiOffset(mjd);
  if (taioffset == null) {
    warn ("MJD too early for TAI, setting TAI = UTC");
    taioffset = 0;
  }
  switch (timetype) {
  case TIME_UTC:
    utc = time;
    ut1 = utc + DAYSPERSEC * dut1;
    break;
  case TIME_UT1:
    ut1 = time;
    utc = ut1 - DAYSPERSEC * dut1;
    break;
  }
  longtime = lon/TWOPI; // Longitude time offset
  lt = utc + times.tz;
  ut0 = ut1 + DAYSPERSEC * dut1ut0;
  ut2 = ut1 + DAYSPERSEC * dut1ut2;
  tai = utc + DAYSPERSEC * taioffset;
  tt  = tai + DAYSPERSEC * TTOFFSET;
  gps = tai + DAYSPERSEC * GPSOFFSET;
  loran = tai + DAYSPERSEC * LORANOFFSET;
  g   = 6.2401 + 0.01720197 * (mjdt - J2000);
  tdb = tt + DAYSPERSEC * (0.001658 * Math.sin(g) + 0.000014 * Math.sin(2 * g));
  earth_mjd(mjd+tdb,pv); // Note: date is TDB
  sun_ra  = get_ra(pv); // Radians
  sun_dec = get_dec(pv); // Radians
  esd  = get_earthdistance(pv);
  eqeq = get_eqeq(mjdt);
  gmst = get_gmst(mjd,ut1);
  gast = gmst + DAYSPERSEC * eqeq;
  lmst = gmst + longtime;
  last = gast + longtime;
  sun_gha = gast - sun_ra/TWOPI;
  st = sun_gha + 0.5;
  lst = st + longtime;
  // Equation of time
  eqt = ut1 - st;
  eqt = putinrange(eqt,-0.5, +0.5, 1.0);
  eqt *= 24 * 60; // In minutes

  times.utc = utc;
  times.ut0 = ut0;
  times.ut1 = ut1;
  times.ut2 = ut2;
  times.tai = tai;
  times.gps = gps;
  times.loran = loran;
  times.tt = tt;
  times.tdb = tdb;
  times.gmst = gmst;
  times.gast = gast;
  times.lmst = lmst;
  times.last = last;
  times.lt = lt;
  times.st = st;
  times.lst = lst;
  times.sun_ra = sun_ra;
  times.sun_dec = sun_dec;
  times.sun_gha = sun_gha;
  times.esd = esd;
  times.eqeq = eqeq;
  times.eqt = eqt;
  times.x = x;
  times.y = y;
}

// Use the global times object
function mjdString(times) { return float_string(times.mjd+times.utc,0,5); }
function ltString(times)  { return day_to_hms_string(times.lt); }
function utcString(times) { return day_to_hms_string(times.utc); }
function ut0String(times) { return day_to_hms_string(times.ut0); }
function ut1String(times) { return day_to_hms_string(times.ut1); }
function ut2String(times) { return day_to_hms_string(times.ut2); }
function taiString(times) { return day_to_hms_string(times.tai); }
function gpsString(times) { return day_to_hms_string(times.gps); }
function loranString(times) { return day_to_hms_string(times.loran); }
function ttString(times)  { return day_to_hms_string(times.tt); }
function tdbString(times) { return day_to_hms_string(times.tdb); }
function gmstString(times){ return day_to_hms_string(times.gmst); }
function gastString(times){ return day_to_hms_string(times.gast); }
function lmstString(times){ return day_to_hms_string(times.lmst); }
function lastString(times){ return day_to_hms_string(times.last); }
function stString(times)  { return day_to_hms_string(times.st); }
function lstString(times) { return day_to_hms_string(times.lst); }
function ghaAriesString(times) { return day_to_dm_string(times.gast); }
function ghaSunString(times)   { return day_to_dm_string(times.sun_gha); }
function raSunString(times)    { return rad_to_hms_string(times.sun_ra); }
function decSunString(times)   { return rad_to_dms_string(times.sun_dec); }
function esdString(times)      { return float_string(times.esd,0,5); }
function eqeqString(times)     { return float_string_signed(times.eqeq,0,3); }
function eqtString(times)      { return signed_sexigesimal_string2(times.eqt,2,0); }
function ntpString(times)      { return make_ntp_string(times.mjd,times.utc); }

function writeTimes(times)
{
  var date = times.date;
  document.writeln (date.toString());
  document.writeln ("MJD: " + mjdString(times));
  document.writeln ("LT:  " + ltString(times));
  document.writeln ("UTC: " + utcString(times));
  document.writeln ("UT0: " + ut0String(times));
  document.writeln ("UT1: " + ut1String(times));
  document.writeln ("UT2: " + ut2String(times));
  document.writeln ("TAI: " + taiString(times));
  document.writeln ("GPS: " + gpsString(times));
  document.writeln ("LORAN: " + loranString(times));
  document.writeln ("TT:  " + ttString(times));
  document.writeln ("TDB: " + tdbString(times));
  document.writeln ("GMST:" + gmstString(times));
  document.writeln ("GAST:" + gastString(times));
  document.writeln ("LMST:" + lmstString(times));
  document.writeln ("LAST:" + lastString(times));
  document.writeln ("ST:" + stString(times));
  document.writeln ("LST:" + lstString(times));
  document.writeln ("NTP:" + ntpString(times));
  document.writeln ("GHA Aries:" + ghaAriesString(times));
  document.writeln ("GHA Sun:  " + ghaSunString(times));
  document.writeln ("Sun RA: " + raSunString(times));
  document.writeln ("Sun Dec:" + decSunString(times));
  document.writeln ("Earth-Sun distance: " + esdString(times));
  document.writeln ("Equation of the equinoxes: " + eqeqString(times));
  document.writeln ("Equation of time: " + eqtString(times));
}