% InversionSpectrumFigure_v1

% Mike Wolovick, 1/17/2023

% This script makes a figure comparing spectral analysis of the inversion
% results for several restricted domains.



clearvars
tic
close all

%% Parameters:

% File names:
% lcurvefile='/albedo/work/user/mwolovic/FRISModeling/LCurves/FRISDomainLCurve_ResTest_SSA_drag_v2_mesh1a.mat';  % include trailing ".mat"
% localpath='/home/mvolovic/Dropbox/FRISModeling/LCurves/';
% %issmscriptpath='/albedo/home/mwolovic/ISSM_installations/trunk-jpl/';
% issmscriptpath='~/Dropbox/ISSM/';
% myscriptpath='~/Dropbox/MATLAB/';
% %figname='/albedo/home/mwolovic/InversionComparisonFigure_v4.png';
% figname='~/Dropbox/FRISModeling/InversionSpectrumFigure_v1a.png';
lcurvefile='/work/ollie/mwolovic/FRISModeling/LCurves/FRISDomainLCurve_ResTest_SSA_drag_v2_mesh1a.mat';  % include trailing ".mat"
localpath='/home/mvolovic/Dropbox/FRISModeling/LCurves/';
%issmscriptpath='/albedo/home/mwolovic/ISSM_installations/trunk-jpl/';
issmscriptpath='/home/ollie/mwolovic/ISSM_installations/trunk-jpl/';
myscriptpath='/home/ollie/mwolovic/MATLAB/';
%figname='/albedo/home/mwolovic/InversionComparisonFigure_v4.png';
figname='/home/ollie/mwolovic/InversionSpectrumFigure_v3.png';

% Server control:
%serverloc='albedo0.dmawi.de';  % can be commented
serverusername='mwolovic';
serverpassword='';  % sorry for saving this in plain text.  
% Maybe if your graphical connection wasn't so damn slow I wouldn't need to run matlab from my home PC, AWI.

% Region to analyze:
zoomboxes=[-600,-600,-466,-466,274,348,348,274;...
    -504,-539,-384,-349,934,1003,1082,1013;...
    -490,-533,-346,-295,753,849,954,872;...
    -1410,-1415,-1282,-1277,90,143,160,106]*1e3;  % [Nx8] m, [x1,x2,x3,x4,y1,y2,y3,y4], points should be listed clockwise from lower left
zoomboxnames={'Academy Glacier';'Slessor Glacier';'Recovery Glacier';'Rutford Ice Stream'};  % Nx1 cell array of strings

% Analysis parameters:
dx_grid=100;                        % m
normbufferdist=.5;                  % unitless, fraction of smaller box dimension
% (smooth buffer region around the margin of each grid box to minimize fft2 artifacts caused by the assumption of periodic BC)

% Ticks and limits:
klims=[.01,1];                            % 1/km

% Other aesthetic settings:
titlefontsize=14;                              % points
axesfontsize=12;                              % points
drivecolor='k';                                 % valid color string or RGB
maxcolor='b';                                  % valid color string or RGB
bestcolor='m';                                  % valid color string or RGB
mincolor='r';                                  % valid color string or RGB

% Subplot controls:
vertaxesbuffer=.08;             % unitless
titlebuffer=.07;                    % unitless
horzaxesbuffer=.06;             % unitless
horznotextbuffer=.02;          % unitless

% Output size and resolution:
pagesize=[20,10];                 % inches (height computed automatically)
resolution=300;                     % dpi

%% Preparation:

% Compute number of zoom boxes:
numzoomboxes=size(zoomboxes,1);

% Compute subplot dimensions:
% Hardcoded subplot arrangement:
numvertplots=2;
numhorzplots=numzoomboxes;
% Compute available space:
availablehorzspace=1-numhorzplots*horznotextbuffer-horzaxesbuffer;
availablevertspace=1-numvertplots*titlebuffer-vertaxesbuffer;
% Compute boxes:
Boxes=cell(numvertplots,numhorzplots);
for ii=1:numvertplots
    for jj=1:numhorzplots
        Boxes{numvertplots-ii+1,jj}=[horzaxesbuffer+(jj-1)*(horznotextbuffer+availablehorzspace/numhorzplots),...
            vertaxesbuffer+(ii-1)*(titlebuffer+availablevertspace/numvertplots),...
            availablehorzspace/numhorzplots,availablevertspace/numvertplots];
    end
end

% Add the ISSM script path:
addpath(genpath(issmscriptpath))

% Add my script path:
addpath(genpath(myscriptpath))

% Check if we're loading remote files:
if exist('serverloc','var')
    % Parse filename:
    slashnum=strfind(lcurvefile,'/');
    remotepath=lcurvefile(1:slashnum(end));
    lcurvename=lcurvefile(slashnum(end)+1:end);
    % Fetch the file:
    if ~exist([localpath,lcurvename],'file')
        ssh2_struct=scp_simple_get(serverloc,serverusername,serverpassword,lcurvename,localpath,remotepath);
    end
    % Load the file:
    load([localpath,lcurvename],'BasicModel','MyInversionParameters','*best*','Lambda')
