I have created my own J2 Orbit Propagator and I am not sure why it is not working properly.

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

% [e-mail]

% 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

% [e-mail]

%

% 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

Answer 1

[Ashok]

Hey @Ryan Goodridge,

Debugging is tricky without knowing the source of the 'Orbit Perturbation' equations mentioned in the code. However, applying a wrap function to the 'u' relation seems to resolve the discontinuity when using an integration time step of 0.1.

u = wrapTo2Pi(w(end) + theta);

Starting from MATLAB R2024a, the Satellite Communication Toolbox allows for plotting a satellite's orbit with the J2 effect. More details can be found here:

https://www.mathworks.com/help/releases/R2024a/aerotbx/ug/hpop-algorithms.html#:~:text=Oblate%20ellipsoid%20(J2).

Here's a snippet that shows how to plot a satellite's orbit including the J2 effect.

startTime = datetime(2020,6,02,8,23,0);
stopTime = startTime + hours(5);
sampleTime = 60;
sc = satelliteScenario(startTime,stopTime,sampleTime,"AutoSimulate", true);
options = numericalPropagator(sc, "GravitationalPotentialModel", "oblate-ellipsoid");
sat = satellite(sc,"threeSatelliteConstellation_1.tle");
ele1 = orbitalElements(sat);
time = datetime(2020,6,02,12,30,0);
pos = states(sat,"CoordinateFrame","inertial");
plot3(pos(1,:), pos(2,:), pos(3,:))
axis equal
grid on
show(sat)

The ‘numericalPropagator’ object is used to set the ‘GravitationalPotentialModel’ to ‘oblate-ellipsoid’, instructing the Satellite Scenario to account for the J2 effect during simulation. For more information, check out this documentation page:

https://www.mathworks.com/help/releases/R2024a/aerotbx/ug/satellitescenario.numericalpropagator.html.

I believe this will assist you!

[Ryan Goodridge]

Hello @Ashok,

I apologize for not closing this question sooner.

The equations in question are sourced from the textbook Orbital Mechanics for Engineering Students by Howard D. Curtis. However, some issues arose during numerical integration due to the increasing complexity of the equations.

The solution turned out to be quite straightforward. Although I can't recall the exact source of the approach, I resolved the issue by separating the rate of change equations into a dedicated function outside the loop, then performed the integration within the loop rather than trying to do it all within the loop.