% InversionComparisonFigure_v1

% Mike Wolovick, 12/9/2022

% This script makes a figure comparing the inversion results for best-value
% lambda, max lambda, and min lambda.  For each lambda value, it shows drag
% coefficient, drag, and velocity misfit.  There are also zoom plots for a
% few key areas.

% v1a:  outlines of zoom boxes can be skewed rectangles, and the zoom boxes
% themselves can be disabled.

% _N_v1:  compares results for different effective pressure fields.  Makes
% a table with stats about the comparison, and makes a figure showing
% comparative L-curves.

% _m_v1:  Compares results for different values of m.  The L-curve figure
% also shows the misfit components.  The map just shows drag coeff, with
% automatically scaled clims

clearvars
tic
close all

%% Parameters:

% File names:
% lcurvefiles={'/albedo/work/user/mwolovic/FRISModeling/LCurves/FRISDomainLCurve_ResTest_SSA_drag_v2_mesh1a.mat';...
%     '/albedo/work/user/mwolovic/FRISModeling/LCurves/FRISDomainLCurve_ResTest_SSA_drag_op_v1_mesh1.mat';...
%     '/albedo/work/user/mwolovic/FRISModeling/LCurves/FRISDomainLCurve_ResTest_SSA_drag_opc_v1_mesh1.mat'};  % cell array of strings, include trailing ".mat", first one is reference
lcurvefiles={'/work/ollie/mwolovic/FRISModeling/LCurves/FRISDomainLCurve_ResTest_SSA_drag_v2_mesh1a.mat';...
    '/work/ollie/mwolovic/FRISModeling/LCurves/FRISDomainLCurve_ResTest_SSA_drag_v2a_m3_mesh1.mat';...
    '/work/ollie/mwolovic/FRISModeling/LCurves/FRISDomainLCurve_ResTest_SSA_drag_v2b_m5_mesh1.mat';...
    '/work/ollie/mwolovic/FRISModeling/LCurves/FRISDomainLCurve_ResTest_SSA_cuas_drag_v1a_mesh1.mat';...
    '/work/ollie/mwolovic/FRISModeling/LCurves/FRISDomainLCurve_ResTest_SSA_cuas_drag_v1a_m3_mesh1.mat';...
    '/work/ollie/mwolovic/FRISModeling/LCurves/FRISDomainLCurve_ResTest_SSA_cuas_drag_v1b_m5_mesh1.mat'};  % cell array of strings, include trailing ".mat", first one is reference
localpath='/home/mvolovic/Dropbox/FRISModeling/LCurves/';
cmapfiles={'/home/ollie/mwolovic/GMA_default.mat';...
    '/home/ollie/mwolovic/GMA_combined.mat';...
    '/home/ollie/mwolovic/GMA_ocean.mat';...
    '/home/ollie/mwolovic/GMA_land.mat'};  % cell array of strings (can be commented)
% cmapfiles={'~/Dropbox/GMA_default.mat';...
%     '~/Dropbox/GMA_combined.mat';...
%     '~/Dropbox/GMA_ocean.mat';...
%     '~/Dropbox/GMA_land.mat'};  % cell array of strings (can be commented)
%issmscriptpath='/albedo/home/mwolovic/ISSM_installations/trunk-jpl/';
% issmscriptpath='~/Dropbox/ISSM/';
% myscriptpath='~/Dropbox/MATLAB/';
issmscriptpath='/home/ollie/mwolovic/ISSM_installations/trunk-jpl/';
myscriptpath='/home/ollie/mwolovic/MATLAB/';
%figname='/albedo/home/mwolovic/InversionComparisonFigure_v4.png';
%figname_lcurves='~/Dropbox/FRISModeling/LCurveComparisonFigure_N_v1.png';
figname_lcurves='/home/ollie/mwolovic/LCurveComparisonFigure_m_v2.png';
figname_maps='/home/ollie/mwolovic/InversionComparisonFigure_m_v2.png';

% N descriptions:
descriptext={'m=1, no N','m=3, no N','m=5, no N','m=1, N_{cuas}','m=3, N_{cuas}','m=5, N_{cuas}'};  % cell array of strings
%descriptext={'m=1','m=3','m=5'};  % cell array of strings

% 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.

% Overall fig-making controls:
makelcurvefig=0;                          % logical
makemapfig=0;                             % logical

% Analysis controls:
velanalysisthrehsold=15;               % m/yr (smoothing weighting btwn 1/2 and 2x this value)
minvelstabilizer=1e-3;                    % m/yr
mintaustabilizer=1;                        % Pa

% Ticks and limits:
jreglims=[1e-5,1e-1];                    % [1x2] unitless
jobslims=[1e-4,1e0];                    % [1x2] unitless
xtick=500;                               % km
ytick=500;                               % km
scalebarlength=50;                % km (km of the zoom plot)
scalebarheight=40;                  % km (km of the main plot)
dragcoeffclims=[1,99];            % [Nx2] percentiles
dragcoeffcmap='GMA_ocean';  % valid colormap string, potentially including names of cmapfiles
flipdragcoeffcmap=0;               % logical
draglims=[0,300];                   % [1x2] kPa
dragtick=50;                             % kPa
dragcmap='GMA_default';         % valid colormap string, potentially including names of cmapfiles
flipdragcmap=0;                        % logical
dragdifflims=[-100,100];               % [1x2] kPa
dragdifftick=25;                         % kPa
dragdiffcmap='FadedJet';           % valid colormap string, potentially including names of cmapfiles

% Zoom plot controls:
showzoombox=1;                          % logical
showzoomboxloc=1;                     % logical
zoomratio=[5,4,5,5];                    % [Nx1] unitless 
zoomboxes=[-504,-539,-384,-349,934,1003,1082,1013;...
    -490,-533,-346,-295,753,849,954,872;...
    -600,-600,-466,-466,274,348,348,274;...
    -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={'Slessor';'Recovery';'Academy';'Rutford'};  % Nx1 cell array of strings

% Other aesthetic settings:
lcurvecolors={'k',[0,.75,0],'m'};               % {1xN} valid color strings
modelmarker='o';                     % valid marker string
modelmarkersize=5;                   % points
curvelinesize=1;                     % points
cornerlinesize=5;                    % points
curvelinestyle='-';                  % valid linestyle string
cornerlinestyle='-';                  % valid linestyle string
bestmarker='d';                      % valid marker string
bestmarkersize=10;                    % points
totcurvecolor='k';                   % valid color string
datacurvecolor='b';                   % valid color string
regcurvecolor='r';                   % valid color string
backplanecolor=[.75,.75,.75];       % valid color string or RGB
titlefontsize=12;                              % points
axesfontsize=10;                              % points
colorbarfontsize=10;                         % points
zoomboxnamefontsize=9;                % points
zoomboxscalebarfontsize=5;                % points
zoomboxnamemargin=.5;                % pixel units
zoomboxnamebuffer=15;                  % km (km of the main plot, offset from upper left corner)
zoomboxscalebarbufferx=[320,370,220,50]; % [1x3] km (km of the main plot, offset from lower right corner)
zoomboxscalebarbuffery=[10,10,10,219]; % [1x3] km (km of the main plot, offset from lower right corner)
outlinestripewidth=25;                          % km (km of the main plot, only for zoom plot 1)

% Subplot controls:
vertaxesbuffer_map=.04;             % unitless
titlebuffer_map=.03;                     % unitless
horzaxesbuffer_map=.04;             % unitless
horznotextbuffer_map=.05;          % unitless
horzcolorbarbuffer_map=.09;        % unitless
vertaxesbuffer_lcurve=.05;          % unitless
horzaxesbuffer_lcurve=.04;          % unitless
titlebuffer_lcurve=.04;             % unitless
horznotextbuffer_lcurve=.02;        % unitless
lcurvefrac=.5;                      % unitless<1

% Output size and resolution:
lcurvepagesize=[20,15];           % [1x2] inches
pagewidth_map=18;                       % inches (height computed automatically)
resolution=300;                     % dpi

%% Preparation:

% Check number of N:
numlcurves=length(lcurvefiles);
if length(descriptext)~=numlcurves
    error('Length of "ndescriptext" did not match length of "lcurvefiles".')
end

% Compute boxes for L-Curve figure:
numsidevertplots=numlcurves;
numsidehorzplots=2;
Boxes_lcurves=cell(numsidevertplots*numsidehorzplots+1,1);
Boxes_lcurves{1}=[horzaxesbuffer_lcurve,vertaxesbuffer_lcurve,lcurvefrac-horzaxesbuffer_lcurve-horznotextbuffer_lcurve,1-titlebuffer_lcurve-vertaxesbuffer_lcurve];
availablehorzspace=1-lcurvefrac-numsidehorzplots*(horzaxesbuffer_lcurve+horznotextbuffer_lcurve);
availablevertspace=1-numsidevertplots*(vertaxesbuffer_lcurve+titlebuffer_lcurve);
for ii=1:numsidevertplots*numsidehorzplots
    d1=ceil(ii/2);
    d2=ii-(d1-1)*2;
    Boxes_lcurves{ii+1}=[lcurvefrac+d2*horzaxesbuffer_lcurve+(d2-1)*(horznotextbuffer_lcurve+availablehorzspace/numsidehorzplots),...
        vertaxesbuffer_lcurve+(numsidevertplots-d1)*(vertaxesbuffer_lcurve+titlebuffer_lcurve+availablevertspace/numsidevertplots),availablehorzspace/numsidehorzplots,availablevertspace/numsidevertplots];
end

% Compute boxes for map figure:
% Hardcoded subplot arrangement:
numvertplots=ceil(numlcurves/3);
numhorzplots=3;
% Compute subplot dimensions:
% Compute available space:
availablehorzspace=1-(numhorzplots)*horznotextbuffer_map-horzaxesbuffer_map;
availablevertspace=1-numvertplots*titlebuffer_map-vertaxesbuffer_map;
% Compute boxes:
Boxes_maps=cell(numvertplots,numhorzplots);
for ii=1:numvertplots
    for jj=1:numhorzplots
        Boxes_maps{numvertplots-ii+1,jj}=[horzaxesbuffer_map+(jj-1)*(horznotextbuffer_map+availablehorzspace/numhorzplots),...
            vertaxesbuffer_map+(ii-1)*(titlebuffer_map+availablevertspace/numvertplots),...
            availablehorzspace/numhorzplots,availablevertspace/numvertplots];
    end
end

% Check number of zoom boxes:
numzoomboxes=size(zoomboxes,1);
if numzoomboxes>4
    error('Sorry I dont know how to place more than 4 zoom boxes')
end

% Load colormaps:
if exist('cmapfiles','var')
    for ii=1:length(cmapfiles)
        % load the colormap:
        load(cmapfiles{ii},'thiscolormap')
        % rename the colormap:
        slashnum=strfind(cmapfiles{ii},'/');
        thiscmapname=cmapfiles{ii}(slashnum(end)+1:end-4);
        eval([thiscmapname,'=thiscolormap;'])
    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(lcurvefiles{1},'/');
    remotepath=lcurvefiles{1}(1:slashnum(end));
    lcurvename=lcurvefiles{1}(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(lcurvefiles{1},'BasicModel','MyInversionParameters','*best*','Lambda')
end

% Rename model structure:
md=BasicModel;

% Unpack inversion parameters:
unpack(MyInversionParameters)

% 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)))];

