/* * * * * * * * * *    SUNRISET SCRIPT    * * * * * * * * * *\
*                                                             *
*     Computes Sun rise/set times, start/end of twilight.     *
*                                                             *
*        Originally written as DAYLEN.C, 1989-08-16           *
*            Modified to SUNRISET.C, 1992-12-01               *
*              (c) Paul Schlyter, 1989, 1992                  *
*                                                             *
*  Released to the public domain by Paul Schlyter, Dec. 1992  *
*                                                             *
*  --  Modified by Bo Johansson  1995-05-12 , 1995-10-14  --  *
*  --  Ported to JavaScript by Bo Johansson  1998-05-22   --  *
*  --  Very minor revision by Bo Johansson  1999-10-25    --  *
*  --   Very minor Y2K fix by Bo Johansson  2000-01-03    --  *
*  --  Minor Opera 4 bug fix by Bo Johansson  2000-10-14  --  *
*                                                             *
*  --  http://w1.545.telia.com/%7Eu54504162/index_e.htm   --  *
\* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */

// Global variables:
  
var twAngle = -6.0;    // For civil twilight, set to -12.0 for
                       // nautical, and -18.0 for astr. twilight

var srAngle = -35.0/60.0;    // For sunrise/sunset

var sRiseT;      // For sunrise/sunset times
var sSetT;
var srStatus;

var twStartT;    // For twilight times
var twEndT;
var twStatus;

var sDIST;       // Solar distance, astronomical units
var sRA;         // Sun's Right Ascension
var sDEC;        // Sun's declination

var sLON;        // True solar longitude


//-------------------------------------------------------------
// A function to compute the number of days elapsed since
// 2000 Jan 0.0  (which is equal to 1999 Dec 31, 0h UT)
//-------------------------------------------------------------
function dayDiff2000(y, m, d)
{
  return (367.0 * y -((7 *(y +((m + 9) / 12))) / 4) +
          ((275 * m)/9) + d - 730530.0);
}


//-------------------------------------------------------------
// Some conversion factors between radians and degrees

var RADEG = 180.0 / Math.PI;
var DEGRAD = Math.PI / 180.0;


//-------------------------------------------------------------
// The trigonometric functions in degrees

function sind(x) { return Math.sin(x * DEGRAD); }
function cosd(x) { return Math.cos(x * DEGRAD); }
function tand(x) { return Math.tan(x * DEGRAD); }

function atand(x) { return (RADEG * Math.atan(x)); }
function asind(x) { return (RADEG * Math.asin(x)); }
function acosd(x) { return (RADEG * Math.acos(x)); }

//function atan2d(y,x) { return (RADEG*Math.atan2(y,x)); }

function atan2d(y,x)          //###  Hack to work in MSIE3 ###
{
  var at2 = (RADEG * Math.atan(y/x));

  if( x < 0 && y < 0)
    at2 -= 180;

  else if( x < 0 && y > 0)
    at2 += 180;

  return at2;
}


//-------------------------------------------------------------
// This function computes times for sunrise/sunset. 
// Sunrise/sunset is considered to occur when the Sun's upper
// limb is 35 arc minutes below the horizon (this accounts for
// the refraction of the Earth's atmosphere).
//-------------------------------------------------------------
function sunRiseSet(year, month, day, lat, lon)
{
  return (sunTimes(year, month, day, lat, lon, srAngle, 1));
}



//-------------------------------------------------------------
// This function computes the start and end times of civil
// twilight. Civil twilight starts/ends when the Sun's center
// is 6 degrees below the horizon.
//-------------------------------------------------------------
function civTwilight(year, month, day, lat, lon)
{
  return (sunTimes( year, month, day, lat, lon, twAngle, 0));
}



