// 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);

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);

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);

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);

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);

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);

// MoonPos calculates the Moon position, based on Meeus chapter 45
function MoonPos(year,month,day,hours) 
{
 var jd, T, T2, T3, T4;
 var LP, M, MP, F;
 var A1, A2, A3, E, E2;
 var Sl, Sr, Eterm, Sb;
 var mglong, mglat, mr, obl, ra, dec;
 var i;

  jd = julian(year,month,day,hours);
  T  = (jd-2451545.0)/36525;
  T2 = T*T;
  T3 = T2*T;
  T4 = T3*T;
  // Moons mean longitude L'
  LP = 218.3164477+481267.88123421*T-0.0015786*T2+T3/538841.0-T4/65194000.0;
  // Moons mean elongation
  D = 297.8501921+445267.1114034*T-0.0018819*T2+T3/545868.0-T4/113065000.0;
  // Suns mean anomaly
  M = 357.5291092+35999.0502909*T-0.0001536*T2+T3/24490000.0;
  // Moons mean anomaly M'
  MP = 134.9633964+477198.8675055*T+0.0087414*T2+T3/69699.0-T4/14712000.0;
  // Moons argument of latitude
  F = 93.2720950+483202.0175233*T-0.0036539*T2-T3/3526000.0+T4/863310000.0;

  // Additional arguments
  A1 = 119.75+131.849*T;
  A2 = 53.09+479264.290*T;
  A3 = 313.45+481266.484*T;
  E  = 1-0.002516*T-0.0000074*T2;
  E2 = E*E;
  // Sums of periodic terms from table 45.A and 45.B
  Sl=0.0; Sr=0.0;  
  for(i=0; i<60; i++) {
    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)); }

  Sb=0.0;
  for(i=0; i<60; i++) {
    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
  mglong = rev(LP+Sl/1000000.0);
  mglat  = rev(Sb/1000000.0);
  if(mglat > 180.0) mglat-=360;
  mr = Math.round(385000.56+Sr/1000.0);
  // Obliquity of Ecliptic
  obl = 23.4393-3.563E-9*(jd-2451543.5);
  // RA and dec
  ra  = rev(atan2d(sind(mglong)*cosd(obl)-tand(mglat)*sind(obl), cosd(mglong)))/15.0;
  dec = rev(asind(sind(mglat)*cosd(obl)+cosd(mglat)*sind(obl)*sind(mglong)));
  if(dec > 180.0) dec-=360;

 return new Array(ra,dec,mr);
}

// returns an array containing rise and set times (-1 if does not occur).
function MoonRise(year,month,day,latitude,longitude) 
{  
 // array to return -1 means not found       
 var riseset = new Array(-1,-1); 
 // elh is the elevation at the hour elhdone is true if elh calculated
 var elh     = new Array();     
 var elhdone = new Array();
 var hours, rad, altaz;
 var rise, found, hfirst, hlast;
 var elfirst, ellast, eld;
 
  for(i=0; i<=24; i++) elhdone[i] = false;     
  // Compute the moon elevation at start and end of day
  // store elevation at the hours in an array elh to save search time
  hours   = 0;
  rad     = MoonPos(year,month,day,hours);
  altaz   = radtoaa(rad[0],rad[1],year,month,day,hours,latitude,longitude);
  elh[0]  = altaz[0]; elhdone[0] = true; hours = 24;
  rad     = MoonPos(year,month,day,hours);
  altaz   = radtoaa(rad[0],rad[1],year,month,day,hours,latitude,longitude);
  elh[24] = altaz[0]; elhdone[24] = true;
  // search for moonrise and set
  for(rise=0; rise<2; rise++) {
     found = false; hfirst = 0; hlast = 24;
     // Try a binary chop on the hours to speed the search
     while(Math.ceil((hlast-hfirst)/2) > 1) {
        hmid = hfirst+Math.round((hlast-hfirst)/2);
        if(!elhdone[hmid]) {
           hours = hmid;
           rad   = MoonPos(year,month,day,hours);
           altaz = radtoaa(rad[0],rad[1],year,month,day,hours,latitude,longitude);
           elh[hmid] = altaz[0]; elhdone[hmid] = true; }
           
        if(((rise == 0) && (elh[hfirst] <= 0.0) && (elh[hmid] >= 0.0)) ||
           ((rise == 1) && (elh[hfirst] >= 0.0) && (elh[hmid] <= 0.0))) {
           hlast = hmid; found = true; continue; }
           
        if(((rise == 0) && (elh[hmid] <= 0.0) && (elh[hlast] >= 0.0)) ||
           ((rise == 1) && (elh[hmid] >= 0.0) && (elh[hlast] <= 0.0))) {
           hfirst = hmid; found = true; continue; }
        break; 
     }
     // If the binary chop did not find a 1 hour interval
     if((hlast-hfirst) > 1) 
        for(i=hfirst; i<hlast; i++) {
           found = false;
           if(!elhdone[i+1]) {
              hours = i+1;
              rad   = MoonPos(year,month,day,hours);
              altaz = radtoaa(rad[0],rad[1],year,month,day,hours,latitude,longitude);
              elh[hours] = altaz[0]; elhdone[hours] = true; }
              
           if(((rise == 0) && (elh[i] <= 0.0) && (elh[i+1] >= 0.0)) ||
              ((rise == 1) && (elh[i] >= 0.0) && (elh[i+1] <= 0.0))) {
              hfirst = i; hlast = i+1; found = true; break; } }
     // simple linear interpolation for the minutes
     if(found) {
       elfirst = elh[hfirst]; ellast = elh[hlast];
       hours = hfirst+0.5;
       rad   = MoonPos(year,month,day,hours);
       altaz = radtoaa(rad[0],rad[1],year,month,day,hours,latitude,longitude);
       if((rise == 0) && (altaz[0] <= 0.0)) { hfirst = hours; elfirst = altaz[0]; }
       if((rise == 0) && (altaz[0] > 0.0))  { hlast  = hours; ellast  = altaz[0]; }
       if((rise == 1) && (altaz[0] <= 0.0)) { hlast  = hours; ellast  = altaz[0]; }
       if((rise == 1) && (altaz[0] > 0.0))  { hfirst = hours; elfirst = altaz[0]; }
       eld = Math.abs(elfirst)+Math.abs(ellast);
       riseset[rise] = hfirst+(hlast-hfirst)*Math.abs(elfirst)/eld; }
  } 
    
 return(riseset);
}

