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

4 views (last 30 days)
Aswathi M on 7 Oct 2015
Answered: Meysam Mahooti on 26 May 2021
how to get the orbital parameters of a satellite formation using at low earth orbit?
James Tursa on 7 Oct 2015
Please be more specific. Do you have position & velocity vectors? Are you trying to do orbit determination? Or what?

Meysam Mahooti on 26 May 2021
%--------------------------------------------------------------------------
%
% 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
% Omega Longitude of the ascending node [rad]
% omega Argument of pericenter [rad]
% at time t_a
%
% Notes:
% The function cannot be used with state vectors describing a circular
% or non-inclined orbit.
%
%
%--------------------------------------------------------------------------
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