Code covered by the BSD License

### Highlights from sun_position.m

4.75

4.8 | 24 ratings Rate this file 113 Downloads (last 30 days) File Size: 1.28 MB File ID: #4605

# sun_position.m

10 Mar 2004 (Updated )

Sun position given observer time/location.

### Editor's Notes:

This file was selected as MATLAB Central Pick of the Week

File Information
Description

sun = sun_position(time, location)

This function computes the sun position (zenith and azimuth angle at the observer location) as a function of the observer local time and position.

It is an implementation of the algorithm presented by Reda et Andreas in:
Reda, I., Andreas, A. (2003) Solar position algorithm for solar radiation application. National Renewable Energy Laboratory (NREL) Technical report NREL/TP-560-34302, Revised January 2008.

This document is avalaible at http://rredc.nrel.gov/solar/codesandalgorithms/spa/ and is included in the .zip file

This algorithm is based on numerical approximation of the exact equations. The authors of the original paper state that this algorithm should be precise at +/- 0.0003 degrees. I have compared it to NOAA solar table (http://www.srrb.noaa.gov/highlights/sunrise/azel.html) and to USNO solar table (http://aa.usno.navy.mil/data/docs/AltAz.html) and found very good correspondance (up to the precision of those tables), except for large zenith angle, where the refraction by the atmosphere is significant (difference of about 1 degree). Note that in this code the correction for refraction in the atmosphere has been implemented for a temperature of 10C (283 kelvins) and a pressure of 1010 mbar. See the subfunction «sun_topocentric_zenith_angle_calculation» for possible modification to explicitely model the effects of temperature and pressure as describe in Reda & Andreas (2003).

Input parameters:
time: a structure that specify the time when the sun position is calculated.
time.year: year. Valid for [-2000, 6000]
time.month: month [1-12]
time.day: calendar day [1-31]
time.hour: local hour [0-23]
time.min: minute [0-59]
time.sec: second [0-59]
time.UTC: offset hour from UTC. Local time = Greenwich time + time.UTC
This input can also be passed using the Matlab time format ('dd-mmm-yyyy HH:MM:SS'). In that case, the time has to be specified as UTC time (time.UTC = 0)

location: a structure that specify the location of the observer
location.latitude: latitude (in degrees, north of equator is positive)
location.longitude: longitude (in degrees, positive for east of Greenwich) location.altitude: altitude above mean sea level (in meters)

Output parameters
sun: a structure with the calculated sun position
sun.zenith = zenith angle in degrees (angle from the vertical)
sun.azimuth = azimuth angle in degrees, eastward from the north.

Only the sun zenith and azimuth angles are returned as output, but a lot of other parameters are calculated that could also be extracted as output of this function. See the documentation in the code.

Exemple of use

location.longitude = -105.1786;
location.latitude = 39.742476;
location.altitude = 1830.14;
time.year = 2003;
time.month = 10;
time.day = 17;
time.hour = 12;
time.min = 30;
time.sec = 30;
time.UTC = -7;

sun = sun_position(time, location);

sun =

zenith: 50.1080438859849
azimuth: 194.341174010338

Acknowledgements

This file inspired Sun Azimuth Data, Mspa, and Get Geometric Location And Other Information.

MATLAB release MATLAB 6.5 (R13)
06 Mar 2014
10 Feb 2014

Hi, I have a quick question. Is there a way to use this function in reverse? I have a known lat/lon position, the year, the day, and the 3 solar elevation angles. Based on those givens, I am trying to determine the time of day. Is that possible with this?
Thanks!
Andy

10 Jan 2014
05 Nov 2013

Hello, I have a problem! Can you help me with it?

///sun_position
Error using sun_position (line 96)
Not enough input arguments.

25 Oct 2013

Hi Nasir,
I (= Jeroen) am just a user of this sun_position function, perhaps your question is meant for Vincent Roy (=author of sun_position).
Still, I'd be glad to help if I can.

Can you post the code that generates your error message ?

23 Oct 2013

Hi Jeroen,
It's giving me error
"Undefined function 'sun_position' for input arguments of type 'struct'.
"
Could you plz check it
thanks,
Nasir

26 Apr 2013
11 Apr 2013

I am attempting to calculate the clear sky solar irradiance for a given latitude with the equation:

I = I0(cozZ)*Tr*Tpg*Tw*Ta

where cozZ is the cosine of the solar zenith angle. which is given by
cosZ = sin(lat)*sin(dec) + cos(lat)*cos(dec)*cos(H),
where H is the hour angle given by:
H = (pi/12)*(tnoon - t)
where t is the local solar time.

So, the zenith returned from this function can be used as Z in the first equation shown here?

16 Nov 2012

the paper reference to a GUI interface, but it's not. How can I get it?

26 Oct 2012

As unlikely as it is, it is the second time I needed this function

26 Oct 2012
10 Jul 2012

Dan,

You can set the altitude with the location.altitude parameter (in units of meter above sea level)

Vince

28 Jun 2012

Any thoughts on how I can add an altitude to this calculation?

(I need to do this for a balloon-based instrument between 0 and 30 km in altitude..)

28 May 2012

Anandhaprabu,

Don't forget to set your UTC time difference (time.UTC).

Vince

20 Apr 2012

If i run the program the error will come like this

??? Reference to non-existent field 'UTC'.Error in ==> sun_position>julian_calculation at 179ut_time = ((time.hour - time.UTC)/24) + (time.min/(60*24)) +(time.sec/(60*60*24)); % time of day in UT time.Error in ==> sun_position at 96julian = julian_calculation(time);.

let me know what to do.

07 Mar 2012

Dear Jeroen,

Thanks, that's excellent, works perfectly!!
I would have never have thought of that either,

Cheers
Cat

06 Mar 2012

@Cat, @Robin: I hope this example will help. One-hour steps, for the first 5 days of the current month.

Regards - Jeroen

location.longitude = 4.3;
location.latitude = 51.6;
location.altitude = 0;

time.year = 2008;
time.month = 5;
time.day = 16;
time.hour = 16;
time.min = 35;
time.sec = 0;
time.UTC = 1;

sun=sun_position(time, location); %One moment, one place.

StepSize=datenum(0,0,0,1,0,0); %one hour
DN_ar=datenum(2012,03,01):StepSize:datenum(2012,3,5);

for i=1:length(DN_ar)
dn=DN_ar(i);
[time.year time.month time.day time.hour time.min time.sec] = datevec(dn);
sun=sun_position(time, location);
Z(i)=sun.zenith;
A(i)=sun.azimuth;
end
plot(Z,'b');
hold on;
plot(A,'r');

06 Mar 2012

Hi, I am asking a similar question to Robin, in that I would like to calculate the azimuth and zenith for every hour of 100 days for the same place. I have made vectors for each of the location and time fields and now have the location and time structures, but the model appears to accept scalars only as an input. I was going to use a for loop to get around this, but I don't know how to get the for loop to run through the vectors in the function. I have tried this so far:
for t=1:1:1896
sun(t)=sun_position(time(t),location(t))
end
??? Undefined function or method 'sun_position' for input
arguments of type 'struct'. I am sure this is really simple but I have never used structures before, and the for loop help doesn't mention them.

Many thanks,
Cat

16 Nov 2011

Hi Robin,
If those matrices are to be multiplied, you have to try dimensions in the order of 2X3*3X2.
Matrix dimensions are No of rows X(by) No of columns.
Try to match No of columns of first matrix with No of rows of the second. The the product will be a matrix of No of Rows of first BY No of columns of second: a 5 BY 7 times matrix CAN be multiplied with a 7 BY 2 and you'll get a 5 BY 2 matrix.
You can also try the dot product by adding a dot in front of the operator as in .*
Hope this helps.

16 Nov 2011

Hello,
Thanks for a great script! My question is whether this will also work to run serial zenith angles for every hour of a day for n days? I am a pretty new Matlab user so my skills are not great yet. To do this will I have to modify the script much? I tried running it once with time.day = 1:5 and time.hour = 0-23 but there is a "matrix dimensions must agree error". Before I go on a hunt to try and resolve it I thought I would check to see if anyone else has done this?

Thanks,
Robin

17 Jun 2011

comment telecharger le fichier .m

08 May 2010

Hi folks,

I was thinking about the rigid approach of Reda and Andreas in respect to delta_t, or the difference between UT and TT < http://en.wikipedia.org/wiki/%CE%94T >

I did a error check for the period between 2010 and 2020 by adopting the proposed in SPA's documentation 64.797s and by a fourier expression of delta_t which I derived using data generated by the HORIZON's software at the JPL(UNI PASADENA) < http://ssd.jpl.nasa.gov/horizons.cgi#top >.

The function generated is inherently periodic and describes the the ephemerides data with a high precision ( st.dev = 1.5e-5).

I am currently analyzing the effects of this deviation on SPA's precision and a noticeable improvement can be observed in the e-5 order of magnitude for azimuth and zenith.

The payback of implementing the function is great, compared to the minor increase in computational cost.

After line number 203-204 of the sun_position function implement the following:
*******************************************************************************
% Compute serial date number since 2010-01-01 00:00:00
JD2010 = julian.day - 2455196.5;
% Compute delta_t using fourier expression derived from HORIZON ephemeris
% PRECISION: rmse = 9.811e-006s; std = 1.5e-5s;
delta_t = 66.18 - 0.0001033*cos(JD2010*0.0172) + 0.001654*sin(JD2010*0.0172);
*********************************************************************************

Regards,
Lubo

16 Jan 2010
15 Jul 2009

23 Jan 2009

Great tool! Good job! However, I wonder, if there will be an update on this script, after reading Jeff Gullett's remarks on the code?

06 Jan 2009

I have found a few other discrepancies in the algorithm. In the 2005 version of the document, the algorithm to compute the topocentric sun declination (eqn #39, line #825) contains the variable 'y' in the denominator, but the code presented later in the pdf document has the variable 'x' in its place. In the 2008 version of the algorithm document, this difference is resolved and the variable should be 'x' (as in the C-code, so it needs to be changed in the m-code presented on this page). Also, when computing the mean obliquity of the ecliptic (eqn #24, line #733), the documented steps have one of the terms as "-4680.93", but again, the C-code from the appendix used the term "-4680.96". Unfortunately, this discrepancy still exists in the 2008 document, so I don't know which term is correct. Using the wrong one will result in an error greater than the maximum error claimed by the algorithm designers.

19 Nov 2008

Congrats guy!!

11 Nov 2008

Has anyone done this same script but in excel? I've been trying it and I still have some formula wrong... I guess I'm confuse with the t ut or the radians and degrees...
Would anyone want to give it a check?

24 Oct 2008

The algorithm has a flaw... The atmospheric refraction correction has a number of vertical asymptotes below -5 degrees true elevation to the sun. For elevation angles that hit these asymptotes, I was getting points that resulted in positive elevations, when the neighboring points would come out with the expected result of "sun well below the horizon". You need to add a conditional prior to applying the atmospheric correction to only add this in if the true elevation is greater than -5 degrees.

As for the code itself, it is somewhat disorganized and difficult to read, but it does work. I had a specific application where I needed to improved the efficiency/readability of the code, so I ended up rewriting the thing from scratch based on the paper, except I store everything in radians instead of degrees. This is still a nice application, though.

17 Oct 2008

Thank you for your perfect program!!!

24 Oct 2007

This m-file is great! Very useful, well written.

03 Aug 2007

Thanks for the program.....The reference from usno.navy is very well accepted and well known for scientist and engineer

08 Feb 2007

Well report

22 Jan 2007

The code works as advertise and is well documented. It even comes with a DOE report on the equations involved. The learning curve is extremely short. Thanks!

29 May 2006
15 Feb 2006
18 Nov 2005
28 Aug 2005

Very high-quality. I have a translation to Delphi

27 Aug 2005

Cool!

10 Jul 2005

Bingo! This is exactly what I am looking for.

07 Oct 2004

It is a grate work and very halpful for me.Thanks a lot.

17 Mar 2004

Very useful script!

15 Mar 2004

Small bugs fix (in julian day calculation)

18 Mar 2004

Small bugs fix (in julian day calculation)

18 Mar 2004

Changed information in the header (help display). location.elevation should have been location.altitude.

13 Apr 2004

Following a suggestion from Jody Klymak (jklymak@ucsd.edu), allowed the 'time' input to be passed as a Matlab time string.

23 Feb 2009

Thanks to the sharp eyes of Jeff Gullett, corrected eq. 39 as per updated Resa and Andreas paper (see line 825) and added a condition for the refractive correction (see lines 854-858, now only applied for sun elevation greater than -5 deg).