Code covered by the BSD License  

Highlights from
The Gravity Perturbed Hohmann Transfer

from The Gravity Perturbed Hohmann Transfer by David Eagle
MATLAB script for solving the Hohmann transfer problem perturbed by non-spherical Earth gravity.

[fid, alt1, alt2, inc1, inc2, raan1, jdutc0, gst0, ...
function [fid, alt1, alt2, inc1, inc2, raan1, jdutc0, gst0, ...
    req, mu, lgrav mgrav, gravfname] = readdata(filename)

% read orbital elements and simulation
% control data file

% required by phohmann.m

% input

%  filename = name of oota data file

% output

%  fid = file id

%  oev1 = initial orbit orbital elements array
%  oev2 = final orbit orbital elements array

% Orbital Mechanics with MATLAB

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

% 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 67 lines of data

for i = 1:1:67
    
    cline = fgetl(fid);
    
    switch i
        
        case 10
            
            % initial UTC calendar date
            
            tl = size(cline);
            
            ci = findstr(cline, ',');
            
            % extract month, day and year
            
            month = str2double(cline(1:ci(1)-1));
            
            day = str2double(cline(ci(1)+1:ci(2)-1));
            
            year = str2double(cline(ci(2)+1:tl(2)));
            
        case 15
            
            % initial UTC epoch
            
            utstr = cline;
            
            tl = size(utstr);
            
            ci = findstr(utstr, ',');
            
            % extract hours, minutes and seconds
            
            utc_hr = str2num(utstr(1:ci(1)-1));
            
            utc_min = str2num(utstr(ci(1)+1:ci(2)-1));
            
            utc_sec = str2num(utstr(ci(2)+1:tl(2)));
            
            % compute utc julian date
            
            day = day + utc_hr / 24.0 + utc_min / 1440.0 + utc_sec / 86400.0;
            
            jdutc0 = julian(month, day, year);
            
            gst0 = gast1 (jdutc0);
            
        case 23
            
            alt1 = str2double(cline);
            
        case 27
            
            inc1 = str2double(cline);
            
        case 31
            
            raan1 = str2double(cline);
            
        case 39
            
            alt2 = str2double(cline);
            
        case 43
            
            inc2 = str2double(cline);
            
        case 51
            
            mu = str2double(cline);
            
        case 55
            
            req = str2double(cline);
            
        case 59
            
            gravfname = cline;
            
        case 63
            
            lgrav = str2double(cline);
            
        case 67
            
            mgrav = str2double(cline);
                 
    end
end

status = fclose(fid);

Contact us