% Script to reproduce plots from Proc Roy Soc submission, September 2024

close all
clear all

% Figure settings:
set(0, 'defaultaxesfontsize', 20);
set(0, 'defaultaxeslinewidth', 2);
set(0, 'defaultlinelinewidth', 2);
set(0, 'defaultpatchlinewidth', 0.7);

%% Figure 3
figure(1)
a = load('figure3data');

set(0, 'DefaultAxesFontSize', 14)

corners = 0.05;
width = 0.75;
starty2 = 0.15;
height = (1-2*starty2-corners)/2;
starty1 = starty2 + height + starty2;
startx = (1-width)/2;

subplot('Position',[startx,starty1,width,height]);
surf(a.X,1-a.Y,a.Z,'edgecolor','none'); view([0,90]);
hold on;
plot3(a.X,1-ones(size(a.X))*a.Y(1*length(a.Y)/5),ones(size(a.X)),'b','LineWidth',0.5);
hold off;
xlim([min(a.X),0.3]);
colormap gray
xlabel('$y/L$','interpreter','latex');
ylabel('$x/L$','interpreter','latex');

insetx1 = 0.192;
insetx2 = 0.207;
xmax = 0.3;

rectangle('Position', [insetx1, 0, insetx2-insetx1, 1], ...
    'Curvature', 0.0, ...
    'FaceColor', [0, 0, 0, 1], ...
    'EdgeColor', [1, 0, 0, 0.7], ...
    'LineWidth',2);

h = axes('Position',[startx+(insetx1/xmax)*width,starty1...
    ,width*(insetx2-insetx1)/xmax,height]);
surf(a.X,1-a.Y,a.Z,'edgecolor','none'); view([0,90]);
hold on;
plot3(a.X,1-ones(size(a.X))*a.Y(1*length(a.Y)/5),ones(size(a.X)),'b','LineWidth',0.5);
hold off;
xlim([insetx1,insetx2]);
colormap gray
clim([-1e-6,1e-6])
xticks([]);
yticks([]);
h.XAxis.Visible = 'off';
h.YAxis.Visible = 'off';

subplot('Position',[startx,starty2,width,height]);
plot(a.X,a.Z(1*length(a.Y)/5,:),'b');
xlabel('$y/L$','interpreter','latex');
ylabel('$w$','interpreter','latex');

%% Figure 4

figure(2)
set(0, 'DefaultAxesFontSize', 17)
set(0, 'DefaultLineLineWidth', 1.5);
cols = parula(6);

width = 0.28;
height = 0.2;
left = 0.07;
h1 = 0.5;
h2 = h1 - 1.4*height;

% Tension 1
c = load('figure4profile1');
d = load('figure4fourier1.mat');
f = d.f; x = d.x; n = double(d.n);

subplot('position',[left h1 width 2*height])
surf(c.X,1+c.Y/c.width,c.Z,'edgecolor','none'); view([0,90]); hold on;
surf(c.X,c.Y/c.width,c.Z*1e5,'edgecolor','none'); hold off;
caxis([min(min(c.Z)),max(max(c.Z))]); colormap gray
xlim([0,c.length]); xticks([]);
ylabel('$y$','interpreter','latex');
title('$T/Eh = 1.65 \times 10^{-5}$','interpreter','latex');

subplot('position',[left h2 width height])
for i = 1:length(d.n)
    F = f{i}; X = x{i};
    plot(X/d.length,n(i)^2*F.^2,'color',cols(1+mod(i-1,length(cols)),:)); hold on;
end

xlim([0,1]); ylim([0,5e-3]);
ax = gca; ax.Clipping = 'off'; set(gca,'TickDir','out');
xlabel('$x/L$','interpreter','latex');
ylabel('$n^2 C_n^2$','interpreter','latex');
hold off;

% Tension 2
left = 0.39;
c = load('figure4profile2');
d = load('figure4fourier2');
f = d.f; x = d.x; n = double(d.n);
subplot('position',[left h1 width 2*height])
surf(c.X,1+c.Y/c.width,c.Z,'edgecolor','none'); view([0,90]); hold on;
surf(c.X,c.Y/c.width,c.Z*1e5,'edgecolor','none'); hold off;
caxis([min(min(c.Z)),max(max(c.Z))]); colormap gray
xlim([0,c.length]); xticks([]); yticks([]);
title('$T/Eh = 3.3 \times 10^{-5}$','interpreter','latex');

subplot('position',[left h2 width height])
for i = 1:length(d.n)
    F = f{i}; X = x{i};
    plot(X/d.length,n(i)^2*F.^2,'color',cols(1+mod(i-1,length(cols)),:)); hold on;
