% MeshSizeFigure_v1

% Mike Wolovick, 1/24/2023

% This script makes a figure demonstrating mesh size for the inversion
% paper.  There are subplots for several example mesh size maps, plus a
% plot showing mesh size and number of elements as a function of mesh
% number, and a plot showing size histograms for the fast-flow grounded
% domain.

% This script automatically identifies all meshes in a size series produced
% by MakeISSMMesh_v4

clearvars
close all
tic

%% Parameters:

% File names and paths:
%generalmeshfile='~/Dropbox/FRISModeling/Meshes/FRISDomainMesh_SizeSeries_v1_mesh';
%gridfile='~/Dropbox/CombinedGrids/AntarcticaCombinedGrids_v4.mat';
%issmscriptpath='~/Dropbox/ISSM/';
%figname='~/Dropbox/FRISModeling/MeshSizeFigure_v2.png';  
generalmeshfile='/work/ollie/mwolovic/FRISModeling/Meshes/FRISDomainMesh_SizeSeries_v1_mesh';
gridfile='/work/ollie/mwolovic/Inputs/AntarcticaCombinedGrids_v4.mat';
velinputfile='/work/ollie/mwolovic/Inputs/AntarcticaMergedVel_MEASURESv2withphase_plus_FRmosaic_velonly_900m.mat'; 
issmscriptpath='/home/ollie/mwolovic/ISSM_installations/trunk-jpl/';
figname='/home/ollie/mwolovic/MeshSizeFigure_v2.png'; 

% Which meshes to show in a map?
meshes2show=[1,5,10];                % [1xN] integer

% Threshold for analyzing fast-flow areas:
velanalysisthrehsold=15;               % m/yr (smoothing weighting btwn 1/2 and 2x this value)

% Map cropping limits: (to speed up velocity interpolation)
xlims_crop=[-1800,1000]*1e3;           % [1x2] m
ylims_crop=[-300,1750]*1e3;            % [1x2] m

% Axes limits for figures:
xtick=500;                             % km
ytick=500;                             % km
meshsizeclims=[1e-1,1e2];              % [1x2] km
mediansizelims=[3e-1,1e2];             % [1x2] km
sizehistlims=[1e-1,1e2];               % [1x2] km
controlfieldlims=[0,4];                % [1x2] controlfield units
numelementlims=[6e3,2e6];              % [1x2] number
numbins=50;                            % integer

% Other aesthetic settings:
sizecmap='parula';                     % valid colormap string
sizecolor='r';                         % valid color string or RGB
numbercolor='b';                       % valid color string or RGB
groundedlinestyle='-';                 % valid linestyle string
fastlinestyle='--';                    % valid linestyle string
slowlinestyle=':';                     % valid linestyle string
groundedmarker='o';                    % valid marker string
fastmarker='s';                        % valid marker string
slowmarker='^';                        % valid marker string
histcolors={[0,0,.75];[.25,1,.25];[.75,0,0]}; % {Nx1} color strings or RGB
axesfontsize=13;                       % points
titlefontsize=14;                          % points
colorbarfontsize=13;                   % points
legendfontsize=12;                    % points

% Subplot controls:
vertaxesbuffer=.04;                    % unitless
titlebuffer=.04;                       % unitless
horzaxesbuffer=.05;                    % unitless
horznotextbuffer=.01;                  % unitless

% Output size and resolution:
pagesize=[20,10];                      % [1x2] inches (single map figures)
resolution=300;                        % dpi

% Duh parameters:
g=9.8;
secondsperyear=60*60*24*365.25;

%% Preparation:

% Communicate:
disp('Doing initial preparation and loading grids')

% Identify mesh files:
meshfiles=dir([generalmeshfile,'*.mat']);
nummeshes=length(meshfiles);

% Check number of mesh files:
if nummeshes==0
    error('Unable to locate mesh files')
elseif max(meshes2show)>nummeshes
    error('max(meshes2show) cannot be greater than number of meshes')
