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.

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