// JavaScript by Peter Hayes
// http://www.aphayes.pwp.blueyonder.co.uk/
// Functions for the moon
// Meeus first edition table 45.A Longitude and distance of the moon
var T45AD=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);
var T45AM=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);
var T45AMP=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);
var T45AF=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);
var T45AL=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,0);
var T45AR=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);
// Meeus table 45B latitude of the moon
var T45BD=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);
var T45BM=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);
var T45BMP=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);
var T45BF=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);
var T45BL=new Array(5128122,280602,277693,173237,55413,46271,32573,17198,9266,8822,8216,4324,4200,-3359,2463,2211,2065,-1870,1828,-1794,-1749,-1565,-1491,-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);
// Support Functions
function getFullYear(now) {
  var year = now.getYear();
  if (year==0) {year=2000};
  if (year<1900) {year=year+1900};
  return year;
}
function rev(angle){return angle-Math.floor(angle/360.0)*360.0;}
function sind(angle){return Math.sin((angle*Math.PI)/180.0);}
function cosd(angle){return Math.cos((angle*Math.PI)/180.0);}
function tand(angle){return Math.tan((angle*Math.PI)/180.0);}
function asind(c){return (180.0/Math.PI)*Math.asin(c);}
function acosd(c){return (180.0/Math.PI)*Math.acos(c);}
function atan2d(y,x){return (180.0/Math.PI)*Math.atan(y/x)-180.0*(x<0);}
function julian(year,month,day,hours) {return  dayno(year,month,day,hours)+2451543.5}
function dayno(year,month,day,hours) {
  // Day number is a modified Julian date,day 0 is 2000 January 0.0
  // which corresponds to a Julian date of 2451543.5
  var d= 367*year-Math.floor(7*(year+Math.floor((month+9)/12))/4)
         +Math.floor((275*month)/9)+day-730530+hours/24;
  return d;
}
function jdtocd(jd) {
  // The calendar date from julian date
  // Returns year, month, day, day of week, s, minutes, seconds
  var Z=Math.floor(jd+0.5); var F=jd+0.5-Z;
  if (Z < 2299161){
    var A=Z;
  } else {
    var alpha=Math.floor((Z-1867216.25)/36524.25);
    var A=Z+1+alpha-Math.floor(alpha/4);
  }
  var B=A+1524; var C=Math.floor((B-122.1)/365.25); var D=Math.floor(365.25*C);
  var E=Math.floor((B-D)/30.6001); var d=B-D-Math.floor(30.6001*E)+F;
  if (E < 14) {var month=E-1;} else {var month=E-13;}
  if ( month>2) {var year=C-4716;} else {var year=C-4715;}
  var day=Math.floor(d); var h=(d-day)*24; var hours=Math.floor(h);
  var m=(h-hours)*60; var minutes=Math.floor(m); var seconds=Math.round((m-minutes)*60);
  if (seconds >= 60) {minutes=minutes+1;seconds=seconds-60;}
  if (minutes >= 60) {hours=hours+1;minutes=0;}
  var dw=Math.floor(jd+1.5)-7*Math.floor((jd+1.5)/7);
  return new Array(year,month,day,dw,hours,minutes,seconds);  
}

function computeJD() {
var mydate=new Date(); var YY=mydate.getFullYear();
if (YY < 1000){YY+=1900}
var DD=mydate.getDay()
var MM=mydate.getMonth()
var DD=mydate.getDate()
    with (Math) {  
      GGG = 1;
      if (YY <= 1585) GGG = 0;
      JD = -1 * floor(7 * (floor((MM + 9) / 12) + YY) / 4);
      S = 1;
      if ((MM - 9)<0) S=-1;
      A = abs(MM - 9);
      J1 = floor(YY + S * floor(A / 7));
      J1 = -1 * floor((floor(J1 / 100) + 1) * 3 / 4);
      JD = JD + floor(275 * MM / 9) + DD + (GGG * J1);
      JD = JD + 1721027 + 2 * GGG + 367 * YY;
    }
      return JD;
}

