function [EL,AZ,SOLDST,DEC,HA] = sunang(YYYY, DAY, HH, lat, long, refraction_flag, varargin)
% sunang - a converted fortran program from Dr. Wiscombe at NASA.
% Calculates the elevation angle, azimuthal angle, and Earth-Sun distance.
% 
% SYNTAX:
%       [EL,AZ,SOLDIST,DEC,HA] = sunang(YYYY,DAY,HH,lat,long)
% 
% INPUTS:
%	YYYY	= 1x1 scalar - Year between 1950 and 2050 
%				(approximations become invalid outside this range)
%	DAY 	= 1x1 scalar - Day of year. Must be between 1 and 366
%	HH 		= 1x1 scalar - (local hour) + (time zone number) + (Daylight
%       		savings Time correction; -1 or 0) where (local hour) range is 0 to 
%       		24, (time zone number) range is -12 to +12, and (Daylight Time 
%       		correction) is -1 if on Daylight Time (summer half of year), 0
%       		otherwise.
%	lat		= 1x1 scalar - 90 to -90 with + in North
%	long	= 1x1 scalar - -180 to 180 with + in East 
%
% OUTPUTS:
%	EL		= 1x1 scalar - Solar elevation angle (0 @ horizon, 90 @ zenith)
%	AZ		= 1x1 scalar - Azimuth angle
%	SOLDST	= 1x1 scalar - Normalized solar distance (D/D_0)
%	DEC		= 1x1 scalar - Declination angle
%	HA		= 1x1 scalar - Solar hour angle (- before noon, + after noon)

% The solar angle/position calculations are based on the sunang program
% from Warren Wiscombe's ftp site: ftp://climate1.gsfc.nasa.gov/wiscombe/Solar_Rad/SunAngles/sunae.f
% 
%    Authors:  Dr. Joe Michalsky (joe@asrc.albany.edu)
%              Dr. Lee Harrison (lee@asrc.albany.edu)
%              Atmospheric Sciences Research Center
%              State University of New York
%              Albany, New York
% 
%    Modified by:  Dr. Warren Wiscombe (wiscombe@climate.gsfc.nasa.gov)
%                  NASA Goddard Space Flight Center
%                  Code 913
%                  Greenbelt, MD 20771
% 
% With reference to:
% 
%      Michalsky, J., 1988: The Astronomical Almanac's algorithm for
%         approximate solar position (1950-2050), Solar Energy 40,
%         227-235 (but the version of this program in the Appendix
%         contains errors and should not be used)
% 
% Converted to Matlab:
%       K. Lapo - 07/12
 
% The below commands can be used to test the code to make sure it is doing
% what you think it is doing.
% -------------------------------------------------------------------------
% NASA Goddard:
% [EL,AZ,SOLDST,DEC,HA] = sunang(1995,352,9.525+5,39,-76.8);
% [EL;AZ;SOLDST;DEC*180/pi;HA*180/pi]
% These values in sunang.m
%     lat = 39.0; long = -76.8; year = 1995; day = 352; zone = 5; lst =
%     9.525
% Should give
%     solar azimuth = 143.20 deg; elev = 18.08 deg; 
%     hour angle = -38.35; declination = -23.38; 
%     air mass = 3.19; solar distance = 0.9840
% -------------------------------------------------------------------------
%  Atmospheric Sci Res Ctr (Albany, NY):
% [EL,AZ,SOLDST,DEC,HA] = sunang(1996,12,14.6+5,42.7,-73.83);
% [EL;AZ;SOLDST;DEC*180/pi;HA*180/pi]
% These values in sunang.m
%     lat = 42.7; long = -73.83; year = 1996; day = 12; zone = 5; lst = 14.5;
% Should give
%     solar azimuth = 216.68 deg; elev = 16.79 deg; 
%     hour angle = 37.99; declination = -21.66; 
%     air mass = 3.43; solar distance = 0.9835


% Initial Checks
if nargin == 6
    if ~islogical(refraction_flag)
        error('The flag for using the refraction calculation must be a logical')
    end
elseif nargin == 5
    refraction_flag = 0;
end

if YYYY < 1950 || YYYY > 2100
    error('bad input variable YEAR')
elseif DAY < 1 || DAY > 366  
     error( 'bad input variable DAY')
elseif HH < -13.0 || HH > 36.0  
     error( 'bad input variable HOUR')
elseif lat < -90.0 || lat > 90.0 
     error( 'bad input variable lat')
