clc
clear all
close all

wa = 300;
ha = 300;

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

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

for a = 1: 1
  set( H.Ax{ 1,a},'xlim',[8 * 0.8, 32 / 0.8],'xtick',[8,16,32],...
    'xscale','log','yscale','log','ylim',[0.0005,0.5],'ytick',[0.001,0.01,0.1])
  if a>1
    set( H.Ax{ 1,a},'yticklabels',[])
  end

  % Empty line objects for legend
  line('parent',H.Ax{ 1,a},'xdata',[],'ydata',[],'color','b','linestyle','--','linewidth',2);
end

xlabel( H.Ax{1,1},'n_z')

ylabel( H.Ax{1,1},'err u')

%% experiment A, L = 160 km

% Read data

% High-resolution solution
u_hi_twopoint   = load('ISMIP_HOM_A_160_linearvisc_twopoint_onesided_nz128.mat','u');
u_hi_threepoint = load('ISMIP_HOM_A_160_linearvisc_threepoint_onesided_nz128.mat','u');
u_hi_ghostpoint = load('ISMIP_HOM_A_160_linearvisc_ghostpoint_nz128.mat','u');
u_hi_twopoint   = max( max( u_hi_twopoint.u( :,:,1)));
u_hi_threepoint = max( max( u_hi_threepoint.u( :,:,1)));
u_hi_ghostpoint = max( max( u_hi_ghostpoint.u( :,:,1)));

nzs = [8,10,12,16,20,24,32];

A_160_twopoint   = zeros( size( nzs));
A_160_threepoint = zeros( size( nzs));
A_160_ghostpoint = zeros( size( nzs));

for nzi = 1: length( nzs)
  
  nz = nzs( nzi);
  
  % Two-point
  filename = ['ISMIP_HOM_A_160_linearvisc_twopoint_onesided_nz' num2str( nz) '.mat'];
  u = load( filename,'u');
  A_160_twopoint( nzi) = abs( (max( max( u.u( :,:,1))) - u_hi_twopoint) / u_hi_twopoint);
  
  % Three-point
  filename = ['ISMIP_HOM_A_160_linearvisc_threepoint_onesided_nz' num2str( nz) '.mat'];
  u = load( filename,'u');
  A_160_threepoint( nzi) = abs( (max( max( u.u( :,:,1))) - u_hi_threepoint) / u_hi_threepoint);
  
  % Ghost point
  filename = ['ISMIP_HOM_A_160_linearvisc_ghostpoint_nz' num2str( nz) '.mat'];
  u = load( filename,'u');
  A_160_ghostpoint( nzi) = abs( (max( max( u.u( :,:,1))) - u_hi_ghostpoint) / u_hi_ghostpoint);
  
end

% Plot results

colors = linspecer( 3);

% Empty line objects for legend
line('parent',H.Ax{ 1,1},'xdata', [],'ydata', [],'color',colors( 1,:),'linestyle','--','linewidth',2)
line('parent',H.Ax{ 1,1},'xdata', [],'ydata', [],'color',colors( 2,:),'linestyle','--','linewidth',2)
line('parent',H.Ax{ 1,1},'xdata', [],'ydata', [],'color',colors( 3,:),'linestyle','--','linewidth',2)

% Raw results
c = colors( 1,:);
line('parent',H.Ax{ 1,1},'xdata',nzs,'ydata',A_160_twopoint,'linestyle','none','marker','o','markerfacecolor',c,...
  'markeredgecolor',c,'markersize',8)
c = colors( 2,:);
line('parent',H.Ax{ 1,1},'xdata',nzs,'ydata',A_160_threepoint,'linestyle','none','marker','o','markerfacecolor',c,...
  'markeredgecolor',c,'markersize',8)
c = colors( 3,:);
line('parent',H.Ax{ 1,1},'xdata',nzs,'ydata',A_160_ghostpoint,'linestyle','none','marker','o','markerfacecolor',c,...
  'markeredgecolor',c,'markersize',8)

% Fit loglinear curves
p_twopoint = polyfit( log( nzs), log( A_160_twopoint), 1);
f_twopoint = exp( polyval( p_twopoint, log( nzs) ));

p_threepoint = polyfit( log( nzs), log( A_160_threepoint), 1);
f_threepoint = exp( polyval( p_threepoint, log( nzs) ));

p_ghostpoint = polyfit( log( nzs), log( A_160_ghostpoint), 1);
f_ghostpoint = exp( polyval( p_ghostpoint, log( nzs) ));

% Plot loglinear fits
legendstr = {};
c = colors( 1,:);
line('parent',H.Ax{ 1,1},'xdata', nzs,'ydata', f_twopoint,'color',c,'linestyle','--','linewidth',2)
legendstr{ end+1} = ['R = ' num2str( round( 100 * -p_twopoint( 1)) / 100 )];

c = colors( 2,:);
line('parent',H.Ax{ 1,1},'xdata', nzs,'ydata', f_threepoint,'color',c,'linestyle','--','linewidth',2)
legendstr{ end+1} = ['R = ' num2str( round( 100 * -p_threepoint( 1)) / 100 )];

c = colors( 3,:);
line('parent',H.Ax{ 1,1},'xdata', nzs,'ydata', f_ghostpoint,'color',c,'linestyle','--','linewidth',2)
legendstr{ end+1} = ['R = ' num2str( round( 100 * -p_ghostpoint( 1)) / 100 )];

legend( H.Ax{ 1,1}, legendstr,'location','southwest')
