%Replicating the buttressing figures of Furst et al. (2016) - The Safety Band of Antarctic Ice Shelves
glacier = 'Ronne';
%load(['Models/Model_' glacier '_CalvingVM_local.mat'])
%Quiver plot
%plotmodel(md,'data',[md.results.StressbalanceSolution.Vx md.results.StressbalanceSolution.Vy])

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if strcmp(glacier,'Amery')
	PSIthresh=0.4;
elseif strcmp(glacier,'Denman')||strcmp(glacier,'Shackleton')
	PSIthresh=0.1;
elseif strcmp(glacier,'Filchner')||strcmp(glacier,'Ronne')
	PSIthresh=0.4;
elseif strcmp(glacier,'LarsenC')
	PSIthresh=0.3;
elseif strcmp(glacier,'PineIsland')
	PSIthresh=0.7;
elseif strcmp(glacier,'Ross')
	PSIthresh=0.5;
elseif strcmp(glacier,'Thwaites')
	PSIthresh=0.4;
elseif strcmp(glacier,'Totten')
	PSIthresh=0.1;
end

%%%% Superimpose final ice front position from transient on buttressing plot %%%%%
if strcmp(glacier,'Amery')
   EC = 50000000;
   VM = 150000;
   sr = 100e3; %scale ruler
elseif strcmp(glacier,'Denman')
   EC = 20000000;
   VM = 400000;
   sr = 50e3;
elseif strcmp(glacier,'Filchner')
   EC = 80000000;
   VM = 190000;
   sr = 200e3;
elseif strcmp(glacier,'LarsenC')
   EC = 30000000000;
   VM = 140000;
   sr = 200e3;
elseif strcmp(glacier,'PineIsland')
   EC = 300000000;
   VM = 275000;
   sr = 50e3;
elseif strcmp(glacier,'Ronne')
   EC = 400000000;
   VM = 120000;
   sr = 200e3;
elseif strcmp(glacier,'Ross')
   EC = 300000000;
   VM = 105000;
   sr = 200e3;
elseif strcmp(glacier,'Shackleton')
   EC = 200000000;
   VM = 280000;
   sr = 50e3;
elseif strcmp(glacier,'Thwaites')
   EC = 200000000;
   VM = 320000;
elseif strcmp(glacier,'Totten')
   EC = 250000000;
   VM = 260000;
   sr = 50e3;
end

load(['Models/Model_' glacier '_CalvingLevermann' num2str(EC) '.mat']);

%find elements that are floating ice
pos_floating = find(min(md.mask.ocean_levelset(md.mesh.elements),[],2)<0 & min(md.mask.ice_levelset(md.mesh.elements),[],2)<0);

%Recover stresses
vx = md.results.StressbalanceSolution.Vx;
vy = md.results.StressbalanceSolution.Vy;
vel = md.results.StressbalanceSolution.Vel;
md = mechanicalproperties(md, vx, vy);

%Constants
g = md.constants.g;
rho_i = md.materials.rho_ice;
rho_w = md.materials.rho_water;

%Element averaged thickness
H = mean(md.geometry.thickness(md.mesh.elements),2);
depth = mean(md.geometry.base(md.mesh.elements),2)-mean(md.geometry.bed(md.mesh.elements),2);
vx = mean(vx(md.mesh.elements),2);
vy = mean(vy(md.mesh.elements),2);
vel = mean(vel(md.mesh.elements),2);

%Allocate Kn
Kn = NaN(length(md.mesh.elements),1);
buttress = NaN(md.mesh.numberofelements,1);

for el=pos_floating'

	%Build stress tensor
   tau_xx = md.results.deviatoricstress.xx(el);
   tau_yy = md.results.deviatoricstress.yy(el);
   tau_xy = md.results.deviatoricstress.xy(el);
   R = [2*tau_xx+tau_yy   tau_xy;tau_xy   2*tau_yy+tau_xx];

	%Decide on outward pointing unit vector
	if 0 %"Flow buttressing" - the supplement
		vector = [vx(el); vy(el)];
	else %"Max buttressing" - the main text
		vector = md.results.deviatoricstress.principalaxis2(el,:)';
	end
	n = vector./sqrt(vector(1)^2+vector(2)^2);
	
	%Buildsigma_nn and then Kn
	N = n'*R*n;
   N0 = 0.5*g*rho_i.*(1-rho_i./rho_w).*H(el);
	Kn(el) = 1 - N/N0;
