Code covered by the BSD License  

Highlights from
Closest Approach Between the Earth and Heliocentric Objects

Closest Approach Between the Earth and Heliocentric Objects

by

 

06 Dec 2012 (Updated )

MATLAB script that predicts closest approach between the Earth and heliocentric objects.

hel_eqm (t, y)
function ydot = hel_eqm (t, y)

% heliocentric equations of motion
% with planetary perturbations

% Battin's f(q) formulation

% input

%  t = current simulation time (days)

% output

%  ydot = first order equations of motion

% Orbital Mechanics with Matlab

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

global xmu jdtcm

% current julian date (relative to tcm event)

jdate = jdtcm + t;

% distance from sun to spacecraft

rsun2sc = norm(y(1:3));

rrsun2sc = -xmu(1) / rsun2sc^3;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% calculate planetary position vectors
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

nplanets = 6;

for i = 1:1:nplanets

    svplanet = jplephem(jdate, i, 11);

    rplanet = svplanet(1:3);

    for j = 1:1:3
        rp(i, j) = rplanet(j);
    end
end

% compute planet-centered position vector of spacecraft

for i = 1:1:nplanets
    rp2sc(i, 1) = y(1) - rp(i, 1);
    
    rp2sc(i, 2) = y(2) - rp(i, 2);
    
    rp2sc(i, 3) = y(3) - rp(i, 3);
end

% compute f(q) functions for each planet

for k = 1:1:nplanets
    q(k) = dot(y(1:3), y(1:3) - 2.0 * rp(k, :)') / dot(rp(k, :), rp(k, :));

    f(k) = q(k) * ((3.0 + 3.0 * q(k) + q(k) * q(k)) / (1.0 + (1.0 + q(k))^1.5));

    d3(k) = norm(rp2sc(k, :)) * norm(rp2sc(k, :)) * norm(rp2sc(k, :));
end

% compute planetary perturbations

for j = 1:1:3
    
    accp(j) = 0.0;

    for k = 1:1:nplanets
        
        accp(j) = accp(j) - xmu(k + 1) * (y(j) + f(k) * rp(k, j)) / d3(k);
        
    end
end

% compute integration vector

ydot = [ y(4)
    y(5)
    y(6)
    accp(1) + y(1) * rrsun2sc
    accp(2) + y(2) * rrsun2sc
    accp(3) + y(3) * rrsun2sc];

Contact us