%% sunrs.m calculates the time of ground sunrise and sunset
%% at the points of interest, using John Titheridges eqns.
%% Note that the slong parameter values -1 corresponds to no
%% night and 1 to no day.
%%
%% Author: Phillip A. Webb
%% Last updated: February 2 2006
%%

function [srise,sset,slong,sundec]=sunrs(latmat,day,month,year)

% latmat = latitude of points of interest
% day = day of point of interest
% month = month of point of interest
% year = year of point of interest
% sunrse = local sunrse (fractional hours)
% sunset = local sunset (fractional hours)
% slong = slong parameter
% sundec = solar declination (rads)

sz=size(latmat);                                                 %size of input

%convert all input to column vectors
latmat=latmat(:);
day=day(:);
month=month(:);
year=year(:);

dim=[31 28 31 30 31 30 31 31 30 31 30 31]';                      %days in month
dpoi=dim(month);                                                 %day in poi of interest month 
lyc=rem(year,400)==0 | (rem(year,4)==0 & rem(year,100)~=0);      %check if leap year
if any(lyc)& month==2                                            %reset leap year February dpoi to 29 
 dpoi(any(lyc)& month==2)=29;
end 
 
amon=month+(day-1)./dpoi;                                        %day of year in months        
sundec=0.41015*sin((amon-3.7)*0.5236);                           %sun declination
slong=-tan(latmat).*tan(sundec);
srise=zeros(size(latmat));                                            %sunrise matrix (0=no night)
srise(slong>=1)=12;                                                   %no day        
srise(slong>-1 & slong<1)=12-acos(slong(slong>-1 & slong<1))*3.82;    %time of ground sunrise

%Resize output to input size
srise=reshape(srise,sz);
sset=24-srise;                                                        %time of ground sunset
slong=reshape(slong,sz);
sundec=reshape(sundec,sz);