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.

caevent(objfunc, prtfunc, ti, tf, dt, dtsml)
```function caevent(objfunc, prtfunc, ti, tf, dt, dtsml)

% predict close approach events

% input

%  objfunc = objective function
%  prtfunc = display results function
%  ti      = initial simulation time
%  tf      = final simulation time
%  dt      = step size used for bounding minima
%  dtsml   = small step size used to determine whether
%            the function is increasing or decreasing

% Orbital Mechanics with MATLAB

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

% convergence tolerance

tolm = 1.0e-8;

% initialization

df = 1;

t = ti;

while (1)

% find where function first starts decreasing

while (1)
if (df <= 0)
break;
end

t = t + dt;

ft = feval(objfunc, t);

ftemp = feval(objfunc, t - dtsml);

df = ft - ftemp;

end

% function decreasing - find where function
% first starts increasing

while (1)

el = t;

t = t + dt;

ft = feval(objfunc, t);

ftemp = feval(objfunc, t - dtsml);

df = ft - ftemp;

if (df > 0)
break;
end

end

er = t;

% calculate minimum using Brent's method

[tmin, fmin] = minima(objfunc, el, er, tolm);

el = er;

% print results at extrema

prtevent(objfunc, prtfunc, tmin);

% check for end of simulation

if (t >= tf)

break;

end

end

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

function prtevent (objfunc, prtfunc, topt)

% display close approach conditions

% input

%  objfunc = objective function
%  prtfunc = display results function
%  topt    = extrema time

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

global xmu aunit jdatei ysaved drsaved camin

if (drsaved > camin)

return;

end

% evaluate objective function

froot = feval(objfunc, topt);

jdate = jdatei + topt;

feval(prtfunc, jdate);

fprintf ('\nheliocentric orbital elements');
fprintf ('\n(Earth mean equator and equinox J2000)');
fprintf ('\n--------------------------------------\n');

for i = 1:1:3

r(i) = ysaved(i);

v(i) = ysaved(i + 3);

end

oev = eci2orb1 (xmu(1), r, v);

oeprint(xmu(1), oev);

svprint(aunit * r, aunit * v / 86400.0);

```