// JavaScript by Peter Hayes
// http://www.aphayes.pwp.blueyonder.co.uk/ or
// http://www.peter-hayes.freeserve.co.uk/
// Copyright 2001-2004
// This code is made freely available but please keep this notice.
// The calculations are approximate but should be good enough for general use,
// I accept no responsibility for errors in astronomy or coding.

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 julian(year,month,day,hours) {
  // Returns JD given the MJD
   return  dayno(year,month,day,hours)+2451543.5
}

function s_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 SunRiseSet(year,month,day,latitude,longitude) {
  // Based on method in sci.astro FAQ by Paul Schlyter
  // returns an array holding rise and set times in UTC hours
  // NOT accurate at high latitudes 
  // latitude = your local latitude: north positive, south negative
  // longitude = your local longitude: east positive, west negative
  // Calculate the Suns position at noon local zone time
  var d=dayno(year,month,day,12.0-longitude/15);
  var oblecl=23.4393-3.563E-7*d;
  var w=282.9404+4.70935E-5*d;
  var M=356.0470+0.9856002585*d;
  var e=0.016709-1.151E-9*d;
  var E=M+e*(180/Math.PI)*sind(M)*(1.0+e*cosd(M));
  var A=cosd(E)-e;
  var B=Math.sqrt(1-e*e)*sind(E);
  var slon=w+atan2d(B,A);
  var sRA=atan2d(sind(slon)*cosd(oblecl),cosd(slon));
  while(sRA<0)sRA+=360; while(sRA>360)sRA-=360; sRA=sRA/15;
  var sDec=asind(sind(oblecl)*sind(slon));
  // Time sun is on the meridian
  var lst=s_local_sidereal(year,month,day,12-longitude/15,longitude);
  var MT=12.0-longitude/15+sRA-lst;
  while(MT<0)MT+=24; while(MT>24)MT-=24;
  // hour angle
  var cHA0=(sind(-0.8333333)-sind(latitude)*sind(sDec))/(cosd(latitude)*cosd(sDec));
  var HA0=acosd(cHA0);
  HA0=rev(HA0)/15;
  // return rise and set times
  return new Array((MT-HA0),(MT+HA0));
}