//-------------------------------------------------------------
// The "workhorse" function for sun rise/set times
//
// year,month,date = calendar date, 1801-2099 only.
// Eastern longitude positive, Western longitude negative
// Northern latitude positive, Southern latitude negative
//
// altit = the altitude which the Sun should cross. Set to
//         -35/60 degrees for rise/set, -6 degrees for civil,
//         -12 degrees for nautical and -18 degrees for
//         astronomical twilight.
//
// sUppLimb: non-zero -> upper limb, zero -> center. Set to
//           non-zero (e.g. 1) when computing rise/set times,
//           and to zero when computing start/end of twilight.
//
// Status:  0 = sun rises/sets this day, times stored in Global
//              variables
//          1 = sun above the specified "horizon" 24 hours.
//              Rising set to time when the sun is at south,
//              minus 12 hours while Setting is set to the
//              south time plus 12 hours.
//         -1 = sun is below the specified "horizon" 24 hours.
//              Rising and Setting are both set to the time
//              when the sun is at south.
//-------------------------------------------------------------
function sunTimes(year, month, day, lat, lon, altit, sUppLimb)
{
  var dayDiff;    // Days since 2000 Jan 0.0 (negative before)
  var sRadius;    // Sun's apparent radius
  var diuArc;     // Diurnal arc
  var sSouthT;    // Time when Sun is at south
  var locSidT;    // Local sidereal time

  var stCode = 0;     // Status code from function - usually 0

      /* Compute dayDiff of 12h local mean solar time */
      
  dayDiff = dayDiff2000(year, month, day) + 0.5 - lon/360.0;

      /* Compute local sideral time of this moment */

  locSidT = revolution( GMST0(dayDiff) + 180.0 + lon );

      /* Compute Sun's RA + Decl at this moment */

  sunRaDec( dayDiff );

      /* Compute time when Sun is at south - in hours UT */
      

  sSouthT = 12.0 - rev180(locSidT - sRA)/15.0;


      /* Compute the Sun's apparent radius, degrees */

  sRadius = 0.2666 / sDIST;


      /* Do correction to upper limb, if necessary */

  if ( sUppLimb != 0)            // For surise/sunset
    altit -= sRadius;


      /* Compute the diurnal arc that the Sun traverses */
      /* to reach the specified altitude altit:         */

    var cost = ( sind(altit) - sind(lat) * sind(sDEC) ) /
               ( cosd(lat) * cosd(sDEC) );

    if ( cost >= 1.0 )
    {
       stCode = -1;
       diuArc = 0.0;       // Sun always below altit
    }

    else if ( cost <= -1.0 )
    {
       stCode = 1;
       diuArc = 12.0;      // Sun always above altit
    }

    else
    {                             // *** Opera bug fix 2000-10-14
                                  // The diurnal arc, hours
       diuArc = acosd(cost)/15.0;
    }


      /* Compute rise and set times - in hours UT */

  if ( sUppLimb != 0)       // For sunrise/sunset
  {
    sRiseT = sSouthT - diuArc;

    if(sRiseT < 0)          // Sunrise day before
      sRiseT += 24;
      
    sSetT  = sSouthT + diuArc;

    if(sSetT > 24)          // Sunset next day
      sSetT -= 24;

    srStatus = stCode;
  }

  else                      // For twilight times
  {
    twStartT = sSouthT - diuArc;

    if(twStartT < 0)
      twStartT += 24;
      
    twEndT  = sSouthT + diuArc;

    if(twEndT > 24)
      twEndT -= 24;

    twStatus = stCode;
  }

} //================== sunTimes() =====================


//-------------------------------------------------------------
// This function computes the sun's spherical coordinates
//-------------------------------------------------------------
function sunRaDec(dayDiff)
{
  var eclObl;   // Obliquity of ecliptic
                // (inclination of Earth's axis)
  var x;
  var y;
  var z;

      /* Compute Sun's ecliptical coordinates */

  sunPos( dayDiff );

      /* Compute ecliptic rectangular coordinates (z=0) */

  x = sDIST * cosd(sLON);

  y = sDIST * sind(sLON);

      /* Compute obliquity of ecliptic */
      /* (inclination of Earth's axis) */

//  eclObl = 23.4393 - 3.563E-7 * dayDiff;

  eclObl = 23.4393 - 3.563/10000000 * dayDiff;  // for Opera

      /* Convert to equatorial rectangular coordinates */
      /* - x is uchanged                               */

  z = y * sind(eclObl);
  y = y * cosd(eclObl);

      /* Convert to spherical coordinates */

  sRA = atan2d( y, x );
  sDEC = atan2d( z, Math.sqrt(x*x + y*y) );

} //================= sunRaDec() =======================



