diff --git a/SunPosition/Contents.m b/SunPosition/Contents.m
new file mode 100644
index 0000000..922b4a6
--- /dev/null
+++ b/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.
diff --git a/SunPosition/EarthEphemeris.m b/SunPosition/EarthEphemeris.m
new file mode 100644
index 0000000..eb00c1b
--- /dev/null
+++ b/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
\ No newline at end of file
diff --git a/SunPosition/Ephem_coefficients.mat b/SunPosition/Ephem_coefficients.mat
deleted file mode 100644
index fceb9af..0000000
Binary files a/SunPosition/Ephem_coefficients.mat and /dev/null differ
diff --git a/SunPosition/Ephemeris.m b/SunPosition/Ephemeris.m
deleted file mode 100644
index dc3741d..0000000
--- a/SunPosition/Ephemeris.m
+++ /dev/null
@@ -1,57 +0,0 @@
-function [ declin, radiusvector, solar_lon, varargout ] = Ephemeris( matdates )
-%Ephemeris: [ declin, radiusvector, solar_lon [,eqTime] ] = Ephemeris( matdates )
-% using calculations from the NOAA Solar Calculator
-%   http://www.esrl.noaa.gov/gmd/grad/solcalc/
-%
-%input
-%   matdates - scalar or vector (or matrix) of MATLAB datenums (UTC)
-%
-%output
-%   declination & solar longitude in degrees
-%   radius vector in AU
-%   optional - equation of time
-
-% matdates to julian centuries
-assert(isnumeric(matdates),'argument must be numeric')
-dt = datetime(matdates(:),'ConvertFrom','datenum');
-G = (juliandate(dt)-2451545)/36525; % Julian century
-
-% calculations - the letters correspond to the collumns in the spreadsheets
-% at 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 = matdates(:)-floor(matdates(:)); % local time
-AB = mod(E*1440+V,1440); % true solar time, min
-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 needed
-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
-if nargout>0
-    varargout{1} = V;
-end
-
-end
\ No newline at end of file
diff --git a/SunPosition/Examples/.MATLABDriveTag b/SunPosition/Examples/.MATLABDriveTag
new file mode 100644
index 0000000..21132d2
--- /dev/null
+++ b/SunPosition/Examples/.MATLABDriveTag
@@ -0,0 +1 @@
+c34090b3-2335-418a-abe4-e4d93e142b01
\ No newline at end of file
diff --git a/SunPosition/Examples/Analemma.mlx b/SunPosition/Examples/Analemma.mlx
new file mode 100644
index 0000000..c500212
Binary files /dev/null and b/SunPosition/Examples/Analemma.mlx differ
diff --git a/SunPosition/Examples/GlobalDaylight.mlx b/SunPosition/Examples/GlobalDaylight.mlx
new file mode 100644
index 0000000..b177fe5
Binary files /dev/null and b/SunPosition/Examples/GlobalDaylight.mlx differ
diff --git a/SunPosition/Examples/PlotSolarPath.mlx b/SunPosition/Examples/PlotSolarPath.mlx
new file mode 100644
index 0000000..d593b99
Binary files /dev/null and b/SunPosition/Examples/PlotSolarPath.mlx differ
diff --git a/SunPosition/Examples/SunPositionNow.mlx b/SunPosition/Examples/SunPositionNow.mlx
new file mode 100644
index 0000000..81a2b89
Binary files /dev/null and b/SunPosition/Examples/SunPositionNow.mlx differ
diff --git a/SunPosition/ExcelSource/.MATLABDriveTag b/SunPosition/ExcelSource/.MATLABDriveTag
new file mode 100644
index 0000000..9ba6433
--- /dev/null
+++ b/SunPosition/ExcelSource/.MATLABDriveTag
@@ -0,0 +1 @@
+3428071b-a9a1-4380-a3f1-b1ad7458042b
\ No newline at end of file
diff --git a/SunPosition/ExcelSource/NOAA_Solar_Calculations_day.xlsx b/SunPosition/ExcelSource/NOAA_Solar_Calculations_day.xlsx
new file mode 100644
index 0000000..31b4656
Binary files /dev/null and b/SunPosition/ExcelSource/NOAA_Solar_Calculations_day.xlsx differ
diff --git a/SunPosition/azimuthPreference.m b/SunPosition/azimuthPreference.m
new file mode 100644
index 0000000..e9bac84
--- /dev/null
+++ b/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
\ No newline at end of file
diff --git a/SunPosition/daylength.m b/SunPosition/daylength.m
deleted file mode 100644
index bfdc8c8..0000000
--- a/SunPosition/daylength.m
+++ /dev/null
@@ -1,27 +0,0 @@
-function hours = daylength(declin, latitude)
-% omega = daylight(declin, latitude)
-%
-% length of day, does not consider refraction
-%
-% Input
-%   declin - solar declination, + in northen hemisphere
-%   latitude - in degrees, + in northern hemisphere
-%   (if both arguments are not scalars, they need to be same size)
-%
-% Output
-%   hours - length of daylight
-
-if ~isscalar(declin) && ~isscalar(latitude)
-    assert(isequal(size(declin),size(latitude)),...
-        'if not scalars, arguments must be of same size')
-end
-
-arg = -tand(latitude).*tand(declin);
-% reset argument for total darkness (>1) or total sunlight(<-1)
-t = abs(arg)>1;
-if nnz(t)
-    arg(t) = sign(arg(t));
-end
-hours = 2*acosd(arg)/15;
-
-end
\ No newline at end of file
diff --git a/SunPosition/daylight.m b/SunPosition/daylight.m
deleted file mode 100644
index 46382a2..0000000
--- a/SunPosition/daylight.m
+++ /dev/null
@@ -1,28 +0,0 @@
-function omega = daylight(doy, latitude)
-% omega = daylight(doy, latitude)
-%
-% length of day, does not consider refraction
-%
-% Input
-%   doy - calendar day of year (1-366)
-%   latitude - in degrees, + in northern hemisphere
-%   (if both arguments are not scalars, they need to be same size)
-%
-% Output
-%   omega - length of daylight in degrees of rotation (i.e. 360 = 24 hrs)
-
-if ~isscalar(doy) && ~isscalar(latitude)
-    assert(isequal(size(doy),size(latitude)),...
-        'if not scalars, arguments must be of same size')
-end
-
-declination = declination_simple(doy);
-arg = -tand(latitude).*tand(declination);
-% reset argument for total darkness (>1) or total sunlight(<-1)
-t = abs(arg)>1;
-if nnz(t)
-    arg(t) = sign(arg(t));
-end
-omega = 2*acosd(arg);
-
-end
\ No newline at end of file
diff --git a/SunPosition/declination_simple.m b/SunPosition/declination_simple.m
deleted file mode 100644
index 3bd6186..0000000
--- a/SunPosition/declination_simple.m
+++ /dev/null
@@ -1,13 +0,0 @@
-% approximate (better than 1/10th degree) solar declination in degrees
-function dec=declination_simple(doy)
-% dec=declination_simple(doy)
-% input doy numerical day of year, starting 1 at 01 January
-A0=0.292033450096867;
-A=[-22.982158535385800 -0.287744385009193 -0.134810057611257]';
-B=[3.751757849830815 0.063793006836552 0.103899662868187]';
-omega=0.017115168775574;
-
-angles=doy*omega;
-dec=A0+A(1)*cos(angles)+B(1)*sin(angles)+A(2)*cos(2*angles)+...
-    B(2)*sin(2*angles)+A(3)*cos(3*angles)+B(3)*sin(3*angles);
-end
\ No newline at end of file
diff --git a/SunPosition/doc/.MATLABDriveTag b/SunPosition/doc/.MATLABDriveTag
new file mode 100644
index 0000000..e23a0bc
--- /dev/null
+++ b/SunPosition/doc/.MATLABDriveTag
@@ -0,0 +1 @@
+9d13a2a4-65cf-4ec6-8096-e8895ca916e0
\ No newline at end of file
diff --git a/SunPosition/doc/GettingStarted.mlx b/SunPosition/doc/GettingStarted.mlx
new file mode 100644
index 0000000..d9ab23f
Binary files /dev/null and b/SunPosition/doc/GettingStarted.mlx differ
diff --git a/SunPosition/kasten.m b/SunPosition/kasten.m
deleted file mode 100644
index 423db6e..0000000
--- a/SunPosition/kasten.m
+++ /dev/null
@@ -1,28 +0,0 @@
-function airmass = kasten( mu0 )
-% airmass = kasten( mu0 )
-%Kasten Relative optical airmass from cosine of the solar zenith angle
-%   Equation from Kasten, F., and A. T. Young (1989), Revised optical air
-%   mass tables and approximation formula, Applied Optics, 28, 4735-4738,
-%   doi: 10.1364/AO.28.004735.
-
-% coefficients
-a = 0.50572;
-b = 6.07995;
-c = 1.6364;
-
-% solar elevation in degrees, from horizon upward
-gam = asind(mu0);
-
-% Kasten-Young equation - set to NaN if below horizon
-t = mu0 < 0;
-if nnz(t)
-    airmass = NaN(size(mu0));
-    airmass(~t) = 1 ./(sind(gam(~t))+a.*(gam(~t)+b).^(-c));
-else
-    airmass = 1 ./(sind(gam)+a.*(gam+b).^(-c));
-end
-
-% slight correction for overhead sun
-airmass(airmass<1) = 1;
-
-end
\ No newline at end of file
diff --git a/SunPosition/radius_vector.m b/SunPosition/radius_vector.m
deleted file mode 100644
index 50cebe1..0000000
--- a/SunPosition/radius_vector.m
+++ /dev/null
@@ -1,20 +0,0 @@
-function rv = radius_vector( doy )
-%radius_vector: approximation for Earth-Sun distance in AU
-%
-% Input
-%   doy - calendar day-of year (can be a vector of days)
-%
-% Output
-%   radius vector
-
-a0 = 1.000128376223074;
-a1 = -0.016662032168947;
-b1 = -0.001013122469254;
-a2 = -1.2783185533524e-4;
-b2 = -2.152377408901758e-5;
-w = 0.017191017439612;
-
-rv = a0 + a1*cos(doy*w) + b1*sin(doy*w) + a2*cos(2*doy*w) + b2*sin(2*doy*w);
-
-end
-
diff --git a/SunPosition/sunRiseSet.m b/SunPosition/sunRiseSet.m
index 25e5650..2d2150f 100644
--- a/SunPosition/sunRiseSet.m
+++ b/SunPosition/sunRiseSet.m
@@ -1,21 +1,30 @@
 function [ sunrise, sunset ] = sunRiseSet( lat,lon,matdates,varargin )
 % [ sunrise, sunset ] = sunRiseSet( lat,lon,matdates [,P,T] )
