0;

widthCM = 13.5;
fig = figure('units','centimeters','position',...
    [5 -5 widthCM widthCM]);

dim = 2;
beta0 = 1.2;

if dim == 2
    lambda = 2pi3;
    omega  = 3;
    machvec = 00.010.1;
    lengthvec = (2.50.510)'lambda;
    n_contours = 15;
else
    lambda = 2pi0.1;
    omega  = 0.1;
    machvec = 00.00250.025;
    lengthvec = (0.50.52)'lambda;
    n_contours = 20;
end
machmach = repmat(machvec,length(lengthvec),1);
lengthlength = repmat(lengthvec,1,length(machvec));
s_hats_hat = zeros(length(lengthvec),length(machvec));
    

for I_L = 1length(lengthvec)
    L = lengthvec(I_L)lambda;
    if L  10
        lengthstring = strcat('0',sprintf('%02.1f',L));
    else
        lengthstring = sprintf('%02.1f',L);
    end
    if dim == 2
        filename = strcat('s_hat_machdata2D_',lengthstring,'_omega3.0');
        load(filename,'s_hat2vec','geomparams')
        X_i = geomparams.X_i;
        s_hats_hat(I_L,) = (abs(s_hat2vec'X_i));
    else
        filename = strcat('s_hat_machdata3D_',lengthstring,'_omega0.1');
        load(filename,'s_hat2vec','geomparams')
        R_i = geomparams.R_i;
        s_hats_hat(I_L,) = (abs(s_hat2vec'R_i));
    end
end

upperclim = max(s_hats_hat,[],'all');
lowerclim = min(s_hats_hat,[],'all');
midpoint = mean(s_hats_hat(,1));
maxdeviation = max([upperclim-midpoint midpoint-lowerclim]);


meandrift = mean(s_hats_hat(,end) - s_hats_hat(,1));

if dim == 2
    contourlevels = (midpoint-maxdeviation)...
        maxdeviationn_contoursmidpoint+maxdeviation;
elseif dim == 3
    lowerclim = min(s_hats_hat,[],'all');
    contourlevels = lowerclim2midpointn_contours2midpoint;
end



set_custom_colormap(1)
[~,c] = contourf(machmach,lengthlength,s_hats_hat);
c.LevelList = contourlevels;

% p = pcolor(machmach,lengthlength,s_hats_hat);
% p.EdgeColor = 'none';

clim([contourlevels(1),contourlevels(length(contourlevels))]);

xlabel('$$M$$')
ylabel('$$Llambda$$')
cbar = colorbar;
cbar.Label.Interpreter = 'latex';
s_sf_min = 1(2pibeta0machvec(end));
if gca().YLim(2)  s_sf_min
    hold on
    minmach = (2pigca().YLim(2)beta0)^-1;
    maxmach = machvec(end);
    newmachvec = minmach(maxmach-minmach)10maxmach;
    plot(newmachvec,1.(2pibeta0newmachvec),'g','LineWidth',3)
end


if dim == 2
    cbar.Label.String = '$$tilde{s}_0^1X_1$$';
elseif dim == 3
    cbar.Label.String = '$$tilde{s}_0^1R_1$$';
end
exportgraphics(gcf,'s_hat_machplot.pdf','ContentType','vector')

