Code covered by the BSD License  

Highlights from
Optimal Finite-burn Interplanetary Injection from Earth Orbit

from Optimal Finite-burn Interplanetary Injection from Earth Orbit by David Eagle
A MATLAB script for optimizing finite-burn interplanetary injection trajectories.

[fid, itarget, nsegments, xmass0, thrmag, xisp, oev_initial, ...
function [fid, itarget, nsegments, xmass0, thrmag, xisp, oev_initial, ...
    c3_target, rla_target, dla_target, tcoast, igrav] = read_escape(filename)

% read finite-burn hyperbolic injection data file

% input

%  filename = name of data file

% output

%  fid = file id

%  itarget        = type of targeting/optimization
%                   (1 = C3 only, 2 = C3 and DLA, 3 = C3, RLA and DLA)
%  nsegments      = number of time segments
%  xmass0         = initial spacecraft mass (kilograms)
%  thrmag         = thrust magnitude (newtons)
%  xisp           = specific impulse (seconds)
%  oev_initial(1) = initial semimajor axis (kilometers)
%  oev_initial(2) = initial orbital eccentricity
%  oev_initial(3) = initial orbital inclination (radians)
%  oev_initial(4) = initial argument of perigee (radians)
%  oev_initial(5) = initial right ascension of the ascending node (radians)
%  oev_initial(6) = initial true anomaly (radians)
%  c3_target      = final c3 orbital energy target (km/sec)^2
%  dla_target     = final dla target (radians)
%  tcoast         = final coast time (seconds)
%  igrav          = j2 gravity perturbation flag

% Orbital Mechanics with MATLAB

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

dtr = pi / 180.0;

% 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!!');
    
    return;
    
end

% read 65 lines of data file

for i = 1:1:65
    
    cline = fgetl(fid);
    
    switch i
        
        case 12
            
            % type of targeting/optimization
            
            itarget = str2double(cline);
            
        case 15
            
            % number of trajectory segments
            
            nsegments = str2double(cline);
            
        case 18
            
            % initial spacecraft mass (kilograms)
            
            xmass0 = str2double(cline);
            
        case 21
            
            % thrust magnitude (newtons)
            
            thrmag = str2double(cline);
            
        case 24
            
            % specific impulse (seconds)
            
            xisp = str2double(cline);
            
        case 31
            
            % initial semimajor axis (kilometers)
            
            oev_initial(1) = str2double(cline);
            
        case 34
            
            % initial orbital eccentricity
            
            oev_initial(2) = str2double(cline);
            
        case 37
            
            % initial orbital inclination (radians)
            
            oev_initial(3) = dtr * str2double(cline);
            
        case 40
            
            % argument of perigee (radians)
            
            oev_initial(4) = dtr * str2double(cline);
            
        case 43
            
            % initial raan (radians)
            
            oev_initial(5) = dtr * str2double(cline);
            
        case 46
            
            % initial true anomaly (radians)
            
            oev_initial(6) = dtr * str2double(cline);
            
        case 53
            
            % c3 orbital energy target (km/sec)^2
            
            c3_target = str2double(cline);
            
        case 56
            
            % RLA target (radians)
            
            rla_target = dtr * str2double(cline);
            
        case 59
            
            % DLA target (radians)
            
            dla_target = dtr * str2double(cline);
            
        case 62
            
            % final coast time (seconds)
            
            tcoast = 60.0 * str2double(cline);
            
        case 65
            
            % j2 gravity perturbation indicator
            
            igrav = str2double(cline);
            
    end
    
end

status = fclose(fid);

Contact us