% written by Allison Shaw, Contact ashaw@umn.edu for assistance.
% started: 6 Oct 2021
% last updated: 28 November 2022
%
% takes 1-2 days to run

clear all
%%---PARAMETERS-----------------------------------------------------------$
lambda = 10;   % time to half change in infection
gamma = 0.8;   % midpoint of survival-v-intensity
alpha_s = 10;   % slope of survival-v-intensity
phi_max = 2;   % maximum per capita fecundity
alpha_f = 0.8; % slope of fecundity-v-intensity
%
ns = 10; % number of stages
pflag = 0; % whether (1) or not (0) to plot time series when ESS not found
%
betavec = [-1 -0.5 0 0.5]; % scalar of infection dynamics at x=1
muvec = linspace(0,0.02,81);   % slope of survival-v-migration
psivec = linspace(0,0.1,81);   % degree of migratory culling (interaction infection/migration)
%%---PARAMETERS-----------------------------------------------------------$

NM1 = NaN(length(psivec),length(muvec),length(betavec));
M1 = NaN(length(psivec),length(muvec),length(betavec));
m_strategy_vec = NaN(ns,length(psivec),length(muvec),length(betavec));
p_by_distance_all = NaN(6,length(psivec),length(muvec),length(betavec));
ESSfound_all = NaN(length(psivec),length(muvec),length(betavec));

tic
for j = 1:length(betavec)
    beta_min = betavec(j);    % scalar of infection dynamics at x=1
    for jj = 1:length(psivec)
        psi = psivec(jj);  % degree of migratory culling (interaction infection/migration)
        for jjj = 1:length(muvec)
            mu = muvec(jjj); % slope of survival-v-migration
            
            [j jj jjj]
            [m_res_strategy,ESSfound,~,~,~,~,p_by_distance,IavgNM,IavgM] = do_search_ESS(lambda,beta_min,gamma,alpha_s,mu,psi,phi_max,alpha_f,ns,pflag);
            
            NM1(jj,jjj,j) = IavgNM(1,:);
            M1(jj,jjj,j) = IavgM(1,:);
            ESSfound_all(jj,jjj,j) = ESSfound;
            
            m_strategy_vec(:,jj,jjj,j) = m_res_strategy;
            p_by_distance_all(:,jj,jjj,j) = p_by_distance;

            
            name = strcat(['sims_sensitivity_lambda=' num2str(lambda) '_gamma=' num2str(gamma) '_alphaS=' num2str(alpha_s) '_phimax=' num2str(phi_max) '_alphaF=' num2str(alpha_f) '.mat']);
            save(name)
        end
    end
end
toc
