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.

cae2ho.m
% cae2ho.m          December 2, 2013

% closest approach conditions between the Earth and an asteroid,
% comet or object in a heliocentric elliptic orbit

% numerical integration using Cowell's
% method and one-dimensional minimization

% DE421 ephemeris

% Orbital Mechanics with MATLAB

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

global emu xmu xmmu jdatei rkcoef aunit csqr

global ysaved tisaved camin

global ephname iephem km

% J2000 equatorial-to-ecliptic transformation matrix

eq2000 = [[1.000000000000000 0  0]; ...
    [0   0.917482062069182   0.397777155931914]; ...
    [0  -0.397777155931914   0.917482062069182]];

% astrodynamic and utility constants

rtd = 180 / pi;               % radians to degrees

dtr = pi / 180;               % degrees to radians

% astronomical unit (kilometers)

aunit = 149597870.691d0;

% speed of light (meters/second; DE421 value)

clight = 299792.458d0;

% speed of light squared (au**2/day**2)

csqr = (clight * 86400.0d0 / aunit)^2;

% earth gravitational constant (km**3/sec**2)

emu = 398600.436233d0;

% moon gravitational constant (km**3/sec**2)

xmmu = 4902.800121244414d0;

xmmu = xmmu * 86400.0d0^2 / aunit^3;

% earth equatorial radius (kilometers)

req = 6378.1363d0;

% **************************************
% gravitational constants (au**3/day**2)
% DE 421 values
% **************************************

% sun

xmu(1) = 0.295912208285591100d-3;

smu = xmu(1);

% mercury

xmu(2) = 0.491254957186794000D-10;

% venus

xmu(3) = 0.724345233269844100D-09;

% earth

xmu(4) = emu * 86400.0d0^2 / aunit^3;

% mars

xmu(5) = 0.954954869562239000D-10;

% jupiter

xmu(6) = 0.282534584085505000D-06;

% saturn

xmu(7) = 0.845970607330847800D-07;

% uranus

xmu(8) = 0.129202482579265000D-07;

% neptune

xmu(9) = 0.152435910924974000D-07;

% pluto

xmu(10) = 0.217844105199052000D-11;

% initialize rkf integrator

rkcoef = 1;

% initialize DE421 algorithm

ephname = 'de421.bin';

iephem = 1;

km = 0;

% initialize deltat-t algorithm

itime = 1;

% begin simulation

clc; home;

fprintf('\n                       program cae2ho\n');

fprintf('\n< closest approach between the Earth and a heliocentric object >\n\n');

% request name of data file

[filename, pathname] = uigetfile('*.in', 'please select an input data file');

[fid, bodname, oev1, iframe, jdatei] = readoe2(filename);

while (1)
    
    fprintf('\nplease input the search duration (days)\n');
    
    ndays = input('? ');
    
    if (ndays > 0)
        
        break;
        
    end
    
end

while(1)
    
    fprintf('\n\nwould you like to enforce a minimum distance constraint (y = yes, n = no)\n');
    
    slct = input('? ', 's');
    
    if (slct == 'y' || slct == 'n')
        
        break;
        
    end
    
end

if (lower(slct) == 'y')
    
    % request minimum distance constraint
    
    while(1)
        
        fprintf('\nplease input the minimum distance constraint (AU)\n');
        
        camin = input('? ');
        
        if (camin >= 0)
            
            break;
            
        end
        
    end
    
else
    
    % no minimum distance constraint
    
    camin = 0;
    
end

% compute initial true anomaly (radians)

ecc = oev1(2);

manom = oev1(6);

[eanom, tanom] = kepler1 (manom, ecc);

oev1(6) = tanom;

% calculate initial heliocentric position and velocity vectors

[rwrk, vwrk] = orb2eci(xmu(1), oev1);

% calculate orbital period

period = sqrt(oev1(1) * oev1(1) * oev1(1) / xmu(1));

% load initial state vector array

if (iframe == 1)
    
    % convert ecliptic to equatorial
    
    rtmp = eq2000' * rwrk;
    
    vtmp = eq2000' * vwrk;
    
    for i = 1:1:3
        
        ysaved(i) = rtmp(i);
        
        ysaved(i + 3) = vtmp(i);
        
    end
    
    oev_wrk = eci2orb1 (xmu(1), rtmp, vtmp);
    
else
    
    % no transformation required
    
    for i = 1:1:3
        
        ysaved(i) = rwrk(i);
        
        ysaved(i + 3) = vwrk(i);
        
    end
   
    oev_wrk = eci2orb1 (xmu(1), rwrk, vwrk);
    
end

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

[cdstr, utstr] = jd2str(jdatei);

fprintf('\ncalendar date        ');

disp(cdstr);

fprintf('\nTDB time             ');

disp(utstr);

fprintf('\nTDB Julian date      %16.8f \n\n', jdatei);

oeprint(xmu(1), oev_wrk);

[r, v] = orb2eci(xmu(1), oev_wrk);

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

% initialize search parameters

ti = 0;

tf = ndays;

tisaved = ti;

dt = period / 4;

dtsml = 0.1;

% find close approach events

caevent('cae2o', 'cae2oprt', ti, tf, dt, dtsml);

Contact us