Code covered by the BSD License
 atan3 (a, b)
four quadrant inverse tangent
 cae2o(x)
Earthtoobject closest approach objective function
 cae2oprt(jdate)
print Earthtoobject closest approach conditions
 caeqm (time, y)
heliocentric equations of motion
 caevent(objfunc, prtfunc,...predict close approach events
 eci2orb1 (mu, r, v)
convert eci state vector to six classical orbital
 gdate (jdate)
convert Julian date to Gregorian (calendar) date
 hel_eqm (t, y)
heliocentric equations of motion
 heqms (t, y)
first order form of the heliocentric equations of motion
 jd2str(jdate)
convert Julian date to string equivalent
 jplephem (et, ntarg, ncent)
reads the jpl planetary ephemeris and gives
 julian (month, day, year)
Julian date
 kepler1 (manom, ecc)
solve Kepler's equation for circular,
 keycheck
pause and request user input
 minima (f, a, b, tolm)
onedimensional minimization
 oeprint(mu, oev)
print six classical orbital elements
 oeprint1(mu, oev, ittype)
print six classical orbital elements
 orb2eci(mu, oev)
convert classical orbital elements to eci state vector
 readoe2(filename)
read orbital elements data file
 rkf78 (deq, neq, ti, tf, ...solve first order system of differential equations
 svprint(r, v)
print position and velocity vectors and magnitudes
 uvector(x)
unit vector
 cae2ho.m

View all files
Closest Approach Between the Earth and Heliocentric Objects
by
David Eagle
06 Dec 2012
(Updated
02 Dec 2013)
MATLAB script that predicts closest approach between the Earth and heliocentric objects.

cae2o(x)

function y = cae2o(x)
% Earthtoobject closest approach objective function
% required by cae2ho.m
% input
% x = simulation time argument (days)
% output
% y = Earthtoobject distance (au)
% Orbital Mechanics with MATLAB
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
global tisaved ysaved jdatei drsaved
% compute object's position vector
ttmp = tisaved;
for i = 1:1:6
ytmp(i) = ysaved(i);
end
% number of differential equations
neq = 6;
% step size guess (days)
h = 0.1;
% integrate from tisaved to requested time = x
tetol = 1.0e15;
ysaved = rkf78('heqms', neq, ttmp, x, h, tetol, ytmp);
% compute object's position vector
for i = 1:1:3
robject(i) = ysaved(i);
end
% compute heliocentric position vector of the Earth
jdate = jdatei + x;
result = jplephem(jdate, 3, 11);
% vector from Earth to object
for i = 1:1:3
re2a(i) = robject(i)  result(i);
end
% Earthtoobject scalar distance (au)
srange = norm(re2a);
drsaved = srange;
% save objective function value and current time
y = srange;
tisaved = x;


Contact us