-%sunRiseSet UTC times of sunrise and sunset
+%sunRiseSet times of sunrise and sunset, as MATLAB datetimes
 %
 % Input
 %   lat,lon - latitudes and longitudes, must be same size if not scalars
-%   matdates - dates in MATLAB datenum format, UTC but part of day ignored,
-%       must be same size as lat,lon if not scalar and they are not scalar
-% Optional input, must have both if used, and causes refraction to be
-% considered
-%   P,T - pressure (kPa) and temperature (K)
+%   matdates - dates in MATLAB datenum or datetime format, hour and
+%       minutes are ignored
+%       if datetime format, specify the TimeZone, otherwise UTC is assumed
+%       must be scalar or same size as lat,lon
+% Optional input, must have both if used, and causes refraction that depends
+%   on pressure and temperature to be considered
+%   P,T - pressure (hPa) and temperature (K)
+%
+% Output
+%   MATLAB datetimes corresponding to the input dates in the time zone
+%       specified in the input
+%   if 24 hr daylight, sunrise and sunset are set to the beginning and
+%       ending of the day
+%   if 24 hr darkness, sunrise and sunset are both set to NaT
 
 % arguments
 narginchk(3,5)
 p = inputParser;
-addRequired(p,'lat',@isnumeric)
-addRequired(p,'lon',@isnumeric)
-addRequired(p,'matdates',@isnumeric)
+addRequired(p,'lat',@(x) isnumeric(x) && all(abs(x(:))<=90))
+addRequired(p,'lon',@(x) isnumeric(x) && all(abs(x(:))<=180))
+addRequired(p,'matdates',@(x) isnumeric(x) || isdatetime(x))
 addOptional(p,'P',[],@isnumeric)
 addOptional(p,'T',[],@isnumeric)
 parse(p,lat,lon,matdates,varargin{:})
