solar position calculator

solarPosition.m

function [angles,projection] = solarPosition(datetime,latitude,longitude,time_zone,rotation,dst)
%SOLARPOSITION Calculate solar position using most basic algorithm
%   This is the most basic algorithm. It is documented in Seinfeld &
%   Pandis, Duffie & Beckman and Wikipedia.
%
% [ANGLES,PROJECTION] = SOLARPOSITION(DATE,TIME,LATITUDE,LONGITUDE,TIME_ZONE)
% returns ZENITH & AZIMUTH for all DATE & TIME pairs at LATITUDE, LONGITUDE.
% ANGLES = [ZENITH,AZIMUTH] and PROJECTION = [PHI_X, PHI_Y]
% PHI_X is projection on x-z plane & PHI_Y is projection on y-z plane.
% DATETIME can be string, vector [YEAR, MONTH, DAY, HOURS, MINUTES, SECONDS],
%   cellstring or matrix N x [YEAR, MONTH, DAY, HOURS, MINUTES, SECONDS] for N
%   times.
% LATITUDE [degrees] and LONGITUDE [degrees] are the coordinates of the site.
% TIME_ZONE [hours] of the site.
% ROTATION [degrees] clockwise rotation of system relative to north.
% DST [logical] flag for daylight savings time, typ. from March to November
%   in the northern hemisphere.
%
% References:
% http://en.wikipedia.org/wiki/Solar_azimuth_angle
% http://en.wikipedia.org/wiki/Solar_elevation_angle
%
% Mark A. Mikofski
% Copyright (c) 2013
%

%% datetime
if iscellstr(datetime) || ~isvector(datetime)
    datetime = datenum(datetime); % [days] dates & times
else
    datetime = datetime(:); % convert datenums to row
end
date = floor(datetime); % [days]
[year,~,~] = datevec(date);
time = datetime - date; % [days]
%% constants
toRadians = @(x)x*pi/180; % convert degrees to radians
toDegrees = @(x)x*180/pi; % convert radians to degrees
%% Equation of time
d_n = mod(date-datenum(year,1,1)+1,365); % day number
B = 2*pi*(d_n-81)/365; % ET parameter
ET = 9.87*sin(2*B)-7.53*cos(B)-1.5*sin(B); % [minutes] equation of time
% approximate solar time
solarTime = ((time*24-double(dst))*60+4*(longitude-time_zone*15)+ET)/60/24;
latitude_rad = toRadians(latitude); % [radians] latitude
rotation_rad = toRadians(rotation); % [radians] field rotation
t_h = (solarTime*24-12)*15; % [degrees] hour angle
t_h_rad = toRadians(t_h); % [radians]
delta = -23.45 * cos(2*pi*(d_n+10)/365); % [degrees] declination
delta_rad = toRadians(delta); % [radians]
theta_rad = acos(sin(latitude_rad)*sin(delta_rad)+ ...
    cos(latitude_rad)*cos(delta_rad).*cos(t_h_rad)); % [radians] zenith
theta = toDegrees(theta_rad); % [degrees] zenith
elevation = 90 - theta; % elevation
day = elevation>0; % day or night?
cos_phi = (cos(theta_rad)*sin(latitude_rad)- ...
    sin(delta_rad))./(sin(theta_rad)*cos(latitude_rad)); % cosine(azimuth)
% azimuth [0, 180], absolute value measured from due south, so east = west = 90,
% south = 0, north = 180
phi_south = acos(min(1,max(-1,cos_phi)));
% azimuth [0, 360], measured clockwise from due north, so east = 90,
% south = 180, and west = 270 degrees
phi_rad = NaN(size(phi_south)); % night azimuth is NaN
% shift from ATAN to ATAN2, IE: use domain from 0 to 360 degrees instead of
% from -180 to 180
phi_rad(day) = pi + sign(t_h(day)).*phi_south(day); % Shift domain to 0-360 deg
% projection of sun angle on x-z plane, measured from z-direction (up)
phi_x = toDegrees(atan2(sin(phi_rad-rotation_rad).*sin(theta_rad), ...
    cos(theta_rad))); % [degrees]
% projection of sun angle on y-z plane, measured from z-direction (up)
phi_y = toDegrees(atan2(cos(phi_rad-rotation_rad).*sin(theta_rad), ...
    cos(theta_rad))); % [degrees]
phi = toDegrees(phi_rad); % [degrees] azimuth
angles = [theta, phi]; % [degrees] zenith, azimuth
projection = [phi_x,phi_y]; % [degrees] x-z plane, y-z plane
end

solarPosition_Example.m

%% initalize workspace
clear('all'),close('all'),clc

%% Solar Position Calculator Examples
% This example uses
% <http://www.mathworks.com/matlabcentral/fileexchange/32359-datefig-figure-handle-class-with-automatically-readjusting-date-ticks datefig>
% from the file exchange and a mex version of <http://www.nrel.gov/midc/solpos/ NREL's C/C++ source for SOLPOS>.

% The following examples show different uses of datetime and DST for
% Golden, Colorado.
lat = 39.5; % [arc-degrees] latitude
long = -103.8; % [arc-degrees] longitude
TZ = -7; % [hrs] offset from UTC, during standard time
rot = 0; % [arc-degrees] rotation clockwise from north

%% Datetime as cell array in local time during DayLight Savings Time

