clear all; clc
% model written by Allison Shaw (contact for assistance: ashaw@umn.edu)
% started April 2021
%   last updated: 10 January 2023
% takes <1 minute to run
%
% generate mechanism figure

tic
%-----PARAMETERS----------------------------------------------------------%
    v = 0.5;      % dispersal variance
    ngens = 2;    % number of generations
    tau = 1;      % number of types to track in the population (from 0.1 to 1)
    lambda = 2;   % density-independent growth rate
    b = 1;        % density-dependence parameter
    sflag = 0;    % without (0) or with (1) spatial sorting
    mu = 0;       % mortality of dispersers
    a_init = 0;   % initial Allee threshold
    a_lo = 0.3;   % Allee threshold (for plotting only)
    a_hi = 0.7;   % Allee threshold (for plotting only)
    p_lo = 0.6;  % test proportion dispersing
    p_hi = 0.9;  % test proportion dispersing
%-----PARAMETERS----------------------------------------------------------%


%-----MODEL SET UP--------------------------------------------------------%
    eps1 = 1e-15;    % threshold for a 'zero' population density
    eps2 = 1e-1;     % threshold for good kernel integration
    ncrit = 0.001;   % threshold for the edge of the population

    nodes = (2^15)+1;      % number of nodes/bins to have in the domain (2^m + 1)
    diameter = 700;        % length of domain/space to simulate across
    radius = diameter/2;

    % set up 1-D spatial domain
    x = linspace(-radius,radius,nodes);          % vector of node locations
    x2 = linspace(-diameter,diameter,2*nodes-1); % vector of node locations, extended domain
    dx = diameter/(nodes-1);                     % distance between nodes

    pvec = linspace(0.1,1,tau); % propensity to disperse
    
    % dispersal kernel
    k = exp(-sqrt((2*(x2).^2)/v))./sqrt(2*v);
    if abs(1-trapz(x2,k)) > eps2; error('bad dispersal kernel integration: increase nodes'); end
    
    % matrix to store population densities across time (dim 1), type (dim 2),
    % and space (dim 3)
    n = zeros(ngens+1,tau,length(x));
    
    % variables to store information
    popedge = NaN(ngens+1,1);
    speed = NaN(ngens,1);
    core_avg = NaN;
    edge_avg = NaN;
%-----MODEL SET UP--------------------------------------------------------%

%set up initial conditions
IC_rad = 0.5; % radius of initial population
IC_dens = 1;  % initial density
temp = find(abs(x) <= IC_rad);
n(1,:,temp) = IC_dens*ones(size(n(1,:,temp))); % set pop with initial conditions
clear temp

Nvia = sum(squeeze(sum(n(1,:,:)))>ncrit); % number viable pop locations
i = 0;

%-----DO SINGLE INTERATION------------------------------------------------%
    i = i+1;
    
    % INITIAL CONDITIONS
    n0 = squeeze(n(i,:,:));  % tau x nodes

    % GROWTH
    ntot = n0; % DON'T SUM THIS SINCE ONLY A SINGLE TYPE
    g = b./(b+ntot); % strength of density dependence
    g(ntot<a_init) = 0;     % account for Allee effect
    
    if sflag==0 % no spatial sorting
        n1 = repmat(lambda*ntot.*g./tau,tau,1); % offspring are evenly distributed
    elseif sflag==1 % with spatial sorting
        n1 = lambda*n0.*repmat(g,tau,1); % offspring inherit parent's type
    else
        error('unrecognized sflag')
    end
    
    % LOW DISPERSAL - USING ONLY DIM OF 1 FOR TAU
    % calculate convolutions
    % only pvec(j) of each type disperse
    % only 1-mu of dispersers survive
    nx = fft_conv(k,(1-mu)*p_lo.*n1);
    % convert to population density
    n2_lo = dx*nx(:,nodes:length(x2))';
    % fix ends
    n2_lo(1) = n2_lo(1)/2;
    n2_lo(nodes) = n2_lo(nodes)/2;
    % add in non-dispersing individuals
    n2_lo = n2_lo + (1-p_lo).*n1;   
    % set any densities that are too low to zero
    temp = find(n2_lo < eps1);
    n2_lo(temp) = zeros(size(n2_lo(temp)));
    clear temp nx

    % HIGH DISPERSAL - USING ONLY DIM OF 1 FOR TAU
    % calculate convolutions
    % only pvec(j) of each type disperse
    % only 1-mu of dispersers survive
    nx = fft_conv(k,(1-mu)*p_hi.*n1);
    % convert to population density
    n2_hi = dx*nx(:,nodes:length(x2))';
    % fix ends
    n2_hi(1) = n2_hi(1)/2;
    n2_hi(nodes) = n2_hi(nodes)/2;
    % add in non-dispersing individuals
    n2_hi = n2_hi + (1-p_hi).*n1;   
    % set any densities that are too low to zero
    temp = find(n2_hi < eps1);
    n2_hi(temp) = zeros(size(n2_hi(temp)));
    clear temp nx