@@ -23,48 +32,81 @@
 T = p.Results.T;
 assert(~xor(isempty(P),isempty(T)),...
     'if P or T are used, both must be specified')
-assert(isequal(size(lat),size(lon)),...
-    'if lat and lon not scalars, must be same size')
+lat = p.Results.lat;
+lon = p.Results.lon;
+if isdatetime(p.Results.matdates)
+    matdates = p.Results.matdates;
+    if isempty(matdates.TimeZone)
+        matdates.TimeZone = 'Etc/GMT';
+    end
+    % eliminate the part of the day
+    dn = datenum(matdates);
+    matdates = datetime(floor(dn),'ConvertFrom','datenum',...
+        'TimeZone',matdates.TimeZone);
+else
+    matdates = datetime(floor(p.Results.matdates),'ConvertFrom','datenum',...
+        'timezone','Etc/GMT');
+end
 assert((isscalar(lat) || isscalar(matdates)) ||...
     isequal(size(lat),size(matdates)),...
     'if lat/lon and matdates not scalars, must be same size')
 
 % make same size if one is a scalar
-if xor(isscalar(lat),isscalar(matdates))
-    if isscalar(lat)
-        lat = ones(size(matdates))*lat;
-        lon = ones(size(matdates))*lon;
-    else
-        matdates = ones(size(lat))*matdates;
-    end
-end
+[lat,lon,matdates] = checkSizes(lat,lon,matdates);
 if ~isempty(P)
     assert(isscalar(P) || isequal(size(P),size(lat)),...
         'if P is not scalar, must be same size as lat/lon or matdates')