else
    % Load basic model:
    load(lcurvefile,'BasicModel','MyInversionParameters','*best*','Lambda')
end

% Rename model structure:
md=BasicModel;

% Unpack inversion parameters:
unpack(MyInversionParameters)

% Identify max/min acceptable models:
bestmodelind_max=find(Lambda<=bestlambda_maxallowable,1,'last');
bestmodelind_min=find(Lambda>=bestlambda_minallowable,1,'first');

% Compute mesh properties:
% Extract 2D mesh:
try
    x2d=md.mesh.x2d;
    y2d=md.mesh.y2d;
    elements2d=md.mesh.elements2d;
    oldnumlayers=md.mesh.numberoflayers;
catch
    x2d=md.mesh.x;
    y2d=md.mesh.y;
    elements2d=md.mesh.elements;
    oldnumlayers=1;
end
numelements2d=size(elements2d,1);
numvertices2d=length(x2d);
% Compute edge lengths:
a=sqrt((x2d(elements2d(:,1))-x2d(elements2d(:,2))).^2+(y2d(elements2d(:,1))-y2d(elements2d(:,2))).^2);
b=sqrt((x2d(elements2d(:,2))-x2d(elements2d(:,3))).^2+(y2d(elements2d(:,2))-y2d(elements2d(:,3))).^2);
c=sqrt((x2d(elements2d(:,3))-x2d(elements2d(:,1))).^2+(y2d(elements2d(:,3))-y2d(elements2d(:,1))).^2);
% Compute area:
s=(a+b+c)/2;
Area=sqrt(s.*(s-a).*(s-b).*(s-c));
% Compute hydraulic diameter:
HydraulicDiameter_elements=2*Area./s;
% Compute element centers:
CentX=mean(x2d(elements2d),2);
CentY=mean(y2d(elements2d),2);
% Find mesh limits:
xlims=[min(x2d),max(x2d)];
ylims=[min(y2d),max(y2d)];
% Extract masks on 2D mesh and truncate:
OceanLevelSet_vertices=max(-1,min(1,md.mask.ocean_levelset(1:numvertices2d)));
IceLevelSet_vertices=max(-1,min(1,md.mask.ice_levelset(1:numvertices2d)));
% Compute ocean and ice area fractions:
OceanAreaFraction_elements=(1-mean(OceanLevelSet_vertices(elements2d),2))/2;
IceAreaFraction_elements=(1-mean(IceLevelSet_vertices(elements2d),2))/2;
% Compute area totals:
totalarea=sum(Area);
totalicearea=sum(Area.*IceAreaFraction_elements);
totalgroundedarea=sum(Area.*IceAreaFraction_elements.*(1-OceanAreaFraction_elements));
% Compute gradient of linear basis functions:
BasisGradX_elements=[(y2d(elements2d(:,2))-y2d(elements2d(:,3)))./(x2d(elements2d(:,1)).*y2d(elements2d(:,2))-x2d(elements2d(:,2)).*y2d(elements2d(:,1))-x2d(elements2d(:,1)).*y2d(elements2d(:,3))...
    +x2d(elements2d(:,3)).*y2d(elements2d(:,1))+x2d(elements2d(:,2)).*y2d(elements2d(:,3))-x2d(elements2d(:,3)).*y2d(elements2d(:,2))),...
    (y2d(elements2d(:,3))-y2d(elements2d(:,1)))./(x2d(elements2d(:,2)).*y2d(elements2d(:,3))-x2d(elements2d(:,3)).*y2d(elements2d(:,2))-x2d(elements2d(:,2)).*y2d(elements2d(:,1))...
    +x2d(elements2d(:,1)).*y2d(elements2d(:,2))+x2d(elements2d(:,3)).*y2d(elements2d(:,1))-x2d(elements2d(:,1)).*y2d(elements2d(:,3))),...
    (y2d(elements2d(:,1))-y2d(elements2d(:,2)))./(x2d(elements2d(:,3)).*y2d(elements2d(:,1))-x2d(elements2d(:,1)).*y2d(elements2d(:,3))-x2d(elements2d(:,3)).*y2d(elements2d(:,2))...
    +x2d(elements2d(:,2)).*y2d(elements2d(:,3))+x2d(elements2d(:,1)).*y2d(elements2d(:,2))-x2d(elements2d(:,2)).*y2d(elements2d(:,1)))];
