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.

caeqm (time, y)
function yddot = caeqm (time, y)
   
% heliocentric equations of motion
% with planetary perturbations

% required by cae2ho.m

% Celestial Computing with MATLAB

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

global aunit xmu jdatei iephtype

% current julian date

jdate = jdatei + time;

% distance from sun to object

rs2sc = sqrt(y(1) * y(1) + y(2) * y(2) + y(3) * y(3));
   
rrs2sc = -xmu(1) / rs2sc^3;
   
% calculate planetary position vectors

% Mercury

if (iephtype == 1)
   result = jplephem(jdate, 2, 11);   
else
   [result, ierr] = slp96 (jdate, 1, 11, 1, 1, 1);
end

for i = 1:1:3
    rmercury(i) = result(i);
end

rpm(1) = sqrt(dot(rmercury, rmercury));

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

% Venus

if (iephtype == 1)
   result = jplephem(jdate, 2, 11);
else
   [result, ierr] = slp96 (jdate, 2, 11, 1, 1, 1);
end

for i = 1:1:3
   rvenus(i) = result(i);
end

rpm(2) = sqrt(dot(rvenus, rvenus));

for i = 1:1:3
    rp(2, i) = rvenus(i);
end

% Earth

if (iephtype == 1)
   result = jplephem(jdate, 3, 11);
else
   [result, ierr] = slp96 (jdate, 12, 11, 1, 1, 1);
end

for i = 1:1:3
    rearth(i) = result(i);
end

rpm(3) = sqrt(dot(rearth, rearth));

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

% Mars

if (iephtype == 1)
   result = jplephem(jdate, 4, 11);
else
   [result, ierr] = slp96 (jdate, 4, 11, 1, 1, 1);
end

for i = 1:1:3
   rmars(i) = result(i);
end

rpm(4) = sqrt(dot(rmars, rmars));

for i = 1:1:3
    rp(4, i) = rmars(i);
end

% Jupiter

if (iephtype == 1)
   result = jplephem(jdate, 5, 11);
else
   [result, ierr] = slp96 (jdate, 5, 11, 1, 1, 1);   
end

for i = 1:1:3
   rjupiter(i) = result(i);
end

rpm(5) = sqrt(dot(rjupiter, rjupiter));

for i = 1:1:3
    rp(5, i) = rjupiter(i);
end

% Saturn

if (iephtype == 1)
   result = jplephem(jdate, 6, 11);
else
   [result, ierr] = slp96 (jdate, 6, 11, 1, 1, 1);   
end

for i = 1:1:3
   rsaturn(i) = result(i);
end

rpm(6) = sqrt(dot(rsaturn, rsaturn));

for i = 1:1:3
    rp(6, i) = rsaturn(i);
end

for i = 1:1:6
    rrpm(i) = xmu(i + 1) / (rpm(i)^3);
end

% compute planetocentric position vectors of spacecraft

for i = 1:1:6
    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 planetocentric distances of spacecraft

for i = 1:1:6
    rp2scm(i) = sqrt(rp2sc(i, 1)^2 + rp2sc(i, 2)^2 + rp2sc(i, 3)^2);
    rrp2scm(i) = xmu(i + 1) / (rp2scm(i)^3);
end

% compute planetary perturbations

for j = 1:1:3
    accp(j) = 0;
    for i = 1:1:6
        accp(j) = accp(j) - rp2sc(i, j) * rrp2scm(i) - rp(i, j) * rrpm(i);
    end
end

% compute integration vector

yddot(1) = y(4);

yddot(2) = y(5);

yddot(3) = y(6);

yddot(4) = accp(1) + y(1) * rrs2sc;
    
yddot(5) = accp(2) + y(2) * rrs2sc;
    
yddot(6) = accp(3) + y(3) * rrs2sc;

Contact us