end


%Final eigencalving buttressing
clear md
load(['Models/Model_' glacier '_CalvingLevermann' num2str(EC) '.mat']);
vx = md.results.TransientSolution(end).Vx;
vy = md.results.TransientSolution(end).Vy;
vel = md.results.TransientSolution(end).Vel;

md = mechanicalproperties(md, vx, vy);
vx = mean(vx(md.mesh.elements),2);
vy = mean(vy(md.mesh.elements),2);
vel = mean(vel(md.mesh.elements),2);
H = mean(md.results.TransientSolution(end).Thickness(md.mesh.elements),2);
%Allocate Kn_ec
Kn_ec = NaN(length(md.mesh.elements),1);
buttress = NaN(md.mesh.numberofelements,1);
pos_floating_ec = find(min(md.results.TransientSolution(end).MaskOceanLevelset(md.mesh.elements),[],2)<0 & min(md.results.TransientSolution(end).MaskIceLevelset(md.mesh.elements),[],2)<0);

for el=pos_floating_ec'

   %Build stress tensor
   tau_xx = md.results.deviatoricstress.xx(el);
   tau_yy = md.results.deviatoricstress.yy(el);
   tau_xy = md.results.deviatoricstress.xy(el);
   R = [2*tau_xx+tau_yy   tau_xy;tau_xy   2*tau_yy+tau_xx];

   %Decide on outward pointing unit vector
   if 0 %"Flow buttressing" - the supplement
      vector = [vx(el); vy(el)];
   else %"Max buttressing" - the main text
      vector = md.results.deviatoricstress.principalaxis2(el,:)';
   end
   n = vector./sqrt(vector(1)^2+vector(2)^2);

   %Buildsigma_nn and then Kn
   N = n'*R*n;
   N0 = 0.5*g*rho_i.*(1-rho_i./rho_w).*H(el);
   Kn_ec(el) = 1 - N/N0;
end

%Final von Mises buttressing
clear md vx vy vel H
load(['Models/Model_' glacier '_CalvingVM' num2str(VM) '.mat']);
vx = md.results.TransientSolution(end).Vx;
vy = md.results.TransientSolution(end).Vy;
vel = md.results.TransientSolution(end).Vel;
md = mechanicalproperties(md, vx, vy);
vx = mean(vx(md.mesh.elements),2);
vy = mean(vy(md.mesh.elements),2);
vel = mean(vel(md.mesh.elements),2);
H = mean(md.results.TransientSolution(end).Thickness(md.mesh.elements),2);
%Allocate Kn_vm
Kn_vm = NaN(length(md.mesh.elements),1);
buttress = NaN(md.mesh.numberofelements,1);
pos_floating_vm = find(min(md.results.TransientSolution(end).MaskOceanLevelset(md.mesh.elements),[],2)<0 & min(md.results.TransientSolution(end).MaskIceLevelset(md.mesh.elements),[],2)<0);

for el=pos_floating_vm'

   %Build stress tensor
   tau_xx = md.results.deviatoricstress.xx(el);
   tau_yy = md.results.deviatoricstress.yy(el);
   tau_xy = md.results.deviatoricstress.xy(el);
   R = [2*tau_xx+tau_yy   tau_xy;tau_xy   2*tau_yy+tau_xx];

   %Decide on outward pointing unit vector
   if 0 %"Flow buttressing" - the supplement
      vector = [vx(el); vy(el)];
   else %"Max buttressing" - the main text
      vector = md.results.deviatoricstress.principalaxis2(el,:)';
   end
   n = vector./sqrt(vector(1)^2+vector(2)^2);

   %Buildsigma_nn and then Kn
   N = n'*R*n;
   N0 = 0.5*g*rho_i.*(1-rho_i./rho_w).*H(el);
   Kn_vm(el) = 1 - N/N0;