function computeDJ(JD) {
  with (Math) {
    IJD = floor(JD);
    L = floor(IJD + 68569);
    N = floor(4 * L / 146097);
    L = L - floor((146097*N + 3)/4);
    I = floor(4000*(L + 1)/1461001);
    L = L - floor(1461 * I / 4) + 31;
    J = floor(80 * L / 2447);
    K = L - floor(2447 * J / 80);
    L = floor(J/11);
    J = J + 2 - 12*L;
    I = 100*(N - 49) + I + L;
  }
  if (JD != 0) {
   return (I+","+J+","+K)
  }
}

function local_sidereal(year,month,day,hours,lon) {
  // Compute local siderial time in degrees
  // year, month, day and hours are the Greenwich date and time
  // lon is the observers longitude
  var d=dayno(year,month,day,hours);
  var lst=(98.9818+0.985647352*d+hours*15+lon);
  return rev(lst)/15;
}

function hmstring(t) {
// takes hours and returns hh:mm
  var hours=t;
  while (hours >= 24) {hours-=24;}
  while (hours < 0) {hours+=24;}
  var minutes=Math.round(60.0*(hours-Math.floor(hours)));
  hours=Math.floor(hours);
  if (minutes >= 60) {hours+=1; minutes-=60;}
  if (hours >= 24) {hours-=24;}
  var hmstr=(t < 0) ? "-" : "";
  hmstr+=((hours < 10) ? "0" : "")+hours;
  hmstr+=((minutes < 10) ? ":0" : ":")+minutes;
  return hmstr;
}

function MoonPos(year,month,day,hours) {
  // julian date
	var jd=julian(year,month,day,hours);
  var T=(jd-2451545.0)/36525; var T2=T*T; var T3=T2*T; var T4=T3*T;
  // Moons mean longitude L'
  var LP=218.3164477+481267.88123421*T-0.0015786*T2+T3/538841.0-T4/65194000.0;
  // Moons mean elongation
  var D=297.8501921+445267.1114034*T-0.0018819*T2+T3/545868.0-T4/113065000.0;
  // Suns mean anomaly
  var M=357.5291092+35999.0502909*T-0.0001536*T2+T3/24490000.0;
  // Moons mean anomaly M'
  var MP=134.9633964+477198.8675055*T+0.0087414*T2+T3/69699.0-T4/14712000.0;
  // Moons argument of latitude
  var F=93.2720950+483202.0175233*T-0.0036539*T2-T3/3526000.0+T4/863310000.0;
  // Additional arguments
  var A1=119.75+131.849*T; var A2=53.09+479264.290*T;
  var A3=313.45+481266.484*T; var E=1-0.002516*T-0.0000074*T2; var E2=E*E;
  // Sums of periodic terms from table 45.A and 45.B
  var Sl=0.0; var Sr=0.0;
  for (var i=0; i<60; i++) {
    var Eterm=1;
    if (Math.abs(T45AM[i])==1) Eterm=E;
    if (Math.abs(T45AM[i])==2) Eterm=E2;
    Sl+=T45AL[i]*Eterm*sind(rev(T45AD[i]*D+T45AM[i]*M+T45AMP[i]*MP+T45AF[i]*F));
    Sr+=T45AR[i]*Eterm*cosd(rev(T45AD[i]*D+T45AM[i]*M+T45AMP[i]*MP+T45AF[i]*F));
  }
  var Sb=0.0;
  for (var i=0; i<60; i++) {
    var Eterm=1;
    if (Math.abs(T45BM[i])==1) Eterm=E;
    if (Math.abs(T45BM[i])==2) Eterm=E2;
    Sb+=T45BL[i]*Eterm*sind(rev(T45BD[i]*D+T45BM[i]*M+T45BMP[i]*MP+T45BF[i]*F));
  }
  // Additional additive terms
  Sl=Sl+3958*sind(rev(A1))+1962*sind(rev(LP-F))+318*sind(rev(A2));
  Sb=Sb-2235*sind(rev(LP))+382*sind(rev(A3))+175*sind(rev(A1-F))+
     175*sind(rev(A1+F))+127*sind(rev(LP-MP))-115*sind(rev(LP+MP));
  // geocentric longitude,latitude and distance
  var mglong=rev(LP+Sl/1000000.0);
  var mglat=rev(Sb/1000000.0);
  if (mglat > 180.0) mglat=mglat-360;
  
   var mr=Math.round(385000.56+Sr/1000.0);
  // Obliquity of Ecliptic
  var obl=23.4393-3.563E-9*(jd-2451543.5);
  // RA and dec
  var ra=rev(atan2d(sind(mglong)*cosd(obl)-tand(mglat)*sind(obl),
                    cosd(mglong)))/15.0;
  var dec=rev(asind(sind(mglat)*cosd(obl)+cosd(mglat)*sind(obl)*sind(mglong)));
  if (dec > 180.0) dec=dec-360;
  return new Array(ra,dec,mr,mglong,mglat,obl);
}


