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.

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