Sun, azimuth and altitude

sun.js

/**
 * http://stjarnhimlen.se/comp/tutorial.html
 * http://www.planetcalc.ru/320/
 */

// day
var year = 2016;
var month = 4;
var day = 18;

// Moscow
var latitude = 55.755830;
var longitude = 37.617780;

var localTime = 19 + (50 / 60); // hour + (minutes / 60)
var timeZone = 3; // Moscow, +3:00
var GMT = localTime - timeZone;

// "day number"
var d = (function (Y, M, D) {
    return Math.round(367 * Y - (7 * (Y + ((M + 9) / 12))) / 4 + (275 * M) / 9 + D - 730530);
})(year, month, day);

var w = 282.9404 + 4.70935E-5 * d;   // longitude of perihelion (degrees)
var e = 0.016709 - 1.151E-9 * d;     // eccentricity
var M = 356.0470 + 0.9856002585 * d; // mean anomaly (degrees)

M = rev(M);

var oblecl = 23.4393 - 3.563E-7 * d; // obliquity of the ecliptic
var L = rev(w + M); // Sun's mean longitude

var E = M + (180 / Math.PI) * e * Math.sin(deg2rad(M)) * (1 + e * Math.cos(deg2rad(M)));

var x = Math.cos(deg2rad(E)) - e;
var y = Math.sin(deg2rad(E)) * Math.sqrt(1 - e * e);

var r = Math.sqrt(x * x + y * y); // Heliocentric distance: the planet's distance from the Sun
var v = rad2deg(Math.atan2(y, x));

var lon = rev(v + w);

x = r * Math.cos(deg2rad(lon));
y = r * Math.sin(deg2rad(lon));

var xequat = x;
var yequat = y * Math.cos(deg2rad(oblecl));
var zequat = y * Math.sin(deg2rad(oblecl));

var RA = rad2deg(atan2(yequat, xequat)) / 15;
var Decl = rad2deg(atan2(zequat, Math.sqrt(xequat * xequat + yequat * yequat)));

var GMST0 = L / 15 + 12;
var SIDTIME = GMST0 + GMT + longitude / 15;
var HA = (SIDTIME - RA) * 15;

x = Math.cos(deg2rad(HA)) * Math.cos(deg2rad(Decl));
y = Math.sin(deg2rad(HA)) * Math.cos(deg2rad(Decl));
var z = Math.sin(deg2rad(Decl));

var xhor = x * Math.sin(deg2rad(latitude)) - z * Math.cos(deg2rad(latitude));
var yhor = y;
var zhor = x * Math.cos(deg2rad(latitude)) + z * Math.sin(deg2rad(latitude));

var azimuth = rad2deg(atan2(yhor, xhor)) + 180;
var altitude = rad2deg(atan2(zhor, Math.sqrt(xhor * xhor + yhor * yhor)));

if (azimuth < 0) {
    azimuth += 360;
}

console.log("azimuth = " + azimuth);
console.log("altitude = " + altitude);

function rev(x) {
    var rv = x - Math.floor(x / 360) * 360;
    if (rv < 0) {
        rv += 360;
    }
    return rv;
}

function rad2deg(x) {
    return 180 * x / Math.PI;
}
function deg2rad(x) {
    return x * Math.PI / 180;
}
function atan2(y, x) {
    return Math.atan(y / x) - Math.PI * (x < 0);
}