% InversionThermalSettingFigure_v1

% Mike Wolovick, 2/28/2023

% This script makes the thermal setting figure for the inversion paper.  

% Subplots:  
% a) surftemp
% b) accum
% c) gflux+shearheat
% d) basaltemp
% e) meltrate
% f) bbar

clearvars
tic
close all

%% Parameters:

% File names:
thermalfile='/work/ollie/mwolovic/Inputs/FRISDomain1DMelt.mat';
boundaryfile='/work/ollie/mwolovic/FRISModeling/Boundaries/BoundaryPickingTest_v3_50km_expanded.exp';
modelfile='/work/ollie/mwolovic/FRISModeling/LCurves/FRISDomainLCurve_ResTest_SSA_drag_v2_mesh1a.mat';
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)
figname='/home/ollie/mwolovic/InversionThermalSettingFigure_v2.png'; 

% Ticks and limits and colorscales:
xtick=500;                               % km
ytick=500;                               % km
surftemplims=[-60,-10];                    % [1x2] deg C
surftemptick=10;                           % deg C
surfcmap='GMA_default';                     % valid colormap string, potentially including names of cmapfiles
flipsurfcmap=0;                                 % logical
accumlims=[0,60];                    % [1x2] cm/yr
accumtick=10;                           % cm/yr
accumcmap='parula';                     % valid colormap string, potentially including names of cmapfiles
flipaccumcmap=0;                                 % logical
heatflowlims=[29.9,1e4];          % [1x2] mW/m^2
heatflowcmap='GMA_ocean';       % valid colormap string, potentially including names of cmapfiles
flipheatflowcmap=0;                                 % logical
basaltemplims=[-25,0];                    % [1x2] deg C
basaltemptick=5;                           % deg C
basalcmap='GMA_default';                     % valid colormap string, potentially including names of cmapfiles
flipbasalcmap=0;                                 % logical
meltlims=[1e-4,1.001];          % [1x2] m/yr
meltcmap='GMA_ocean';       % valid colormap string, potentially including names of cmapfiles
flipmeltcmap=0;                                 % logical
bbarlims=[0,2];                    % [1x2] MPa*yr^-1/3
bbartick=0.5;                           % MPa*yr^-1/3
bbarcmap='GMA_default';                     % valid colormap string, potentially including names of cmapfiles
flipbbarcmap=0;                                 % logical

% Other aesthetic settings:
boundarylinewidth=2;                % points
boundarylinecolor='r';                % valid color string or RGB
masklinewidth=1;                       % points
maskcolor='w';                            % valid color string or RGB
axesfontsize=13;                       % points
titlefontsize=14;                          % points
colorbarfontsize=13;                   % points

% Subplot controls:
vertaxesbuffer=.09;             % unitless
titlebuffer=.05;                     % unitless
horzaxesbuffer=.03;             % unitless
horznotextbuffer=0;          % unitless
horzcolorbarbuffer=.09;        % unitless

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

% Seconds per year:
secondsperyear=60*60*24*365.25;

%% Preparation:

% Hardcoded subplot arrangement:
numvertplots=2;
numhorzplots=3;
% Compute subplot dimensions:
% Compute available space:
availablehorzspace=1-(numhorzplots-1)*horznotextbuffer-horzaxesbuffer-numhorzplots*horzcolorbarbuffer;
availablevertspace=1-numvertplots*titlebuffer-vertaxesbuffer;
% Compute boxes:
Boxes=cell(numvertplots,numhorzplots);
for ii=1:numvertplots
    for jj=1:numhorzplots
        Boxes{numvertplots-ii+1,jj}=[horzaxesbuffer+(jj-1)*(horznotextbuffer+horzcolorbarbuffer+availablehorzspace/numhorzplots),...
            vertaxesbuffer+(ii-1)*(titlebuffer+availablevertspace/numvertplots),...
            availablehorzspace/numhorzplots,availablevertspace/numvertplots];
    end
end
% Modify boxes for colorbars:
for ii=1:numvertplots*numhorzplots
    Boxes{ii}=Boxes{ii}+[0,0,horzcolorbarbuffer,0];
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

% Load domain boundary:
DomainBoundary=expread(boundaryfile);

% Compute axes limits:
xlims=[min(DomainBoundary.x/1000),max(DomainBoundary.x/1000)];
ylims=[min(DomainBoundary.y/1000),max(DomainBoundary.y/1000)];

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

