File Exchange

## Vectorized Solar Azimuth and Elevation Estimation

version 1.1.0.0 (3.58 KB) by
Predict the topocentric solar position defined by geodetic lat, lon, Alt, and a universal time

Updated 18 Apr 2013

Predict the azimuth and elevation of the Sun within +/- 1 degree at any geodetic latitude, longitude and altitude. Due to popular demand, this routine has been vectorized for speed.

Function Call: [Az El] = SolarAzEl('2008/02/18 13:10:00',60,15,0)

Input List:
UTC Date and Time - Use format YYYY/MM/DD hh:mm:ss or MATLAB date vector dimensions can be [N x 1]
Latitude - Site Latitude in degrees -90:90 -> S(-) N(+) dimensions can be [N x 1]
Longitude - Site Longitude in degrees -180:180 W(-) E(+) dimensions can be [N x 1]
Altitude - Site Altitude in km dimensions can be [N x 1]

Output List:
Az - Solar Azimuth angle in degrees [N x 1]
El - Solar Elevation/Altitude Angle in degrees [N x 1]

### Cite As

Darin Koblick (2020). Vectorized Solar Azimuth and Elevation Estimation (https://www.mathworks.com/matlabcentral/fileexchange/23051-vectorized-solar-azimuth-and-elevation-estimation), MATLAB Central File Exchange. Retrieved .

Khalid CHOUA

Serge

No correction for atmospheric refraction, which is around 0.5 deg when sun is on the horizon. So this routine give the Sun's true azel not apparent azel. Otherwise the accuracy appears to be much better then +-1 deg.

Ned Gulley

Thanks for the file! I used it to write this blog entry: https://blogs.mathworks.com/community/2018/01/24/when-is-noon/

Serge

Here is a metricised version that can be 80x faster in some cases.

function [Az,El] = SolarAzElq(UTC,Lat,Lon,Alt)
%Sun possition from time and location (matricised)
% [Az,El] = SolarAzEl(UTC,Lat,Lon,Alt)
%UTC: UTC Time (MatLab's datenum or 'yyyy-mm-dd HH:MM:SS' cellstr or char)
%Lat: Latitude [-90 90] (deg)
%Lon: Longitude [-180 180] (deg)
%Alt: Altitude above sea level, optional (km)
%Az: Azimuth location of the sun (deg)
%El: Elevation location of the sun (deg)
%
%Example:
% [Az,El] = SolarAzElq('1991-05-19 13:00:00',50,10,0)
%
%Example:
% [UTC,Lat,Lon] = ndgrid(730486:1/24:730487,-90:5:90,-180:5:180);
% tic,[Az,El] = SolarAzElq(UTC,Lat,Lon);toc
%
%References:
% http://stjarnhimlen.se/comp/tutorial.html#5
% http://www.stargazing.net/kepler/altaz.html RA,DEC to Az,Alt

%History:
% Darin C. Koblick 02/17/2009 Authos
% Darin C. Koblick 04/16/2013 Vectorized
% Serge Kharabash 09/02/2016 Metricised

%Test:
% [UTC,Lat,Lon] = ndgrid(730486:10:730852,-90:5:90,-180:5:180);
% tic,[Az1,El1] = SolarAzElq(UTC(:),Lat(:),Lon(:));toc,t1=toc;
% tic,[Az2,El2] = SolarAzEl(UTC(:),Lat(:),Lon(:),0);toc,t2=toc;
% max(abs(Az1-Az2)),max(abs(El1-El2)),t2/t1

%defaults
if nargin<4 || isempty(Alt), Alt = 0; end
d2r = pi/180; %radiance to degrees conversion factor
r2d = 180/pi; %radiance to degrees conversion factor

%checks
if ischar(UTC)
UTC = cellstr(UTC);
end
if iscell(UTC)
UTC = reshape(datenum(UTC(:),'yyyy-mm-dd HH:MM:SS'),size(UTC));
end

%julian date
[year,month,day,hour,min,sec] = datevec(UTC);
if ndims(UTC)>2 %#ok<ISMAT>
year = reshape(year ,size(UTC));
month = reshape(month,size(UTC));
day = reshape(day ,size(UTC));
hour = reshape(hour ,size(UTC));
min = reshape(min ,size(UTC));
sec = reshape(sec ,size(UTC));
end
[jd,UTH] = juliandate(year,month,day,hour,min,sec);
day = jd - 2451543.5;

%Keplerian elements for the Sun (geocentric)
w = 282.9404 + 4.70935e-5 * day; %longitude of perihelion degrees
e = 0.016709 - 1.151e-9 * day; %eccentricity
M = mod(356.0470 + 0.9856002585 * day, 360); %mean anomaly degrees
L = w + M; %Sun's mean longitude degrees
oblecl = (23.4393 - 3.563e-7 * day)*d2r; %Sun's obliquity of the ecliptic, rad

%auxiliary angle
E = M + r2d*e.*sin(M*d2r).*(1+e.*cos(M*d2r));

%rectangular coordinates in the plane of the ecliptic (x toward perhilion)
x = cos(E*d2r)-e;
year = sin(E*d2r).*sqrt(1-e.^2);

%distance and true anomaly
r = sqrt(x.^2 + year.^2);
v = atan2(year,x)*r2d;

%longitude of the sun
lon = v + w;

%ecliptic rectangular coordinates
xeclip = r.*cos(lon*d2r);
yeclip = r.*sin(lon*d2r);
zeclip = 0;

%rotate to equitorial rectangular coordinates
xequat = xeclip;
yequat = yeclip.*cos(oblecl) + zeclip*sin(oblecl);
zequat = yeclip.*sin(0.409115648642983) + zeclip*cos(oblecl);

%convert to RA and Dec
r = sqrt(xequat.^2 + yequat.^2 + zequat.^2) - (Alt/149598000); %roll up the altitude correction

%local siderial time
GMST0 = mod(L+180,360)/15;
SIDTIME = GMST0 + UTH + Lon/15;

%replace RA with hour angle HA
HA = 15*SIDTIME - RA * r2d;

%convert to rectangular coordinate system
x = cos(HA*d2r).*cos(delta);
year = sin(HA*d2r).*cos(delta);
z = sin(delta);

%rotate along an axis going east-west
xhor = x.*cos((90-Lat)*d2r) - z.*sin((90-Lat)*d2r);
yhor = year;
zhor = x.*sin((90-Lat)*d2r) + z.*cos((90-Lat)*d2r);

%find Az and El
Az = atan2(yhor,xhor) * r2d + 180;
El = asin(zhor) * r2d;

function [jd,UTH] = juliandate(year,month,day,hour,min,sec)
%calculate julian date & J2000 value
UTH = hour + min/60 + sec/3600; %J2000
idx = month <= 2;
year(idx) = year(idx) - 1;
month(idx) = month(idx) + 12;
jd = floor(365.25*(year+4716)) + floor(30.6001*(month+1)) + 2 - ...
floor(year/100) + floor(floor(year/100)/4) + day - 1524.5 + ...
UTH/24;

Tunahan Liman

I need a code something like that. How can I get the codes?

Superb. This function is well written and well documented. Thanks for sharing!

L. He

very nice program. I compare it with
http://www.mathworks.com/matlabcentral/fileexchange/4605-sun-position-m

the result is almost same. That one is more accurate,but not Vectorized yet.

jd = juliandate(datestr([y,mo,d,h,mi,s],'yyyy/mm/dd HH:MM:SS'),'yyyy/mm/dd HH:MM:SS');
do not use datestr but use date_num directly.
this will boost the speed, a lot

Per

Very useful function. With a few updates it can handle vector time input. I would prefer the use of matlab UTC time input in order to speed up.

Phil

Ropey when using vector times: (line 36 generates a vector eccentricity: line 42 then requires an edit to force array multiply not matrix multiply). Still unable to get vector time version to agree with loop version....
As Mr. Picky, I would prefer time argin to be Matlab datenum, not string.
HOWEVER, this is the only code I've found that gives Azimuth round the full 360: most are 0-180 and it's up to you to find if its in the east or west... due to using code like
http://en.wikipedia.org/wiki/Solar_azimuth_angle
Thanks

Anthony Kendall

Excellent function, it's fast, compact, and easily modified for my particular needs. Thank you very much! BTW, I compared it with sun position tables, (http://www.srrb.noaa.gov/highlights/sunrise/azel.html) and it does very well.