{ RCSid $Id: sun2.cal,v 1.2 2019/10/18 23:42:52 greg Exp $ }
{
	Improved solar angle calculations based on "The Astronomical Almanac's
	Algorithm for Approximate Solar Position (1950-2050)" by Joseph J.
	Michalsky, published in 1988 in Solar Energy, Vol. 40, No. 3

	Translated by G.Ward, October 2019

	Inputs:
		LAT =		Site latitude (degrees north of Equator)
		LON =		Site longitude (degrees west of Greenwich)
		SM =		Standard meridian (degrees west of Greenwich)
		TIME =		Standard time (fractional hours 0-24)
		DAY =		Day of the month
		MONTH =		Month (January == 1)
		YEAR =		Calendar year AD (1950-2050)

	Outputs:
		SAZI =		Solar altitude (degrees above horizon)
		SALT =		Solar azimuth (degrees west of south)
}
				{ Defaults for San Francisco, CA }
LAT = 37.8152145;
LON = 122.04;
SM = 120;
DEG : PI/180;
				{ Constants and functions }
normcir(a,c) : if(-a, normcir(a+c, c), if(c-a, a, normcir(a-c, c)));
and(a,b) : if(a, b, a);
leap(yr) : floor(yr/4 + .125) - yr/4 + .125;
				{ Julian date }
delta = YEAR - 1949;
subjd = DAY + select(MONTH,0,31,59,90,120,151,181,212,243,273,304,334) +
		if(and(MONTH-2.5,leap(YEAR)), 1, 0);
UTIME = TIME + SM/15;
jd = 2432916.5 + delta*365 + floor(delta/4) + subjd + UTIME/24;
				{ Ecliptic coordinates }
n = jd - 2451545;
L = normcir(280.460 + 0.9856474*n, 360);
g = normcir(357.528 + 0.9856003*n, 360);
l = normcir(L + 1.915*sin(g*DEG) + 0.020*sin(2*DEG*g), 360);
ep = 23.439 - 4e-7*n;
				{ Reused values }
sin_l = sin(l*DEG);
sin_LAT = sin(LAT*DEG);
				{ Right ascension and declination angles }
ra = normcir(atan2(sin_l*cos(ep*DEG), cos(l*DEG))/DEG, 360);
sin_dec = sin(ep*DEG)*sin_l;
cos_dec = sqrt(1 - sin_dec*sin_dec);
dec = asin(sin_dec) / DEG;
				{ Greenwich mean sidereal time }
gmst = normcir(6.697375 + 0.0657098242*n + UTIME, 24);
lmst = normcir(gmst - LON/15, 24);
ha = 15*(normcir(lmst - ra/15 + 12, 24) - 12);
				{ Final solar angles }
sin_SALT = sin_dec*sin_LAT + cos_dec*cos(LAT*DEG)*cos(ha*DEG);
SALT = asin(sin_SALT) / DEG;
az = asin(-cos_dec*sin(ha*DEG)/cos(SALT*DEG)) / DEG;
SAZI = normcir(if(sin_dec-sin_SALT*sin_LAT, az, 180-az), 360) - 180;