%% sundec.m calculates the solar declination angle using Titheridge's equations.
%%

function sdec=sundec(day,month,year)

% day = day of interest
% month = month of interest
% year = year of interest
% sdec = solar declination (rads)

sz=size(day);                                                    %size of input

%convert all input to column vectors
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        
sdec=0.41015*sin((amon-3.7)*0.5236);                             %sun's declination

%Resize output to input size
sdec=reshape(sdec,sz);
