close all

constants
mu = G*M_E;

h = 300e3;
r1 = R_E+h;
step = 0.1;
start = 0.1;
R = start:step:100;
r2 = r1*R;

v1 = visviva_circ(r1);

%Plot Hohmann w.r.t R
deltaVtot = sqrt(mu/r1).*(sqrt(2*r2./(r1+r2)) - 1) + sqrt(mu./r2).*(1 - sqrt(2*r1./(r1+r2))) % Total deltaV hohmann.
deltaVnorm = deltaVtot/v1;

plot(R, deltaVnorm,'k')

grid on
xlabel('$R = r_2/r_1$','Interpreter','latex')
ylabel('$\Delta v/v_1$','Interpreter','latex')

legend("Hohmann")
legend('show')

hold on


% Bi-elliptical w.r.t R for range of alpha
alpha_arr = [5 11.94 15.58 20 30 50 100]
for alpha = alpha_arr
    rb = r1*alpha;

    a1 = (r1+rb)/2;
    a2 = (r2+rb)./2;

    deltaV1 = sqrt(2*mu/r1 - mu/a1) - sqrt(mu/r1);
    deltaV2 = sqrt(2*mu/rb - mu./a2) - sqrt(2*mu/rb - mu/a1);
    deltaV3 = sqrt(2.*mu./r2 - mu./a2) - sqrt(mu./r2);

    deltaVtot = deltaV1 + deltaV2

    for i = 1:length(R)
        if R(i) <= alpha
            deltaVtot(i) = deltaVtot(i) + deltaV3(i);
        else
            deltaVtot(i) = deltaVtot(i) - deltaV3(i);
        end
    end

    deltaVnorm = deltaVtot/v1

    plot(R, deltaVnorm)
    ylim([0 1])
end

% Plot a=infinity theoretical limit.
deltaVinf = sqrt(mu./r1).*(sqrt(2)-1).*(1+sqrt(r1./r2));
plot(R, deltaVinf./v1)

legend(['Hohmann' string(alpha_arr) 'infty'])