end

xlim([0,1]); ylim([0,5e-3]); yticks([]);
ax = gca; ax.Clipping = 'off'; set(gca,'TickDir','out');
xlabel('$x/L$','interpreter','latex');
hold off;

% Tension 3
left = 0.71;
c = load('figure4profile3.mat');
d = load('figure4fourier3.mat');
f = d.f; x = d.x; n = double(d.n);
subplot('position',[left h1 width 2*height])
surf(c.X,1+c.Y/c.width,c.Z,'edgecolor','none'); view([0,90]); hold on;
surf(c.X,c.Y/c.width,c.Z*1e5,'edgecolor','none'); hold off;
caxis([min(min(c.Z)),max(max(c.Z))]); colormap gray
xlim([0,c.length]); xticks([]); yticks([]);
title('$T/Eh = 6.6 \times 10^{-5}$','interpreter','latex');

subplot('position',[left h2 width height])
for i = 1:length(d.n)
    F = f{i}; X = x{i};
    plot(X/d.length,n(i)^2*F.^2,'color',cols(1+mod(i-1,length(cols)),:)); hold on;
end

xlim([0,1]); ylim([0,5e-3]); yticks([]);
ax = gca; ax.Clipping = 'off'; set(gca,'TickDir','out');
xlabel('$x/L$','interpreter','latex');
hold off;

%% Figure 5a

figure(3)

set(0, 'DefaultAxesFontSize', 17)
set(0, 'DefaultLineLineWidth', 1.5);

% Tension 1
d = load('figure4fourier1.mat');
f = d.f; x = d.x; n = double(d.n); L = d.length;
X = linspace(0,d.length,1000);

xvals = zeros(length(d.n)-2,1);
for i = 1:length(d.n)-2
    F0 = spline(x{i},f{i},X);
    F1 = spline(x{i+1},f{i+1},X);
    xpts = find((n(i+1)*F1>n(i)*F0).*(n(i)^2*F0.^2>1e-3));
    xvals(i) = X(xpts(1));
end

Npref = @(X) 1./((1/3)*(1-(7/8)*(X/L)));
lambda = @(X) d.width./Npref(X);
xdiffs = xvals(2:end)-xvals(1:end-1);
lambdas = lambda(0.5*(xvals(1:end-1)+xvals(2:end)));
loglog(lambdas(2:end)/5.2,xdiffs(2:end)/L,'-x'); hold on;

% Tension 2
d = load('figure4fourier2');
f = d.f; x = d.x; n = double(d.n);

xvals = zeros(length(d.n)-2,1);
for i = 1:length(d.n)-2
    F0 = spline(x{i},f{i},X);
    F1 = spline(x{i+1},f{i+1},X);
    xpts = find((n(i+1)*F1>n(i)*F0).*(n(i)^2*F0.^2>1e-3));
    xvals(i) = X(xpts(1));
end

Npref = @(X) 1./((1/3)*(1-(7/8)*(X/L)));
lambda = @(X) d.width./Npref(X);
xdiffs = xvals(2:end)-xvals(1:end-1);
lambdas = lambda(0.5*(xvals(1:end-1)+xvals(2:end)));
loglog(lambdas(2:end)/5.2,xdiffs(2:end)/L,'-x');

% Tension 3
d = load('figure4fourier3.mat');
f = d.f; x = d.x; n = double(d.n);

xvals = zeros(length(d.n)-2,1);
for i = 1:length(d.n)-2
    F0 = spline(x{i},f{i},X);
    F1 = spline(x{i+1},f{i+1},X);
    xpts = find((n(i+1)*F1>n(i)*F0).*(n(i)^2*F0.^2>1e-3));
    xvals(i) = X(xpts(1));
end

Npref = @(X) 1./((1/3)*(1-(7/8)*(X/L)));
lambda = @(X) d.width./Npref(X);
xdiffs = xvals(2:end)-xvals(1:end-1);
lambdas = lambda(0.5*(xvals(1:end-1)+xvals(2:end)));
loglog(lambdas(2:end)/5.2,xdiffs(2:end)/L,'-x');

pre = 2.35;
X = linspace(0.5,0.8)/5.2;
Y = pre*X.^2;
plot(X,Y,'k-');
plot(X,min(Y)*ones(size(X)),'k-');
plot(max(X)*[1,1],[min(Y),max(Y)],'k-');
hold off;

legend('low tension','medium tension','high tension','location','south east');
xlim([2e-1,2e-0]/5.2); ylim([2e-1,2e1]/80);

text(0.65/5.2,1.5/L,'1');
text(0.85/5.2,2.5/L,'2');

