/
sunlon2time.m
83 lines (74 loc) · 3.54 KB
/
sunlon2time.m
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
function [time eot] = sunlon2time(sunlon,ecc,lpe,tottime,obl)
% [sday eot] = sunlon2time(sunlon,ecc,lpe,tottime,obl)
%
% Given a particular eccentricity and longitude of perihelion, get time of tropical year
% associated with a particular geocentric solar longitude, i.e. by accounting for
% conservation of angular momentum during orbit (Kepler 2nd Law).
%
% Input
% =====
% sunlon = Keplerian geocentric solar longitude in degrees ('lambda', i.e. 'v' relative to NH spring equinox)
% 0 = NH Spring, 90 = NH Summer, 180 = NH Autumn, 270 = NH Winter
% Either 1 value (used as constant if other inputs are vector), or a vector of values.
% ecc = eccentricity, array
% lpe = heliocentric longitude of perihelion (a.k.a omega-bar), array (in radians)
% tottime = total time in the year, single value, any constant time unit you want. Use empty, [], for 365.24.
% obl = optional fifth input variable, for when calculating eot. Obliquity (in radians).
%
% Apologies that sunlon is in degrees and lpe in radians, it seemed like a good idea at the time.
% Input can be vectorised in various ways, but double check output.
%
% Output
% ======
% time = Time interval of tropical year (where interval 0 is northern spring equinox).
% Same dims as sunlon.
% eot = equation of time (minutes). Requires fifth input variable obl.
% Same dims as sunlon. Returns empty if obl not supplied.
%
% B.C. Lougheed, June 2020, Matlab 2019a
% Updated April 2023 to include eot.
%
% See following for background, as well as comments in the script:
% Meeus, J., (1998). Astronomical Algorithms, 2nd ed. Willmann-Bell, Inc., Richmond, Virginia. (specifically Chapter 30).
% https://dr-phill-edwards.eu/Science/EOT.html (for equation of time)
if isempty(tottime) == 1
tottime = 365.24;
end
% change lpe from heliocentric to geocentric
omega = lpe+pi; % add 180
omega(omega>=2*pi) = omega(omega>=2*pi) - 2*pi; % wrap to 360
% First, get day of anchor day (dz) relative to perihelion
vz = 2*pi - omega; % v of spring equinox relative to perihelion
vz(vz>2*pi) = vz(vz>2*pi) - 2*pi;
Ez = 2 * atan( tan(vz/2) .* sqrt((1-ecc)./(1+ecc)) ); % Meeus (1998) page 195, solve for E
Mz = Ez-ecc.*sin(Ez); % Meeus page 195, solve for M (Kepler equation). M is the circular orbit equivalent of v
Mz(Mz<0) = pi + (pi-Mz(Mz<0)*-1); % inbound to perihelion
dz = Mz/(2*pi) .* tottime;
% Second, get day of target day (dx) relative to perihelion
vx = vz + deg2rad(sunlon);
vx(vx>2*pi) = vx(vx>2*pi) - 2*pi;
Ex = 2 * atan( tan(vx/2) .* sqrt((1-ecc)./(1+ecc)) ); % Meeus (1998) page 195, solve for E
Mx = Ex-ecc.*sin(Ex); % Meeus page 195, solve for M (Kepler equation). M is the circular orbit equivalent of v
Mx(Mx<0) = pi + (pi-Mx(Mx<0)*-1); % inbound to perihelion (probably not necessary)
dx = Mx/(2*pi) .* tottime;
% Finally, get day of target day (dx) relative to day of anchor day (dz)
dx(dx<dz) = dx(dx<dz) + tottime; % for dz in next orbital period relative to perihelion, keep in same orbital period relative to NH spring equinox
time = dx - dz;
% eliminate rounding errors at zero
time(sunlon==0) = 0;
time(sunlon==360) = 0;
% calculate eot
if nargin > 4 || ~isempty(obl)
sunlon = deg2rad(sunlon);
% https://dr-phill-edwards.eu/Science/EOT.html (somebody who actually explains it clearly)
% eccentricity component
dtecc = rad2deg(Mx-vx) * 4; % four minutes per degree
% obliquity component
alpha = atan2(sin(sunlon) * cos(obl), cos(sunlon));
alpha(alpha<0) = alpha(alpha<0) + 2*pi;
dtobl = rad2deg((sunlon-alpha)) * 4;
eot = dtecc + dtobl;
else
eot = [];
end
end % end main function