%% soco2.m is a reduced version of my MATLAB version of the soco
%% subroutine from IRI95. It only calculates the solar declination
%% angle (rads).
%%

function dec=soco2(lonmat,timmat,dno)

% lonmat = geographic longitude of points of interest (rads).
% timmat = univerisal time (UT) of points of interest. 
% dno = day number matrix.
%
% dec = solar declination (rads)

lonmat2=rad2deg(longch2(lonmat));                               %convert lonmat value to 0 to 360        
timl=timut2l(timmat,lonmat);                                    %local time

%required constants
p1=0.017203534;
p2=0.034407068;
p3=0.051610602;
p4=0.068814136;
p6=0.103221204;

wlon=360-lonmat2;                                               %convert to WEST longitude

%time of equinox for 1980
tmc=(timl+wlon/15)/24;
tmc(tmc>=1)=tmc(tmc>=1)-1;
td=dno+tmc;                                       
te=td+0.9369;

%declination of the sun
dca=23.256*sin(p1*(te-82.242))+0.381*sin(p2*(te-44.855));
dcb=0.167*sin(p3*(te-23.355))-0.013*sin(p4*(te+11.97));
dcc=0.011*sin(p6*(te-10.41))+0.339137;
dec=deg2rad(dca+dcb+dcc);
