Moon.js

//
// This script will NOT work without errors in Netscape 2. This
// is not only to do with the array.length property in Moon().
// Developed using Navigator 3.04, tested on Navigator 4.6 and
// MSIE 5. I do want a Javascript 1.0 version.
//
// This main function gets the UT (strictly TDT) date and time
// from form, does some basic input checking, calls the Moon
// function, the nutation function and the ecliptic function, and
// writes the resulting outputs to the form in the page.
//
// Speed: 100000 multiples in 25 seconds on a P75 - 67 seconds
// 46 seconds for Math.sin in radians.
//
 
function moondrawf(y8, m8, d8, h8, min8) {
var days;
days=day2000(y8, m8, d8, h8 + min8/60);
var EcMoon;
var EqMoon;
var Nutat;
EcMoonNut = new Array(0,0,0,0);
 
EcMoon = Moon(days);
nutat = nutation(days);
EcMoonNut[1] = EcMoon[1] + nutat[1];
EcMoonNut[2] = EcMoon[2] + nutat[2];
EcMoonNut[3] = EcMoon[3] + nutat[3];
EqMoon = equatorial(days, EcMoonNut);
EqMoon [4]=EcMoon[4]
 
return EqMoon;
}
 
 
 
 
 
 
 
 
function doCalcs(form) {
var days;
 
//
// Get date and time code from user, isolate the year, month, day and hours
// and minutes, and do some basic error checking! This only works for AD years
//
var g, y, m, d, bit, h, min, bk;
g = form.Date.value;
y = Math.floor(g / 10000);
m = Math.floor( (g - y*10000) / 100);
d = Math.floor(g - y*10000 - m*100 );
bit = (g - Math.floor(g))*100;
h = Math.floor(bit);
min = Math.floor(bit*100 - h * 100 + 0.5);
//
// primative error checking - accounting for right number of
// days per month including leap years. Using bk variable to
// prevent multiple alerts. See functions isleap(y)
// and goodmonthday(y, m, d).
//
bk = 0;
if(g < 14000000) {	
bk = 1;
alert("Routines are not accurate enough to work back that" +
" far\nErrors larger than 1 arcmin");
}
if(g > 25000000) {
bk = 1;	
alert("Routines are not accurate enough to work that far into the future\n" +
"Errors larger than 1 arcmin");
}	
if( ((m<1) || (m>12)) && (bk !=1)) {
bk = 1;
alert("Months are not right - type date again");
}
if((goodmonthday(y, m, d) ==0) && (bk !=1)) {
bk = 1;
alert("Wrong number of days for the month or not a leap year - type date again");
}
if ((h>23) && (bk != 1)) {
bk = 1;
alert("Hours are not right - type date again");
}
if((min>59) && (bk !=1)) {
alert("Minutes are not right - type date again");
}
//
// Get the number of days since J2000.0 using day2000() function
//
days = day2000(y, m, d, h + min/60);
 
//
// Proceed to calculate by calling functions
//
 
EcMoon = new Array(0,0,0,0);
EqMoon = new Array(0,0,0,0);
Nutat = new Array(0,0,0,0);
EcMoonNut = new Array(0,0,0,0);
 
EcMoon = Moon(days);
nutat = nutation(days);
EcMoonNut[1] = EcMoon[1] + nutat[1];
EcMoonNut[2] = EcMoon[2] + nutat[2];
EcMoonNut[3] = EcMoon[3] + nutat[3];
EqMoon = equatorial(days, EcMoonNut);
 
//
// write to the textboxes in the page
//
 
form.MoonLong.value = round(EcMoonNut[1] * RAD2DEG, 3);
form.MoonLat.value = round(EcMoonNut[2] * RAD2DEG, 3);
form.MoonDist.value = round(EcMoonNut[3], 0);
form.MoonRA.value = round(EqMoon[1]/15 *RAD2DEG, 4);
form.MoonDec.value = round(EqMoon[2] * RAD2DEG, 3);
}
 
