Code covered by the BSD License

### Highlights from sun_position.m

4.75
4.8 | 24 ratings Rate this file 238 Downloads (last 30 days) File Size: 1.27 MB File ID: #4605 Version: 1.1

# sun_position.m

### Vincent Roy (view profile)

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
MATLAB release MATLAB 6.5 (R13)
24 Jul 2015 Aziz KHIDACH

### Aziz KHIDACH (view profile)

Hi guys !
i need some help here !
i download the program from here , and when i try to run it it give me this message :
"" Undefined function 'sun_position' for input arguments of type 'struct' ""

Thanks

Comment only
17 Apr 2015 studentU

### studentU (view profile)

xxw

Comment only
09 Mar 2015 John Wood

### John Wood (view profile)

Reda & Andreas have issued a set of corrections for their algorithm at http://www.sciencedirect.com/science/article/pii/S0038092X07000059

Comment only
27 Feb 2015 Eric

### Eric (view profile)

It seems there's a slight bug in the implementation for the calculation of the Julian Day. The reference describes the INT operation in Equation 4 as the Integer of the calculation terms. The authors state that INT(-8.7) = -8.

However, the paper implements this with the floor function, but floor(-8.7) = -9. I believe the fix function in Matlab should be used instead. This will only affect years prior to -4716.

Comment only
28 Nov 2014 Clara

### Clara (view profile)

My program is telling me that I have an error in line 47, this code.

Undefined function 'sun_position_gha_dec' for input arguments of type 'struct'.

Error in plot_orbit_v3 (line 47)
sunp = sun_position_gha_dec(time);
Line 47:
sunp = sun_position_gha_dec(time);

table = [
fix(sunp.gha) abs(sunp.gha-fix(sunp.gha))*60;
fix(sunp.dec) abs(sunp.dec-fix(sunp.dec))*60]';

Comment only
06 Mar 2014 David Young

### David Young (view profile)

10 Feb 2014 Ander

### Ander (view profile)

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

Comment only
10 Jan 2014 Willem

### Willem (view profile)

05 Nov 2013 Xander

### Xander (view profile)

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

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

Comment only
25 Oct 2013 jhh

### jhh (view profile)

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 ?

Comment only

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

Comment only
26 Apr 2013 Leah

### Leah (view profile)

11 Apr 2013 Richard

### Richard (view profile)

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?

Comment only
16 Nov 2012 marilia15

### marilia15 (view profile)

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

26 Oct 2012 Jaroslaw Tuszynski

### Jaroslaw Tuszynski (view profile)

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

Comment only
26 Oct 2012 Jaroslaw Tuszynski

### Jaroslaw Tuszynski (view profile)

10 Jul 2012 Vincent Roy

### Vincent Roy (view profile)

Dan,

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

Vince

Comment only
28 Jun 2012 Dan Weaver

### Dan Weaver (view profile)

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 Vincent Roy

### Vincent Roy (view profile)

Anandhaprabu,

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

Vince

Comment only
20 Apr 2012 Anandhaprabu Dandabany

### Anandhaprabu Dandabany (view profile)

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.

Comment only
07 Mar 2012 Catriona

### Catriona (view profile)

Dear Jeroen,

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

Cheers
Cat

Comment only
06 Mar 2012 jhh

### jhh (view profile)

@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');

Comment only
06 Mar 2012 Catriona

### Catriona (view profile)

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

Comment only
16 Nov 2011 Rafael

### Rafael (view profile)

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.

Comment only
16 Nov 2011 Robin

### Robin (view profile)

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

Comment only
17 Jun 2011 automat

### automat (view profile)

comment telecharger le fichier .m

Comment only
08 May 2010 Lyuboslav

### Lyuboslav (view profile)

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

Comment only
16 Jan 2010 Lyuboslav

15 Jul 2009 K E

### K E (view profile)

23 Jan 2009 MathDelight

### MathDelight (view profile)

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 Jeff Gullett

### Jeff Gullett (view profile)

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.

Comment only
19 Nov 2008 EDGAR RAMOS

### EDGAR RAMOS (view profile)

Congrats guy!!

11 Nov 2008 Lu Salinas

### Lu Salinas (view profile)

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?

Comment only
24 Oct 2008 Jeff Gullett

### Jeff Gullett (view profile)

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 Lin Xinan

Thank you for your perfect program!!!

24 Oct 2007 P Johnstone

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

03 Aug 2007 Andora Suwisto

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

08 Feb 2007 Son Nguyen

Well report

22 Jan 2007 George Young

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!

15 Feb 2006 DUNG TRAN THANH
18 Nov 2005 Steve Wales
28 Aug 2005 frederik Steensen

Very high-quality. I have a translation to Delphi

27 Aug 2005 Mat Saenger

Cool!

10 Jul 2005 Xin Jin

Bingo! This is exactly what I am looking for.

07 Oct 2004 lu shengwang

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

17 Mar 2004 Andrea Bergami

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 1.1

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).