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 ISMIP-HOM ensemble results
foldername ='D:\Documents\ProjectPlioceen\UU spullen\momentum_balance_numerics\ismip_all\';

model_types = {...
  'aas1','FS';...
  'aas2','FS';...
  'ahu1','HO';...
  'ahu2','HO';...
  'bds1','HO';...
  'cma1','FS';...
  'cma2','HO';...
  'dpo1','HO';...
  'fpa1','HO';...
  'fpa2','FS';...
  'fsa1','HO';...
  'ghg1','FS';...
  'jvj1','FS';...
  'lpe1','HO';...
  'mbr1','HO';...
  'mmr1','FS';...
  'mtk1','HO';...
  'oga1','FS';...
  'oso1','HO';...
  'rhi1','FS';...
  'rhi2','HO';...
  'rhi3','FS';...
  'rhi4','HO';...
  'rhi5','HO';...
  'spr1','FS';...
  'ssu1','FS';...
  'tpa1','HO';...
  'yko1','FS'};

models = dir(foldername);
models = models(3:end);

% ensemble.a160 = [];
% ensemble.a080 = [];
% ensemble.a040 = [];
% ensemble.a020 = [];
% ensemble.a010 = [];
% ensemble.a005 = [];
% 
% for mi = 1: length(models)
%   
%   modeldata = dir([foldername models(mi).name]);
%   modeldata = modeldata(3:end);
%   
%   % Go over all experiments, check if this model has them.
%   flds = fields(ensemble);
%   for xi = 1:length(flds)
%     ex = flds{xi};
%     
%     for di = 1:length(modeldata)
%       mdname = modeldata(di).name;
%       str = [ex '.txt'];
%       if length(mdname) >= length(str)
%         if strcmpi(mdname(end-length(str)+1:end),str)
%           % This is the experiment from this model
%           
%           disp(['Reading data from model ' models(mi).name ', experiment ' ex])
%           
%           fid = fopen([foldername models(mi).name '/' mdname]);
%           temp = textscan(fid,'%s','delimiter','\n','MultipleDelimsAsOne',1); temp = temp{1};
%           fclose(fid);
%           
%           n = length(temp);
%           nx = sqrt(n);
%           if nx-floor(nx)>0
%             error('whaa!')
%           end
%           x_vec = zeros(n,1);
%           y_vec = zeros(n,1);
%           u_vec = zeros(n,1);
%           v_vec = zeros(n,1);
%           w_vec = zeros(n,1);
%           for i = 1:n
%             temp2 = textscan(temp{i},'%f %f %f %f %f %f %f %f');
%             x_vec(i) = temp2{1};
%             y_vec(i) = temp2{2};
%             u_vec(i) = temp2{3};
%             v_vec(i) = temp2{4};
%             w_vec(i) = temp2{5};
%           end
%           
%           ensemble.(ex).(models(mi).name).x = reshape(x_vec,[nx,nx]);
%           ensemble.(ex).(models(mi).name).y = reshape(y_vec,[nx,nx]);
%           ensemble.(ex).(models(mi).name).u = reshape(u_vec,[nx,nx]);
%           ensemble.(ex).(models(mi).name).v = reshape(v_vec,[nx,nx]);
%           ensemble.(ex).(models(mi).name).w = reshape(w_vec,[nx,nx]);
%           
%           if (ensemble.(ex).(models(mi).name).x(1,1) == ensemble.(ex).(models(mi).name).x(end,1))
%             ensemble.(ex).(models(mi).name).x = ensemble.(ex).(models(mi).name).x';
%             ensemble.(ex).(models(mi).name).y = ensemble.(ex).(models(mi).name).y';
%             ensemble.(ex).(models(mi).name).u = ensemble.(ex).(models(mi).name).u';
%             ensemble.(ex).(models(mi).name).v = ensemble.(ex).(models(mi).name).v';
%             ensemble.(ex).(models(mi).name).w = ensemble.(ex).(models(mi).name).w';
%           end
%           
%         end
%       end
%     end
%   end
%   
% end
% 
% save('ISMIP_HOM_A_ensemble.mat','ensemble');
load('ISMIP_HOM_A_ensemble.mat')