//
// Define some global constants
//
 
var DEG2RAD = Math.PI / 180;
var RAD2DEG = 180 / Math.PI;
var TwoPi = Math.PI * 2;
 
//
// Moon returns mean ecliptic long and lat of Moon to 10 arc sec (long)
// and 3 arc sec (lat) for a few centuries either side of J2000.0.
// This function is taken from Arkkana Peck's Javascript code at
// http://www.shallowsky.com/moon/hitchhiker.html
// but with the addition of the Moon distance series
//
 
function Moon(day)	{
// Time measured in Julian centuries from epoch J2000.0:
var T = day / 36525;
var T2 = T*T;
var T3 = T2*T;
var T4 = T3*T;
ecliptic = new Array(0, 0, 0, 0, 0);
// Mean elongation of the moon:
var D = angle( 297.8502042
+ 445267.1115168 * T
- 0.0016300 * T2
+ T3 / 545868
+ T4 / 113065000 );
// Sun's mean anomaly:
var M = angle
( 357.5291092
+ 35999.0502909 * T
- 0.0001536 * T2
+ T3 / 24490000 );
// Moon's mean anomaly:
var Mprime = angle
( 134.9634114
+ 477198.8676313 * T
+ 0.0089970 * T2
- T3 / 3536000
+ T4 / 14712000 );
 
 
// Done getting the phase angle; now calculate the mean
// longitude and latitude.
 
// Moon's mean longitude:
var Lprime = angle(218.3164591 + 481267.88134236 * T
- .0013268 * T2 + T3 / 538841 - T4 / 65194000);
 
// Moon's argument of latitude (mean distance from ascending node):
var F = angle(93.2720993 + 483202.0175273 * T
- .0034029 * T2 - T3 / 3526000 + T4 / 863310000);
 
// Now, the fearsome neverending tables! (KPB - These are available
// online in machine format!)
 
DcA = new Array(0, 2, 2, 0, 0, 0, 2, 2, 2, 2, 0, 1, 0, 2, 0, 0,
4, 0, 4, 2, 2, 1, 1, 2, 2, 4, 2, 0, 2, 2, 1, 2,
0, 0, 2, 2, 2, 4, 0, 3, 2, 4, 0, 2, 2, 2, 4, 0,
4, 1, 2, 0, 1, 3, 4, 2, 0, 1, 2, 2);
 
McA = new Array(0, 0, 0, 0, 1, 0, 0, -1, 0, -1, 1, 0, 1, 0, 0, 0,
0, 0, 0, 1, 1, 0, 1, -1, 0, 0, 0, 1, 0, -1, 0, -2,
1, 2, -2, 0, 0, -1, 0, 0, 1, -1, 2, 2, 1, -1, 0,
0, -1, 0, 1, 0, 1, 0, 0, -1, 2, 1, 0, 0);
 
MpcA = new Array(1, -1, 0, 2, 0, 0, -2, -1, 1, 0, -1, 0, 1, 0, 1,
1, -1, 3, -2, -1, 0, -1, 0, 1, 2, 0, -3, -2,
-1, -2, 1, 0,
2, 0, -1, 1, 0, -1, 2, -1, 1, -2, -1, -1, -2, 0,
1, 4, 0, -2, 0, 2, 1, -2, -3, 2, 1, -1, 3, -1);
 
FcA = new Array(0, 0, 0, 0, 0, 2, 0, 0, 0, 0, 0, 0, 0, -2, 2, -2,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0, 0, 0,
0, 0, 0, -2, 2, 0, 2, 0, 0, 0, 0, 0, 0, -2, 0, 0,
0, 0, -2, -2, 0, 0, 0, 0, 0, 0, 0, -2);
 
 
// Longitude coefficients from Meeus' table A
 
ScA = new Array(6288774, 1274027, 658314, 213618, -185116, -114332,
58793, 57066, 53322, 45758, -40923, -34720, -30383,
15327, -12528, 10980, 10675, 10034, 8548, -7888,
-6766, -5163, 4987, 4036, 3994, 3861, 3665, -2689,
-2602, 2390, -2348, 2236,
-2120, -2069, 2048, -1773, -1595, 1215, -1110, -892,
-810, 759, -713, -700, 691, 596, 549, 537, 520,
-487, -399, -381, 351, -340, 330, 327, -323, 299, 294);
 
// Radius vector coefficients from Meeus' Table A
 
ScAR = new Array(-20905355, -3699111, -2955968, -569925,
48888, -3149, 246158, -152138, -170733, -204586, -129620,
108743, 104755, 10321, 0, 79661, -34782, -23210, -21636,
24208, 30824, -8379, -16675, -12831, -10445, -11650, 14403,
-7003, 0, 10056, 6322, -9884, 5751, 0, -4950, 4130, 0,
-3958, 0, 3258, 2616, -1897, -2117, 2354, 0, 0, -1423,
-1117, -1571, -1739, 0, -4421, 0, 0, 0, 0, 1165, 0, 0,
8752);
 
DcB = new Array(0, 0, 0, 2, 2, 2, 2, 0, 2, 0, 2, 2, 2, 2, 2, 2, 2, 0,
4, 0, 0, 0, 1, 0, 0, 0, 1, 0, 4, 4,
0, 4, 2, 2, 2, 2, 0, 2, 2, 2, 2, 4, 2, 2, 0, 2, 1, 1,
0, 2, 1, 2, 0, 4, 4, 1, 4, 1, 4, 2);
 
McB = new Array(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 1, -1, -1, -1,
1, 0, 1, 0, 1, 0, 1, 1, 1, 0, 0, 0, 0,
0, 0, 0, 0, -1, 0, 0, 0, 0, 1, 1, 0, -1, -2, 0, 1, 1,
1, 1, 1, 0, -1, 1, 0, -1, 0, 0, 0, -1, -2);
 
MpcB = new Array(0, 1, 1, 0, -1, -1, 0, 2, 1, 2, 0, -2, 1, 0, -1, 0,
-1, -1, -1, 0, 0, -1, 0, 1, 1, 0, 0, 3, 0, -1,
1, -2, 0, 2, 1, -2, 3, 2, -3, -1, 0, 0, 1, 0, 1, 1, 0,
0, -2, -1, 1, -2, 2, -2, -1, 1, 1, -1, 0, 0);
 
FcB = new Array(1, 1, -1, -1, 1, -1, 1, 1, -1, -1, -1, -1, 1, -1, 1,
1, -1, -1, -1, 1, 3, 1, 1, 1, -1, -1, -1, 1, -1, 1,
-3, 1, -3, -1, -1, 1, -1, 1, -1, 1, 1, 1, 1, -1, 3, -1,
-1, 1, -1, -1, 1, -1, 1, -1, -1, -1, -1, -1, -1, 1);
 
ScB = new Array(5128122, 280602, 277693, 173237, 55413, 46271, 32573,
17198, 9266, 8822, 8216, 4324, 4200, -3359, 2463, 2211,
2065, -1870, 1828, -1794, -1749, -1565, -1492, -1475,
-1410, -1344, -1335, 1107, 1021, 833,
777, 671, 607, 596, 491, -451, 439, 422, 421, -366,
-351, 331, 315, 302, -283, -229, 223, 223, -220, -220,
-185, 181, -177, 176, 166, -164, 132, -119, 115, 107);
 
 
// term involving the decreasing eccentricity of the earth's orbit:
var E = 1 - .002516 * T - .0000074 * T2;
var E2 = E*E;
 
var sigma_l = 0.;
var sigma_b = 0.;
var sigma_r = 0.;
// Sum over table A coeffs - this bit doesn't like Javascript 1.0
//
 
for (var i=0; i < ScA.length; ++i)
{
var ME = 1;
if (MpcA[i] == 1 || MpcA[i] == -1)
ME = E;
else if (MpcA[i] == 2 || MpcA[i] == -2)
ME = E2;
sigma_l += ScA[i] * Math.sin(DcA[i] * D + McA[i] * ME * M
+ MpcA[i] * Mprime + FcA[i] * F);
sigma_r += ScAR[i] *
Math.cos(DcA[i] * D + McA[i] * ME * M + MpcA[i] * Mprime
+ FcA[i] * F);
 
}
 
for (var i=0; i < ScB.length ; ++i)
{
var ME = 1;
if (MpcB[i] == 1 || MpcB[i] == -1)
ME = E;
else if (MpcB[i] == 2 || MpcB[i] == -2)
ME = E2;
sigma_b += ScB[i] *
Math.sin(DcB[i] * D + McB[i] * ME * M + MpcB[i] * Mprime
+ FcB[i] * F);
}
 
// Three intermediate arguments, in radians:
A1 = angle(119.75 + 131.849*T); // effect of Venus
A2 = angle(53.09 + 479264.29*T); // effect of Jupiter
A3 = angle(313.45 + 481266*T); // ??
 
 
sigma_l += 3958 * Math.sin(A1) + 1962 * Math.sin(Lprime - F)
+ 318 * Math.sin(A2);
sigma_b += 382 * Math.sin(A3) - 2235 * Math.sin(Lprime)
+ 175 * Math.sin(A1 - F) + 175 * Math.sin(A1 + F)
+ 127 * Math.sin(Lprime-Mprime) - 115 * Math.sin(Lprime+Mprime);
// Note: sigmas are still in DEGREES (actually millions of degrees)
 
// Finally, we can calculate the coordinates of the moon (in radians):
ecliptic[1] = modrad(Lprime + angle(sigma_l / 1000000));
 
ecliptic[2] = angle(sigma_b / 1000000);
ecliptic[3] = sigma_r / 1000 + 385000.56;
var phase=Math.abs((M-Mprime)/180)
if (phase>1){
phase=2-phase;
}
ecliptic[4] = D;
 
// alert(ecliptic[1]+' '+ecliptic[2]);
return ecliptic;
 
}
 
