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];
)