end

%Calculate front-averaged parameters
%figure(1)
%plotmodel(md,'data',Kn,'caxis',[0 1],'overlay',1,'colormap',jet(12),'smooth',1);

pos_front = find(min(md.mask.ice_levelset(md.mesh.elements),[],2)<0 & max(md.mask.ice_levelset(md.mesh.elements),[],2)>0);

maxCFthick = max(H(pos_front));
maxCFvel = max(vel(pos_front));
maxCFbuttress = max(Kn(pos_front));
maxCFdepth = max(depth(pos_front));
minCFbuttress = min(Kn(pos_front));
minCFdepth = min(depth(pos_front));
meanCFthick = mean(H(pos_front),'omitnan');
meanCFvel = mean(vel(pos_front),'omitnan');
meanCFbuttress = mean(Kn(pos_front),'omitnan');
meanCFdepth = mean(depth(pos_front));
medCFthick = median(H(pos_front),'omitnan');
medCFvel = median(vel(pos_front),'omitnan');
medCFbuttress = median(Kn(pos_front),'omitnan');
medCFdepth = median(depth(pos_front));

meanFLOATINGbuttress = mean(Kn(pos_floating));
medFLOATINGbuttress = median(Kn(pos_floating));
%propFLOATINGbuttressLessThan1 = length(nonzeros(Kn(pos_floating)<1))/length(nonzeros(Kn(pos_floating)));
areas=GetAreas(md.mesh.elements,md.mesh.x,md.mesh.y);
%propareaslessPSI = sum(areas(Kn(pos_floating)<PSIthresh))/sum(areas(pos_floating)); %summed area of elements less than PSI
%propareaslessPSI = sum(areas(Kn(pos_floating)<PSIthresh))/sum(areas(Kn(pos_floating)<1|Kn(pos_floating)>=1)); %summed area of elements less than PSI
%Kn_noNaN = Kn(isnan(Kn)==0)
floatingareas = areas(pos_floating);
floatingareas_lessPSI = floatingareas(Kn(pos_floating)<PSIthresh);
floatingareas_less01 = floatingareas(Kn(pos_floating)<0.1);
floatingareas_less02 = floatingareas(Kn(pos_floating)<0.2);
floatingareas_less03 = floatingareas(Kn(pos_floating)<0.3);
floatingareas_less04 = floatingareas(Kn(pos_floating)<0.4);
floatingareas_less05 = floatingareas(Kn(pos_floating)<0.5);
floatingareas_less06 = floatingareas(Kn(pos_floating)<0.6);
floatingareas_less07 = floatingareas(Kn(pos_floating)<0.7);
floatingareas_less08 = floatingareas(Kn(pos_floating)<0.8);
floatingareas_less09 = floatingareas(Kn(pos_floating)<0.9);
floatingareas_lessONE = floatingareas(Kn(pos_floating)<1.0);
propareaslessPSI = sum(floatingareas_lessPSI)/sum(floatingareas);
propareasless01 = sum(floatingareas_less01)/sum(floatingareas);
propareasless02 = sum(floatingareas_less02)/sum(floatingareas);
propareasless03 = sum(floatingareas_less03)/sum(floatingareas);
propareasless04 = sum(floatingareas_less04)/sum(floatingareas);
propareasless05 = sum(floatingareas_less05)/sum(floatingareas);
propareasless06 = sum(floatingareas_less06)/sum(floatingareas);
propareasless07 = sum(floatingareas_less07)/sum(floatingareas);
propareasless08 = sum(floatingareas_less08)/sum(floatingareas);
propareasless09 = sum(floatingareas_less09)/sum(floatingareas);
propareaslessONE = sum(floatingareas_lessONE)/sum(floatingareas);

