1 |
{ RCSid $Id: sun2.cal,v 1.3 2019/10/19 00:01:33 greg Exp $ } |
2 |
{ |
3 |
Improved solar angle calculations based on "The Astronomical Almanac's |
4 |
Algorithm for Approximate Solar Position (1950-2050)" by Joseph J. |
5 |
Michalsky, 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 AD (1950-2050) |
17 |
|
18 |
Outputs: |
19 |
STIME = Solar time (hours past midnight) |
20 |
SAZI = Solar altitude (degrees above horizon) |
21 |
SALT = Solar azimuth (degrees west of south) |
22 |
} |
23 |
{ Defaults for San Francisco, CA } |
24 |
LAT = 37.8152145; |
25 |
LON = 122.04; |
26 |
SM = 120; |
27 |
DEG : PI/180; |
28 |
{ Constants and functions } |
29 |
normcir(a,c) : if(-a, normcir(a+c, c), if(c-a, a, normcir(a-c, c))); |
30 |
and(a,b) : if(a, b, a); |
31 |
leap(yr) : floor(yr/4 + .125) - yr/4 + .125; |
32 |
{ Julian date } |
33 |
delta = YEAR - 1949; |
34 |
subjd = DAY + select(MONTH,0,31,59,90,120,151,181,212,243,273,304,334) + |
35 |
if(and(MONTH-2.5,leap(YEAR)), 1, 0); |
36 |
UTIME = TIME + SM/15; |
37 |
jd = 2432916.5 + delta*365 + floor(delta/4) + subjd + UTIME/24; |
38 |
{ Ecliptic coordinates } |
39 |
n = jd - 2451545; |
40 |
L = normcir(280.460 + 0.9856474*n, 360); |
41 |
g = normcir(357.528 + 0.9856003*n, 360); |
42 |
l = normcir(L + 1.915*sin(g*DEG) + 0.020*sin(2*DEG*g), 360); |
43 |
ep = 23.439 - 4e-7*n; |
44 |
{ Reused values } |
45 |
sin_l = sin(l*DEG); |
46 |
sin_LAT = sin(LAT*DEG); |
47 |
{ Right ascension and declination angles } |
48 |
ra = normcir(atan2(sin_l*cos(ep*DEG), cos(l*DEG))/DEG, 360); |
49 |
sin_dec = sin(ep*DEG)*sin_l; |
50 |
cos_dec = sqrt(1 - sin_dec*sin_dec); |
51 |
{ Greenwich mean sidereal time } |
52 |
gmst = normcir(6.697375 + 0.0657098242*n + UTIME, 24); |
53 |
lmst = normcir(gmst - LON/15, 24); |
54 |
STIME = normcir(lmst - ra/15 + 12, 24); |
55 |
ha = 15*(STIME - 12); |
56 |
{ Final solar angles } |
57 |
sin_SALT = sin_dec*sin_LAT + cos_dec*cos(LAT*DEG)*cos(ha*DEG); |
58 |
SALT = asin(sin_SALT) / DEG; |
59 |
az = asin(-cos_dec*sin(ha*DEG)/cos(SALT*DEG)) / DEG; |
60 |
SAZI = normcir(if(sin_dec-sin_SALT*sin_LAT, az, 180-az), 360) - 180; |