//-------------------------------------------------------------
// Computes the Sun's ecliptic longitude and distance
// at an instant given in dayDiff, number of days since
// 2000 Jan 0.0.  The Sun's ecliptic latitude is not
// computed, since it's always very near 0. 
//-------------------------------------------------------------
function sunPos(dayDiff)
{
  var M;       // Mean anomaly of the Sun
  var w;       // Mean longitude of perihelion
               // Note: Sun's mean longitude = M + w
  var e;       // Eccentricity of Earth's orbit
  var eAN;     // Eccentric anomaly
  var x;       // x, y coordinates in orbit
  var y;
  var v;       // True anomaly

      /* Compute mean elements */

  M = revolution( 356.0470 + 0.9856002585 * dayDiff );

//  w = 282.9404 + 4.70935E-5 * dayDiff;
//  e = 0.016709 - 1.151E-9 * dayDiff;

  w = 282.9404 + 4.70935/100000 * dayDiff;    // for Opera
  e = 0.016709 - 1.151/1000000000 * dayDiff;  // for Opera

      /* Compute true longitude and radius vector */

  eAN = M + e * RADEG * sind(M) * ( 1.0 + e * cosd(M) );

  x = cosd(eAN) - e;
  y = Math.sqrt( 1.0 - e*e ) * sind(eAN);

  sDIST = Math.sqrt( x*x + y*y );    // Solar distance

  v = atan2d( y, x );                // True anomaly

  sLON = v + w;                      // True solar longitude

  if ( sLON >= 360.0 )
    sLON -= 360.0;                   // Make it 0..360 degrees

} //=================== sunPos() =============================


var INV360 = 1.0 / 360.0;


//-------------------------------------------------------------
// Reduce angle to within 0..360 degrees
//-------------------------------------------------------------
function revolution( x )
{
  return (x - 360.0 * Math.floor( x * INV360 ));
}


//-------------------------------------------------------------
// Reduce angle to within -180..+180 degrees
//-------------------------------------------------------------
function rev180( x )
{
  return ( x - 360.0 * Math.floor( x * INV360 + 0.5 ) );
}


//-------------------------------------------------------------
// This function computes GMST0, the Greenwhich Mean Sidereal
// Time at 0h UT (i.e. the sidereal time at the Greenwhich
// meridian at 0h UT).
// GMST is then the sidereal time at Greenwich at any time of
// the day.  GMST0 is generalized as well, and is defined as:
//
//  GMST0 = GMST - UT
//
// This allows GMST0 to be computed at other times than 0h UT
// as well.  While this sounds somewhat contradictory, it is
// very practical:
// Instead of computing GMST like:
//
//  GMST = (GMST0) + UT * (366.2422/365.2422)
//
// where (GMST0) is the GMST last time UT was 0 hours, one simply
// computes:
//
//  GMST = GMST0 + UT
//
// where GMST0 is the GMST "at 0h UT" but at the current moment!
// Defined in this way, GMST0 will increase with about 4 min a
// day.  It also happens that GMST0 (in degrees, 1 hr = 15 degr)
// is equal to the Sun's mean longitude plus/minus 180 degrees!
// (if we neglect aberration, which amounts to 20 seconds of arc
// or 1.33 seconds of time)
//-------------------------------------------------------------
function GMST0( dayDiff )
{
  var const1 = 180.0 + 356.0470 + 282.9404;

//  var const2 = 0.9856002585 + 4.70935E-5;

  var const2 = 0.9856002585 + 4.70935/100000;  // for Opera


  return ( revolution( const1 + const2 * dayDiff ) );
      
} //=================== GMST0() =========================


//* * * * * * * *   SUNRISET SCRIPT - END -  * * * * * * * * *