//
// Find the obliquity of the ecliptic in rads
//
 
function obliquity(day)	{
var t = day / 36525;
return angle((84381.448 + t*(-46.8150 + t*(-0.00059 + t*0.001813)))/ 3600);
}
 
//
// Nutation (good to 0.5 arcsec in dphi and 0.1 arcsec in depsilon)
// returns dphi and depsilon in an array in radians so the
// array can be added to ecliptic coords
//
 
function nutation(day) {
var t = day /36525;
var Om = angle(125.04452 - 1934.136261 * t);
var Lm = angle(218.3165 + 481267.8813 * t);
var Ls = angle(280.4665 + 36000.7698 * t);	
var dphi = -17.20 * Math.sin(Om) - 1.32 * Math.sin(2 * Ls)
- 0.23 * Math.sin(2 * Lm) + 0.21 * Math.sin(2 * Om);
var deps = 9.20 * Math.cos(Om) + 0.57 * Math.cos(2 * Ls)
+ 0.10 * Math.cos(2 * Lm) - 0.09 * Math.cos(2*Om);
nut = new Array(0, 0, 0, 0);
nut[1] = dphi * DEG2RAD /3600;
nut[2] = deps * DEG2RAD / 3600;
return nut
}
 
//
// Equatorial - converts ecliptic geocentrics to equatorial geocentrics
// uses obliquity of date, so if using J2000.0 coords, put day = 0
// array ecl has longitude, latitude and distance in [1], [2], [3]
//
 
