function [S,GST] = sunGEI(yyyy,ddd,seconds)
% function [S,GST] = sunGEI(yyyy,ddd,seconds)
% returns the unit vector position of the sun in GEI coordinates
% It also returns grenwich mean sidereal time (GST) in radians
% yyyy,ddd,seconds are integers in UT (GMT)
% alternate calling: S = sunGEI(decimalDOY)
% valid for yyyy = 1901-2099

if nargin == 1,
UT = yyyy;
yyyy = 1999;
ddd = floor(UT);
seconds = floor(24*60*60*(UT-ddd));
end
rad = 57.29578;

if (yyyy < 1901) | (yyyy > 2099),
error(sprintf('sunGEI: Year %d out of range',yyyy));
end

fday = seconds/86400;
dj = 365*(yyyy-1900)+(yyyy-1901)/4 + ddd + fday -0.5;
t = dj/36525;
vl = mod(279.696678 + 0.9856473354*dj,360.0);
GST = mod(279.696678 + 0.9856473354*dj+360*fday+180,360.0)/rad;
g = mod(358.475845 + 0.985600267*dj,360.0)/rad;
slong = vl+(1.91946-0.004789*t)*sin(g) + 0.020094*sin(2*g);
obliq = (23.45229-0.0130125*t)/rad;
slp = (slong-0.005686)/rad;
sind = sin(obliq)*sin(slp);
cosd = sqrt(1-sind^2);
SDEC = rad*atan(sind/cosd);
SRASN = 180-rad*atan2(1/tan(obliq)*sind/cosd,-cos(slp)/cosd);

X = cos(SRASN/rad)*cos(SDEC/rad);
Y = sin(SRASN/rad)*cos(SDEC/rad);
Z = sin(SDEC/rad);
S = [X;Y;Z];