% plotting R_0 and I, I/N values for voluntary emigration model

clc

%% parameters

beta1 = 0.37;
beta2 = 0.37;
gamma = 0.14;
Lambda1 = 3.58;
Lambda2 = 4.76*10^3;
lambda1 = 0.5;
lambda2 = 0.5;
mu1 = 3.58*10^(-5);
mu2 = 0.05;
phi = 0.97;

step = 0.0001;
omega = 0:step:1;

%% compute R_0

R0 = 1/mu2*sqrt(Lambda2*beta1*beta2*(mu1+omega)*lambda1*lambda2/ ...
    (Lambda1*(gamma+mu1)*(lambda1+mu1)*(lambda2+mu2)));

%plot(omega,R0)

%% compute EE value of E compartment

a = -mu2*omega*(lambda1+mu1)^2*(lambda2+mu2).*(beta2*mu1*lambda1*...
    (mu1+omega)+mu2*omega*(lambda1+mu1)^2*(lambda2+mu2)*(gamma+mu1));

b = -Lambda2*beta1*beta2*mu1^2*lambda1*lambda2*(lambda1+mu1)*...
    (mu1+omega)-2*Lambda2*mu1*mu2^2*omega*(lambda1+mu1)^2*...
    (lambda2+mu2)*(gamma+mu1)-Lambda1*beta2*mu1^2*mu2*lambda1*...
    (lambda1+mu1)*(lambda2+mu2)*(mu1+omega);

c = Lambda1*Lambda2*beta1*beta2*mu1^2*lambda1*lambda2*(mu2+omega)-...
    Lambda1^2*mu1^2*mu2^2*(lambda1+mu1)*(lambda2+mu2)*(gamma+mu1);

D = b.^2-4*a.*c;

Estar = (-b-sqrt(D))./(2*a);

% separate computation when omega=0 because then a=0
Estar(1) = (Lambda1*Lambda2*beta1*beta2*mu1*lambda1*lambda2 - ...
    Lambda1^2*mu2^2*(lambda1+mu1)*(lambda2+mu2)*(gamma+mu1))/...
    (Lambda2*beta1*beta2*mu1*lambda1*lambda2 + ...
    Lambda1*beta2*mu1*mu2*lambda1*(lambda1+mu1)*(lambda2+mu2));

%% compute EE values of other relevant compartments

Nstar = (Lambda1*mu1+(lambda1+mu1)*omega.*Estar)./(mu1*(mu1+omega));

Istar = phi*lambda1*Estar/(gamma+mu1);

Zstar = Lambda2*beta2*mu1*lambda1*lambda2*(mu1+omega).*Estar./...
    (mu2*(lambda2+mu2)*(beta2*mu1*lambda1*(mu1+omega).*Estar+...
    Lambda1*mu1*mu2*(gamma+mu1)+...
    mu2*omega*(lambda1+mu1)*(gamma+mu1).*Estar));

%% graphs of I* and I*/N*

figure(4);
plot(omega,Istar,'k','LineWidth',1.5,'MarkerSize',3)
xlim([0.0001 0.01]);
set(gcf,'Units','Inches','Position',[0, 0, 3, 2.5]);
ax = gca;
xlabel('migration rate $\omega$','Interpreter','latex')
set(gca,'FontSize',10)
ylabel('infectious individuals $I^*$','Interpreter','latex')
set(gca,'TickLabelInterpreter', 'latex');

Iproportion = Istar./Nstar;
figure(5);
plot(omega,Iproportion,'k','LineWidth',1.5,'MarkerSize',3)
xlim([0.001 1]);
set(gcf,'Units','Inches','Position',[0, 0, 3, 2.5]);
ax = gca;
xlabel('migration rate $\omega$','Interpreter','latex')
set(gca,'FontSize',10)
ylabel('proportion infectious $I^*/N^*$','Interpreter','latex')
set(gca,'TickLabelInterpreter', 'latex');


%% compute probability of getting infected

f1 = beta1*Zstar./Nstar;

pgi = f1./(mu1+f1)*phi*lambda1/(mu1+lambda1);

figure(1);
plot(omega,pgi,'k','LineWidth',1.5,'MarkerSize',3)
axis([0 0.02 0.86 0.98]);
set(gcf,'Units','Inches','Position',[0, 0, 3, 2.5]);
ax = gca;
xlabel('migration rate $\omega$','Interpreter','latex')
set(gca,'FontSize',10)
ylabel('probability of infection','Interpreter','latex')
set(gca,'TickLabelInterpreter', 'latex');


%% look at emigration conditions

Cs = 10;
Cb = 0:step:1;
points = 1/step+1; % the number of sample points
heatmap = zeros(points); % will show if migrate or not
for cost=1:points
    for mr=1:points
        if Cb(cost)+Cs*omega(mr)<pgi(mr)
            % it's better to emigrate
            heatmap(mr,cost)=1;
        end
    end
end

figure(2);
Emigration = imagesc(heatmap);
ax = gca;
nullclinescolors = [0.5 0.5 0.5; % not emigrate
    1 1 1]; % emigrate
colormap(nullclinescolors)
colormap(ax);
set(gcf,'Units','Inches','Position',[0, 0, 3, 2.5]);
ax = gca;
set(gca,'FontSize',10);
ax.TickLabelInterpreter = 'latex';
ax.XTick = [0 0.2*points 0.4*points 0.6*points 0.8*points 1*points];
ax.XTickLabel = {'0', '0.2', '0.4', '0.6', '0.8', '1'};
ax.YDir = 'normal'; % flips y-axis
ax.YTick = [0 0.2*points 0.4*points 0.6*points 0.8*points 1*points];
ax.YTickLabel = {'0', '0.2', '0.4', '0.6', '0.8', '1'};
hl = xlabel('relative base cost $\tilde{C}_b$');
set(hl,'interpreter','latex');
hl = ylabel('migration rate $\omega$');
set(hl,'interpreter','latex');

%% mandatory emigration

Ndiff = -(Nstar-Nstar(1));
Ndiff(2:4) = 0;
Cme = 0.00024;

% payoff in the mandatory migration
Eme = -Istar - Cme*Ndiff;
Eme(1:4) = Eme(5);

% omega_me = omega(1:501);
% Eme_me = Eme(1:501);

figure(3);
plot(omega,Eme,'k','LineWidth',1.5,'MarkerSize',3)
xlim([0.0005 0.01]);
ylim([-24 -23.6]);
set(gcf,'Units','Inches','Position',[0, 0, 3, 2.5]);
ax = gca;
xlabel('migration rate $\omega$','Interpreter','latex')
set(gca,'FontSize',10)
ylabel('mandatory emigration payoff','Interpreter','latex')
set(gca,'TickLabelInterpreter', 'latex');




%