%% orbtrsf.m transforms between the orbital plane
%% and the geographical co-ordinates or vice versa,
%% depending on the what is inputted as one of the
%% parameters. Note that orbtrsf.m converts between
%% cartesian co-ordinates of these two.
%% The orbital plane position is defined
%% relative to the geographical co-ordinates by the
%% inclination of the orbit and the longitude at which
%% the orbit crosses the equator ie by two angular
%% parameters, which orbtrsf.m also requires as input.
%% The inclination and longitude of crosses are given
%% in degrees. Input and output co-ordinates are given
%% as rows.

function [xco,yco,zco]=orbtrsf(xci,yci,zci,incl,lone,tft)

% xci = x co-ordinates in initial co-ordinate system
% yci = y co-ordinates in initial co-ordinate system
% zci = z co-ordinates in initial co-ordinate system
% incl = inclination of orbit in geographic frame (degrees)
% lone = longitude of orbit's asending equatorial crossing (degrees)
% tft = transform type required (1 = geographical to orbit, 2 = orbit to geographic)
% xco = x co-ordinates in final co-ordinate system
% yco = y co-ordinates in final co-ordinate system
% zco = z co-ordinates in final co-ordinate system

incl=deg2rad(incl);           %convert from degrees to rads
lone=deg2rad(lone);           %convert from degrees to rads

% First calculate the elements of the transformation matrix.
% This is based on the theory in Appendix 3 in Introduction
% to Space Physics. What I have done is calculate the axis
% of the orbital co-ordinate system in terms of the geographical
% co-ordinate system. (Phil to do this, you consider each axis
% in the new co-ordinate system one at time. You assume that it
% has unit length ie 1 and calculate that new axis x, y and z component
% values in the old sysytem. Thus x in the new system has components
% (x1,x2,x3) in the old system.

x1=cos(lone);
x2=sin(lone);
x3=0;

y1=-cos(incl)*sin(lone);
y2=cos(incl)*cos(lone);
y3=sin(incl);

z1=sin(incl)*sin(lone);
z2=-sin(incl)*cos(lone);
z3=cos(incl);

trmat=[x1 x2 x3;y1 y2 y3;z1 z2 z3];   %transformation matrix

if tft==1                             %transformation from geographical to orbit
 res=trmat*[xci;yci;zci];             %required results
 xco=res(1,:);                        %x co-ordinate in orbital co-ordinates
 yco=res(2,:);                        %y co-ordinate in orbital co-ordinates
 zco=res(3,:);                        %z co-ordinate in orbital co-ordinates
else
 res=trmat'*[xci;yci;zci];            %required results
 xco=res(1,:);                        %x co-ordinate in geographical co-ordinates
 yco=res(2,:);                        %y co-ordinate in geographical co-ordinates
 zco=res(3,:);                        %z co-ordinate in geographical co-ordinates
end 