elseif length(meshes2show)~=length(histcolors)
    error('length(meshes2show) must be equal to length(histcolors)')
end

% Ensure that meshes are sorted correctly:
meshnumber=zeros(nummeshes,1);
for ii=1:nummeshes
    thismeshnum=strfind(meshfiles(ii).name,'mesh');
    thisdotnum=strfind(meshfiles(ii).name,'.');
    meshnumber(ii)=str2double(meshfiles(ii).name(thismeshnum+4:thisdotnum-1));
end
[~,sortedind]=sort(meshnumber);
meshfiles=meshfiles(sortedind);

% Load meshing parameters:
load([meshfiles(1).folder,'/',meshfiles(1).name],'MeshingParameters')
thisgridfile=gridfile;
thisvelfile=velinputfile;
unpack(MeshingParameters)
gridfile=thisgridfile;
velinputfile=thisvelfile;

% Compute subplot dimensions:
% Hardcoded subplot arrangement:
numvertplots=2;
numhorzplots=ceil((length(meshes2show)+2)/2);
% Compute available space:
availablehorzspace=1-(numhorzplots+1)*(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)*(horzaxesbuffer+availablehorzspace/numhorzplots),...
            vertaxesbuffer+(ii-1)*(vertaxesbuffer+titlebuffer+availablevertspace/numvertplots),...
            availablehorzspace/numhorzplots,availablevertspace/numvertplots];
    end
end

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

% Load gridded input file:
load(gridfile,'X','Y','Mask','IceThick','BedElev','SurfElev')
if dissweight~=0
    load(gridfile,"DrivingStress_alongflow")
end
if ~exist('velinputfile','var')
    load(gridfile,'U_obs','V_obs','VelMag_obs')
end

% Crop input grids:
ind1=find(Y>=ylims_crop(1),1,'first');
ind2=find(Y<=ylims_crop(2),1,'last');
ind3=find(X>=xlims_crop(1),1,'first');
ind4=find(X<=xlims_crop(2),1,'last');
X_grid=X(ind3:ind4);
Y_grid=Y(ind1:ind2);
Mask_grid=Mask(ind1:ind2,ind3:ind4);
BedElev_grid=BedElev(ind1:ind2,ind3:ind4);
SurfElev_grid=SurfElev(ind1:ind2,ind3:ind4);
IceThick_grid=IceThick(ind1:ind2,ind3:ind4);
if ~exist('velinputfile','var')
    U_obs_grid=U_obs(ind1:ind2,ind3:ind4);
    V_obs_grid=V_obs(ind1:ind2,ind3:ind4);
    VelMag_obs_grid=VelMag_obs(ind1:ind2,ind3:ind4);
end
if dissweight~=0
    DrivingStress_grid=DrivingStress_alongflow(ind1:ind2,ind3:ind4);
end
clear X Y Mask U_obs V_obs VelMag_obs BedElev IceThick SurfElev DrivingStress_alongflow

% Evaluate grid size:
[ysize,xsize]=size(Mask_grid);
dx_grid=X_grid(2)-X_grid(1);

% Reorder the mask:
ReorderedMask_grid=Mask_grid;  % 0 (ocean) and 2 (grounded ice) remain the same
ReorderedMask_grid(Mask_grid==3)=1; % shelf
ReorderedMask_grid(Mask_grid==1)=3; % rock

