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.

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