xlabel('$\lambda_{\mathrm{loc}}/W$','interpreter','latex');
ylabel('$d/L$','interpreter','latex');
set(gca,'LineWidth',2)

%% Figure 5b

figure(4)

string = ["figure4fourier1.mat",...
    "figure4fourier2.mat",...
    "figure4fourier3.mat"];
L = 80; nx = 1e4;

for S = 1:length(string)
    a = load(string(S));     
    x = a.x; f = a.f;
    X = linspace(x{1}(1),x{1}(end),nx);
    F = zeros(length(x),length(X));
    nvals = 2 + (1:length(x));
    wrlength = zeros(size(X));

    for i = 1:length(x)
        F(i,:) = spline(x{i},f{i},X);
        wrlength = wrlength + (nvals(i)*F(i,:)).^2;
    end

    N = zeros(size(X));

    for j = 1:length(X)
        [~,m] = max(nvals'.*F(:,j));
        N(j) = nvals(m);
    end

    plot(X/L,N,'LineWidth',1); hold on;

end

Npref = 1./((1/3)*(1-(7/8)*(X/max(X))));
plot(X/L,Npref,'k--','LineWidth',1);
hold off;

xlabel('$x$','interpreter','latex');
ylabel('$n$','interpreter','latex');
set(gca,'fontsize',16)



%% Figure 6

alphas = [0.1 0.3 0.5 1 4.5 10];

for i = 1:length(alphas)

    figure

table = load(strcat('profilespaper/profilespaper',num2str(alphas(i)),'.mat'));
table = table.tab;

subplot('Position',[0.11 0.78 0.78 0.15])
        plot(table.table1.X,table.table2.kappa);
        xlim([-5,5]); ylim([-2,2]);
        %xlabel('$x$','interpreter','latex');
        ylabel('$\kappa$','interpreter','latex');
        set(gca,'xticklabel',[])

        title(strcat('$\alpha = ',num2str(alphas(i)),'$'),'interpreter','latex');

subplot('Position',[0.11 0.57 0.78 0.15])
        plot(table.table2.X,table.table2.Z);
        
        ylabel('$Z_k$','interpreter','latex');
        set(gca,'xticklabel',[])
        %title(['$\alpha$ = ',num2str(alphas(i))],'interpreter','latex');
xlim([-5,5])

y = table.table3.Y;
subplot('Position',[0.11 0.13 0.78 0.39])
        surf(table.table3.X,table.table3.Y,(table.table3.w)','EdgeColor','none'); view([0,90]);
        %xlabel('$x$','interpreter','latex');
        xlabel('$\bar{X}$','interpreter','latex');
        ylabel('$y/W$','interpreter','latex');
        
        set(gca,'FontSize',15)
        ylim([min(min(y)),max(max(y))]);
    xlim([-5,5])

end

%% Figure 8

figure

set(0, 'defaultaxeslinewidth', 2);
set(0, 'defaultlinelinewidth', 2);
set(0, 'defaultpatchlinewidth', 0.7);
set(0, 'defaultaxesfontsize', 14);

load('Alphas.dat');
load('amps5.dat');
load('energy2.dat');
load('energy3.dat');
load('energy5.dat');

figure(1);
subplot(2,1,1);
semilogx(Alphas,amps5(:,2),'b',Alphas,amps5(:,3),'r',Alphas,amps5(:,1),'g'); shg

hold on;
yvals = linspace(0,1);
plot(4*ones(size(yvals)),yvals,'k--','LineWidth',1)
plot(8*ones(size(yvals)),yvals,'k--','LineWidth',1)
hold off;

obj = legend('$j=\pm 1$ modes','$j = 0$ mode','$j = \pm 2$ modes');
set(obj,'interpreter','latex','Location','SouthWest','FontSize',14);

xlim([1e-1,1e+2]);
ylim([0,1]);

xticks([1e-1,1,4,8, 100]);
set(gca,'Layer','top'); 
set(gca,'LineWidth',1)
box off;

xlabel('$\alpha$','interpreter','latex');
ylabel('$\max\limits_{x} |Z_j|$','interpreter','latex');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

subplot(2,1,2);
cols = parula(4);
semilogx(Alphas,energy5,'color',cols(1,:));
hold on;
semilogx(Alphas,energy3,'color',cols(2,:));
semilogx(Alphas,energy2,'color',cols(3,:));
hold off;

Ylim = 40;

xlim([1e-1,1e+2]);
ylim([0,Ylim]);

hold on;
yvals = Ylim*linspace(0,1);
plot(4*ones(size(yvals)),yvals,'k--','LineWidth',1)
plot(8*ones(size(yvals)),yvals,'k--','LineWidth',1)
hold off;

hold on;
plot(Alphas,8*ones(size(Alphas)),'k--','LineWidth',1);
hold off;

obj = legend('5-mode','3-mode','2-mode');
set(obj,'interpreter','latex','Location','NorthEast','FontSize',14);

xticks([1e-1,1,4,8, 100]);
set(gca,'Layer','top'); 
set(gca,'LineWidth',1)
box off;

xlabel('$\alpha$','interpreter','latex');
ylabel('energy','interpreter','latex');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% reran code used above with higher resolution for alpha near 8 for inset
load('Alphas2.dat');
load('phiamp5_2.dat');

figure(1);
subplot(2,1,1);
a = axes('position', [0.72 0.76 0.18 0.15]);
loglog(8-Alphas2,sin(phiamp5_2),'r.','MarkerSize',8);
hold on;
loglog(8-Alphas2,sqrt(8-Alphas2)*sqrt(pi/2)*0.248,'k','LineWidth',1.5);
hold off
xlim([1e-2,8]);

xlh = xlabel('$8-\alpha$','interpreter','latex');
xlh.Position(2) = xlh.Position(2) +0.005;
ylh = ylabel('$\max\limits_x |Z_0|$','interpreter','latex');
ylh.Position(1) = ylh.Position(1) +0.0008;
ylh.Position(2) = ylh.Position(2) -0.06;

a.XRuler.TickLabelGapOffset = -2; 
set(gca,'FontSize',10)

set(gca,'LineWidth',1)

%% Figure S1

figure

n0 = 5:20;
cols = parula(length(n0)+4);


for i = 1:length(n0)

    table = load(strcat('supp/',num2str(i),'.mat'));
    table = table.table;
    
    plot(table.X,table.Y1,'--','color',cols(i+1,:));
    hold on;
    plot(table.X,table.Y2,'--','color',cols(i+2,:));
    plot(table.X,table.Y3,'--','color',cols(i+3,:));
    plot(table.X,table.Y4,'--','color',cols(i+4,:));


    plot(table.Xb,table.Yb1,'--','color',cols(i+3,:));
    plot(table.Xb,table.Yb2,'--','color',cols(i+2,:));
    plot(table.Xb,table.Yb3,'--','color',cols(i+1,:));
    plot(table.Xb,table.Yb4,'--','color',cols(i,:));  

    plot(table.XSE, table.YSE,'color',cols(i+2,:));    

end

hold off;
xlim([40,75]);
xlabel('$x$','interpreter','latex');
ylabel('$k_i A_i$','interpreter','latex');


%% Figure S2
s2 = load('SuppFig2.mat');

table1 = s2.s2.table1;
table2 = s2.s2.table2;

figure

subplot(3,2,1);
plot(table1.x,table1.phi,'c');
xlim([-5,5])
ylim([0 pi/2])
xlabel('$\bar{X}$','interpreter','latex')
ylabel('$\phi$','interpreter','latex')

subplot(3,2,3);
plot(table1.x,table1.theta,'c');
xlim([-5,5]);
ylim([0,pi/2]);
xlabel('$\bar{X}$','interpreter','latex')
ylabel('$\theta$','interpreter','latex')

subplot(3,2,5)
plot(table1.x,table1.X1,'b',table1.x,table1.X3,'r',table1.x,table1.X5,'g')
xlim([-5,5])
xlabel('$\bar{X}$','interpreter','latex')
ylabel('$Z_{-2},Z_{0},Z_{+2}$','interpreter','latex')

subplot(3,2,2);
plot(table2.x,table2.phi,'m');
xlim([-5,5])
ylim([0 pi/2])
xlabel('$\bar{X}$','interpreter','latex')
ylabel('$\phi$','interpreter','latex')

subplot(3,2,4);
plot(table2.x,table2.theta,'m');
xlim([-5,5]);
ylim([0,pi/2]);
xlabel('$\bar{X}$','interpreter','latex')
ylabel('$\theta$','interpreter','latex')

subplot(3,2,6)
plot(table2.x,table2.X1,'b',table2.x,table2.X3,'r',table2.x,table2.X5,'g')
xlim([-5,5])
xlabel('$\bar{X}$','interpreter','latex')
ylabel('$Z_{-2},Z_{0},Z_{+2}$','interpreter','latex')

%% Figure S3

a = load('SI_Figure3.dat');

figure
plot(a(:,1),a(:,2),'b'); 
xlabel('$X$','interpreter','latex');
ylabel('$\Theta$','interpreter','latex');
xlim([min(a(:,1)),max(a(:,1))])














