/*

Code from http://www.srrb.noaa.gov/highlights/sunrise/sunrise.html

*/

function isLeapYear(yr) 
{
  return ((yr % 4 == 0 && yr % 100 != 0) || yr % 400 == 0);
}

function radToDeg(angleRad) 
{
  return (180.0 * angleRad / Math.PI);
}

function degToRad(angleDeg) 
{
  return (Math.PI * angleDeg / 180.0);
}

function calcDayOfYear(mn, dy, lpyr) 
{
  var k = (lpyr ? 1 : 2);
  var doy = Math.floor((275 * mn)/9) - k * Math.floor((mn + 9)/12) + dy -30;
  return doy;
}

function calcDayOfWeek(juld)
{
  return (Math.floor(juld+0.5)%7)+1;
}

function calcJD(year, month, day)
{
  if (month <= 2) 
  {
    year -= 1;
    month += 12;
  }
  var A = Math.floor(year/100);
  var B = 2 - A + Math.floor(A/4);

  var JD = Math.floor(365.25*(year + 4716)) + Math.floor(30.6001*(month+1)) + day + B - 1524.5;
  return JD;
}

function calcTimeJulianCent(jd)
{
  var T = (jd - 2451545.0)/36525.0;
  return T;
}

function calcJDFromJulianCent(t)
{
  var JD = t * 36525.0 + 2451545.0;
  return JD;
}

function calcGeomMeanLongSun(t)
{
  var L0 = 280.46646 + t * (36000.76983 + 0.0003032 * t);
  while(L0 > 360.0)
  {
    L0 -= 360.0;
  }
  while(L0 < 0.0)
  {
    L0 += 360.0;
  }
  return L0;		// in degrees
}

function calcGeomMeanAnomalySun(t)
{
  var M = 357.52911 + t * (35999.05029 - 0.0001537 * t);
  return M;		// in degrees
}

function calcEccentricityEarthOrbit(t)
{
  var e = 0.016708634 - t * (0.000042037 + 0.0000001267 * t);
  return e;		// unitless
}

function calcSunEqOfCenter(t)
{
  var m = calcGeomMeanAnomalySun(t);

  var mrad = degToRad(m);
  var sinm = Math.sin(mrad);
  var sin2m = Math.sin(mrad+mrad);
  var sin3m = Math.sin(mrad+mrad+mrad);

  var C = sinm * (1.914602 - t * (0.004817 + 0.000014 * t)) + sin2m * (0.019993 - 0.000101 * t) + sin3m * 0.000289;
  return C;		// in degrees
}

function calcSunTrueLong(t)
{
  var l0 = calcGeomMeanLongSun(t);
  var c = calcSunEqOfCenter(t);

  var O = l0 + c;
  return O;		// in degrees
}

function calcSunTrueAnomaly(t)
{
  var m = calcGeomMeanAnomalySun(t);
  var c = calcSunEqOfCenter(t);

  var v = m + c;
  return v;		// in degrees
}

function calcSunRadVector(t)
{
  var v = calcSunTrueAnomaly(t);
  var e = calcEccentricityEarthOrbit(t);
 
  var R = (1.000001018 * (1 - e * e)) / (1 + e * Math.cos(degToRad(v)));
  return R;		// in AUs
}

function calcSunApparentLong(t)
{
  var o = calcSunTrueLong(t);

  var omega = 125.04 - 1934.136 * t;
  var lambda = o - 0.00569 - 0.00478 * Math.sin(degToRad(omega));
  return lambda;		// in degrees
}

function calcMeanObliquityOfEcliptic(t)
{
  var seconds = 21.448 - t*(46.8150 + t*(0.00059 - t*(0.001813)));
  var e0 = 23.0 + (26.0 + (seconds/60.0))/60.0;
  return e0;		// in degrees
}

function calcObliquityCorrection(t)
{
  var e0 = calcMeanObliquityOfEcliptic(t);

  var omega = 125.04 - 1934.136 * t;
  var e = e0 + 0.00256 * Math.cos(degToRad(omega));
  return e;		// in degrees
}