BasisGradY_elements=[-(x2d(elements2d(:,2))-x2d(elements2d(:,3)))./(x2d(elements2d(:,1)).*y2d(elements2d(:,2))-x2d(elements2d(:,2)).*y2d(elements2d(:,1))-x2d(elements2d(:,1)).*y2d(elements2d(:,3))...
    +x2d(elements2d(:,3)).*y2d(elements2d(:,1))+x2d(elements2d(:,2)).*y2d(elements2d(:,3))-x2d(elements2d(:,3)).*y2d(elements2d(:,2))),...
    -(x2d(elements2d(:,3))-x2d(elements2d(:,1)))./(x2d(elements2d(:,2)).*y2d(elements2d(:,3))-x2d(elements2d(:,3)).*y2d(elements2d(:,2))-x2d(elements2d(:,2)).*y2d(elements2d(:,1))...
    +x2d(elements2d(:,1)).*y2d(elements2d(:,2))+x2d(elements2d(:,3)).*y2d(elements2d(:,1))-x2d(elements2d(:,1)).*y2d(elements2d(:,3))),...
    -(x2d(elements2d(:,1))-x2d(elements2d(:,2)))./(x2d(elements2d(:,3)).*y2d(elements2d(:,1))-x2d(elements2d(:,1)).*y2d(elements2d(:,3))-x2d(elements2d(:,3)).*y2d(elements2d(:,2))...
    +x2d(elements2d(:,2)).*y2d(elements2d(:,3))+x2d(elements2d(:,1)).*y2d(elements2d(:,2))-x2d(elements2d(:,2)).*y2d(elements2d(:,1)))];

% Compute driving stress:
% Compute surface gradients on the elements:
SurfGradX_elements=sum(BasisGradX_elements.*md.geometry.surface(elements2d),2);
SurfGradY_elements=sum(BasisGradY_elements.*md.geometry.surface(elements2d),2);
% Average surface gradient from elements onto vertices:
% pre-allocate:
SurfGradX_vertices=zeros(numvertices2d,1);
SurfGradY_vertices=zeros(numvertices2d,1);
Denom_vertices=zeros(numvertices2d,1);
% loop through elements:
for jj=1:numelements2d
    % Add to the running sums:
    SurfGradX_vertices(elements2d(jj,:))=SurfGradX_vertices(elements2d(jj,:))+SurfGradX_elements(jj);
    SurfGradY_vertices(elements2d(jj,:))=SurfGradY_vertices(elements2d(jj,:))+SurfGradY_elements(jj);
    Denom_vertices(elements2d(jj,:))=Denom_vertices(elements2d(jj,:))+1;
end
% finish computing means:
SurfGradX_vertices=SurfGradX_vertices./Denom_vertices;
SurfGradY_vertices=SurfGradY_vertices./Denom_vertices;
% Compute driving stress:
DrivingStressX_vertices=-md.materials.rho_ice*md.constants.g*SurfGradX_vertices.*md.geometry.thickness(1:numvertices2d);
DrivingStressY_vertices=-md.materials.rho_ice*md.constants.g*SurfGradY_vertices.*md.geometry.thickness(1:numvertices2d);
DrivingStressMag_vertices=sqrt(DrivingStressX_vertices.^2+DrivingStressY_vertices.^2);
% Resolve driving stress in the flow direction:
ObsVx=md.inversion.vx_obs(end-numvertices2d+1:end);
ObsVy=md.inversion.vy_obs(end-numvertices2d+1:end);
ObsVel=sqrt(ObsVx.^2+ObsVy.^2);
DrivingStressFlowDir_vertices=(DrivingStressX_vertices.*ObsVx+DrivingStressY_vertices.*ObsVy)./ObsVel;
% Set to zero under floating shelf and in ice-free areas:
DrivingStressFlowDir_vertices(OceanLevelSet_vertices<0|IceLevelSet_vertices>0)=0;

% Generate the active mask:
ActiveMask_drag=IceLevelSet_vertices<=0&OceanLevelSet_vertices>=0;
ActiveMask_elements=mean(ActiveMask_drag(elements2d),2);

% Unpack the good data mask:
GoodData_obs_vertices=md.miscellaneous.dummy(1:numvertices2d);

% Ensure exact rectangles:
% Compute side lengths:
Dist=zeros(numzoomboxes,4);
for ii=1:numzoomboxes
    Dist(ii,:)=sqrt([diff(zoomboxes(ii,1:4)).^2+diff(zoomboxes(ii,5:8)).^2,(zoomboxes(ii,1)-zoomboxes(ii,4))^2+(zoomboxes(ii,5)-zoomboxes(ii,8))^2]);
end
% Compute mean width and length:
MeanWidth=.5*(Dist(:,1)+Dist(:,3));
MeanLength=.5*(Dist(:,2)+Dist(:,4));
% Compute heading vectors: (real=x, imag=y)
% set up to point in consistent directions:
% 1->2, 2->3, 4->3, 1->4
HeadingVector=[(zoomboxes(:,2)-zoomboxes(:,1)+1i*(zoomboxes(:,6)-zoomboxes(:,5)))./Dist(:,1),(zoomboxes(:,3)-zoomboxes(:,2)+1i*(zoomboxes(:,7)-zoomboxes(:,6)))./Dist(:,2),...
    (zoomboxes(:,3)-zoomboxes(:,4)+1i*(zoomboxes(:,7)-zoomboxes(:,8)))./Dist(:,3),(zoomboxes(:,4)-zoomboxes(:,1)+1i*(zoomboxes(:,8)-zoomboxes(:,5)))./Dist(:,4)];
