%Produces some game theory plots and pertinant information for the
%migration model.
clear
clc

pop = 100000;
beta1 = 0.37;
beta2 = 0.37;
mu1 = 1/(365.25*76.4);
mu2 = 1/20;
alpha1 = pop*mu1;
phi = 0.97;
R0 = 4.072727272727272727272727272727;
lambda1 = 0.5;
lambda2 = 0.5;
ggamma = 1/7;
alpha2 = (R0^2)*(mu2^2)*alpha1*(ggamma+mu1)*(lambda1+mu1)*(lambda2+mu2)/(mu1*beta1*beta2*lambda1*lambda2);

%Use a range of omega values to compute other values:
omega = 0:0.0001:1;
C = zeros(length(omega),1);
R0var = zeros(length(omega),1);
I2var = zeros(length(omega),1);
I2N2var = zeros(length(omega),1);
for i=1:length(omega)
    %EE:
    a = -mu2*omega(i)*(lambda1+mu1)^2*(lambda2+mu2)*(beta2*mu1*lambda1*(mu1+omega(i))+mu2*omega(i)*(lambda1+mu1)*(ggamma+mu1));
    b = -alpha2*beta1*beta2*mu1^2*lambda1*lambda2*(lambda1+mu1)*(mu1+omega(i))-alpha1*mu1*mu2*(lambda1+mu1)*(lambda2+mu2)*(beta2*mu1*lambda1*(mu1+omega(i))+mu2*omega(i)*(lambda1+mu1)*(ggamma+mu1))-alpha1*mu1*mu2^2*omega(i)*(lambda1+mu1)^2*(lambda2+mu2)*(ggamma+mu1);
    c = alpha1*alpha2*beta1*beta2*mu1^2*lambda1*lambda2*(mu1+omega(i))-alpha1^2*mu1^2*mu2^2*(lambda1+mu1)*(lambda2+mu2)*(ggamma+mu1);
    E2 = (-b - sqrt(b^2-4*a*c))/(2*a);
    S2 = (1/(mu1+omega(i)))*(alpha1-(lambda1+mu1)*E2);
    I2 = phi*lambda1*E2/(ggamma+mu1);
    Ia2 = (1-phi)*lambda1*E2/(ggamma+mu1);
    R2 = (ggamma/mu1)*(I2+Ia2);
    N2 = (alpha1-omega(i)*S2)/mu1;
    X2 = alpha2*(ggamma+mu1)*N2/(beta2*lambda1*E2+mu2*(ggamma+mu1)*N2);
    Y2 = (alpha2-mu2*X2)/(lambda2+mu2);
    Z2 = (lambda2/mu2)*Y2;

    R0 = (1/mu2)*sqrt(alpha2*beta1*beta2*lambda1*lambda2*(mu1+omega(i))/(alpha1*(lambda1+mu1)*(lambda2+mu2)*(ggamma+mu1)));
    if R0 < 1
        %Adjust equilibrium to DFE if necessary:
        Z2 = alpha2/mu2;
        N2 = alpha1/(mu1+omega(i));
        Rtracker = [Rtracker; omega(i)];
    end
    f = beta1*Z2/N2;
    %Sets C=pi
    C(i) = f/(mu1+omega(i)+f);
    %Keep track of R0, I*, and I*/N*:
    R0var(i) = R0;
    I2var(i) = I2;
    I2N2var(i) = I2/N2;
end

figure
box on
hold on
set(gcf,'Units','Inches','Position',[0, 0, 3, 2.5]);
plot(C,omega,'k','LineWidth',1.5,'MarkerSize',3)
%plot(linspace(0,1),rHI*ones(length(linspace(0,1))),'k--','LineWidth',1,'MarkerSize',1)
xlabel('relative cost of protection $C$','Interpreter','latex')
%xticks([0 0.2 0.4 0.6 0.8 1])
%xticklabels({'$0$','$0.2$','$0.4$','$0.6$','$0.8$','$1$'})
set(gca, 'FontSize', 10)
ylabel('migration rate $\omega$','Interpreter','latex')
%yticks([0 0.2 0.4 0.6 0.8 1])
%yticklabels({'$0$','$0.2$','$0.4$','$0.6$','$0.8$','$1$'})
set(gca,'TickLabelInterpreter', 'latex');
%plot([max(max(C)) 1],[0 0],'k','LineWidth',1.5,'MarkerSize',3)


figure
box on
%hold on
set(gcf,'Units','Inches','Position',[0, 0, 3, 2.5]);
semilogy(omega,R0var,'k','LineWidth',1.5,'MarkerSize',3)
axis([0 1 0 1000]);
%plot(linspace(0,1),rHI*ones(length(linspace(0,1))),'k--','LineWidth',1,'MarkerSize',1)
xlabel('migration rate $\omega$','Interpreter','latex')
%xticks([0 0.2 0.4 0.6 0.8 1])
%xticklabels({'$0$','$0.2$','$0.4$','$0.6$','$0.8$','$1$'})
set(gca, 'FontSize', 10)
ylabel('basic reproduction number $R_0$','Interpreter','latex')
%yticks([0 0.2 0.4 0.6 0.8 1])
%yticklabels({'$0$','$0.2$','$0.4$','$0.6$','$0.8$','$1$'})
set(gca,'TickLabelInterpreter', 'latex');
%plot([max(max(C)) 1],[0 0],'k','LineWidth',1.5,'MarkerSize',3)

figure
box on
hold on
set(gcf,'Units','Inches','Position',[0, 0, 3, 2.5]);
plot(omega,I2var,'k','LineWidth',1.5,'MarkerSize',3)
%plot(linspace(0,1),rHI*ones(length(linspace(0,1))),'k--','LineWidth',1,'MarkerSize',1)
xlabel('migration rate $\omega$','Interpreter','latex')
%xticks([0 0.2 0.4 0.6 0.8 1])
%xticklabels({'$0$','$0.2$','$0.4$','$0.6$','$0.8$','$1$'})
set(gca, 'FontSize', 10)
ylabel('infected humans $I^*$','Interpreter','latex')
%yticks([0 0.2 0.4 0.6 0.8 1])
%yticklabels({'$0$','$0.2$','$0.4$','$0.6$','$0.8$','$1$'})
set(gca,'TickLabelInterpreter', 'latex');
%plot([max(max(C)) 1],[0 0],'k','LineWidth',1.5,'MarkerSize',3)

figure
box on
hold on
set(gcf,'Units','Inches','Position',[0, 0, 3, 2.5]);
plot(omega,I2N2var,'k','LineWidth',1.5,'MarkerSize',3)
%plot(linspace(0,1),rHI*ones(length(linspace(0,1))),'k--','LineWidth',1,'MarkerSize',1)
xlabel('migration rate $\omega$','Interpreter','latex')
%xticks([0 0.2 0.4 0.6 0.8 1])
%xticklabels({'$0$','$0.2$','$0.4$','$0.6$','$0.8$','$1$'})
set(gca, 'FontSize', 10)
ylabel('\% of infected humans $I^*/N^*$','Interpreter','latex')
%yticks([0 0.2 0.4 0.6 0.8 1])
%yticklabels({'$0$','$0.2$','$0.4$','$0.6$','$0.8$','$1$'})
set(gca,'TickLabelInterpreter', 'latex');
%plot([max(max(C)) 1],[0 0],'k','LineWidth',1.5,'MarkerSize',3)