Code covered by the BSD License

# Bi-elliptic Transfer Between Coplanar Circular Orbits

### David Eagle (view profile)

02 Nov 2012 (Updated )

Bi-elliptic orbital transfer between two coplanar circular Earth orbits.

bielliptic.m
```% bielliptic.m       July 9, 2013

% three impulse, outer bi-elliptic transfer
% between coplanar circular orbits

% includes three-dimensional orbit graphics
% and graphical primer vector analysis

% Orbital Mechanics with MATLAB

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

clear all;

global mu req ri rf pvi pvdi

dtr = pi / 180.0;

% read astrodynamic and utility constants

om_constants;

% request inputs

clc; home;

fprintf('\nBi-elliptic Orbit Transfer Analysis\n');

while (1)

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

alti = input('? ');

if (alti > 0)
break;
end

end

while (1)

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

altf = input('? ');

if (altf > 0)
break;
end

end

while(1)

fprintf('\n\ntype of intermediate altitude computation \n');

fprintf('\n  <1> optimal\n');

fprintf('\n  <2> user-defined');

fprintf('\n\n selection (1 or 2)\n');

iac_type = input('? ');

if (iac_type == 1 || iac_type == 2)
break;
end

end

if (iac_type == 2)

while (1)

fprintf('\n\nplease input the bi-elliptic altitude (kilometers)\n');

altb = input('? ');

if (altb > 0)
break;
end

end

end

% initial and final orbit radii (kilometers)

ri = alti + req;

rf = altf + req;

if (iac_type == 1)

xmin = rf;

xmax = 100.0 * rf;

options = optimset('TolFun', 1.0e-6, 'TolX', 1.0e-6);

[x, fx, exitflag] = fminbnd('befunc', xmin, xmax, options);

% intermediate orbit radius and altitude (kilometers)

rb = x;

altb = rb - req;

else

% user-defined intermediate altitude

rb = altb + req;

end

% semimajor axes of the elliptical transfer orbits (kilometers)

smat1 = (ri + rb) / 2.0;

smat2 = (rb + rf) / 2.0;

% compute first ellipse radii (kilometers) and eccentricity (non-dimensional)

rpt1 = req + alti;

rat1 = req + altb;

ecct1 = (rat1 - rpt1) / (rat1 + rpt1);

% compute second ellipse radii (kilometers) and eccentricity (non-dimensional)

rpt2 = req + altf;

rat2 = req + altb;

ecct2 = (rat2 - rpt2) / (rat2 + rpt2);

% initial circular orbit velocity (kilometers/second)

vi = sqrt(mu / ri);

% first ellipse perigee velocity (kilometers/second)

vt1a = sqrt((2.0 * mu / ri) - (mu / smat1));

% first ellipse apogee velocity (kilometers/second)

vt1b = sqrt((2.0 * mu / rb) - (mu / smat1));

% second ellipse apogee velocity (kilometers/second)

vt2b = sqrt((2.0 * mu / rb) - (mu / smat2));

% second ellipse perigee velocity (kilometers/second)

vt2c = sqrt((2.0 * mu / rf) - (mu / smat2));

% final circular orbit velocity (kilometers/second)

vf = sqrt(mu / rf);

% compute delta-v contributions (kilometers/second)

dva = abs(vt1a - vi);

dvb = abs(vt2b - vt1b);

dvc = abs(vf - vt2c);

% compute individual and total transfer time (seconds)

tof1 = pi * sqrt(smat1^3 / mu);

tof2 = pi * sqrt(smat2^3 / mu);

tof = tof1 + tof2;

% print results

clc; home;

fprintf('\n\nBi-elliptic Orbit Transfer Analysis');
fprintf('\n-----------------------------------\n\n');

fprintf('initial orbit altitude           %14.4f kilometers \n\n', alti);

fprintf('initial orbit radius             %14.4f kilometers \n\n', req + alti);

fprintf('initial orbit velocity           %14.4f meters/second \n\n\n', 1000.0 * vi);

fprintf('first ellipse perigee altitude   %14.4f kilometers \n\n', alti);

fprintf('first ellipse perigee radius     %14.4f kilometers \n\n', rpt1);

fprintf('first ellipse apogee altitude    %14.4f kilometers \n\n', altb);

fprintf('first ellipse apogee radius      %14.4f kilometers \n\n', rat1);

fprintf('first ellipse perigee velocity   %14.4f meters/second \n\n', 1000.0 * vt1a);

fprintf('first ellipse apogee velocity    %14.4f meters/second \n\n', 1000.0 * vt1b);

fprintf('first ellipse eccentricity       %14.8f \n\n\n', ecct1);

fprintf('second ellipse perigee altitude  %14.4f kilometers \n\n', altf);

fprintf('second ellipse perigee radius    %14.4f kilometers \n\n', rpt2);

fprintf('second ellipse apogee altitude   %14.4f kilometers \n\n', altb);

fprintf('second ellipse apogee radius     %14.4f kilometers \n\n', rat2);

fprintf('second ellipse perigee velocity  %14.4f meters/second \n\n', 1000.0 * vt2c);

fprintf('second ellipse apogee velocity   %14.4f meters/second \n\n', 1000.0 * vt2b);

fprintf('second ellipse eccentricity      %14.8f \n\n\n', ecct2);

fprintf('final orbit altitude             %14.4f kilometers \n\n', altf);

fprintf('final orbit radius               %14.4f kilometers \n\n', req + altf);

fprintf('final orbit velocity             %14.4f meters/second \n\n\n', 1000.0 * vf);

fprintf('first delta-v                    %14.4f meters/second \n\n', 1000.0 * dva);

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

fprintf('third delta-v                    %14.4f meters/second \n\n', 1000.0 * dvc);

fprintf('total delta-v                    %14.4f meters/second \n\n\n', ...
1000.0 * (dva + dvb + dvc));

fprintf('first ellipse transfer time      %14.4f hours \n', tof1 / 3600.0);

fprintf('                                 %14.4f days \n\n', tof1 / 86400.0);

fprintf('second ellipse transfer time     %14.4f hours \n', tof2 / 3600.0);

fprintf('                                 %14.4f days \n\n', tof2 / 86400.0);

fprintf('total transfer time              %14.4f hours \n', tof / 3600.0);

fprintf('                                 %14.4f days \n\n\n', tof / 86400.0);

%%%%%%%%%%%%%%%%%%%%%%%%%%
% create primer graphics %
%%%%%%%%%%%%%%%%%%%%%%%%%%

% orbital elements and state vector of initial circular orbit

oevi(1) = ri;
oev1(2) = 0.0;
oevi(3) = 0.0;
oevi(4) = 0.0;
oevi(5) = 0.0;
oevi(6) = 0.0;

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

% orbital elements and state vector at perigee of the first ellipse

oevti1(1) = smat1;
oevti1(2) = ecct1;
oevti1(3) = 0.0;
oevti1(4) = 0.0;
oevti1(5) = 0.0;
oevti1(6) = 0.0;

[rtp1, vtp1] = orb2eci(mu, oevti1);

% orbital elements and state vector at apogee of the first ellipse

oevtf2(1) = smat1;
oevtf2(2) = ecct1;
oevtf2(3) = 0.0;
oevtf2(4) = 0.0;
oevtf2(5) = 0.0;
oevtf2(6) = 180.0 * dtr;

[rta1, vta1] = orb2eci(mu, oevtf2);

% orbital elements and state vector at apogee of the second ellipse

oevti2(1) = smat2;
oevti2(2) = ecct2;
oevti2(3) = 0.0;
oevti2(4) = 0.0;
oevti2(5) = 0.0;
oevti2(6) = 180.0 * dtr;

[rta2, vta2] = orb2eci(mu, oevti2);

% orbital elements and state vector at perigee of the second ellipse

oevtf2(1) = smat2;
oevtf2(2) = ecct2;
oevtf2(3) = 0.0;
oevtf2(4) = 0.0;
oevtf2(5) = 0.0;
oevtf2(6) = 0.0;

[rtp2, vtp2] = orb2eci(mu, oevtf2);

% orbital elements and state vector of final circular orbit

oevf(1) = rf;
oevf(2) = 0.0;
oevf(3) = 0.0;
oevf(4) = 0.0;
oevf(5) = 0.0;
oevf(6) = 0.0;

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

% compute orbital periods (seconds)

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

period2 = 2.0 * pi * oevti1(1) * sqrt(oevti1(1) / mu);

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

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

deltat1 = period1 / 300;

simtime1 = -deltat1;

deltat2 = 0.5 * period2 / 300;

simtime2 = -deltat2;

deltat3 = 0.5 * period3 / 300;

simtime3 = -deltat3;

deltat4 = period4 / 300;

simtime4 = -deltat4;

for i = 1:1:301

simtime1 = simtime1 + deltat1;

simtime2 = simtime2 + deltat2;

simtime3 = simtime3 + deltat3;

simtime4 = simtime4 + deltat4;

% initial orbit position vector (ER)

[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;

% first transfer orbit position vector (ER)

[rwrk, vwrk] = twobody2 (mu, simtime2, rtp1, vtp1);

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

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

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

% second transfer orbit position vector (ER)

[rwrk, vwrk] = twobody2 (mu, simtime3, rta2, vta2);

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

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

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

% final orbit position vector (ER)

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

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

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

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

end

% 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];

figure (1);

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 first transfer orbit

plot3(rp2_x, rp2_y, rp2_z, '-b', 'LineWidth', 1.5);

plot3(rp2_x(end), rp2_y(end), rp2_z(end), 'ob');

% plot second transfer orbit

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

plot3(rp3_x(end), rp3_y(end), rp3_z(end), 'ob');

% plot final orbit

plot3(rp4_x, rp4_y, rp4_z, '-g', 'LineWidth', 1.5);

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

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

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

title('Bi-elliptic Transfer: Initial, Transfer and Final Orbits', 'FontSize', 16);

axis equal;

view(50, 40);

rotate3d on;

print -depsc -tiff -r300 bielliptic1.eps

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% process first transfer ellipse
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% initial deltav

dvi = (vtp1 - vi)';

% final deltav

dvf = (vta2 - vta1)';

% perform primer vector initialization

pviniz(tof1, rtp1, vtp1, dvi, dvf);

% number of graphic data points

npts = 300;

% create primer vector and derivative data

dt = tof1 / npts;

for i = 1:1:npts + 1

t = (i - 1) * dt;

if (t == 0)

% initial value of primer magnitude and derivative

pvm = norm(pvi);

pvdm = dot(pvi, pvdi) / pvm;

else

% primer vector and derivative magnitudes at time t

[pvm, pvdm] = pvector(rtp1, vtp1, t);

end

time_wrk(i) = t / 86400.0;

pvm_wrk(i) = pvm;

pvdm_wrk(i) = pvdm;

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% process second transfer ellipse
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% initial deltav

dvi = (vta2 - vta1)';

% final deltav

dvf = (vf - vtp2)';

% perform primer vector initialization

pviniz(tof2, rta2, vta2, dvi, dvf);

% number of graphic data points

npts = 300;

% create primer vector and derivative data

dt = tof2 / npts;

for i = 1:1:npts + 1

t = (i - 1) * dt;

if (t == 0)

% initial value of primer magnitude and derivative

pvm = norm(pvi);

pvdm = dot(pvi, pvdi) / pvm;

else

% primer vector and derivative magnitudes at time t

[pvm, pvdm] = pvector(rta2, vta2, t);

end

time_wrk(i + 301) = (tof1 + t) / 86400.0;

pvm_wrk(i + 301) = pvm;

pvdm_wrk(i + 301) = pvdm;

end

figure(2);

hold on;

plot(time_wrk, pvm_wrk, '-r', 'LineWidth', 1.5);

plot(time_wrk(1), pvm_wrk(1), 'or');

plot(time_wrk(301), pvm_wrk(301), 'or');

plot(time_wrk(end), pvm_wrk(end), 'or');

title('Primer Vector Analysis of the Coplanar Bi-elliptic Transfer', 'FontSize', 16);

xlabel('simulation time (days)', 'FontSize', 12);

ylabel('primer vector magnitude', 'FontSize', 12);

grid;

% create eps graphics file with tiff preview

print -depsc -tiff -r300 primer.eps;

% plot behavior of magnitude of primer derivative

figure(3);

hold on;

plot(time_wrk, pvdm_wrk, '-r', 'LineWidth', 1.5);

plot(time_wrk(1), pvdm_wrk(1), 'or');

plot(time_wrk(301), pvdm_wrk(301), 'or');

plot(time_wrk(end), pvdm_wrk(end), 'or');

title('Primer Vector Analysis of the Coplanar Bi-elliptic Transfer', 'FontSize', 16);

xlabel('simulation time (days)', 'FontSize', 12);

ylabel('primer derivative magnitude', 'FontSize', 12);

grid;

% create eps graphics file with tiff preview

print -depsc -tiff -r300 primer_der.eps;
```