% Compute mean (distance-weighted) heading vector:
% mean vector points in same direction as 2->3 and 1->4
MeanHeadingVector=sum([imag(HeadingVector(:,1))-1i*real(HeadingVector(:,1)),HeadingVector(:,2),imag(HeadingVector(:,3))-1i*real(HeadingVector(:,3)),HeadingVector(:,4)].*Dist,2)./sum(Dist,2);
% Construct new zoom box based on mean dimensions and heading, fixed
% center:
ZoomBoxCentX=mean(zoomboxes(:,1:4),2);
ZoomBoxCentY=mean(zoomboxes(:,5:8),2);
zoomboxes=[repmat(ZoomBoxCentX,[1,4])+.5*[-MeanLength.*real(MeanHeadingVector)+MeanWidth.*imag(MeanHeadingVector),-MeanLength.*real(MeanHeadingVector)-MeanWidth.*imag(MeanHeadingVector),...
    MeanLength.*real(MeanHeadingVector)-MeanWidth.*imag(MeanHeadingVector),MeanLength.*real(MeanHeadingVector)+MeanWidth.*imag(MeanHeadingVector)],...
    repmat(ZoomBoxCentY,[1,4])+.5*[-MeanLength.*imag(MeanHeadingVector)-MeanWidth.*real(MeanHeadingVector),-MeanLength.*imag(MeanHeadingVector)+MeanWidth.*real(MeanHeadingVector),...
    MeanLength.*imag(MeanHeadingVector)+MeanWidth.*real(MeanHeadingVector),MeanLength.*imag(MeanHeadingVector)-MeanWidth.*real(MeanHeadingVector)]];

% Compute buffer distances:
BufferDist=normbufferdist*min(MeanLength,MeanWidth);

% Compute rotated coordinates relative to corners:
x2d_rot=repmat(real(MeanHeadingVector'),[numvertices2d,1]).*(repmat(x2d,[1,numzoomboxes])-repmat(zoomboxes(:,1)',[numvertices2d,1]))+...
    repmat(-imag(MeanHeadingVector'),[numvertices2d,1]).*(repmat(y2d,[1,numzoomboxes])-repmat(zoomboxes(:,5)',[numvertices2d,1]));
y2d_rot=repmat(imag(MeanHeadingVector'),[numvertices2d,1]).*(repmat(x2d,[1,numzoomboxes])-repmat(zoomboxes(:,1)',[numvertices2d,1]))+...
    repmat(real(MeanHeadingVector'),[numvertices2d,1]).*(repmat(y2d,[1,numzoomboxes])-repmat(zoomboxes(:,5)',[numvertices2d,1]));

% Identify vertices within the zoom boxes:
WithinZoomBox_vertices=x2d_rot>=0&x2d_rot<=repmat(MeanLength',[numvertices2d,1])&y2d_rot>=0&y2d_rot<=repmat(MeanWidth',[numvertices2d,1]);

% Identify elements within the zoom boxes:
WithinZoomBoxes_elements=zeros(numelements2d,numzoomboxes);
for ii=1:numzoomboxes
    WithinZoomBoxes_elements(:,ii)=mean(reshape(WithinZoomBox_vertices(elements2d,ii),[],3),2);
end

% Compute mean ice thicknesses:
Icethick_elements=mean(md.geometry.thickness(elements2d),2);
boxesmeanicethick=zeros(numzoomboxes,1);
for ii=1:numzoomboxes
    boxesmeanicethick(ii)=sum(Area.*Icethick_elements.*WithinZoomBoxes_elements(:,ii))/sum(Area.*WithinZoomBoxes_elements(:,ii));
end

% Load all model drags:
% pre-allocate:
Drag=zeros(numvertices2d,4);
MeanDrag=zeros(4,numzoomboxes);
% first one is driving stress:
Drag(:,1)=DrivingStressFlowDir_vertices;
% Smoothly fill driving stress in ice-free areas:
Drag(IceLevelSet_vertices>0,1)=NaN;
Drag(:,1)=LaplacianInterpolation_mesh(x2d,y2d,elements2d,Drag(:,1));
% Compute mean within each box:
ThisDrag_elements=mean(reshape(Drag(elements2d,1),[numelements2d,3]),2);
for ii=1:numzoomboxes
    MeanDrag(1,ii)=sum(Area.*ThisDrag_elements.*WithinZoomBoxes_elements(:,ii))/sum(Area.*WithinZoomBoxes_elements(:,ii));
