4.8

4.8 | 15 ratings Rate this file 250 downloads (last 30 days) File Size: 875.78 KB File ID: #4605

sun_position.m

by Vincent Roy

 

10 Mar 2004 (Updated 23 Aug 2005)

Sun position given observer time/location.

I am interested in collaboration

Download Now | Watch this File

File Information
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)
Zip File Content  
Other Files sun_position.m,
34302.pdf
Tags for This File  
Everyone's Tags
Tags I've Applied
Add New Tags Please login to tag files.
Comments and Ratings (17)
17 Mar 2004 Andrea Bergami Very useful script!  
07 Oct 2004 lu shengwang It is a grate work and very halpful for me.Thanks a lot.
10 Jul 2005 Xin Jin Bingo! This is exactly what I am looking for.
27 Aug 2005 Mat Saenger Cool!
28 Aug 2005 frederik Steensen Very high-quality. I have a translation to Delphi
18 Nov 2005 Steve Wales  
15 Feb 2006 DUNG TRAN THANH  
29 May 2006 nadeem nazar  
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!
08 Feb 2007 Son Nguyen Well report
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
24 Oct 2007 P Johnstone This m-file is great! Very useful, well written.
17 Oct 2008 Lin Xinan Thank you for your perfect program!!!
24 Oct 2008 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.
11 Nov 2008 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?
19 Nov 2008 EDGAR RAMOS Congrats guy!!
06 Jan 2009 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.
Please login to add a comment or rating.
Updates
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 Aug 2005 Modifications to the julian_calculation function. Bug was tranforming a string (time) to a structure directly. This was accepted until R14SP2. Script should now be compliant with future releases of Matlab...

Public Submission Policy

NOTICE: Any content you submit to MATLAB Central, including personal information, is not subject to the protections which may be afforded information collected under other sections of The MathWorks, Inc. Web site. You are entirely responsible for all content that you upload, post, e-mail, transmit or otherwise make available via MATLAB Central. The MathWorks does not control the content posted by visitors to MATLAB Central and, does not guarantee the accuracy, integrity, or quality of such content. Under no circumstances will The MathWorks be liable in any way for any content not authored by The MathWorks, or any loss or damage of any kind incurred as a result of the use of any content posted, e-mailed, transmitted or otherwise made available via MATLAB Central. Read the complete Disclaimer prior to use.

Contact us at files@mathworks.com