Code covered by the BSD License

# Closest Approach Between the Earth and Heliocentric Objects

### 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;

```