-    if isscalar(P) && ~isscalar(lat)
-        P = ones(size(lat))*P;
-        T = ones(size(lat))*T;
-    end
 end
 % hold size in case we vectorize
 N = size(lat);
 lat = lat(:);
 lon = lon(:);
 matdates = matdates(:);
+
+% values just based on declination at time zone at noon (not accounting for
+% change in declination during the day)
+[declin,~,~,eqT] = EarthEphemeris(matdates+hours(12));
+% refraction correction
 refrac = ~isempty(P);
 if refrac
     P = P(:);
     T = T(:);
 end
-
-r = zeros(size(lat));
-s = zeros(size(r));
-for k=1:length(lat)
-    if refrac
-        [r(k),s(k)] = sunUpDown(lat(k),lon(k),matdates(k),P(k),T(k));
+if refrac
+    if isscalar(P)
+        muR = fzero(@xref,0);
     else
-        [r(k),s(k)] = sunUpDown(lat(k),lon(k),matdates(k));
+        assert(isequal(size(P),size(lat)),'P and lat must be of same sizes')
+        muR = zeros(size(lat));
+        % P & T are passed to xref internally, so must do one by one
+        holdP = P;
+        holdT = T;
+        for k=1:length(holdP)
+            P = holdP(k);
+            T = holdT(k);
+            muR(k) = fzero(@xref,0);
+        end
     end
+else
+    muR = 0;
+end
+arg = muR./(cosd(lat).*cosd(declin))-tand(lat).*tand(declin);
+r = NaT(size(lat),'TimeZone',matdates.TimeZone);
+s = NaT(size(lat),'TimeZone',matdates.TimeZone);
+ts = arg<-1; %total sunlight
+td = arg>1; %total darkness, leave as NaT
+t = ~ts & ~td; %all else
+tz = hoursFromTimeZone(matdates);
+solarNoon = (720-4*lon-eqT+tz*60)/1440;
+if nnz(t)
+    HA = acosd(arg(t));
+    r(t) = matdates(t)+days(solarNoon(t)-4*HA/1440);
+    s(t) = matdates(t)+days(solarNoon(t)+4*HA/1440);
+end
+if nnz(ts)
+    HA = 180;
+    r(ts) = matdates(ts)+days(solarNoon(ts)-4*HA/1440);
+    s(ts) = r(ts)+days(1);
 end
 
 % resize if necessary
