File Exchange

## Vectorized Solar Azimuth and Elevation Estimation

version 1.1 (3.58 KB) by

Predict the topocentric solar position defined by geodetic lat, lon, Alt, and a universal time

4.33333
6 Ratings

Updated

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]

Serge

### Serge (view profile)

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

### Tunahan Liman (view profile)

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

### L. He (view profile)

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

### Per (view profile)

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

### Phil (view profile)

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

### Anthony Kendall (view profile)

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.