%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')
%%%%%%%%%%%%%%%%%%%%%%%%%%%
