% Parameter values
beta1 = 1;
beta2 = 0.4;
h1 = 0.5;
h2 = 0.5;
a1 = 1;
a2 = 2;
c1 = 20*h1;
c2 = 20*h2;

% Functional response of Model 2

F = @(N1, N2) a1*beta1*(1+c1*beta1*N1).*N1.^2 ./ ...
    ( a1*N1.*(1+c1*beta1*N1).*(1+beta1*h1*N1) + a2*N2.*(1+c2*beta2*N2).*(1+beta2*h2*N2) );

n = 100;
N1max = 15;
N2max = 15;
N1 = linspace(1e-9, N1max, n);
N2 = linspace(1e-9, N2max, n);

[NN1, NN2] = meshgrid(N1, N2);

%%%%%%%%%%%%%%%%%%%%%%%%%
% 3D figure, full range %
%%%%%%%%%%%%%%%%%%%%%%%%%

figure(1)
clf
ax = axes;
h = surf(NN1, NN2, F(NN1, NN2));
ax.YDir = 'reverse';
h.EdgeColor = 'none';

hold on
spacing = 10;
for i = 1:spacing:length(NN1(:,1))
    plot3(NN1(:, i), NN2(:, i), F(NN1(:, i), NN2(:, i)), ':k');
    plot3(NN1(i, :), NN2(i, :), F(NN1(i, :), NN2(i, :)), ':k');
end

h.FaceColor = 'none';
h.EdgeAlpha = '1';
h.LineStyle = ':';
set(gca, 'xtick',[], 'ytick', [],'ztick',[]);
ax.LineWidth = 2;
ax.GridColor = 'none';

hold on
plot3([0 N1max N1max 0 0], [0 0 N2max N2max 0], [0 0 0 0 0], '-k', 'LineWidth', 2);
h1 = plot3(N1, N2, F(N1, N2), '-k', 'LineWidth', 1.5);
h2 = plot3(N1, N2max - N2, F(N1, N2max - N2), '-.k', 'LineWidth', 1.5);
h3 = plot3(N1, repmat(2, size(N1)), F(N1, 2), '--k', 'LineWidth', 1.5);

view([-10 10])
exportgraphics(ax,'M2_surf_A.png','Resolution',200) 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%      2D figure, full range     %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

figure(2)
clf
ax = axes;
h = surf(NN1, NN2, F(NN1, NN2));
h.EdgeColor = 'none';
ax.YDir = 'reverse';
ax.LineWidth = 2;
ax.GridColor = 'none';

hold on
spacing = 10;
for i = 1:spacing:length(NN1(:,1))
    plot3(NN1(:, i), NN2(:, i), F(NN1(:, i), NN2(:, i)), ':k');
    plot3(NN1(i, :), NN2(i, :), F(NN1(i, :), NN2(i, :)), ':k');
end

h.FaceColor = 'none';
h.EdgeAlpha = '1';
h.LineStyle = ':';

hold on
plot3([0 N1max N1max 0 0], [0 0 N2max N2max 0], [0 0 0 0 0], '-k', 'LineWidth', 2);
h1 = plot3(N1, N2, F(N1, N2), '-k', 'LineWidth', 1.5);
h2 = plot3(N1, N2max - N2, F(N1, N2max - N2), '-.k', 'LineWidth', 1.5);
h3 = plot3(N1, repmat(2, size(N1)), F(N1, 2), '--k', 'LineWidth', 1.5);

box on
set(gca,'xtick',[],'ytick',[],'ztick',[]);
view([0 0])

exportgraphics(ax,'M2_surf_B_full.png','Resolution',200) 

%%%%%%%%%%%%%%%%%%%%%
% 2D, Zoomed figure %
%%%%%%%%%%%%%%%%%%%%%

N1max = 4;
N2max = 15;
N1 = linspace(1e-9, N1max, n);
N2 = linspace(1e-9, N2max, n);
[NN1, NN2] = meshgrid(N1, N2);

figure(3)
clf
ax = axes;
h = surf(NN1, NN2, F(NN1, NN2));
ax.YDir = 'reverse';
h.EdgeColor = 'none';

hold on
spacing = 10;
for i = 1:spacing:length(NN1(:,1))
    plot3(NN1(:, i), NN2(:, i), F(NN1(:, i), NN2(:, i)), ':k');
    plot3(NN1(i, :), NN2(i, :), F(NN1(i, :), NN2(i, :)), ':k');
end

h.FaceColor = 'none';
h.EdgeAlpha = '1';
h.LineStyle = ':';

ax.LineWidth = 2;
ax.GridColor = 'none';

set(gca, 'xtick',[], 'ytick', [],'ztick',[]);
hold on
plot3([0 N1max N1max 0 0], [0 0 N2max N2max 0], [0 0 0 0 0], '-k', 'LineWidth', 2);
h1 = plot3(N1, N1, F(N1, N1), '-k', 'LineWidth', 1.5);
h2 = plot3(N1, N2max - N1, F(N1, N2max - N1), '-.k', 'LineWidth', 1.5);
h3 = plot3(N1, repmat(2, size(N1)), F(N1, 2), '--k', 'LineWidth', 1.5);

view([0 0])

exportgraphics(ax, 'M2_surf_B_zoom.png', 'Resolution', 200) 