function [R, Rt] = geo2gsm(dn)
% Subroutine to calculate the transformation matrix from GEO to GSM coordinates (and back)
%	Created October 1, 2000 (In Japan)
%	function [R, Rt] = geo2gsm(dn);
%		dn		=		datenum of event
%		R		=		the forward matrix (old coord is GEO and new is GSM)
%		Rt		=		the inverse transformation (from GSM to GEO)
%					(Both matrices are returned as 9-col arrays at the times dn)

%		S(GEO)	=		the Sun vector in GEO coordinates [sungei,sungeo] = sun(iyr,iday,dut)
%		mu(1971)=		The north end of the magnetic dipole axis mu = [.06859,-.18602,.98015]

%		Get the three unit vectors of GSM coordinates. R is a matrix with rows [Xgsm;Ygsm;Zgsm]
%		Note that mu is contained in X-Z plane of GSM

% DEFINE THE NORTH DIPOLE
mu 	= [.06859,-.18602,.98015];					%GEO coordinates of north dipole

% 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;					%Decimal Hour of day

% GET THE UNIT SUN VECTOR IN GEO COORDINATES (Xgsm = sungeo)
[sungei,sungeo,gst] = sun(iyr,iday,dut);		%These are 3-col arrays

% CREATE THE Ygsm VECTOR PERPENDIULAR TO BOTH sungeo and mu (then normalized)
Yx		=   mu(2)*sungeo(:,3)-mu(3)*sungeo(:,2);%Cross product elements
Yy		= -[mu(1)*sungeo(:,3)-mu(3)*sungeo(:,1)];
Yz		=   mu(1)*sungeo(:,2)-mu(2)*sungeo(:,1);
Y		= sqrt(Yx.^2+Yy.^2+Yz.^2);								%Length of Y
Yx		= Yx./Y;								%Make a unit vector
Yy		= Yy./Y;
Yz		= Yz./Y;

% CREATE THE Z COMPONENT (ANOTHER CROSS PRODUCT OF TWO ORTHOGONAL UNIT VECTORS) Z =  X x Y
Zx		=   sungeo(:,2).*Yz-sungeo(:,3).*Yy;
Zy		= -[sungeo(:,1).*Yz-sungeo(:,3).*Yx];
Zz		=   sungeo(:,1).*Yy-sungeo(:,2).*Yx;

% CREATE THE FORWARD TRANSFORMATION FROM GEO TO GSM
R 		= [sungeo(:,1), sungeo(:,2), sungeo(:,3), Yx, Yy, Yz, Zx, Zy, Zz];
Rt		= [sungeo(:,1), Yx, Zx, sungeo(:,2), Yy, Zy, sungeo(:,3), Yz, Zz];