@@ -75,24 +117,21 @@
     sunrise = r;
     sunset = s;
 end
-
+    function mu_Refracted = xref(mu0)
+        mu_Refracted = refracted(mu0,P,T);
+    end
 end
 
-function [r,s] = sunUpDown(lat,lon,d,varargin)
-if nargin>3
-    P = varargin{1};
-    T = varargin{2};
-    refrac = true;
+function hz=hoursFromTimeZone(d)
+assert(isdatetime(d),'argument must be a datetime')
+if isempty(d.TimeZone)
+    hz = 0;
 else
-    refrac = false;
-end
-
-[declin,~,~,eqT] = Ephemeris(d);
-solarNoon = (720-4*lon-eqT)/1440;
-arg = cosd(90.833)./(cosd(lat).*cosd(declin))-tand(lat).*tand(declin);
-arg(arg<-1) = -1;
-arg(arg>1) = 1;
-HA = acosd(arg); % hour angle, sunrise
-r = solarNoon-HA*4/1440;
-s = solarNoon+HA*4/1440;
+    d2 = d;
+    d2.TimeZone = 'Etc/UTC';
+    [yy,mm,dd,hh,m,s] = datevec(d);
+    tempD = datetime(yy,mm,dd,hh,m,s,'TimeZone',d2.TimeZone);
+    [H,M,~] = hms(tempD-d2);
+    hz = H+M/60;
 end
+end
\ No newline at end of file
diff --git a/SunPosition/sunang.m b/SunPosition/sunang.m
index 503afde..47df56f 100644
--- a/SunPosition/sunang.m
+++ b/SunPosition/sunang.m
@@ -1,7 +1,6 @@
 function [mu0, phi0, airmass] = sunang( lat, lon, declin, omega, varargin )
 % [mu0, phi0, airmass] = sunang( lat, lon, declin, omega, ... )
-%Sun angles on horizontal surface, with optional correction for
-%   refraction
+%Sun angles on horizontal surface, with optional correction for refraction
 %   (uses the distance function from the Mapping Toolbox)
 %
 %  Input (in degrees)
@@ -10,43 +9,57 @@
 %   declin, solar declination
 %   omega, longitude at which sun is vertical
 %
-%  Optional arguments that go after omega, in the following order:
-%   zeroflag - if true, sets output negative cosines to zero and their
+%  Optional arguments as name-value pairs that go after omega, case-
+%   insensitive, in any order:
+%   'negzero' - if true, sets output negative cosines to zero and their
 %       corresponding azimuths to NaN, default true
 %       if false, returns negative cosines and their azimuths
-%   (next optional arguments are a pair, supply neither or both)
-%   P - pressure, to correct for atmospheric path length and
-%       refraction, default ignore refraction
-%   T - temperature in Kelvin, to correct for atmospheric path length and
-%   	refraction, default ignore refraction
+%   (next optional arguments are a pair to correct for refraction, supply
+%   neither or both, default is to ignore refraction)
+%   'P' - pressure in hPa (same as mb)
+%   'T' - temperature in Kelvin
 %  Output:
 %   mu0, cosine of solar zenith angle
-%   phi0, solar azimuth (degrees, from south, + ccw)
+%   phi0, solar azimuth (degrees counter-clockwise from south or degrees
+%       clocksie from north, depending on how azimuthPreference is set)
 %   airmass - relative atmospheric path length, where 1.0 is the path
 %       length at solar zenith angle = 0
