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

how to get the orbital parameters of a satellite formation using at low earth orbit?

[James Tursa]

Please be more specific. Do you have position & velocity vectors? Are you trying to do orbit determination? Or what?

Answer 1

[Meysam Mahooti]

%--------------------------------------------------------------------------

%

% 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