% Load separate velocity input file:
if exist('velinputfile','var')
    % Load:
    load(velinputfile,'X','Y','VelMag_merged','U_merged','V_merged')
    % Crop:
    ind1=find(Y>=ylims_crop(1),1,'first');
    ind2=find(Y<=ylims_crop(2),1,'last');
    ind3=find(X>=xlims_crop(1),1,'first');
    ind4=find(X<=xlims_crop(2),1,'last');
    X=X(ind3:ind4);
    Y=Y(ind1:ind2);
    U_merged=U_merged(ind1:ind2,ind3:ind4);
    V_merged=V_merged(ind1:ind2,ind3:ind4);
    VelMag_merged=VelMag_merged(ind1:ind2,ind3:ind4);
    % Smooth before interpolation:
    dx_vel=X(2)-X(1);
    if dx_grid/dx_vel>sqrt(2)
        U_merged=ThreeByThree(U_merged,round((.5*dx_grid/dx_vel)^2));
        V_merged=ThreeByThree(V_merged,round((.5*dx_grid/dx_vel)^2));
        VelMag_merged=ThreeByThree(VelMag_merged,round((.5*dx_grid/dx_vel)^2));
    end
    % Interpolate:
    U_obs_grid=interp2(X,Y,U_merged,repmat(X_grid,[ysize,1]),repmat(Y_grid,[1,xsize]));
    V_obs_grid=interp2(X,Y,V_merged,repmat(X_grid,[ysize,1]),repmat(Y_grid,[1,xsize]));
    VelMag_obs_grid=interp2(X,Y,VelMag_merged,repmat(X_grid,[ysize,1]),repmat(Y_grid,[1,xsize]));
    % Clear extraneous variables:
    clear X Y U_merged V_merged VelMag_merged
end

%% Velocity Processing:

% Communicate:
disp('Processing velocity field...')

% Record original velocity grid:
VelMag_obs_grid_original=VelMag_obs_grid;

% Smoothly fill in holes in observed velocity with neighboring values:
% identify bad values:
BadValues_grid=isnan(VelMag_obs_grid)|VelMag_obs_grid<=mingoodvalue4velfill;
% expand the bad values region by 2 in all directions:
BadValues_grid(2:end-1,2:end-1)=BadValues_grid(1:end-2,1:end-2)|BadValues_grid(2:end-1,1:end-2)|BadValues_grid(3:end,1:end-2)|...
    BadValues_grid(1:end-2,2:end-1)|BadValues_grid(2:end-1,2:end-1)|BadValues_grid(3:end,2:end-1)|...
    BadValues_grid(1:end-2,3:end)|BadValues_grid(2:end-1,3:end)|BadValues_grid(3:end,3:end);
BadValues_grid(2:end-1,2:end-1)=BadValues_grid(1:end-2,1:end-2)|BadValues_grid(2:end-1,1:end-2)|BadValues_grid(3:end,1:end-2)|...
    BadValues_grid(1:end-2,2:end-1)|BadValues_grid(2:end-1,2:end-1)|BadValues_grid(3:end,2:end-1)|...
    BadValues_grid(1:end-2,3:end)|BadValues_grid(2:end-1,3:end)|BadValues_grid(3:end,3:end);
% Flag bad values and interpolate:
U_obs_grid(BadValues_grid)=NaN;
U_obs_grid=LaplacianInterpolation(U_obs_grid,9);
V_obs_grid(BadValues_grid)=NaN;
V_obs_grid=LaplacianInterpolation(V_obs_grid,9);
VelMag_obs_grid(BadValues_grid)=NaN;
VelMag_obs_grid=LaplacianInterpolation(VelMag_obs_grid,9);

% Smooth out velocity:
if velwavelength>0
    U_obs_grid=ThreeByThree(U_obs_grid,round((.5*velwavelength/dx_grid)^2));
    V_obs_grid=ThreeByThree(V_obs_grid,round((.5*velwavelength/dx_grid)^2));
    VelMag_obs_grid=ThreeByThree(VelMag_obs_grid,round((.5*velwavelength/dx_grid)^2));
end


%% Generate control field:

% Communicate:
disp('Generating control field...')

% Initialize control field:
ControlField_grid=zeros(ysize,xsize);

% Add flux:
if fluxweight~=0
    % Compute component:
    Flux_grid=VelMag_obs_grid.*IceThick_grid;
    % Apply floor:
    Flux_grid=sqrt(Flux_grid.^2+fluxfloor^2);
    % Apply mask constraints:
    Flux_grid(ReorderedMask_grid<.5|ReorderedMask_grid>2.5)=0; % 0 outside the ice sheet or shelf
    % Add to control field:
    ControlField_grid=ControlField_grid+fluxweight*log10(Flux_grid/fluxfloor);