+%
+%Examples
+%   [declin,~,solar_lon] = EarthEphemeris(datetime('2020-06-30 5:30','TimeZone','Etc/GMT+8'))
+%   [mu0,azm,airmass] = sunang(38,-119,declin,solar_lon) % w/o refraction
+%   [mu0,azm,airmass] = sunang(38,-119,declin,solar_lon,'P',1000,'T',288) % w refraction
+%   [mu0,azm,airmass] = sunang(30:5:50,-119,declin,solar_lon) % multiple latitudes
+%   sun below horizon
+%   [declin,~,solar_lon] = EarthEphemeris(datetime('2020-06-30 2:30','TimeZone','Etc/GMT+8'))
+%   [mu0,azm,airmass] = sunang(38,-119,declin,solar_lon) % cosine set to zero
+%   [mu0,azm,airmass] = sunang(38,-119,declin,solar_lon,'negzero',false) % cosine negative
 
-% check for, and if necessary process, optional arguments
-optargin=size(varargin,2);
-assert( nargin-optargin >= 3,...
-    'need 3 required arguments: latitude, declination, omega')
-if optargin >= 1
-    zeroflag = varargin{1};
-    if optargin == 2
-        error('if you specify pressure, you must specify temperature')
-    elseif optargin == 3
-        P = varargin{2};
-        T = varargin{3};
-    elseif optargin ~= 1
-        error('too many optional arguments: can be 1 or 3')
-    end
-else
-    % default is to set negative cosines to zero
-    zeroflag = true;
-end
+p = inputParser;
+addRequired(p,'lat',@(x) isnumeric(x) && all(abs(x(:))<=90))
+addRequired(p,'lon',@(x) isnumeric(x) && all(abs(x(:))<=180))
+addRequired(p,'declin',@(x) isnumeric(x) && all(abs(x(:))<=23.6))
+addRequired(p,'omega',@(x) isnumeric(x) && all(abs(x(:))<=180))
+% default is to set negative cosines to zero
+addParameter(p,'negzero',true,@(x) isscalar(x) && (isnumeric(x) || islogical(x)))
+% default is to not account for refraction
+addParameter(p,'P',[],@(x) isnumeric(x) && all(x(:)>0))
+addParameter(p,'T',[],@(x) isnumeric(x) && all(x(:)>0))
+parse(p,lat,lon,declin,omega,varargin{:});
+
+negzero = logical(p.Results.negzero);
+convertAzm = azimuthPreference;
+
+assert(~xor(isempty(p.Results.P),isempty(p.Results.T)),...
+    'if you specify pressure or temperature, you must specify both')
 
 if ~isscalar(lat)
-    assert(isequal(size(lat),size(lon)),...
-        'if not scalars, lat & lon must be same size')
+    if ~isscalar(lon)
+        assert(isequal(size(lat),size(lon)),...
+            'if not scalars, lat & lon must be same size')
+    else
+        [lat,lon] = checkSizes(lat,lon);
+    end
     if ~isscalar(declin)
         assert(isequal(size(declin),size(omega),size(lat)),...
             'if not scalars, decline & omega must be same size as lat/lon')
@@ -55,24 +68,55 @@
 
 [arclen, phi0] = distance(lat, lon, declin, omega);
 mu0 = cosd(arclen);
-% translate so that 0 is south, positive counter-clockwise
-phi0 = 180-phi0;
+% if convertAzm is true, translate so that 0 is south, positive counter-clockwise
+% otherwise MATLAB convention is clockwise from north
+if convertAzm
+    phi0 = 180-phi0;
+end
 
 % atmospheric refraction
-if optargin == 3
-    mu0 = refracted(mu0, P/10, T);
+if ~isempty(p.Results.P)
+    mu0 = refracted(mu0, p.Results.P, p.Results.T);
 end
 
 % relative airmass
 airmass = kasten(mu0);
 
 % set negative cosines to zero
-if zeroflag
+if negzero
     t = mu0 < 0;
     if nnz(t)
         mu0(t) = 0;
         phi0(t) = NaN;
     end
 end
