%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Name: draworbit.m
% Author: J. Kristoffer Hultgren
%
% Description:
%
% Draws the orbit with radius r (m), inclination 
% inc (rad) and right ascension omega (rad), returns the x,y,z-coordinates
% as a matrix M, and exports the matrix M as a space separated textfile
% named output.
%
% Input:  r     - Radius of orbit (m)
%         inc   - Inclination of orbit (wrt the Sun line)
%         omega - Argument of right ascension
%
% Output: [r inc omega]
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function elem = draworbit(r,inc,omega)

elem = [r inc omega];
v = 0:2*pi/359:2*pi;
x = r.*cos(v);
y = r.*sin(v);
z = zeros(1,length(x));

% Defines the rotation matrix
R = [cos(omega) -sin(omega)*cos(inc) sin(omega)*sin(inc); 
        sin(omega) cos(omega)*cos(inc) -cos(omega)*sin(inc); 
        0 sin(inc) cos(inc);];
    
% Performs the rotation
Result = R*[x;y;z];

x2 = Result(1,:);
y2 = Result(2,:);
z2 = Result(3,:);

plot3(x2,y2,z2)
hold on

%Plot starting point
a = linspace(0,0,360);
xS = r.*cos(a);
yS = r.*sin(a);
P = R*[xS;yS;z];
xxS = P(1,:);
yyS = P(2,:);
zzS = P(3,:);
plot3(xxS,yyS,zzS,'h','Color','r')

axis([-r r -r r -r r])
grid
box
xlabel('x - Sunwards')
ylabel('y')
zlabel('z - Ecliptic North')

% Writes the resulting matrix to a text file
M = [x2' y2' z2'];
dlmwrite('orbit.dat',M,' ')

% Inserts text at the beginning of the text file
!sed -i 1i"#START" orbit.dat

