4 views (last 30 days)

Show older comments

Hello,

I am trying to perform plot a of sun synchronous orbit using GMAT (a mission design software which has a script format). I found this code for matlab:

clc;

clear all;

mu = 398600.440; % Earth’s gravitational parameter [km^3/s^2]

Re = 6378; % Earth radius [km]

J2 = 0.0010826269; % Second zonal gravity harmonic of the Earth

we = 1.99106e-7; % Mean motion of the Earth in its orbit around the Sun [rad/s]

% Input

Alt = 250:5:1000; % Altitude,Low Earth orbit (LEO)

a = Alt + Re; % Mean semimajor axis [km]

e = 0.0; % Eccentricity

h = a*(1 - e^2); % [km]

n = (mu./a.^3).^0.5; % Mean motion [s-1]

tol = 1e-10; % Error tolerance

% Initial guess for the orbital inclination

i0 = 180/pi*acos(-2/3*(h/Re).^2*we./(n*J2));

err = 1e1;

while(err >= tol )

% J2 perturbed mean motion

np = n.*(1 + 1.5*J2*(Re./h).^2.*(1 - e^2)^0.5.*(1 - 3/2*sind(i0).^2));

i = 180/pi*acos(-2/3*(h/Re).^2*we./(np*J2));

err = abs(i - i0);

i0 = i;

end

plot(Alt,i,'.b');

grid on;hold on;

xlabel('Altitude,Low Earth orbit (LEO)');

ylabel('Mean orbital inclination');

title('Sun-Synchronous Circular Orbit,Inclination vs Altitude(LEO,J2 perturbed)');

hold off;

However, GMAT does not interpret Alt = 250:5:1000. And so I've been trying to modify it such as using the for loop on matlab first

for Alt = 250:5:1000

a = Alt + Re; % Mean semimajor axis [km]

e = 0.0; % Eccentricity

h = a*(1 - e^2); % [km]

n = (mu./a.^3).^0.5; % Mean motion [s-1]

end

But it's not working...

This is the graph Im trying to achieve with matlab if the modification is done right and before applying it on GMAT

Thank you

Alan Stevens
on 30 May 2021

When I run the program you've listed above it produces the graph you show!

Alan Stevens
on 30 May 2021

Ok. Like this then:

mu = 398600.440; % Earth’s gravitational parameter [km^3/s^2]

Re = 6378; % Earth radius [km]

J2 = 0.0010826269; % Second zonal gravity harmonic of the Earth

we = 1.99106e-7; % Mean motion of the Earth in its orbit around the Sun [rad/s]

% Input

Alt = 250:5:1000; % Altitude,Low Earth orbit (LEO)

i = zeros(1,numel(Alt)); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for k = 1:numel(Alt) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

a = Alt(k) + Re; % Mean semimajor axis [km] %%%%%%%%%%%%%%%%%%

e = 0.0; % Eccentricity

h = a*(1 - e^2); % [km]

n = (mu./a.^3).^0.5; % Mean motion [s-1]

tol = 1e-10; % Error tolerance

% Initial guess for the orbital inclination

i0 = 180/pi*acos(-2/3*(h/Re).^2*we./(n*J2));

err = 1e1;

while(err >= tol )

% J2 perturbed mean motion

np = n.*(1 + 1.5*J2*(Re./h).^2.*(1 - e^2)^0.5.*(1 - 3/2*sind(i0).^2));

i(k) = 180/pi*acos(-2/3*(h/Re).^2*we./(np*J2)); %%%%%%%%%%%%%%%%%%

err = abs(i(k) - i0); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

i0 = i(k); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

end

end

plot(Alt,i,'.b');

grid on;hold on;

xlabel('Altitude,Low Earth orbit (LEO)');

ylabel('Mean orbital inclination');

title('Sun-Synchronous Circular Orbit,Inclination vs Altitude(LEO,J2 perturbed)');

hold off;

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!