function equatorial(day, ecl) {
var t = day /36525;
var obl = obliquity(day);
var co = Math.cos(obl);
var so = Math.sin(obl);
var sl = Math.sin(ecl[1]);
var Y = sl * co - Math.tan(ecl[2]) * so;
eq = new Array(0 ,0 , 0, 0);
eq[1] = ratan2(Y, Math.cos(ecl[1]));
var sd = Math.sin(ecl[2]) * co + Math.cos(ecl[2]) * so * sl;
eq[2] = Math.asin(sd);
eq[3] = ecl[3];
return eq;
}
 
//
// Some trig and astronomical functions - these functions
// all use local variables and do not read or write to
// form. They take arguments (arrays, numbers, strings)
// and return values (arrays, numbers, strings)
//
 
function ratan2(y, x) {
var a55;
if ((x == 0) && (y == 0)) {
return 0;
}
else	{
a55 = Math.atan(y / x);
if (x < 0) {
a55 = a55 + Math.PI;
}
if (y < 0 && x > 0) {
a55 = a55 + 2 * Math.PI;
}
return a55;
}
}
 
// convert degrees to a valid angle in radians:
function angle(deg)
{
return range(deg) * DEG2RAD;
}
 
function modrad(x)
{
var a9, b
b = x / TwoPi;
a9 = TwoPi * (b - ipart(b));
if (a9 < 0 ) {
a9 = a9 + TwoPi
}
return a9
}
 