end
% Loop through models:
for jj=1:3
    % Get this model ind:
    if jj==1
        thismodelind=bestmodelind_min;
    elseif jj==2
        thismodelind=bestmodelind;
    else
        thismodelind=bestmodelind_max;
    end
    % Construct this file name:
    slashnum=strfind(lcurvefile,'/');
    thisbase=10^floor(log10(Lambda(thismodelind)));
    thisname=[lcurvefile(slashnum(end)+1:end-4),'_Lambda_',num2str(thisbase*.01*round(100*Lambda(thismodelind)/thisbase)),'.mat'];
    % Check if we're loading remote files:
    if exist('serverloc','var')
        % Fetch the file:
        if ~exist([localpath,thisname],'file')
            ssh2_struct=scp_simple_get(serverloc,serverusername,serverpassword,thisname,localpath,remotepath);
        end
        % Load the file:
        load([localpath,thisname],'md')
    else
        % Load this model:
        load([lcurvefile(1:slashnum(end)),thisname],'md')
    end
    % Extract variables we need:
     if numlayers>1 && ~strcmp(flowequation,'SSA')
        ModelBasalVel=md.results.StressbalanceSolution.Vel(md.mesh.vertexonbase==1);
        ModelDragCoeff=md.results.StressbalanceSolution.FrictionCoefficient(md.mesh.vertexonbase==1);
        if ~strcmp(nsource,'none')
            EffectivePressure_vertices=md.friction.effective_pressure(md.mesh.vertexonbase==1);
        end
     else
        ModelBasalVel=md.results.StressbalanceSolution.Vel;
        ModelDragCoeff=md.results.StressbalanceSolution.FrictionCoefficient;
        if ~strcmp(nsource,'none')
            EffectivePressure_vertices=md.friction.effective_pressure;
        end
     end
     % Compute drag:
     Drag(:,jj+1)=(ModelDragCoeff.^2).*((ModelBasalVel/secondsperyear).^(1/m));
     if ~strcmp(nsource,'none')
         Drag(:,jj+1)=Drag(:,jj+1).*(EffectivePressure_vertices.^r);
     end
     % Set drag to zero under floating shelf:
     Drag(OceanLevelSet_vertices<0&IceLevelSet_vertices<0,jj+1)=0;
     % Smoothly fill drag in ice-free areas:
     Drag(IceLevelSet_vertices>0,jj+1)=NaN;
     Drag(:,jj+1)=LaplacianInterpolation_mesh(x2d,y2d,elements2d,Drag(:,jj+1));
     % Compute mean within each box:
     ThisDrag_elements=mean(reshape(Drag(elements2d,jj+1),[numelements2d,3]),2);
     for ii=1:numzoomboxes
         MeanDrag(jj+1,ii)=sum(Area.*ThisDrag_elements.*WithinZoomBoxes_elements(:,ii))/sum(Area.*WithinZoomBoxes_elements(:,ii));
     end
end

%% Spectral analysis:

% Generate spatial grids:
% pre-allocate:
SpatialGrids=struct('X_grid',[],'Y_grid',[],'xsize',[],'ysize',[],'Drag_grid',[]);
% Loop through zoom boxes:
for ii=1:numzoomboxes
    % generate first grid:
    %     SpatialGrids(ii).X_grid=dx_grid*[0:1:ceil(MeanLength(ii)/dx_grid)];
    %     SpatialGrids(ii).Y_grid=dx_grid*[0:1:ceil(MeanWidth(ii)/dx_grid)]';
    SpatialGrids(ii).X_grid=dx_grid*[floor(-BufferDist(ii)/dx_grid):1:ceil((MeanLength(ii)+BufferDist(ii))/dx_grid)];
    SpatialGrids(ii).Y_grid=dx_grid*[floor(-BufferDist(ii)/dx_grid):1:ceil((MeanWidth(ii)+BufferDist(ii))/dx_grid)]';
    % get grid size:
    SpatialGrids(ii).ysize=length(SpatialGrids(ii).Y_grid);
    SpatialGrids(ii).xsize=length(SpatialGrids(ii).X_grid);
    % Truncate to even number of grid samples:
    if rem(SpatialGrids(ii).xsize,2)==1
        SpatialGrids(ii).X_grid=SpatialGrids(ii).X_grid(1:end-1);
        SpatialGrids(ii).xsize=SpatialGrids(ii).xsize-1;
    end
    if rem(SpatialGrids(ii).ysize,2)==1
        SpatialGrids(ii).Y_grid=SpatialGrids(ii).Y_grid(1:end-1);
        SpatialGrids(ii).ysize=SpatialGrids(ii).ysize-1;
    end
end

