| 1 |
greg |
1.1 |
{ RCSid $Id$ } |
| 2 |
|
|
{ |
| 3 |
|
|
Improved solar angle calculations based on "The Astronomical Almanac's |
| 4 |
|
|
Algorithm for Approximate Solar Position (1950-2050)" by Joseph J. Michalsky, |
| 5 |
|
|
published in 1988 in Solar Energy, Vol. 40, No. 3 |
| 6 |
|
|
|
| 7 |
|
|
Translated by G.Ward, October 2019 |
| 8 |
|
|
|
| 9 |
|
|
Inputs: |
| 10 |
|
|
LAT = Site latitude (degrees north of Equator) |
| 11 |
|
|
LON = Site longitude (degrees west of Greenwich) |
| 12 |
|
|
SM = Standard meridian (degrees west of Greenwich) |
| 13 |
|
|
TIME = Standard time (fractional hours 0-24) |
| 14 |
|
|
DAY = Day of the month |
| 15 |
|
|
MONTH = Month (January == 1) |
| 16 |
|
|
YEAR = Calendar year |
| 17 |
|
|
|
| 18 |
|
|
Outputs: |
| 19 |
|
|
SAZI = Solar altitude (degrees above horizon) |
| 20 |
|
|
SALT = Solar azimuth (degrees west of south) |
| 21 |
|
|
} |
| 22 |
|
|
DEG : PI/180; |
| 23 |
|
|
normcir(a,c) : if(-a, normcir(a+c, c), if(c-a, a, normcir(a-c, c))); |
| 24 |
|
|
and(a,b) : if(a, b, a); |
| 25 |
|
|
mod(n,d) : n - floor(n/d)*d; |
| 26 |
|
|
leap(yr) : .2 - mod(yr+.1, 4); |
| 27 |
|
|
|
| 28 |
|
|
LAT = 37.8152145; { Defaults for San Francisco, CA } |
| 29 |
|
|
LON = 122.04; |
| 30 |
|
|
SM = 120; |
| 31 |
|
|
{ Julian date } |
| 32 |
|
|
delta = YEAR - 1949; |
| 33 |
|
|
subjd = DAY + select(MONTH,0,31,59,90,120,151,181,212,243,273,304,334) + |
| 34 |
|
|
if(and(MONTH-2.5,leap(YEAR)), 1, 0); |
| 35 |
|
|
UTIME = TIME + SM/15; |
| 36 |
|
|
jd = 2432916.5 + delta*365 + floor(delta/4) + subjd + UTIME/24; |
| 37 |
|
|
{ Ecliptic coordinates } |
| 38 |
|
|
n = jd - 2451545; |
| 39 |
|
|
L = normcir(280.460 + 0.9856474*n, 360); |
| 40 |
|
|
g = normcir(357.528 + 0.9856003*n, 360); |
| 41 |
|
|
l = normcir(L + 1.915*sin(g*DEG) + 0.020*sin(2*DEG*g), 360); |
| 42 |
|
|
ep = 23.439 - 4e-7*n; |
| 43 |
|
|
{ Reused values } |
| 44 |
|
|
sin_l = sin(l*DEG); |
| 45 |
|
|
sin_LAT = sin(LAT*DEG); |
| 46 |
|
|
cos_dec = cos(dec*DEG); |
| 47 |
|
|
{ Right ascension and declination angles } |
| 48 |
|
|
ra = normcir(atan2(sin_l*cos(ep*DEG), cos(l*DEG))/DEG, 360); |
| 49 |
|
|
dec = asin(sin(ep*DEG)*sin_l) / DEG; |
| 50 |
|
|
{ Greenwich mean sidereal time } |
| 51 |
|
|
gmst = normcir(6.697375 + 0.0657098242*n + UTIME, 24); |
| 52 |
|
|
lmst = normcir(gmst - LON/15, 24); |
| 53 |
|
|
ha = 15*(normcir(lmst - ra/15 + 12, 24) - 12); |
| 54 |
|
|
{ Final solar angles } |
| 55 |
|
|
SALT = asin(sin(dec*DEG)*sin_LAT + |
| 56 |
|
|
cos_dec*cos(LAT*DEG)*cos(ha*DEG)) / DEG; |
| 57 |
|
|
elc = asin(sin(dec*DEG)/sin_LAT) / DEG; |
| 58 |
|
|
az = asin(-cos_dec*sin(ha*DEG)/cos(SALT*DEG)) / DEG; |
| 59 |
|
|
SAZI = normcir(if(SALT-elc, 180-az, az), 360) - 180; |
