%plotting normalized deltaVs for Hohmann and Bi-elliptic. alpha is rc/ra,
%beta is rb/ra. deltaVs are normalized by initial circular orbit velocity
clc; close all; clear all;
%CASE beta greater than alpha
alpha = linspace(1,70,100);
beta = linspace(70,1000,100);
for j = 1:100
for i = 1:100
DeltaV_H(i) = 1/sqrt(alpha(i)) - sqrt(2)*((1-alpha(i))/(sqrt(alpha(i)*(1+alpha(i)))))-1;
DeltaV_BE(i,j) = sqrt(2*(alpha(i)+beta(j))/(alpha(i)*beta(j)))-(1+sqrt(alpha(i)))/sqrt(alpha(i))-(1-beta(j))*(sqrt(2/(beta(j)*(1+beta(j)))));
end
end
figure(1)
plot(alpha, DeltaV_H,'b','LineWidth', 3)
xlabel('\alpha')
ylabel('normalized \DeltaV')
grid on
hold on
plot(alpha, DeltaV_BE(:,1),'r', 'LineWidth', 3)
title('\beta > \alpha')
for j = 2:100
plot(alpha, DeltaV_BE(:,j),'g')
end
legend('Hohmann')
%%%%%%%%%%%%%%%%%%%%%%%%%%%
%CASE beta less than alpha
alpha = linspace(10,70,100);
beta = linspace(1,10,100);
for j = 1:100
for i = 1:100
DeltaV_H(i) = 1/sqrt(alpha(i)) - sqrt(2)*((1-alpha(i))/(sqrt(alpha(i)*(1+alpha(i)))))-1;
DeltaV_BE(i,j) = sqrt(2*(alpha(i)+beta(j))/(alpha(i)*beta(j)))+(1-sqrt(alpha(i)))/sqrt(alpha(i))-(1-beta(j))*(sqrt(2/(beta(j)*(1+beta(j)))));
end
end
figure(2)
plot(alpha, DeltaV_H,'b','LineWidth', 3)
xlabel('\alpha')
ylabel('normalized \DeltaV')
grid on
hold on
plot(alpha, DeltaV_BE(:,1),'r', 'LineWidth', 3)
title('\beta < \alpha')
for j = 2:100
plot(alpha, DeltaV_BE(:,j),'g')
end
legend('Hohmann')
%%%%%%%%%%%%%%%%%%%%%%%%%%%