+end
+
+function airmass = kasten( mu0 )
+% airmass = kasten( mu0 )
+%Kasten Relative optical airmass from cosine of the solar zenith angle
+%   Equation from Kasten, F., and A. T. Young (1989), Revised optical air
+%   mass tables and approximation formula, Applied Optics, 28, 4735-4738,
+%   doi: 10.1364/AO.28.004735.
+
+% coefficients
+a = 0.50572;
+b = 6.07995;
+c = 1.6364;
+
+% solar elevation in degrees, from horizon upward
+gam = asind(mu0);
+
+% Kasten-Young equation - set to NaN if below horizon
+t = mu0 < 0;
+if nnz(t)
+    airmass = NaN(size(mu0));
+    airmass(~t) = 1 ./(sind(gam(~t))+a.*(gam(~t)+b).^(-c));
+else
+    airmass = 1 ./(sind(gam)+a.*(gam+b).^(-c));
+end
+
+% slight correction for overhead sun
+airmass(airmass<1) = 1;
 
 end
\ No newline at end of file
diff --git a/SunPosition/sunhorz.m b/SunPosition/sunhorz.m
deleted file mode 100644
index ef2a148..0000000
--- a/SunPosition/sunhorz.m
+++ /dev/null
@@ -1,49 +0,0 @@
-function omega = sunhorz( latitude, longitude, declination, varargin )
-% omega = sunhorz( lat, lon, declin, [P,T] )
-%sunhorz longitudes of sun at sunrise & sunset with option to account for
-%       refraction
-%  Input, all in degrees, all scalar
-%   latitude
-%   longitude
-%   declination
-%
-%  Variable scalar input arguments, must be 2 to account for refraction
-%   P - atmospheric pressure, kPa
-%   T - temperature, Kelvin
-%
-%  Output:
-%   omega, vector of length 2, solar longitudes at sunrise & sunset
-%   (east of location at sunrise, west at sunset)
-
-global LAT DECLIN PRESS TEMP ZEROFLAG
-
-assert(isscalar(latitude) && isscalar(longitude) && isscalar(declination),...
-    'all inputs must be scalars')
-arg = -tand(latitude).*tand(declination); % now scalars, but may incorporate vectors later
-% reset argument for total darkness (>1) or total sunlight(<-1)
-if abs(arg)>1
-    arg = sign(arg);
-end
-
-optarg = size(varargin,2);
-assert(optarg == 0 || optarg ==2, 'must have 0 or 2 optional arguments')
-if optarg == 0
-    % standard method, ignoring atmospheric refraction
-    o = acosd(arg);
-else
-    % numerical method accounting for atmospheric refraction
-    PRESS = varargin{1};
-    TEMP = varargin{2};
-    ZEROFLAG = false;
-    LAT = latitude;
-    DECLIN = declination;  
-    o = fzero(@psunangr,acosd(arg));
-end
-
-omega = [longitude+o longitude-o];
-t = omega>180;
-omega(t) = omega(t)-180;
-t = omega<-180;
-omega(t) = abs(180+omega(t));
-
-end
\ No newline at end of file
diff --git a/SunPosition/sunslope.m b/SunPosition/sunslope.m
index 4079d97..09325f9 100644
--- a/SunPosition/sunslope.m
+++ b/SunPosition/sunslope.m
@@ -1,14 +1,28 @@
 function mu = sunslope( mu0,phi0,S,A )
-%sunslope cosine of illumination angle on slope
-%
 % mu = sunslope( mu0,phi0,S,A )
+%sunslope cosine of illumination angle on slope
 %
 % Input (angles in degrees)
 %   mu0, cosine of sun angle on flat surface
-%   phi0, sun azimuth, degrees, +ccw from south
+%   phi0, sun azimuth, degrees, direction set by the azimuthPreference function
 %   S, slope angle, degrees, from horizontal
-%   A, slope azimuth, degrees, +ccw from south
+%   A, slope azimuth, degrees, direction set by the azimuthPreference function
+%
+% Examples, the values for A are [-135 45] if your azimuthPreference is set
+%   to convert to counter-clockwise from south, + east, - west.
+%   If your azimuth preference is set to the MATLAB convention, clockwise from
+%   0° north, the A should be [315 145]
+%   Northern Hemisphere summer morning, Mt Blanc, NW and SE slopes
+%   [declin,~,sol_lon] = EarthEphemeris(datetime('2020-07-04 10:00','TimeZone','Europe/Paris'))
+%   [mu0,phi0] = sunang(45.8328,6.865,declin,sol_lon)
+%   mu = sunslope(mu0,phi0,35,A)
+%   Northern Hemisphere winter morning, Mt Blanc, NW and SE slopes
+%   [declin,~,sol_lon] = EarthEphemeris(datetime('2020-01-04 10:00','TimeZone','Europe/Paris'))
+%   [mu0,phi0] = sunang(45.8328,6.865,declin,sol_lon)
+%   mu = sunslope(mu0,phi0,35,A)
+
 
