forked from sandialabs/MATLAB_PV_LIB
-
Notifications
You must be signed in to change notification settings - Fork 0
/
pvl_ephemeris.m
199 lines (168 loc) · 8.46 KB
/
pvl_ephemeris.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
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
function [SunAz, SunEl, ApparentSunEl, SolarTime]= pvl_ephemeris(Time, Location, varargin)
% PVL_EPHEMERIS Calculates the position of the sun given time, location, and optionally pressure and temperature
%
% Syntax
% [SunAz, SunEl, ApparentSunEl, SolarTime]=pvl_ephemeris(Time, Location)
% [SunAz, SunEl, ApparentSunEl, SolarTime]=pvl_ephemeris(Time, Location, Pressure)
% [SunAz, SunEl, ApparentSunEl, SolarTime]=pvl_ephemeris(Time, Location, Pressure, Temperature)
% [SunAz, SunEl, ApparentSunEl, SolarTime]=pvl_ephemeris(Time, Location, 'temperature', Temperature)
%
% Description
% [SunAz, SunEl, ApparentSunEl, SolarTime]=pvl_ephemeris(Time, Location)
% Uses the given time and location structs to give sun positions with
% the pressure assumed to be 1 atm (101325 Pa) and the temperature
% assumed to be 12 C.
% [SunAz, SunEl, ApparentSunEl, SolarTime]=pvl_ephemeris(Time, Location, Pressure)
% Uses the given time and location structs with the given pressure to
% determine sun positions. The temperature is assumed to be 12C.
% Pressure must be given in Pascals (1atm = 101325 Pa). If site pressure
% is unknown but site altitude is known, a conversion function may be
% used.
% [SunAz, SunEl, ApparentSunEl, SolarTime]=pvl_ephemeris(Time, Location, Pressure, Temperature)
% Uses the given time and location structs with the given pressure and
% temperature to determine sun positions. Pressure must be given in
% Pascals, and temperature must be given in C.
% [SunAz, SunEl, ApparentSunEl, SolarTime]=pvl_ephemeris(Time, Location, 'temperature', Temperature)
% Uses the given time and location structs with the given temperature
% (in C) to determine sun positions. Default pressure is 101325 Pa.
%
% Input Parameters:
% Time is a struct with the following elements, note that all elements
% can be column vectors, but they must all be the same length. Time is
% entered as a struct which must include a value for the offset from
% UTC. Time is absolutely specified by the date, time, and the number
% of hours of offset from that date and time to UTC. For example, if
% you live in Boston, USA, and your data is timestamped in local standard
% time, your UTC offset should be -5.
%
% Time.year = The year in the gregorian calendar
% Time.month = the month of the year (January = 1 to December = 12)
% Time.day = the day of the month
% Time.hour = the hour of the day
% Time.minute = the minute of the hour
% Time.second = the second of the minute
% Time.UTCOffset = the UTC offset code, using the convention
% that a positive UTC offset is for time zones east of the prime meridian
% (e.g. EST = -5)
%
% Location is a struct with the following elements
% Location.latitude = vector or scalar latitude in decimal degrees (positive is
% northern hemisphere)
% Location.longitude = vector or scalar longitude in decimal degrees (positive is
% east of prime meridian)
% Location.altitude = an optional component of the Location struct, not
% used in the ephemeris code directly, but it may be used to calculate
% standard site pressure (see pvl_alt2pres function)
%
% Output Parameters:
% SunAz = Azimuth of the sun in decimal degrees from North. 0 = North to 270 = West
% SunEl = Actual elevation (not accounting for refraction)of the sun
% in decimal degrees, 0 = on horizon. The complement of the True Zenith
% Angle.
% ApparentSunEl = Apparent sun elevation accounting for atmospheric
% refraction. This is the complement of the Apparent Zenith Angle.
% SolarTime = solar time in decimal hours (solar noon is 12.00).
%
% References
% Grover Hughes' class and related class materials on Engineering
% Astronomy at Sandia National Laboratories, 1985.
%
% See also PVL_MAKETIMESTRUCT PVL_MAKELOCATIONSTRUCT PVL_ALT2PRES
% PVL_GETAOI PVL_SPA
p = inputParser;
p.addRequired('Time',@isstruct);
p.addRequired('Location',@isstruct);
p.addOptional('pressure',101325, @(x) all(isvector(x) & isnumeric(x) & x>=0));
p.addOptional('temperature',12, @(x) all(isvector(x) & isnumeric(x) & x>=-273.15));
p.parse(Time, Location, varargin{:});
Latitude = p.Results.Location.latitude(:);
% the inversion of longitude is due to the fact that this code was
% originally written for the convention that positive longitude were for
% locations west of the prime meridian. However, the correct convention (as
% of 2009) is to use negative longitudes for locations west of the prime
% meridian. Therefore, the user should input longitude values under the
% correct convention (e.g. Albuquerque is at -106 longitude), but it needs
% to be inverted for use in the code.
Longitude = -1 * p.Results.Location.longitude(:);
Year = p.Results.Time.year(:);
Month = p.Results.Time.month(:);
Day = p.Results.Time.day(:);
Hour = p.Results.Time.hour(:);
Minute = p.Results.Time.minute(:);
Second = p.Results.Time.second(:);
% the inversion of UTC offset is due to the fact that this code was
% originally written for the convention that positive offset values were
% positive for locations west of the prime merdian. However, the correct
% convention (as of 2009) is to use negative offset codes for locaitons
% west of the prime merdian. Therefore, the user should input offset values
% under the correct convention (e.g. EST = -5), but it needs to be inverted
% for use in the following code.
TZone = -1 * p.Results.Time.UTCOffset(:);
if isscalar(p.Results.pressure)
pressure = p.Results.pressure*ones(size(Year));
else
pressure = p.Results.pressure(:);
end
if isscalar(p.Results.temperature)
temperature = p.Results.temperature*ones(size(Year));
else
temperature = p.Results.temperature(:);
end
if ~(numel(pressure)==numel(temperature)) && (numel(pressure)==numel(Time.Year)) && ((numel(pressure)==numel(Latitude)) || isscalar(Latitude))
error(['An error has occurred in ',mfilename,'. One of the input'...
' variables is not of correct size.']);
end
DayOfYear = pvl_date2doy(Year, Month, Day);
DecHours = Hour + Minute./60 + Second./3600;
RadtoDeg = 180 / pi;
DegtoRad = pi / 180;
Abber = 20/3600;
LatR = Latitude * DegtoRad;
UnivDate = DayOfYear + floor((DecHours + TZone)/24);
UnivHr = mod((DecHours + TZone), 24);
Yr = Year-1900;
YrBegin = 365 * Yr + floor((Yr-1)/4)-0.5;
Ezero = YrBegin + UnivDate;
T = Ezero / 36525;
GMST0 = 6/24 +38/1440 + (45.836 + 8640184.542 * T + 0.0929 * T.^2)/86400;
GMST0 = 360 * (GMST0 - floor(GMST0));
GMSTi = mod(GMST0 + 360*(1.0027379093 * UnivHr / 24),360);
LocAST = mod((360 + GMSTi - Longitude), 360);
EpochDate = Ezero + UnivHr / 24;
T1 = EpochDate / 36525;
ObliquityR = DegtoRad * (23.452294 - 0.0130125 * T1 - 0.00000164 * T1.^2 ...
+ 0.000000503 * T1.^3);
MlPerigee = 281.22083 + 0.0000470684 * EpochDate + 0.000453 * T1 .^ 2 + ...
0.000003 * T1 .^ 3;
MeanAnom = mod((358.47583 + 0.985600267 * EpochDate - 0.00015 * T1 .^ 2 - ...
0.000003 * T1 .^ 3), 360);
Eccen = 0.01675104 - 0.0000418 * T1 - 0.000000126 * T1 .^ 2;
EccenAnom = MeanAnom;
E=0;
while max(abs(EccenAnom - E)) > 0.0001;
E = EccenAnom;
EccenAnom = MeanAnom + RadtoDeg .* Eccen .* sin(DegtoRad .* E);
end
TrueAnom = 2 * mod(RadtoDeg * atan2(((1 + Eccen) ./ (1 - Eccen)).^ 0.5 .* tan(DegtoRad * EccenAnom / 2), 1), 360) ;
EcLon = mod(MlPerigee + TrueAnom, 360) - Abber ;
EcLonR = DegtoRad * EcLon;
DecR = asin(sin(ObliquityR) .* sin(EcLonR));
Dec = RadtoDeg * DecR;
RtAscen = RadtoDeg * atan2(cos(ObliquityR).*(sin(EcLonR)),cos(EcLonR));
HrAngle = LocAST - RtAscen ;
HrAngleR = DegtoRad .* HrAngle ;
HrAngle = HrAngle - (360 .* sign(HrAngle) .* (abs(HrAngle) > 180));
SunAz = RadtoDeg .* atan2(-1 * sin(HrAngleR), cos(LatR) .* tan(DecR) - sin(LatR) .* cos(HrAngleR));
SunAz = SunAz + (SunAz < 0) * 360; %shift from range of [-180,180] to [0,360]
SunEl = asind(cos(LatR) .* cos(DecR) .* cos(HrAngleR) + sin(LatR) .* sin(DecR));
SolarTime = (180 + HrAngle) / 15;
% Calculate the refraction of the sun until the actual center of the sun is
% 1 degree below the horizon.
TanEl = tan(DegtoRad * SunEl);
Refract = zeros(length(SunEl),1) + ...
(and(SunEl > 5, SunEl <= 85) .* (58.1 ./ TanEl - 0.07 ./ (TanEl.^3) + .000086 ./ (TanEl.^5))) + ...
(and(SunEl > -0.575, SunEl <=5) .* (SunEl .* (-518.2 + SunEl .* (103.4 + SunEl .* (-12.79 + SunEl .* 0.711))) +1735)) + ...
(and(SunEl > -1 ,SunEl <= -0.575) .* ((-20.774 ./ TanEl)));
Refract = Refract .* (283 ./ (273 + temperature)) .* pressure ./ 101325 ./ 3600;
% Generate apparent sun elevation including refraction
ApparentSunEl = SunEl + Refract;