hohmann transfer in matlab

main.m
r1 = 7000; % initial orbit radius [km]
r2 = 21000; % final orbit radius [km]
mu = 398600; % gravitational parameter [km^3/s^2]
v1 = sqrt(mu/r1); % initial orbit velocity [km/s]
v2 = sqrt(mu/r2); % final orbit velocity [km/s]
225 chars
6 lines

main.m
a_trans = (r1 + r2)/2; % transfer orbit semi-major axis [km]
delta_v1 = sqrt(mu/r1)*(sqrt(2*r2/(r1+r2)) - 1); % velocity increment at initial orbit [km/s]
delta_v2 = sqrt(mu/r2)*(1 - sqrt(2*r1/(r1+r2))); % velocity increment at final orbit [km/s]
247 chars
4 lines

main.m
T_trans = 2*pi*sqrt(a_trans^3/mu); % transfer orbit period [s]
TOF = T_trans/2; % time of flight [s]
101 chars
3 lines

main.m
theta1 = 0; % true anomaly at initial point [rad]
theta2 = pi; % true anomaly at final point [rad]
99 chars
3 lines

main.m
syms E
M = sqrt(mu/a_trans^3)*E; % mean anomaly [rad]
f = @(E) E - a_trans/r1*sin(E) - M;
E0 = fsolve(f, pi/2); % initial eccentric anomaly [rad]
nu0 = 2*atan(sqrt((1 + e)/2)/sqrt((1 - e)/2)*tan(E0/2)); % initial true anomaly [rad]
E1 = fsolve(f, pi/2); % final eccentric anomaly [rad]
nu1 = 2*atan(sqrt((1 + e)/2)/sqrt((1 - e)/2)*tan(E1/2)); % final true anomaly [rad]

figure('Color', 'w')
hold on
% plot initial orbit
t = linspace(0, 2*pi);
x = r1*cos(t);
y = r1*sin(t);
plot(x, y, 'k--')
% plot final orbit
x = r2*cos(t);
y = r2*sin(t);
plot(x, y, 'k--')
% plot transfer orbit
e = sqrt(1 - (r1+r2)/(2*a_trans));
theta = linspace(0,2*pi,10000);
r = a_trans*(1 - e^2)./(1 + e*cos(theta-nu0));
[x,y] = pol2cart(theta,r);
plot(x,y,'b')
axis equal
legend('Initial orbit', 'Final orbit', 'Transfer orbit')
xlabel('x [km]')
ylabel('y [km]')
title('Hohmann transfer')
864 chars
31 lines