function [R,Rt] = gei2geo(dn)
% Calculate the transformation from GEI (geocentric equatorial inertial) to GEO (geographic) coordinates
%	Created Oct 1, 2000 (Japan)
%		dn		=		MATLAB datenum
%		R		=		Forward transformation from gei to geo (9-col array with rows for each time dn)
%		Rt		=		Inverse transformation

dr	= pi/180;
Z	= zeros(size(dn));
O	= ones(size(dn));

%Time processing
[iyr,m,d,h,m,s]	= datevec(dn);					%Get the elements of the date
dn0				= datenum(iyr,1,1,0,0,0);		%First instant of the year in dn
iday			= fix(dn-dn0+1);						%Day of year
dut				= rem(dn,1)*24;					%Hour of day

% GET gst (Greenwich Mean Sidereal Time) from the SUN routine
[sungei,sungeo,gst] = sun(iyr,iday,dut);						%Sun vector in two systems and gst

% CREATE THE ELEMENTS OF THE FORWARD TRANSFORMATION (GEI TO GEO)
R	= [cos(dr*gst), sin(dr*gst),Z, -sin(dr*gst), cos(dr*gst),Z,Z,Z,O];		%Forward transform from gei to geo

% CREATE THE ELEMENTS OF THE INVERSE TRANSFORMATION (GEO TO GEI)
Rt	= [cos(dr*gst), -sin(dr*gst),Z, sin(dr*gst), cos(dr*gst),Z,Z,Z,O];		%Inverse trans from geo to gei
return