% Generate the active mask:
ActiveMask_drag=IceLevelSet_vertices<=0&OceanLevelSet_vertices>=0&md.mesh.vertexonboundary(1:numvertices2d)==0;
ActiveMask_elements=mean(ActiveMask_drag(elements2d),2);

% Unpack the good data mask:
GoodData_obs_vertices=md.miscellaneous.dummy(1:numvertices2d);
GoodDataWeight_elements=mean(GoodData_obs_vertices(elements2d).^2,2);

% Extract observed velocity:
U_obs_vertices=md.inversion.vx_obs(end-numvertices2d+1:end);
V_obs_vertices=md.inversion.vy_obs(end-numvertices2d+1:end);
VelMag_obs_vertices=md.inversion.vel_obs(end-numvertices2d+1:end);

% Compute weighting function for fast-flowing areas:
FastFlowWeight_vertices=ones(numvertices2d,1);
FastFlowWeight_vertices(VelMag_obs_vertices<.5*velanalysisthrehsold)=0;
FastFlowWeight_vertices(VelMag_obs_vertices>=.5*velanalysisthrehsold&VelMag_obs_vertices<2*velanalysisthrehsold)=...
    .5*(1+sin(.5*pi*log2(VelMag_obs_vertices(VelMag_obs_vertices>=.5*velanalysisthrehsold&VelMag_obs_vertices<2*velanalysisthrehsold)/velanalysisthrehsold)));
FastFlowWeight_elements=mean(FastFlowWeight_vertices(elements2d),2);

% Compute page height:
pageheight_map=pagewidth_map*((ylims(2)-ylims(1))/(xlims(2)-xlims(1)))*(availablehorzspace/availablevertspace)*(numvertplots/numhorzplots);

% Ensure exact rectangles:
if showzoombox || showzoomboxloc
    % 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)]];
    % fixed corner:
    %     zoomboxes=[repmat(zoomboxes(:,1),[1,4])+[zeros(numzoomboxes,1),MeanWidth.*(-imag(MeanHeadingVector)),MeanWidth.*(-imag(MeanHeadingVector))+MeanLength.*real(MeanHeadingVector),MeanLength.*real(MeanHeadingVector)],...
    %         repmat(zoomboxes(:,5),[1,4])+[zeros(numzoomboxes,1),MeanWidth.*real(MeanHeadingVector),MeanWidth.*real(MeanHeadingVector)+MeanLength.*imag(MeanHeadingVector),MeanLength.*imag(MeanHeadingVector)]];
end

% Check whether to show the zoom boxes themselves:
if showzoombox
    % 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 zoom box placement:
    zoomboxes_big=zeros(numzoomboxes,4);
    zoomboxes_big(1,:)=[xlims(1),ylims(2)-zoomratio(1)*MeanWidth(1),xlims(1)+zoomratio(1)*MeanLength(1),ylims(2)];
    if numzoomboxes>=2
        zoomboxes_big(2,:)=[xlims(2)-zoomratio(2)*MeanLength(2),ylims(2)-zoomratio(2)*MeanWidth(2),xlims(2),ylims(2)];
    end
    if numzoomboxes>=3
        zoomboxes_big(3,:)=[xlims(2)-zoomratio(3)*MeanLength(3),ylims(1),xlims(2),ylims(1)+zoomratio(3)*MeanWidth(3)];
    end
    if numzoomboxes>=4
        zoomboxes_big(4,:)=[xlims(1),zoomboxes_big(1,2)-zoomratio(4)*MeanWidth(4),xlims(1)+zoomratio(4)*MeanLength(4),zoomboxes_big(1,2)];
    end
    % Note: zoomboxes_big are [xmin,ymin,xmax,ymax], because the big boxes
    % are rectangles aligned with the coordinate axes.
end

%% Collect Data:

% Pre-allocate:
EffectivePressure_vertices=ones(numvertices2d,numlcurves);
EffectivePressure_elements=ones(numelements2d,numlcurves);
DragCoeff_vertices=zeros(numvertices2d,numlcurves);
DragCoeff_elements=zeros(numelements2d,numlcurves);
Drag_vertices=zeros(numvertices2d,numlcurves);
Drag_elements=zeros(numelements2d,numlcurves);
VelMag_vertices=zeros(numvertices2d,numlcurves);   
VelMag_elements=zeros(numelements2d,numlcurves);  
VelMagMisfit_vertices=zeros(numvertices2d,numlcurves);   
VelMagMisfit_elements=zeros(numelements2d,numlcurves);  
ObsMisfit_vertices=zeros(numvertices2d,numlcurves);   
ObsMisfit_elements=zeros(numelements2d,numlcurves);   
RegMisfit_elements=zeros(numelements2d,numlcurves);  % no interpolation from vertices to elements required for reg term
BestLambda=zeros(numlcurves,1);
MaxLambda=zeros(numlcurves,1);
MinLambda=zeros(numlcurves,1);
DeltaLambdaRatio=zeros(numlcurves,1);
DragCoeffScale=zeros(numlcurves,1);
VelScale=zeros(numlcurves,1);
InversionParameters=cell(numlcurves,1);
LCurveData=struct('LogLambda_points',[],'LogRegCost_points',[],'LogDataCost_points',[],'IsOutlier',[],'LogRegCost_curve',[],'LogDataCost_curve',[],...
    'LogDataCurvature_curve',[],'LogRegCurvature_curve',[],'TotalCurvature_curve',[],...
    'bestlambdaind_maxallowable',[],'bestlambdaind_minallowable',[],'bestlambdaind',[]);

