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.

readoe2(filename)
function [fid, bodname, oev, iframe, jdate] = readoe2(filename)

% read orbital elements data file

% required by cae2ho.m

% input

%  filename = name of orbital elements data file

% output

%  fid = file id

%  oev(1) = semimajor axis (au)
%  oev(2) = orbital eccentricity
%  oev(3) = orbital inclination
%  oev(4) = argument of perigee
%  oev(5) = right ascension of the ascending node
%  oev(6) = mean anomaly
%  jdate  = TDB julian date
%  iframe = fundamental frame (1 = ecliptic, 2 = equator)

% NOTE: all angular elements are returned in radians

% Orbital Mechanics with MATLAB

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

dtr = pi / 180;

% open data file

fid = fopen(filename, 'r');

% check for file open error

if (fid == -1)
    
    clc; home;
    
    fprintf('\n\n error: cannot find this file!!');
    
    keycheck;
    
    return;
    
end

% read 54 lines of data file

for i = 1:1:54
    
    cline = fgetl(fid);
    
    switch i
        
        case 9
            
            bodname = cline;
            
        case 14
            
            tl = size(cline);
            
            ci = findstr(cline, ',');
            
            % extract month, day and year
            
            month = str2num(cline(1:ci(1)-1));
            
            day = str2num(cline(ci(1)+1:ci(2)-1));
            
            year = str2num(cline(ci(2)+1:tl(2)));
            
        case 19
            
            tl = size(cline);
            
            ci = findstr(cline, ',');
            
            % extract hours, minutes and seconds
            
            uthr = str2num(cline(1:ci(1)-1));
            
            utmin = str2num(cline(ci(1)+1:ci(2)-1));
            
            utsec = str2num(cline(ci(2)+1:tl(2)));
            
            % compute julian date
            
            jdate = julian(month, day, year) ...
                + uthr / 24 + utmin / 1440 + utsec / 86400;
            
        case 24
            
            oev(1) = str2num(cline);
            
        case 29
            
            oev(2) = str2num(cline);
            
        case 34
            
            oev(3) = dtr * str2num(cline);
            
        case 39
            
            oev(4) = dtr * str2num(cline);
            
        case 44
            
            oev(5) = dtr * str2num(cline);
            
        case 49
            
            oev(6) = dtr * str2num(cline);
            
        case 54
            
            iframe = str2num(cline);
    end
    
end

status = fclose(fid);

Contact us