% cell array of local date-times during DST
datetimes = {'3/31/2009 10:00 AM','3/31/2009 11:00 AM'};
DST = true; % local time is daylight savings time
% calculate solar position
a = solarPosition(datetimes,lat,long,TZ,rot,DST); % [arc-degrees]
% compare to NREL SOLPOS calculator http://www.nrel.gov/midc/solpos/
a_nrel = a; 
for n = 1:size(a,1)
    dt_nrel = datevec(datetimes{n});dt_nrel(:,4) = dt_nrel(:,4)-DST;
    [a_nrel(n,:),~] = solpos([lat,long,TZ],dt_nrel,[1013,25]);
    night = a_nrel(:,1)>90;a_nrel(night,:) = NaN(sum(night),2);
end
fprintf('%-20s %-3s %-7s %-8s %-7s %-8s %-7s %-8s\n','date/time','DST', ...
    'Zenith','Azimuth','Ze_NREL','Az_NREL','Ze_Err %','Az_Err %')
fprintf('%-20s|%-3s|%-7s|%-8s|%-7s|%-8s|%-7s|%-8s|\n',repmat('-',1,20), ...
    '---','-------','--------','-------','--------','-------','--------')
for n = 1:size(a,1)
  fprintf('%20s %-3s %7.4f %8.4f %7.4f %8.4f %7.4f %8.4f\n',datetimes{n}, ...
      DST*'yes'+~DST*'no ',a(n,1),a(n,2),a_nrel(n,1),a_nrel(n,2), ...
      100*(a(n,1)-a_nrel(n,1))./a_nrel(n,1), ...
      100*(a(n,2)-a_nrel(n,2))./a_nrel(n,2))
end

fprintf('%-20s|%-3s|%-7s|%-8s|%-7s|%-8s|%-7s|%-8s|\n',repmat('-',1,20), ...
    '---','-------','--------','-------','--------','-------','--------')
fprintf('\n')

%% Datetime as matrix of date-vectors with DST removed.

% array of standard date-times, DST removed
datetimes = [2009,3,31,9,0,0;2009,3,31,10,0,0];
DST = false; % daylight savings time has been removed
% calculate solar position
a = solarPosition(datetimes,lat,long,TZ,rot,DST); % [arc-degrees]
% compare to NREL SOLPOS calculator http://www.nrel.gov/midc/solpos/
a_nrel = a; 
for n = 1:size(a,1)
    dt_nrel = datetimes(n,:);dt_nrel(:,4) = dt_nrel(:,4)-DST;
    [a_nrel(n,:),~] = solpos([lat,long,TZ],dt_nrel,[1013,25]);
    night = a_nrel(:,1)>90;a_nrel(night,:) = NaN(sum(night),2);
end

fprintf('%-20s %-3s %-7s %-8s %-7s %-8s %-7s %-8s\n','date/time','DST', ...
    'Zenith','Azimuth','Ze_NREL','Az_NREL','Ze_Err %','Az_Err %')
fprintf('%-20s|%-3s|%-7s|%-8s|%-7s|%-8s|%-7s|%-8s|\n',repmat('-',1,20), ...
    '---','-------','--------','-------','--------','-------','--------')
for n = 1:size(a,1)
  fprintf('%20s %-3s %7.4f %8.4f %7.4f %8.4f %7.4f %8.4f\n', ...
      datestr(datetimes(n,:)),DST*'yes'+~DST*'no ',a(n,1),a(n,2), ...
      a_nrel(n,1),a_nrel(n,2),100*(a(n,1)-a_nrel(n,1))./a_nrel(n,1), ...
      100*(a(n,2)-a_nrel(n,2))./a_nrel(n,2))
end
fprintf('%-20s|%-3s|%-7s|%-8s|%-7s|%-8s|%-7s|%-8s|\n',repmat('-',1,20), ...
    '---','-------','--------','-------','--------','-------','--------')
fprintf('\n')

%% Datetime as vector of date-numbers with DST removed for plotting zenith & azimuth

% linear sequence of `datenum` for one day at hourly intervals
datetimes = linspace(datenum([2009,3,31]),datenum([2009,4,1]),25); % [days]
DST = false; % daylight savings time has been removed
% calculate solar position
a = solarPosition(datetimes,lat,long,TZ,rot,DST); % [arc-degrees]
% compare to NREL SOLPOS calculator http://www.nrel.gov/midc/solpos/
a_nrel = a; 
for n = 1:size(a,1)
    dt_nrel = datevec(datetimes(n));dt_nrel(:,4) = dt_nrel(:,4)-DST;
    [a_nrel(n,:),~] = solpos([lat,long,TZ],dt_nrel,[1013,25]);
    night = a_nrel(:,1)>90;a_nrel(night,:) = NaN(sum(night),2);
end
ax = datefig.plotyy(datetimes,[a,a_nrel], ...
    datetimes,(a-a_nrel)./a_nrel*100);
title('Solar Position - 3/31/2009, Golden, CO [39.5\circ, -103.8\circ, -7.0]')
ylabel(ax(1),'angle [\circ]')
ylabel(ax(2),'rel. err. [%]')
legend('zenith','azimuth','ze_{NREL}','az_{NREL}','ze_{err}','az_{err}')

%% Errors
% Comparison of the solar position calculator with NREL's SOLPOS, shows
% good overall agreement. NREL's SOLPOS considers atmospheric refraction in
% the calculation of zenith, which accounts for some error in this solar
% calculator. Max relative error in zenith, including the error due to
% refraction is near than 2%. Error in azimuth is near zero.