Code covered by the BSD License  

Highlights from
Optimal Impulsive Orbital Transfer

Optimal Impulsive Orbital Transfer

by

 

27 Nov 2012 (Updated )

MATLAB script for the solution of the one and two impulse orbit transfer between two Earth orbits.

oota_matlab.m
% oota_matlab.m        July 15, 2013

% optimal impulsive orbital transfer analysis

% true anomaly grid search, three-dimensional trajectory
% graphics and primer vector optimality analysis

% Orbital Mechanics with MATLAB

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

clear all;

% global information

global dtr rtd mu req output

global top2eci soln filename fname

global sma ecc inc argper raan

global r1 r2 vt1 vt2 dv1 dv2

% initialize solution number

soln = 0;

% angular conversion factors

dtr = pi / 180.0;

rtd = 180.0 / pi;

% begin simulation

clc; home;

fprintf('\n         program oota_matlab\n');

fprintf('\n< optimal orbital transfer analysis >\n\n');

% request name of user-defined input data file

[filename, pathname] = uigetfile('*.in', 'please select an input data file');

% read simulation definition data file

[fid, mu, req, oev1, oev2, phi01, phi02, ...
    dphi1, dphi2, nphi1, nphi2] = read_oota(filename);

% extract file name up to first period

l1 = strfind(filename, '.');

fname = filename(1:l1 -1);

% load classical orbital elements of the initial orbit

sma(1) = oev1(1);
ecc(1) = oev1(2);
inc(1) = oev1(3) * dtr;
argper(1) = oev1(4) * dtr;
raan(1) = oev1(5) * dtr;

% load classical orbital elements of the final orbit

sma(2) = oev2(1);
ecc(2) = oev2(2);
inc(2) = oev2(3) * dtr;
argper(2) = oev2(4) * dtr;
raan(2) = oev2(5) * dtr;

% perform algorithm initialization

oota_iniz;

% convert initial search angles to radians

phi01 = phi01 * dtr;

phi02 = phi02 * dtr;

% convert initial search angle increments to radians

dphi1 = dphi1 * dtr;

dphi2 = dphi2 * dtr;

tol = 1.0e-4;

% perform two-dimensional grid search

x0 = zeros(2, 1);

for i = 1:1:nphi1
    
    for j = 1:1:nphi2
        
        x0(1) = phi01 + dphi1 * (i - 1);
        
        x0(2) = phi02 + dphi2 * (j - 1);
        
        % save current values
        
        xs1 = x0(1);
        
        xs2 = x0(2);
        
        dx = abs(x0(2) - x0(1));
        
        if ((abs(dx - pi) > tol) && (dx > tol) && (abs(dx - 2.0 * pi) > tol))
            
            options = optimset('MaxFunEvals', 500, 'TolFun', 1.0e-8, 'TolX', 1.0e-8);
            
            % perform 2-dimensional minimization
            
            [xmin, fval, exitflag, output] = fminsearch('ootafun1', x0, options);
            
            % evaluate the current solution
            
            f = ootafun1(xmin);
            
            % compute eci delta-v vectors
            
            dv1eci = top2eci * dv1';
            
            dv2eci = top2eci * dv2';
            
            % compute eci delta-v magnitudes
            
            dv1ecim = norm(dv1eci);
            
            dv2ecim = norm(dv2eci);
            
            % transfer orbit - first impulse
            
            rt1eci = top2eci * r1';
            
            vt1eci = top2eci * vt1';
            
            oevt1 = eci2orb1 (mu, rt1eci, vt1eci);
            
            % transfer orbit - second impulse
            
            rt2eci = top2eci * r2';
            
            vt2eci = top2eci * vt2';
            
            oevt2 = eci2orb1 (mu, rt2eci, vt2eci);
            
            % transfer time-of-flight (seconds)
            
            dtof = tof1(mu, oevt1(1), oevt1(2), oevt1(6), oevt2(6));
            
            if (dv1ecim > 0.001 && dv2ecim > 0.001)
                
                % --------------------
                % two impulse solution
                % --------------------
                
                % analyze primer behavior for optimality
                
                iopt = primer(dtof, rt1eci, vt1eci, dv1eci, dv2eci);
                
            else
                
                % --------------------
                % one impulse solution
                % --------------------
                
                soln = soln + 1;
                
                make_textfile(soln);
                
                tgraphics(soln, dtof);
                
            end
            
        end
        
    end
    
end

close('all');

fclose('all');

Contact us