% Generate wavenumber grids:
% pre-allocate:
SpectralGrids=struct('Kx',[],'Ky',[],'dkx',[],'dky',[],'Kx_mesh',[],'Ky_mesh',[],'Kmag',[],'DragSpectra',[]);
% compute sampling frequency:
ks=1/dx_grid;
% Loop through zoom boxes:
for ii=1:numzoomboxes
    % Generate positive wavenumber grid:
    SpectralGrids(ii).Kx=ks*(0:(SpatialGrids(ii).xsize/2))/SpatialGrids(ii).xsize;
    SpectralGrids(ii).Ky=ks*(0:(SpatialGrids(ii).ysize/2))/SpatialGrids(ii).ysize;
    SpectralGrids(ii).dkx=SpectralGrids(ii).Kx(2)-SpectralGrids(ii).Kx(1);
    SpectralGrids(ii).dky=SpectralGrids(ii).Ky(2)-SpectralGrids(ii).Ky(1);
    % Generate full wavenumber grid:
    if rem(SpatialGrids(ii).xsize,2)==1
        SpectralGrids(ii).Kx=[SpectralGrids(ii).Kx,fliplr(-SpectralGrids(ii).Kx(2:end))];
    else
        SpectralGrids(ii).Kx=[SpectralGrids(ii).Kx,fliplr(-SpectralGrids(ii).Kx(2:end-1))];
    end
    if rem(SpatialGrids(ii).ysize,2)==1
        SpectralGrids(ii).Ky=[SpectralGrids(ii).Ky,fliplr(-SpectralGrids(ii).Ky(2:end))]';
    else
        SpectralGrids(ii).Ky=[SpectralGrids(ii).Ky,fliplr(-SpectralGrids(ii).Ky(2:end-1))]';
    end
    % Make meshgrid:
    [SpectralGrids(ii).Kx_mesh,SpectralGrids(ii).Ky_mesh]=meshgrid(SpectralGrids(ii).Kx,SpectralGrids(ii).Ky);
    % Compute spectral magnitude:
    SpectralGrids(ii).Kmag=sqrt(SpectralGrids(ii).Kx_mesh.^2+SpectralGrids(ii).Ky_mesh.^2);
end

% Generate 1D spectral sampling:
Kmag_1d=ks*(0:(median([SpatialGrids(:).xsize,SpatialGrids(:).ysize])/2))/median([SpatialGrids(:).xsize,SpatialGrids(:).ysize]);
Kmag_1d=Kmag_1d(Kmag_1d<=1.05*klims(2)/1000);
num1dk=length(Kmag_1d);
dk=Kmag_1d(2)-Kmag_1d(1);

% Interpolate maps to spatial grids:
% Loop through boxes:
for ii=1:numzoomboxes
    % Pre-allocate cell:
    SpatialGrids(ii).Drag_grid=cell(4,1);
    % Loop through drag maps:
    for jj=1:4
        % Basic interpolation:
        SpatialGrids(ii).Drag_grid{jj}=InterpFromMeshToGrid(elements2d,x2d_rot(:,ii),y2d_rot(:,ii),Drag(:,jj),SpatialGrids(ii).X_grid,SpatialGrids(ii).Y_grid,0);
        % Check if we have a buffer:
        if normbufferdist>0
            % Set areas outside of the analysis box to NaN:
            SpatialGrids(ii).Drag_grid{jj}(repmat(SpatialGrids(ii).X_grid<0|SpatialGrids(ii).X_grid>MeanLength(ii),[SpatialGrids(ii).ysize,1])|...
                repmat(SpatialGrids(ii).Y_grid<0|SpatialGrids(ii).Y_grid>MeanWidth(ii),[1,SpatialGrids(ii).xsize]))=NaN;
            % Set margins to mean value:
            SpatialGrids(ii).Drag_grid{jj}(1,:)=MeanDrag(jj,ii);
            SpatialGrids(ii).Drag_grid{jj}(end,:)=MeanDrag(jj,ii);
            SpatialGrids(ii).Drag_grid{jj}(:,1)=MeanDrag(jj,ii);
            SpatialGrids(ii).Drag_grid{jj}(:,end)=MeanDrag(jj,ii);
            % Smoothly interpolate within buffer zone:
            SpatialGrids(ii).Drag_grid{jj}=LaplacianInterpolation(SpatialGrids(ii).Drag_grid{jj},9);
        end
    end
end

% Take 2D spectra:
for ii=1:numzoomboxes
    SpectralGrids(ii).DragSpectra=cell(4,1);
    for jj=1:4
        SpectralGrids(ii).DragSpectra{jj}=fft2(SpatialGrids(ii).Drag_grid{jj}/(ks*SpatialGrids(ii).xsize*ks*SpatialGrids(ii).ysize)); % denominator normalizes to 2D power spectral density; units= Pa^2/(m^-2)
    end
end