end

% Add dissipation:
if dissweight~=0
    % Compute component:
    Dissipation_grid=VelMag_obs_grid.*DrivingStress_grid/secondsperyear;
    % Apply mask constraints:
    Dissipation_grid(ReorderedMask_grid<1.5|ReorderedMask_grid>2.5)=0; % zero outside of grounded ice sheet
    % Apply floor:
    Dissipation_grid=sqrt(Dissipation_grid.^2+dissfloor^2);
    % Add to control field:
    ControlField_grid=ControlField_grid+dissweight*log10(Dissipation_grid/dissfloor);
end

% Add strain rate:
if strainweight~=0
    % Compute component:
    [StrainRateXX_grid,DUDY_grid]=gradient(U_obs_grid,dx_grid,dx_grid);
    [DVDX_grid,StrainRateYY_grid]=gradient(V_obs_grid,dx_grid,dx_grid);
    StrainRateXY_grid=.5*(DUDY_grid+DVDX_grid);
    StrainRateMag_grid=sqrt(StrainRateXY_grid.^2+.5*(StrainRateXX_grid.^2+StrainRateYY_grid.^2));
    % Apply mask constraints:
    StrainRateMag_grid(ReorderedMask_grid<.5|ReorderedMask_grid>2.5)=0; % 0 outside the ice sheet or shelf
    % Apply floor:
    StrainRateMag_grid=sqrt(StrainRateMag_grid.^2+strainfloor^2);
    % Add to control field:
    ControlField_grid(ReorderedMask_grid==1)=ControlField_grid(ReorderedMask_grid==1)+strainweight*log10(StrainRateMag_grid(ReorderedMask_grid==1)/strainfloor);
end

% Add velocity:
if velweight~=0
    % Apply floor:
    VelMag_obs_grid=sqrt(VelMag_obs_grid.^2+velfloor^2);
    % Apply mask constraints:
    VelMag_obs_grid(ReorderedMask_grid<.5|ReorderedMask_grid>2.5)=0; % 0 outside the ice sheet or shelf
    % Add to control field:
    ControlField_grid=ControlField_grid+velweight*log10(VelMag_obs_grid/velfloor);
end

% Add mask:
if maskweight~=0
    ControlField_grid(ReorderedMask_grid>0)=ControlField_grid(ReorderedMask_grid>0)+maskweight*(ReorderedMask_grid(ReorderedMask_grid>0)+frontmaskshift);
end

% Smooth out control field:
if controlfieldwavelength>0
    ControlField_grid=ThreeByThree(ControlField_grid,round((.5*controlfieldwavelength/dx_grid)^2));
end

%% Load model structures:

% Communicate:
disp('Loading and processing meshes...')

% Initialize model structure:
md(1:nummeshes)=model();

% Pre-allocate things:
HydraulicDiameter=cell(nummeshes,1);
Edges=linspace(log(sizehistlims(1)*1e3),log(sizehistlims(2)*1e3),numbins+1);
Counts_grounded=zeros(numbins,nummeshes);
Counts_fast=zeros(numbins,nummeshes);
Counts_slow=zeros(numbins,nummeshes);
MedianSize_grounded=zeros(nummeshes,1);
MedianSize_fast=zeros(nummeshes,1);
MedianSize_slow=zeros(nummeshes,1);
Numelements_grounded=zeros(nummeshes,1);
Numelements_fast=zeros(nummeshes,1);
Numelements_slow=zeros(nummeshes,1);