// the illuminated percentage from Meeus chapter 46
function MoonPhase(year,month,day) 
{  
 var j, T, T2, T3, T4;
 var D, MP, M, pa, il;
        
  j  = dayno(year,month,day,12)+2451543.5;
  T  = (j-2451545.0)/36525;
  T2 = T*T; T3 = T2*T; T4 = T3*T;
  // Moons mean elongation Meeus second edition
  D = 297.8501921+445267.1114034*T-0.0018819*T2+T3/545868.0-T4/113065000.0; 
  // Moons mean anomaly M' Meeus second edition
  MP = 134.9633964+477198.8675055*T+0.0087414*T2+T3/69699.0-T4/14712000.0;
  // Suns mean anomaly
  M = 357.5291092+35999.0502909*T-0.0001536*T2+T3/24490000.0;
  // phase angle
  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);
  
  if(pa > 180.0) 
     pa = Math.floor(pa/360.0)*360.0 - pa;
  else if(pa < 0.0) 
     pa -= Math.floor(pa/360.0)*360.0;
         
  if(pa > 180.0)
     pa-=360.0;  
  // illumination percent
  il = 100.0*(1.0 + cosd(pa)) / 2.0;
  //il = 100.0*(1.0+cosd(pa)) / 2.0;
     
 return new Array(pa,il);
}

// Returns the dates for New and Full Moons.
// Converted from Basic Sky & Telescope, March, 1985
function moons(year) 
{
 var moondates=new Array();
 var R1, Y, K0, T, T2, T3;
 var J0, F0, M0, M1, B1, K9;
 var J, F, M5, M6, B6;
 
  R1 = Math.PI/180;
  Y  = year;
  K0 = Math.floor((Y-1900)*12.3685);
  T  = (Y-1899.5)/100; T2 = T*T; T3 = T*T*T;
  J0 = 2415020+29*K0;
  F0 = 0.0001178*T2-0.000000155*T3;
  F0 = F0+0.75933+0.53058868*K0;
  F0 = F0-0.000837*T-0.000335*T2;
  M0 = K0*0.08084821133;
  M0 = 360*(M0-Math.floor(M0))+359.2242;
  M0 -= 0.0000333*T2;
  M0 -= 0.00000347*T3;
  M1 = K0*0.07171366128;
  M1 = 360*(M1-Math.floor(M1))+306.0253;
  M1 += 0.0107306*T2;
  M1 += 0.00001236*T3;
  B1 = K0*0.08519585128;
  B1 = 360*(B1-Math.floor(B1))+21.2964;
  B1 -= 0.0016528*T2;
  B1 -= 0.00000239*T3;
  
  for(K9=0; K9<=28; K9++) {
    J  = J0+14*K9;
    F  = F0+0.765294*K9;
    K  = K9/2;
    M5 = (M0+K*29.10535608)*R1;
    M6 = (M1+K*385.81691806)*R1;
    B6 = (B1+K*390.67050646)*R1;
    F -= 0.4068*Math.sin(M6);
    F += (0.1734-0.000393*T)*Math.sin(M5);
    F += 0.0161*Math.sin(2*M6);
    F += 0.0104*Math.sin(2*B6);
    F -= 0.0074*Math.sin(M5-M6);
    F -= 0.0051*Math.sin(M5+M6);
    F += 0.0021*Math.sin(2*M5);
    F += 0.0010*Math.sin(2*B6-M6);
    J  = J+F;
    moondates[K9] = J;
  }

 return moondates;
}

-->