Sun Synchronous with For Loop

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

Answer 1

[Alan Stevens]

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

[Daynah Rodriguez]

yes, the code runs fine. But I want to modify it using a for loop with matlab first. GMAT does not interpret Alt = 250:5:1000

[Alan Stevens]

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;