%% Set up GUI

wa = 350;
ha = 300;

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

wf = wa*3 + sum( margins_hor);
hf = ha   + sum( margins_ver);

H.Fig = figure('color','w','position',[200,200,wf,hf]);
for a = 1: 3
  x = (a-1)*wa+sum(margins_hor(1:a));
  y = margins_ver(2);
  H.Ax( 1,a) = axes('parent',H.Fig,'units','pixels','position',[x,y,wa,ha],...
    'fontsize',20,'xlim',[0,1],'ylim',[0,125],'xtick',[0,0.2,0.4,0.6,0.8,1],'ytick',[0,25,50,75,100,125],...
    'xgrid','on','ygrid','on');
end

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 - ISMIP-HOM ensemble
c_FS  = [0.2,0.5,1.0];
c_HO  = [0.0,0.0,0.0];
% patch( 'parent',H.Ax( 1,2),'vertices',[],'faces',[],'facecolor',c_FS,'edgecolor','none','facealpha',0.7)
patch( 'parent',H.Ax( 1,2),'vertices',[],'faces',[],'facecolor',c_HO,'edgecolor','none','facealpha',0.3)
% line(  'parent',H.Ax( 1,2),'xdata',[],'ydata',[],'color',c_FS,'linewidth',3);
line(  'parent',H.Ax( 1,2),'xdata',[],'ydata',[],'color',c_HO,'linewidth',3,'linestyle','--');