%-----DO SINGLE INTERATION------------------------------------------------%


fs1 = 10;  % axes labels
fs3 = 09;  % axis numbering
lw2 = 1; % fig edges
lw1 = 1.5; % figure lines

width = 8;
height = 6;
xpos = 3;
ypos = 2;
sx = 0.13;
sy = 0.18;
w = 0.84;
he = 0.8;
dy = 0.12;

c1 = [0.2 0.2 0.2]; % color 1
c2 = [0.6 0.6 0.6]; % color 2

figure(3); clf
hh = gcf;
set(hh,'PaperUnits','centimeters');
set(hh,'Units','centimeters');
set(gcf,'Position',[xpos ypos width height])

axes('position',[sx sy w he])
    % plot initial population distribution
    plot(x,squeeze(n(1,1,:)),'k--','LineWidth',lw1);
    hold on
    % plot two alternative next-year population distributions
    plot(x,n2_hi,'color',c1,'LineWidth',lw1);
    plot(x,n2_lo,'color',c2,'LineWidth',lw1);
    % find the population edge under each Allee threshold
    ind_lo = find(n2_hi>a_lo,1,'last');
    ind_hi = find(n2_lo>a_hi,1,'last');
    e_lo = x(ind_lo);
    e_hi = x(ind_hi);
    % plot two potential Allee thresholds
    line([0 e_lo],[a_lo a_lo],'color','k','LineStyle',':')
    line([0 e_hi],[a_hi a_hi],'color','k','LineStyle',':')
    % plot the population edge under each threshold
    line([e_lo e_lo],[0 n2_hi(ind_lo)],'color','k','LineStyle',':')
    line([e_hi e_hi],[0 n2_lo(ind_hi)],'color','k','LineStyle',':')
    ylabel('density','fontsize',fs1)
    xlabel('space','fontsize',fs1)
    axis([0 1.5 0 1])
    set(gca,'XTick',[e_hi e_lo])
    set(gca,'XTickLabel',{'e_{hi}','e_{lo}'})
    set(gca,'YTick',[a_lo a_hi])
    set(gca,'YTickLabel',{'a_{lo}','a_{hi}'})
    legend('n(x,0)','n(x,1) with p=0.9','n(x,1) with p=0.6')
    set(gca,'FontSize',fs3,'LineWidth',lw2,'Fontname', 'Arial');


% Backup previous settings
prePaperType = get(hh,'PaperType');
prePaperPosition = get(hh,'PaperPosition');
prePaperSize = get(hh,'PaperSize');

% Make changing paper type possible
set(hh,'PaperType','<custom>');
% Set the page size and position to match the figure's dimensions
position = get(hh,'Position');
set(hh,'PaperPosition',[0,0,position(3:4)]);
set(hh,'PaperSize',position(3:4));

%print -dpdf -r600 Figure1.pdf
print -djpeg -r600 fig1.jpg
print -depsc fig1.eps

% Restore the previous settings
set(hh,'PaperType',prePaperType);
set(hh,'PaperPosition',prePaperPosition);
set(hh,'PaperSize',prePaperSize);

toc