Code covered by the BSD License

# The Gravity Perturbed Hohmann Transfer

### David Eagle (view profile)

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

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;

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);
```