elseif long < -180.0 || long > 180.0 
     error( 'bad input variable long')
end
     
RPD = pi / 180;

% current Julian date (actually add 2,400,000 
% for true JD);  LEAP = leap days since 1949;
% 32916.5 is midnite 0 jan 1949 minus 2.4e6
DELTA = YYYY - 1949;
LEAP = DELTA / 4;
JD = 32916.5 + (DELTA*365 + LEAP + DAY) + HH / 24;

% last yr of century not leap yr unless divisible
% by 400 (not executed for the allowed YEAR range,
% but left in so our successors can adapt this for 
% the following 100 years)
if rem(YYYY, 100) == 0 && rem(YYYY, 400) ~= 0
    JD = JD - 1;
end

% ecliptic coordinates
% 51545.0 + 2.4e6 = noon 1 jan 2000
TIME  = JD - 51545.0;

% force mean longitude between 0 and 360 degs
MNLONG = 280.460 + 0.9856474*TIME;
MNLONG = rem(MNLONG, 360);
if MNLONG < 0
    MNLONG = MNLONG + 360;
end

% mean anomaly in radians between 0 and 2*pi
MNANOM = 357.528 + 0.9856003*TIME;
MNANOM = rem(MNANOM, 360);
if MNANOM < 0
    MNANOM = MNANOM + 360;
end
% Convert to radians
MNANOM = MNANOM*RPD;

% ecliptic longitude and obliquity 
% of ecliptic in radians
ECLONG = MNLONG + 1.915*sin(MNANOM) + 0.020*sin(2.*MNANOM);
ECLONG = rem(ECLONG, 360);
if ECLONG < 0
    ECLONG = ECLONG + 360;
end
OBLQEC = 23.439 - 0.0000004*TIME;
ECLONG = ECLONG*RPD;
OBLQEC = OBLQEC*RPD;

% right ascension
NUM = cos(OBLQEC)*sin(ECLONG);
DEN = cos(ECLONG);
RA = atan(NUM/DEN);

% Force right ascension between 0 and 2*pi
if DEN < 0
    RA = RA + pi;
elseif NUM < 0
    RA = RA + 2*pi;
end

% declination
DEC = asin(sin(OBLQEC)*sin(ECLONG));

% Greenwich mean sidereal time in hours
GMST = 6.697375 + 0.0657098242*TIME + HH;

% Hour not changed to sidereal time since 
% 'time' includes the fractional day
GMST = rem(GMST,24);
if GMST < 0
    GMST = GMST + 24;
end
    
% local mean sidereal time in radians
LMST = GMST + long / 15;
LMST = rem(LMST, 24);
if LMST < 0
    LMST = LMST + 24;
end 
LMST = LMST*15.*RPD;

% hour angle in radians between -pi and pi
HA  = LMST - RA;

if HA < -pi
    HA  = HA + 2*pi;
elseif HA > pi   
    HA  = HA - 2*pi;
end

% solar azimuth and elevation
EL  = asin(sin(DEC)*sin(lat*RPD) + cos(DEC)*cos(lat*RPD)*cos(HA));
AZ  = asin(-cos(DEC)*sin(HA)/cos(EL));

% Put azimuth between 0 and 2*pi radians
if sin(DEC)-sin(EL)*sin(lat*RPD) >= 0
    if sin(AZ) < 0
        AZ  = AZ + 2*pi;
    end
else
    AZ  = pi - AZ;
end

% Convert elevation and azimuth to degrees
EL = EL / RPD;
AZ = AZ / RPD;

%   ======== Refraction correction for U.S. Standard Atmos. ==========
%      (assumes elevation in degs) (3.51823=1013.25 mb/288 K)
if refraction_flag
    
    if EL >= 19.225
        REFRAC = 0.00452*3.51823 / tan(EL*RPD);
    elseif EL > -0.766 && EL < 19.225
        REFRAC = 3.51823 * ( 0.1594 + EL*(0.0196 + 0.00002*EL) ) / ( 1. + EL*(0.505 + 0.0845*EL) );
    elseif EL <= -0.766
        REFRAC = 0;
    end
    EL = EL + REFRAC;
    
end
%  ===================================================================

% distance to sun in A.U.
SOLDST = 1.00014 - 0.01671*cos(MNANOM) - 0.00014*cos(2.*MNANOM);

% checks
if EL < -90.0 || EL > 90.0
    error('output argument EL out of range')
end
if AZ < 0.0 || AZ > 360.0
    error('output argument AZ out of range')
end

end
