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

 

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