MATLAB Answers

Sun Synchronous with For Loop

4 views (last 30 days)
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

Accepted Answer

Alan Stevens
Alan Stevens on 30 May 2021
When I run the program you've listed above it produces the graph you show!
  2 Comments
Alan Stevens
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;

Sign in to comment.

More Answers (0)

Community Treasure Hunt

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

Start Hunting!