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.

rkf78 (deq, neq, ti, tf, h, tetol, x)
function xout = rkf78 (deq, neq, ti, tf, h, tetol, x)

% solve first order system of differential equations

% Runge-Kutta-Fehlberg 7(8) method

% input

%  deq   = name of function which defines the
%          system of differential equations
%  neq   = number of differential equations
%  ti    = initial simulation time
%  tf    = final simulation time
%  h     = initial guess for integration step size
%  tetol = truncation error tolerance (non-dimensional)
%  x     = integration vector at time = ti

% output

%  xout  = integration vector at time = tf

% Celestial Computing with MATLAB

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

global rkcoef ch alph beta

if (rkcoef == 1)

   % allocate arrays

   ch = zeros(13, 1);
   alph = zeros(13, 1);
   beta = zeros(13, 12);

   % define integration coefficients

   ch(6) = 34 / 105;
   ch(7) = 9 / 35;
   ch(8) = ch(7);
   ch(9) = 9 / 280;
   ch(10) = ch(9);
   ch(12) = 41 / 840;
   ch(13) = ch(12);
 
   alph(2) = 2 / 27;
   alph(3) = 1 / 9;
   alph(4) = 1 / 6;
   alph(5) = 5 / 12;
   alph(6) = 0.5;
   alph(7) = 5 / 6;
   alph(8) = 1 / 6;
   alph(9) = 2 / 3;
   alph(10) = 1 / 3;
   alph(11) = 1;
   alph(13) = 1;

   beta(2, 1) = 2 / 27;
   beta(3, 1) = 1 / 36;
   beta(4, 1) = 1 / 24;
   beta(5, 1) = 5 / 12;
   beta(6, 1) = 0.05;
   beta(7, 1) = -25 / 108;
   beta(8, 1) = 31 / 300;
   beta(9, 1) = 2;
   beta(10, 1) = -91 / 108;
   beta(11, 1) = 2383 / 4100;
   beta(12, 1) = 3 / 205;
   beta(13, 1) = -1777 / 4100;
   beta(3, 2) = 1 / 12;
   beta(4, 3) = 1 / 8;
   beta(5, 3) = -25 / 16;
   beta(5, 4) = -beta(5, 3);
   beta(6, 4) = 0.25;
   beta(7, 4) = 125 / 108;
   beta(9, 4) = -53 / 6;
   beta(10, 4) = 23 / 108;
   beta(11, 4) = -341 / 164;
   beta(13, 4) = beta(11, 4);
   beta(6, 5) = 0.2;
   beta(7, 5) = -65 / 27;
   beta(8, 5) = 61 / 225;
   beta(9, 5) = 704 / 45;
   beta(10, 5) = -976 / 135;
   beta(11, 5) = 4496 / 1025;
   beta(13, 5) = beta(11, 5);
   beta(7, 6) = 125 / 54;
   beta(8, 6) = -2 / 9;
   beta(9, 6) = -107 / 9;
   beta(10, 6) = 311 / 54;
   beta(11, 6) = -301 / 82;
   beta(12, 6) = -6 / 41;
   beta(13, 6) = -289 / 82;
   beta(8, 7) = 13 / 900;
   beta(9, 7) = 67 / 90;
   beta(10, 7) = -19 / 60;
   beta(11, 7) = 2133 / 4100;
   beta(12, 7) = -3 / 205;
   beta(13, 7) = 2193 / 4100;
   beta(9, 8) = 3;
   beta(10, 8) = 17 / 6;
   beta(11, 8) = 45 / 82;
   beta(12, 8) = -3 / 41;
   beta(13, 8) = 51 / 82;
   beta(10, 9) = -1 / 12;
   beta(11, 9) = 45 / 164;
   beta(12, 9) = 3 / 41;
   beta(13, 9) = 33 / 164;
   beta(11, 10) = 18 / 41;
   beta(12, 10) = 6 / 41;
   beta(13, 10) = 12 / 41;
   beta(13, 12) = 1;

   % set initialization indicator

   rkcoef = 0;
end

% compute integration "direction"

sdt = sign(tf - ti);

dt = abs(h) * sdt;

while (1)

   % load "working" time and integration vector

   twrk = ti;
       
   xwrk = x;

   % check for last dt

   if (abs(dt) > abs(tf - ti)) 
      dt = tf - ti;
   end

   % check for end of integration period

   if (abs(ti - tf) < 0.00000001)
      xout = x;
      return;
   end

   % evaluate equations of motion

   xdot = feval(deq, ti, x);
   
   for i = 1:1:neq
       f(i, 1) = xdot(i);
   end

   % compute solution

   for k = 2:1:13
       kk = k - 1;
       for i = 1:1:neq
           temp = 0;
           for j = 1:1:kk
               temp = temp + beta(k, j) * f(i, j);
           end
           x(i) = xwrk(i) + dt * temp;
       end
       
       ti = twrk + alph(k) * dt;

       xdot = feval(deq, ti, x);

       for j = 1:1:neq
           f(j, k) = xdot(j);
       end

   end
       
   for i = 1:1:neq
       temp = 0;

       for l = 1:1:13
           temp = temp + ch(l) * f(i, l);
       end

       x(i) = xwrk(i) + dt * temp;
   end

   % truncation error calculations

   xerr = tetol;

   for i = 1:1:neq
       ter = abs((f(i, 1) + f(i, 11) - f(i, 12) - f(i, 13)) * ch(12) * dt);
       tol = abs(x(i)) * tetol + tetol;
       tconst = ter / tol;
       if (tconst > xerr)
          xerr = tconst;
       end
   end

   % compute new step size

   dt = 0.8 * dt * (1 / xerr) ^ (1 / 8);
   
   if (xerr > 1)
      % reject current step
      ti = twrk;
      x = xwrk;
   else
      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      % accept current step              %
      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      % perform graphics, additional     %
      % calculations, etc. at this point %
      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   end
end



Contact us