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.

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

Contact us