sun_position.m

Sun position given observer time/location.

Description

sun = sun_position(time, location)  
 
This function compute 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.  
 
This document is avalaible at www.osti.gov/bridge 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 as been implemented for a temperature of 10C (283 kelvins) and a pressure of 1010 mbar. See the subfunction «sun_topocentric_zenith_angle_calculation» for a possible modification to explicitely model the effect 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 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 submission has inspired the following:; Sun Azimuth Data
• MATLAB release: MATLAB 6.5 (R13)
• Other Files: sun_position.m, ; 34302.pdf

[Andrea Bergami]

Very useful script!  

[lu shengwang]

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

[Xin Jin]

Bingo! This is exactly what I am looking for.

[Mat Saenger]

Cool!

[frederik Steensen]

Very high-quality. I have a translation to Delphi

[Steve Wales]

 

[DUNG TRAN THANH]

 

[nadeem nazar]

 

[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!

[Son Nguyen]

Well report

[Andora Suwisto]

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

[P Johnstone]

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

[Lin Xinan]

Thank you for your perfect program!!!

[Jeff Gullett]

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.

[Lu Salinas]

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?

[EDGAR RAMOS]

Congrats guy!!

[Jeff Gullett]

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.