% Convert from 2D spectra to 1D spectral sampling:
% Pre-allocate:
DragPowerSpectra_1d=zeros(num1dk,4,numzoomboxes);
DragCoherence_1d=zeros(num1dk,3,numzoomboxes);
% Loop through zoom boxes:
for ii=1:numzoomboxes
    % Loop through 1D samples:
    for kk=1:num1dk

        % Identify nearby spectral points:
        TheseK=abs(SpectralGrids(ii).Kmag-Kmag_1d(kk))<3*dk;

        % Compute weighting function:
        TheseWeights=SpectralGrids(ii).dkx*SpectralGrids(ii).dky*exp(-.5*(((SpectralGrids(ii).Kmag(TheseK)-Kmag_1d(kk))/dk).^2))/(dk*sqrt(2*pi));

        % Explanation of weighting function:
        % dkx*dky integrates power in kx,ky space;
        % exp() restricts integral to thin region around a circular arc of constant kmag;
        % denominator normalizes to unit integral, evaluated along a 1D radial line in kx,ky space.

        % Perform weighted sum for 1D power spectral density:
        for jj=1:4
            DragPowerSpectra_1d(kk,jj,ii)=sum(TheseWeights.*(abs(SpectralGrids(ii).DragSpectra{jj}(TheseK)).^2));
        end

        % Perform weighted average for coherence:
        for jj=1:3
            DragCoherence_1d(kk,jj,ii)=sum(TheseWeights.*(real(SpectralGrids(ii).DragSpectra{jj+1}(TheseK)).*real(SpectralGrids(ii).DragSpectra{1}(TheseK))+...
                imag(SpectralGrids(ii).DragSpectra{jj+1}(TheseK)).*imag(SpectralGrids(ii).DragSpectra{1}(TheseK))))/...
                sum(TheseWeights.*abs(SpectralGrids(ii).DragSpectra{jj+1}(TheseK)).*abs(SpectralGrids(ii).DragSpectra{1}(TheseK)));
        end

    end
end

% Note on units and meaning of 1D spectra:
% The 1D spectra integrate around a circle in kx,ky space.  The gaussian
% weighting function is only there to account for the fact that the 2D
% sampling of kx,ky space is on a rectangular grid, not along circles of
% constant kmag.  The guassian weighting function is therefore normalized
% to unit RADIAL (kmag) integral.  However, it should be noted that this is not
% the same as being normalized to unit integral in 2D kx,ky space.  This is
% important, because the product we want to end up with is total power per
% differential 1D wavenumber.  Higher wavenumbers integrate along a larger
% circle in kx,ky space, and that is okay.  They have more kx,ky
% contributing power and that should increase their total.  Thus, the
% weighting function first integrates all power in the 2D grid (sum of 2D spectra
% multiplied by dkx*dky), then multiplies by a radial gaussian that has
% been normalized to unit integral along the kmag radial line.

% Thus, the units of the 2D spectra in kx,ky space are:  Pa^2/(m^-2), but
% the units of the 1D spectra in radial kmag space are:  Pa^2/(m^-1).

%% Make Figure:

% Make the figure:
figure(1)

% Compute power ylims:
firstkmagind=find(Kmag_1d*1000>=klims(1),1,'first');
lastkmagind=find(Kmag_1d*1000<=klims(2),1,'last');
powerylims=(10.^[floor(log10(min(min(min(DragPowerSpectra_1d(firstkmagind:lastkmagind,:,:)))))),ceil(log10(max(max(max(DragPowerSpectra_1d(firstkmagind:lastkmagind,:,:))))))])*1e-9;
powerylims=powerylims.*[.999,1.001];