//
// Calendar functions
// this is the usual days since J2000 function
// Gregorian only!!
//
 
function day2000(y8, m, d, h) {
var d1, b, c, greg;
greg = y8*10000 + m*100 + d;
if (m == 1 || m == 2) {
y8 = y8 - 1;
m = m + 12;
}
a = Math.floor(y8/100);
b = 2 - a + Math.floor(a/4)
c = Math.floor(365.25 * y8);
d1 = Math.floor(30.6001 * (m + 1));
return (b + c + d1 -730550.5 + d + h/24);
}
 
//
// Leap year detecting function (gregorian calendar)
// returns 1 for leap year and 0 for non-leap year
//
 
function isleap(y8) {
var a9, g;	
 
// assume not a leap year...
a9 = 0;
 
// ...flag leap year candidates...
if (y8 % 4 == 0) a9 = 1;
 
// ...if year is a century year then not leap...
if (y8 % 100 ==0 ) a9 = 0;
 
// ...except if century year divisible by 400...
if (y8 % 400 == 0) a9 = 1;
 
// ...and so done according to Pope Gregory's wishes.
return a9;
}
 
//
// Month and day number checking function
// This will work OK for Julian or Gregorian
// providing isleap() is defined appropriately
// Returns 1 if Month and Day combination OK,
// and 0 if month and day combination impossible
//
 
function goodmonthday(y, m, d) {
var a9, leap;
leap = isleap(y);
 
// assume OK
a9 = 1;
 
// first deal with zero day number!
if (d == 0) a9 = 0;
 
// Sort Feburary next
if ((m==2) && (leap ==1) && (d > 29)) a9= 0;
if ((m==2) && (d > 28) && (leap ==0))	a9 = 0;
 
// then the rest of the months - 30 days hath...
if(((m==4) || (m == 6) || (m == 9) || (m==11)) && d > 30) a9 = 0;
 
// ...31 days hath all the rest...
if (d > 31) a9 = 0;
 
// ...and so done
return a9;
}	
 
//
// Toolbox type things
//
 
 
function ipart(x) {
var a9;
if (x> 0) {
a9 = Math.floor(x);
}
else {
a9 = Math.ceil(x);
}	
return a9;
}
 
function fpart(x) {
return x % 1;
}
 
function range(x) {
var a9, b9
b9 = x / 360;
a9 = 360 * (b9 - ipart(b9));
if (a9 < 0 ) {
a9 = a9 + 360
}
return a9
}
 
//
// round rounds the number num to dp decimal places
// the second line is some C like jiggery pokery I
// found in an O'Reilly book which means if dp is null
// you get 2 decimal places.
//
function round(num, dp) {
// dp = (!dp ? 2: dp);
return Math.round (num * Math.pow(10, dp)) / Math.pow(10, dp);
}
 
//
// The close HTML comment below stops the script hiding, and then
// the close script and close head HTML tags return you to a normal
// web page.
//