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.

tpbvp(x)
```function [f, g] = tpbvp(x)

% two point boundary value objective function and
% and modified equinoctial element constraints

% input

%  x = current delta-v vectors
%      x(1:3) = first delta-v vector
%      x(4:6) = second delta-v vector

% output

%  f(1) = objective function (delta-v magnitude)
%  f(2) = semiparameter equality constraint
%  f(3) = f element equality constraint
%  f(4) = g element equality constraint
%  f(5) = h^2 + k^2 equality constraint (non-equatorial mission orbit)
%       = h equality constraint (equatorial mission orbit)
%  f(6) = k equality constraint (equatorial mission orbit)

% Orbital Mechanics with MATLAB

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

global mu ri vi xi tof itarget

global tevent yevent mee_target mee_final

% load current state vector of transfer orbit after first impulse

xi(1) = ri(1);
xi(2) = ri(2);
xi(3) = ri(3);

xi(4) = vi(1) + x(1);
xi(5) = vi(2) + x(2);
xi(6) = vi(3) + x(3);

% set up options for ode45

options = odeset('RelTol', 1.0e-12, 'AbsTol', 1.0e-12, 'Events', @nc_event);

% solve for nodal crossing condition

rwrk = xi(1:3);

vwrk = xi(4:6);

% maximum search duration = 102% of two-body transfer time (seconds)

tend = 1.02 * tof;

[t, ysol, tevent, yevent, ie] = ode45(@ceqm1, [0 tend], [rwrk vwrk], options);

% objective function (total delta-v magnitude)

f(1) = norm(x(1:3)) + norm(x(4:6));

% current position vector at nodal crossing

rf = yevent(1:3);

% current velocity vector after second delta-v

vf(1) = yevent(4) + x(4);

vf(2) = yevent(5) + x(5);

vf(3) = yevent(6) + x(6);

% compute current modified equinoctial orbital elements of mission orbit

mee_final = eci2mee(mu, rf, vf);

% enforce semiparameter and circular orbit constraints (p, f = g = 0)

f(2) = mee_final(1) - mee_target(1);

f(3) = mee_final(2) - mee_target(2);

f(4) = mee_final(3) - mee_target(3);

% enforce mission orbit inclination constraint(s)

if (itarget <= 1.0d-8)

% equatorial orbit (h = k = 0)

f(5) = mee_final(4);

f(6) = mee_final(5);

else

% non-equatorial mission orbit (h^2 + k^2 = 0)

f(5) = (mee_final(4)^2 + mee_final(5)^2) - (mee_target(4)^2 + ...
mee_target(5)^2);

end

% transpose objective function/constraints vector

f = f';

% no derivatives

g = [];
```