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.

hohmann.m
% hohmann.m       September 11, 2012

% Hohmann two impulse orbit transfer between
% planar and non-coplanar circular orbits

% Orbital Mechanics with MATLAB

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

global rtd dtr

global mu req hn1 hn2 hn3 dinc

% astrodynamic and utility constants

om_constants;

% Brent root-finding tolerance

rtol = 0.00000001;

% request inputs

clc; home;

fprintf('\n   Hohmann Orbit Transfer Analysis\n');

while (1)
    fprintf('\n\nplease input the initial altitude (kilometers)\n');

    alt1 = input('? ');

    if (alt1 > 0.0)
        break;
    end
end

while (1)
    fprintf('\n\nplease input the final altitude (kilometers)\n');

    alt2 = input('? ');

    if (alt2 > 0.0)
        break;
    end
end

while (1)
    fprintf('\n\nplease input the initial orbital inclination (degrees)');
    fprintf('\n(0 <= inclination <= 180)\n');

    inc1 = input('? ');

    if (inc1 >= 0.0 && inc1 <= 180.0)
        break;
    end
end

while (1)
    fprintf('\n\nplease input the final orbital inclination (degrees)');
    fprintf('\n(0 <= inclination <= 180)\n');

    inc2 = input('? ');

    if (inc2 >= 0.0 && inc2 <= 180.0)
        break;
    end
end

% convert inclinations to radians

inc1 = inc1 * dtr;

inc2 = inc2 * dtr;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% solve the orbit transfer problem %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% calculate total inclination change (radians)

dinc = abs(inc2 - inc1);

% compute geocentric radii of initial and final orbits (kilometers)

r1 = req + alt1;

r2 = req + alt2;

% compute "normalized" radii

hn1 = sqrt(2.0 * r2 / (r2 + r1));

hn2 = sqrt(r1 / r2);

hn3 = sqrt(2.0 * r1 / (r2 + r1));

% compute "local circular velocity" of initial and final orbits (km/sec)

v1 = sqrt(mu / r1);

v2 = sqrt(mu / r2);

% compute transfer orbit semimajor axis (kilometers)

smat = 0.5 * (r1 + r2);

% compute transfer orbit eccentricity (non-dimensional)

ecct = (max(r1, r2) - min(r1, r2)) / (r1 + r2);

% compute transfer orbit perigee and apogee radii and velocities

rp = smat * (1.0 - ecct);

ra = smat * (1.0 + ecct);

vt1 = sqrt(2.0 * mu * ra / (rp * (rp + ra)));

vt2 = sqrt(2.0 * mu * rp / (ra * (rp + ra)));

% compute transfer orbit period (seconds)

taut = 2.0 * pi * sqrt(smat^3 / mu);

if (abs(dinc) == 0)
    
    % coplanar orbit transfer

    if (r2 > r1)
       
       % higher-to-lower transfer
        
       dv1 = vt1 - v1;
    
       dv2 = v2 - vt2;
       
    else
        
       % lower-to-higher transfer
       
       dv1 = v1 - vt2;
       
       dv2 = vt1 - v2;
       
    end
    
    dinc1 = 0;
    
    dinc2 = 0;
    
else
    
    % non-coplanar orbit transfer

    [xroot, froot] = brent('hohmfunc', 0, dinc, rtol);

    % calculate delta-v's

    dinc1 = xroot;

    dinc2 = dinc - dinc1;

    dv1 = v1 * sqrt(1 + hn1 * hn1 - 2.0 * hn1 * cos(dinc1));

    dv2 = v1 * sqrt(hn2 * hn2 * hn3 * hn3 + hn2 * hn2 ...
        - 2.0 * hn2 * hn2 * hn3 * cos(dinc2));
end

% print results

clc; home;

fprintf('\n   Hohmann Orbit Transfer Analysis \n\n\n');

fprintf('initial orbit altitude            %10.4f kilometers \n\n', alt1);

fprintf('initial orbit inclination         %10.4f degrees \n\n', inc1 * rtd);

fprintf('initial orbit velocity            %10.4f meters/second \n\n\n', 1000 * v1);

fprintf('final orbit altitude              %10.4f kilometers \n\n', alt2);

fprintf('final orbit inclination           %10.4f degrees \n\n', inc2 * rtd);

fprintf('final orbit velocity              %10.4f meters/second \n', 1000 * v2);

fprintf('\n\nfirst inclination change          %10.4f degrees\n\n', dinc1 * rtd);

fprintf('second inclination change         %10.4f degrees\n\n', dinc2 * rtd);

fprintf('total inclination change          %10.4f degrees\n\n\n', rtd * (dinc1 + dinc2));

fprintf('first delta-v                     %10.4f meters/second \n\n', 1000 * dv1);

fprintf('second delta-v                    %10.4f meters/second \n\n', 1000 * dv2);

fprintf('total delta-v                     %10.4f meters/second \n\n\n', 1000 * (dv1 + dv2));

fprintf('transfer orbit semimajor axis     %10.4f kilometers \n\n', smat);

fprintf('transfer orbit eccentricity       %10.8f \n\n', ecct);

fprintf('transfer orbit perigee velocity   %10.4f meters/second \n\n', 1000 * vt1);

fprintf('transfer orbit apogee velocity    %10.4f meters/second \n\n', 1000 * vt2);

fprintf('transfer orbit coast time         %10.4f seconds \n', 0.5 * taut);

fprintf('                                  %10.4f minutes \n', 0.5 * taut / 60.0);

fprintf('                                  %10.4f hours \n\n', 0.5 * taut / 3600.0);

Contact us