Sunrise/Sunset calculations

See http://williams.best.vwh.net/sunrise_sunset_algorithm.htm

by medmunds

HTML

<form>
    <label>Lat: <input id='lat' value='37.761'></label>
    <label>Lon: <input id='lon' value='-122.402'></label>
    <label>Date: <input type='date' id='date'></label>
    <button>Calculate</button>
</form>
<dl id='output'></dl>

CSS

label, button {
    display: block;
}
label[type='datetime'] {
    width: 140px;
}

dt {
    float: left;
    width: 4em;
}

JavaScript

function solar_event(date, latitude, longitude, rising, zenith) {
    var year = date.getUTCFullYear(),
        month = date.getUTCMonth() + 1,
        day = date.getUTCDate();

    var floor = Math.floor,
        degtorad = function(deg) {
            return Math.PI * deg / 180;
        },
        radtodeg = function(rad) {
            return 180 * rad / Math.PI;
        },
        sin = function(deg) {
            return Math.sin(degtorad(deg));
        },
        cos = function(deg) {
            return Math.cos(degtorad(deg));
        },
        tan = function(deg) {
            return Math.tan(degtorad(deg));
        },
        asin = function(x) {
            return radtodeg(Math.asin(x));
        },
        acos = function(x) {
            return radtodeg(Math.acos(x));
        },
        atan = function(x) {
            return radtodeg(Math.atan(x));
        },
        modpos = function(x, m) {
            return ((x % m) + m) % m;
        };

    // 1. first calculate the day of the year
    var N1 = floor(275 * month / 9),
        N2 = floor((month + 9) / 12),
        N3 = (1 + floor((year - 4 * floor(year / 4) + 2) / 3)),
        N = N1 - (N2 * N3) + day - 30;

    // 2. convert the longitude to hour value and calculate an approximate time
    var lngHour = longitude / 15,
        t = N + (((rising ? 6 : 18) - lngHour) / 24);

    // 3. calculate the Sun's mean anomaly
    var M = (0.9856 * t) - 3.289;

    // 4. calculate the Sun's true longitude
    var L = M + (1.916 * sin(M)) + (0.020 * sin(2 * M)) + 282.634;
    L = modpos(L, 360); // NOTE: L potentially needs to be adjusted into the range [0,360) by adding/subtracting 360
    // 5a. calculate the Sun's right ascension
    var RA = atan(0.91764 * tan(L));
    RA = modpos(RA, 360); // NOTE: RA potentially needs to be adjusted into the range [0,360) by adding/subtracting 360
    // 5b. right ascension value needs to be in the same quadrant as L
    var Lquadrant = (floor(L / 90)) * 90,
       ...