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.

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

% first order form of the heliocentric equations of motion
% in the EME2000 coordinate system

% JPL ephemeris and Battin's f(q) formulation

% input

%  time = elapsed time since jdate0 (days)
%  y    = current state vector (au and au/day)

% output

%  ydot = equations of motion

% Orbital Mechanics with Matlab

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

global xmu xmmu jdatei csqr

accp = zeros(3, 1);

accm = zeros(3, 1);

accr = zeros(3, 1);

rp = zeros(9, 3);

rp2sc = zeros(9, 3);

q = zeros(9, 1);

f = zeros(9, 1);

d3 = zeros(9, 1);

% current TDB julian date

jdate = jdatei + t;

% distance from sun to spacecraft (kilometers)

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

% point-mass solar gravity multiplier

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

% calculate planetary position vectors

for i = 1:1:9
    
    rrd = jplephem (jdate, i, 11);
    
    for j = 1:1:3
        
        rp(i, j) = rrd(j);
        
    end
    
end

% compute planet-centered position vectors of spacecraft

for i = 1:1:9
    
    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:9
    
    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
    
    for k = 1:1:9
        
        accp(j) = accp(j) - xmu(k + 1) * (y(j) + f(k) * rp(k, j)) / d3(k);
        
    end
    
end

% compute lunar perturbation

rrd = jplephem (jdate, 10, 11);

rmoon = rrd(1:3);

% compute selenocentric position vector of the spacecraft

for i = 1:1:3
    
    rm2sc(i) = y(i) - rmoon(i);
    
end

% f(q) formulation

for i = 1:1:3
    
    vtmp(i) = y(i) - 2.0d0 * rmoon(i);
    
end

dot1 = dot(y(1:3), vtmp);

dot2 = dot(rmoon, rmoon);

qmoon = dot1 / dot2;

fmoon = qmoon * ((3.0d0 + 3.0d0 * qmoon + qmoon * qmoon) ...
    / (1.0d0 + (1.0d0 + qmoon)^1.5d0));

d3moon = norm(rm2sc)^3;

% point-mass gravity of the moon

for i = 1:1:3
    
    accm(i) = -xmmu * (y(i) + fmoon * rmoon(i)) / d3moon;
    
end

% compute relativistic perturbation

rs2sc = sqrt(y(1) * y(1) + y(2) * y(2) + y(3) * y(3));

rrs2sc = -xmu(1) / rs2sc^3;

vasqrm = y(4) * y(4) + y(5) * y(5) + y(6) * y(6);

rdotv = y(1) * y(4) + y(2) * y(5) + y(3) * y(6);

for i = 1:1:3
    
    accr(i) = (-rrs2sc / csqr) ...
        * (((4.0d0 * xmu(1) / rs2sc) - vasqrm) ...
        * y(i) + 4.0d0 * rdotv * y(i + 3));
    
end

% compute integration vector

ydot(1) = y(4);
ydot(2) = y(5);
ydot(3) = y(6);

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

Contact us