% Load and process all meshes:
for ii=1:nummeshes

    % Load this mesh:
    thisfile=load([meshfiles(ii).folder,'/',meshfiles(ii).name],'md');
    md(ii)=thisfile.md;

    % Compute mesh size:
    % Compute edge lengths:
    a=sqrt((md(ii).mesh.x(md(ii).mesh.elements(:,1))-md(ii).mesh.x(md(ii).mesh.elements(:,2))).^2+(md(ii).mesh.y(md(ii).mesh.elements(:,1))-md(ii).mesh.y(md(ii).mesh.elements(:,2))).^2);
    b=sqrt((md(ii).mesh.x(md(ii).mesh.elements(:,2))-md(ii).mesh.x(md(ii).mesh.elements(:,3))).^2+(md(ii).mesh.y(md(ii).mesh.elements(:,2))-md(ii).mesh.y(md(ii).mesh.elements(:,3))).^2);
    c=sqrt((md(ii).mesh.x(md(ii).mesh.elements(:,3))-md(ii).mesh.x(md(ii).mesh.elements(:,1))).^2+(md(ii).mesh.y(md(ii).mesh.elements(:,3))-md(ii).mesh.y(md(ii).mesh.elements(:,1))).^2);
    % Compute area:
    s=(a+b+c)/2;
    Area=sqrt(s.*(s-a).*(s-b).*(s-c));
    % Compute hydraulic diameter:
    HydraulicDiameter{ii}=2*Area./s;

    % Compute grounded area fraction:
    GroundedFrac_vertices=md(ii).mask.ice_levelset<0&md(ii).mask.ocean_levelset>0;
    GroundedFrac_elements=mean(GroundedFrac_vertices(md(ii).mesh.elements),2);

    % Interpolate velocity onto mesh:
    VelMag_vertices=MultiWavelengthInterpolator(X_grid,Y_grid,VelMag_obs_grid_original,md(ii).mesh.x,md(ii).mesh.y,md(ii).mesh.elements,'vertices');

    % Compute weighting function for fast-flowing areas:
    FastFlowWeight_vertices=ones(md(ii).mesh.numberofvertices,1);
    FastFlowWeight_vertices(VelMag_vertices<.5*velanalysisthrehsold)=0;
    FastFlowWeight_vertices(VelMag_vertices>=.5*velanalysisthrehsold&VelMag_vertices<2*velanalysisthrehsold)=...
        .5*(1+sin(.5*pi*log2(VelMag_vertices(VelMag_vertices>=.5*velanalysisthrehsold&VelMag_vertices<2*velanalysisthrehsold)/velanalysisthrehsold)));
    FastFlowWeight_elements=mean(FastFlowWeight_vertices(md(ii).mesh.elements),2);

    % Compute histograms:
    % Compute log of size:
    LogHydraulicDiameter=log(HydraulicDiameter{ii});
    % Loop through bins:
    for jj=1:numbins
        % Generate area counts:
        Counts_grounded(jj,ii)=sum(Area(LogHydraulicDiameter>=Edges(jj)&LogHydraulicDiameter<Edges(jj+1)).*...
            GroundedFrac_elements(LogHydraulicDiameter>=Edges(jj)&LogHydraulicDiameter<Edges(jj+1)));
        Counts_fast(jj,ii)=sum(Area(LogHydraulicDiameter>=Edges(jj)&LogHydraulicDiameter<Edges(jj+1)).*...
            FastFlowWeight_elements(LogHydraulicDiameter>=Edges(jj)&LogHydraulicDiameter<Edges(jj+1)).*GroundedFrac_elements(LogHydraulicDiameter>=Edges(jj)&LogHydraulicDiameter<Edges(jj+1)));
        Counts_slow(jj,ii)=sum(Area(LogHydraulicDiameter>=Edges(jj)&LogHydraulicDiameter<Edges(jj+1)).*...
            FastFlowWeight_elements(LogHydraulicDiameter>=Edges(jj)&LogHydraulicDiameter<Edges(jj+1)).*(1-GroundedFrac_elements(LogHydraulicDiameter>=Edges(jj)&LogHydraulicDiameter<Edges(jj+1))));
    end
    % Normalize area counts to probability density on log scale:
    Counts_grounded(:,ii)=Counts_grounded(:,ii)/(sum(Area.*GroundedFrac_elements)*(Edges(2)-Edges(1)));
    Counts_fast(:,ii)=Counts_fast(:,ii)/(sum(Area.*GroundedFrac_elements.*FastFlowWeight_elements)*(Edges(2)-Edges(1)));
    Counts_slow(:,ii)=Counts_slow(:,ii)/(sum(Area.*GroundedFrac_elements.*(1-FastFlowWeight_elements))*(Edges(2)-Edges(1)));

    % Compute median mesh sizes:
    % sort by size:
    [SortedSize_grounded,SortedInd_grounded]=sort(HydraulicDiameter{ii}(GroundedFrac_elements>1e-3));
    [SortedSize_fast,SortedInd_fast]=sort(HydraulicDiameter{ii}(GroundedFrac_elements>1e-3&FastFlowWeight_elements>1e-3));
    [SortedSize_slow,SortedInd_slow]=sort(HydraulicDiameter{ii}(GroundedFrac_elements>1e-3&FastFlowWeight_elements<1-1e-3));
    % Construct sorted and weighted area lists:
    SortedArea_grounded=Area(GroundedFrac_elements>1e-3).*GroundedFrac_elements(GroundedFrac_elements>1e-3);
    SortedArea_grounded=SortedArea_grounded(SortedInd_grounded);
    SortedArea_fast=Area(GroundedFrac_elements>1e-3&FastFlowWeight_elements>1e-3).*...
        GroundedFrac_elements(GroundedFrac_elements>1e-3&FastFlowWeight_elements>1e-3).*FastFlowWeight_elements(GroundedFrac_elements>1e-3&FastFlowWeight_elements>1e-3);
    SortedArea_fast=SortedArea_fast(SortedInd_fast);
    SortedArea_slow=Area(GroundedFrac_elements>1e-3&FastFlowWeight_elements<1-1e-3).*...
        GroundedFrac_elements(GroundedFrac_elements>1e-3&FastFlowWeight_elements<1-1e-3).*(1-FastFlowWeight_elements(GroundedFrac_elements>1e-3&FastFlowWeight_elements<1-1e-3));
    SortedArea_slow=SortedArea_slow(SortedInd_slow);
    % Compute cumulative area:
    CumulativeArea_grounded=cumsum(SortedArea_grounded);
    CumulativeArea_fast=cumsum(SortedArea_fast);
    CumulativeArea_slow=cumsum(SortedArea_slow);
    % Interpolate medians:
    MedianSize_grounded(ii)=interp1(CumulativeArea_grounded/CumulativeArea_grounded(end),SortedSize_grounded,.5,'linear');
    MedianSize_fast(ii)=interp1(CumulativeArea_fast/CumulativeArea_fast(end),SortedSize_fast,.5,'linear');
    MedianSize_slow(ii)=interp1(CumulativeArea_slow/CumulativeArea_slow(end),SortedSize_slow,.5,'linear');

    % Compute number of elements:
    Numelements_grounded(ii)=sum(GroundedFrac_elements);
    Numelements_fast(ii)=sum(GroundedFrac_elements.*FastFlowWeight_elements);
    Numelements_slow(ii)=sum(GroundedFrac_elements.*(1-FastFlowWeight_elements));

    % Compute things for the first mesh:
    if ii==1
        % Compute mean ice thickness:
        meanicethick=sum(Area.*GroundedFrac_elements.*mean(md(ii).geometry.thickness(md(ii).mesh.elements),2))/sum(Area.*GroundedFrac_elements);
        % Compute axes limits:
        xlims=[min(md(ii).mesh.x),max(md(ii).mesh.x)]/1000;
        ylims=[min(md(ii).mesh.y),max(md(ii).mesh.y)]/1000;
    end

