Code covered by the BSD License  

Highlights from
The Gravity Perturbed Hohmann Transfer

The Gravity Perturbed Hohmann Transfer

by

 

28 Feb 2013 (Updated )

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