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.

kepler1 (manom, ecc)
function [eanom, tanom] = kepler1 (manom, ecc)
     
% solve Kepler's equation for circular,
% elliptic and hyperbolic orbits
   
% Danby's method

% input

%  manom = mean anomaly (radians)
%  ecc   = orbital eccentricity (non-dimensional)

% output

%  eanom = eccentric anomaly (radians)
%  tanom = true anomaly (radians)

% Celestial Computing with MATLAB

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

global niter

% define convergence criterion

ktol = 1.0e-10;

pi2 = 2 * pi;

xma = manom - pi2 * fix(manom/pi2);

% initial guess

if (ecc == 0)
   % circular orbit
   
   tanom = xma;
   eanom = xma;
   return;
elseif (ecc < 1)
   % elliptic orbit
   
   eanom = xma + 0.85 * sign(sin(xma)) * ecc;   
   
else
   % hyperbolic orbit
   
   eanom = log(2 * xma / ecc + 1.8);
end

% perform iterations

niter = 0;

while(1)
   if (ecc < 1)
   
      % elliptic orbit
      
      s = ecc * sin(eanom);
      c = ecc * cos(eanom);
      
      f = eanom - s - xma;  
      fp = 1 - c;
      fpp = s;
      fppp = c;
   else
   
      % hyperbolic orbit
      
      s = ecc * sinh(eanom);
      c = ecc * cosh(eanom);
      
      f = s - eanom - xma;
      fp = c - 1;
      fpp = s;
      fppp = c;
   end
   
   niter = niter + 1;
   
   % check for convergence

   if (abs(f) <= ktol | niter > 20)
      break;
   end
   
   % update eccentric anomaly
         
   delta = -f / fp;
   
   deltastar = -f / (fp + 0.5 * delta * fpp);
   
   deltak = -f / (fp + 0.5 * deltastar * fpp ...
            + deltastar * deltastar * fppp / 6);
         
   eanom = eanom + deltak;
end

if (niter > 20)
   clc; home;
   fprintf('\n\n   more than 20 iterations in kepler1 \n\n');
   keycheck;
end

% compute true anomaly

if (ecc < 1)
   % elliptic orbit
   
   sta = sqrt(1 - ecc * ecc) * sin(eanom);
   cta = cos(eanom) - ecc;
else
   % hyperbolic orbit
   
   sta = sqrt(ecc * ecc - 1) * sinh(eanom);
   cta = ecc - cosh(eanom);
end

tanom = atan3(sta, cta);


Contact us