From c0c46f450d6159d5dc5d8cdfb7d2b7fe43949d93 Mon Sep 17 00:00:00 2001 From: Riley Date: Tue, 14 Jun 2016 09:44:23 -0600 Subject: [PATCH 1/3] Update to pvl_getISDdata to use Haversine formula and include input parsing --- pvl_getISDdata.m | 61 ++++++++++++++++------- pvl_readtmy3.m | 125 ++++++++++++++++++++++++++++++++++++++++++++++- 2 files changed, 166 insertions(+), 20 deletions(-) diff --git a/pvl_getISDdata.m b/pvl_getISDdata.m index d622d37..2764f8a 100644 --- a/pvl_getISDdata.m +++ b/pvl_getISDdata.m @@ -1,4 +1,4 @@ -function outfilen = pvl_getISDdata(lat,lon,yr,archive) +function outfilen = pvl_getISDdata(latitude,longitude,year,archive) % % getISDdata: fetch data for one year of integrated surface data for location % @@ -6,14 +6,16 @@ % data = getISDdata(latitude,longitude,year,archive) % % Inputs: -% latitude - latitude of location of interest -% longitude - longitude of location of interest +% latitude - latitude of location of interest in decimal degrees +% longitude - longitude of location of interest in decimal degrees % year - desired year % archive - string specifying path to local ISD archive. Current % directory is used if archive is not specified. % % Output: -% outfilen - name for the returned data file +% outfilen - name for the returned data file in the form "USAF-WBAN-yyyy" +% where USAF is the USAF ID of the station and WBAN is the WBAN ID of +% the system, as provided by the ISD. % % Notes: % This function copies up to two files from NOAA's ftp site to the @@ -40,18 +42,30 @@ mod = vers(pos+6); if versyr<2013 || ((versyr==2013) && mod=='a') display(['Matlab ' vers ' detected']) - display('Function pvl_getISDdata requires R2013b or later') - return + error('Function pvl_getISDdata requires R2013b or later') end - -%% Get index data -% Change pwd to a local directory if you want to use an archive -if nargin < 4 +%% Check the input data +p = inputParser; +p.addRequired('latitude',@(x) all(isscalar(x) & isnumeric(x) & x<=90 & x>=-90)); +p.addRequired('longitude', @(x) all(isscalar(x) & isnumeric(x) & x<=180 & x>=-180)); +p.addRequired('year',@(x) all(isscalar(x) & isnumeric(x) & (x >= 1901))); +p.addOptional('archive', @(x) ischar(x)); +p.parse(latitude, longitude, year, varargin{:}); + +defaultchecker = {'archive'}; +if ~any(strcmp(defaultchecker,p.UsingDefaults)) archive = pwd; end +lat = p.Results.latitude; +lon = p.Results.longitude; +yr = p.Results.year; + +%% Get index data +% Change pwd to a local directory if you want to use an archive + isdfile = 'isd-history.csv'; -fname = [archive '\' isdfile]; +fname = [archive filesep isdfile]; if ~exist(fname,'file') % Get index file ftphandle = ftp('ftp.ncdc.noaa.gov'); @@ -69,17 +83,23 @@ index.Properties.VariableNames{'ELEV_M_'} = 'ELEV'; index.USAF = char(index.USAF); index.WBAN = char(index.WBAN); -index.LAT(strcmp('',index.LAT)) = {'99999'}; -index.LON(strcmp('',index.LON)) = {'99999'}; -index.ELEV(strcmp('',index.ELEV)) = {'99999'}; +index.LAT(strcmp('',index.LAT)) = {'NaN'}; +index.LON(strcmp('',index.LON)) = {'NaN'}; +index.ELEV(strcmp('',index.ELEV)) = {'NaN'}; index.LAT = str2num(char(index.LAT)); %#ok<*ST2NM> index.LON = str2num(char(index.LON)); index.ELEV = str2num(char(index.ELEV)); index.BEGIN = datenum(char(index.BEGIN),'yyyymmdd'); index.END = datenum(char(index.END),'yyyymmdd'); -%% Find row for closest location (might iterate through several locations) -dist = hypot(index.LAT-lat,index.LON-lon); +%% Find row for closest location (might iterate through several locations) using haversine formula +% calculate the argument for the haversine formula, make sure the argument +% is between -1 and 1 +arg = sqrt((sind((index.LAT - lat) ./ 2)).^2 + cosd(index.LAT).*cosd(lat).*(sind((index.LON-lon)./2)).^2); +arg(arg < -1) = -1; +arg(arg > 1) = 1; +% compute the distance in meters between the provided lat/lon and the each lat/lon from the isd. +dist = 2 .* 6371008.8 .* asin(arg); % mean earth radius is 6371.0088 km [~,idx] = sort(dist); % Find first row that contains year (might iterate through several years) @@ -99,10 +119,15 @@ % % get file from ftp server if not in archive -fname = [archive '\' outfilen]; +fname = [archive filesep outfilen]; if ~exist(fname,'file') addr = sprintf('ftp://ftp.ncdc.noaa.gov/pub/data/noaa/%4d/',yr); - gunzip([addr '/' outfilen '.gz'],archive); + try + gunzip([addr '/' outfilen '.gz'],archive); + catch + error(['Unable to reach the ftp site for the integrated surface database.' ... + ' The data file cannot be accessed.']); + end %% Alt approach if gunzip doesn't seem to be working % % cd(ftphandle,num2str(yr)); % % gunzip(mget(ftphandle,[fname '.gz'])); diff --git a/pvl_readtmy3.m b/pvl_readtmy3.m index cc3d595..3c72ff5 100644 --- a/pvl_readtmy3.m +++ b/pvl_readtmy3.m @@ -15,6 +15,9 @@ % If a FileName is not provided, the user will be prompted to browse to % an appropriate TMY3 file. % +% TMY3 format was changed slightly on January 19, 2015 ([2]) to include +% Present weather variables along with changes to other variables. +% % Input % FileName - an optional argument which allows the user to select which % TMY3 format file should be read. A file path may also be necessary if @@ -106,12 +109,23 @@ % TMYData.Lprecipquantity - The period of accumulation for the liquid precipitation depth field, hour % TMYData.LprecipSource - See [1], Table 1-5, 8760x1 cell array of strings % TMYData.LprecipUncertainty - See [1], Table 1-6 +% TMYData.PresWth - Present weather code, see [2]. Only available in +% TMY3 files created after the January 2015 TMY3 update. +% TMYData.PresWthSource - Present weather code source, see [2]. Only +% available in TMY3 files created after the Jan. 2015 TMY3 update. +% TMYData.PresWthUncertainty - Present weather code uncertainty, see +% [2]. Only available in TMY3 files after the Jan. 2015 update. % % Reference % [1] Wilcox, S and Marion, W., 2008. Users Manual for TMY3 Data Sets, % NREL/TP-581-43156, National Renewable Energy Laboratory. Available at % . % +% [2] National Solar Radiation Data Base, 1991-2005 Update: Typical +% Meteorological Year 3. +% http://rredc.nrel.gov/solar/old_data/nsrdb/1991-2005/tmy3/ retrieved June +% 6, 2016. +% % % See also % DATEVEC PVL_MAKELOCATIONSTRUCT PVL_MAKETIMESTRUCT PVL_READTMY2 @@ -136,7 +150,8 @@ Header2 = textscan(FileID, '%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s',1,'Delimiter',','); DataLines = textscan(FileID, '%s%s%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%s%f%f%s%f%f%s%f%f%s%f%f%s%f%f%s%f%f%s%f%f%s%f%f%s%f%f%s%f%f%s%f%f%s%f%f%s%f%f%f%s%f',8760,'Delimiter',','); ST = fclose(FileID); -%% Read in the data from the struct created from the textscan + +% Parse out the header information since it is common to both formats TMYData.SiteID = Header1{1}; TMYData.StationName = Header1{2}; TMYData.StationState = Header1{3}; @@ -144,7 +159,109 @@ TMYData.SiteLatitude = Header1{5}; TMYData.SiteLongitude = Header1{6}; TMYData.SiteElevation = Header1{7}; +%% Determine if the file is from pre-January 19, 2015 or post-20150119 +if isempty(strfind(char(DataLines{1,1}{1,1}), 'PresWth')) + % The string 'PresWth' was not found in DataLines{1,1}{1,1}, so we have a + % pre-2015 TMY3 file + TMYData = ParsePre2015TMY(DataLines); +else + % The string 'PresWth' was found in DataLines{1,1}{1,1}, so we have a + % post-2015 TMY3 file. We'll need to re-import the file to accomodate + % the additional columns. + FileID = fopen(FilePathAndNameAndExt); + Header1 = textscan(FileID, '%f%q%q%f%f%f%f', 1,'Delimiter',','); + Header2 = textscan(FileID, '%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s',1,'Delimiter',','); + DataLines = textscan(FileID, '%s%s%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%s%f%f%s%f%f%s%f%f%s%f%f%s%f%f%s%f%f%s%f%f%s%f%f%s%f%f%s%f%f%s%f%f%s%f%f%s%f%f%f%s%f%f%s%f',8760,'Delimiter',','); + ST = fclose(FileID); + + TMYData = ParsePost2015TMY(DataLines); +end +end +function TMYData = ParsePre2015TMY(DataLines) + +%% Read in the data from the struct created from the textscan for TMY3 +% This fuction will read pre-2015 TMY3 files and sort them into column +% vectors. A different format is used for TMY3 files after the January 2015 +% update, and thus requires a different function. + +TMYData.DateString = DataLines{1}; +TMYData.TimeString = DataLines{2}; +TMYData.ETR = DataLines{3}; +TMYData.ETRN = DataLines{4}; +TMYData.GHI = DataLines{5}; +TMYData.GHISource = DataLines{6}; +TMYData.GHIUncertainty = DataLines{7}; +TMYData.DNI = DataLines{8}; +TMYData.DNISource = DataLines{9}; +TMYData.DNIUncertainty = DataLines{10}; +TMYData.DHI = DataLines{11}; +TMYData.DHISource = DataLines{12}; +TMYData.DHIUncertainty = DataLines{13}; +TMYData.GHillum = DataLines{14}; +TMYData.GHillumSource = DataLines{15}; +TMYData.GHillumUncertainty = DataLines{16}; +TMYData.DNillum = DataLines{17}; +TMYData.DNillumSource = DataLines{18}; +TMYData.DNillumUncertainty = DataLines{19}; +TMYData.DHillum = DataLines{20}; +TMYData.DHillumSource = DataLines{21}; +TMYData.DHillumUncertainty = DataLines{22}; +TMYData.Zenithlum = DataLines{23}; +TMYData.ZenithlumSource = DataLines{24}; +TMYData.ZenithlumUncertainty = DataLines{25}; +TMYData.TotCld = DataLines{26}; +TMYData.TotCldSource = DataLines{27}; +TMYData.TotCldUnertainty = DataLines{28}; +TMYData.OpqCld = DataLines{29}; +TMYData.OpqCldSource = DataLines{30}; +TMYData.OpqCldUncertainty = DataLines{31}; +TMYData.DryBulb = DataLines{32}; +TMYData.DryBulbSource = DataLines{33}; +TMYData.DryBulbUncertainty = DataLines{34}; +TMYData.DewPoint = DataLines{35}; +TMYData.DewPointSource = DataLines{36}; +TMYData.DewPointUncertainty = DataLines{37}; +TMYData.RHum = DataLines{38}; +TMYData.RHumSource = DataLines{39}; +TMYData.RHumUncertainty = DataLines{40}; +TMYData.Pressure = DataLines{41}; +TMYData.PressureSource = DataLines{42}; +TMYData.PressureUncertainty = DataLines{43}; +TMYData.Wdir = DataLines{44}; +TMYData.WdirSource = DataLines{45}; +TMYData.WdirUncertainty = DataLines{46}; +TMYData.Wspd = DataLines{47}; +TMYData.WspdSource = DataLines{48}; +TMYData.WspdUncertainty = DataLines{49}; +TMYData.Hvis = DataLines{50}; +TMYData.HvisSource = DataLines{51}; +TMYData.HvisUncertainty = DataLines{52}; +TMYData.CeilHgt = DataLines{53}; +TMYData.CeilHgtSource = DataLines{54}; +TMYData.CeilHgtUncertainty = DataLines{55}; +TMYData.Pwat = DataLines{56}; +TMYData.PwatSource = DataLines{57}; +TMYData.PwatUncertainty = DataLines{58}; +TMYData.AOD = DataLines{59}; +TMYData.AODSource = DataLines{60}; +TMYData.AODUncertainty = DataLines{61}; +TMYData.Alb = DataLines{62}; +TMYData.AlbSource = DataLines{63}; +TMYData.AlbUncertainty = DataLines{64}; +TMYData.Lprecipdepth = DataLines{65}; +TMYData.Lprecipquantity = DataLines{66}; +TMYData.LprecipSource = DataLines{67}; +TMYData.LprecipUncertainty = DataLines{68}; +%% Create a MATLAB datenum from the string date and time +TMYData.DateNumber = datenum(strcat(TMYData.DateString,'-',TMYData.TimeString),'mm/dd/yyyy-HH:MM'); %MATLAB datenum format +end + + + +function TMYData = ParsePost2015TMY(DataLines) + +%% Read in the data from the struct created from the textscan TMYData.DateString = DataLines{1}; TMYData.TimeString = DataLines{2}; TMYData.ETR = DataLines{3}; @@ -213,5 +330,9 @@ TMYData.Lprecipquantity = DataLines{66}; TMYData.LprecipSource = DataLines{67}; TMYData.LprecipUncertainty = DataLines{68}; +TMYData.PresWth = DataLines{69}; +TMYData.PresWthSource = DataLines{70}; +TMYData.PresWthUncertainty = DataLines{71}; %% Create a MATLAB datenum from the string date and time -TMYData.DateNumber = datenum(strcat(TMYData.DateString,'-',TMYData.TimeString),'mm/dd/yyyy-HH:MM'); %MATLAB datenum format \ No newline at end of file +TMYData.DateNumber = datenum(strcat(TMYData.DateString,'-',TMYData.TimeString),'mm/dd/yyyy-HH:MM'); %MATLAB datenum format +end \ No newline at end of file From 4bd7bfe78a9eef5b648a59ba2423c1922a015ac3 Mon Sep 17 00:00:00 2001 From: Riley Date: Thu, 23 Jun 2016 10:29:27 -0600 Subject: [PATCH 2/3] Update to the default delta T calculation within the SPA code. --- pvl_spa.m | 121 +++++++++++++++++++++++++++++++++++++++++++++++++++--- 1 file changed, 116 insertions(+), 5 deletions(-) diff --git a/pvl_spa.m b/pvl_spa.m index 9304a6c..0f31024 100644 --- a/pvl_spa.m +++ b/pvl_spa.m @@ -60,8 +60,8 @@ % % delta_t = The current difference, in seconds, between international atomic time % (TAI) and UT1. This value is obtained by observation (see [2]). If -% omitted, the default value is 66.3 + 0.6175*(Year-2012). This -% default prediction can be changed as necessary. +% omitted, the default value is provided according to [3]. This default +% can be changed if necessary, but the default is preferred. % % Output Parameters: % SunAz = Azimuth of the sun in decimal degrees from North. 0 = North to 270 = West @@ -76,8 +76,10 @@ % Radiaiton Applications. NREL Report No. TP-560-34302. Revised % January 2008. http://www.nrel.gov/docs/fy08osti/34302.pdf as viewed % 2012-05-07 -% [2] DeltaT projections and current values from US Naval Observatory, -% http://asa.usno.navy.mil/SecK/DeltaT.html as viewed 2012-05-07 +% [2] Delta T and Universal Time, Fred Spenak, retrieved from +% http://eclipse.gsfc.nasa.gov/SEhelp/deltaT.html on 2016-06-22 +% [3] Polynomial Expressions for Delta T, retrieved from +% http://eclipse.gsfc.nasa.gov/SEhelp/deltatpoly2004.html on 2016-06-22 % % See also PVL_MAKETIMESTRUCT PVL_MAKELOCATIONSTRUCT PVL_ALT2PRES % PVL_GETAOI PVL_EPHEMERIS @@ -87,7 +89,7 @@ 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.addOptional('delta_t', 66.3 + .6175*(Time.year-2012), @(x) all(isvector(x) & isnumeric(x) & x>-8000 & x<8000)); +p.addOptional('delta_t', calculate_deltaT(Time.year, Time.month), @(x) all(isvector(x) & isnumeric(x) & x>-8000 & x<8000)); p.parse(Time, Location, varargin{:}); % Read in optional values from the struct p @@ -441,6 +443,115 @@ end + +function [DeltaT]=calculate_deltaT(year, month) + + % The equations aren't meant to be used outside of the range of years + % [-1999, 3000] + if any(year < -1999 | year > 3000) + warning(['Delta T is unknown for years before -1999 and after 3000.' ... + ' DeltaT values will be calculated, but the calculations ' ... + 'are not intended to be used for these years.']) + end + + y = year + (month + 0.5) ./ 12; + + DeltaT = zeros(size(year)); % create a starting vector + + % Implement the piecewise polynomial expressions for delta T as provided + % by http://eclipse.gsfc.nasa.gov/SEhelp/deltatpoly2004.html on + % June 22, 2016. + + % It is curious that in the years < -500 and >= 2150 the equation uses + % the year itself rather than the decimal year as used by all other + % time periods, but this is what the NASA eclipse site says... + DeltaT(year < -500) = ... + -20 ... + + 32 .* (((year(year < -500) - 1820) ./ 100).^2); + DeltaT(year >= -500 & year < 500) = ... + 10583.6 ... + - 1014.41 .* (y(year >= -500 & year < 500) ./ 100) ... + + 33.78311 .* ((y(year >= -500 & year < 500) ./ 100)).^2 ... + - 5.952053 .* ((y(year >= -500 & year < 500) ./ 100)).^3 ... + - 0.1798452 .* ((y(year >= -500 & year < 500) ./ 100)).^4 ... + + 0.022174192 .* ((y(year >= -500 & year < 500) ./ 100)).^5 ... + + 0.0090316521 .* ((y(year >= -500 & year < 500) ./ 100)).^6; + DeltaT(year >= 500 & year < 1600) = ... + 1574.2 ... + - 556.01 .* ((y(year >= 500 & year < 1600) - 1000) ./ 100) ... + + 71.23472 .* ((y(year >= 500 & year < 1600) - 1000) ./ 100).^2 ... + + 0.319781 .* ((y(year >= 500 & year < 1600) - 1000) ./ 100).^3 ... + - 0.8503463 .* ((y(year >= 500 & year < 1600) - 1000) ./ 100).^4 ... + - 0.005050998 .* ((y(year >= 500 & year < 1600) - 1000) ./ 100).^5 ... + + 0.0083572073 .* ((y(year >= 500 & year < 1600) - 1000) ./ 100).^6; + DeltaT(year >= 1600 & year < 1700) = ... + 120 ... + - 0.9808 .* (y(year >= 1600 & year < 1700) - 1600) ... + - 0.01532 .* (y(year >= 1600 & year < 1700) - 1600).^2 ... + + ((y(year >= 1600 & year < 1700) - 1600).^3) ./ 7129; + DeltaT(year >= 1700 & year < 1800) = ... + 8.83 ... + + 0.1603 .* (y(year >= 1700 & year < 1800) - 1700) ... + - 0.0059285 .* (y(year >= 1700 & year < 1800) - 1700).^2 ... + + 0.00013336 .* (y(year >= 1700 & year < 1800) - 1700).^3 ... + - ((y(year >= 1700 & year < 1800) - 1700).^4) ./ 1174000; + DeltaT(year >= 1800 & year < 1860) = ... + 13.72 ... + - 0.332447 .* (y(year >= 1800 & year < 1860) - 1800) ... + + 0.0068612 .* (y(year >= 1800 & year < 1860) - 1800).^2 ... + + 0.0041116 .* (y(year >= 1800 & year < 1860) - 1800).^3 ... + - 0.00037436 .* (y(year >= 1800 & year < 1860) - 1800).^4 ... + + 0.0000121272 .* (y(year >= 1800 & year < 1860) - 1800).^5 ... + - 0.0000001699 .* (y(year >= 1800 & year < 1860) - 1800).^6 ... + + 0.000000000875 .* (y(year >= 1800 & year < 1860) - 1800).^7; + DeltaT(year >= 1860 & year < 1900) = ... + 7.62 ... + + 0.5737 .* (y(year >= 1860 & year < 1900) - 1860) ... + - 0.251754 .* (y(year >= 1860 & year < 1900) - 1860).^2 ... + + 0.01680668 .* (y(year >= 1860 & year < 1900) - 1860).^3 ... + - 0.0004473624 .* (y(year >= 1860 & year < 1900) - 1860).^4 ... + + ((y(year >= 1860 & year < 1900) - 1860).^5) ./ 233174; + DeltaT(year >= 1900 & year < 1920) = ... + -2.79 ... + + 1.494119 .* (y(year >= 1900 & year < 1920) - 1900) ... + - 0.0598939 .* (y(year >= 1900 & year < 1920) - 1900).^2 ... + + 0.0061966 .* (y(year >= 1900 & year < 1920) - 1900).^3 ... + - 0.000197 .* (y(year >= 1900 & year < 1920) - 1900).^4; + DeltaT(year >= 1920 & year < 1941) = ... + 21.20 ... + + 0.84493 .* (y(year >= 1920 & year < 1941) - 1920) ... + - 0.076100 .* (y(year >= 1920 & year < 1941) - 1920).^2 ... + + 0.0020936 .* (y(year >= 1920 & year < 1941) - 1920).^3; + DeltaT(year >= 1941 & year < 1961) = ... + 29.07 ... + + 0.407 .* (y(year >= 1941 & year < 1961) - 1950) ... + - (y(year >= 1941 & year < 1961) - 1950).^2 ./ 233 ... + + (y(year >= 1941 & year < 1961) - 1950).^3 ./ 2547; + DeltaT(year >= 1961 & year < 1986) = ... + 45.45 ... + + 1.067 .* (y(year >= 1961 & year < 1986) - 1975) ... + - (y(year >= 1961 & year < 1986) - 1975).^2 ./ 260 ... + - (y(year >= 1961 & year < 1986) - 1975).^3 ./ 718; + DeltaT(year >= 1986 & year < 2005) = ... + 63.86 ... + + 0.3345 .* (y(year >= 1986 & year < 2005) - 2000) ... + - 0.060374 .* (y(year >= 1986 & year < 2005) - 2000).^2 ... + + 0.0017275 .* (y(year >= 1986 & year < 2005) - 2000).^3 ... + + 0.000651814 .* (y(year >= 1986 & year < 2005) - 2000).^4 ... + + 0.00002373599 .* (y(year >= 1986 & year < 2005) - 2000).^5; + DeltaT(year >= 2005 & year < 2050) = ... + 62.92 ... + + 0.32217 .* (y(year >= 2005 & year < 2050) - 2000) ... + + 0.005589 .* (y(year >= 2005 & year < 2050) - 2000).^2; + DeltaT(year >= 2050 & year < 2150) = ... + -20 ... + + 32 .* ((y(year >= 2050 & year < 2150) - 1820) ./ 100).^2 ... + - 0.5628 .* (2150 - y(year >= 2050 & year < 2150)); + DeltaT(year >= 2150) = ... % Same equation as <-500, provided separately for clarity and expansion + -20 ... + + 32 .* (((year(year >= 2150) - 1820) ./ 100).^2); +end + %% % This function just loads up the tables necessary function [L0table, L1table, L2table, L3table, L4table, L5table,... From 3abd593cf1ac5f25463a4a3d0e051879e1b7547f Mon Sep 17 00:00:00 2001 From: DanRiley Date: Thu, 23 Jun 2016 11:06:23 -0600 Subject: [PATCH 3/3] Update pvl_spa.m --- pvl_spa.m | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/pvl_spa.m b/pvl_spa.m index 0f31024..7b9059e 100644 --- a/pvl_spa.m +++ b/pvl_spa.m @@ -445,7 +445,9 @@ end function [DeltaT]=calculate_deltaT(year, month) - +% This function generates delta T based upon year and month using +% piecewise polynomial expressions. + % The equations aren't meant to be used outside of the range of years % [-1999, 3000] if any(year < -1999 | year > 3000) @@ -972,4 +974,4 @@ 0 0 0 0]; -end \ No newline at end of file +end