Skip to content

Commit 4bd7bfe

Browse files
DanRileyDanRiley
committed
Update to the default delta T calculation within the SPA code.
1 parent c0c46f4 commit 4bd7bfe

File tree

1 file changed

+116
-5
lines changed

1 file changed

+116
-5
lines changed

‎pvl_spa.m

Lines changed: 116 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -60,8 +60,8 @@
6060
%
6161
% delta_t = The current difference, in seconds, between international atomic time
6262
% (TAI) and UT1. This value is obtained by observation (see [2]). If
63-
% omitted, the default value is 66.3 + 0.6175*(Year-2012). This
64-
% default prediction can be changed as necessary.
63+
% omitted, the default value is provided according to [3]. This default
64+
% can be changed if necessary, but the default is preferred.
6565
%
6666
% Output Parameters:
6767
% SunAz = Azimuth of the sun in decimal degrees from North. 0 = North to 270 = West
@@ -76,8 +76,10 @@
7676
% Radiaiton Applications. NREL Report No. TP-560-34302. Revised
7777
% January 2008. http://www.nrel.gov/docs/fy08osti/34302.pdf as viewed
7878
% 2012-05-07
79-
% [2] DeltaT projections and current values from US Naval Observatory,
80-
% http://asa.usno.navy.mil/SecK/DeltaT.html as viewed 2012-05-07
79+
% [2] Delta T and Universal Time, Fred Spenak, retrieved from
80+
% http://eclipse.gsfc.nasa.gov/SEhelp/deltaT.html on 2016-06-22
81+
% [3] Polynomial Expressions for Delta T, retrieved from
82+
% http://eclipse.gsfc.nasa.gov/SEhelp/deltatpoly2004.html on 2016-06-22
8183
%
8284
% See also PVL_MAKETIMESTRUCT PVL_MAKELOCATIONSTRUCT PVL_ALT2PRES
8385
% PVL_GETAOI PVL_EPHEMERIS
@@ -87,7 +89,7 @@
8789
p.addRequired('Location',@isstruct);
8890
p.addOptional('pressure',101325, @(x) all(isvector(x) & isnumeric(x) & x>=0));
8991
p.addOptional('temperature',12, @(x) all(isvector(x) & isnumeric(x) & x>=-273.15));
90-
p.addOptional('delta_t', 66.3 + .6175*(Time.year-2012), @(x) all(isvector(x) & isnumeric(x) & x>-8000 & x<8000));
92+
p.addOptional('delta_t', calculate_deltaT(Time.year, Time.month), @(x) all(isvector(x) & isnumeric(x) & x>-8000 & x<8000));
9193
p.parse(Time, Location, varargin{:});
9294

9395
% Read in optional values from the struct p
@@ -441,6 +443,115 @@
441443

442444