% Loop through zoom boxes:
for ii=1:numzoomboxes

    % Make power spectral plot:
    % Make subplot:
    subplot('Position',Boxes{1,ii})
    % Plot spectra:
    plot(Kmag_1d(1:end-1)*1000,squeeze(DragPowerSpectra_1d(1:end-1,1,ii))*1e-9,drivecolor)
    hold on
    plot(Kmag_1d(1:end-1)*1000,squeeze(DragPowerSpectra_1d(1:end-1,2,ii))*1e-9,mincolor)
    plot(Kmag_1d(1:end-1)*1000,squeeze(DragPowerSpectra_1d(1:end-1,3,ii))*1e-9,bestcolor)
    plot(Kmag_1d(1:end-1)*1000,squeeze(DragPowerSpectra_1d(1:end-1,4,ii))*1e-9,maxcolor)
    % Set to log/log scale:
    set(gca,'XScale','log')
    set(gca,'YScale','log')
    % Plot and label select wavelengths:
    theseylims=powerylims;
    plot((1000/boxesmeanicethick(ii))*[1,1],theseylims,'--k')
    text(1000/boxesmeanicethick(ii),theseylims(1),'H','HorizontalAlignment','left','VerticalAlignment','bottom','FontSize',axesfontsize)
    plot((1000/(2*boxesmeanicethick(ii)))*[1,1],theseylims,'--k')
    text(1000/(2*boxesmeanicethick(ii)),theseylims(1),'2H','HorizontalAlignment','left','VerticalAlignment','bottom','FontSize',axesfontsize)
    %plot((1000/(3*boxesmeanicethick(ii)))*[1,1],theseylims,'--k')
    %text(1000/(3*boxesmeanicethick(ii)),theseylims(1),'3H','HorizontalAlignment','left','VerticalAlignment','bottom')
    plot((1000/(5*boxesmeanicethick(ii)))*[1,1],theseylims,'--k')
    text(1000/(5*boxesmeanicethick(ii)),theseylims(1),'5H','HorizontalAlignment','left','VerticalAlignment','bottom','FontSize',axesfontsize)
    plot((1000/(10*boxesmeanicethick(ii)))*[1,1],theseylims,'--k')
    text(1000/(10*boxesmeanicethick(ii)),theseylims(1),'10H','HorizontalAlignment','left','VerticalAlignment','bottom','FontSize',axesfontsize)
    % Set up remaining:
    set(gca,'TickDir','out')
    xlim(klims)
    ylim(theseylims)
    %xlabel('Wavenumber (km^{-1})')
    set(gca,'XTickLabel',[])
    set(gca,'YTick',10.^[floor(log10(powerylims(1))):ceil(log10(powerylims(2)))])
    if ii==1
        ylabel('Power spectral density (kPa^2(km^{-1})^{-1})','FontSize',axesfontsize)
        %set(gca,'YTickLabel',10.^[floor(log10(powerylims(1))):ceil(log10(powerylims(2)))])
    else
        set(gca,'YTickLabel',[])
    end
    set(gca,'FontSize',axesfontsize)
    legend('\tau_d','\tau_b, min \lambda','\tau_b, best \lambda','\tau_b, max \lambda','location','SouthWest','FontSize',axesfontsize)
    title({[char(96+ii),') Power spectra of basal stress,'];zoomboxnames{ii}},'FontSize',titlefontsize)

    % Make coherence plot:
    % Make subplot:
    subplot('Position',Boxes{2,ii})
    % Plot coherence spectra:
    plot(Kmag_1d(1:end-1)*1000,squeeze(DragCoherence_1d(1:end-1,1,ii)),mincolor)
    hold on
    plot(Kmag_1d(1:end-1)*1000,squeeze(DragCoherence_1d(1:end-1,2,ii)),bestcolor)
    plot(Kmag_1d(1:end-1)*1000,squeeze(DragCoherence_1d(1:end-1,3,ii)),maxcolor)
    % Plot and label select wavelengths:
    theseylims=[0,1];
    plot((1000/boxesmeanicethick(ii))*[1,1],theseylims,'--k')
    text(1000/boxesmeanicethick(ii),theseylims(1),'H','HorizontalAlignment','left','VerticalAlignment','bottom','FontSize',axesfontsize)
    plot((1000/(2*boxesmeanicethick(ii)))*[1,1],theseylims,'--k')
    text(1000/(2*boxesmeanicethick(ii)),theseylims(1),'2H','HorizontalAlignment','left','VerticalAlignment','bottom','FontSize',axesfontsize)
    %plot((1000/(3*boxesmeanicethick(ii)))*[1,1],theseylims,'--k')
    %text(1000/(3*boxesmeanicethick(ii)),theseylims(1),'3H','HorizontalAlignment','left','VerticalAlignment','bottom')
    plot((1000/(5*boxesmeanicethick(ii)))*[1,1],theseylims,'--k')
    text(1000/(5*boxesmeanicethick(ii)),theseylims(1),'5H','HorizontalAlignment','left','VerticalAlignment','bottom','FontSize',axesfontsize)
    plot((1000/(10*boxesmeanicethick(ii)))*[1,1],theseylims,'--k')
    text(1000/(10*boxesmeanicethick(ii)),theseylims(1),'10H','HorizontalAlignment','left','VerticalAlignment','bottom','FontSize',axesfontsize)
    % Set up remaining:
    set(gca,'TickDir','out')
    set(gca,'XScale','log')
    xlim(klims)
    ylim([0,1])
    xlabel('Wavenumber (km^{-1})','FontSize',axesfontsize)
    if ii==1
        ylabel('Coherence (unitless)','FontSize',axesfontsize)
    else
        set(gca,'YTickLabel',[])
    end
    set(gca,'FontSize',axesfontsize)
    legend('min \lambda','best \lambda','max \lambda','location','SouthWest','FontSize',axesfontsize)
    title({[char(96+numzoomboxes+ii),') Spectral coherence of \tau_b with \tau_d,'];zoomboxnames{ii}},'FontSize',titlefontsize)

end

% Save figure:
set(gcf,'PaperSize',pagesize)
set(gcf,'PaperPosition',[0,0,pagesize])
print('-dpng',figname,['-r',num2str(resolution)])

% Final display:
disp('Done!')
toc