function [sungei,sungeo,gst] = sun(iyr,iday,dut)
% M-FILE TO CREATE A GEOCENTRIC EQUATORIAL INERTIAL AND GEOGRAPHIC SUN VECTOR
% [sungei,sungeo] = sun(iyr,iday,dut);
% Derived from subroutine SUN in Russel paper on coord transform Mods Sept 2000
% iyr = Integer year of desired time
% iday = Integer day of year (Jan 1 = 1)
% dut = Universal time in decimal hours
% sungei = Unit Sun vector as row of a 3-col array for each input time in GEI coords
% sungeo = Unit Sun vector in GEO coords, also row of a 3-column array
% gst = Greenwich Mean Sideral Time (degrees)

degrad = pi/180.;
if(iyr<1901 | iyr>2099)
disp('Year out of range of subroutine SUN')
return
end
fday = dut/24;
dj = 365*(iyr-1900)+fix((iyr-1901)/4)+iday+fday-0.5;
t = dj/36525;
vl = rem(279.696678+0.9856473354*dj,360);
gst = rem(279.690983+0.9856473354*dj+360*fday+180.,360);
g = rem(358.475845+0.985600267*dj,360)*degrad;
slong = vl+(1.91946-0.004789*t).*sin(g)+0.020094*sin(2*g);
obliq = (23.45229-0.0130125*t)*degrad;
slp = (slong-.005686)*degrad;
sind = sin(obliq).*sin(slp);
cosd = sqrt(1-sind.^2);
sdec = atan(sind./cosd)/degrad;
srasn = 180.0-(1/degrad)*atan2((1./tan(obliq)).*sind./cosd,-cos(slp)./cosd);
%
%TRANSFORMATION FROM GEI TO GEO COORDINATES
zcol = zeros(size(iyr)); %Column of zeros
ocol = ones(size(iyr)); %Column of ones
% NOTE: I am writing the transformation matrix [Row1;Row2;Row3] as [Row1, Row2, Row3] {a 9 col array}
%rgei_geo= [cos(degrad*gst),sin(degrad*gst),zcol;-sin(degrad*gst),cos(degrad*gst)zcol;zcol,zcol,ocol];
Row1 = [cos(degrad*gst),sin(degrad*gst),zcol];
Row2 = [-sin(degrad*gst),cos(degrad*gst),zcol];
Row3 = [ zcol,zcol,ocol];

%SUN VECTOR IN GEI COORDINATES
xsun = cos(degrad*srasn).*cos(degrad*sdec);
ysun = sin(degrad*srasn).*cos(degrad*sdec);
zsun = sin(degrad*sdec);
%NOTE: I am writing the col sun vector Rsun as a three col array [xsun, ysun, zsun]
%sungei = [xsun;ysun;zsun]; %As a column vectorf
sungei = [xsun,ysun,zsun];

%TRANSFORM SUN VECTOR TO GEO COORDINATES
% I don't know how to do this with multidimensional arrays so old fashioned element by element
%sungeo = rgei_geo*sungei; This won't work because both (R) and Vect are time dependent
xsungeo = sum((Row1.*sungei)')';
ysungeo = sum((Row2.*sungei)')';
zsungeo = sum((Row3.*sungei)')';
sungeo = [xsungeo,ysungeo,zsungeo];