fprintf('Max calving front thickness is %f\n', maxCFthick);
fprintf('Max calving front velocity is %f\n', maxCFvel);
fprintf('Max calving front buttressing is %f\n', maxCFbuttress);
fprintf('Min calving front buttressing is %f\n', minCFbuttress);
fprintf('Mean calving front thickness is %f\n', meanCFthick);
fprintf('Mean calving front velocity is %f\n', meanCFvel);
fprintf('Mean calving front buttressing is %f\n', meanCFbuttress);
fprintf('Median calving front thickness is %f\n', medCFthick);
fprintf('Median calving front velocity is %f\n', medCFvel);
fprintf('Median calving front buttressing is %f\n', medCFbuttress);
disp('%%%% DEPTH %%%%')
fprintf('Max calving front depth is %f\n', maxCFdepth);
fprintf('Min calving front depth is %f\n', minCFdepth);
fprintf('Mean calving front depth is %f\n', meanCFdepth);
fprintf('Median calving front depth is %f\n', medCFdepth);
disp('%%%% FLOATING ELEMENTS %%%%')
fprintf('Mean floating ice buttressing is %f\n', meanFLOATINGbuttress);
fprintf('Median floating ice buttressing is %f\n', medFLOATINGbuttress);
fprintf('Proportion floating ice less than 1 is %f\n', propareaslessONE);
fprintf('Proportion floating ice less than PSI is %f\n', propareaslessPSI);
disp('%%%% FLOATING ELEMENTS, RATIO EXPERIMENT  %%%%')
fprintf('%f\n', propareasless01);
fprintf('%f\n', propareasless02);
fprintf('%f\n', propareasless03);
fprintf('%f\n', propareasless04);
fprintf('%f\n', propareasless05);
fprintf('%f\n', propareasless06);
fprintf('%f\n', propareasless07);
fprintf('%f\n', propareasless08);
fprintf('%f\n', propareasless09);
fprintf('%f\n', propareaslessONE);
%fprintf('Centered calving front thickness is %f\n', C_CFthick);
%fprintf('Centered calving front velocity is %f\n', C_CFvel);
%fprintf('Centered calving front buttressing is %f\n', C_CFbuttress);


fig = figure(1);
md=radarpower(md);
[X Y]=meshgrid(md.radaroverlay.x,md.radaroverlay.y);
vmin=0; vmax = 1;
firstfront=isoline(md, md.mask.ice_levelset,'output','matrix');
options = pairoptions('tightsubplot',1,'gap',[.005 .005],'vmargin',[.05 .05],'hmargin',[.02 .02]);
set(gcf,'color','w','Position',[600 200 1100 700]);clf;
cm = getcolormap(pairoptions('colormap',copper(12)));
colorpalette=lines(5); colorpalette(2,:) = [1 0 0];

md1 = loadmodel(['Models/Model_' glacier '_CalvingLevermann' num2str(EC) '.mat']);
md2 = loadmodel(['Models/Model_' glacier '_CalvingVM' num2str(VM) '.mat']);

switch glacier
	case 'Amery'
		percent_x_max=80;
	otherwise
		percent_x_max=100;