function calcSunRtAscension(t)
{
  var e = calcObliquityCorrection(t);
  var lambda = calcSunApparentLong(t);
 
  var tananum = (Math.cos(degToRad(e)) * Math.sin(degToRad(lambda)));
  var tanadenom = (Math.cos(degToRad(lambda)));
  var alpha = radToDeg(Math.atan2(tananum, tanadenom));
  return alpha;		// in degrees
}

function calcSunDeclination(t)
{
  var e = calcObliquityCorrection(t);
  var lambda = calcSunApparentLong(t);

  var sint = Math.sin(degToRad(e)) * Math.sin(degToRad(lambda));
  var theta = radToDeg(Math.asin(sint));
  return theta;		// in degrees
}

function calcEquationOfTime(t)
{
  var epsilon = calcObliquityCorrection(t);
  var l0 = calcGeomMeanLongSun(t);
  var e = calcEccentricityEarthOrbit(t);
  var m = calcGeomMeanAnomalySun(t);

  var y = Math.tan(degToRad(epsilon)/2.0);
  y *= y;

  var sin2l0 = Math.sin(2.0 * degToRad(l0));
  var sinm   = Math.sin(degToRad(m));
  var cos2l0 = Math.cos(2.0 * degToRad(l0));
  var sin4l0 = Math.sin(4.0 * degToRad(l0));
  var sin2m  = Math.sin(2.0 * degToRad(m));

  var Etime = y * sin2l0 - 2.0 * e * sinm + 4.0 * e * y * sinm * cos2l0
  		- 0.5 * y * y * sin4l0 - 1.25 * e * e * sin2m;

  return radToDeg(Etime)*4.0;	// in minutes of time
}

function calcHourAngleSunrise(lat, solarDec)
{
  var latRad = degToRad(lat);
  var sdRad  = degToRad(solarDec)

  var HAarg = (Math.cos(degToRad(90.833))/(Math.cos(latRad)*Math.cos(sdRad))-Math.tan(latRad) * Math.tan(sdRad));

  var HA = (Math.acos(Math.cos(degToRad(90.833))/(Math.cos(latRad)*Math.cos(sdRad))-Math.tan(latRad) * Math.tan(sdRad)));

  return HA;		// in radians
}

function calcHourAngleSunset(lat, solarDec)
{
  var latRad = degToRad(lat);
  var sdRad  = degToRad(solarDec)

  var HAarg = (Math.cos(degToRad(90.833))/(Math.cos(latRad)*Math.cos(sdRad))-Math.tan(latRad) * Math.tan(sdRad));

  var HA = (Math.acos(Math.cos(degToRad(90.833))/(Math.cos(latRad)*Math.cos(sdRad))-Math.tan(latRad) * Math.tan(sdRad)));

  return -HA;		// in radians
}

function calcSunriseUTC(JD, latitude, longitude)
{
  var t = calcTimeJulianCent(JD);

  // *** Find the time of solar noon at the location, and use
        //     that declination. This is better than start of the 
        //     Julian day

  var noonmin = calcSolNoonUTC(t, longitude);
  var tnoon = calcTimeJulianCent (JD+noonmin/1440.0);

  // *** First pass to approximate sunrise (using solar noon)

  var eqTime = calcEquationOfTime(tnoon);
  var solarDec = calcSunDeclination(tnoon);
  var hourAngle = calcHourAngleSunrise(latitude, solarDec);

  var delta = longitude - radToDeg(hourAngle);
  var timeDiff = 4 * delta;	// in minutes of time
  var timeUTC = 720 + timeDiff - eqTime;	// in minutes

  // alert("eqTime = " + eqTime + "\nsolarDec = " + solarDec + "\ntimeUTC = " + timeUTC);

  // *** Second pass includes fractional jday in gamma calc

  var newt = calcTimeJulianCent(calcJDFromJulianCent(t) + timeUTC/1440.0); 
  eqTime = calcEquationOfTime(newt);
  solarDec = calcSunDeclination(newt);
  hourAngle = calcHourAngleSunrise(latitude, solarDec);
  delta = longitude - radToDeg(hourAngle);
  timeDiff = 4 * delta;
  timeUTC = 720 + timeDiff - eqTime; // in minutes

  // alert("eqTime = " + eqTime + "\nsolarDec = " + solarDec + "\ntimeUTC = " + timeUTC);

  return timeUTC;
}

