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

David Eagle (view profile)

 

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