/* JavaScript functions for positional astronomy by John Walker -- September, MIM http://www.fourmilab.ch/ This program is in the public domain. */ // Frequently-used constants var J2000 = 2451545.0, /* Julian day of J2000 epoch */ JulianCentury = 36525.0, /* Days in Julian century */ AstronomicalUnit = 149597870.0; /* Astronomical unit in kilometres */ /* ASTOR -- Arc-seconds to radians. */ function astor(a) { return a * (Math.PI / (180.0 * 3600.0)); } /* DTR -- Degrees to radians. */ function dtr(d) { return (d * Math.PI) / 180.0; } /* RTD -- Radians to degrees. */ function rtd(r) { return (r * 180.0) / Math.PI; } /* FIXANGLE -- Range reduce angle in degrees. */ function fixangle(a) { return a - 360.0 * (Math.floor(a / 360.0)); } /* FIXANGR -- Range reduce angle in radians. */ function fixangr(a) { return a - (2 * Math.PI) * (Math.floor(a / (2 * Math.PI))); } /* UTCTOJ -- Convert GMT date and time to astronomical Julian time (i.e. Julian date plus day fraction, expressed as a double). Note that the mon (month) argument is in Unix index form: 0 = January, 11 = December. */ function utctoj(year, mon, mday, hour, min, sec) { /* Algorithm as given in Meeus, Astronomical Algorithms, Chapter 7, page 61 */ var a, b, m, y; m = mon + 1; y = year; if (m <= 2) { y--; m += 12; } /* Determine whether date is in Julian or Gregorian calendar based on canonical date of calendar reform. */ if ((year < 1582) || ((year == 1582) && ((mon < 9) || (mon == 9 && mday < 5)))) { b = 0; } else { a = Math.floor(y / 100); b = 2 - a + Math.floor(a / 4); } return ((Math.floor((365.25 * (y + 4716))) + Math.floor((30.6001 * (m + 1))) + mday + b) - 1524.5) + ((sec + 60 * (min + 60 * hour)) / 86400.0); } /* JYEAR -- Convert Julian date to year, month, day, which are returned as an Array. */ function jyear(td) { var z, f, a, alpha, b, c, d, e, mm; td += 0.5; z = Math.floor(td); f = td - z; if (z < 2299161.0) { a = z; } else { alpha = Math.floor((z - 1867216.25) / 36524.25); a = z + 1 + alpha - Math.floor(alpha / 4); } b = a + 1524; c = Math.floor((b - 122.1) / 365.25); d = Math.floor(365.25 * c); e = Math.floor((b - d) / 30.6001); mm = Math.floor((e < 14) ? (e - 1) : (e - 13)); return new Array( Math.floor((mm > 2) ? (c - 4716) : (c - 4715)), mm, Math.floor(b - d - Math.floor(30.6001 * e) + f) ); } /* JHMS -- Convert Julian time to hour, minutes, and seconds, returned as a three-element array. */ function jhms(j) { var ij; j += 0.5; /* Astronomical to civil */ ij = (j - Math.floor(j)) * 86400.0; return new Array( Math.floor(ij / 3600), Math.floor((ij / 60) % 60), Math.floor(ij % 60)); } /* GMST -- Calculate Greenwich Mean Sidereal Time for a given instant expressed as a Julian date and fraction. */ function gmst(jd) { var t, theta0; /* Time, in Julian centuries of 36525 ephemeris days, measured from the epoch 1900 January 0.5 ET. */ t = ((Math.floor(jd + 0.5) - 0.5) - 2415020.0) / JulianCentury; theta0 = 6.6460656 + 2400.051262 * t + 0.00002581 * t * t; t = (jd + 0.5) - (Math.floor(jd + 0.5)); theta0 += (t * 24.0) * 1.002737908; theta0 = (theta0 - 24.0 * (Math.floor(theta0 / 24.0))); return theta0; } /* OBLIQEQ -- Calculate the obliquity of the ecliptic for a given Julian date. This uses Laskar's tenth-degree polynomial fit (J. Laskar, Astronomy and Astrophysics, Vol. 157, page 68 [1986]) which is accurate to within 0.01 arc second between AD 1000 and AD 3000, and within a few seconds of arc for +/-10000 years around AD 2000. If we're outside the range in which this fit is valid (deep time) we simply return the J2000 value of the obliquity, which happens to be almost precisely the mean. */ var oterms = new Array ( -4680.93, -1.55, 1999.25, -51.38, -249.67, -39.05, 7.12, 27.87, 5.79, 2.45 ); function obliqeq(jd) { var eps = 23 + (26 / 60.0) + (21.448 / 3600.0), u, v, i; v = u = (jd - J2000) / (JulianCentury * 100); if (Math.abs(u) < 1.0) { for (i = 0; i < 10; i++) { eps += (oterms[i] / 3600.0) * v; v *= u; } } return eps; } /* Periodic terms for nutation in longiude (delta \Psi) and obliquity (delta \Epsilon) as given in table 21.A of Meeus, "Astronomical Algorithms", first edition. */ var nutArgMult = new Array( 0, 0, 0, 0, 1, -2, 0, 0, 2, 2, 0, 0, 0, 2, 2, 0, 0, 0, 0, 2, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, -2, 1, 0, 2, 2, 0, 0, 0, 2, 1, 0, 0, 1, 2, 2, -2, -1, 0, 2, 2, -2, 0, 1, 0, 0, -2, 0, 0, 2, 1, 0, 0, -1, 2, 2, 2, 0, 0, 0, 0, 0, 0, 1, 0, 1, 2, 0, -1, 2, 2, 0, 0, -1, 0, 1, 0, 0, 1, 2, 1, -2, 0, 2, 0, 0, 0, 0, -2, 2, 1, 2, 0, 0, 2, 2, 0, 0, 2, 2, 2, 0, 0, 2, 0, 0, -2, 0, 1, 2, 2, 0, 0, 0, 2, 0, -2, 0, 0, 2, 0, 0, 0, -1, 2, 1, 0, 2, 0, 0, 0, 2, 0, -1, 0, 1, -2, 2, 0, 2, 2, 0, 1, 0, 0, 1, -2, 0, 1, 0, 1, 0, -1, 0, 0, 1, 0, 0, 2, -2, 0, 2, 0, -1, 2, 1, 2, 0, 1, 2, 2, 0, 1, 0, 2, 2, -2, 1, 1, 0, 0, 0, -1, 0, 2, 2, 2, 0, 0, 2, 1, 2, 0, 1, 0, 0, -2, 0, 2, 2, 2, -2, 0, 1, 2, 1, 2, 0, -2, 0, 1, 2, 0, 0, 0, 1, 0, -1, 1, 0, 0, -2, -1, 0, 2, 1, -2, 0, 0, 0, 1, 0, 0, 2, 2, 1, -2, 0, 2, 0, 1, -2, 1, 0, 2, 1, 0, 0, 1, -2, 0, -1, 0, 1, 0, 0, -2, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 2, 0, -1, -1, 1, 0, 0, 0, 1, 1, 0, 0, 0, -1, 1, 2, 2, 2, -1, -1, 2, 2, 0, 0, -2, 2, 2, 0, 0, 3, 2, 2, 2, -1, 0, 2, 2 ); var nutArgCoeff = new Array( -171996, -1742, 92095, 89, /* 0, 0, 0, 0, 1 */ -13187, -16, 5736, -31, /* -2, 0, 0, 2, 2 */ -2274, -2, 977, -5, /* 0, 0, 0, 2, 2 */ 2062, 2, -895, 5, /* 0, 0, 0, 0, 2 */ 1426, -34, 54, -1, /* 0, 1, 0, 0, 0 */ 712, 1, -7, 0, /* 0, 0, 1, 0, 0 */ -517, 12, 224, -6, /* -2, 1, 0, 2, 2 */ -386, -4, 200, 0, /* 0, 0, 0, 2, 1 */ -301, 0, 129, -1, /* 0, 0, 1, 2, 2 */ 217, -5, -95, 3, /* -2, -1, 0, 2, 2 */ -158, 0, 0, 0, /* -2, 0, 1, 0, 0 */ 129, 1, -70, 0, /* -2, 0, 0, 2, 1 */ 123, 0, -53, 0, /* 0, 0, -1, 2, 2 */ 63, 0, 0, 0, /* 2, 0, 0, 0, 0 */ 63, 1, -33, 0, /* 0, 0, 1, 0, 1 */ -59, 0, 26, 0, /* 2, 0, -1, 2, 2 */ -58, -1, 32, 0, /* 0, 0, -1, 0, 1 */ -51, 0, 27, 0, /* 0, 0, 1, 2, 1 */ 48, 0, 0, 0, /* -2, 0, 2, 0, 0 */ 46, 0, -24, 0, /* 0, 0, -2, 2, 1 */ -38, 0, 16, 0, /* 2, 0, 0, 2, 2 */ -31, 0, 13, 0, /* 0, 0, 2, 2, 2 */ 29, 0, 0, 0, /* 0, 0, 2, 0, 0 */ 29, 0, -12, 0, /* -2, 0, 1, 2, 2 */ 26, 0, 0, 0, /* 0, 0, 0, 2, 0 */ -22, 0, 0, 0, /* -2, 0, 0, 2, 0 */ 21, 0, -10, 0, /* 0, 0, -1, 2, 1 */ 17, -1, 0, 0, /* 0, 2, 0, 0, 0 */ 16, 0, -8, 0, /* 2, 0, -1, 0, 1 */ -16, 1, 7, 0, /* -2, 2, 0, 2, 2 */ -15, 0, 9, 0, /* 0, 1, 0, 0, 1 */ -13, 0, 7, 0, /* -2, 0, 1, 0, 1 */ -12, 0, 6, 0, /* 0, -1, 0, 0, 1 */ 11, 0, 0, 0, /* 0, 0, 2, -2, 0 */ -10, 0, 5, 0, /* 2, 0, -1, 2, 1 */ -8, 0, 3, 0, /* 2, 0, 1, 2, 2 */ 7, 0, -3, 0, /* 0, 1, 0, 2, 2 */ -7, 0, 0, 0, /* -2, 1, 1, 0, 0 */ -7, 0, 3, 0, /* 0, -1, 0, 2, 2 */ -7, 0, 3, 0, /* 2, 0, 0, 2, 1 */ 6, 0, 0, 0, /* 2, 0, 1, 0, 0 */ 6, 0, -3, 0, /* -2, 0, 2, 2, 2 */ 6, 0, -3, 0, /* -2, 0, 1, 2, 1 */ -6, 0, 3, 0, /* 2, 0, -2, 0, 1 */ -6, 0, 3, 0, /* 2, 0, 0, 0, 1 */ 5, 0, 0, 0, /* 0, -1, 1, 0, 0 */ -5, 0, 3, 0, /* -2, -1, 0, 2, 1 */ -5, 0, 3, 0, /* -2, 0, 0, 0, 1 */ -5, 0, 3, 0, /* 0, 0, 2, 2, 1 */ 4, 0, 0, 0, /* -2, 0, 2, 0, 1 */ 4, 0, 0, 0, /* -2, 1, 0, 2, 1 */ 4, 0, 0, 0, /* 0, 0, 1, -2, 0 */ -4, 0, 0, 0, /* -1, 0, 1, 0, 0 */ -4, 0, 0, 0, /* -2, 1, 0, 0, 0 */ -4, 0, 0, 0, /* 1, 0, 0, 0, 0 */ 3, 0, 0, 0, /* 0, 0, 1, 2, 0 */ -3, 0, 0, 0, /* -1, -1, 1, 0, 0 */ -3, 0, 0, 0, /* 0, 1, 1, 0, 0 */ -3, 0, 0, 0, /* 0, -1, 1, 2, 2 */ -3, 0, 0, 0, /* 2, -1, -1, 2, 2 */ -3, 0, 0, 0, /* 0, 0, -2, 2, 2 */ -3, 0, 0, 0, /* 0, 0, 3, 2, 2 */ -3, 0, 0, 0 /* 2, -1, 0, 2, 2 */ ); /* NUTATION -- Calculate the nutation in longitude, deltaPsi, and obliquity, deltaEpsilon for a given Julian date jd. Results are returned as a two element Array giving (deltaPsi, deltaEpsilon). */ function nutation(jd) { var deltaPsi, deltaEpsilon, i, j, t = (jd - 2451545.0) / 36525.0, t2, t3, to10, ta = new Array, dp = 0, de = 0, ang; t3 = t * (t2 = t * t); /* Calculate angles. The correspondence between the elements of our array and the terms cited in Meeus are: ta[0] = D ta[0] = M ta[2] = M' ta[3] = F ta[4] = \Omega */ ta[0] = dtr(297.850363 + 445267.11148 * t - 0.0019142 * t2 + t3 / 189474.0); ta[1] = dtr(357.52772 + 35999.05034 * t - 0.0001603 * t2 - t3 / 300000.0); ta[2] = dtr(134.96298 + 477198.867398 * t + 0.0086972 * t2 + t3 / 56250.0); ta[3] = dtr(93.27191 + 483202.017538 * t - 0.0036825 * t2 + t3 / 327270); ta[4] = dtr(125.04452 - 1934.136261 * t + 0.0020708 * t2 + t3 / 450000.0); /* Range reduce the angles in case the sine and cosine functions don't do it as accurately or quickly. */ for (i = 0; i < 5; i++) { ta[i] = fixangr(ta[i]); } to10 = t / 10.0; for (i = 0; i < 63; i++) { ang = 0; for (j = 0; j < 5; j++) { if (nutArgMult[(i * 5) + j] != 0) { ang += nutArgMult[(i * 5) + j] * ta[j]; } } dp += (nutArgCoeff[(i * 4) + 0] + nutArgCoeff[(i * 4) + 1] * to10) * Math.sin(ang); de += (nutArgCoeff[(i * 4) + 2] + nutArgCoeff[(i * 4) + 3] * to10) * Math.cos(ang); } /* Return the result, converting from ten thousandths of arc seconds to radians in the process. */ deltaPsi = dtr(dp / (3600.0 * 10000.0)); deltaEpsilon = dtr(de / (3600.0 * 10000.0)); return new Array(deltaPsi, deltaEpsilon); } /* ECLIPTOEQ -- Convert celestial (ecliptical) longitude and latitude into right ascension (in degrees) and declination. We must supply the time of the conversion in order to compensate correctly for the varying obliquity of the ecliptic over time. The right ascension and declination are returned as a two-element Array in that order. */ function ecliptoeq(jd, Lambda, Beta) { var eps, Ra, Dec; /* Obliquity of the ecliptic. */ eps = dtr(obliqeq(jd)); log += "Obliquity: " + rtd(eps) + "\n"; Ra = rtd(Math.atan2((Math.cos(eps) * Math.sin(dtr(Lambda)) - (Math.tan(dtr(Beta)) * Math.sin(eps))), Math.cos(dtr(Lambda)))); log += "RA = " + Ra + "\n"; Ra = fixangle(rtd(Math.atan2((Math.cos(eps) * Math.sin(dtr(Lambda)) - (Math.tan(dtr(Beta)) * Math.sin(eps))), Math.cos(dtr(Lambda))))); Dec = rtd(Math.asin((Math.sin(eps) * Math.sin(dtr(Lambda)) * Math.cos(dtr(Beta))) + (Math.sin(dtr(Beta)) * Math.cos(eps)))); return new Array(Ra, Dec); } /* LOCALPOS -- Calculate local hour angle, altitude, and azimuth and return in an Array as follows: 0: Local hour angle (radians) 1: Altitude (degrees) 2: Azimuth (degrees) We return the local hour angle in radians because you're probably going to pass it to a trig function which expects radians anyway. */ function localpos(jd, ra, dec, siteLat, siteLong) { var igmst, latsin, latcos, lha, alt, az; igmst = gmst(jd); latsin = Math.sin(dtr(siteLat)); latcos = Math.cos(dtr(siteLat)); lha = dtr(fixangle((igmst * 15) - siteLong - ra)); az = rtd(Math.atan2(Math.sin(lha), Math.cos(lha) * latsin - Math.tan(dtr(dec)) * latcos)); alt = rtd(Math.asin(latsin * Math.sin(dtr(dec)) + latcos * Math.cos(dtr(dec)) * Math.cos(lha))); return new Array(lha, alt, az); }