clc
clear all
close all

choice_surf_BC_schemes = {
  'twopoint_onesided'
  'threepoint_onesided'
  'ghostpoint'
  };
legend_strings = {
  'Two-point, one-sided'
  'Three-point, one-sided'
  'Ghost point'
  };

%% Read and plot own results

wa = 350;
ha = 300;

margins_hor = [90,25,25,25];
margins_ver = [50,90];

H = setup_multipanel_figure( wa, ha, margins_hor, margins_ver);

set( H.Ax{ 1,1},'xlim',[0,1],'ylim',[0,450],'xtick',[0,0.2,0.4,0.6,0.8,1]);
set( H.Ax{ 1,2},'xlim',[0,1],'ylim',[0,450],'xtick',[0,0.2,0.4,0.6,0.8,1]);
set( H.Ax{ 1,3},'xlim',[0,1],'ylim',[0,450],'xtick',[0,0.2,0.4,0.6,0.8,1]);

set( H.Ax{ 1,2},'yticklabels','');
set( H.Ax{ 1,3},'yticklabels','');

xlabel( H.Ax{ 1,1},'x / L');
xlabel( H.Ax{ 1,2},'x / L');
xlabel( H.Ax{ 1,3},'x / L');

ylabel( H.Ax{ 1,1},'u (m/yr)');

title( H.Ax{ 1,1},'A');
title( H.Ax{ 1,2},'B');
title( H.Ax{ 1,3},'C');

% Empty objects for legend - nz lines
nzs = [8,10,12,16,20,24,32];
nzs_legendstrings = {'nz = 8','nz = 10','nz = 12','nz = 16','nz = 20','nz = 24','nz = 32'};
linecolors = linspecer( length( nzs),'sequential');
for nzi = 1: length( nzs)
  line('parent',H.Ax{ 1,1},'xdata',[],'ydata',[],'color',linecolors( nzi,:),'linewidth',2);
end

% Zoom axeses
for a = 1:3
  x = (a-1)*wa+sum(margins_hor(1:a));
  y = margins_ver(2);
  % Zoom axes
  xs = x+wa*0.7;
  ys = y+ha*0.15;
  was = wa*0.3;
  has = ha*0.3;
  xlims = [0.2,0.3];
  ylims = [37 ,50 ];
  H.Axs( 1,a) = axes('parent',H.Fig,'units','pixels','position',[xs,ys,was,has],'box','on',...
    'fontsize',24,'xlim',xlims,'ylim',ylims,'xtick',xlims,'ytick',ylims);
  % zoom box
  xb = [xlims(1), xlims(2), xlims(2), xlims(1), xlims(1)];
  yb = [ylims(1), ylims(1), ylims(2), ylims(2), ylims(1)];
  line('parent',H.Ax{ 1,a},'xdata',xb,'ydata',yb);
  % zoom lines
  line('parent',H.Ax{ 1,a},'xdata',[xlims(1),(xs-x    )/wa],'ydata',[ylims(2),((ys-y+has)/ha)*250],'linestyle',':')
  line('parent',H.Ax{ 1,a},'xdata',[xlims(2),(xs-x+was)/wa],'ydata',[ylims(1),((ys-y    )/ha)*250],'linestyle',':')
end
  
for si = 1: length( choice_surf_BC_schemes)

  choice_surf_BC_scheme = choice_surf_BC_schemes{ si};

  for nzi = 1: length( nzs)

    nz = nzs( nzi);

    run_name = ['ISMIP_HOM_A_160_Glenreg_' choice_surf_BC_scheme '_nz' num2str( nz)];
    filename = [run_name '.mat'];
    load( filename);

    j = ceil( grid.ny / 4);

    x = 0.5 + grid.x / 160e3;
    u_surf_b = u( :,:,1);
    u_surf_a = map_cx_ax_2D( grid, map_cy_ay_2D( grid, u_surf_b));
    u_trans = u_surf_a( :,j);
    
    x = [x; 1+x];
    u_trans = [u_trans; u_trans];

    line('parent',H.Ax{  1,si},'xdata',x,'ydata',u_trans,'linewidth',2,'color',linecolors( nzi,:));
    line('parent',H.Axs( 1,si),'xdata',x,'ydata',u_trans,'linewidth',2,'color',linecolors( nzi,:));

  end
  
end

legend( H.Ax{ 1,1},nzs_legendstrings,'location','northwest');