443445
end
446+
447+
function [DeltaT]=calculate_deltaT(year, month)
448+
449+
% The equations aren't meant to be used outside of the range of years
450+
% [-1999, 3000]
451+
if any(year < -1999 | year > 3000)
452+
warning(['Delta T is unknown for years before -1999 and after 3000.' ...
453+
' DeltaT values will be calculated, but the calculations ' ...
454+
'are not intended to be used for these years.'])
455+
end
456+
457+
y = year + (month + 0.5) ./ 12;
458+
459+
DeltaT = zeros(size(year)); % create a starting vector
460+
461+
% Implement the piecewise polynomial expressions for delta T as provided
462+
% by http://eclipse.gsfc.nasa.gov/SEhelp/deltatpoly2004.html on
463+
% June 22, 2016.
464+
465+
% It is curious that in the years < -500 and >= 2150 the equation uses
466+
% the year itself rather than the decimal year as used by all other
467+
% time periods, but this is what the NASA eclipse site says...
468+
DeltaT(year < -500) = ...
469+
-20 ...
470+
+ 32 .* (((year(year < -500) - 1820) ./ 100).^2);
471+
DeltaT(year >= -500 & year < 500) = ...
472+
10583.6 ...
473+
- 1014.41 .* (y(year >= -500 & year < 500) ./ 100) ...
474+
+ 33.78311 .* ((y(year >= -500 & year < 500) ./ 100)).^2 ...
475+
- 5.952053 .* ((y(year >= -500 & year < 500) ./ 100)).^3 ...
476+
- 0.1798452 .* ((y(year >= -500 & year < 500) ./ 100)).^4 ...
477+
+ 0.022174192 .* ((y(year >= -500 & year < 500) ./ 100)).^5 ...
478+
+ 0.0090316521 .* ((y(year >= -500 & year < 500) ./ 100)).^6;
479+
DeltaT(year >= 500 & year < 1600) = ...
480+
1574.2 ...
481+
- 556.01 .* ((y(year >= 500 & year < 1600) - 1000) ./ 100) ...
482+
+ 71.23472 .* ((y(year >= 500 & year < 1600) - 1000) ./ 100).^2 ...
483+
+ 0.319781 .* ((y(year >= 500 & year < 1600) - 1000) ./ 100).^3 ...
484+
- 0.8503463 .* ((y(year >= 500 & year < 1600) - 1000) ./ 100).^4 ...
485+
- 0.005050998 .* ((y(year >= 500 & year < 1600) - 1000) ./ 100).^5 ...
486+
+ 0.0083572073 .* ((y(year >= 500 & year < 1600) - 1000) ./ 100).^6;
487+
DeltaT(year >= 1600 & year < 1700) = ...
488+
120 ...
489+
- 0.9808 .* (y(year >= 1600 & year < 1700) - 1600) ...
490+
- 0.01532 .* (y(year >= 1600 & year < 1700) - 1600).^2 ...
491+
+ ((y(year >= 1600 & year < 1700) - 1600).^3) ./ 7129;
492+
DeltaT(year >= 1700 & year < 1800) = ...
493+
8.83 ...
494+
+ 0.1603 .* (y(year >= 1700 & year < 1800) - 1700) ...
495+
- 0.0059285 .* (y(year >= 1700 & year < 1800) - 1700).^2 ...
496+
+ 0.00013336 .* (y(year >= 1700 & year < 1800) - 1700).^3 ...
497+
- ((y(year >= 1700 & year < 1800) - 1700).^4) ./ 1174000;
498+
DeltaT(year >= 1800 & year < 1860) = ...
499+
13.72 ...
500+
- 0.332447 .* (y(year >= 1800 & year < 1860) - 1800) ...
501+
+ 0.0068612 .* (y(year >= 1800 & year < 1860) - 1800).^2 ...
502+
+ 0.0041116 .* (y(year >= 1800 & year < 1860) - 1800).^3 ...
503+
- 0.00037436 .* (y(year >= 1800 & year < 1860) - 1800).^4 ...
504+
+ 0.0000121272 .* (y(year >= 1800 & year < 1860) - 1800).^5 ...
505+
- 0.0000001699 .* (y(year >= 1800 & year < 1860) - 1800).^6 ...
506+
+ 0.000000000875 .* (y(year >= 1800 & year < 1860) - 1800).^7;
507+
DeltaT(year >= 1860 & year < 1900) = ...
508+
7.62 ...
509+
+ 0.5737 .* (y(year >= 1860 & year < 1900) - 1860) ...
510+
- 0.251754 .* (y(year >= 1860 & year < 1900) - 1860).^2 ...
511+
+ 0.01680668 .* (y(year >= 1860 & year < 1900) - 1860).^3 ...
512+
- 0.0004473624 .* (y(year >= 1860 & year < 1900) - 1860).^4 ...
513+
+ ((y(year >= 1860 & year < 1900) - 1860).^5) ./ 233174;
514+
DeltaT(year >= 1900 & year < 1920) = ...
515+
-2.79 ...
516+
+ 1.494119 .* (y(year >= 1900 & year < 1920) - 1900) ...
517+
- 0.0598939 .* (y(year >= 1900 & year < 1920) - 1900).^2 ...
518+
+ 0.0061966 .* (y(year >= 1900 & year < 1920) - 1900).^3 ...
519+
- 0.000197 .* (y(year >= 1900 & year < 1920) - 1900).^4;
520+
DeltaT(year >= 1920 & year < 1941) = ...
521+
21.20 ...
522+
+ 0.84493 .* (y(year >= 1920 & year < 1941) - 1920) ...
523+
- 0.076100 .* (y(year >= 1920 & year < 1941) - 1920).^2 ...
524+
+ 0.0020936 .* (y(year >= 1920 & year < 1941) - 1920).^3;
525+
DeltaT(year >= 1941 & year < 1961) = ...
526+
29.07 ...
527+
+ 0.407 .* (y(year >= 1941 & year < 1961) - 1950) ...
528+
- (y(year >= 1941 & year < 1961) - 1950).^2 ./ 233 ...
529+
+ (y(year >= 1941 & year < 1961) - 1950).^3 ./ 2547;
530+
DeltaT(year >= 1961 & year < 1986) = ...
531+
45.45 ...
532+
+ 1.067 .* (y(year >= 1961 & year < 1986) - 1975) ...
533+
- (y(year >= 1961 & year < 1986) - 1975).^2 ./ 260 ...
534+
- (y(year >= 1961 & year < 1986) - 1975).^3 ./ 718;
535+
DeltaT(year >= 1986 & year < 2005) = ...
536+
63.86 ...
537+
+ 0.3345 .* (y(year >= 1986 & year < 2005) - 2000) ...
538+
- 0.060374 .* (y(year >= 1986 & year < 2005) - 2000).^2 ...
539+
+ 0.0017275 .* (y(year >= 1986 & year < 2005) - 2000).^3 ...
540+
+ 0.000651814 .* (y(year >= 1986 & year < 2005) - 2000).^4 ...
541+
+ 0.00002373599 .* (y(year >= 1986 & year < 2005) - 2000).^5;
542+
DeltaT(year >= 2005 & year < 2050) = ...
543+
62.92 ...
544+
+ 0.32217 .* (y(year >= 2005 & year < 2050) - 2000) ...
545+
+ 0.005589 .* (y(year >= 2005 & year < 2050) - 2000).^2;
546+
DeltaT(year >= 2050 & year < 2150) = ...
547+
-20 ...
548+
+ 32 .* ((y(year >= 2050 & year < 2150) - 1820) ./ 100).^2 ...
549+
- 0.5628 .* (2150 - y(year >= 2050 & year < 2150));
550+
DeltaT(year >= 2150) = ... % Same equation as <-500, provided separately for clarity and expansion
551+
-20 ...
552+
+ 32 .* (((year(year >= 2150) - 1820) ./ 100).^2);
553+
end
554+
444555
%%
445556
% This function just loads up the tables necessary
446557
function [L0table, L1table, L2table, L3table, L4table, L5table,...

0 commit comments

Comments
 (0)