% Loop through models:
for thismodel=1:numlcurves

    % Check if we're loading remote files:
    if exist('serverloc','var')
        % Parse filename:
        slashnum=strfind(lcurvefiles{thismodel},'/');
        remotepath=lcurvefiles{thismodel}(1:slashnum(end));
        lcurvename=lcurvefiles{thismodel}(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 basic file:
        load([localpath,lcurvename],'Log*Cost_bestcurve','*Curvature_bestcurve','Lambda','DataCost','RegularizationCost','IsOutlier','bestlambda*','bestmodelind','MyInversionParameters')
    else
        % Load the basic file:
        load(lcurvefiles{thismodel},'Log*Cost_bestcurve','*Curvature_bestcurve','Lambda','DataCost','RegularizationCost','IsOutlier','bestlambda*','bestmodelind','MyInversionParameters')
    end
    
    % Unpack parameters:
    unpack(MyInversionParameters)
    
    % Record parameters:
    InversionParameters{thismodel}=MyInversionParameters;

    % Construct this file name:
    slashnum=strfind(lcurvefiles{thismodel},'/');
    thisbase=10^floor(log10(Lambda(bestmodelind)));
    thisname=[lcurvefiles{thismodel}(slashnum(end)+1:end-4),'_Lambda_',num2str(thisbase*.01*round(100*Lambda(bestmodelind)/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([lcurvefiles{thismodel}(1:slashnum(end)),thisname],'md')
    end

    % Assign simple things:
    BestLambda(thismodel)=bestlambda;
    MaxLambda(thismodel)=bestlambda_maxallowable;
    MinLambda(thismodel)=bestlambda_minallowable;
    DeltaLambdaRatio(thismodel)=exp(bestlambda_loguncertainty);
    DragCoeffScale(thismodel)=dragcoeffstdev;
    
    % Assign velocity scale:
     if strcmp(velabsscale,'mean')
         VelScale(thismodel)=velmean;
    elseif strcmp(velabsscale,'median')
        VelScale(thismodel)=velmedian;
    elseif strcmp(velabsscale,'stdev')
        VelScale(thismodel)=velstdev;
    elseif strcmp(velabsscale,'rms')
        VelScale(thismodel)=velrms;
     else
        error('I do not recognize the value of "velabsscale".')
    end
    
    % Assign lcurve data:
    LCurveData(thismodel).LogLambda_points=log(Lambda);
    LCurveData(thismodel).LogRegCost_points=log(RegularizationCost);
    LCurveData(thismodel).LogDataCost_points=log(DataCost);
    LCurveData(thismodel).IsOutlier=IsOutlier;
    LCurveData(thismodel).LogRegCost_curve=LogRegCost_bestcurve;
    LCurveData(thismodel).LogDataCost_curve=LogDataCost_bestcurve;
    LCurveData(thismodel).LogDataCurvature_curve=LogDataCurvature_bestcurve;
    LCurveData(thismodel).LogRegCurvature_curve=LogRegCurvature_bestcurve;
    LCurveData(thismodel).TotalCurvature_curve=TotalCurvature_bestcurve;
    LCurveData(thismodel).bestlambdaind=bestlambdaind;
    LCurveData(thismodel).bestlambdaind_maxallowable=bestlambdaind_maxallowable;
    LCurveData(thismodel).bestlambdaind_minallowable=bestlambdaind_minallowable;

    % Check mesh size:
    if numlayers>1 && ~strcmp(flowequation,'SSA')
        if md.mesh.numberofvertices2d~=numvertices2d || md.mesh.numberofelements2d~=numelements2d
            error(['Mesh number ',num2str(thismodel),' does not match reference mesh.'])
        end
    else
        if md.mesh.numberofvertices~=numvertices2d || md.mesh.numberofelements~=numelements2d
            error(['Mesh number ',num2str(thismodel),' does not match reference mesh.'])
        end
    end

    % Extract variables we need:
    if numlayers>1 && ~strcmp(flowequation,'SSA')
        ModelSurfU=md.results.StressbalanceSolution.Vx(md.mesh.vertexonsurface==1);
        ModelSurfV=md.results.StressbalanceSolution.Vy(md.mesh.vertexonsurface==1);
        ModelSurfVel=md.results.StressbalanceSolution.Vel(md.mesh.vertexonsurface==1);
        ModelBasalVel=md.results.StressbalanceSolution.Vel(md.mesh.vertexonbase==1);
        DragCoeff_vertices(:,thismodel)=md.results.StressbalanceSolution.FrictionCoefficient(md.mesh.vertexonbase==1);
        if ~strcmp(nsource,'none')
            EffectivePressure_vertices(:,thismodel)=md.friction.effective_pressure(md.mesh.vertexonbase==1);
        else
            EffectivePressure_vertices(:,thismodel)=ones(numvertices2d,1);
        end
    else
        ModelBasalVel=md.results.StressbalanceSolution.Vel;
        ModelSurfVel=ModelBasalVel;
        ModelSurfU=md.results.StressbalanceSolution.Vx;
        ModelSurfV=md.results.StressbalanceSolution.Vy;
        DragCoeff_vertices(:,thismodel)=md.results.StressbalanceSolution.FrictionCoefficient;
        if ~strcmp(nsource,'none')
            EffectivePressure_vertices(:,thismodel)=md.friction.effective_pressure;
        else
            EffectivePressure_vertices(:,thismodel)=ones(numvertices2d,1);
        end
    end
    
    % Interpolate cost function coefficients to elements:
    AreReg=costfunctions>500;
    CostFunctionCoeffs_elements=zeros(numelements2d,length(costfunctions));
    for jj=1:length(costfunctions)
        CostFunctionCoeffs_elements(:,jj)=mean(reshape(md.inversion.cost_functions_coefficients(elements2d,jj),[numelements2d,3]),2);
        if AreReg(jj)
            CostFunctionCoeffs_elements(:,jj)=CostFunctionCoeffs_elements(:,jj)/Lambda(bestmodelind);
        end
    end
    
    % Compute regularization term:
    % Compute gradients:
    DragCoeffGradX_elements=sum(BasisGradX_elements.*reshape(DragCoeff_vertices(elements2d,thismodel),[numelements2d,3]),2);
    DragCoeffGradY_elements=sum(BasisGradY_elements.*reshape(DragCoeff_vertices(elements2d,thismodel),[numelements2d,3]),2);
    % Compute (squared) magnitude and apply coefficients:
    RegMisfit_elements(:,thismodel)=.5*(DragCoeffGradX_elements.^2+DragCoeffGradY_elements.^2).*CostFunctionCoeffs_elements(:,AreReg);

    % Replace zeros:
    DragCoeff_vertices(:,thismodel)=max(mindragcoeff,DragCoeff_vertices(:,thismodel));

    % Compute drag:
    Drag_vertices(:,thismodel)=(DragCoeff_vertices(:,thismodel).^2).*((ModelBasalVel/secondsperyear).^(1/m));
    if ~strcmp(nsource,'none')
        Drag_vertices(:,thismodel)=Drag_vertices(:,thismodel).*(EffectivePressure_vertices(:,thismodel).^r);
    end
    % Set drag to zero under floating shelf and in ice-free areas:
    Drag_vertices(OceanLevelSet_vertices<0|IceLevelSet_vertices>0,thismodel)=0;

    % Compute velocity magnitude misfit:
    VelMagMisfit_vertices(:,thismodel)=ModelSurfVel-VelMag_obs_vertices;
    
    % Compute observational misfit:
    ObsMisfit_vertices(:,thismodel)=.5*((ModelSurfU-U_obs_vertices).^2+(ModelSurfV-V_obs_vertices).^2)/(secondsperyear^2);
    
    % Assign (sliding) velocity:
    VelMag_vertices(:,thismodel)=ModelBasalVel;

    % Interpolate to elements:
    EffectivePressure_elements(:,thismodel)=mean(reshape(EffectivePressure_vertices(elements2d,thismodel),[numelements2d,3]),2);
    DragCoeff_elements(:,thismodel)=mean(reshape(DragCoeff_vertices(elements2d,thismodel),[numelements2d,3]),2);
    Drag_elements(:,thismodel)=mean(reshape(Drag_vertices(elements2d,thismodel),[numelements2d,3]),2);
    VelMag_elements(:,thismodel)=mean(reshape(VelMag_vertices(elements2d,thismodel),[numelements2d,3]),2);
    VelMagMisfit_elements(:,thismodel)=mean(reshape(VelMagMisfit_vertices(elements2d,thismodel),[numelements2d,3]),2);
    ObsMisfit_elements(:,thismodel)=mean(reshape(ObsMisfit_vertices(elements2d,thismodel),[numelements2d,3]),2).*CostFunctionCoeffs_elements(:,~AreReg);  % apply coefficients as well

end

% Convert to ln(C) and ln(N):
LnC_elements=2*log(DragCoeff_elements);
LnN_elements=log(EffectivePressure_elements);
LnU_elements=log(max(VelMag_elements,minvelstabilizer));
LnTau_elements=log(max(Drag_elements,mintaustabilizer));

%% Compute stats:

% Compute total area for stats computation:
totalstatsarea=[sum(Area.*ActiveMask_elements.*GoodDataWeight_elements),sum(Area.*ActiveMask_elements.*FastFlowWeight_elements.*GoodDataWeight_elements)];

% Pre-allocate:
meanc=zeros(numlcurves,2);
rmsc=zeros(numlcurves,2);
meanlnc=zeros(numlcurves,2);
rmslnc=zeros(numlcurves,2);
meann=zeros(numlcurves,2);
rmsn=zeros(numlcurves,2);
meanlnn=zeros(numlcurves,2);
rmslnn=zeros(numlcurves,2);
meantau=zeros(numlcurves,2);
rmstau=zeros(numlcurves,2);
meanlntau=zeros(numlcurves,2);
rmslntau=zeros(numlcurves,2);
meanlnu=zeros(numlcurves,2);
rmslnu=zeros(numlcurves,2);
rlnclnc=zeros(numlcurves,2);
rnc=zeros(numlcurves,2);
rlnnlnc=zeros(numlcurves,2);
rtautau=zeros(numlcurves,2);
rntau=zeros(numlcurves,2);
rmslnclnc=zeros(numlcurves,2);
rmstautau=zeros(numlcurves,2);
rmsvelmisfit=zeros(numlcurves,2);
jobs=zeros(numlcurves,2);
jreg=zeros(numlcurves,2);
jobsvelequiv=zeros(numlcurves,2);

% for all stats variables, first column is whole domain and second column
% is fast-flow domain.

% Loop through models:
for ii=1:numlcurves
    
    % Compute basic means and variances:
    % linear C:  (coorected from k to C)
    meanc(ii,1)=sum(Area.*ActiveMask_elements.*GoodDataWeight_elements.*(DragCoeff_elements(:,ii).^2))/totalstatsarea(1);
    rmsc(ii,1)=sqrt(sum(Area.*ActiveMask_elements.*GoodDataWeight_elements.*(((DragCoeff_elements(:,ii).^2)-meanc(ii,1)).^2))/totalstatsarea(1));
    meanc(ii,2)=sum(Area.*ActiveMask_elements.*FastFlowWeight_elements.*GoodDataWeight_elements.*(DragCoeff_elements(:,ii).^2))/totalstatsarea(2);
    rmsc(ii,2)=sqrt(sum(Area.*ActiveMask_elements.*FastFlowWeight_elements.*GoodDataWeight_elements.*(((DragCoeff_elements(:,ii).^2)-meanc(ii,2)).^2))/totalstatsarea(2));
    % lnC:
    meanlnc(ii,1)=sum(Area.*ActiveMask_elements.*GoodDataWeight_elements.*LnC_elements(:,ii))/totalstatsarea(1);
    rmslnc(ii,1)=sqrt(sum(Area.*ActiveMask_elements.*GoodDataWeight_elements.*((LnC_elements(:,ii)-meanlnc(ii,1)).^2))/totalstatsarea(1));
    meanlnc(ii,2)=sum(Area.*ActiveMask_elements.*FastFlowWeight_elements.*GoodDataWeight_elements.*LnC_elements(:,ii))/totalstatsarea(2);
    rmslnc(ii,2)=sqrt(sum(Area.*ActiveMask_elements.*FastFlowWeight_elements.*GoodDataWeight_elements.*((LnC_elements(:,ii)-meanlnc(ii,2)).^2))/totalstatsarea(2));
    % linear N:
    meann(ii,1)=sum(Area.*ActiveMask_elements.*GoodDataWeight_elements.*EffectivePressure_elements(:,ii))/totalstatsarea(1);
    rmsn(ii,1)=sqrt(sum(Area.*ActiveMask_elements.*GoodDataWeight_elements.*((EffectivePressure_elements(:,ii)-meann(ii,1)).^2))/totalstatsarea(1));
    meann(ii,2)=sum(Area.*ActiveMask_elements.*FastFlowWeight_elements.*GoodDataWeight_elements.*EffectivePressure_elements(:,ii))/totalstatsarea(2);
    rmsn(ii,2)=sqrt(sum(Area.*ActiveMask_elements.*FastFlowWeight_elements.*GoodDataWeight_elements.*((EffectivePressure_elements(:,ii)-meann(ii,2)).^2))/totalstatsarea(2));
    % lnN:
    meanlnn(ii,1)=sum(Area.*ActiveMask_elements.*GoodDataWeight_elements.*LnN_elements(:,ii))/totalstatsarea(1);
    rmslnn(ii,1)=sqrt(sum(Area.*ActiveMask_elements.*GoodDataWeight_elements.*((LnN_elements(:,ii)-meanlnn(ii,1)).^2))/totalstatsarea(1));
    meanlnn(ii,2)=sum(Area.*ActiveMask_elements.*FastFlowWeight_elements.*GoodDataWeight_elements.*LnN_elements(:,ii))/totalstatsarea(2);
    rmslnn(ii,2)=sqrt(sum(Area.*ActiveMask_elements.*FastFlowWeight_elements.*GoodDataWeight_elements.*((LnN_elements(:,ii)-meanlnn(ii,2)).^2))/totalstatsarea(2));
    % linear tau:
    meantau(ii,1)=sum(Area.*ActiveMask_elements.*GoodDataWeight_elements.*Drag_elements(:,ii))/totalstatsarea(1);
    rmstau(ii,1)=sqrt(sum(Area.*ActiveMask_elements.*GoodDataWeight_elements.*((Drag_elements(:,ii)-meantau(ii,1)).^2))/totalstatsarea(1));
    meantau(ii,2)=sum(Area.*ActiveMask_elements.*FastFlowWeight_elements.*GoodDataWeight_elements.*Drag_elements(:,ii))/totalstatsarea(2);
    rmstau(ii,2)=sqrt(sum(Area.*ActiveMask_elements.*FastFlowWeight_elements.*GoodDataWeight_elements.*((Drag_elements(:,ii)-meantau(ii,2)).^2))/totalstatsarea(2));
    % ln tau:
    meanlntau(ii,1)=sum(Area.*ActiveMask_elements.*GoodDataWeight_elements.*LnTau_elements(:,ii))/totalstatsarea(1);
    rmslntau(ii,1)=sqrt(sum(Area.*ActiveMask_elements.*GoodDataWeight_elements.*((LnTau_elements(:,ii)-meanlntau(ii,1)).^2))/totalstatsarea(1));
    meanlntau(ii,2)=sum(Area.*ActiveMask_elements.*FastFlowWeight_elements.*GoodDataWeight_elements.*LnTau_elements(:,ii))/totalstatsarea(2);
    rmslntau(ii,2)=sqrt(sum(Area.*ActiveMask_elements.*FastFlowWeight_elements.*GoodDataWeight_elements.*((LnTau_elements(:,ii)-meanlntau(ii,2)).^2))/totalstatsarea(2));
    % ln U:
    meanlnu(ii,1)=sum(Area.*ActiveMask_elements.*GoodDataWeight_elements.*LnU_elements(:,ii))/totalstatsarea(1);
    rmslnu(ii,1)=sqrt(sum(Area.*ActiveMask_elements.*GoodDataWeight_elements.*((LnU_elements(:,ii)-meanlnu(ii,1)).^2))/totalstatsarea(1));
    meanlnu(ii,2)=sum(Area.*ActiveMask_elements.*FastFlowWeight_elements.*GoodDataWeight_elements.*LnU_elements(:,ii))/totalstatsarea(2);
    rmslnu(ii,2)=sqrt(sum(Area.*ActiveMask_elements.*FastFlowWeight_elements.*GoodDataWeight_elements.*((LnU_elements(:,ii)-meanlnu(ii,2)).^2))/totalstatsarea(2));
    
    % Compute correlations:
    % R(N,Cw):
    if numlcurves>4 && ii>=4
        rnc(ii,1)=sum(Area.*ActiveMask_elements.*GoodDataWeight_elements.*(EffectivePressure_elements(:,ii)-meann(ii,1)).*((DragCoeff_elements(:,ii-3).^2)-meanc(ii-3,1)))/(totalstatsarea(1)*rmsn(ii,1)*rmsc(ii-3,1));
        rnc(ii,2)=sum(Area.*ActiveMask_elements.*FastFlowWeight_elements.*GoodDataWeight_elements.*(EffectivePressure_elements(:,ii)-meann(ii,2)).*((DragCoeff_elements(:,ii-3).^2)-meanc(ii-3,2)))/(totalstatsarea(2)*rmsn(ii,2)*rmsc(ii-3,2));
    else
        rnc(ii,1)=sum(Area.*ActiveMask_elements.*GoodDataWeight_elements.*(EffectivePressure_elements(:,ii)-meann(ii,1)).*((DragCoeff_elements(:,1).^2)-meanc(1,1)))/(totalstatsarea(1)*rmsn(ii,1)*rmsc(1,1));
        rnc(ii,2)=sum(Area.*ActiveMask_elements.*FastFlowWeight_elements.*GoodDataWeight_elements.*(EffectivePressure_elements(:,ii)-meann(ii,2)).*((DragCoeff_elements(:,1).^2)-meanc(1,2)))/(totalstatsarea(2)*rmsn(ii,2)*rmsc(1,2));
    end
    % R(ln(N),ln(Cw)):
    rlnnlnc(ii,1)=sum(Area.*ActiveMask_elements.*GoodDataWeight_elements.*(LnN_elements(:,ii)-meanlnn(ii,1)).*(LnC_elements(:,1)-meanlnc(1,1)))/(totalstatsarea(1)*rmslnn(ii,1)*rmslnc(1,1));
    rlnnlnc(ii,2)=sum(Area.*ActiveMask_elements.*FastFlowWeight_elements.*GoodDataWeight_elements.*(LnN_elements(:,ii)-meanlnn(ii,2)).*(LnC_elements(:,1)-meanlnc(1,2)))/(totalstatsarea(2)*rmslnn(ii,2)*rmslnc(1,2));
    % R(ln(C),ln(Cw)):
    rlnclnc(ii,1)=sum(Area.*ActiveMask_elements.*GoodDataWeight_elements.*(LnC_elements(:,ii)-meanlnc(ii,1)).*(LnC_elements(:,1)-meanlnc(1,1)))/(totalstatsarea(1)*rmslnc(ii,1)*rmslnc(1,1));
    rlnclnc(ii,2)=sum(Area.*ActiveMask_elements.*FastFlowWeight_elements.*GoodDataWeight_elements.*(LnC_elements(:,ii)-meanlnc(ii,2)).*(LnC_elements(:,1)-meanlnc(1,2)))/(totalstatsarea(2)*rmslnc(ii,2)*rmslnc(1,2));
    % R(tau,tau_w):
    rtautau(ii,1)=sum(Area.*ActiveMask_elements.*GoodDataWeight_elements.*(Drag_elements(:,ii)-meantau(ii,1)).*(Drag_elements(:,1)-meantau(1,1)))/(totalstatsarea(1)*rmstau(ii,1)*rmstau(1,1));
    rtautau(ii,2)=sum(Area.*ActiveMask_elements.*FastFlowWeight_elements.*GoodDataWeight_elements.*(Drag_elements(:,ii)-meantau(ii,2)).*(Drag_elements(:,1)-meantau(1,2)))/(totalstatsarea(2)*rmstau(ii,2)*rmstau(1,2));
    % R(N,tau_w):
    rntau(ii,1)=sum(Area.*ActiveMask_elements.*GoodDataWeight_elements.*(EffectivePressure_elements(:,ii)-meann(ii,1)).*(Drag_elements(:,1)-meantau(1,1)))/(totalstatsarea(1)*rmsn(ii,1)*rmstau(1,1));
    rntau(ii,2)=sum(Area.*ActiveMask_elements.*FastFlowWeight_elements.*GoodDataWeight_elements.*(EffectivePressure_elements(:,ii)-meann(ii,2)).*(Drag_elements(:,1)-meantau(1,2)))/(totalstatsarea(2)*rmsn(ii,2)*rmstau(1,2));
    
    % Compute RMS differences:
    rmslnclnc(ii,1)=sqrt(sum(Area.*ActiveMask_elements.*GoodDataWeight_elements.*((LnC_elements(:,ii)-meanlnc(ii,1)-(LnC_elements(:,1)-meanlnc(1,1))).^2))/totalstatsarea(1));
    rmslnclnc(ii,2)=sqrt(sum(Area.*ActiveMask_elements.*FastFlowWeight_elements.*GoodDataWeight_elements.*((LnC_elements(:,ii)-meanlnc(ii,2)-(LnC_elements(:,1)-meanlnc(1,2))).^2))/totalstatsarea(2));
    rmstautau(ii,1)=sqrt(sum(Area.*ActiveMask_elements.*GoodDataWeight_elements.*((Drag_elements(:,ii)-Drag_elements(:,1)).^2))/totalstatsarea(1));
    rmstautau(ii,2)=sqrt(sum(Area.*ActiveMask_elements.*FastFlowWeight_elements.*GoodDataWeight_elements.*((Drag_elements(:,ii)-Drag_elements(:,1)).^2))/totalstatsarea(2));
    
    % Compute RMS vel misfit:
    rmsvelmisfit(ii,1)=sqrt(sum(Area.*ActiveMask_elements.*GoodDataWeight_elements.*(VelMagMisfit_elements(:,ii).^2))/totalstatsarea(1));
    rmsvelmisfit(ii,2)=sqrt(sum(Area.*ActiveMask_elements.*FastFlowWeight_elements.*GoodDataWeight_elements.*(VelMagMisfit_elements(:,ii).^2))/totalstatsarea(2));
    
    % Compute cost terms:
    jobs(ii,1)=sum(Area.*ObsMisfit_elements(:,ii));
    jobs(ii,2)=sum(Area.*FastFlowWeight_elements.*ObsMisfit_elements(:,ii))*(totalstatsarea(1)/totalstatsarea(2));
    jreg(ii,1)=sum(Area.*RegMisfit_elements(:,ii));
    jreg(ii,2)=sum(Area.*FastFlowWeight_elements.*RegMisfit_elements(:,ii))*(totalstatsarea(1)/totalstatsarea(2));
    
    % Compute velocity equivalent of obs misfit:
    jobsvelequiv(ii,:)=sqrt(jobs(ii,:))*VelScale(ii);
    
end

% Compute squared quantities:
varlnc=rmslnc.^2;
varlnn=rmslnn.^2;
varlntau=rmslntau.^2;
varlnu=rmslnu.^2;
r2nc=rnc.^2;
r2lnnlnc=rlnnlnc.^2;
r2lnclnc=rlnclnc.^2;
r2tautau=rtautau.^2;
r2ntau=rntau.^2;

% compute variance ratios:
varratiolnclnc=varlnc./repmat(varlnc(1,:),[numlcurves,1]);
varratiolnclntau=varlnc./varlntau;
varratiolnclnu=varlnc./varlnu;
varratiojobsjobs=jobs./repmat(jobs(1,:),[numlcurves,1]);
varratiototal=varratiolnclnc.*varratiojobsjobs;

%% Communicate Stats:


% Display stats, whole domain:
string1='RMS(ln(C))=';
string2='var(ln(C))=';
string3='Variance ratio (ln(C)/ln(C1))=';
string4='R(N,C1)=';
string5='R2(N,C1)=';
string6='R(ln(N),ln(C1))=';
string7='R2(ln(N),ln(C1))=';
string8='R(ln(C),ln(C1))=';
string9='R2(ln(C),ln(C1))=';
string10='R(tau,tau1)=';
string11='R2(tau,tau1)=';
string12='R(N,tau1)=';
string13='R2(N,tau1)=';
string14='RMS(ln(C)-ln(C1))=';
string15='RMS(tau-tau1)=';
string16='RMS(u_model-u_obs)=';
string17='Jobs=';
string18='Jreg=';
string19='RMS(ln(N))=';
string20='var(ln(N))=';
string21='Jobs, equiv du=';
string22='RMS(ln(tau))=';
string23='var(ln(tau))=';
string24='RMS(ln(Ub))=';
string25='var(ln(Ub))=';
string26='Variance ratio (ln(C)/ln(tau))=';
string27='Variance ratio (ln(C)/ln(u))=';
string28='Variance ratio (Jobs/Jobs(1))=';
string29='Variance ratio total (Jobsratio*lnCratio)=';
for ii=1:numlcurves
    string1=[string1,num2str(.01*round(100*rmslnc(ii,1))),', '];
    string2=[string2,num2str(.01*round(100*varlnc(ii,1))),', '];
    string3=[string3,num2str(.01*round(100*varratiolnclnc(ii,1))),', '];
    string4=[string4,num2str(.01*round(100*rnc(ii,1))),', '];
    string5=[string5,num2str(.01*round(100*r2nc(ii,1))),', '];
    string6=[string6,num2str(.01*round(100*rlnnlnc(ii,1))),', '];
    string7=[string7,num2str(.01*round(100*r2lnnlnc(ii,1))),', '];
    string8=[string8,num2str(.01*round(100*rlnclnc(ii,1))),', '];
    string9=[string9,num2str(.01*round(100*r2lnclnc(ii,1))),', '];
    string10=[string10,num2str(.01*round(100*rtautau(ii,1))),', '];
    string11=[string11,num2str(.01*round(100*r2tautau(ii,1))),', '];
    string12=[string12,num2str(.01*round(100*rntau(ii,1))),', '];
    string13=[string13,num2str(.01*round(100*r2ntau(ii,1))),', '];
    string14=[string14,num2str(.01*round(100*rmslnclnc(ii,1))),', '];
    string15=[string15,num2str(.01*round(100*rmstautau(ii,1)/1000)),', '];
    string16=[string16,num2str(.01*round(100*rmsvelmisfit(ii,1))),', '];
    string17=[string17,num2str(.01*round(100*jobs(ii,1)*10000)),', '];
    string18=[string18,num2str(.01*round(100*jreg(ii,1)*10000)),', '];
    string19=[string19,num2str(.01*round(100*rmslnn(ii,1))),', '];
    string20=[string20,num2str(.01*round(100*varlnn(ii,1))),', '];
    string21=[string21,num2str(.01*round(100*jobsvelequiv(ii,1))),', '];
    string22=[string22,num2str(.01*round(100*rmslntau(ii,1))),', '];
    string23=[string23,num2str(.01*round(100*varlntau(ii,1))),', '];
    string24=[string24,num2str(.01*round(100*rmslnu(ii,1))),', '];
    string25=[string25,num2str(.01*round(100*varlnu(ii,1))),', '];
    string26=[string26,num2str(.01*round(100*varratiolnclntau(ii,1))),', '];
    string27=[string27,num2str(.01*round(100*varratiolnclnu(ii,1))),', '];
    string28=[string28,num2str(.01*round(100*varratiojobsjobs(ii,1))),', '];
    string29=[string29,num2str(.01*round(100*varratiototal(ii,1))),', '];
end
disp('...')
disp('Whole Domain Stats:')
disp(string1(1:end-1))
disp(string2(1:end-1))
disp(string3(1:end-1))
disp(string4(1:end-1))
disp(string5(1:end-1))
disp(string6(1:end-1))
disp(string7(1:end-1))
disp(string8(1:end-1))
disp(string9(1:end-1))
disp(string10(1:end-1))
disp(string11(1:end-1))
disp(string12(1:end-1))
disp(string13(1:end-1))
disp(string14(1:end-1))
disp([string15(1:end-1),' kPa'])
disp([string16(1:end-1),' m/yr'])
disp([string17(1:end-1),' x10^-4'])
disp([string18(1:end-1),' x10^-4'])
disp(string19(1:end-1))
disp(string20(1:end-1))
disp([string21(1:end-1),' m/yr'])
disp(string22(1:end-1))
disp(string23(1:end-1))
disp(string24(1:end-1))
disp(string25(1:end-1))
disp(string26(1:end-1))
disp(string27(1:end-1))
disp(string28(1:end-1))
disp(string29(1:end-1))

% Display stats, fast flow:
string1='RMS(ln(C))=';
string2='var(ln(C))=';
string3='Variance ratio (ln(C)/ln(C1))=';
string4='R(N,C1)=';
string5='R2(N,C1)=';
string6='R(ln(N),ln(C1))=';
string7='R2(ln(N),ln(C1))=';
string8='R(ln(C),ln(C1))=';
string9='R2(ln(C),ln(C1))=';
string10='R(tau,tau1)=';
string11='R2(tau,tau1)=';
string12='R(N,tau1)=';
string13='R2(N,tau1)=';
string14='RMS(ln(C)-ln(C1))=';
string15='RMS(tau-tau1)=';
string16='RMS(u_model-u_obs)=';
string17='Jobs=';
string18='Jreg=';
string19='RMS(ln(N))=';
string20='var(ln(N))=';
string21='Jobs, equiv du=';
string22='RMS(ln(tau))=';
string23='var(ln(tau))=';
string24='RMS(ln(Ub))=';
string25='var(ln(Ub))=';
string26='Variance ratio (ln(C)/ln(tau))=';
string27='Variance ratio (ln(C)/ln(u))=';
string28='Variance ratio (Jobs/Jobs(1))=';
string29='Variance ratio total (Jobsratio*lnCratio)=';
for ii=1:numlcurves
    string1=[string1,num2str(.01*round(100*rmslnc(ii,2))),', '];
    string2=[string2,num2str(.01*round(100*varlnc(ii,2))),', '];
    string3=[string3,num2str(.01*round(100*varratiolnclnc(ii,2))),', '];
    string4=[string4,num2str(.01*round(100*rnc(ii,2))),', '];
    string5=[string5,num2str(.01*round(100*r2nc(ii,2))),', '];
    string6=[string6,num2str(.01*round(100*rlnnlnc(ii,2))),', '];
    string7=[string7,num2str(.01*round(100*r2lnnlnc(ii,2))),', '];
    string8=[string8,num2str(.01*round(100*rlnclnc(ii,2))),', '];
    string9=[string9,num2str(.01*round(100*r2lnclnc(ii,2))),', '];
    string10=[string10,num2str(.01*round(100*rtautau(ii,2))),', '];
    string11=[string11,num2str(.01*round(100*r2tautau(ii,2))),', '];
    string12=[string12,num2str(.01*round(100*rntau(ii,2))),', '];
    string13=[string13,num2str(.01*round(100*r2ntau(ii,2))),', '];
    string14=[string14,num2str(.01*round(100*rmslnclnc(ii,2))),', '];
    string15=[string15,num2str(.01*round(100*rmstautau(ii,2)/1000)),', '];
    string16=[string16,num2str(.01*round(100*rmsvelmisfit(ii,2))),', '];
    string17=[string17,num2str(.01*round(100*jobs(ii,2)*10000)),', '];
    string18=[string18,num2str(.01*round(100*jreg(ii,2)*10000)),', '];
    string19=[string19,num2str(.01*round(100*rmslnn(ii,2))),', '];
    string20=[string20,num2str(.01*round(100*varlnn(ii,2))),', '];
    string21=[string21,num2str(.01*round(100*jobsvelequiv(ii,2))),', '];
    string22=[string22,num2str(.01*round(100*rmslntau(ii,2))),', '];
    string23=[string23,num2str(.01*round(100*varlntau(ii,2))),', '];
    string24=[string24,num2str(.01*round(100*rmslnu(ii,2))),', '];
    string25=[string25,num2str(.01*round(100*varlnu(ii,2))),', '];
    string26=[string26,num2str(.01*round(100*varratiolnclntau(ii,2))),', '];
    string27=[string27,num2str(.01*round(100*varratiolnclnu(ii,2))),', '];
    string28=[string28,num2str(.01*round(100*varratiojobsjobs(ii,2))),', '];
    string29=[string29,num2str(.01*round(100*varratiototal(ii,2))),', '];
end
disp('...')
disp('Fast-Flow Stats:')
disp(string1(1:end-1))
disp(string2(1:end-1))
disp(string3(1:end-1))
disp(string4(1:end-1))
disp(string5(1:end-1))
disp(string6(1:end-1))
disp(string7(1:end-1))
disp(string8(1:end-1))
disp(string9(1:end-1))
disp(string10(1:end-1))
disp(string11(1:end-1))
disp(string12(1:end-1))
disp(string13(1:end-1))
disp(string14(1:end-1))
disp([string15(1:end-1),' kPa'])
disp([string16(1:end-1),' m/yr'])
disp([string17(1:end-1),' x10^-4'])
disp([string18(1:end-1),' x10^-4'])
disp(string19(1:end-1))
disp(string20(1:end-1))
disp([string21(1:end-1),' m/yr'])
disp(string22(1:end-1))
disp(string23(1:end-1))
disp(string24(1:end-1))
disp(string25(1:end-1))
disp(string26(1:end-1))
disp(string27(1:end-1))
disp(string28(1:end-1))
disp(string29(1:end-1))
disp('...')

%toc 
%return


%% Make L-Curves Figure:

if makelcurvefig
    
    % Make the figure:
    figure
    
    % Plot cost vs cost, log scale:
    % Make subplot:
    subplot('Position',Boxes_lcurves{1})
    % pre-allocate handles:
    thesehandles=zeros(numlcurves,1);
    % loop through L-curves:
    for thismodel=1:numlcurves
        % Plot model costs:
        plot(exp(LCurveData(thismodel).LogRegCost_points(~LCurveData(thismodel).IsOutlier)),exp(LCurveData(thismodel).LogDataCost_points(~LCurveData(thismodel).IsOutlier)),...
            'Marker',modelmarker,'MarkerFaceColor',lcurvecolors{thismodel},'Color',lcurvecolors{thismodel},'LineStyle','none','MarkerSize',modelmarkersize)
        hold on
        plot(exp(LCurveData(thismodel).LogRegCost_points(LCurveData(thismodel).IsOutlier)),exp(LCurveData(thismodel).LogDataCost_points(LCurveData(thismodel).IsOutlier)),...
            'Marker',modelmarker,'MarkerFaceColor','none','Color',lcurvecolors{thismodel},'LineStyle','none','MarkerSize',modelmarkersize)
        % Plot smooth costs:
        plot(exp(LCurveData(thismodel).LogRegCost_curve),exp(LCurveData(thismodel).LogDataCost_curve),...
            'Marker','none','Color',lcurvecolors{thismodel},'LineStyle',curvelinestyle,'LineWidth',curvelinesize);
        % Plot corner region:
        thesehandles(thismodel)=plot(exp(LCurveData(thismodel).LogRegCost_curve(LCurveData(thismodel).bestlambdaind_minallowable:LCurveData(thismodel).bestlambdaind_maxallowable)),...
            exp(LCurveData(thismodel).LogDataCost_curve(LCurveData(thismodel).bestlambdaind_minallowable:LCurveData(thismodel).bestlambdaind_maxallowable)),...
            'Marker','none','Color',lcurvecolors{thismodel},'LineStyle',cornerlinestyle,'LineWidth',cornerlinesize);
        % Plot best lambda value:
        plot(exp(LCurveData(thismodel).LogRegCost_curve(LCurveData(thismodel).bestlambdaind)),exp(LCurveData(thismodel).LogDataCost_curve(LCurveData(thismodel).bestlambdaind)),...
            'Marker',bestmarker,'MarkerFaceColor',lcurvecolors{thismodel},'Color',lcurvecolors{thismodel},'LineStyle','none','MarkerSize',bestmarkersize)
    end
    % Note outliers:
    %     for thismodel=1:numlcurves
    %         if sum(LCurveData(thismodel).IsOutlier)>0
    %             numoutliers=sum(LCurveData(thismodel).IsOutlier);
    %             minoutlierjobs=exp(min(LCurveData(thismodel).LogDataCost_points(LCurveData(thismodel).IsOutlier)));
    %             minoutlierpow10=floor(log10(minoutlierjobs));
    %             text(jreglims(2)*.9,jobslims(2),['\uparrow ',num2str(sum(LCurveData(thismodel).IsOutlier)),' outliers for ',descriptext{thismodel},', with J_{obs}>10^{',num2str(minoutlierpow10),'} \uparrow'],...
    %                 'HorizontalAlignment','right','VerticalAlignment','top','Color',lcurvecolors{thismodel})
    %         end
    %     end
    % Set up axes etc:
    xlim(jreglims.*[.999,1.001])
    ylim(jobslims.*[.999,1.001])
    set(gca,'XScale','log')
    set(gca,'YScale','log')
    set(gca,'TickDir','out')
    set(gca,'PlotBoxAspectRatio',[(log(jreglims(2))-log(jreglims(1)))/(log(jobslims(2))-log(jobslims(1))),1,1])
    xlabel('Regularization cost, J_{reg} (unitless)')
    ylabel('Data cost, J_{obs} (unitless)')
    legend(thesehandles,descriptext,'location','SouthWest')
    title('a) L-Curve Plots (cost vs cost)')
    
    
    % Loop through L-curves:
    for thismodel=1:numlcurves
        
        % Compute lambda subsampling:
        LogLambda_path=linspace(LCurveData(thismodel).LogLambda_points(find(~LCurveData(thismodel).IsOutlier,1,'first')),...
            LCurveData(thismodel).LogLambda_points(find(~LCurveData(thismodel).IsOutlier,1,'last')),length(LCurveData(thismodel).LogDataCost_curve))';
        
        % Plot cost components:
        % Make subplot:
        subplot('Position',Boxes_lcurves{2*thismodel})
        % pre-allocate handles:
        thesehandles=zeros(4,1);
        % Plot good model points:
        thesehandles(1)=plot(exp(LCurveData(thismodel).LogLambda_points(~LCurveData(thismodel).IsOutlier)),exp(LCurveData(thismodel).LogDataCost_points(~LCurveData(thismodel).IsOutlier)),'Marker','o','MarkerFaceColor','b','Color','b','LineStyle','none');
        hold on
        thesehandles(2)=plot(exp(LCurveData(thismodel).LogLambda_points(~LCurveData(thismodel).IsOutlier)),exp(LCurveData(thismodel).LogRegCost_points(~LCurveData(thismodel).IsOutlier)),'Marker','o','MarkerFaceColor','r','Color','r','LineStyle','none');
        % Plot bad model points:
        if sum(LCurveData(thismodel).IsOutlier)>0
            plot(exp(LCurveData(thismodel).LogLambda_points(LCurveData(thismodel).IsOutlier)),exp(LCurveData(thismodel).LogDataCost_points(LCurveData(thismodel).IsOutlier)),'Marker','o','MarkerFaceColor','none','Color','b','LineStyle','none')
            plot(exp(LCurveData(thismodel).LogLambda_points(LCurveData(thismodel).IsOutlier)),exp(LCurveData(thismodel).LogRegCost_points(LCurveData(thismodel).IsOutlier)),'Marker','o','MarkerFaceColor','none','Color','r','LineStyle','none')
        end
        % Plot smooth curve:
        thesehandles(3)=plot(exp(LogLambda_path),exp(LCurveData(thismodel).LogDataCost_curve),'b','LineWidth',1);
        thesehandles(4)=plot(exp(LogLambda_path),exp(LCurveData(thismodel).LogRegCost_curve),'r','LineWidth',1);
        % Compute axes limits:
        thesexlims=10.^[floor(log10(lambdamin)),ceil(log10(lambdamax))];
        theseylims=[min(jreglims(1),jobslims(1)),max(jreglims(2),jobslims(2))];
        % Set up axes etc:
        set(gca,'XScale','log')
        set(gca,'YScale','log')
        xlim(thesexlims.*[.999,1.001])
        ylim(theseylims.*[.999,1.001])
        set(gca,'XTick',10.^(log10(thesexlims(1)):log10(thesexlims(2))))
        if thismodel==numlcurves
            xlabel('\lambda value (unitless)')
        else
            set(gca,'XTickLabel',[])
        end
        ylabel('Cost Components J (unitless)')
        legend(thesehandles,'J_{obs} (individual models)','J_{reg} (individual models)','J_{obs} (smooth curve)','J_{reg} (smooth curve)','location','SouthWest')
        %         if sum(LCurveData(thismodel).IsOutlier)>0
        %             legend('J_{obs} (individual models)','J_{reg} (individual models)','J_{obs} (outlier models)','J_{reg} (outlier models)','J_{obs} (smooth curve)','J_{reg} (smooth curve)','location','SouthWest')
        %         else
        %             legend('J_{obs} (individual models)','J_{reg} (individual models)','J_{obs} (smooth curve)','J_{reg} (smooth curve)','location','SouthWest')
        %         end
        title([char(96+2*thismodel),') Cost components for ',descriptext{thismodel}])
        
        % Plot curvatures:
        % Make subplot:
        subplot('Position',Boxes_lcurves{2*thismodel+1})
        % Identify max curvature:
        maxcurve=LCurveData(thismodel).TotalCurvature_curve(LCurveData(thismodel).bestlambdaind-1);
        % Plot best curve curvatures:
        plot(exp(LogLambda_path(2:end-1)),LCurveData(thismodel).LogDataCurvature_curve,'Color',datacurvecolor,'LineWidth',1)
        hold on
        plot(exp(LogLambda_path(2:end-1)),LCurveData(thismodel).LogRegCurvature_curve,'Color',regcurvecolor,'LineWidth',1)
        plot(exp(LogLambda_path(2:end-1)),LCurveData(thismodel).TotalCurvature_curve,'Color',totcurvecolor,'LineWidth',1)
        % Show best value and uncertainty:
        plot(BestLambda(thismodel),maxcurve,'Color','k','MarkerFaceColor','k','Marker','o')
        text(BestLambda(thismodel),maxcurve*1.025,{'Best';'value'},'HorizontalAlignment','center','VerticalAlignment','bottom')
        plot(BestLambda(thismodel)*[1/DeltaLambdaRatio(thismodel),DeltaLambdaRatio(thismodel)],maxcurve*[1,1],'k')
        plot((BestLambda(thismodel)/DeltaLambdaRatio(thismodel))*[1,1],maxcurve*[.975,1.025],'k')
        plot(BestLambda(thismodel)*DeltaLambdaRatio(thismodel)*[1,1],maxcurve*[.975,1.025],'k')
        % Show max/min lambda and uncertainties:
        %plot(BestLambda(thismodel)*[1,1],maxcurve*[.5,1],'--k')
        %plot([LCurveData(thismodel).bestlambda_minallowable,LCurveData(thismodel).bestlambda_maxallowable],.5*maxcurve*[1,1],'--k')
        plot(MinLambda(thismodel),.5*maxcurve,'Color','k','MarkerFaceColor','k','Marker','o')
        plot(MinLambda(thismodel)*[1/DeltaLambdaRatio(thismodel),DeltaLambdaRatio(thismodel)],.5*maxcurve*[1,1],'k')
        plot((MinLambda(thismodel)/DeltaLambdaRatio(thismodel))*[1,1],.5*maxcurve*[.975,1.025],'k')
        plot(MinLambda(thismodel)*DeltaLambdaRatio(thismodel)*[1,1],.5*maxcurve*[.975,1.025],'k')
        plot(MinLambda(thismodel)*[1,.5],.5*maxcurve+[0,.01],'k')
        text(MinLambda(thismodel)*.5,.5*maxcurve+.01,{'Min';'value'},'HorizontalAlignment','right','VerticalAlignment','bottom','Color','k')
        plot(MaxLambda(thismodel),.5*maxcurve,'Color','k','MarkerFaceColor','k','Marker','o')
        plot(MaxLambda(thismodel)*[1/DeltaLambdaRatio(thismodel),DeltaLambdaRatio(thismodel)],.5*maxcurve*[1,1],'k')
        plot((MaxLambda(thismodel)/DeltaLambdaRatio(thismodel))*[1,1],.5*maxcurve*[.975,1.025],'k')
        plot(MaxLambda(thismodel)*DeltaLambdaRatio(thismodel)*[1,1],.5*maxcurve*[.975,1.025],'k')
        plot(MaxLambda(thismodel)*[1,2],.5*maxcurve+[0,.01],'k')
        text(MaxLambda(thismodel)*2,.5*maxcurve+.01,{'Max';'value'},'HorizontalAlignment','left','VerticalAlignment','bottom','Color','k')
        % Determine axes limits:
        thesexlims=10.^[floor(log10(exp(LCurveData(thismodel).LogLambda_points(1)))),ceil(log10(exp(LCurveData(thismodel).LogLambda_points(end))))];
        theseylims=maxcurve*[-.25,1.25];
        % Plot zero line:
        plot(thesexlims,[0,0],'--k')
        % Set up axes etc:
        xlim(thesexlims.*[.999,1.001])
        ylim(theseylims)
        set(gca,'XScale','log')
        set(gca,'XTick',10.^(log10(thesexlims(1)):log10(thesexlims(2))))
        if thismodel==numlcurves
            xlabel('\lambda value (unitless)')
        else
            set(gca,'XTickLabel',[])
        end
        ylabel('Cost curvature, $\frac{d^{2}(ln(J))}{d(ln(\lambda))^{2}}$','interpreter','latex')
        title([char(96+2*thismodel+1),') Cost curvature for ',descriptext{thismodel}])
        if log(BestLambda(thismodel))>mean([LCurveData(thismodel).LogLambda_points(1),LCurveData(thismodel).LogLambda_points(end)])
            locationstring='NorthWest';
        else
            locationstring='NorthEast';
        end
        legend('J_{obs}','J_{reg}','Total','location',locationstring)
    end
    
    
    % Save figure:
    set(gcf,'PaperSize',lcurvepagesize)
    set(gcf,'PaperPosition',[0,0,lcurvepagesize])
    print('-dpng',figname_lcurves,['-r',num2str(resolution)])
    
end
%toc
%return


%% Make Map Figure:

if makemapfig
    
    % Make the figure:
    figure
    
    % Make faded jet:
    MakeFadedJet;
    
    % Loop through models:
    for thismodel=1:numlcurves
        
        % Generate drag coeff units string:
        if ~strcmp(InversionParameters{thismodel}.nsource,'none')
            if InversionParameters{thismodel}.r==1
                if InversionParameters{thismodel}.m==1
                    cunitsstring='(m/a)^{-1}';
                else
                    cunitsstring=['(m/a)^{-1/',num2str(InversionParameters{thismodel}.m),'}'];
                end
            else
                if InversionParameters{thismodel}.m==1
                    cunitsstring=['Pa^{',num2str(1-InversionParameters{thismodel}.r),'}(m/a)^{-1}'];
                else
                    cunitsstring=['Pa^{',num2str(1-InversionParameters{thismodel}.r),'}(m/a)^{-1/',num2str(InversionParameters{thismodel}.m),'}'];
                end
            end
        else
            if InversionParameters{thismodel}.m==1
                cunitsstring='Pa(m/a)^{-1}';
            else
                cunitsstring=['Pa(m/a)^{-1/',num2str(InversionParameters{thismodel}.m),'}'];
            end
        end
        
        % Compute drag coeff ticks and limits:
        % Sort elements by drag coefficient:
        [SortedDragCoeff_elements,SortedInd]=sort(2*log10(DragCoeff_elements(ActiveMask_elements>.99,thismodel)/(secondsperyear^(1/(2*InversionParameters{thismodel}.m)))));
        % Compute cumulative sorted area:
        CumulativeSortedArea=Area(ActiveMask_elements>.99);
        CumulativeSortedArea=cumsum(CumulativeSortedArea(SortedInd));
        % Interpolate percentiles:
        cclims=interp1(CumulativeSortedArea/CumulativeSortedArea(end),SortedDragCoeff_elements,dragcoeffclims/100);
        %         if exist('dragcoeffclims','var')
        %             cclims=log10(dragcoeffclims(thismodel,:));
        %         else
        %             % Compute limits from percentiles:
        %             cclims=2*log10(prctile(DragCoeff_vertices(ActiveMask_drag,thismodel)/(secondsperyear^(1/(2*InversionParameters{thismodel}.m))),[.1,99.9]));
        %             % Round up/down to powers of ten:
        %             cclims=[floor(cclims(1)),ceil(cclims(2))];
        %         end
        % Compute ticks and labels for color bar:
        % Compute first base and tick for color bar:
        thiscbase=10^floor(cclims(1));
        thistick=thiscbase;
        % Create all ticks and labels for color bar:
        ccticks=zeros(100,1);
        cclabels=cell(100,1);
        for ii=1:100
            % Record this tick:
            ccticks(ii)=log10(thistick);
            % Check if we hit the next power of ten:
            if abs(thistick-10*thiscbase)/thiscbase<1e-12 || (ii==1 && abs(thistick-thiscbase)/thiscbase<1e-12)
                % Generate label:
                cclabels{ii}=['10^{',num2str(round(log10(thistick))),'}'];
                % Advance base:
                thiscbase=thistick;
            end
            % Advance tick:
            thistick=thistick+thiscbase;
            % Break from loop and truncate:
            if thistick>(1+1e-12)*(10^cclims(2))
                ccticks=ccticks(1:ii);
                cclabels=cclabels(1:ii);
                break
            end
        end
        
        
        % Plot drag coefficient:
        % Make subplot:
        subplot('Position',Boxes_maps{ceil(thismodel/3),thismodel-3*floor((thismodel-1)/3)})
        % Draw back plane:
        patch('Vertices',[[xlims(1);xlims(2);xlims(2);xlims(1);xlims(1)],...
            [ylims(1);ylims(1);ylims(2);ylims(2);ylims(1)]]/1e3,'Faces',1:5,'EdgeColor','none','FaceColor',backplanecolor)
        hold on
        % Draw main figure:
        patch('Vertices',[x2d/1000,y2d/1000],'Faces',elements2d(IceAreaFraction_elements>.5&OceanAreaFraction_elements<.5,:),'EdgeColor','none','FaceColor','interp','FaceVertexCData',2*log10(DragCoeff_vertices(:,thismodel)/(secondsperyear^(1/(2*InversionParameters{thismodel}.m)))))
        % Draw zoom boxes:
        for ii=1:numzoomboxes
            % Check whether to show zoom box location:
            if showzoomboxloc
                % draw small outline on main map:
                plot([zoomboxes(ii,1:4),zoomboxes(ii,1)]/1e3,[zoomboxes(ii,5:8),zoomboxes(ii,5)]/1e3,'k','LineWidth',1)
            end
            % Check whether to do the zoom box itself:
            if showzoombox
                % drag backbround patch:
                patch('Vertices',[[zoomboxes_big(ii,1);zoomboxes_big(ii,3);zoomboxes_big(ii,3);zoomboxes_big(ii,1);zoomboxes_big(ii,1)],...
                    [zoomboxes_big(ii,2);zoomboxes_big(ii,2);zoomboxes_big(ii,4);zoomboxes_big(ii,4);zoomboxes_big(ii,2)]]/1e3,'Faces',1:5,'EdgeColor','none','FaceColor',backplanecolor)
                % draw zoom patch:
                patch('Vertices',[(zoomboxes_big(ii,1)+zoomratio(ii)*x2d_rot(:,ii))/1000,(zoomboxes_big(ii,2)+zoomratio(ii)*y2d_rot(:,ii))/1000],...
                    'Faces',elements2d(IceAreaFraction_elements>.5&OceanAreaFraction_elements<.5&WithinZoomBoxes_elements(:,ii)>.01,:),'EdgeColor','none','FaceColor','interp','FaceVertexCData',2*log10(DragCoeff_vertices(:,thismodel)/(secondsperyear^(1/(2*InversionParameters{thismodel}.m)))))
                % draw covering patch:
                if ii==1 || ii==4
                    patch('Vertices',[[zoomboxes_big(ii,1);zoomboxes_big(ii,3);zoomboxes_big(ii,3);zoomboxes_big(ii,3)+outlinestripewidth*1e3;zoomboxes_big(ii,3)+outlinestripewidth*1e3;zoomboxes_big(ii,1);zoomboxes_big(ii,1)],...
                        [zoomboxes_big(ii,2);zoomboxes_big(ii,2);zoomboxes_big(ii,4);zoomboxes_big(ii,4);zoomboxes_big(ii,2)-outlinestripewidth*1e3;zoomboxes_big(ii,2)-outlinestripewidth*1e3;zoomboxes_big(ii,2)]]/1e3,'Faces',1:7,'EdgeColor','none','FaceColor',backplanecolor)
                end
                % drag big outline:
                plot([zoomboxes_big(ii,1),zoomboxes_big(ii,3),zoomboxes_big(ii,3),zoomboxes_big(ii,1),zoomboxes_big(ii,1)]/1e3,...
                    [zoomboxes_big(ii,2),zoomboxes_big(ii,2),zoomboxes_big(ii,4),zoomboxes_big(ii,4),zoomboxes_big(ii,2)]/1e3,'k','LineWidth',1)
                % Label:
                if ii==1
                    text(zoomboxes_big(ii,1)/1e3+zoomboxnamebuffer,zoomboxes_big(ii,4)/1e3-zoomboxnamebuffer,zoomboxnames{ii},...
                        'HorizontalAlignment','left','VerticalAlignment','top','Color','k','BackgroundColor','w','FontSize',zoomboxnamefontsize,'Margin',zoomboxnamemargin)
                elseif ii==2
                    text(zoomboxes_big(ii,3)/1e3-zoomboxnamebuffer,zoomboxes_big(ii,2)/1e3-zoomboxnamebuffer,zoomboxnames{ii},...
                        'HorizontalAlignment','right','VerticalAlignment','top','Color','k','BackgroundColor','w','FontSize',zoomboxnamefontsize,'Margin',zoomboxnamemargin)
                elseif ii==3
                    text(zoomboxes_big(ii,3)/1e3-zoomboxnamebuffer,zoomboxes_big(ii,4)/1e3+zoomboxnamebuffer,zoomboxnames{ii},...
                        'HorizontalAlignment','right','VerticalAlignment','bottom','Color','k','BackgroundColor','w','FontSize',zoomboxnamefontsize,'Margin',zoomboxnamemargin)
                else
                    text(zoomboxes_big(ii,3)/1e3-zoomboxnamebuffer,zoomboxes_big(ii,2)/1e3-zoomboxnamebuffer,zoomboxnames{ii},...
                        'HorizontalAlignment','right','VerticalAlignment','top','Color','k','BackgroundColor','w','FontSize',zoomboxnamefontsize,'Margin',zoomboxnamemargin)
                end
                % scale bar:
                patch('Vertices',[zoomboxes_big(ii,3)/1e3-zoomboxscalebarbufferx(ii)-zoomratio(ii)*scalebarlength*[1;0;0;1;1],zoomboxes_big(ii,2)/1e3+zoomboxscalebarbuffery(ii)+scalebarheight*[0;0;1;1;0]],'Faces',1:5,'EdgeColor','k','FaceColor','w')
                text(zoomboxes_big(ii,3)/1e3-zoomboxscalebarbufferx(ii)-.5*zoomratio(ii)*scalebarlength,zoomboxes_big(ii,2)/1e3+zoomboxscalebarbuffery(ii),[num2str(scalebarlength),' km'],...
                    'HorizontalAlignment','center','VerticalAlignment','bottom','FontSize',zoomboxscalebarfontsize)
            end
        end
        % Set up axes:
        xlim(xlims/1000)
        ylim(ylims/1000)
        set(gca,'XTick',xtick*[ceil(xlims(1)/(1e3*xtick)):1:floor(xlims(2)/(1e3*xtick))])
        set(gca,'YTick',ytick*[ceil(ylims(1)/(1e3*ytick)):1:floor(ylims(2)/(1e3*ytick))])
        if rem(thismodel,3)==1
            ylabel('Northing (km)','FontSize',axesfontsize)
        else
            set(gca,'YTickLabel',[])
        end
        if numlcurves-thismodel<3
            xlabel('Easting (km)')
        else
            set(gca,'XTickLabel',[])
        end
        set(gca,'PlotBoxAspectRatio',[(xlims(2)-xlims(1))/(ylims(2)-ylims(1)),1,1])
        set(gca,'TickDir','out')
        set(gca,'FontSize',axesfontsize)
        %set(gca,'PositionConstraint','innerposition')
        % Set up color scale:
        eval(['thiscmap=colormap(gca,',dragcoeffcmap,');'])
        if flipdragcoeffcmap
            colormap(gca,flipud(thiscmap))
        else
            colormap(gca,thiscmap)
        end
        caxis(gca,cclims)
        hc=colorbar;
        set(hc,'YTick',ccticks)
        set(hc,'YTickLabel',cclabels)
        set(get(hc,'Ylabel'),'String',cunitsstring,'FontSize',colorbarfontsize)
        set(hc,'TickDir','out')
        set(hc,'FontSize',colorbarfontsize)
        %thisposition=get(gca,'Position');
        %set(gca,'Position',[thisposition(1),Boxes{1,thismodel}(2),thisposition(3:4)])
        % Make title:
        title([char(96+thismodel),') Drag coefficient (C) for ',descriptext{thismodel}],'FontSize',titlefontsize)
        
        
    end
    
    % Save figure:
    set(gcf,'PaperSize',[pagewidth_map,pageheight_map])
    set(gcf,'PaperPosition',[0,0,pagewidth_map,pageheight_map])
    print('-dpng',figname_maps,['-r',num2str(resolution)])
    
end

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

%% Cutting Room Floor:


% % Plot drag:
% % Make subplot:
% subplot('Position',Boxes_maps{2,thismodel})
% % Draw back plane:
% patch('Vertices',[[xlims(1);xlims(2);xlims(2);xlims(1);xlims(1)],...
%     [ylims(1);ylims(1);ylims(2);ylims(2);ylims(1)]]/1e3,'Faces',1:5,'EdgeColor','none','FaceColor',backplanecolor)
% hold on
% % Draw main figure:
% if thismodel==1
%     patch('Vertices',[x2d/1000,y2d/1000],'Faces',elements2d(IceAreaFraction_elements>.5&OceanAreaFraction_elements<.5,:),'EdgeColor','none','FaceColor','interp','FaceVertexCData',Drag_vertices(:,thismodel)/1e3)
% else
%     patch('Vertices',[x2d/1000,y2d/1000],'Faces',elements2d(IceAreaFraction_elements>.5&OceanAreaFraction_elements<.5,:),'EdgeColor','none','FaceColor','interp','FaceVertexCData',(Drag_vertices(:,thismodel)-Drag_vertices(:,1))/1e3)
% end
% % Draw zoom boxes:
% for ii=1:numzoomboxes
%     % Check whether to show zoom box location:
%     if showzoomboxloc
%         % draw small outline on main map:
%         plot([zoomboxes(ii,1:4),zoomboxes(ii,1)]/1e3,[zoomboxes(ii,5:8),zoomboxes(ii,5)]/1e3,'k','LineWidth',1)
%     end
%     % Check whether to do the zoom box itself:
%     if showzoombox
%         % drag backbround patch:
%         patch('Vertices',[[zoomboxes_big(ii,1);zoomboxes_big(ii,3);zoomboxes_big(ii,3);zoomboxes_big(ii,1);zoomboxes_big(ii,1)],...
%             [zoomboxes_big(ii,2);zoomboxes_big(ii,2);zoomboxes_big(ii,4);zoomboxes_big(ii,4);zoomboxes_big(ii,2)]]/1e3,'Faces',1:5,'EdgeColor','none','FaceColor',backplanecolor)
%         % draw zoom patch:
%         if thismodel==1
%             patch('Vertices',[(zoomboxes_big(ii,1)+zoomratio(ii)*x2d_rot(:,ii))/1000,(zoomboxes_big(ii,2)+zoomratio(ii)*y2d_rot(:,ii))/1000],...
%                 'Faces',elements2d(IceAreaFraction_elements>.5&OceanAreaFraction_elements<.5&WithinZoomBoxes_elements(:,ii)>.01,:),'EdgeColor','none','FaceColor','interp','FaceVertexCData',Drag_vertices(:,thismodel)/1e3)
%         else
%             patch('Vertices',[(zoomboxes_big(ii,1)+zoomratio(ii)*x2d_rot(:,ii))/1000,(zoomboxes_big(ii,2)+zoomratio(ii)*y2d_rot(:,ii))/1000],...
%                 'Faces',elements2d(IceAreaFraction_elements>.5&OceanAreaFraction_elements<.5&WithinZoomBoxes_elements(:,ii)>.01,:),'EdgeColor','none','FaceColor','interp','FaceVertexCData',(Drag_vertices(:,thismodel)-Drag_vertices(:,1))/1e3)
%         end
%         % draw covering patch:
%         if ii==1 || ii==4
%             patch('Vertices',[[zoomboxes_big(ii,1);zoomboxes_big(ii,3);zoomboxes_big(ii,3);zoomboxes_big(ii,3)+outlinestripewidth*1e3;zoomboxes_big(ii,3)+outlinestripewidth*1e3;zoomboxes_big(ii,1);zoomboxes_big(ii,1)],...
%                 [zoomboxes_big(ii,2);zoomboxes_big(ii,2);zoomboxes_big(ii,4);zoomboxes_big(ii,4);zoomboxes_big(ii,2)-outlinestripewidth*1e3;zoomboxes_big(ii,2)-outlinestripewidth*1e3;zoomboxes_big(ii,2)]]/1e3,'Faces',1:7,'EdgeColor','none','FaceColor',backplanecolor)
%         end
%         % drag big outline:
%         plot([zoomboxes_big(ii,1),zoomboxes_big(ii,3),zoomboxes_big(ii,3),zoomboxes_big(ii,1),zoomboxes_big(ii,1)]/1e3,...
%             [zoomboxes_big(ii,2),zoomboxes_big(ii,2),zoomboxes_big(ii,4),zoomboxes_big(ii,4),zoomboxes_big(ii,2)]/1e3,'k','LineWidth',1)
%         % label:
%         if ii==1
%             text(zoomboxes_big(ii,1)/1e3+zoomboxnamebuffer,zoomboxes_big(ii,4)/1e3-zoomboxnamebuffer,zoomboxnames{ii},...
%                 'HorizontalAlignment','left','VerticalAlignment','top','Color','k','BackgroundColor','w','FontSize',zoomboxnamefontsize,'Margin',zoomboxnamemargin)
%         elseif ii==2
%             text(zoomboxes_big(ii,3)/1e3-zoomboxnamebuffer,zoomboxes_big(ii,2)/1e3-zoomboxnamebuffer,zoomboxnames{ii},...
%                 'HorizontalAlignment','right','VerticalAlignment','top','Color','k','BackgroundColor','w','FontSize',zoomboxnamefontsize,'Margin',zoomboxnamemargin)
%         elseif ii==3
%             text(zoomboxes_big(ii,3)/1e3-zoomboxnamebuffer,zoomboxes_big(ii,4)/1e3+zoomboxnamebuffer,zoomboxnames{ii},...
%                 'HorizontalAlignment','right','VerticalAlignment','bottom','Color','k','BackgroundColor','w','FontSize',zoomboxnamefontsize,'Margin',zoomboxnamemargin)
%         else
%             text(zoomboxes_big(ii,3)/1e3-zoomboxnamebuffer,zoomboxes_big(ii,2)/1e3-zoomboxnamebuffer,zoomboxnames{ii},...
%                 'HorizontalAlignment','right','VerticalAlignment','top','Color','k','BackgroundColor','w','FontSize',zoomboxnamefontsize,'Margin',zoomboxnamemargin)
%         end
%         % scale bar:
%         patch('Vertices',[zoomboxes_big(ii,3)/1e3-zoomboxscalebarbufferx(ii)-zoomratio(ii)*scalebarlength*[1;0;0;1;1],zoomboxes_big(ii,2)/1e3+zoomboxscalebarbuffery(ii)+scalebarheight*[0;0;1;1;0]],'Faces',1:5,'EdgeColor','k','FaceColor','w')
%         text(zoomboxes_big(ii,3)/1e3-zoomboxscalebarbufferx(ii)-.5*zoomratio(ii)*scalebarlength,zoomboxes_big(ii,2)/1e3+zoomboxscalebarbuffery(ii),[num2str(scalebarlength),' km'],...
%             'HorizontalAlignment','center','VerticalAlignment','bottom','FontSize',zoomboxscalebarfontsize)
%     end
% end
% % Set up axes:
% xlim(xlims/1000)
% ylim(ylims/1000)
% set(gca,'XTick',xtick*[ceil(xlims(1)/(1e3*xtick)):1:floor(xlims(2)/(1e3*xtick))])
% set(gca,'YTick',ytick*[ceil(ylims(1)/(1e3*ytick)):1:floor(ylims(2)/(1e3*ytick))])
% if thismodel==1
%     ylabel('Northing (km)','FontSize',axesfontsize)
% else
%     set(gca,'YTickLabel',[])
% end
% xlabel('Northing (km)')
% set(gca,'PlotBoxAspectRatio',[(xlims(2)-xlims(1))/(ylims(2)-ylims(1)),1,1])
% set(gca,'TickDir','out')
% set(gca,'FontSize',axesfontsize)
% %set(gca,'PositionConstraint','innerposition')
% % Set up color scale:
% if thismodel==1
%     eval(['thiscmap=colormap(gca,',dragcmap,');'])
%     if flipdragcmap
%         colormap(gca,flipud(thiscmap))
%     else
%         colormap(gca,thiscmap)
%     end
%     caxis(gca,draglims)
% else
%     eval(['thiscmap=colormap(gca,',dragdiffcmap,');'])
%     colormap(gca,thiscmap)
%     caxis(gca,dragdifflims)
% end
% hc=colorbar;
% if thismodel==1
%     set(get(hc,'Ylabel'),'String','kPa','FontSize',colorbarfontsize)
%     set(hc,'YTick',dragtick*[ceil(draglims(1)/dragtick):1:floor(draglims(2)/dragtick)])
% else
%     set(get(hc,'Ylabel'),'String','kPa','FontSize',colorbarfontsize)
%     set(hc,'YTick',dragdifftick*[ceil(dragdifflims(1)/dragdifftick):1:floor(dragdifflims(2)/dragdifftick)])
% end
% set(hc,'TickDir','out')
% set(hc,'FontSize',colorbarfontsize)
% % Make title:
% if thismodel==1
%     title([char(96+numlcurves+thismodel),') Drag (\tau_b) for ',descriptext{thismodel}],'FontSize',titlefontsize)
% else
%     title([char(96+numlcurves+thismodel),') Drag Difference (\tau_{b}(m)-\tau_{b}(m=1)) for ',descriptext{thismodel}],'FontSize',titlefontsize)
% end