+[mu0,phi0,S,A] = checkSizes(mu0,phi0,S,A);
 mu = mu0.*cosd(S) + sqrt((1-mu0).*(1+mu0)).*sind(S).*cosd(phi0-A);
 if isscalar(mu) && mu<0
     mu = 0;
diff --git a/SunPosition/util/.MATLABDriveTag b/SunPosition/util/.MATLABDriveTag
new file mode 100644
index 0000000..74b9294
--- /dev/null
+++ b/SunPosition/util/.MATLABDriveTag
@@ -0,0 +1 @@
+d16f2717-c99c-4825-94a0-51ac6a46ef01
\ No newline at end of file
diff --git a/SunPosition/refracted.m b/SunPosition/util/refracted.m
similarity index 67%
rename from SunPosition/refracted.m
rename to SunPosition/util/refracted.m
index 5913622..e3199d2 100644
--- a/SunPosition/refracted.m
+++ b/SunPosition/util/refracted.m
@@ -5,11 +5,12 @@
 %   algorithm for solar radiation applications (revised),
 %   NREL/TP-560-34302, 56 pp, National Renewable Energy Laboratory, Golden,
 %   CO, doi: 10.2172/15003974.
+%   (in their equation, T is in deg C)
 %
 % Input:    mu0 - cosine of solar zenith angle, unrefracted
-%           P - pressure in kPa
+%           P - pressure in hPa (same as mb)
 %           T - temperature in Kelvin
-% (OK to use annual averages for P & T)
+% (OK to use estimates for P & T as the sensitivity is small)
 %
 % Output: cosine of refracted solar zenith angle
 
@@ -19,25 +20,25 @@
 if t==0
     assert(isequal(size(mu0),size(P),size(T)),...
         'all inputs vectors, so they must be the same size')
-elseif t==1
+elseif t==1 % 2 of the 3 inputs are not scalar
     if isscalar(T)
         assert(isequal(size(mu0),size(P)),...
             'mu0 and P are both vectors but of different sizes')
     elseif isscalar(P)
         assert(isequal(size(mu0),size(T)),...
             'mu0 and T are both vectors but of different sizes')
-    else if isscalar(mu0)
-            assert(isequal(size(T),size(P)),...
-                'P and T are both vectors but of different sizes')
-        end
+    elseif isscalar(mu0)
+        assert(isequal(size(T),size(P)),...
+            'P and T are both vectors but of different sizes')
     end
 end
 
 % check P & T within range
-tp = P > 200;
-assert(nnz(tp)==0, 'pressure must be in kPa, max value is %g', max(P));
-tt = T<50;
-assert(nnz(tt)==0, 'Temp in Kelvin, min value is %g', min(T));
+tp = P>1200 | P<200;
+assert(nnz(tp)==0, 'pressure must be in hPa, current range is [%g %g]',...
+    min(P(:)), max(P(:)))
+tt = T<30;
+assert(nnz(tt)==0, 'Temp in Kelvin, current minimum value is %g', min(T(:)))
 
 % Code following takes advantage of the sin/cos relationship,
 % i.e. sind(90-x) = cosd(x) etc.
@@ -46,7 +47,7 @@
 e0 = asind(mu0);
 
 % refraction - equation (42) from Reda-Andreas
-delE = (P/101).*(283./T).*(1.02./(60*tand(e0 + 10.3./(e0+5.11))));
+delE = (P/1010).*(283./T).*(1.02./(60*tand(e0 + 10.3./(e0+5.11))));
 
 % cosine of refracted solar zenith angle
 muR = sind(e0+delE);