function calcSolNoonUTC(t, longitude)
{
  // First pass uses approximate solar noon to calculate eqtime
  var tnoon = calcTimeJulianCent(calcJDFromJulianCent(t) + longitude/360.0);
  var eqTime = calcEquationOfTime(tnoon);
  var solNoonUTC = 720 + (longitude * 4) - eqTime; // min

  var newt = calcTimeJulianCent(calcJDFromJulianCent(t) -0.5 + solNoonUTC/1440.0); 

  eqTime = calcEquationOfTime(newt);
  // var solarNoonDec = calcSunDeclination(newt);
  solNoonUTC = 720 + (longitude * 4) - eqTime; // min
  
  return solNoonUTC;
}

function calcSunsetUTC(JD, latitude, longitude)
{
  var t = calcTimeJulianCent(JD);

  // *** Find the time of solar noon at the location, and use
        //     that declination. This is better than start of the 
        //     Julian day

  var noonmin = calcSolNoonUTC(t, longitude);
  var tnoon = calcTimeJulianCent (JD+noonmin/1440.0);

  // First calculates sunrise and approx length of day

  var eqTime = calcEquationOfTime(tnoon);
  var solarDec = calcSunDeclination(tnoon);
  var hourAngle = calcHourAngleSunset(latitude, solarDec);

  var delta = longitude - radToDeg(hourAngle);
  var timeDiff = 4 * delta;
  var timeUTC = 720 + timeDiff - eqTime;

  // first pass used to include fractional day in gamma calc

  var newt = calcTimeJulianCent(calcJDFromJulianCent(t) + timeUTC/1440.0); 
  eqTime = calcEquationOfTime(newt);
  solarDec = calcSunDeclination(newt);
  hourAngle = calcHourAngleSunset(latitude, solarDec);

  delta = longitude - radToDeg(hourAngle);
  timeDiff = 4 * delta;
  timeUTC = 720 + timeDiff - eqTime; // in minutes

  return timeUTC;
}

function isDST(year, month, day)
{
  if(month<3 || month>10)
  {
    return 0;
  }
  else if(month>3 && month<10)
  {
    return 1;
  }
  else
  {
    lastSunday=31-(calcDayOfWeek(calcJD(year, month, 31))%7);
    return (day>=lastSunday)==(month==3);
  }
}

function writeTime(time)
{
  min=Math.floor(time+0.5);
  hour=Math.floor(min/60);
  min=min%60;
  f="";
  f += Math.floor(hour/10);
  f += hour%10;
  f += ":";
  f += Math.floor(min/10);
  f += min%10;
  return f;
}

function sunRiseAndSet(lang)
{
  longitude = -12.4832;
  latitude = 55.7606;

  UTCoffset = 60;
  DSToffset = 120;

  var d = new Date();
  year = d.getFullYear();
  month = d.getMonth()+1;
  day = d.getDate();

  JD = calcJD(year,month,day);
  if(isDST(year,month,day))
  {
    timeOffset = DSToffset;
  }
  else
  {
    timeOffset = UTCoffset;
  }

  document.getElementById('riseandset').innerHTML = (lang==0?'Sol: ':'Sun: ')+writeTime(calcSunriseUTC(JD, latitude, longitude)+timeOffset)+
                 (lang==0?' til: ':' to: ')+writeTime(calcSunsetUTC(JD, latitude, longitude)+timeOffset);
}