how to determine the classical orbital parameter of satellite using matlab code?

21 次查看(过去 30 天)
how to get the orbital parameters of a satellite formation using at low earth orbit?

回答(1 个)

Meysam Mahooti
Meysam Mahooti 2021-5-26
%--------------------------------------------------------------------------
%
% Elements: Computes orbital elements from two given position vectors and
% associated times
%
% Inputs:
% GM Gravitational coefficient
% (gravitational constant * mass of central body)
% Mjd_a Time t_a (Modified Julian Date)
% Mjd_b Time t_b (Modified Julian Date)
% r_a Position vector at time t_a
% r_b Position vector at time t_b
% Outputs:
% Keplerian elements (a,e,i,Omega,omega,M)
% a Semimajor axis
% e Eccentricity
% i Inclination [rad]
% Omega Longitude of the ascending node [rad]
% omega Argument of pericenter [rad]
% M Mean anomaly [rad]
% at time t_a
%
% Notes:
% The function cannot be used with state vectors describing a circular
% or non-inclined orbit.
%
% Last modified: 2018/01/27 M. Mahooti
%
%--------------------------------------------------------------------------
function [a,e,i,Omega,omega,M] = Elements(GM,Mjd_a,Mjd_b,r_a,r_b)
% Calculate vector r_0 (fraction of r_b perpendicular to r_a) and the
% magnitudes of r_a,r_b and r_0
pi2 = 2*pi;
s_a = norm(r_a);
e_a = r_a/s_a;
s_b = norm(r_b);
fac = dot(r_b,e_a);
r_0 = r_b-fac*e_a;
s_0 = norm(r_0);
e_0 = r_0/s_0;
% Inclination and ascending node
W = cross(e_a,e_0);
Omega = atan2(W(1),-W(2)); % Long. ascend. node
Omega = mod(Omega,pi2);
i = atan2(sqrt(W(1)^2+W(2)^2),W(3)); % Inclination
if (i==0)
u = atan2(r_a(2),r_a(1));
else
u = atan2(+e_a(3),(-e_a(1)*W(2)+e_a(2)*W(1)));
end
% Semilatus rectum
tau = sqrt(GM)*86400*abs(Mjd_b-Mjd_a);
eta = FindEta(r_a,r_b,tau);
p = (s_a*s_0*eta/tau)^2;
% Eccentricity, true anomaly and argument of perihelion
cos_dnu = fac/s_b;
sin_dnu = s_0/s_b;
ecos_nu = p/s_a-1;
esin_nu = (ecos_nu*cos_dnu-(p/s_b-1))/sin_dnu;
e = sqrt(ecos_nu^2+esin_nu^2);
nu = atan2(esin_nu,ecos_nu);
omega = mod(u-nu,pi2);
% Perihelion distance, semimajor axis and mean motion
a = p/(1-e^2);
n = sqrt(GM/abs(a^3));
% Mean anomaly and time of perihelion passage
if (e<1)
E = atan2(sqrt((1-e)*(1+e))*esin_nu,ecos_nu+e^2);
M = mod(E-e*sin(E),pi2);
else
sinhH = sqrt((e-1)*(e+1))*esin_nu/(e+e*ecos_nu);
M = e*sinhH-log(sinhH+sqrt(1+sinhH^2));
end

类别

Help CenterFile Exchange 中查找有关 Satellite and Orbital Mechanics 的更多信息

Community Treasure Hunt

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

Start Hunting!

Translated by