% Load thermal file grids:
load(thermalfile,'X','Y','Accum','SurfTemp','BasalTemp','Gflux','ShearHeat','MeltRate','Mask','MaterialParameters')

% Evaluate grid size:
[ysize,xsize]=size(Mask);

% Generate contour structure for mask:
C_main=contourc(X/1000,Y/1000,Mask,[.5,1.5,2.5]);

% Load model:
load(modelfile,'BasicModel')

%% Make Figure:

% Make a figure:
figure(1)

% Make faded jet:
MakeFadedJet;


% Surface temperature plot:
% Make subplot:
subplot('Position',Boxes{1,1})
% Generate RGB:
% pre-allocate:
RGB=zeros(ysize,xsize,3);
% get this colormap:
eval(['thiscmap=colormap(gca,',surfcmap,');'])
if flipsurfcmap
    thiscmap=flipud(thiscmap);
end
% loop through RGB:
for ii=1:3
    RGB(:,:,ii)=interp1(linspace(0,1,size(thiscmap,1))',thiscmap(:,ii),max(0,min(1,(SurfTemp-surftemplims(1))/(surftemplims(2)-surftemplims(1)))),'linear');
end
% set color outside of the ice:
RGB(repmat(Mask==0,[1,1,3]))=.75;
% Draw variable:
image(X/1000,Y/1000,RGB)
set(gca,'ydir','normal')
hold on
% Plot grounding line and ice front:
startind=1;
done=0;
while done==0
    stopind=startind+C_main(2,startind);
    plot(C_main(1,startind+1:stopind),C_main(2,startind+1:stopind),'Color',maskcolor,'LineWidth',masklinewidth)
    startind=stopind+1;
    if startind>size(C_main,2)
        done=1;
    end
end
% Plot domain outline:
plot(DomainBoundary.x/1000,DomainBoundary.y/1000,'Color',boundarylinecolor,'LineWidth',boundarylinewidth)
% Set up axes:
xlim(xlims)
ylim(ylims)
set(gca,'XTick',xtick*[ceil(xlims(1)/xtick):1:floor(xlims(2)/xtick)])
set(gca,'YTick',ytick*[ceil(ylims(1)/ytick):1:floor(ylims(2)/ytick)])
ylabel('Northing (km)','FontSize',axesfontsize)
set(gca,'XTickLabel',[])
set(gca,'PlotBoxAspectRatio',[(xlims(2)-xlims(1))/(ylims(2)-ylims(1)),1,1])
set(gca,'TickDir','out')
set(gca,'FontSize',axesfontsize)
% Set up color scale and colorbar:
colormap(gca,thiscmap)
caxis(surftemplims)
hc=colorbar;
set(hc,'YTick',surftemptick*[ceil(surftemplims(1)/surftemptick):1:floor(surftemplims(2)/surftemptick)])
set(get(hc,'Ylabel'),'String','^{\circ}C','FontSize',colorbarfontsize)
set(hc,'TickDir','out')
set(hc,'FontSize',colorbarfontsize)
% Make title:
title('a) Surface temperature','FontSize',titlefontsize)


% accumulation plot:
% Make subplot:
subplot('Position',Boxes{1,2})
% Generate RGB:
% pre-allocate:
RGB=zeros(ysize,xsize,3);
% get this colormap:
eval(['thiscmap=colormap(gca,',accumcmap,');'])
if flipaccumcmap
    thiscmap=flipud(thiscmap);
end
% loop through RGB:
for ii=1:3
    RGB(:,:,ii)=interp1(linspace(0,1,size(thiscmap,1))',thiscmap(:,ii),max(0,min(1,(Accum*100-accumlims(1))/(accumlims(2)-accumlims(1)))),'linear');
end
% set color outside of the ice:
RGB(repmat(Mask==0|Mask==3,[1,1,3]))=.75;
% Draw variable:
image(X/1000,Y/1000,RGB)
set(gca,'ydir','normal')
hold on
% Plot grounding line and ice front:
startind=1;
done=0;
while done==0
    stopind=startind+C_main(2,startind);
    plot(C_main(1,startind+1:stopind),C_main(2,startind+1:stopind),'Color',maskcolor,'LineWidth',masklinewidth)
    startind=stopind+1;
    if startind>size(C_main,2)
        done=1;
    end
end
% Plot domain outline:
plot(DomainBoundary.x/1000,DomainBoundary.y/1000,'Color',boundarylinecolor,'LineWidth',boundarylinewidth)
% Set up axes:
xlim(xlims)
ylim(ylims)
set(gca,'XTick',xtick*[ceil(xlims(1)/xtick):1:floor(xlims(2)/xtick)])
set(gca,'YTick',ytick*[ceil(ylims(1)/ytick):1:floor(ylims(2)/ytick)])
set(gca,'YTickLabel',[])
set(gca,'XTickLabel',[])
set(gca,'PlotBoxAspectRatio',[(xlims(2)-xlims(1))/(ylims(2)-ylims(1)),1,1])
set(gca,'TickDir','out')
set(gca,'FontSize',axesfontsize)
% Set up color scale and colorbar:
colormap(gca,thiscmap)
caxis(accumlims)
hc=colorbar;
set(hc,'YTick',accumtick*[ceil(accumlims(1)/accumtick):1:floor(accumlims(2)/accumtick)])
set(get(hc,'Ylabel'),'String','cm/yr','FontSize',colorbarfontsize)
set(hc,'TickDir','out')
set(hc,'FontSize',colorbarfontsize)
% Make title:
title('b) Accumulation rate','FontSize',titlefontsize)


% heat flow plot:
% Make subplot:
subplot('Position',Boxes{1,3})
% Generate RGB:
% pre-allocate:
RGB=zeros(ysize,xsize,3);
% get this colormap:
eval(['thiscmap=colormap(gca,',heatflowcmap,');'])
if flipheatflowcmap
    thiscmap=flipud(thiscmap);
end
% loop through RGB:
for ii=1:3
    RGB(:,:,ii)=interp1(linspace(0,1,size(thiscmap,1))',thiscmap(:,ii),max(0,min(1,(log((Gflux+ShearHeat)*1000)-log(heatflowlims(1)))/(log(heatflowlims(2))-log(heatflowlims(1))))),'linear');
end
% set color outside of the grounded ice:
RGB(repmat(Mask~=2,[1,1,3]))=.75;
% Draw image:
image(X/1000,Y/1000,RGB)
set(gca,'ydir','normal')
hold on
% Plot grounding line and ice front:
startind=1;
done=0;
while done==0
    stopind=startind+C_main(2,startind);
    plot(C_main(1,startind+1:stopind),C_main(2,startind+1:stopind),'Color',maskcolor,'LineWidth',masklinewidth)
    startind=stopind+1;
    if startind>size(C_main,2)
        done=1;
    end
end
% Plot domain outline:
plot(DomainBoundary.x/1000,DomainBoundary.y/1000,'Color',boundarylinecolor,'LineWidth',boundarylinewidth)
% Set up axes:
set(gca,'ydir','normal')
xlim(xlims)
ylim(ylims)
set(gca,'XTick',xtick*[ceil(xlims(1)/xtick):1:floor(xlims(2)/xtick)])
set(gca,'YTick',ytick*[ceil(ylims(1)/ytick):1:floor(ylims(2)/ytick)])
set(gca,'XTickLabel',[])
set(gca,'YTickLabel',[])
set(gca,'PlotBoxAspectRatio',[(xlims(2)-xlims(1))/(ylims(2)-ylims(1)),1,1])
set(gca,'TickDir','out')
set(gca,'FontSize',axesfontsize)
% Set up color scale:
colormap(gca,thiscmap)
caxis(log10(heatflowlims))
% Compute first base and tick for color bar:
thisbase=10^floor(log10(heatflowlims(1)));
thistick=thisbase*ceil(heatflowlims(1)/thisbase);
% Create all ticks and labels for color bar:
heatflowticks=zeros(100,1);
heatflowticklabels=cell(100,1);
for ii=1:100
    % Record this tick:
    heatflowticks(ii)=log10(thistick);
    % Check if we hit the next power of ten:
    if abs(thistick-10*thisbase)/thisbase<1e-12 || (ii==1 && abs(thistick-thisbase)/thisbase<1e-12)
        % Generate label:
        heatflowticklabels{ii}=['10^{',num2str(round(log10(thistick))),'}'];
        % Advance base:
        thisbase=thistick;
    end
    % Advance tick:
    thistick=thistick+thisbase;
    % Break from loop and truncate:
    if thistick>(1+1e-12)*heatflowlims(2)
        heatflowticks=heatflowticks(1:ii);
        heatflowticklabels=heatflowticklabels(1:ii);
        break
    end
end
% Make colorbar:
hc=colorbar;
set(hc,'YTick',heatflowticks)
set(hc,'YTickLabel',heatflowticklabels)
set(get(hc,'Ylabel'),'String','mW/m^2','FontSize',colorbarfontsize)
set(hc,'TickDir','out')
set(hc,'FontSize',colorbarfontsize)
% Make title:
title('c) Combined heat flow','FontSize',titlefontsize)


% basal temperature plot:
% Make subplot:
subplot('Position',Boxes{2,1})
% Generate RGB:
% pre-allocate:
RGB=zeros(ysize,xsize,3);
% get this colormap:
eval(['thiscmap=colormap(gca,',basalcmap,');'])
if flipbasalcmap
    thiscmap=flipud(thiscmap);
end
% loop through RGB:
for ii=1:3
    RGB(:,:,ii)=interp1(linspace(0,1,size(thiscmap,1))',thiscmap(:,ii),max(0,min(1,(BasalTemp-basaltemplims(1))/(basaltemplims(2)-basaltemplims(1)))),'linear');
end
% set color outside of the ice:
RGB(repmat(Mask==0|Mask==3,[1,1,3]))=.75;
% Draw variable:
image(X/1000,Y/1000,RGB)
set(gca,'ydir','normal')
hold on
% Plot grounding line and ice front:
startind=1;
done=0;
while done==0
    stopind=startind+C_main(2,startind);
    plot(C_main(1,startind+1:stopind),C_main(2,startind+1:stopind),'Color',maskcolor,'LineWidth',masklinewidth)
    startind=stopind+1;
    if startind>size(C_main,2)
        done=1;
    end
end
% Plot domain outline:
plot(DomainBoundary.x/1000,DomainBoundary.y/1000,'Color',boundarylinecolor,'LineWidth',boundarylinewidth)
% Set up axes:
xlim(xlims)
ylim(ylims)
set(gca,'XTick',xtick*[ceil(xlims(1)/xtick):1:floor(xlims(2)/xtick)])
set(gca,'YTick',ytick*[ceil(ylims(1)/ytick):1:floor(ylims(2)/ytick)])
ylabel('Northing (km)','FontSize',axesfontsize)
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 up color scale and colorbar:
colormap(gca,thiscmap)
caxis(basaltemplims)
hc=colorbar;
set(hc,'YTick',basaltemptick*[ceil(basaltemplims(1)/basaltemptick):1:floor(basaltemplims(2)/basaltemptick)])
set(get(hc,'Ylabel'),'String','^{\circ}C','FontSize',colorbarfontsize)
set(hc,'TickDir','out')
set(hc,'FontSize',colorbarfontsize)
% Make title:
title('d) Basal temperature','FontSize',titlefontsize)


% melt rate plot:
% Make subplot:
subplot('Position',Boxes{2,2})
% Generate RGB:
% pre-allocate:
RGB=zeros(ysize,xsize,3);
% get this colormap:
eval(['thiscmap=colormap(gca,',meltcmap,');'])
if flipmeltcmap
    thiscmap=flipud(thiscmap);
end
% loop through RGB:
for ii=1:3
    RGB(:,:,ii)=interp1(linspace(0,1,size(thiscmap,1))',thiscmap(:,ii),max(0,min(1,(log(max(0,MeltRate))-log(meltlims(1)))/(log(meltlims(2))-log(meltlims(1))))),'linear');
end
% set color outside of the grounded ice:
RGB(repmat(Mask~=2,[1,1,3]))=.75;
% Draw image:
image(X/1000,Y/1000,RGB)
set(gca,'ydir','normal')
hold on
% Plot grounding line and ice front:
startind=1;
done=0;
while done==0
    stopind=startind+C_main(2,startind);
    plot(C_main(1,startind+1:stopind),C_main(2,startind+1:stopind),'Color',maskcolor,'LineWidth',masklinewidth)
    startind=stopind+1;
    if startind>size(C_main,2)
        done=1;
    end
end
% Plot domain outline:
plot(DomainBoundary.x/1000,DomainBoundary.y/1000,'Color',boundarylinecolor,'LineWidth',boundarylinewidth)
% Set up axes:
set(gca,'ydir','normal')
xlim(xlims)
ylim(ylims)
set(gca,'XTick',xtick*[ceil(xlims(1)/xtick):1:floor(xlims(2)/xtick)])
set(gca,'YTick',ytick*[ceil(ylims(1)/ytick):1:floor(ylims(2)/ytick)])
xlabel('Easting (km)','FontSize',axesfontsize)
set(gca,'YTickLabel',[])
set(gca,'PlotBoxAspectRatio',[(xlims(2)-xlims(1))/(ylims(2)-ylims(1)),1,1])
set(gca,'TickDir','out')
set(gca,'FontSize',axesfontsize)
% Set up color scale:
colormap(gca,thiscmap)
caxis(log10(meltlims))
% Compute first base and tick for color bar:
thisbase=10^floor(log10(meltlims(1)));
thistick=thisbase*ceil(meltlims(1)/thisbase);
% Create all ticks and labels for color bar:
meltticks=zeros(100,1);
meltticklabels=cell(100,1);
for ii=1:100
    % Record this tick:
    meltticks(ii)=log10(thistick);
    % Check if we hit the next power of ten:
    if abs(thistick-10*thisbase)/thisbase<1e-12 || (ii==1 && abs(thistick-thisbase)/thisbase<1e-12)
        % Generate label:
        meltticklabels{ii}=['10^{',num2str(round(log10(thistick))),'}'];
        % Advance base:
        thisbase=thistick;
    end
    % Advance tick:
    thistick=thistick+thisbase;
    % Break from loop and truncate:
    if thistick>(1+1e-12)*meltlims(2)
        meltticks=meltticks(1:ii);
        meltticklabels=meltticklabels(1:ii);
        break
    end
end
% Make colorbar:
hc=colorbar;
set(hc,'YTick',meltticks)
set(hc,'YTickLabel',meltticklabels)
set(get(hc,'Ylabel'),'String','m/yr','FontSize',colorbarfontsize)
set(hc,'TickDir','out')
set(hc,'FontSize',colorbarfontsize)
% Make title:
title('e) Basal melt rate','FontSize',titlefontsize)


% rheology B
% Make subplot:
subplot('Position',Boxes{2,3})
% get this colormap:
eval(['thiscmap=colormap(gca,',bbarcmap,');'])
if flipbbarcmap
    thiscmap=flipud(thiscmap);
end
% Make back plane patch:
patch('Vertices',[[xlims(1);xlims(2);xlims(2);xlims(1);xlims(1)],[ylims(1);ylims(1);ylims(2);ylims(2);ylims(1)]],'Faces',1:5,'FaceColor',[.75,.75,.75],'EdgeColor','none')
% Make data patch:
IceFraction_elements=mean(BasicModel.mask.ice_levelset(BasicModel.mesh.elements)<0,2);
patch('Vertices',[BasicModel.mesh.x/1000,BasicModel.mesh.y/1000],'Faces',BasicModel.mesh.elements(IceFraction_elements>.5,:),'EdgeColor','none','FaceColor','interp','FaceVertexCData',BasicModel.materials.rheology_B*1e-6*(secondsperyear^-(1/3)))
hold on
% Plot grounding line and ice front:
startind=1;
done=0;
while done==0
    stopind=startind+C_main(2,startind);
    plot(C_main(1,startind+1:stopind),C_main(2,startind+1:stopind),'Color',maskcolor,'LineWidth',masklinewidth)
    startind=stopind+1;
    if startind>size(C_main,2)
        done=1;
    end
end
% Set up axes:
xlim(xlims)
ylim(ylims)
set(gca,'XTick',xtick*[ceil(xlims(1)/xtick):1:floor(xlims(2)/xtick)])
set(gca,'YTick',ytick*[ceil(ylims(1)/ytick):1:floor(ylims(2)/ytick)])
xlabel('Easting (km)','FontSize',axesfontsize)
set(gca,'YTickLabel',[])
set(gca,'PlotBoxAspectRatio',[(xlims(2)-xlims(1))/(ylims(2)-ylims(1)),1,1])
set(gca,'TickDir','out')
set(gca,'FontSize',axesfontsize)
% Set up color scale and colorbar:
colormap(gca,thiscmap)
caxis(bbarlims)
hc=colorbar;
set(hc,'YTick',bbartick*[ceil(bbarlims(1)/bbartick):1:floor(bbarlims(2)/bbartick)])
set(get(hc,'Ylabel'),'String','MPa yr^{-1/3}','FontSize',colorbarfontsize)
set(hc,'TickDir','out')
set(hc,'FontSize',colorbarfontsize)
% Make title:
title('f) Rheological stiffness (B)','FontSize',titlefontsize)


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

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