function radtoaa(ra,dec,year,month,day,hours,lat,lon) {
  // convert ra and dec to altitude and azimuth
  // year,month, day and hours are the Greenwich date and time
  // lat and lon are the observers latitude and longitude
  var lst=local_sidereal(year,month,day,hours,lon);
  var x=cosd(15.0*(lst-ra))*cosd(dec);
  var y=sind(15.0*(lst-ra))*cosd(dec); var z=sind(dec);
 
  // rotate so z is the local zenith
  var xhor=x*sind(lat)-z*cosd(lat); var yhor=y; var zhor=x*cosd(lat)+z*sind(lat);
  var azimuth=rev(atan2d(yhor,xhor)+180.000); // so 0 degrees is north
  var altitude=atan2d(zhor,Math.sqrt(xhor*xhor+yhor*yhor));
  var lstg=local_sidereal(year,month,day,hours,0);
  var xg=cosd(15.00*(lstg-ra))*cosd(dec); 
  var yg=sind(15.00*(lstg-ra))*cosd(dec); var zg=sind(dec);
  var xhorg=xg*sind(0)-zg*cosd(0); var yhorg=yg; var zhorg=xg*cosd(0.000)+zg*sind(0.000);
  var azimuthg=rev(atan2d(yhorg,xhorg)+180.000); // so 0 degrees is north
  var xg=acosd(cosd(15.000*(lstg-ra)));
   if (azimuthg<360.000&&azimuthg>180.000){xg=-xg}; 
  return new Array(altitude,azimuth,xg);
}
function MoonPhase(year,month,day,hour) {
  // the illuminated percentage from Meeus chapter 46
  var j=dayno(year,month,day,hour)+2451543.5;
  var T=(j-2451545.0)/36525;
  var T2=T*T;
  var T3=T2*T;
  var T4=T3*T;
  // Moons mean elongation Meeus first edition
  // var D=297.8502042+445267.1115168*T-0.0016300*T2+T3/545868.0-T4/113065000.0;
  // Moons mean elongation Meeus second edition
  var D=297.8501921+445267.1114034*T-0.0018819*T2+T3/545868.0-T4/113065000.0; 
  // Moons mean anomaly M' Meeus first edition
  // var MP=134.9634114+477198.8676313*T+0.0089970*T2+T3/69699.0-T4/14712000.0;
  // Moons mean anomaly M' Meeus second edition
  var MP=134.9633964+477198.8675055*T+0.0087414*T2+T3/69699.0-T4/14712000.0;
  // Suns mean anomaly
  var M=357.5291092+35999.0502909*T-0.0001536*T2+T3/24490000.0;
  // phase angle
  var pa=180.0-D-6.289*sind(MP)+2.1*sind(M)-1.274*sind(2*D-MP)
         -0.658*sind(2*D)-0.214*sind(2*MP)-0.11*sind(D);
  return(rev(pa));
}