I am writing a code to plot an orbit with considering the effects of the earths oblateness and I am having trouble getting it to work properly. I am not sure the problem exactly, but the x coordinates seem to jump at around the 3000 second mark. Also, when plotted in 3 dimensions there are lines that are bisecting the orbit indicating that the x or y values are not continues however plotted in 2 dimensions they look continues. I am fairly confused with what I did so any help would be appreciated. The code is below!
% Orbit calculator for j2 effects
% Ryan Goodridge
% 101147736
% ryangoodridge@cmail.carleton.ca
% Date performed: 03/10/2023
clc
close all;
clear all;
% Useful constants
Re = 6378; % radius of earth
u = 398600; % MU of earth it goes moo moo
d2r = pi/180; % rad/deg
J2 = 0.0010836;
% user defined orbital elements
a = Re + 750; % km
e = 0.0000001;
i = 80*d2r; % rad
tp = 0; % sec
RAAN = 270*d2r; % rad
w = 0*d2r; % rad
% Calculating the Orbit radius vector, velocity magnitude, and time
Orbit1 = Orbitcalc(a,e,i,w,RAAN,u,J2,Re); % Using orbit function to calculate the position of a spacecraft relitive to the time
r_P = [Orbit1(1,:); Orbit1(2,:); Orbit1(3,:)]; % retreving data from function
% first ECI point
C_IP = C3(Orbit1(7,1))*C1(Orbit1(8,1))*C3(Orbit1(9,1));
r_I = C_IP*r_P(:,1);
% creating a rotation matrix for perifocal to ECI
sz = size(Orbit1);
sz = sz(1,2);
for j = 2:1:sz
C_IP = C3(Orbit1(7,j))*C1(Orbit1(8,j))*C3(Orbit1(9,j));
r_I(:,end+1) = C_IP*r_P(:,j);
end
R = sqrt(Orbit1(1,:).^2 + Orbit1(2,:).^2);
figure
plot(Orbit1(10,:),R)
hold on
plot(Orbit1(10,:),Orbit1(11,:))
hold off
% Plot of orbit in parifocal
figure
plot(Orbit1(1,:),Orbit1(2,:))
hold on
plot(0,0,'g--o')
title('Orbit in Parifocal');
ylabel('Distance [m]');
xlabel('Distance [m]');
% testing to see what is wrong
figure
plot(Orbit1(10,:),Orbit1(1,:))
title('x v time')
figure
plot(Orbit1(10,:),Orbit1(2,:))
title('y v time')
% Plot of orbit in ECI
figure
plot3(r_I(1,:),r_I(2,:),r_I(3,:),'k','linewidth', 2)
grid on
hold on
% [xs,ys,zs]=sphere(50);
% xs = xs*Re;
% ys = ys*Re;
% zs = zs*Re;
% obj = surf(xs,ys,zs);
set(gca,'FontSize',9,'FontName', 'Times')
xlabel('rx_I (km)')
ylabel('ry_I (km)')
zlabel('rz_I (km)')
axis equal
% Orbit Propogator
function Orbitdata = Orbitcalc(a,e,i,w,RAAN,mu,J2,Re)
% initial calculations
T = (2*pi*a^(3/2))/(sqrt(mu)); % Period of orbit
Rp = a*(1 - e); % Radius of parigee
Ra = a + a*e; % Radius of Apsis
R = Ra; % inital position magnitude
% inital position components
R_x = Rp;
R_y = 0;
tol = 10^-8; %rad %% Tolerance of the orbit
V = sqrt(2*mu/R - mu/a); % Inital velocity
for t = 0:0.1:T-0.1
% inital calculations
p = a(end)*sqrt(1 - e(end)^2);
% Orbit calculations
M = 2*pi*t/T; % mean anomaly
E = CalcEA(M,e(end),tol); % eccentric anomaly
theta = 2*atan(sqrt((1+e(end))/(1-e(end)))*tan(E/2)); % true anomaly
R_mag = a(end)*(1-e(end)^2)/(1+e(end)*cos(theta)); % radius of orbit
% Orbital Perturbations
u = w(end) + theta;
f_r = (-3*mu*J2*(Re^2)*(1 - 3*(sin(i(1,end))*sin(i(1,end)))*(sin(u)*sin(u))))/(2*R_mag^4); % Eq 3.66
f_theta = ((-3*u*J2*(Re^2))/(2*(R_mag)^4))*((sin(i(1,end))*sin(i(1,end)))* (sin(2*u))); % Eq 3.67
f_Z = ((-3*u*J2*(Re^2))/(2*(R_mag)^4))*((sin(2*i(end)))*(sin(u))); % Eq 3.68
a_dot = ((2*a(end)^2)/(sqrt(p*mu)))*(e(end)*sin(theta)*f_r + (1 + e(end)*cos(theta))*f_theta); % Eq 3.15
e_dot = sqrt(p/mu)*(sin(theta)*f_r + (2*cos(theta) + e(end)*(1 + (cos(theta))^2))*f_theta/(1 + cos(theta))); % Eq 3.16
i_dot = sqrt(p/mu)*(sin(u)*f_Z/(1 + e(end)*cos(theta))); % eq 3.17
RAAN_dot = sqrt(p/mu)*(sin(u)*f_Z/(sin(i(end))*(1 + e(end)*cos(theta)))); % Eq 3.18
w_dot = sqrt(p/mu)*(-cos(theta)*f_r/e(end) + (2 + e(end)*cos(theta)*sin(theta))*f_theta/(e(end)*(1 + e(end)*cos(theta))) - sin(u)*f_Z/((1 + e(end)*cos(theta))*tan(i(end)))); % Eq 3.19
% Store Values
a(end+1) = a(end) + a_dot*0.1; % new values of a
e(end+1) = e(end) + e_dot*0.1; % new values of e
i(end+1) = i(end) + i_dot*0.1; % new values of i
RAAN(end+1) = RAAN(end) + RAAN_dot*0.1; % new values of RAAN
w(end+1) = w(end) + w_dot*0.1; % new values of w
R(end+1) = R_mag; % magnitude of position vector
R_x(end+1) = R_mag*cos(theta); % x compoent of radius
R_y(end+1) = R_mag*sin(theta); % y compoent of radius
V(end+1) = sqrt(2*((-mu/(2*a(end)))+(mu/R_mag))); % velocity of spacecraft
end
time = 0:0.1:T; % time
sz = size(R_x);
Z = zeros(sz); % z compoent of radius
Orbitdata = [R_x; R_y; Z; V; a; e; i; RAAN; w; time; R]; % output
end
% rotation matrix
function A = C1(x)
if nargin > 1
error('Too many inputs. See help file')
end
A = [ 1 0 0;
0 cos(x) -sin(x);
0 sin(x) cos(x)];
end
% rotation matrix
function A = C2(x)
if nargin > 1
error('Too many inputs. See help file')
end
A = [cos(x) 0 -sin(x);
0 1 0;
sin(x) 0 cos(x)];
end
% rotation matrix
function A = C3(x)
if nargin > 1
error('Too many inputs. See help file')
end
A = [cos(x) -sin(x) 0;
sin(x) cos(x) 0;
0 0 1];
end
% Orbit eccentric anomaly, Kepler's equation keplers equation
% Richard Rieber
% 1/23/2005
% rrieber@gmail.com
%
% Revision 8/21/07: Fixed typo in line 38 if statement
% Added H1 line for lookfor functionality
%
% function E = CalcEA(M,ecc,tol)
%
% Purpose: Solves for eccentric anomaly, E, from a given mean anomaly, M,
% and eccentricty, ecc. Performs a simple Newton-Raphson iteration
%
% Inputs: o M - Mean anomaly in radians.
% o ecc - Eccentricity of the orbit.
% o tol - a tolerance at which to terminate iterations; Default
% is 10^-8 radians. [OPTIONAL]
%
% Outputs: o E - The eccentric anomaly in radians.
%
% E = CalcEA(M,ecc) uses default tolerances
%
% E = CalcEA(M,ecc,tol) will use a user specified tolerance, tol
%
function E = CalcEA(M,ecc,tol)
%Checking for user inputed tolerance
if nargin == 2
%using default value
tol = 10^-8;
elseif nargin > 3
error('Too many inputs. See help CalcEA')
elseif nargin < 2
error('Too few inputs. See help CalcEA')
end
if (M > -pi && M < 0) || M > pi
E = M - ecc;
else
E = M + ecc;
end
Etemp = E + (M - E + ecc*sin(E))/(1-ecc*cos(E));
Etemp2 = E;
while abs(Etemp - Etemp2) > tol
Etemp = Etemp2;
Etemp2 = Etemp + (M - Etemp + ecc*sin(Etemp))/(1-ecc*cos(Etemp));
end
E = Etemp2;
end