end
percent_x_min = 0; percent_y_min = 0; percent_y_max = 100;
x_mesh_min = (percent_x_min * (max(md.mesh.x) - min(md.mesh.x)) / 100) + min(md.mesh.x);
x_mesh_max = (percent_x_max * (max(md.mesh.x) - min(md.mesh.x)) / 100) + min(md.mesh.x);
idx_x_max = md.mesh.x(:)<=x_mesh_max;
idx_x_min = md.mesh.x(:)>=x_mesh_min;
idx_X_min = X(1,:)>=x_mesh_min;
idx_X_max = X(1,:)<=x_mesh_max;
X_min = find(idx_X_min,1,'first');
X_max = find(idx_X_max,1,'last');
mesh_x_subset = intersect(md.mesh.x(idx_x_max),md.mesh.x(idx_x_min),'stable');
%data_in_subset = intersect(data_in(idx_x_max),data_in(idx_x_min),'stable');
y_mesh_min = (percent_y_min * (max(md.mesh.y) - min(md.mesh.y)) / 100) + min(md.mesh.y);
y_mesh_max = (percent_y_max * (max(md.mesh.y) - min(md.mesh.y)) / 100) + min(md.mesh.y);
idx_y_max = md.mesh.y(:)<=y_mesh_max;
idx_y_min = md.mesh.y(:)>=y_mesh_min;
idx_Y_min = Y(:,1)>=y_mesh_min;
idx_Y_max = Y(:,1)<=y_mesh_max;
Y_min = find(idx_Y_min,1,'first');
Y_max = find(idx_Y_max,1,'last');
mesh_y_subset = intersect(md.mesh.y(idx_y_max),md.mesh.y(idx_y_min),'stable');
%data = InterpFromMeshToGrid(md.mesh.elements,md.mesh.x,md.mesh.y,[data_in],X(1,1:end),Y(1:end,1),NaN);
%mask = InterpFromMeshToGrid(md.mesh.elements,md.mesh.x,md.mesh.y,[mask_in],X(1,1:end),Y(1:end,1),NaN);


%data(mask>0) = NaN;


for i=1:3
   switch(i)
      case 1; md=md1; titlestr='t=0 yrs'; data_in = averaging(md,Kn,0); mask_in = md.results.TransientSolution(1).MaskIceLevelset;
      case 2; md=md1; titlestr='t=200 yrs (Eigencalving)'; data_in = averaging(md,Kn_ec,0); mask_in = md.results.TransientSolution(end).MaskIceLevelset;
      case 3; md=md2; titlestr='t=200 yrs (Von Mises)'; data_in = averaging(md,Kn_vm,0); mask_in = md.results.TransientSolution(end).MaskIceLevelset;
      otherwise error('not supported');
   end
   subplotmodel(3,1,i,options);
   %data_in = md.results.TransientSolution(end).Vel; %data_in(data_in==0)=NaN;
	%mask_in = md.results.TransientSolution(end).MaskIceLevelset; %data_in(data_in==0)=NaN;
   data = InterpFromMeshToGrid(md.mesh.elements,md.mesh.x,md.mesh.y,[data_in],X(1,1:end),Y(1:end,1),NaN);
   mask = InterpFromMeshToGrid(md.mesh.elements,md.mesh.x,md.mesh.y,[mask_in],X(1,1:end),Y(1:end,1),NaN);
   data(mask>0) = NaN;
   overlay(X(1,X_min:X_max),Y(Y_min:Y_max,1), data(Y_min:Y_max,X_min:X_max), 'caxis',[vmin vmax],'colormap',cm,'hemisphere',-1);
   hold on
   lastfront=isoline(md1, md1.results.TransientSolution(end).MaskIceLevelset,'output','matrix');
	lastfrontvm = isoline(md2, md2.results.TransientSolution(end).MaskIceLevelset,'output','matrix');
   if i==1
      h2=plot(lastfront(:,1),lastfront(:,2),'-','Color', colorpalette(1,:),'LineWidth',2.5);
      hold on
		h3=plot(lastfrontvm(:,1),lastfrontvm(:,2),'-','Color',colorpalette(2,:),'LineWidth',2.5);
		scaleruler(pairoptions('scaleruler',[min(md.mesh.x)+(max(md.mesh.x)-min(md.mesh.x))/21 max(md.mesh.y)-(max(md.mesh.y)-min(md.mesh.y))/9 sr (max(md.mesh.y)-min(md.mesh.y))/30 1],'scaleruler_fontsize',12,'FontColor','k'));
		legend('EC (200 yrs)','VM (200 yrs)','Location','northeast');
		hold off
		hreg = showregion(min(md.mesh.x),min(md.mesh.y),max(md.mesh.x),max(md.mesh.y),[.55 .05 .3 .3],'color',[1 0 0],'offset2',0*10^3,'LineWidth',1,'hemisphere',-1,'marker',1);
		axes(hreg);
	end
   hold off
end

cb = manualcb(vmin,vmax,cm,'Position',[.75 .3 .02 .6],'orientation','vertical');
