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

% NComparisonFigure:  just makes a figure showing the different maps of N

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_op_v1_mesh1.mat';...
    '/work/ollie/mwolovic/FRISModeling/LCurves/FRISDomainLCurve_ResTest_SSA_drag_opc_v1_mesh1.mat';...
    '/work/ollie/mwolovic/FRISModeling/LCurves/FRISDomainLCurve_ResTest_SSA_cuas_drag_v1a_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='/home/ollie/mwolovic/NComparisonFigure_v2.png';

% N descriptions:
ndescriptext={'N_{op}','N_{opc}','N_{cuas}'};  % 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.

% Ticks and limits:
xtick=500;                               % km
ytick=500;                               % km
scalebarlength=50;                % km (km of the zoom plot)
scalebarheight=40;                  % km (km of the main plot)
nlims=[0,40;0,40;0,10];                   % [Nx2] MPa
ntick=[10,10,2];                             % [1xN] MPa
ncmap='GMA_default';         % valid colormap string, potentially including names of cmapfiles

% Zoom plot controls:
showzoombox=0;                          % logical
showzoomboxloc=0;                     % 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:
backplanecolor=[.75,.75,.75];       % valid color string or RGB
titlefontsize=14;                              % points
axesfontsize=13;                              % points
colorbarfontsize=13;                         % 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=.10;             % unitless
titlebuffer_map=.04;                     % unitless
horzaxesbuffer_map=.05;             % unitless
horznotextbuffer_map=.01;          % unitless
horzcolorbarbuffer_map=.02;        % unitless

% Output size and resolution:
pagewidth_map=18;                       % inches (height computed automatically)
resolution=300;                     % dpi

%% Preparation:

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

% Hardcoded subplot arrangement:
numvertplots=1;
numhorzplots=numn;
% Compute subplot dimensions:
% Compute available space:
availablehorzspace=1-(numhorzplots-1)*horznotextbuffer_map-horzaxesbuffer_map-numhorzplots*horzcolorbarbuffer_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+horzcolorbarbuffer_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);

% Compute page height:
pageheight=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,numn);

% Loop through models:
for thisn=1:numn

    % Check if we're loading remote files:
    if exist('serverloc','var')
        % Parse filename:
        slashnum=strfind(lcurvefiles{thisn},'/');
        remotepath=lcurvefiles{thisn}(1:slashnum(end));
        lcurvename=lcurvefiles{thisn}(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],'BasicModel')
    else
        % Load the basic file:
        load(lcurvefiles{thisn},'BasicModel')
    end

    % Check mesh size:
    if numlayers>1 && ~strcmp(flowequation,'SSA')
        if BasicModel.mesh.numberofvertices2d~=numvertices2d || BasicModel.mesh.numberofelements2d~=numelements2d
            error(['Mesh number ',num2str(thisn),' does not match reference mesh.'])
        end
    else
        if BasicModel.mesh.numberofvertices~=numvertices2d || BasicModel.mesh.numberofelements~=numelements2d
            error(['Mesh number ',num2str(thisn),' does not match reference mesh.'])
        end
    end
    
    % Extract effective pressure:
    EffectivePressure_vertices(:,thisn)=BasicModel.friction.effective_pressure(1:numvertices2d);

end


%% Make Map Figure:

% Make the figure:
figure

% Make faded jet:
MakeFadedJet;

% Loop through models:
for thisn=1:numn
    
    % Plot effective pressure:
    % Make subplot:
    subplot('Position',Boxes_maps{1,thisn})
    % 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',EffectivePressure_vertices(:,thisn)/1e6)
    % 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',EffectivePressure_vertices(:,thisn)/1e6)
            % 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 thisn==1
        ylabel('Northing (km)','FontSize',axesfontsize)
    else
        set(gca,'YTickLabel',[])
    end
    xlabel('Easting (km)','FontSize',axesfontsize)
    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,',ncmap,');'])
    colormap(gca,thiscmap)
    caxis(gca,nlims(thisn,:))
    hc=colorbar;
    set(get(hc,'Ylabel'),'String','MPa','FontSize',colorbarfontsize)
    set(hc,'YTick',ntick(thisn)*[ceil(nlims(thisn,1)/ntick(thisn)):1:floor(nlims(thisn,2)/ntick(thisn))])
    set(hc,'TickDir','out')
    set(hc,'FontSize',colorbarfontsize)
    % Make title:
    title([char(96+thisn),') Effective Pressure for ',ndescriptext{thisn}],'FontSize',titlefontsize)
    
end

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


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

%% Cutting Room Floor:

% %% 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);
% 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);
% LCurveData=struct('LogLambda_points',[],'LogRegCost_points',[],'LogDataCost_points',[],'IsOutlier',[],'LogRegCost_curve',[],'LogDataCost_curve',[],...
%     'LogDataCurvature_curve',[],'LogRegCurvature_curve',[],'TotalCurvature_curve',[],...
%     'bestlambdaind_maxallowable',[],'bestlambdaind_minallowable',[],'bestlambdaind',[]);
% 
% % Record reference parameters:
% ReferenceInversionParameters=MyInversionParameters;
% 
% % 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)
% 
%     % 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);
%         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;
%         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);
% 
%     % 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);
%     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);
% 
% %% 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);
% 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:
%     meanc(ii,1)=sum(Area.*ActiveMask_elements.*GoodDataWeight_elements.*DragCoeff_elements(:,ii))/totalstatsarea(1);
%     rmsc(ii,1)=sqrt(sum(Area.*ActiveMask_elements.*GoodDataWeight_elements.*((DragCoeff_elements(:,ii)-meanc(ii,1)).^2))/totalstatsarea(1));
%     meanc(ii,2)=sum(Area.*ActiveMask_elements.*FastFlowWeight_elements.*GoodDataWeight_elements.*DragCoeff_elements(:,ii))/totalstatsarea(2);
%     rmsc(ii,2)=sqrt(sum(Area.*ActiveMask_elements.*FastFlowWeight_elements.*GoodDataWeight_elements.*((DragCoeff_elements(:,ii)-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));
%     
%     % Compute correlations:
%     % R(N,Cw):
%     rnc(ii,1)=sum(Area.*ActiveMask_elements.*GoodDataWeight_elements.*(EffectivePressure_elements(:,ii)-meann(ii,1)).*(DragCoeff_elements(:,1)-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)-meanc(1,2)))/(totalstatsarea(2)*rmsn(ii,2)*rmsc(1,2));
%     % 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;
% r2nc=rnc.^2;
% r2lnnlnc=rlnnlnc.^2;
% r2lnclnc=rlnclnc.^2;
% r2tautau=rtautau.^2;
% r2ntau=rntau.^2;
% 
% % compute variance ratio:
% varratiolnclnc=varlnc./repmat(varlnc(1,:),[numlcurves,1]);
% 
% %% Communicate Stats:
% 
% 
% % Display stats, whole domain:
% string1='RMS(ln(C))=';
% string2='var(ln(C))=';
% string3='Variance ratio=';
% 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=';
% 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))),', '];
% 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'])
% 
% % Display stats, fast flow:
% string1='RMS(ln(C))=';
% string2='var(ln(C))=';
% string3='Variance ratio=';
% 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=';
% 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))),', '];
% 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('...')
% 
% toc 
% return
