clc
clear all
close all

wa = 300;
ha = 300;

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

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

for a = 1: 3
  set( H.Ax{ 1,a},'xlim',[8 * 0.8, 32 / 0.8],'xtick',[8,16,32],...
    'xscale','log','yscale','log','ylim',[0.01,1.0])
  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

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

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

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

%% Read data

% High-resolution solution
u_hi = load('ISMIP_HOM_C_5_Glen_ghostpoint_nz128.mat','u');
u_hi = max( max( u_hi.u( :,:,1)));

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

twopoint   = zeros( size( nzs));
threepoint = zeros( size( nzs));
ghostpoint = zeros( size( nzs));

for nzi = 1: length( nzs)
  
  nz = nzs( nzi);
  
  % Two-point
  filename = ['ISMIP_HOM_C_5_Glen_twopoint_onesided_nz' num2str( nz) '.mat'];
  u = load( filename,'u');
  twopoint( nzi) = abs( (max( max( u.u( :,:,1))) - u_hi) / u_hi);
  
  % Three-point
  filename = ['ISMIP_HOM_C_5_Glen_threepoint_onesided_nz' num2str( nz) '.mat'];
  u = load( filename,'u');
  threepoint( nzi) = abs( (max( max( u.u( :,:,1))) - u_hi) / u_hi);
  
  % Ghost point
  filename = ['ISMIP_HOM_C_5_Glen_ghostpoint_nz' num2str( nz) '.mat'];
  u = load( filename,'u');
  ghostpoint( nzi) = abs( (max( max( u.u( :,:,1))) - u_hi) / u_hi);
  
end

%% Plot results

% Raw results
line('parent',H.Ax{ 1,1},'xdata',nzs,'ydata',twopoint,'linestyle','none','marker','o','markerfacecolor','r',...
  'markeredgecolor','k','markersize',8)
line('parent',H.Ax{ 1,2},'xdata',nzs,'ydata',threepoint,'linestyle','none','marker','o','markerfacecolor','r',...
  'markeredgecolor','k','markersize',8)
line('parent',H.Ax{ 1,3},'xdata',nzs,'ydata',ghostpoint,'linestyle','none','marker','o','markerfacecolor','r',...
  'markeredgecolor','k','markersize',8)

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

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

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

% Plot loglinear fits
line('parent',H.Ax{ 1,1},'xdata', nzs,'ydata', f_twopoint,'color','b','linestyle','--','linewidth',2)
legend( H.Ax{ 1,1},['R = ' num2str( -p_twopoint( 1))])

line('parent',H.Ax{ 1,2},'xdata', nzs,'ydata', f_threepoint,'color','b','linestyle','--','linewidth',2)
legend( H.Ax{ 1,2},['R = ' num2str( -p_threepoint( 1))])

line('parent',H.Ax{ 1,3},'xdata', nzs,'ydata', f_ghostpoint,'color','b','linestyle','--','linewidth',2)
legend( H.Ax{ 1,3},['R = ' num2str( -p_ghostpoint( 1))])