SUN_ANGLES Finds solar conditions. [sza,saa] = sun_angles(mjd,lat,lon) calculates the Sun's zenith and azimuth angles for a given location and time, from Astronomical Almanac (see Michalsky, 1988). Refraction is not considered. See also: mjd2localtime, date2mjd
0001 function [sza,saa] = sun_angles(mjd,lat,lon) 0002 % SUN_ANGLES Finds solar conditions. 0003 % [sza,saa] = sun_angles(mjd,lat,lon) calculates the Sun's zenith and 0004 % azimuth angles for a given location and time, from Astronomical Almanac 0005 % (see Michalsky, 1988). 0006 % 0007 % Refraction is not considered. 0008 % 0009 % See also: mjd2localtime, date2mjd 0010 0011 %---------------------------------------------------------------- 0012 % Craig Haley 20/09/01 0013 % 02-10-05 CSH vectorized 0014 %---------------------------------------------------------------- 0015 0016 DEGREE = pi/180.0; 0017 RADIAN = 180.0/pi; 0018 0019 %calculate position of sun in celestial 0020 [ra,dec] = sun_celest(mjd); 0021 0022 %calculate the Greenwich mean siderial time in hours 0023 gmst = mjd2gmst(mjd); 0024 0025 %calculate local mean siderial time in radians 0026 lmst = gmst+lon/15.0; 0027 lmst = mod(lmst,24.0); 0028 lmst(lmst < 0) = lmst(lmst < 0)+24; 0029 lmst = lmst*15.0*DEGREE; 0030 0031 %calculate hour angle in radians between -PI and PI 0032 %here we ignore the equation of the equinoxes factor 0033 ha = lmst-ra; 0034 ha(ha < -pi) = ha(ha < -pi)+2*pi; 0035 ha(ha > pi) = ha(ha > pi)-2*pi; 0036 0037 %calculate altitude angle in radians 0038 a = asin(sin(dec).*sin(lat*DEGREE)+cos(dec).*cos(ha).*cos(lat*DEGREE)); 0039 0040 %calculate azimuth in radians between 0 and 2*PI 0041 az = atan2(sin(ha),cos(ha).*sin(lat*DEGREE)-tan(dec).*cos(lat*DEGREE))+pi; 0042 0043 % %calculate the refraction correction factor 0044 % P = 1016.0; %pressure in millibars 0045 % T = 23.5; %temperature in C 0046 % if (a < 15.0) 0047 % R = 0.00452*P/((273.0+T)*tan(a)); 0048 % else 0049 % R = P*(0.1594+0.0196*a+0.00002*a*a)/((273.0+T)*(1.0+0.505*a+0.0845*a*a)); 0050 % end 0051 0052 %calculate zenith and azimuth angles in degrees 0053 sza = 90.0-(a*RADIAN); %add R to correct for refraction 0054 saa = az*RADIAN;