% 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,3),'xdata',[],'ydata',[],'color',linecolors( nzi,:),'linewidth',3);
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.12;
  ys = y+ha*0.35;
  was = wa*0.3;
  has = ha*0.3;
  xlims = [0.2,0.3];
  ylims = [1  ,5  ];
  H.Axs( 1,a) = axes('parent',H.Fig,'units','pixels','position',[xs,ys,was,has],'box','on',...
    'fontsize',20,'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(1),((ys-y    )/ha)*125],'linestyle',':')
  line('parent',H.Ax( 1,a),'xdata',[xlims(2),(xs-x+was)/wa],'ydata',[ylims(1),((ys-y    )/ha)*125],'linestyle',':')
end

%% Plot results
  
for si = 1: length( choice_surf_BC_schemes)

  choice_surf_BC_scheme = choice_surf_BC_schemes{ si};

  %% Plot Pattyn2008 ensemble
  for a = 1

    if a == 1
      ex = 'a160';
    elseif a == 2
      ex = 'a080';
    elseif a == 3
      ex = 'a040';
    elseif a == 4
      ex = 'a020';
    elseif a == 5
      ex = 'a010';
    elseif a == 6
      ex = 'a005';
    end

    x_FS = linspace(0,1,200)';
    u_FS = [];
    u_HO = [];

    patch_HO  = patch( 'parent',H.Ax(  1,si),'xdata',[],'ydata',[],'facecolor',c_HO,'edgecolor','none','facealpha',0.3);
%     patch_FS  = patch( 'parent',H.Ax(  1,si),'xdata',[],'ydata',[],'facecolor',c_FS,'edgecolor','none','facealpha',0.7);
    line_HO   = line(  'parent',H.Ax(  1,si),'xdata',[],'ydata',[],'color',c_HO,'linewidth',3,'linestyle','--');
%     line_FS   = line(  'parent',H.Ax(  1,si),'xdata',[],'ydata',[],'color',c_FS,'linewidth',3);

    patch_HOs = patch( 'parent',H.Axs( 1,si),'xdata',[],'ydata',[],'facecolor',c_HO,'edgecolor','none','facealpha',0.3);
%     patch_FSs = patch( 'parent',H.Axs( 1,si),'xdata',[],'ydata',[],'facecolor',c_FS,'edgecolor','none','facealpha',0.7);
    line_HOs  = line(  'parent',H.Axs( 1,si),'xdata',[],'ydata',[],'color',c_HO,'linewidth',3,'linestyle','--');
%     line_FSs  = line(  'parent',H.Axs( 1,si),'xdata',[],'ydata',[],'color',c_FS,'linewidth',3);

    flds = fields(ensemble.(ex));
    for mi = 1:length(flds)
      m = flds{mi};

      x = ensemble.(ex).(m).x(:,1);
      y = ensemble.(ex).(m).y(1,:)';
      yi = round(0.25*length(y));
      u = ensemble.(ex).(m).u(:,yi);

      % Determine if this model is FS or HO
      FS = false;
      HO = false;
      for mii = 1:size(model_types,1)
        if strcmpi(model_types{mii,1},m)
          if strcmpi(model_types{mii,2},'FS')
            FS = true;
          else
            HO = true;
          end
        end
      end
      if ~(FS || HO)
        for mii = 1:size(model_types,1)
          if strcmpi(model_types{mii,1}(1:3),m(1:3))
            if strcmpi(model_types{mii,2},'FS')
              FS = true;
            else
              HO = true;
            end
          end
        end
      end
      if ~(FS || HO)
        % Unknown model?
        continue
      end

      % Add to data ranges for HO/FS models
      up = interp1(x,u,x_FS);
      if FS
        u_FS(:,end+1) = up;
      else
        u_HO(:,end+1) = up;
      end

    end

    % Missing data points
    m = true(size(x_FS));
    for i = 1:length(x_FS)
      if sum(isnan(u_FS(i,:)))+sum(isnan(u_HO(i,:)))>0
        m(i) = false;
      end
    end
    u_FS = u_FS(m,:);
    u_HO = u_HO(m,:);
    x_FS = x_FS(m);

    % ISMIP-HOM ensemble data
    uav_FS = mean(u_FS,2);
    uav_HO = mean(u_HO,2);
    sigma_FS = zeros(size(x_FS));
    sigma_HO = zeros(size(x_FS));
    for i = 1:size(u_FS,1)
      sigma_FS(i) = std(u_FS(i,:));
      sigma_HO(i) = std(u_HO(i,:));
      if isnan(sigma_FS(i)); sigma_FS(i) = 0; end
      if isnan(sigma_HO(i)); sigma_HO(i) = 0; end
    end
    umin_FS = uav_FS - sigma_FS;
    umax_FS = uav_FS + sigma_FS;
    umin_HO = uav_HO - sigma_HO;
    umax_HO = uav_HO + sigma_HO;

    xdata = [x_FS;flipud(x_FS)];
    ydata = [umin_FS;flipud(umax_FS)];
%     set(patch_FS ,'xdata',xdata,'ydata',ydata)
%     set(patch_FSs,'xdata',xdata,'ydata',ydata)

    xdata = [x_FS;flipud(x_FS)];
    ydata = [umin_HO;flipud(umax_HO)];
    set(patch_HO ,'xdata',xdata,'ydata',ydata)
    set(patch_HOs,'xdata',xdata,'ydata',ydata)

    set(line_HO ,'xdata',x_FS,'ydata',uav_HO)
    set(line_HOs,'xdata',x_FS,'ydata',uav_HO)
%     set(line_FS ,'xdata',x_FS,'ydata',uav_FS)
%     set(line_FSs,'xdata',x_FS,'ydata',uav_FS)

  end

  %% Read and plot own results

  for nzi = 1: length( nzs)

    nz = nzs( nzi);

    run_name = ['ISMIP_HOM_A_160_Glen_' 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),'Full-Stokes','Higher-Order','FS mean','HO mean','location','northwest')
legend( H.Ax( 1,2),'Pattyn2008 ensemble','Pattyn2008 mean','location','southeast')

pos = get( H.Ax( 1,3),'position');
legend( H.Ax( 1,3),nzs_legendstrings,'location','eastoutside');
set( H.Ax( 1,3),'position',pos);