1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
|
/*
* Parts of this file have been taken from the Redshift project:
* https://github.com/jonls/redshift/
*
* Copyright 2010 Jon Lund Steffensen
* Copyright 2023 Robert Tari
*
* This program is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation; version 3.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU Lesser General Public License for more details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*
* Authors:
* Jon Lund Steffensen <jonlst@gmail.com>
* Robert Tari <robert@tari.in>
*/
#include <math.h>
#include "solar.h"
#define RAD(x) ((x)*(M_PI/180))
#define DEG(x) ((x)*(180/M_PI))
/* Julian centuries since J2000.0 from Julian day */
static double
jcent_from_jd(double jd)
{
return (jd - 2451545.0) / 36525.0;
}
/* Angular elevation at the location for the given hour angle.
lat: Latitude of location in degrees
decl: Declination in radians
ha: Hour angle in radians
Return: Angular elevation in radians */
static double
elevation_from_hour_angle(double lat, double decl, double ha)
{
return asin(cos(ha)*cos(RAD(lat))*cos(decl) +
sin(RAD(lat))*sin(decl));
}
/* Geometric mean anomaly of the sun.
t: Julian centuries since J2000.0
Return: Geometric mean anomaly in radians. */
static double
sun_geom_mean_anomaly(double t)
{
return RAD(357.52911 + t*(35999.05029 - t*0.0001537));
}
/* Equation of center of the sun.
t: Julian centuries since J2000.0
Return: Center(?) in radians */
static double
sun_equation_of_center(double t)
{
/* Use the first three terms of the equation. */
double m = sun_geom_mean_anomaly(t);
double c = sin(m)*(1.914602 - t*(0.004817 + 0.000014*t)) +
sin(2*m)*(0.019993 - 0.000101*t) +
sin(3*m)*0.000289;
return RAD(c);
}
/* Geometric mean longitude of the sun.
t: Julian centuries since J2000.0
Return: Geometric mean logitude in radians. */
static double
sun_geom_mean_lon(double t)
{
/* FIXME returned value should always be positive */
return RAD(fmod(280.46646 + t*(36000.76983 + t*0.0003032), 360));
}
/* True longitude of the sun.
t: Julian centuries since J2000.0
Return: True longitude in radians */
static double
sun_true_lon(double t)
{
double l_0 = sun_geom_mean_lon(t);
double c = sun_equation_of_center(t);
return l_0 + c;
}
/* Apparent longitude of the sun. (Right ascension).
t: Julian centuries since J2000.0
Return: Apparent longitude in radians */
static double
sun_apparent_lon(double t)
{
double o = sun_true_lon(t);
return RAD(DEG(o) - 0.00569 - 0.00478*sin(RAD(125.04 - 1934.136*t)));
}
/* Mean obliquity of the ecliptic
t: Julian centuries since J2000.0
Return: Mean obliquity in radians */
static double
mean_ecliptic_obliquity(double t)
{
double sec = 21.448 - t*(46.815 + t*(0.00059 - t*0.001813));
return RAD(23.0 + (26.0 + (sec/60.0))/60.0);
}
/* Corrected obliquity of the ecliptic.
t: Julian centuries since J2000.0
Return: Currected obliquity in radians */
static double
obliquity_corr(double t)
{
double e_0 = mean_ecliptic_obliquity(t);
double omega = 125.04 - t*1934.136;
return RAD(DEG(e_0) + 0.00256*cos(RAD(omega)));
}
/* Declination of the sun.
t: Julian centuries since J2000.0
Return: Declination in radians */
static double
solar_declination(double t)
{
double e = obliquity_corr(t);
double lambda = sun_apparent_lon(t);
return asin(sin(e)*sin(lambda));
}
/* Eccentricity of earth orbit.
t: Julian centuries since J2000.0
Return: Eccentricity (unitless). */
static double
earth_orbit_eccentricity(double t)
{
return 0.016708634 - t*(0.000042037 + t*0.0000001267);
}
/* Difference between true solar time and mean solar time.
t: Julian centuries since J2000.0
Return: Difference in minutes */
static double
equation_of_time(double t)
{
double epsilon = obliquity_corr(t);
double l_0 = sun_geom_mean_lon(t);
double e = earth_orbit_eccentricity(t);
double m = sun_geom_mean_anomaly(t);
double y = pow(tan(epsilon/2.0), 2.0);
double eq_time = y*sin(2*l_0) - 2*e*sin(m) +
4*e*y*sin(m)*cos(2*l_0) -
0.5*y*y*sin(4*l_0) -
1.25*e*e*sin(2*m);
return 4*DEG(eq_time);
}
/* Julian day from Julian centuries since J2000.0 */
static double
jd_from_jcent(double t)
{
return 36525.0*t + 2451545.0;
}
/* Solar angular elevation at the given location and time.
t: Julian centuries since J2000.0
lat: Latitude of location
lon: Longitude of location
Return: Solar angular elevation in radians */
static double
solar_elevation_from_time(double t, double lat, double lon)
{
/* Minutes from midnight */
double jd = jd_from_jcent(t);
double offset = (jd - round(jd) - 0.5)*1440.0;
double eq_time = equation_of_time(t);
double ha = RAD((720 - offset - eq_time)/4 - lon);
double decl = solar_declination(t);
return elevation_from_hour_angle(lat, decl, ha);
}
/* Julian day from unix epoch */
static double
jd_from_epoch(double t)
{
return (t / 86400.0) + 2440587.5;
}
/* Solar angular elevation at the given location and time.
date: Seconds since unix epoch
lat: Latitude of location
lon: Longitude of location
Return: Solar angular elevation in degrees */
double
solar_elevation(double date, double lat, double lon)
{
double jd = jd_from_epoch(date);
return DEG(solar_elevation_from_time(jcent_from_jd(jd), lat, lon));
}
|
