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.

phohmann.m
```% phohmann.m       July 21, 2013

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

% Orbital Mechanics with MATLAB

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

clear all;

global rtd dtr jdutc0 gst0 rkcoef

global mu req omega hn1 hn2 hn3 dinc

global lgrav mgrav ccoef scoef j2 j4 tof

global ri vi xi itarget

global tevent yevent mee_target mee_final

dtr = pi / 180.0;

rtd = 180.0 / pi;

% initialize rkf78 method

rkcoef = 1;

% earth inertial rotation rate (radians/second)

omega = 7.292115486e-5;

% Brent root-finding tolerance

rtol = 1.0d-8;

% request inputs

clc; home;

fprintf('\nGravity-perturbed Hohmann Orbit Transfer Analysis\n');

% request simulation data file

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

[fid, alt1, alt2, inc1, inc2, raan1, jdutc0, gst0, req, mu, lgrav mgrav, ...

% read earth gravity data file

% extract j2 and j4 zonal gravity coefficients

j2 = -ccoef(3, 1);

j4 = -ccoef(5, 1);

inc1 = inc1 * dtr;

inc2 = inc2 * dtr;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% solve the two-body orbit Hohmann 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('\ntwo-body Hohmann transfer solution\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.0 * 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.0 * 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.0 * dv1);

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

fprintf('total delta-v                     %10.4f meters/second \n\n\n', 1000.0 * (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.0 * vt1);

fprintf('transfer orbit apogee velocity    %10.4f meters/second \n\n', 1000.0 * 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);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% state vectors and orbital elements of the two-body Hohmann transfer
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

ttransfer = 0.5 * taut;

tof = ttransfer;

% state vector prior to initial impulse

oevi(1) = req + alt1;
oevi(2) = 0.0;
oevi(3) = inc1;
oevi(4) = 0.0;
oevi(5) = dtr * raan1;
oevi(6) = 0.0;

[ri, vi] = orb2eci(mu, oevi);

fprintf('\norbital elements and state vector of the initial orbit');
fprintf('\n------------------------------------------------------\n');

oeprint1(mu, oevi, 1);

svprint(ri, vi);

% state vector after initial impulse

oevti(1) = smat;
oevti(2) = ecct;
oevti(3) = inc1 - dinc1;
oevti(4) = 0.0;
oevti(5) = dtr * raan1;
oevti(6) = 0.0;

[rti, vti] = orb2eci(mu, oevti);

fprintf('\norbital elements and state vector of the transfer orbit after the first impulse');
fprintf('\n-------------------------------------------------------------------------------\n');

oeprint1(mu, oevti, 1);

svprint(rti, vti);

% state vector prior to second impulse

oevtf(1) = smat;
oevtf(2) = ecct;
oevtf(3) = inc1 - dinc1;
oevtf(4) = 0.0;
oevtf(5) = dtr * raan1;
oevtf(6) = pi;

[rtf, vtf] = orb2eci(mu, oevtf);

fprintf('\norbital elements and state vector of the transfer orbit prior to second impulse');
fprintf('\n-------------------------------------------------------------------------------\n');

oeprint1(mu, oevtf, 1);

svprint(rtf, vtf);

% state vector after second impulse

oevf(1) = req + alt2;
oevf(2) = 0.0;
oevf(3) = inc2;
oevf(4) = 0.0;
oevf(5) = dtr * raan1;
oevf(6) = pi;

[rf, vf] = orb2eci(mu, oevf);

fprintf('\norbital elements and state vector of the final orbit');
fprintf('\n----------------------------------------------------\n');

oeprint1(mu, oevf, 1);

svprint(rf, vf);

dvi = (vti - vi)';

dvf = (vf - vtf)';

upeci = dvi / norm(dvi);

[pitch1, yaw1] = ueci2angles(ri, vi, upeci);

upeci = dvf / norm(dvf);

[pitch2, yaw2] = ueci2angles(rtf, vtf, upeci);

fprintf('\ninitial delta-v vector, magnitude and steering angles');
fprintf('\n-----------------------------------------------------\n');

fprintf('\nx-component of delta-v      %12.6f  meters/second', 1000.0 * dvi(1));

fprintf('\ny-component of delta-v      %12.6f  meters/second', 1000.0 * dvi(2));

fprintf('\nz-component of delta-v      %12.6f  meters/second', 1000.0 * dvi(3));

fprintf('\n\ndelta-v magnitude           %12.6f  meters/second\n', 1000.0 * norm(dvi));

fprintf('\npitch angle                 %12.6f  degrees', rtd * pitch1);

fprintf('\nyaw angle                   %12.6f  degrees\n', rtd * yaw1);

fprintf('\nfinal delta-v vector, magnitude and steering angles');
fprintf('\n---------------------------------------------------\n');

fprintf('\nx-component of delta-v      %12.6f  meters/second', 1000.0 * dvf(1));

fprintf('\ny-component of delta-v      %12.6f  meters/second', 1000.0 * dvf(2));

fprintf('\nz-component of delta-v      %12.6f  meters/second', 1000.0 * dvf(3));

fprintf('\n\ndelta-v magnitude           %12.6f  meters/second\n', 1000.0 * norm(dvf));

fprintf('\npitch angle                 %12.6f  degrees', rtd * pitch2);

fprintf('\nyaw angle                   %12.6f  degrees\n', rtd * yaw2);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% solve the gravity-perturbed Hohmann transfer
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% set "target" values for tpbvp

itarget = inc2;

mee_target = eci2mee(mu, rf, vf);

% initial guess for components of initial delta-v

xg(1) = vti(1) - vi(1);
xg(2) = vti(2) - vi(2);
xg(3) = vti(3) - vi(3);

% initial guess for components of final delta-v

xg(4) = vf(1) - vtf(1);
xg(5) = vf(2) - vtf(2);
xg(6) = vf(3) - vtf(3);

xg = xg';

% define lower and upper bounds for components of delta-v vectors (kilometers/second)

dvm = norm(xg(1:3));

for i = 1:1:3

xlwr(i) = xg(i) - 0.01;

xupr(i) = xg(i) + 0.01;

end

dvm = norm(xg(4:6));

for i = 4:1:6

xlwr(i) = xg(i) - 0.01;

xupr(i) = xg(i) + 0.01;

end

xlwr = xlwr';

xupr = xupr';

% bounds on objective function

flow(1) = 0.0d0;

fupp(1) = +Inf;

% enforce final modified equinoctial equality constraints

flow(2) = 0.0d0;
fupp(2) = 0.0d0;

flow(3) = 0.0d0;
fupp(3) = 0.0d0;

flow(4) = 0.0d0;
fupp(4) = 0.0d0;

flow(5) = 0.0d0;
fupp(5) = 0.0d0;

if (itarget <= 1.0d-8)

% equatorial orbit constraint

flow(6) = 0.0d0;
fupp(6) = 0.0d0;

end

flow = flow';

fupp = fupp';

snspec('snopt_specs.txt');

% solve the orbital TPBVP using SNOPT

snscreen on;

[x, f, inform, xmul, fmul] = snopt(xg, xlwr, xupr, flow, fupp, 'tpbvp');

snprintfile('off');

snsummary('off');

fprintf('\ngravity-perturbed Hohmann transfer solution\n');

fprintf('\norbital elements and state vector of the initial orbit');
fprintf('\n------------------------------------------------------\n');

oeprint1(mu, oevi, 1);

svprint(ri, vi);

fprintf('\norbital elements and state vector of the transfer orbit after the initial delta-v');
fprintf('\n---------------------------------------------------------------------------------\n');

oevti = eci2orb1 (mu, xi(1:3), xi(4:6));

oeprint1(mu, oevti, 1);

svprint(xi(1:3), xi(4:6));

fprintf('\norbital elements and state vector of the transfer orbit prior to the final delta-v');
fprintf('\n----------------------------------------------------------------------------------\n');

oevtf = eci2orb1 (mu, yevent(1:3), yevent(4:6));

oeprint1(mu, oevtf, 1);

svprint(yevent(1:3), yevent(4:6));

fprintf('\norbital elements and state vector of the final orbit');
fprintf('\n----------------------------------------------------\n');

[rf, vf] = mee2eci(mu, mee_final);

oevf = eci2orb1 (mu, rf, vf);

oeprint1(mu, oevf, 1);

svprint(rf, vf);

upeci = x(1:3) / norm(x(1:3));

[pitch1, yaw1] = ueci2angles(ri, vi, upeci);

upeci = x(4:6) / norm(x(4:6));

[pitch2, yaw2] = ueci2angles(yevent(1:3), yevent(4:6), upeci);

fprintf('\ninitial delta-v vector, magnitude and steering angles');
fprintf('\n-----------------------------------------------------\n');

fprintf('\nx-component of delta-v      %12.6f  meters/second', 1000.0 * x(1));

fprintf('\ny-component of delta-v      %12.6f  meters/second', 1000.0 * x(2));

fprintf('\nz-component of delta-v      %12.6f  meters/second', 1000.0 * x(3));

fprintf('\n\ndelta-v magnitude           %12.6f  meters/second\n', 1000.0 * norm(x(1:3)));

fprintf('\npitch angle                 %12.6f  degrees', rtd * pitch1);

fprintf('\nyaw angle                   %12.6f  degrees\n', rtd * yaw1);

fprintf('\nfinal delta-v vector, magnitude and steering angles');
fprintf('\n---------------------------------------------------\n');

fprintf('\nx-component of delta-v      %12.6f  meters/second', 1000.0 * x(4));

fprintf('\ny-component of delta-v      %12.6f  meters/second', 1000.0 * x(5));

fprintf('\nz-component of delta-v      %12.6f  meters/second', 1000.0 * x(6));

fprintf('\n\ndelta-v magnitude           %12.6f  meters/second\n', 1000.0 * norm(x(4:6)));

fprintf('\npitch angle                 %12.6f  degrees', rtd * pitch2);

fprintf('\nyaw angle                   %12.6f  degrees\n', rtd * yaw2);

fprintf('\ntotal delta-v               %12.6f  meters/second\n', ...
1000.0 * (norm(x(1:3)) + norm(x(4:6))));

fprintf('\ntransfer time               %12.6f  seconds', tevent);
fprintf('\n                            %12.6f  minutes', tevent / 60.0);
fprintf('\n                            %12.6f  hours\n', tevent / 3600.0);

fprintf('\ndegree of gravity model    %3i', lgrav);

fprintf('\norder of gravity model     %3i \n\n', mgrav);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% create trajectory graphics %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% compute orbital periods of the initial and final orbits

period1 = 2.0 * pi * oevi(1) * sqrt(oevi(1) / mu);

period3 = 2.0 * pi * oevf(1) * sqrt(oevf(1) / mu);

deltat1 = period1 / 300;

simtime1 = -deltat1;

deltat3 = period3 / 300;

simtime3 = -deltat3;

for i = 1:1:301

simtime1 = simtime1 + deltat1;

simtime3 = simtime3 + deltat3;

% compute initial orbit "normalized" position vector

[rwrk, vwrk] = twobody2 (mu, simtime1, ri, vi);

rp1_x(i) = rwrk(1) / req;

rp1_y(i) = rwrk(2) / req;

rp1_z(i) = rwrk(3) / req;

% compute final orbit position vector

[rwrk, vwrk] = twobody2 (mu, simtime3, rf, vf);

rp3_x(i) = rwrk(1) / req;

rp3_y(i) = rwrk(2) / req;

rp3_z(i) = rwrk(3) / req;

end

% define ode45 options

options = odeset('RelTol', 1.0e-8, 'AbsTol', 1.0e-10);

% create graphics data for transfer trajectory

[twrk, ysol] = ode45(@ceqm1, [0 tevent], [xi(1:3), xi(4:6)], options);

rwrk = ysol(:, 1:3) / req;

figure(1);

% create axes vectors

xaxisx = [1 1.5];
xaxisy = [0 0];
xaxisz = [0 0];

yaxisx = [0 0];
yaxisy = [1 1.5];
yaxisz = [0 0];

zaxisx = [0 0];
zaxisy = [0 0];
zaxisz = [1 1.5];

hold on;

grid on;

% plot earth

[x y z] = sphere(24);

h = surf(x, y, z);

colormap([127/255 1 222/255]);

set (h, 'edgecolor', [1 1 1]);

% plot coordinate system axes

plot3(xaxisx, xaxisy, xaxisz, '-g', 'LineWidth', 1);

plot3(yaxisx, yaxisy, yaxisz, '-r', 'LineWidth', 1);

plot3(zaxisx, zaxisy, zaxisz, '-b', 'LineWidth', 1);

% plot initial orbit

plot3(rp1_x, rp1_y, rp1_z, '-r', 'LineWidth', 1.5);

plot3(rp1_x(1), rp1_y(1), rp1_z(1), 'ob');

% plot transfer trajectory

plot3(rwrk(:, 1), rwrk(:, 2), rwrk(:, 3), '-b', 'LineWidth', 1.5);

plot3(rwrk(end, 1), rwrk(end, 2), rwrk(end, 3), 'ob');

% plot final orbit

plot3(rp3_x, rp3_y, rp3_z, '-g', 'LineWidth', 1.5);

xlabel('X coordinate (ER)', 'FontSize', 12);

ylabel('Y coordinate (ER)', 'FontSize', 12);

zlabel('Z coordinate (ER)', 'FontSize', 12);

title('Hohmann Transfer: Initial, Transfer and Final Orbits', 'FontSize', 16);

axis equal;

view(160, 30);

rotate3d on;

print('-depsc', 'phohmann1.eps');

saveas(h, 'phohmann1.fig');```