end

%% Make Figure:

% Communicate:
disp('Making figure...')

% Make figure:
figure(1)

% Plot control field:
% Make subplot:
subplot('Position',Boxes{1,1})
imagesc(X_grid/1000,Y_grid/1000,ControlField_grid)
set(gca,'ydir','normal')
caxis(controlfieldlims)
hc=colorbar;
set(get(hc,'Ylabel'),'String','Arbitrary Units','FontSize',colorbarfontsize)
set(hc,'YTick',ceil(controlfieldlims(1)):floor(controlfieldlims(2)))
xlim(xlims)
ylim(ylims)
set(gca,'XTick',xtick*[ceil(xlims(1)/xtick):1:floor(xlims(2)/xtick)])
set(gca,'YTick',ytick*[ceil(ylims(1)/ytick):1:floor(ylims(2)/ytick)])
set(gca,'PlotBoxAspectRatio',[(xlims(2)-xlims(1))/(ylims(2)-ylims(1)),1,1])
xlabel('Easting (km)','FontSize',axesfontsize)
ylabel('Northing (km)','FontSize',axesfontsize)
set(gca,'FontSize',axesfontsize)
title('a) Control field for mesh refinement','FontSize',titlefontsize)


% Plot size line plots:
% make subplot:
subplot('Position',Boxes{1,2})
% Do first y-axis:
yyaxis left
% Plot sizes:
plot(linspace(1,nummeshes,nummeshes)',MedianSize_grounded/1e3,'Color',sizecolor,'LineStyle',groundedlinestyle,'Marker',groundedmarker,'MarkerFaceColor',sizecolor)
hold on
plot(linspace(1,nummeshes,nummeshes)',MedianSize_fast/1e3,'Color',sizecolor,'LineStyle',fastlinestyle,'Marker',fastmarker,'MarkerFaceColor',sizecolor)
plot(linspace(1,nummeshes,nummeshes)',MedianSize_slow/1e3,'Color',sizecolor,'LineStyle',slowlinestyle,'Marker',slowmarker,'MarkerFaceColor',sizecolor)
% plot ice thick:
%plot(linspace(1,nummeshes,nummeshes)',meanicethick*ones(nummeshes,1)/1e3,'Color',sizecolor,'LineStyle',':')
% Set up axes etc:
set(gca,'Yscale','log')
xlim([1,nummeshes])
ylim(mediansizelims)
xlabel('Mesh Number','FontSize',axesfontsize)
ylabel('Median element size (km)','FontSize',axesfontsize)
set(gca,'FontSize',axesfontsize)
set(gca,'xcolor','k')
set(gca,'ycolor',sizecolor)
% Make title:
title('b) Mesh size metrics','FontSize',titlefontsize)
% Do second y-axis:
yyaxis right
% Plot sizes:
plot(linspace(1,nummeshes,nummeshes)',Numelements_grounded,'Color',numbercolor,'LineStyle',groundedlinestyle,'Marker',groundedmarker,'MarkerFaceColor',numbercolor)
hold on
plot(linspace(1,nummeshes,nummeshes)',Numelements_fast,'Color',numbercolor,'LineStyle',fastlinestyle,'Marker',fastmarker,'MarkerFaceColor',numbercolor)
plot(linspace(1,nummeshes,nummeshes)',Numelements_slow,'Color',numbercolor,'LineStyle',slowlinestyle,'Marker',slowmarker,'MarkerFaceColor',numbercolor)
% Set up axes etc:
set(gca,'Yscale','log')
xlim([1,nummeshes])
ylim(numelementlims)
xlabel('Mesh Number','FontSize',axesfontsize)
ylabel('Number of elements','FontSize',axesfontsize)
set(gca,'xcolor','k')
set(gca,'ycolor',numbercolor)
set(gca,'FontSize',axesfontsize)
% Make legend:
h1=plot(-9999*[1,1],9999*[1,1],'Color','k','LineStyle',groundedlinestyle,'Marker',groundedmarker,'MarkerFaceColor','k');
h2=plot(-9999*[1,1],9999*[1,1],'Color','k','LineStyle',fastlinestyle,'Marker',fastmarker,'MarkerFaceColor','k');
h3=plot(-9999*[1,1],9999*[1,1],'Color','k','LineStyle',slowlinestyle,'Marker',slowmarker,'MarkerFaceColor','k');
%h3=plot(-9999*[1,1],9999*[1,1],'Color',sizecolor,'LineStyle',':');
legend([h1;h2;h3],'whole domain','fast-flow','slow-flow','location','SouthWest','FontSize',legendfontsize)


% Make histograms:
% make subplot:
subplot('Position',Boxes{1,3})
% Loop through meshes to show:
legendtext=cell(length(meshes2show),1);
for ii=1:length(meshes2show)
    % Make histogram:
    histogram('BinEdges',exp(Edges)/1e3,'BinCounts',Counts_fast(:,meshes2show(ii))','DisplayStyle','stairs','EdgeColor',histcolors{ii})
    hold on
    % Generate legend text:
    legendtext{ii}=['Mesh ',num2str(meshes2show(ii))];
end
% Set up axes etc:
set(gca,'YAxisLocation','right')
set(gca,'XScale','log')
xlim(sizehistlims.*[.999,1.001])
xlabel('Hydraulic diameter (km)','FontSize',axesfontsize)
ylabel('Probability Density','FontSize',axesfontsize)
set(gca,'FontSize',axesfontsize)
legend(legendtext,'location','NorthEast','FontSize',legendfontsize)
title('c) Select element size distributions','FontSize',titlefontsize)



% Loop through meshes to display:
for ii=1:length(meshes2show)
    % Make this subplot:
    subplot('Position',Boxes{2,ii})
    % Make the patch:
    patch('Vertices',[md(meshes2show(ii)).mesh.x/1000,md(meshes2show(ii)).mesh.y/1000],'Faces',md(meshes2show(ii)).mesh.elements,'EdgeColor','none','FaceColor','flat',...
        'FaceVertexCData',log10(HydraulicDiameter{meshes2show(ii)}/1000))
    theseclims=log10(meshsizeclims);
    caxis(theseclims)
    colormap(sizecmap)
    
    % Make colorbar:
    hc=colorbar;
    set(get(hc,'Ylabel'),'String','Hydraulic Diameter (km)','FontSize',colorbarfontsize)
    % Compute first tick for color bar:
    thisbase=10^floor(theseclims(1));
    thistick=thisbase*ceil((10^theseclims(1))/thisbase);
    % Create all ticks and labels for color bar:
    cticks=zeros(100,1);
    clabels=cell(100,1);
    for jj=1:100
        % Record this tick:
        cticks(jj)=log10(thistick);
        % Check if we hit the next power of ten:
        if abs(thistick-10*thisbase)/thisbase<1e-12 || (jj==1 && abs(thistick-thisbase)/thisbase<1e-12)
            % Generate label:
            clabels{jj}=['10^{',num2str(round(log10(thistick))),'}'];
            % Advance base:
            thisbase=thistick;
        end
        % Advance tick:
        thistick=thistick+thisbase;
        % Break from loop and truncate:
        if thistick>(1+1e-12)*(10^theseclims(2))
            cticks=cticks(1:jj);
            clabels=clabels(1:jj);
            break
        end
    end
    % Set the colorbar ticks and labels:
    set(hc,'YTick',cticks)
    set(hc,'YTickLabel',clabels,'FontSize',colorbarfontsize)
    set(hc,'TickDir','out')
    
    % Fix up the rest:
    xlim(xlims)
    ylim(ylims)
    set(gca,'XTick',xtick*[ceil(xlims(1)/xtick):1:floor(xlims(2)/xtick)])
    set(gca,'YTick',ytick*[ceil(ylims(1)/ytick):1:floor(ylims(2)/ytick)])
    set(gca,'PlotBoxAspectRatio',[(xlims(2)-xlims(1))/(ylims(2)-ylims(1)),1,1])
    xlabel('Easting (km)','FontSize',axesfontsize)
    if ii==1
        ylabel('Northing (km)','FontSize',axesfontsize)
    else
        set(gca,'YTickLabel',[])
    end
    set(gca,'FontSize',axesfontsize)
    title([char(96+3+ii),') Element size map for mesh number ',num2str(meshes2show(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

