updated SunPosition

SunPosition/Contents.m

@@ -0,0 +1,26 @@
+% SUNPOSITION
+% Functions
+%   EarthEphemeris - Returns the solar declination, Earth-Sun radius vector
+%       (distance in astronomical units), solar longitude, and optionally the
+%       equation of time.
+%   sunang - Returns cosine of sun angle and its azimuth in degrees, either
+%       clockwise from north or counter-clockwise from south, and airmass
+%       (path length through atmosphere accounting for Earth curvature).
+%       Optionally accounts for refraction.
+%   sunRiseSet - Times of sunrise and sunset.
+%   sunslope - Cosine of illumination angle on slope
+%
+%About Orientation for Azimuths for sun angles and slopes
+%   Years ago when I started coding for solar geometry in mountains, my go-to
+%   text was William Sellers’ Physical Climatology (1965). Based on his conventions,
+%   I adopted 0° toward South, positive to the East and negative to the West,
+%   putting North at ±180°. This convention is consistent with a right-hand
+%   coordinate system, which we conventionally use by making latitude positive
+%   to the North and longitude positive to the East.
+%   I realize that another common convention puts 0° toward North, with other
+%   azimuths clockwise from north, so I have added this option to sunang.
+%   Set your preference in the function azimuthPreference. Because this preference
+%   is likely persistent, you can set it in the function and not worry about
+%   it again.
+%   The sunslope.m code will work without modification as long as the slope
+%   azimuth uses the same convention as the solar azimuth.

SunPosition/EarthEphemeris.m

@@ -0,0 +1,83 @@
+function [ declin, radiusvector, solar_lon, varargout ] = EarthEphemeris( matdates )
+%EarthEphemeris: [ declin, radiusvector, solar_lon [,eqTime] ] = EarthEphemeris( matdates )
+% using calculations from the NOAA Solar Calculator
+%   http://www.esrl.noaa.gov/gmd/grad/solcalc/
+% calculates the subsolar location and distance from sun for input dates/times
+%
+%input
+%   matdates - scalar or vector or matrix, either of MATLAB datenums (UTC)
+%       or of MATLAB datetimes
+%       If datenum, time zone is assumed to be UTC.
+%       If datetime, time zone should be specified and will be converted to
+%       UTC internally, but is assumed to be UTC if not specified.
+%
+%output
+%   declination & solar longitude in degrees
+%   radius vector in AU
+%   optional - equation of time
+%
+%Examples
+%   [declin,rv,sollon,eqT] = EarthEphemeris(now) % assumes UTC
+%   [declin,rv,sollon,eqT] = EarthEphemeris(datetime(now,'ConvertFrom','datenum','TimeZone','local'))
+%
+%See also
+%   sunang, sunRiseSet, sunslope
+
+% matdates to julian centuries
+p = inputParser;
+addRequired(p,'matdates',@(x) isnumeric(x) || isdatetime(x))
+parse(p,matdates)
+if isnumeric(matdates)
+    dt = datetime(matdates(:),'ConvertFrom','datenum');
+    dt.TimeZone = 'Etc/GMT';
+else
+    dt = matdates(:);
+end
+if isempty(dt.TimeZone)
+    warning('Time Zone not specified, assuming GMT (numerically same as UTC)')
+    dt.TimeZone = 'Etc/GMT';
+else
+    % convert to GMT
+    dt.TimeZone = 'Etc/GMT';
+end
+G = (juliandate(dt)-2451545)/36525; % Julian century
+
+% calculations - the letters correspond to the columns in the spreadsheets
+% in the NOAA Solar Calculator
+I = mod(280.46646+G.*(36000.76983 + G*0.0003032),360); % geom mean long Sun, degrees
+J = 357.52911+G.*(35999.05029 - 0.0001537*G); % geom mean Sun anomaly, degrees
+K = 0.016708634-G.*(0.000042037+0.0000001267*G); % Earth orbit eccentricity
+L = sind(J).*(1.914602-G.*(0.004817+0.000014*G))+sind(2*J).*...
+    (0.019993-0.000101*G)+sind(3*J)*0.000289; % Sun eq of Ctr
+M = I+L; % Sun true longitude
+N = J+L; % Sun true anomaly
+radiusvector = (1.000001018*(1-K.^2))./(1+K.*cosd(N));
+P = M-0.00569-0.00478*sind(125.04-1934.136*G); % Sun apparent long
+Q = 23+(26+((21.448-G.*(46.815+G.*(0.00059-G*0.001813))))/60)/60; % mean obliq ecliptic
+R = Q+0.00256*cosd(125.04-1934.136*G); % obliquity correction
+T = asind(sind(R).*sind(P)); % declination
+U = tand(R/2).^2; % var y
+V = 4*rad2deg(U.*sind(2*I)-2*K.*sind(J)+4*K.*U.*sind(J).*cosd(2*I)-...
+    0.5*U.^2.*sind(4*I)-1.25*K.^2.*sind(2*J)); % eq of time
+E = datenum(dt(:))-floor(datenum(dt(:))); % UTC time in fraction of a day
+AB = mod(E*1440+V,1440); % true solar time, minutes, corrected for equation of time
+t = AB<0;
+solar_lon = AB/4;
+solar_lon(t) = -(solar_lon(t)+180);
+solar_lon(~t) = -(solar_lon(~t)-180);
+declin = T;
+
+% revert to matrix, if input is a matrix
+if ~isequal(size(matdates),size(declin))
+    declin = reshape(declin,size(matdates));
+    solar_lon = reshape(solar_lon,size(matdates));
+    radiusvector = reshape(radiusvector,size(matdates));
+    V = reshape(V,size(matdates));
+end
+
+% equation of time, optional output
+if nargout>3
+    varargout{1} = V;
+end
+
+end

SunPosition/Examples/.MATLABDriveTag

@@ -0,0 +1 @@
+c34090b3-2335-418a-abe4-e4d93e142b01

SunPosition/ExcelSource/.MATLABDriveTag

@@ -0,0 +1 @@
+3428071b-a9a1-4380-a3f1-b1ad7458042b

SunPosition/azimuthPreference.m

@@ -0,0 +1,13 @@
+function convertAzm = azimuthPreference()
+%MATLAB functions that return an azimuth use the convention clockwise from
+%north, i.e. 0° north, east is 90°, west is 270°
+%(but atan2d returns angles in the +/-180 degree range, so ...)
+%If you want to convert MATLAB azimuths to counter-clockwise from
+%south, i.e. 0° south, + to east, - to west, then set convertAzm to true.
+%If you want to use the MATLAB convertion for gradientm, set convertAzm to false.
+
+%comment out the option you don't want to use
+convertAzm = true; % convert azimuths to 0° south, + to east, - to west
+% convertAzm = false; % leave azimuths as clockwise from 0° north
+
+end

SunPosition/doc/.MATLABDriveTag

@@ -0,0 +1 @@
+9d13a2a4-65cf-4ec6-8096-e8895ca916e0