%Paper_figures.m
%
%Charles Brunette - 2020/11/09
% Preparing the figures for the manuscript : A new parameterization for the
% free drift of sea ice (The Cryosphere)
%% Common
addpath('/Users/charlesbrunette/Dropbox/FREEDRIFT2020/Code/Matlab/') % -> MACBOOK
TablePath = '/Users/charlesbrunette/Dropbox/FREEDRIFT2020/Analysis/Freedrift/Data/';
FAITPaths_LOCAL_CB

Make_EASE_grid_CB
%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Figure 1 - Maps of spatial distribution of alpha side by side with piomas
% thickness climatology:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 1.1) Load the cell structure:
TablePath = '/Users/charlesbrunette/Dropbox/FREEDRIFT2020/Analysis/Freedrift/Data/';
load([TablePath 'XDATA_YDATA_cells.mat']);
% -> that cell structure is safe wrt ice concentration - i.e. IABP points
% where A = 0 have been filtered out...

Alpha_lf = NaN(361*361,1);
hi = Alpha_lf;
cacw = Alpha_lf; % also calculate the atmospheric to oceanic drag ratio
rhoa = 1.3; % air density
rhow = 1026; % Arctic water density 

tic
f = waitbar(0,'Browsing Cell structure');  
for i = 19680:361*361
    
    hi(i) = mean(XDATA_cell{i,4});
    
    UiC = complex(YDATA_cell{i,1},YDATA_cell{i,2});
    UaC = complex(XDATA_cell{i,1},XDATA_cell{i,2});
    
     if numel(YDATA_cell{i,1})>=10
        p = polyfit(UaC,UiC,1);
        [theta,rho]=cart2pol(real(p(1)),imag(p(1)));
        Alpha_lf(i) = rho;
        
        cacw(i) = (rho*1e-2)^2*rhow/rhoa; % from alpha = sqrt(rhoaCa/rhowCw)
     end
    
 waitbar(i/(361*361),f,['Browsing Cell structure: ' int2str(i) '/' int2str(361*361) ', TOC=' int2str(toc) 's']);     
    
end

%%
% 1.2) Make the side by side plot
Alpha361=reshape(Alpha_lf,[361,361]);
hi361=reshape(hi,[361,361]);
cacw361=reshape(cacw,[361,361]);

fig = figure(1);
set(fig,'units','normalized','outerposition',[0 0 1 1])

%tlo = tiledlayout(3,4,'TileSpacing','compact','Padding','compact'); % Create Tiles structure
tlo = tiledlayout(4,4,'TileSpacing','compact','Padding','compact'); % Create Tiles structure

nexttile([2 2]); % -> Alpha map

levels = 0:0.25:2.5;

ncpolarm('lat',66) % Create Arctic Mask
box on; hold on;
%pcolorm(EASElat,EASElon,Alpha361');
contourfm(EASElat,EASElon,Alpha361',levels)%,'LineStyle','none');
cb = contourcbar;
%cb.Title.String = '\alpha [%]';
cb.Location = 'westoutside';
%xlabel(cb,'\alpha [%]');
cb.Title.String = '\alpha [%]';
caxis([0 2.5]);                                              
title({'Wind-ice transfer coefficient'},'FontSize', 20);
textm(60,-135,'a)','Fontsize',25);

set(gca,'Fontsize',15);

n2 = nexttile([2 2]); % -> hi map

levels = 0:0.4:4;

ncpolarm('lat',66) % Create Arctic Mask
box on; hold on;
%pcolorm(EASElat,EASElon,hi361');
contourfm(EASElat,EASElon,hi361',levels)
cb2 = contourcbar;
%cb2.Title.String = 'h_i [m]';
cb2.Location = 'eastoutside';
%xlabel(cb2,'h_i [m]');
cb2.Title.String = 'h_i [m]';
caxis([0 4]);
title({'PIOMAS ice thickness'},'FontSize', 20);
colormap(n2,flipud(parula));
textm(60,-135,'b)','Fontsize',25);

set(gca,'Fontsize',15);

%set(tlo.Children,'XTick',[], 'YTick', []); % all in one

% 1.3) scatter plot + linfit of Alpha vs. h_i
% This justifies usinga linear relationship between alpha & u_i

iquick = find(~isnan(Alpha_lf) & ~isnan(hi));
p = deming(hi(iquick),Alpha_lf(iquick)); % deming = total least squares function
xh = [0:0.25:4]';

    %%%% Calculate goodness of fit
    % ypred = hi(iquick).*p(2) + p(1);
    % yres = Alpha_lf(iquick) - ypred;
    % SSresid = sum(yres.^2);
    % SStotal = (length(Alpha_lf(iquick))-1) * var(Alpha_lf(iquick));
    % rsq = 1 - SSresid/SStotal;
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%

%fig12 = figure(12); hold on;
%set(fig12,'units','normalized','outerposition',[0 0 1 1])

nexttile([1 4]) % plot with scatter of of local hi and Alpha coeff;

hold on;
s12 = scatter(hi(iquick),Alpha_lf(iquick),'.');
lf12 = plot([0:0.25:4]', p(1) + p(2).*xh,'-r','LineWidth',3);
ylim([0 3.5]); ylabel('\alpha [%]');
xlim([0 4]); xlabel(['h_i [m]']);
%title('Climatological local wind-ice transfer coefficient vs PIOMAS thickness');

%legend([lf12],{['\alpha = ' num2str(p(1),2) num2str(p(2),2) '*h']});
legend([lf12],{['\alpha = \alpha_{FD} - \beta h']});
text(0.02,0.98,'c)','Units', 'Normalized', 'VerticalAlignment', 'Top','Fontsize',25);

grid on; box on;

set(gca, 'Fontsize', 15);


nexttile([1 4]) % plot with the distribution of the Ca/Cw ratio 
hold on; grid on; box on;
h1 = histogram(cacw,'BinEdges',[0:0.01:1],'Normalization','PDF');
xlim([0 1])
xlabel('C_a/C_w');
ylabel('PDF');
%legend(h1,{'1979-2018'},'location','northeast'); 
text(0.02,0.98,'d)','Units', 'Normalized', 'VerticalAlignment', 'Top','Fontsize',25);
set(gca,'Fontsize',15);


%% Other version of Figure 1.2
% This one shows the ca/cw ratio - seeing if it compares to the classical
% estimate of ~0.2
%CAlculating the ca/cw ratio: (wioth rho_w = 1026 kg/m3 and rho_a = 1.3 kg/m3)
cacw = (Alpha361./100).^2*1026/1.3;

%                                 fig6 = figure(6);
% 
%                                 set(fig6,'units','normalized','outerposition',[0 0 .5 1])
%                                 tlo = tiledlayout(4,3,'TileSpacing','compact','Padding','compact');
% 
%                                 %%%% MAP OF CV
%                                 nexttile([3 3]);
% 
%                                 levels = 0:0.05:0.6;
% 
%                                 ncpolarm('lat',66) % Create Arctic Mask
%                                 box on; hold on;
%                                 %pcolorm(EASElat,EASElon,Alpha361');
%                                 contourfm(EASElat,EASElon,cacw',levels)%,'LineStyle','none');
%                                 cb = contourcbar;
%                                 cb.Location = 'westoutside';
%                                 xlabel(cb,'C_a/C_w');
%                                 caxis([0 0.6]);                                              
%                                 title({'C_a/C_w ratio'},'FontSize', 25);
%                                 textm(60,-135,'a)','Fontsize',25);
%                                 set(gca,'Fontsize',20);
% 
%                                 nexttile([1 3]); %%%%%%%% Distribution of ca/cw speed
%                                 hold on; grid on; box on;
%                                 h1 = histogram(cacw,'BinEdges',[0:0.01:1],'Normalization','PDF');
%                                 xlim([0 1])
%                                 xlabel('C_a/C_w');
%                                 ylabel('PDF');
%                                 text(0.02,0.98,'b)','Units', 'Normalized', 'VerticalAlignment', 'Top','Fontsize',25);
%                                 set(gca,'Fontsize',20);
                                
                                fig6 = figure(6);
                                iquick = find(~isnan(cacw) & ~isnan(hi361));
                                
                                hold on;
                                s12 = scatter(hi(iquick),cacw(iquick),'.');
                                %lf12 = plot([0:0.25:4]', p(1) + p(2).*xh,'-r','LineWidth',3);
                                ylim([0 3.5]); ylabel('C_a/C_w');
                                xlim([0 4]); xlabel(['h_i [m]']);
                                title('C_a/C_w vs PIOMAS thickness');
                              
                                grid on; box on;

                                set(gca, 'Fontsize', 20);
                                
                                

%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% FIGURE 2 & 3 - MAP OF U_w
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%After running ITERATIVE_SOLVER_alpha_uo_CB.m
% Created 3 Uo fields for different periods: uwall,vwall ; uwpr,vwpr ;
% uwpo,vwpro

%%%%%% CAREFUL!!!!! %%%%%
%Problem is solved in the EASE grid with 0° pointting right
%We need to take the transpose of the matlab EASE lat lon grids!!

uocn = NaN(361*361,1); vocn=uocn;
uocnpre = uocn; vocnpre=uocn;
uocnpo = uocn; vocnpo=uocn;

uiclimo = NaN(361*361,1); viclimo = uiclimo;

%%%%% ALL YEARS:
uocn(Ind_isdata) = UO(:,end); %uo(ArcticOceanMask4IceEdge ~=1)=NaN;
vocn(Ind_isdata) = VO(:,end); %vo(ArcticOceanMask4IceEdge ~=1)=NaN;
uocn=reshape(uocn,[361,361]);
vocn=reshape(vocn,[361,361]);
%Rotate uo vo so they point along LAT/LON lines
Uo_ease = uocn + 1i*vocn;
Uo_ll_cms = Uo_ease.*exp(1i*deg2rad(EASElon'));
uomag = sqrt(UO(:,end).^2+VO(:,end).^2);
gridmag = sqrt(uocn.^2+vocn.^2);

    % Compare to mean U_i_buoy over the same period to see  how similar it
    % is...
    load([TablePath 'XDATA_YDATA_cells.mat']); Ind_isdata = find(Datapoint_numbers>0);
    ui = YDATA_cell(:,1); vi = YDATA_cell(:,2);
    uiclimo = cellfun(@(x) mean(x, 'all'), ui);
    viclimo = cellfun(@(x) mean(x, 'all'), vi);
    uiclimmag = sqrt(uiclimo(Ind_isdata).^2+viclimo(Ind_isdata).^2);
    viclimo=reshape(viclimo,[361,361]);
    uiclimo=reshape(uiclimo,[361,361]);
    Ui_ease = uiclimo + 1i*viclimo;
    Ui_ll_cms = Ui_ease.*exp(1i*deg2rad(EASElon'));
    gridmagbuoy = sqrt(uiclimo.^2+viclimo.^2);
        
    %Useful couple of lines:
    %ratioUoUbuoy = gridmag./gridmagbuoy; ratioUoUbuoy = ratioUoUbuoy(:); ratioUoUbuoy(isnan(ratioUoUbuoy))=[];
    %Ypc = prctile(ratioUoUbuoy,95);
    %meanUoUbratio = mean(ratioUoUbuoy(ratioUoUbuoy<=Ypc));
    
%%%% PRE-2000:
 load([TablePath 'XDATA_YDATA_cells_Pre2000.mat']); Ind_isdataPre2000 = find(Datapoint_numbers>0); 
 
uocnpre(Ind_isdataPre2000) = uwpr; %uo(ArcticOceanMask4IceEdge ~=1)=NaN;
vocnpre(Ind_isdataPre2000) = vwpr; %vo(ArcticOceanMask4IceEdge ~=1)=NaN;
uocnpre=reshape(uocnpre,[361,361]);
vocnpre=reshape(vocnpre,[361,361]);
%Rotate uo vo so they point along LAT/LON lines
Uo_easepre = uocnpre + 1i*vocnpre;
Uo_ll_cmspre = Uo_easepre.*exp(1i*deg2rad(EASElon'));
uomagpre = sqrt(uwpr.^2+vwpr.^2);
gridmagpre = sqrt(uocnpre.^2+vocnpre.^2);

    %Calculating the gridded magnitude of Ua of h_i for the different time
    %periods:
    uiv = YDATA_cell(:,1); viv =YDATA_cell(:,2);
    uispeed = cellfun(@(x1,x2) mean(sqrt(x1.^2+x2.^2), 'all'), uiv,viv);
    uav = XDATA_cell(:,1); vav =XDATA_cell(:,2); 
    htot = XDATA_cell(:,4);
    uaspeed = cellfun(@(x1,x2) mean(sqrt(x1.^2+x2.^2), 'all'), uav,vav);
    htotal = cellfun(@(x) mean(x, 'all'), htot);
    uigridpre = NaN(361*361,1); uigridpre(Ind_isdataPre2000)= uispeed(Ind_isdataPre2000); uigridpre=reshape(uigridpre,[361,361]);
    uagridpre = NaN(361*361,1); uagridpre(Ind_isdataPre2000)= uaspeed(Ind_isdataPre2000); uagridpre=reshape(uagridpre,[361,361]);
    hgridpre = NaN(361*361,1); hgridpre(Ind_isdataPre2000)= htotal(Ind_isdataPre2000);  hgridpre=reshape( hgridpre,[361,361]);
    
%%%%% Post-2000:
load([TablePath 'XDATA_YDATA_cells_Post2000.mat']); Ind_isdataPost2000 = find(Datapoint_numbers>0); 

uocnpo(Ind_isdataPost2000) = uwpo; %uo(ArcticOceanMask4IceEdge ~=1)=NaN;
vocnpo(Ind_isdataPost2000) = vwpo; %vo(ArcticOceanMask4IceEdge ~=1)=NaN;
uocnpo=reshape(uocnpo,[361,361]);
vocnpo=reshape(vocnpo,[361,361]);
%Rotate uo vo so they point along LAT/LON lines
Uo_easepo = uocnpo + 1i*vocnpo;
Uo_ll_cmspo = Uo_easepo.*exp(1i*deg2rad(EASElon'));
uomagpo = sqrt(uwpo.^2+vwpo.^2);
gridmagpo = sqrt(uocnpo.^2+vocnpo.^2);

    %Calculating the gridded magnitude of Ua of h_i for the different time
    %periods:
    uiv = YDATA_cell(:,1); viv =YDATA_cell(:,2);
    uispeed = cellfun(@(x1,x2) mean(sqrt(x1.^2+x2.^2), 'all'), uiv,viv);
    uav = XDATA_cell(:,1); vav =XDATA_cell(:,2); 
    htot = XDATA_cell(:,4);
    uaspeed = cellfun(@(x1,x2) mean(sqrt(x1.^2+x2.^2), 'all'), uav,vav);
    htotal = cellfun(@(x) mean(x, 'all'), htot);
    uigridpo = NaN(361*361,1); uigridpo(Ind_isdataPost2000)= uispeed(Ind_isdataPost2000); uigridpo=reshape(uigridpo,[361,361]);
    uagridpo = NaN(361*361,1); uagridpo(Ind_isdataPost2000)= uaspeed(Ind_isdataPost2000); uagridpo=reshape(uagridpo,[361,361]);
    hgridpo = NaN(361*361,1); hgridpo(Ind_isdataPost2000)= htotal(Ind_isdataPost2000);  hgridpo=reshape( hgridpo,[361,361]);

%% Figure itself:

fig = figure(2);   %Make figure handle
set(fig,'units','normalized','outerposition',[0 0 1 1])

tlo = tiledlayout(5,11,'TileSpacing','normal','Padding','normal');

nexttile([5 5]); %%%%%%%% LARGE FIGURE WITH CLIMO CURRENT

ncpolarm('lat',66,'grid','label') % Create Arctic Mask
hold on;

%%%%%%%%% QUIVERMC %%%%%%%%%%%%%%
%'density': % of data points plotted on the map
%'reference': scaling of the vectors %2
[~,cb] = quivermc(EASElat',EASElon',real(Uo_ll_cms),-imag(Uo_ll_cms),'density',40,'reference', 2, 'units','cm/s',...
                                                         'colormap', parula, 'colorbar','southoutside',...
                                                         'linewidth',1,'Fontsize',15);                                                                                                      
xlabel(cb,'Current speed [cm/s]','Fontsize',25)
cb.Ticks = [0:1:10];
caxis([0 10]);
cb.TickLabels = {'0','1','2','3','4','5','6','7','8','9','10'};

%gridm on;                                              
title({'Surface ocean current estimates'},'FontSize', 15);
textm(60,-135,'a)','Fontsize',25);
box on;
set(gca,'Fontsize',25);


% nexttile([3 3]); %%%%%%%% SMALL FIGURE PRE2000
% ncpolarm('lat',66,'grid','label')
% hold on;
% qm2 = quivermc(EASElat',EASElon',real(Uo_ll_cmspre),-imag(Uo_ll_cmspre),'density',40,'reference', 2, 'units','cm/s','colormap', parula,'linewidth',1,'Fontsize',15);
% 
% %xlabel(cb2,'Current speed [cm/s]','Fontsize',15)
% %cb2.Ticks = [0:1:10];
% %caxis([0 10]);
% %cb2.TickLabels = {'0','1','2','3','4','5','6','7','8','9','10'};
% title({'Pre-2000'});
% textm(60,-135,'b)','Fontsize',25);
% box on;
% set(gca,'Fontsize',20);
% 
% 
% nexttile([3 3]); %%%%%%%% SMALL FIGURE POST2000
% ncpolarm('lat',66,'grid','label');
% hold on;
% quivermc(EASElat',EASElon',real(Uo_ll_cmspo),-imag(Uo_ll_cmspo),'density',40,'reference', 2, 'units','cm/s','colormap', parula,'linewidth',1,'Fontsize',15);
% 
% %xlabel(cb3,'Current speed [cm/s]','Fontsize',15)
% %cb3.Ticks = [0:1:10];
% %caxis([0 10]);
% %cb3.TickLabels = {'0','1','2','3','4','5','6','7','8','9','10'};
% title({'Post-2000'});
% textm(60,-135,'c)','Fontsize',25);
% box on;
% set(gca,'Fontsize',20);

%Same figure, but instead of the vector plots for pre / post 2000; showing
%the anomaly with respect ot the whole period:

        n2 =nexttile([3 3]); %%%%%%%% SMALL FIGURE PRE2000
        ncpolarm('lat',66)
        hold on;
        gb = pcolorm(EASElat',EASElon',gridmagpre-gridmag);
        cmap = redblue;
        set(n2,'colormap', cmap);
        caxis([-3 3]);
        title({'Pre-2000 anomaly'});
        textm(60,-135,'b)','Fontsize',25);
        box on;
        set(gca,'Fontsize',20);
        

        n3 = nexttile([3 3]); %%%%%%%% SMALL FIGURE POST2000
        ncpolarm('lat',66);
        hold on;
        gb2 = pcolorm(EASElat',EASElon',gridmagpo-gridmag);
        cba = colorbar('eastoutside');
        set(n3,'colormap', cmap);
        caxis([-3 3]);
        xlabel(cba,'Speed anomaly [cm/s]','Fontsize',25)
        title({'Post-2000 anomaly'});
        textm(60,-135,'c)','Fontsize',25);
        box on;
        set(gca,'Fontsize',20);


nexttile([2 6]); %%%%%%%% Histogramm with the 3:
hold on; grid on; box on;
 %h1 = histogram(uomag,'BinEdges',[0:0.25:30],'Normalization','PDF');
  h2 = histogram(uomagpre,'BinEdges',[0:0.1:8],'Normalization','PDF');
   h3 = histogram(uomagpo,'BinEdges',[0:0.1:8],'Normalization','PDF');

xlim([0 8])
xlabel('Current speed [cm/s]');
ylabel('PDF');
%legend([h1 h2 h3],{'1979-2018','<2000', '>2000'},'location','northeast'); 
legend([h2 h3],{'<2000', '>2000'},'location','northeast'); 
text(0.02,0.98,'d)','Units', 'Normalized', 'VerticalAlignment', 'Top','Fontsize',25);
set(gca,'Fontsize',20);
set(gca,'YAxisLocation','right')

%% Other version of Figure 2
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%% JUST CLIMO U_w by itself:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fig21 = figure(21);   %Make figure handle
set(fig21,'units','normalized','outerposition',[0 0 .5 1])
tlo = tiledlayout(4,3,'TileSpacing','compact','Padding','compact');

nexttile([3 3]); %%%%%%%% LARGE FIGURE WITH CLIMO CURRENT

%ncpolarm('lat',66,'grid','label') % Create Arctic Mask
ncpolarm('lat',66) % -> colormap doesn't like the ncpolarm grid....
hold on;

%%%%%%%%% QUIVERMC %%%%%%%%%%%%%%
%'density': % of data points plotted on the map

%'reference': scaling of the vectors %2
% First approach: Vector plots with length + color indicating magnitude of
% vector
%[~,cb] = quivermc(EASElat',EASElon',real(Uo_ll_cms),-imag(Uo_ll_cms),'density',40,'reference', 1, 'units','cm/s',...
%                                                         'colormap', parula, 'colorbar','westoutside',...
%                                                         'linewidth',1,'Fontsize',15);

%Second approach: A pcorlormap for the magnitude of the vector, then all
%vectors have the same scale
% Colormap sequence
f21col = pcolorm(EASElat',EASElon',gridmag);
cb = colorbar('Westoutside');
colormap(fig21,'parula')
% Vector sequence:
quivermc(EASElat',EASElon',real(Uo_ll_cms)./gridmag,-imag(Uo_ll_cms)./gridmag,'density',30,'reference', 1,'linewidth',1,'Fontsize',15);

xlabel(cb,'U_w [cm/s]','Fontsize',25)
cb.Ticks = [0:1:5]; % vector version was going to 10 cm/s
caxis([0 5]);
cb.TickLabels = {'0','1','2','3','4','5'};

%gridm on;                                              
title({'Under-ice surface ocean currents'},'FontSize', 15);
textm(60,-135,'a)','Fontsize',25);
box on;
set(gca,'Fontsize',25);

nexttile([1 3]); %%%%%%%% Distribution of current speed
hold on; grid on; box on;
h1 = histogram(uomag,'BinEdges',[0:0.25:30],'Normalization','PDF');
xlim([0 10])
xlabel('U_w [cm/s]');
ylabel('PDF');
legend(h1,{'1979-2019'},'location','northeast'); 
text(0.02,0.98,'b)','Units', 'Normalized', 'VerticalAlignment', 'Top','Fontsize',25);
set(gca,'Fontsize',20);
%set(gca,'YAxisLocation','right')
                               
%% 
%%%%% COMPARING Uw CLIMO TO Ui CLIMO!!!!!
        fig21 = figure(21);
        tlo = tiledlayout(1,3,'TileSpacing','normal','Padding','normal');

        nexttile %%% 1) Uw
        ncpolarm('lat',66); hold on;
        %colormap:
        f21col = pcolorm(EASElat',EASElon',gridmag);
        cb = colorbar('Westoutside');
        colormap(fig21,'parula')
        % Vector sequence:
        quivermc(EASElat',EASElon',real(Uo_ll_cms)./gridmag,-imag(Uo_ll_cms)./gridmag,'density',30,'reference', 1,'linewidth',1,'Fontsize',15);
        xlabel(cb,'Current speed [cm/s]','Fontsize',25)
        cb.Ticks = [0:1:5]; % vector version was going to 10 cm/s
        caxis([0 5]);
        cb.TickLabels = {'0','1','2','3','4','5'};                                             
        title({'Under-ice surface ocean current estimates'},'FontSize', 15);
        textm(60,-135,'a)','Fontsize',25);
        box on;
        set(gca,'Fontsize',20);

        nexttile %%% 2) Ui_buoy
        ncpolarm('lat',66); hold on;
        f21col = pcolorm(EASElat',EASElon',gridmag);
        cb = colorbar('Westoutside');
        colormap(fig21,'parula')
        % Vector sequence:
        quivermc(EASElat',EASElon',real(Ui_ll_cms)./gridmagbuoy,-imag(Ui_ll_cms)./gridmagbuoy,'density',30,'reference', 1,'linewidth',1,'Fontsize',15);
        xlabel(cb,'Buoy drift speed [cm/s]','Fontsize',25)
        cb.Ticks = [0:1:5]; % vector version was going to 10 cm/s
        caxis([0 5]);
        cb.TickLabels = {'0','1','2','3','4','5'};                                             
        title({'Climo Buoy drift speed'},'FontSize', 15);
        textm(60,-135,'b)','Fontsize',25);
        box on;
        set(gca,'Fontsize',20);

        n3 = nexttile; %%% 3)Difference
        ncpolarm('lat',66);
        hold on;
        gb3 = pcolorm(EASElat',EASElon',gridmagbuoy-gridmag);
        cb3 = colorbar('westoutside');
        set(n3,'colormap', cmap);
        cb3.Ticks = [-5:1:5];
        caxis([-5 5]);
        cb3.TickLabels = {'-5','-4','-3','-2','-1','0','1','2','3','4','5'};
        xlabel(cb3,'\Delta U_w [cm/s]','Fontsize',25)
        title({'Speed difference'});
        textm(60,-135,'c)','Fontsize',25);
        box on;
        set(gca,'Fontsize',20);

%% %%%%%%%%%%%%%%
%%%%%% Figure showing in depth Pre/Post 2000 difference
%%%%%%%%%%%%%%%%

cmap = redblue; % Creating a red (+) to blue (-) colormap; 

fig22 = figure(22);   %Make figure handle
tlo = tiledlayout(14,8,'TileSpacing','none','Padding','compact');

nexttile([4 4]); %%%%%%%% Pre-2000 U_w

% First approach - vectors with color:
% ncpolarm('lat',66,'grid','label')
% hold on;
% [~,cb1] = quivermc(EASElat',EASElon',real(Uo_ll_cmspre),-imag(Uo_ll_cmspre),'density',40,'reference', 2, 'units','cm/s','colormap', parula, 'colorbar','westoutside','linewidth',1,'Fontsize',15);
% xlabel(cb1,'Current speed [cm/s]','Fontsize',15)
% cb1.Ticks = [0:1:10];
% caxis([0 10]);
% cb1.TickLabels = {'0','1','2','3','4','5','6','7','8','9','10'};
% title({'Currents: Pre-2000'});
% textm(60,-135,'a)','Fontsize',25);
% box on;
% set(gca,'Fontsize',15);

%Second approach - Colormap with unit scale vectors:
    ncpolarm('lat',66); hold on;
    f22acol = pcolorm(EASElat',EASElon',gridmagpre);
    cb = colorbar('Westoutside');
    colormap(gcf,'parula')
    % Vector sequence:
    quivermc(EASElat',EASElon',real(Uo_ll_cmspre)./gridmagpre,-imag(Uo_ll_cmspre)./gridmagpre,'density',30,'reference', 1,'linewidth',1,'Fontsize',15);
    xlabel(cb,'U_w [cm/s]','Fontsize',15)
    cb.Ticks = [0:1:5]; % vector version was going to 10 cm/s
    caxis([0 5]);
    cb.TickLabels = {'0','1','2','3','4','5'};

    %gridm on;                                              
    title({'Ocean current estimates, Pre-2000'},'FontSize', 15);
    textm(60,-135,'a)','Fontsize',25);
    box on;
    set(gca,'Fontsize',15);

nexttile([4 4]); %%%%%%%% Post-2000 U_w

% ncpolarm('lat',66,'grid','label')
% hold on;
% [~,cb1] = quivermc(EASElat',EASElon',real(Uo_ll_cmspo),-imag(Uo_ll_cmspo),'density',40,'reference', 2, 'units','cm/s','colormap', parula, 'colorbar','eastoutside','linewidth',1,'Fontsize',15);
% xlabel(cb1,'Current speed [cm/s]','Fontsize',15)
% cb1.Ticks = [0:1:10];
% caxis([0 10]);
% cb1.TickLabels = {'0','1','2','3','4','5','6','7','8','9','10'};
% title({'Currents: Post-2000'});
% textm(60,-135,'b)','Fontsize',25);
% box on;
% set(gca,'Fontsize',15);

%Second approach - Colormap with unit scale vectors:
    ncpolarm('lat',66); hold on;
    f22bcol = pcolorm(EASElat',EASElon',gridmagpo);
    cb = colorbar('Eastoutside');
    colormap(gcf,'parula')
    % Vector sequence:
    quivermc(EASElat',EASElon',real(Uo_ll_cmspo)./gridmagpo,-imag(Uo_ll_cmspo)./gridmagpo,'density',30,'reference', 1,'linewidth',1,'Fontsize',15);
    xlabel(cb,'U_w [cm/s]','Fontsize',15)
    cb.Ticks = [0:1:5]; % vector version was going to 10 cm/s
    caxis([0 5]);
    cb.TickLabels = {'0','1','2','3','4','5'};
    %gridm on;                                              
    title({'Ocean current estimates, Post-2000'},'FontSize', 15);
    textm(60,-135,'b)','Fontsize',25);
    box on;
    set(gca,'Fontsize',15);

n3 = nexttile([4 4]); %%%%%%%% DELTA U_w

ncpolarm('lat',66);
hold on;
gb3 = pcolorm(EASElat',EASElon',gridmagpo-gridmagpre);
cb3 = colorbar('westoutside');
set(n3,'colormap', cmap);
cb3.Ticks = [-5:1:5];
caxis([-5 5]);
cb3.TickLabels = {'-5','-4','-3','-2','-1','0','1','2','3','4','5'};
xlabel(cb3,'\Delta U_w [cm/s]','Fontsize',15)
title({'Current speed difference'});
textm(60,-135,'c)','Fontsize',25);
box on;
set(gca,'Fontsize',15);

n4 = nexttile([4 4]); %%%%%%%% DELTA U_i_buoy

ncpolarm('lat',66);
hold on;
gb4 = pcolorm(EASElat',EASElon',uigridpo-uigridpre);
cb4 = colorbar('eastoutside');
set(n4,'colormap', cmap);
cb4.Ticks = [-10:2:10];
caxis([-10 10]);
cb4.TickLabels = {'-10','-8','-6','-4','-2','0','2','4','6','8','10'};
xlabel(cb4,'\Delta U_{buoy} [cm/s]','Fontsize',15)
title({'Buoy speed difference'});
textm(60,-135,'d)','Fontsize',25);
box on;
set(gca,'Fontsize',15);

n5 = nexttile([4 4]); %%%%%%%% DELTA hi

ncpolarm('lat',66);
hold on;
gb5 = pcolorm(EASElat',EASElon', hgridpo- hgridpre);
cb5 = colorbar('westoutside');
set(n5,'colormap', cmap);
cb5.Ticks = [-2:0.5:2];
caxis([-2 2]);
cb5.TickLabels = {'-2','-1.5','-1','-0.5','0','0.5','1','1.5','2'};
xlabel(cb5,'\Delta h_i  [m]','Fontsize',15)
title({'Thickness difference'});
textm(60,-135,'e)','Fontsize',25);
box on;
set(gca,'Fontsize',15);

n6 = nexttile([4 4]); %%%%%%%% DELTA U_a

ncpolarm('lat',66);
hold on;
gb6 = pcolorm(EASElat',EASElon',uagridpo-uagridpre);
cb6 = colorbar('eastoutside');
set(n6,'colormap', cmap);
cb6.Ticks = [-5:1:5];
caxis([-5 5]);
cb6.TickLabels = {'-5','-4','-3','-2','-1','0','1','2','3','4','5'};
xlabel(cb6,'\Delta U_a [m/s]','Fontsize',15)
title({'Wind speed difference'});
textm(60,-135,'f)','Fontsize',25);
box on;
set(gca,'Fontsize',15);

n7 = nexttile([2 8]); %%%%%%%% U_w distribution
hold on; grid on; box on;
 h2 = histogram(uomagpre,'BinEdges',[0:0.1:10],'Normalization','PDF','FaceColor','#0000FF','FaceAlpha',0.6);
 h3 = histogram(uomagpo,'BinEdges',[0:0.1:10],'Normalization','PDF','FaceColor','#FF0000','FaceAlpha',0.6);

xlim([0 10])
xlabel('U_w [cm/s]');
ylabel('PDF'); 
legend([h2 h3],{'<2000', '>2000'},'location','northeast'); 
text(0.02,0.98,'g)','Units', 'Normalized', 'VerticalAlignment', 'Top','Fontsize',25);
set(gca,'Fontsize',15);


%% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% FIGURE SHOWING THE SUMMER/WINTER DIFFERENCE
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

        uocnpre = NaN(361*361,1); vocnpre=NaN(361*361,1);
        uocnpo = NaN(361*361,1); vocnpo=NaN(361*361,1);

                %%%% SUMMMER:
                 load([TablePath 'SUMMER_JASOND.mat']);
                 load([TablePath 'XDATA_YDATA_cells_SUMMER.mat']); Ind_summer = find(Datapoint_numbers>0); 

                uocnpre(Ind_summer) = UO(:,end); %uo(ArcticOceanMask4IceEdge ~=1)=NaN;
                vocnpre(Ind_summer) = VO(:,end); %vo(ArcticOceanMask4IceEdge ~=1)=NaN;
                uocnpre=reshape(uocnpre,[361,361]);
                vocnpre=reshape(vocnpre,[361,361]);
                %Rotate uo vo so they point along LAT/LON lines
                Uo_easepre = uocnpre + 1i*vocnpre;
                Uo_ll_cmspre = Uo_easepre.*exp(1i*deg2rad(EASElon'));
                uomagpre = sqrt(UO(:,end).^2+VO(:,end).^2);
                gridmagpre = sqrt(uocnpre.^2+vocnpre.^2);

                    %Calculating the gridded magnitude of Ua of h_i for the different time
                    %periods:
                    uiv = YDATA_cell(:,1); viv =YDATA_cell(:,2);
                    uispeed = cellfun(@(x1,x2) mean(sqrt(x1.^2+x2.^2), 'all'), uiv,viv);
                    uav = XDATA_cell(:,1); vav =XDATA_cell(:,2); 
                    htot = XDATA_cell(:,4);
                    uaspeed = cellfun(@(x1,x2) mean(sqrt(x1.^2+x2.^2), 'all'), uav,vav);
                    htotal = cellfun(@(x) mean(x, 'all'), htot);
                    uigridpre = NaN(361*361,1); uigridpre(Ind_summer)= uispeed(Ind_summer); uigridpre=reshape(uigridpre,[361,361]);
                    uagridpre = NaN(361*361,1); uagridpre(Ind_summer)= uaspeed(Ind_summer); uagridpre=reshape(uagridpre,[361,361]);
                    hgridpre = NaN(361*361,1); hgridpre(Ind_summer)= htotal(Ind_summer);  hgridpre=reshape(hgridpre,[361,361]);

                %%%%% WINTER:
                load([TablePath 'WINTER_JFMAMJ.mat']);
                 load([TablePath 'XDATA_YDATA_cells_WINTER.mat']); Ind_winter = find(Datapoint_numbers>0); 

                uocnpo(Ind_winter) = UO(:,end); %uo(ArcticOceanMask4IceEdge ~=1)=NaN;
                vocnpo(Ind_winter) = VO(:,end); %vo(ArcticOceanMask4IceEdge ~=1)=NaN;
                uocnpo=reshape(uocnpo,[361,361]);
                vocnpo=reshape(vocnpo,[361,361]);
                %Rotate uo vo so they point along LAT/LON lines
                Uo_easepo = uocnpo + 1i*vocnpo;
                Uo_ll_cmspo = Uo_easepo.*exp(1i*deg2rad(EASElon'));
                uomagpo = sqrt(UO(:,end).^2+VO(:,end).^2);
                gridmagpo = sqrt(uocnpo.^2+vocnpo.^2);

                    %Calculating the gridded magnitude of Ua of h_i for the different time
                    %periods:
                    uiv = YDATA_cell(:,1); viv =YDATA_cell(:,2);
                    uispeed = cellfun(@(x1,x2) mean(sqrt(x1.^2+x2.^2), 'all'), uiv,viv);
                    uav = XDATA_cell(:,1); vav =XDATA_cell(:,2); 
                    htot = XDATA_cell(:,4);
                    uaspeed = cellfun(@(x1,x2) mean(sqrt(x1.^2+x2.^2), 'all'), uav,vav);
                    htotal = cellfun(@(x) mean(x, 'all'), htot);
                    uigridpo = NaN(361*361,1); uigridpo(Ind_winter)= uispeed(Ind_winter); uigridpo=reshape(uigridpo,[361,361]);
                    uagridpo = NaN(361*361,1); uagridpo(Ind_winter)= uaspeed(Ind_winter); uagridpo=reshape(uagridpo,[361,361]);
                    hgridpo = NaN(361*361,1); hgridpo(Ind_winter)= htotal(Ind_winter);  hgridpo=reshape( hgridpo,[361,361]);


        %%%%%%% FIGURE SUMMER - WINTER !!!!!! %%%%%%%%%%%%%%

                 cmap = redblue; % Creating a red (+) to blue (-) colormap; 

                fig22 = figure(22);   %Make figure handle
                tlo = tiledlayout(14,8,'TileSpacing','none','Padding','compact');

                nexttile([4 4]); %%%%%%%% WINTER Uw

                %- Colormap with unit scale vectors:
                    ncpolarm('lat',66); hold on;
                    f22acol = pcolorm(EASElat',EASElon',gridmagpo);
                    cb = colorbar('Westoutside');
                    colormap(gcf,'parula')
                    % Vector sequence:
                    quivermc(EASElat',EASElon',real(Uo_ll_cmspo)./gridmagpo,-imag(Uo_ll_cmspo)./gridmagpo,'density',30,'reference', 1,'linewidth',1,'Fontsize',15);
                    xlabel(cb,'U_w [cm/s]','Fontsize',15)
                    cb.Ticks = [0:1:5]; % vector version was going to 10 cm/s
                    caxis([0 5]);
                    cb.TickLabels = {'0','1','2','3','4','5'};

                    %gridm on;                                              
                    title({'Uw estimates, Winter (JFMAMJ)'},'FontSize', 15);
                    textm(60,-135,'a)','Fontsize',25);
                    box on;
                    set(gca,'Fontsize',15);

                nexttile([4 4]); %%%%%%%% SUMMER Uw

                %Second approach - Colormap with unit scale vectors:
                    ncpolarm('lat',66); hold on;
                    f22bcol = pcolorm(EASElat',EASElon',gridmagpre);
                    cb = colorbar('Eastoutside');
                    colormap(gcf,'parula')
                    % Vector sequence:
                    quivermc(EASElat',EASElon',real(Uo_ll_cmspre)./gridmagpre,-imag(Uo_ll_cmspre)./gridmagpre,'density',30,'reference', 1,'linewidth',1,'Fontsize',15);
                    xlabel(cb,'U_w [cm/s]','Fontsize',15)
                    cb.Ticks = [0:1:5]; % vector version was going to 10 cm/s
                    caxis([0 5]);
                    cb.TickLabels = {'0','1','2','3','4','5'};
                    %gridm on;                                              
                    title({'Uw estimates, Summer (JASOND)'},'FontSize', 15);
                    textm(60,-135,'b)','Fontsize',25);
                    box on;
                    set(gca,'Fontsize',15);

                n3 = nexttile([4 4]); %%%%%%%% DELTA U_w

                ncpolarm('lat',66);
                hold on;
                gb3 = pcolorm(EASElat',EASElon',gridmagpre-gridmagpo);
                cb3 = colorbar('westoutside');
                set(n3,'colormap', cmap);
                cb3.Ticks = [-5:1:5];
                caxis([-5 5]);
                cb3.TickLabels = {'-5','-4','-3','-2','-1','0','1','2','3','4','5'};
                xlabel(cb3,'\Delta U_w [cm/s]','Fontsize',15)
                title({'Current speed difference'});
                textm(60,-135,'c)','Fontsize',25);
                box on;
                set(gca,'Fontsize',15);

                n4 = nexttile([4 4]); %%%%%%%% DELTA U_i_buoy

                ncpolarm('lat',66);
                hold on;
                gb4 = pcolorm(EASElat',EASElon',uigridpre-uigridpo);
                cb4 = colorbar('eastoutside');
                set(n4,'colormap', cmap);
                cb4.Ticks = [-10:2:10];
                caxis([-10 10]);
                cb4.TickLabels = {'-10','-8','-6','-4','-2','0','2','4','6','8','10'};
                xlabel(cb4,'\Delta U_{buoy} [cm/s]','Fontsize',15)
                title({'Buoy speed difference'});
                textm(60,-135,'d)','Fontsize',25);
                box on;
                set(gca,'Fontsize',15);

                n5 = nexttile([4 4]); %%%%%%%% DELTA hi

                ncpolarm('lat',66);
                hold on;
                gb5 = pcolorm(EASElat',EASElon', hgridpre- hgridpo);
                cb5 = colorbar('westoutside');
                set(n5,'colormap', cmap);
                cb5.Ticks = [-2:0.5:2];
                caxis([-2 2]);
                cb5.TickLabels = {'-2','-1.5','-1','-0.5','0','0.5','1','1.5','2'};
                xlabel(cb5,'\Delta h_i  [m]','Fontsize',15)
                title({'Thickness difference'});
                textm(60,-135,'e)','Fontsize',25);
                box on;
                set(gca,'Fontsize',15);

                n6 = nexttile([4 4]); %%%%%%%% DELTA U_a

                ncpolarm('lat',66);
                hold on;
                gb6 = pcolorm(EASElat',EASElon',uagridpre-uagridpo);
                cb6 = colorbar('eastoutside');
                set(n6,'colormap', cmap);
                cb6.Ticks = [-5:1:5];
                caxis([-5 5]);
                cb6.TickLabels = {'-5','-4','-3','-2','-1','0','1','2','3','4','5'};
                xlabel(cb6,'\Delta U_a [m/s]','Fontsize',15)
                title({'Wind speed difference'});
                textm(60,-135,'f)','Fontsize',25);
                box on;
                set(gca,'Fontsize',15);

                n7 = nexttile([2 8]); %%%%%%%% U_w distribution
                hold on; grid on; box on;
                 h2 = histogram(uomagpo,'BinEdges',[0:0.1:10],'Normalization','PDF','FaceColor','#0000FF','FaceAlpha',0.6);
                 h3 = histogram(uomagpre,'BinEdges',[0:0.1:10],'Normalization','PDF','FaceColor','#FF0000','FaceAlpha',0.6);

                xlim([0 10])
                xlabel('U_w [cm/s]');
                ylabel('PDF'); 
                legend([h2 h3],{'Summer', 'Winter'},'location','northeast'); 
                text(0.02,0.98,'g)','Units', 'Normalized', 'VerticalAlignment', 'Top','Fontsize',25);
                set(gca,'Fontsize',15);
                
                
                
%% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% FIGURE COMPARING Uo estimates TO DIFFERENT REANALYSIS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

                                % 1) ECCO
                                load('/Users/charlesbrunette/Dropbox/FREEDRIFT2020/Analysis/Freedrift/Data/ECCO_EASE_1992_2017.mat');
   
                                uECCO = uocnEASE';
                                vECCO = vocnEASE';
   
                                %Rotate uo vo so they point along LAT/LON lines
                                UoECCO = uECCO + 1i*vECCO;
                                Uo_ll_cmsECCO = UoECCO.*exp(1i*deg2rad(EASElon'));
                                gridmagECCO = sqrt(uECCO.^2+vECCO.^2);

                               %1) Estimates over ECCO period
                                load([TablePath 'XDATA_YDATA_cells_19922017.mat']); Ind_isdata = find(Datapoint_numbers>0);
                                load('/Users/charlesbrunette/Dropbox/FREEDRIFT2020/Analysis/Freedrift/Data/1992_2017_ThicknessRun.mat')
                                uocn = NaN(361*361,1); vocn=NaN(361*361,1);
                                uocn(Ind_isdata) = UO(:,end); %uo(ArcticOceanMask4IceEdge ~=1)=NaN;
                                vocn(Ind_isdata) = VO(:,end); %vo(ArcticOceanMask4IceEdge ~=1)=NaN;
                                uocn=reshape(uocn,[361,361]);
                                vocn=reshape(vocn,[361,361]);
                                %Rotate uo vo so they point along LAT/LON lines
                                Uo_ease = uocn + 1i*vocn;
                                Uo_ll_cms = Uo_ease.*exp(1i*deg2rad(EASElon'));
                                uomag = sqrt(UO(:,end).^2+VO(:,end).^2);
                                gridmag = sqrt(uocn.^2+vocn.^2);
                                
                                % 2) GLORYS
                                load('/Users/charlesbrunette/Dropbox/FREEDRIFT2020/Analysis/Freedrift/Data/GLORYS_EASE_1993_2018.mat');
   
                                uGLORYS = uocnEASE';
                                vGLORYS = vocnEASE';
   
                                %Rotate uo vo so they point along LAT/LON lines
                                UoGLORYS = uGLORYS + 1i*vGLORYS;
                                Uo_ll_cmsGLORYS = UoGLORYS.*exp(1i*deg2rad(EASElon'));
                                gridmagGLORYS = sqrt(uGLORYS.^2+vGLORYS.^2);

                               %2) Estimates over ECCO period
                                load([TablePath 'XDATA_YDATA_cells_19932018.mat']); Ind_isdata = find(Datapoint_numbers>0);
                                load('/Users/charlesbrunette/Dropbox/FREEDRIFT2020/Analysis/Freedrift/Data/1993_2018_ThicknessRun.mat')
                                uocng = NaN(361*361,1); vocng=NaN(361*361,1);
                                uocng(Ind_isdata) = UO(:,end); %uo(ArcticOceanMask4IceEdge ~=1)=NaN;
                                vocng(Ind_isdata) = VO(:,end); %vo(ArcticOceanMask4IceEdge ~=1)=NaN;
                                uocng=reshape(uocng,[361,361]);
                                vocng=reshape(vocng,[361,361]);
                                %Rotate uo vo so they point along LAT/LON lines
                                Uo_easeg = uocng + 1i*vocng;
                                Uo_ll_cmsg = Uo_easeg.*exp(1i*deg2rad(EASElon'));
                                uomagg = sqrt(UO(:,end).^2+VO(:,end).^2);
                                gridmagg = sqrt(uocng.^2+vocng.^2);
                                
                                
                                 % 3) CPOM/DOT
                                load('/Users/charlesbrunette/Dropbox/FREEDRIFT2020/Analysis/Freedrift/Data/CPOMDOT_EASE_2003_2014.mat');
   
                                ucpom = uocnEASE';
                                vcpom = vocnEASE';
   
                                %Rotate uo vo so they point along LAT/LON lines
                                Uocpom = ucpom + 1i*vcpom;
                                Uo_ll_cmscpom = Uocpom.*exp(1i*deg2rad(EASElon'));
                                gridmagcpom = sqrt(ucpom.^2+vcpom.^2);

                               %2) Estimates over CPOM period
                                load([TablePath 'XDATA_YDATA_cells_20032014.mat']); Ind_isdata = find(Datapoint_numbers>0);
                                load('/Users/charlesbrunette/Dropbox/FREEDRIFT2020/Analysis/Freedrift/Data/2003_2014_ThicknessRun.mat')
                                uocnc = NaN(361*361,1); vocnc=NaN(361*361,1);
                                uocnc(Ind_isdata) = UO(:,end); %uo(ArcticOceanMask4IceEdge ~=1)=NaN;
                                vocnc(Ind_isdata) = VO(:,end); %vo(ArcticOceanMask4IceEdge ~=1)=NaN;
                                uocnc=reshape(uocnc,[361,361]);
                                vocnc=reshape(vocnc,[361,361]);
                                %Rotate uo vo so they point along LAT/LON lines
                                Uo_easec = uocnc + 1i*vocnc;
                                Uo_ll_cmsc = Uo_easec.*exp(1i*deg2rad(EASElon'));
                                uomagc = sqrt(UO(:,end).^2+VO(:,end).^2);
                                gridmagc = sqrt(uocnc.^2+vocnc.^2);
                                
                                
                                
                                % 3) PIOMAS
                                load('/Users/charlesbrunette/Dropbox/FREEDRIFT2020/Analysis/Freedrift/Data/PIOMAS_UO_EASE_1979_2015.mat');
   
                                ucpio = uocnEASE';
                                vcpio = vocnEASE';
   
                                %Rotate uo vo so they point along LAT/LON lines
                                Uocpio = ucpio + 1i*vcpio;
                                Uo_ll_cmscpio = Uocpio.*exp(1i*deg2rad(EASElon'));
                                gridmagcpio = sqrt(ucpio.^2+vcpio.^2);
                                     iz = find(gridmagcpio==0); % Getting rid of zeroes on landmasses
                                     gridmagcpio(iz)=NaN;
                                        
                                        
                                %2) Estimates over PIOMAS period
                                load([TablePath 'XDATA_YDATA_cells_1979_2015.mat']); Ind_isdata = find(Datapoint_numbers>0);
                                load('/Users/charlesbrunette/Dropbox/FREEDRIFT2020/Analysis/Freedrift/Data/1979_2015_ThicknessRun.mat')
                                uocnpio = NaN(361*361,1); vocnpio=NaN(361*361,1);
                                uocnpio(Ind_isdata) = UO(:,end); %uo(ArcticOceanMask4IceEdge ~=1)=NaN;
                                vocnpio(Ind_isdata) = VO(:,end); %vo(ArcticOceanMask4IceEdge ~=1)=NaN;
                                uocnpio=reshape(uocnpio,[361,361]);
                                vocnpio=reshape(vocnpio,[361,361]);
                                %Rotate uo vo so they point along LAT/LON lines
                                Uo_easepio = uocnpio + 1i*vocnpio;
                                Uo_ll_cmspio = Uo_easepio.*exp(1i*deg2rad(EASElon'));
                                uomagpio = sqrt(UO(:,end).^2+VO(:,end).^2);
                                gridmpio = sqrt(uocnpio.^2+vocnpio.^2);       
%%

                        %%%%%%% FIGURE Ocean current comparison %%%%%%%%%%%%%%

                                cmap = redblue; % Creating a red (+) to blue (-) colormap; 

                                fig23 = figure(23);   %Make figure handle
                                tlo = tiledlayout(3,4,'TileSpacing','compact','Padding','compact');

                                nexttile; %%%%%%%% ECCO

                                %- Colormap with unit scale vectors:
                                    ncpolarm('lat',66); hold on;
                                    f23acol = pcolorm(EASElat',EASElon',gridmagECCO);
                                    cb = colorbar('Westoutside');
                                    colormap(gcf,'parula')
                                    % Vector sequence:
                                    quivermc(EASElat',EASElon',real(Uo_ll_cmsECCO)./gridmagECCO,-imag(Uo_ll_cmsECCO)./gridmagECCO,'density',30,'reference', 1,'linewidth',1,'Fontsize',15);
                                    xlabel(cb,'U_w [cm/s]','Fontsize',15)
                                    cb.Ticks = [0:1:5]; % vector version was going to 10 cm/s
                                    caxis([0 5]);
                                    cb.TickLabels = {'0','1','2','3','4','5'};

                                    %gridm on;                                              
                                    title({'ECCO (1992-2017)'},'FontSize', 15);
                                    textm(60,-135,'a)','Fontsize',25);
                                    box on;
                                    set(gca,'Fontsize',15);
                                    
                                    
                                 
                                n4 = nexttile; %%%%%%%% GLORYS

                                %- Colormap with unit scale vectors:
                                    ncpolarm('lat',66); hold on;
                                    f23dcol = pcolorm(EASElat',EASElon',gridmagGLORYS);
                                    %cb = colorbar('Southoutside');
                                    colormap(n4,'parula')
                                    % Vector sequence:
                                    quivermc(EASElat',EASElon',real(Uo_ll_cmsGLORYS)./gridmagGLORYS,-imag(Uo_ll_cmsGLORYS)./gridmagGLORYS,'density',30,'reference', 1,'linewidth',1,'Fontsize',15);
                                    %xlabel(cb,'U_w [cm/s]','Fontsize',15)
                                    %cb.Ticks = [0:1:5]; % vector version was going to 10 cm/s
                                    caxis([0 5]);
                                    %cb.TickLabels = {'0','1','2','3','4','5'};

                                    %gridm on;                                              
                                    title({'GLORYS (1993-2018)'},'FontSize', 15);
                                    textm(60,-135,'b)','Fontsize',25);
                                    box on;
                                    set(gca,'Fontsize',15);
                                    
                                    
                                   n7 = nexttile; %%%%%%%% CPOM/DOT

                                %- Colormap with unit scale vectors:
                                    ncpolarm('lat',66); hold on;
                                    f23gcol = pcolorm(EASElat',EASElon',gridmagcpom);
                                    %cb = colorbar('Southoutside');
                                    colormap(n7,'parula')
                                    % Vector sequence:
                                    quivermc(EASElat',EASElon',real(Uo_ll_cmscpom)./gridmagcpom,-imag(Uo_ll_cmscpom)./gridmagcpom,'density',30,'reference', 1,'linewidth',1,'Fontsize',15);
                                    %xlabel(cb,'U_w [cm/s]','Fontsize',15)
                                    %cb.Ticks = [0:1:5]; % vector version was going to 10 cm/s
                                    caxis([0 5]);
                                    %cb.TickLabels = {'0','1','2','3','4','5'};

                                    %gridm on;                                              
                                    title({'CPOM/DOT (2003-2014)'},'FontSize', 15);
                                    textm(60,-135,'c)','Fontsize',25);
                                    box on;
                                    set(gca,'Fontsize',15); 
                                    
                                    
                                    
                                n10 = nexttile; %%%%%%%% PIOMAS

                                %- Colormap with unit scale vectors:
                                    ncpolarm('lat',66); hold on;
                                    f23jcol = pcolorm(EASElat',EASElon',gridmagcpio);
                                    %cb = colorbar('Southoutside');
                                    colormap(n10,'parula')
                                    % Vector sequence:
                                    quivermc(EASElat',EASElon',real(Uo_ll_cmscpio)./gridmagcpio,-imag(Uo_ll_cmscpio)./gridmagcpio,'density',30,'reference', 1,'linewidth',1,'Fontsize',15);
                                    %xlabel(cb,'U_w [cm/s]','Fontsize',15)
                                    %cb.Ticks = [0:1:5]; % vector version was going to 10 cm/s
                                    caxis([0 5]);
                                    %cb.TickLabels = {'0','1','2','3','4','5'};

                                    %gridm on;                                              
                                    title({'PIOMAS (1979-2015)'},'FontSize', 15);
                                    textm(60,-135,'d)','Fontsize',25);
                                    box on;
                                    set(gca,'Fontsize',15);    
                                    
                                    
                                 nexttile; %%%%%%%% BUOY DERIVED OVER ECCO PERIOD

                                %- Colormap with unit scale vectors:
                                    ncpolarm('lat',66); hold on;
                                    f23bcol = pcolorm(EASElat',EASElon',gridmag);
                                    cb = colorbar('Westoutside');
                                    colormap(gcf,'parula')
                                    % Vector sequence:
                                    quivermc(EASElat',EASElon',real(Uo_ll_cms)./gridmag,-imag(Uo_ll_cms)./gridmag,'density',30,'reference', 1,'linewidth',1,'Fontsize',15);
                                    xlabel(cb,'U_w [cm/s]','Fontsize',15)
                                    cb.Ticks = [0:1:5]; % vector version was going to 10 cm/s
                                    caxis([0 5]);
                                    cb.TickLabels = {'0','1','2','3','4','5'};

                                    %gridm on;                                              
                                    title({'Clim (1992-2017)'},'FontSize', 15);
                                    textm(60,-135,'e)','Fontsize',25);
                                    box on;
                                    set(gca,'Fontsize',15);
                                    
                                    
                                    n5 = nexttile; %%%%%%%% BUOY DERIVED OVER GLORYS PERIOD

                                %- Colormap with unit scale vectors:
                                    ncpolarm('lat',66); hold on;
                                    f23ecol = pcolorm(EASElat',EASElon',gridmagg);
                                    %cb = colorbar('Southoutside');
                                    colormap(n5,'parula')
                                    % Vector sequence:
                                    quivermc(EASElat',EASElon',real(Uo_ll_cmsg)./gridmagg,-imag(Uo_ll_cmsg)./gridmagg,'density',30,'reference', 1,'linewidth',1,'Fontsize',15);
                                    %xlabel(cb,'U_w [cm/s]','Fontsize',15)
                                    %cb.Ticks = [0:1:5]; % vector version was going to 10 cm/s
                                    caxis([0 5]);
                                    %cb.TickLabels = {'0','1','2','3','4','5'};

                                    %gridm on;                                              
                                    title({'Clim (1993-2018)'},'FontSize', 15);
                                    textm(60,-135,'f)','Fontsize',25);
                                    box on;
                                    set(gca,'Fontsize',15);  
                                    
                                    
                                     n8 = nexttile; %%%%%%%% BUOY DERIVED OVER DOT PERIOD

                                %- Colormap with unit scale vectors:
                                    ncpolarm('lat',66); hold on;
                                    f23hcol = pcolorm(EASElat',EASElon',gridmagc);
                                    %cb = colorbar('Southoutside');
                                    colormap(n8,'parula')
                                    % Vector sequence:
                                    quivermc(EASElat',EASElon',real(Uo_ll_cmsc)./gridmagc,-imag(Uo_ll_cmsc)./gridmagc,'density',30,'reference', 1,'linewidth',1,'Fontsize',15);
                                    %xlabel(cb,'U_w [cm/s]','Fontsize',15)
                                    %cb.Ticks = [0:1:5]; % vector version was going to 10 cm/s
                                    caxis([0 5]);
                                    %cb.TickLabels = {'0','1','2','3','4','5'};

                                    %gridm on;                                              
                                    title({'Clim (2003-2014)'},'FontSize', 15);
                                    textm(60,-135,'g)','Fontsize',25);
                                    box on;
                                    set(gca,'Fontsize',15);
                                    
                                  nx =   nexttile; %%%%%%%% BUOY DERIVED OVER PIOMAS period

                                %- Colormap with unit scale vectors:
                                    ncpolarm('lat',66); hold on;
                                    f23hcol = pcolorm(EASElat',EASElon',gridmpio);
                                    %cb = colorbar('Southoutside');
                                    colormap(nx,'parula')
                                    % Vector sequence:
                                    quivermc(EASElat',EASElon',real(Uo_ll_cmspio)./gridmpio,-imag(Uo_ll_cmspio)./gridmpio,'density',30,'reference', 1,'linewidth',1,'Fontsize',15);
                                    %xlabel(cb,'U_w [cm/s]','Fontsize',15)
                                    %cb.Ticks = [0:1:5]; % vector version was going to 10 cm/s
                                    caxis([0 5]);
                                    %cb.TickLabels = {'0','1','2','3','4','5'};

                                    %gridm on;                                              
                                    title({'Clim (1979-2015)'},'FontSize', 15);
                                    textm(60,-135,'h)','Fontsize',25);
                                    box on;
                                    set(gca,'Fontsize',15);
                                       
                                    
                                n3 = nexttile; %%%%%%%% DELTA Estimates - ECCO

                                ncpolarm('lat',66);
                                hold on;
                                gb3 = pcolorm(EASElat',EASElon',gridmag-gridmagECCO);
                                cb3 = colorbar('Westoutside');
                                set(n3,'colormap', cmap);
                                cb3.Ticks = [-5:1:5];
                                caxis([-5 5]);
                                cb3.TickLabels = {'-5','-4','-3','-2','-1','0','1','2','3','4','5'};
                                xlabel(cb3,'\Delta U_w [cm/s]','Fontsize',15)
                                title({'Clim - ECCO'});
                                textm(60,-135,'i)','Fontsize',25);
                                box on;
                                set(gca,'Fontsize',15);
  
                                    
                              n6 = nexttile; %%%%%%%% DELTA Estimates - GLORYS

                                ncpolarm('lat',66);
                                hold on;
                                gb3 = pcolorm(EASElat',EASElon',gridmagg-gridmagGLORYS);
                                %cb3 = colorbar('Southoutside');
                                set(n6,'colormap', cmap);
                                %cb3.Ticks = [-5:1:5];
                                caxis([-5 5]);
                                %cb3.TickLabels = {'-5','-4','-3','-2','-1','0','1','2','3','4','5'};
                                %xlabel(cb3,'\Delta U_w [cm/s]','Fontsize',15)
                                title({'Clim - GLORYS'});
                                textm(60,-135,'j)','Fontsize',25);
                                box on;
                                set(gca,'Fontsize',15);
                                

                                n9 = nexttile; %%%%%%%% DELTA Estimates - CPOM/DOT

                                ncpolarm('lat',66);
                                hold on;
                                gb3 = pcolorm(EASElat',EASElon',gridmagc-gridmagcpom);
                                %cb3 = colorbar('Southoutside');
                                set(n9,'colormap', cmap);
                                %cb3.Ticks = [-5:1:5];
                                caxis([-5 5]);
                                %cb3.TickLabels = {'-5','-4','-3','-2','-1','0','1','2','3','4','5'};
                                %xlabel(cb3,'\Delta U_w [cm/s]','Fontsize',15)
                                title({'Clim - CPOM/DOT'});
                                textm(60,-135,'k)','Fontsize',25);
                                box on;
                                set(gca,'Fontsize',15);
                                
                                     
                          n12 = nexttile;%%%%%%%% DELTA Estimates - PIOMAS

                                ncpolarm('lat',66);
                                hold on;
                                gb3 = pcolorm(EASElat',EASElon',gridmpio-gridmagcpio);
                                %cb3 = colorbar('Southoutside');
                                set(n12,'colormap', cmap);
                                %cb3.Ticks = [-5:1:5];
                                caxis([-5 5]);
                                %cb3.TickLabels = {'-5','-4','-3','-2','-1','0','1','2','3','4','5'};
                                %xlabel(cb3,'\Delta U_w [cm/s]','Fontsize',15)
                                title({'Clim - PIOMAS'});
                                textm(60,-135,'l)','Fontsize',25);
                                box on;
                                set(gca,'Fontsize',15);

                                               

%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% FIGURE 4 - Error figures
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%After running ITERATIVE_SOLVER_alpha_uo_CB.m (Once for each of the parameterization wanted)
%Figures for 3 different takes at the error: Interannual / Seasonal / Spatial:

% First, calculate interannual variability and time series of speed and
% error:

    ubuoy = DATA.u_buoy(Ind_common);  %%%%% -> Ind_common is important here! Ind_common is the estimation of free drift at all buoy points where we have Uw estimates
    vbuoy = DATA.v_buoy(Ind_common);
        speedbuoy = sqrt(ubuoy.^2 + vbuoy.^2);
        angbuoy = atan2d(vbuoy,ubuoy);

    ufd = Ui_fit;
    vfd = Vi_fit;
        speedfd = sqrt(ufd.^2 + vfd.^2);
        angfd = atan2d(vfd,ufd);
        dthetafd = angfd-angbuoy;
        u180 = find(dthetafd<180); dthetafd(u180)= dthetafd(u180)+360; % -> flipping angles over (o) or under (u) 180 degree in the [-180 180] range
        o180 = find(dthetafd>180); dthetafd(o180)= dthetafd(o180)-360;
    
    %%%%%% CHECK PPF U/V direction - might be an error    
    uppf = DATA.u_PPFv4Wind(Ind_common);
    vppf = DATA.v_PPFv4Wind(Ind_common);
        speedppf = sqrt(uppf.^2 + vppf.^2);
        angppf = atan2d(vppf,uppf);
        dthetappf = angppf-angbuoy;
        u180 = find(dthetappf<180); dthetappf(u180)= dthetappf(u180)+360; % -> flipping angles over (o) or under (u) 180 degree in the [-180 180] range
        o180 = find(dthetappf>180); dthetappf(o180)= dthetappf(o180)-360;
        
% Alternative curve: e.g.: constant coefficient
    ufd2 = Ui_fit2;
    vfd2 = Vi_fit2;
        speedfd2 = sqrt(ufd2.^2 + vfd2.^2);        
        angfd2 = atan2d(vfd2,ufd2);
        dthetafd2 = angfd2-angbuoy;
        u180 = find(dthetafd2<180); dthetafd2(u180)= dthetafd2(u180)+360; % -> flipping angles over (o) or under (u) 180 degree in the [-180 180] range
        o180 = find(dthetafd2>180); dthetafd2(o180)= dthetafd2(o180)-360;
        
         ua = DATA.u_ERA5(Ind_common);
        va = DATA.v_ERA5(Ind_common);
        speeda = sqrt(ua.^2 + va.^2);
 
    year =  DATA.year(Ind_common);
    month = DATA.month(Ind_common);
    
    %%% For Figure 5 - Delta U/V:
    % H-param
    dufd = ufd - ubuoy;
    dvfd = vfd - vbuoy;
    % CST-param
    dufd2 = ufd2 - ubuoy;
    dvfd2 = vfd2 - vbuoy;
     % PPF
    duppf = uppf - ubuoy;
    dvppf = vppf - vbuoy;
    

    % 12 month seasonal cycle of drift speed
    CycleTable = NaN(12,8); 
    % [MeanSpeedBuoy MeanSpeedFD STDFD MeanSpeedFD2 STDFD2 MeanSpeedPPF STDPPF]

    for i = 1:12
      im = find(month==i)  ;
      CycleTable(i,1)= nanmean(speedbuoy(im));
      CycleTable(i,2)= nanmean(speedfd(im));
      CycleTable(i,3)= sqrt(nanmean((speedfd(im)-speedbuoy(im)).^2));
      CycleTable(i,4)= nanmean(speedfd2(im));
      CycleTable(i,5)= sqrt(nanmean((speedfd2(im)-speedbuoy(im)).^2));
      CycleTable(i,6)= nanmean(speedppf(im));
      CycleTable(i,7)= sqrt(nanmean((speedppf(im)-speedbuoy(im)).^2)); 
      CycleTable(i,8)= nanmean(speeda(im));
    end

    TSTable = NaN(max(year)-min(year)+1,7); % TimeSeriesTAlbe: TSTable
    % [MeanSpeedBuoy MeanSpeedFD STDFD MeanSpeedFD2 STDFD2 MeanSpeedPPF STDPPF]

    for i = 1:max(year)-min(year)+1
      iy = find(year==1978+i)  ;
      TSTable(i,1)= nanmean(speedbuoy(iy));
      TSTable(i,2)= nanmean(speedfd(iy));
      TSTable(i,3)= sqrt(nanmean((speedfd(iy)-speedbuoy(iy)).^2));
      TSTable(i,4)= nanmean(speedfd2(iy));
      TSTable(i,5)= sqrt(nanmean((speedfd2(iy)-speedbuoy(iy)).^2));
      TSTable(i,6)= nanmean(speedppf(iy));
      TSTable(i,7)= sqrt(nanmean((speedppf(iy)-speedbuoy(iy)).^2));
      TSTable(i,8)= nanmean(speeda(iy));
    end

%% Figures 4,5,6 

%Figure 4 - Panel with seasonality and time series, signal + coefficient of
%variation:
lw = 5; %-> line width
xmonth = [1:12]';
xyear = [1979:2019]';

fig41 = figure(41);   %Make figure handle
tlo = tiledlayout(2,4,'TileSpacing','normal','Padding','normal');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
n1 = nexttile([1 2]); % Seasonal cycle magnitude
hold on; box on;
ahb = plot(xmonth, CycleTable(:,1),'k','LineWidth',lw); % BUOY
ahppf = plot(xmonth, CycleTable(:,6),'-k','LineWidth',lw); ahppf.Color(4) = 0.5; % Polarpathfinder 
ahfd = plot(xmonth, CycleTable(:,2),'b','LineWidth',lw); % Free drift h
ahfd2 = plot(xmonth, CycleTable(:,4),'-.b','LineWidth',lw); % Free drift cst

grid on;
xlim([1 12]);
ylim([2 12]);
ylabel('|U_i| [cm/s]');
set(gca,'FontSize',25);

xticks(xmonth);
xticklabels({'Jan','Feb','Mar','Apr','May','Jun','Jul','Aug','Sep','Oct',...
             'Nov','Dec'});
xlabel('Month');
yticks([2:1:12]');
yticklabels({'2','','4','','6','','8','','10','','12'});

legend([ahb ahppf ahfd ahfd2],{'IABP buoy', 'PPF wind', 'FD thickness','FD constant'},'location','southeast','FontSize',15);

title('Sea ice drift seasonality')   
text(0.02,0.98,'a)','Units', 'Normalized', 'VerticalAlignment', 'Top','Fontsize',25);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
n2 = nexttile([1 2]); % Time series of annual average
hold on; box on;
bhb = plot(xyear, TSTable(:,1),'k','LineWidth',lw); % BUOY
bhppf = plot(xyear, TSTable(:,6),'-k','LineWidth',lw); bhppf.Color(4) = 0.5; % Polarpathfinder 
bhfd = plot(xyear, TSTable(:,2),'b','LineWidth',lw); % Free drift h
bhfd2 = plot(xyear, TSTable(:,4),'-.b','LineWidth',lw); % Free drift cst

grid on;
xlim([1979 2019]);
ylim([2 12]);
ylabel('|U_i| [cm/s]');
set(gca,'FontSize',25);

xticks(1979:1:2019);
xticklabels({'','','1980','','','','','1985','','','','',...
                   '1990','','','','','1995','','','','',...
                   '2000','','','','','2005','','','','',...
                   '2010','','','','','2015','','','',''});
xlabel('Year');
yticks([2:1:12]');
yticklabels({'2','','4','','6','','8','','10','','12'});

legend([bhb bhppf bhfd bhfd2],{'IABP buoy', 'PPF wind', 'FD thickness','FD constant'},'location','southeast','FontSize',15);

title('Annual average drift')   
text(0.02,0.98,'b)','Units', 'Normalized', 'VerticalAlignment', 'Top','Fontsize',25);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
n3 = nexttile([1 2]); % Seasonal cycle coefficient of variation
hold on; box on;
%ahb = plot(xmonth, CycleTable(:,1),'k','LineWidth',lw); % BUOY
chppf = plot(xmonth, CycleTable(:,7)./CycleTable(:,1),'-k','LineWidth',lw); chppf.Color(4) = 0.5; % Polarpathfinder
chfd = plot(xmonth, CycleTable(:,3)./CycleTable(:,1),'b','LineWidth',lw); % Free drift h
chfd2 = plot(xmonth, CycleTable(:,5)./CycleTable(:,1),'-.b','LineWidth',lw); % Free drift cst

grid on;
xlim([1 12]);
ylim([0 1.5]);
ylabel('RMSE/mean');
set(gca,'FontSize',25);

xticks(xmonth);
xticklabels({'Jan','Feb','Mar','Apr','May','Jun','Jul','Aug','Sep','Oct',...
             'Nov','Dec'});
xlabel('Month');
yticks(0:0.25:1.5);
yticklabels({'0','','0.5','','1','','1.5'});

legend([chppf chfd chfd2],{'PPF wind', 'FD thickness','FD constant'},'location','southeast','FontSize',15);

title('Coefficient of variation')   
text(0.02,0.98,'c)','Units', 'Normalized', 'VerticalAlignment', 'Top','Fontsize',25);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
n4 = nexttile([1 2]); % Time series of the coefficient of variation
hold on; box on;
%bhb = plot(xyear, TSTable(:,1),'k','LineWidth',lw); % BUOY
dhppf = plot(xyear, TSTable(:,7)./TSTable(:,1),'-k','LineWidth',lw); dhppf.Color(4) = 0.5; % Polarpathfinder 
dhfd = plot(xyear, TSTable(:,3)./TSTable(:,1),'b','LineWidth',lw); % Free drift h
dhfd2 = plot(xyear, TSTable(:,5)./TSTable(:,1),'-.b','LineWidth',lw); % Free drift cst

grid on;
xlim([1979 2019]);
ylim([0 1.5]);
ylabel('RMSE/mean');
set(gca,'FontSize',25);

xticks(1979:1:2018);
xticklabels({'','','1980','','','','','1985','','','','',...
                   '1990','','','','','1995','','','','',...
                   '2000','','','','','2005','','','','',...
                   '2010','','','','','2015','','','',''});
xlabel('Year');
yticks(0:0.25:1.5);
yticklabels({'0','','0.5','','1','','1.5'});

legend([dhppf dhfd dhfd2],{'PPF wind', 'FD thickness','FD constant'},'location','southeast','FontSize',15);

title('Coefficient of variation')   
text(0.02,0.98,'d)','Units', 'Normalized', 'VerticalAlignment', 'Top','Fontsize',25);


%% Figure 5 - Scatter plot & Histograms of speed & Delta Theta distribution for H-param; cst-param; PPF-wind

spacing = 25; % set to 5 for final figure
deltascale = 40;
xtext = 0.05;
ytext = 0.98;
formatSpec = '%.1f';

fig5 = figure(5);   %Make figure handle
set(fig5,'units','pixel','outerposition',[1281 1 1040 1080]) %[1281 1 960 1080]
tlo = tiledlayout(3,3,'TileSpacing','compact','Padding','compact');

%%%%%%%%%%%%%%%%%%% Scatter plot h-param
nexttile;
hold on; box on; grid on;
title('Thickness parameterization','Fontsize',15);

%First set of axes - scatter plot
scatter_kde(dufd(1:spacing:end),dvfd(1:spacing:end), '.');
ax1 = gca;
ax1.XLabel.String = '\Delta U_i [cm/s]';
ax1.YLabel.String = '\Delta V_i [cm/s]';
set(ax1, 'XLim', [-deltascale deltascale],...
         'YLim', [-deltascale deltascale], ...
         'FontSize', 12);
    
% Plotting two circles with radius equivalent to mean speed:
c1 = viscircles(ax1,[0 0],mean(speedbuoy),'Color','r','Linestyle','-','EnhanceVisibility',false); % Buoy speed circle

legend([c1],{'Buoy mean'});
text(xtext,ytext,'a)','Units', 'Normalized', 'VerticalAlignment', 'Top','Fontsize',25);
set(gca,'FontSize',15);

%%%%%%%%%%%%%%%%%%% Scatter plot constant
nexttile;
hold on; box on; grid on;
title('Constant \alpha','Fontsize',15);

%First set of axes - scatter plot
scatter_kde(dufd2(1:spacing:end),dvfd2(1:spacing:end), '.');
ax2 = gca;
ax2.XLabel.String = '\Delta U_i [cm/s]';
ax2.YLabel.String = '\Delta V_i [cm/s]';
set(ax2, 'XLim', [-deltascale deltascale],...
         'YLim', [-deltascale deltascale], ...
         'FontSize', 12);
    
% Plotting two circles with radius equivalent to mean speed:
c1 = viscircles(ax2,[0 0],mean(speedbuoy),'Color','r','Linestyle','-','EnhanceVisibility',false); % Buoy speed circle

%legend([c1],{'Buoy mean'});
text(xtext,ytext,'b)','Units', 'Normalized', 'VerticalAlignment', 'Top','Fontsize',25);
set(gca,'FontSize',15);


%%%%%%%%%%%%%%%%%%% Scatter plot PPF
nexttile;
hold on; box on; grid on;
title('Polar Pathfinder','Fontsize',15);

%First set of axes - scatter plot
scatter_kde(duppf(1:spacing:end),dvppf(1:spacing:end), '.');
ax3 = gca;
ax3.XLabel.String = '\Delta U_i [cm/s]';
ax3.YLabel.String = '\Delta V_i [cm/s]';
set(ax3, 'XLim', [-deltascale deltascale],...
         'YLim', [-deltascale deltascale], ...
         'FontSize', 12);
    
% Plotting two circles with radius equivalent to mean speed:
c1 = viscircles(ax3,[0 0],mean(speedbuoy),'Color','r','Linestyle','-','EnhanceVisibility',false); % Buoy speed circle

%legend([c1],{'Buoy mean'});
text(xtext,ytext,'c)','Units', 'Normalized', 'VerticalAlignment', 'Top','Fontsize',25);
set(gca,'FontSize',15);


%%%%%%%%%%%%%%%%%%%%%%%%% SPEED DISTRIBUTIONS

%%%%% H-PARAM distribution
nexttile; 
hold on; grid on; box on;
 h1b = histogram(speedbuoy,'BinEdges',[0:.5:25],'Normalization','PDF','FaceColor','#0000FF','FaceAlpha',0.6);
 h1fd = histogram(speedfd,'BinEdges',[0:.5:25],'Normalization','PDF','FaceColor','#FF0000','FaceAlpha',0.6);
xlim([0 25])
xlabel('|U_i| [cm/s]');
ylabel('PDF'); 
legend([h1b h1fd],{'Buoy', '\alpha (h_i)'},'location','northeast'); 
text(xtext,ytext,'d)','Units', 'Normalized', 'VerticalAlignment', 'Top','Fontsize',25);
set(gca,'Fontsize',15);

%%%%% H-constant distribution
nexttile; 
hold on; grid on; box on;
 h2b = histogram(speedbuoy,'BinEdges',[0:.5:25],'Normalization','PDF','FaceColor','#0000FF','FaceAlpha',0.6);
 h2fd = histogram(speedfd2,'BinEdges',[0:.5:25],'Normalization','PDF','FaceColor','#FF0000','FaceAlpha',0.6);
xlim([0 25])
xlabel('|U_i| [cm/s]');
ylabel('PDF'); 
legend([h2b h2fd],{'Buoy', '\alpha (constant)'},'location','northeast'); 
text(xtext,ytext,'e)','Units', 'Normalized', 'VerticalAlignment', 'Top','Fontsize',25);
set(gca,'Fontsize',15);

%%%%% PPF distribution
nexttile; 
hold on; grid on; box on;
 h3b = histogram(speedbuoy,'BinEdges',[0:.5:25],'Normalization','PDF','FaceColor','#0000FF','FaceAlpha',0.6);
 h3fd = histogram(speedppf,'BinEdges',[0:.5:25],'Normalization','PDF','FaceColor','#FF0000','FaceAlpha',0.6);
xlim([0 25])
xlabel('|U_i| [cm/s]');
ylabel('PDF'); 
legend([h3b h3fd],{'Buoy', 'Polar Pathfinder'},'location','northeast'); 
text(xtext,ytext,'f)','Units', 'Normalized', 'VerticalAlignment', 'Top','Fontsize',25);
set(gca,'Fontsize',15);


%%%%%%%%%%%%%%%%%%%%%%%%% Angle error distribution

%%%%% H-PARAM distribution
nexttile; 
hold on; grid on; box on;
 ha1 = histogram(dthetafd,'BinEdges',[-180:5:180],'Normalization','PDF','FaceColor','#0000FF','FaceAlpha',0.6);
 legdum1 = plot(0,-1,'.w');
xlim([-180 180]); ylim([0 0.02]);
xlabel('\Delta \theta [°]');
ylabel('PDF'); 
legend([legdum1],{['$\overline{\Delta \theta}=$' num2str(mean(dthetafd),formatSpec)]},'interpreter','Latex','Fontsize',18);
text(xtext,ytext,'g)','Units', 'Normalized', 'VerticalAlignment', 'Top','Fontsize',25);
set(gca,'Fontsize',15);

%%%%% Constant alpha distribution
nexttile; 
hold on; grid on; box on;
 ha2 = histogram(dthetafd2,'BinEdges',[-180:5:180],'Normalization','PDF','FaceColor','#0000FF','FaceAlpha',0.6);
 legdum2 = plot(0,-1,'.w');
xlim([-180 180]); ylim([0 0.02]);
xlabel('\Delta \theta [°]');
ylabel('PDF'); 
legend([legdum2],{['$\overline{\Delta \theta}=$' num2str(mean(dthetafd2),formatSpec)]},'interpreter','Latex','Fontsize',18);
text(xtext,ytext,'h)','Units', 'Normalized', 'VerticalAlignment', 'Top','Fontsize',25);
set(gca,'Fontsize',15);

%%%%% PPF distribution
nexttile; 
hold on; grid on; box on;
 ha3 = histogram(dthetappf,'BinEdges',[-180:5:180],'Normalization','PDF','FaceColor','#0000FF','FaceAlpha',0.6);
legdum3 = plot(0,-1,'.w');
xlim([-180 180]); ylim([0 0.02]);
xlabel('\Delta \theta [°]');
ylabel('PDF'); 
legend([legdum3],{['$\overline{\Delta \theta}=$' num2str(nanmean(dthetappf),formatSpec)]},'interpreter','Latex','Fontsize',18);
text(xtext,ytext,'i)','Units', 'Normalized', 'VerticalAlignment', 'Top','Fontsize',25);
set(gca,'Fontsize',15);

%%%%%%%%%%%% DENSITY CURVES ON SECONDARY AXES
%Second set of axes - U density curves
%%%%% CAREFUL ABOUT THE SET(AXIS XLIM YLIM) command - use the smae scale in
%%%%% between the different figures so that the scale is the same in
%%%%% between the three gaussian distributions in the three plots

%%%%%%%% FOR H PARAM
ax1_pos = ax1.Position; % position of first axes
ax12 = axes('Position',ax1_pos,...
    'XAxisLocation','top',...
    'YAxisLocation','right',...
    'Color','none',...
    'XColor','none',...
    'YColor','none');

[fu,xiu] = ksdensity(dufd,[-100:0.5:100]); 
line(xiu,fu,'Parent',ax12,'Color','k','Linewidth',2);

set(ax12, 'XLim', [-deltascale deltascale],...
         'YLim', [0 max(fu)*4]);
     
%Third set of axes - V density curves
ax13 = axes('Position',ax1_pos,...
    'XAxisLocation','top',...
    'YAxisLocation','right',...
    'Color','none',...
    'XColor','none',...
    'YColor','none');

[fv,xiv] = ksdensity(dvfd,[-100:0.5:100]); 
line(fv,xiv,'Parent',ax13,'Color','k','Linewidth',2);

set(ax13, 'YLim', [-deltascale deltascale],...
         'XLim', [0 max(fv)*4]);
     
    %%%%%%%% FOR Constant PARAM
    ax2_pos = ax2.Position; % position of first axes
    ax22 = axes('Position',ax2_pos,...
        'XAxisLocation','top',...
        'YAxisLocation','right',...
        'Color','none',...
        'XColor','none',...
        'YColor','none');

    [fu2,xiu2] = ksdensity(dufd2,[-100:0.5:100]); 
    line(xiu2,fu2,'Parent',ax22,'Color','k','Linewidth',2);

    set(ax22, 'XLim', [-deltascale deltascale],...
             'YLim', [0 max(fu)*4]);

    %Third set of axes - V density curves
    ax23 = axes('Position',ax2_pos,...
        'XAxisLocation','top',...
        'YAxisLocation','right',...
        'Color','none',...
        'XColor','none',...
        'YColor','none');

    [fv2,xiv2] = ksdensity(dvfd2,[-100:0.5:100]); 
    line(fv2,xiv2,'Parent',ax23,'Color','k','Linewidth',2);

    set(ax23, 'YLim', [-deltascale deltascale],...
             'XLim', [0 max(fv)*4]);
     
        %%%%%%%% FOR PPF
        ax3_pos = ax3.Position; % position of first axes
        ax32 = axes('Position',ax3_pos,...
            'XAxisLocation','top',...
            'YAxisLocation','right',...
            'Color','none',...
            'XColor','none',...
            'YColor','none');

        [fu3,xiu3] = ksdensity(duppf,[-100:0.5:100]); 
        line(xiu3,fu3,'Parent',ax32,'Color','k','Linewidth',2);

        set(ax32, 'XLim', [-deltascale deltascale],...
                 'YLim', [0 max(fu)*4]);

        %Third set of axes - V density curves
        ax33 = axes('Position',ax3_pos,...
            'XAxisLocation','top',...
            'YAxisLocation','right',...
            'Color','none',...
            'XColor','none',...
            'YColor','none');

        [fv3,xiv3] = ksdensity(dvppf,[-100:0.5:100]); 
        line(fv3,xiv3,'Parent',ax33,'Color','k','Linewidth',2);

        set(ax33, 'YLim', [-deltascale deltascale],...
                 'XLim', [0 max(fv)*4]);
             
%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%% Error on the angle as a function of the drift speed %%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
spacing = 3; %-> uses every x points in the computation of the scatter_kde (otherwise it takes forever to compute the kde density)

% Delta angles were: dthetafd, dtherappf, dthetafd2
             
fig52 = figure(52);   %Make figure handle
set(fig52,'units','normalized','outerposition',[0 0 .5 1])
hold on; box on; grid on;

 scatter_kde(speedfd(1:spacing:end),dthetafd(1:spacing:end), '.');
 
xlim([0 25]); xlabel('|U_i| [cm/s]');
ylim([-180 180]); ylabel('\Delta \theta [°]');
title('Angle error vs drift speed');

set(gca,'Fontsize',20);
c=colorbar;
xlabel(c,'kernel density estimate','Fontsize',15);
      
%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%% MAP OF THE RELATIVE ERROR FOR THE H_PARAM %%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% 1) Need to calculate  LOCAL RMSE AND MEAN BUOY
TablePath = '/Users/charlesbrunette/Dropbox/FREEDRIFT2020/Analysis/Freedrift/Data/';
load([TablePath 'XDATA_YDATA_cells.mat']);
UW = UO(:,end);
VW = VO(:,end);

% Requires the output of ITERATIVE_SOLVER_alpha_uo_CB.m
ALPHAFD = ALPHA(1,end);
BETA = ALPHA(2,end);
THETA = ALPHA(4,end);

rmsemap = NaN(361,361);
meanmap = NaN(361,361);
tic
    for i=1:length(Ind_isdata)
        
        j= Ind_isdata(i);
        UA = XDATA_cell{j,1};
        VA = XDATA_cell{j,2};
        H =  XDATA_cell{j,4};
        UB = YDATA_cell{j,1}; 
        VB = YDATA_cell{j,2}; 
        
        %Thickness dependence parameterization
        UFD = ALPHAFD*(1 - BETA.*H).*(UA.*cos(THETA) - VA.*sin(THETA)) + UW(i);
        VFD = ALPHAFD*(1 - BETA.*H).*(VA.*cos(THETA) + UA.*sin(THETA)) + VW(i);
        
        FDspeed = sqrt(UFD.^2+VFD.^2);
        Bspeed =  sqrt(UB.^2+VB.^2);
        
        rmsemap(j) = sqrt(mean((FDspeed-Bspeed).^2));
        meanmap(j) = mean(Bspeed);
        disp(['Calculating...' int2str(i) '/' int2str(length(Ind_isdata))]);
    end
toc

%% 2) Doing it for a different parameterization:

UW = UO(:,end);
VW = VO(:,end);

% Requires the output of ITERATIVE_SOLVER_alpha_uo_CB.m
ALPHAFD = ALPHA(1,end);
THETA = ALPHA(4,end);

rmsemap2 = NaN(361,361);
meanmap2 = NaN(361,361);
tic
    for i=1:length(Ind_isdata)
        
        j= Ind_isdata(i);
        UA = XDATA_cell{j,1};
        VA = XDATA_cell{j,2};
        UB = YDATA_cell{j,1}; 
        VB = YDATA_cell{j,2}; 
        
        %Thickness dependence parameterization
        UFD = ALPHAFD.*(UA.*cos(THETA) - VA.*sin(THETA)) + UW(i);
        VFD = ALPHAFD.*(VA.*cos(THETA) + UA.*sin(THETA)) + VW(i);
        
        FDspeed = sqrt(UFD.^2+VFD.^2);
        Bspeed =  sqrt(UB.^2+VB.^2);
        
        rmsemap2(j) = sqrt(mean((FDspeed-Bspeed).^2));
        meanmap2(j) = mean(Bspeed);
        disp(['Calculating...' int2str(i) '/' int2str(length(Ind_isdata))]);
    end
toc

% and for polarpathfinder wind component:

rmsemap3 = NaN(361,361);
meanmap3 = NaN(361,361);
tic
    for i=1:length(Ind_isdata)
        
        j= Ind_isdata(i);
        UA = PPFw_cell{j,1};
        VA = PPFw_cell{j,2};
        UB = YDATA_cell{j,1}; 
        VB = YDATA_cell{j,2}; 
        
        %Thickness dependence parameterization
        UFD = UA;
        VFD = VA;
        
        FDspeed = sqrt(UFD.^2+VFD.^2);
        Bspeed =  sqrt(UB.^2+VB.^2);
        
        rmsemap3(j) = sqrt(mean((FDspeed-Bspeed).^2));
        meanmap3(j) = mean(Bspeed);
        disp(['Calculating...' int2str(i) '/' int2str(length(Ind_isdata))]);
    end
toc

%% Figure for only one error map
    
fig6 = figure(6);

set(fig6,'units','normalized','outerposition',[0 0 .5 1])
tlo = tiledlayout(4,3,'TileSpacing','compact','Padding','compact');

%%%% MAP OF CV
nexttile([3 3]);

levels = 0:0.1:1;

ncpolarm('lat',66) % Create Arctic Mask
box on; hold on;
%pcolorm(EASElat,EASElon,Alpha361');
contourfm(EASElat,EASElon,rmsemap'./meanmap',levels)%,'LineStyle','none');
cb = contourcbar;
cb.Location = 'westoutside';
%xlabel(cb,'$C_v=\frac{\sigma_\epsilon}{\mu}$','Interpreter','latex');
xlabel(cb,'RMSE/mean');
caxis([0 1]);                                              
title({'Relative error'},'FontSize', 25);
textm(60,-135,'a)','Fontsize',25);
set(gca,'Fontsize',20);


nexttile([1 3]); %%%%%%%% Distribution of current speed
hold on; grid on; box on;
h1 = histogram(rmsemap'./meanmap','BinEdges',[0:0.05:1],'Normalization','PDF');
xlim([0 1]); ylim([0 0.35])
xlabel('RMSE/mean');
ylabel('PDF');
text(0.02,0.98,'b)','Units', 'Normalized', 'VerticalAlignment', 'Top','Fontsize',25);
set(gca,'Fontsize',20);
%set(gca,'YAxisLocation','right')

%% Figure for multiple error maps
    
fig6 = figure(6);

set(fig6,'units','normalized','outerposition',[0 0 .5 1])
tlo = tiledlayout(6,6,'TileSpacing','compact','Padding','compact');

%%%% MAP OF CV1
nexttile([3 3]);

levels = 0:0.1:1;

ncpolarm('lat',66) % Create Arctic Mask
box on; hold on;
%pcolorm(EASElat,EASElon,Alpha361');
contourfm(EASElat,EASElon,rmsemap3'./meanmap3',levels)%,'LineStyle','none');
%cb = contourcbar;
%cb.Location = 'westoutside';
%xlabel(cb,'$C_v=\frac{\sigma_\epsilon}{\mu}$','Interpreter','latex');
%xlabel(cb,'RMSE/mean');
caxis([min(levels) max(levels)]);                                              
title({'Relative error: PPF'},'FontSize', 15);
textm(60,-135,'a)','Fontsize',25);
set(gca,'Fontsize',15);

%%%% MAP OF CV1
nexttile([3 3]);

levels = 0:0.1:1;

ncpolarm('lat',66) % Create Arctic Mask
box on; hold on;
%pcolorm(EASElat,EASElon,Alpha361');
contourfm(EASElat,EASElon,rmsemap2'./meanmap2',levels)%,'LineStyle','none');
%cb = contourcbar;
%cb.Location = 'westoutside';
%xlabel(cb,'$C_v=\frac{\sigma_\epsilon}{\mu}$','Interpreter','latex');
%xlabel(cb,'RMSE/mean');
caxis([min(levels) max(levels)]);                                              
title({'Relative error: \alpha_{ow}'},'FontSize', 15);
textm(60,-135,'b)','Fontsize',25);
set(gca,'Fontsize',15);

%%%% MAP OF CV2
nexttile([3 3]);

ncpolarm('lat',66) % Create Arctic Mask
box on; hold on;
%pcolorm(EASElat,EASElon,Alpha361');
contourfm(EASElat,EASElon,rmsemap'./meanmap',levels)%,'LineStyle','none');
cb = contourcbar;
cb.Location = 'southoutside';
%xlabel(cb,'$C_v=\frac{\sigma_\epsilon}{\mu}$','Interpreter','latex');
xlabel(cb,'RMSE/mean');
caxis([min(levels) max(levels)]);                                             
title({'Relative error: \alpha_h'},'FontSize', 15);
textm(60,-135,'c)','Fontsize',25);
set(gca,'Fontsize',15);

nexttile([1 3]); %%%%%%%% Distribution of current speed
hold on; grid on; box on;
h1 = histogram(rmsemap3'./meanmap3','BinEdges',[0:0.05:1],'Normalization','PDF');
xlim([0 1]); ylim([0 0.35])
%xlabel('RMSE/mean');
ylabel('PDF');
text(0.02,0.98,'d) PPF','Units', 'Normalized', 'VerticalAlignment', 'Top','Fontsize',25);
set(gca,'Fontsize',15);
%set(gca,'YAxisLocation','right')

nexttile([1 3]); %%%%%%%% Distribution of current speed
hold on; grid on; box on;
h1 = histogram(rmsemap2'./meanmap2','BinEdges',[0:0.05:1],'Normalization','PDF');
xlim([0 1]); ylim([0 0.35])
%xlabel('RMSE/mean');
ylabel('PDF');
text(0.02,0.98,'e) \alpha_{ow}','Units', 'Normalized', 'VerticalAlignment', 'Top','Fontsize',25);
set(gca,'Fontsize',15);
%set(gca,'YAxisLocation','right')

nexttile([1 3]); %%%%%%%% Distribution of current speed
hold on; grid on; box on;
h1 = histogram(rmsemap'./meanmap','BinEdges',[0:0.05:1],'Normalization','PDF');
xlim([0 1]); ylim([0 0.35])
xlabel('RMSE/mean');
ylabel('PDF');
text(0.02,0.98,'f) \alpha_h','Units', 'Normalized', 'VerticalAlignment', 'Top','Fontsize',25);
set(gca,'Fontsize',15);
%set(gca,'YAxisLocation','right')



%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%% MAP OF THE NUMBER OF POINTS GOING IN THE EVALUATIONof UW  %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

load([TablePath 'XDATA_YDATA_cells.mat']);
Npoints = reshape(Datapoint_numbers,[361 361]); Npoints(Npoints<1)=NaN; % Matrix containing number of points going in the evaluation of each Uw_i,j

figappx = figure(99);
set(figappx,'units','normalized','outerposition',[0 0 .5 1])

tlo = tiledlayout(4,2,'TileSpacing','compact','Padding','compact'); % Create Tiles structure

nexttile([3 2]); % -> Npoints map

levels = 0:100:800;

ncpolarm('lat',66) % Create Arctic Mask
box on; hold on;
%pcolorm(EASElat,EASElon,Alpha361');
contourfm(EASElat,EASElon,Npoints',levels)%,'LineStyle','none');
cb = contourcbar;
colormap parula;
cb.Location = 'westoutside';
xlabel(cb,'#points');
caxis([0 800]);                                              
title({'Density of data points'},'FontSize', 20);
textm(60,-135,'a)','Fontsize',25);

set(gca,'Fontsize',20);


nexttile([1 2]); %%%%%%%% Distribution of current speed
hold on; grid on; box on;
h1 = histogram(Npoints,'BinEdges',[0:10:1000]);
xlim([0 1000])
xlabel('Npoints');
ylabel('Count');
text(0.02,0.98,'b)','Units', 'Normalized', 'VerticalAlignment', 'Top','Fontsize',25);
set(gca,'Fontsize',20);
%set(gca,'YAxisLocation','right')





%% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% FIGURE COMPARING Uo estimates TO DIFFERENT REANALYSIS - Second version
% 2003-2014
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%) 0 %2) Estimates over 2003-2014 period
                                load([TablePath 'XDATA_YDATA_cells_20032014.mat']); Ind_isdata = find(Datapoint_numbers>0);
                                load('/Users/charlesbrunette/Dropbox/FREEDRIFT2020/Analysis/Freedrift/Data/2003_2014_ThicknessRun.mat')
                                uocnc = NaN(361*361,1); vocnc=NaN(361*361,1);
                                uocnc(Ind_isdata) = UO(:,end); %uo(ArcticOceanMask4IceEdge ~=1)=NaN;
                                vocnc(Ind_isdata) = VO(:,end); %vo(ArcticOceanMask4IceEdge ~=1)=NaN;
                                uocnc=reshape(uocnc,[361,361]);
                                vocnc=reshape(vocnc,[361,361]);
                                %Rotate uo vo so they point along LAT/LON lines
                                Uo_ease = uocnc + 1i*vocnc;
                                Uo_ll_cms = Uo_ease.*exp(1i*deg2rad(EASElon'));
                                uomag = sqrt(UO(:,end).^2+VO(:,end).^2);
                                gridmag = sqrt(uocnc.^2+vocnc.^2);

                                % 1) ECCO
                                load('/Users/charlesbrunette/Dropbox/FREEDRIFT2020/Analysis/Freedrift/Data/ECCO_EASE_2003_2014.mat');
   
                                uECCO = uocnEASE';
                                vECCO = vocnEASE';
   
                                %Rotate uo vo so they point along LAT/LON lines
                                UoECCO = uECCO + 1i*vECCO;
                                Uo_ll_cmsECCO = UoECCO.*exp(1i*deg2rad(EASElon'));
                                gridmagECCO = sqrt(uECCO.^2+vECCO.^2);

                                
                                % 2) GLORYS
                                load('/Users/charlesbrunette/Dropbox/FREEDRIFT2020/Analysis/Freedrift/Data/GLORYS_EASE_2003_2014.mat');
   
                                uGLORYS = uocnEASE';
                                vGLORYS = vocnEASE';
   
                                %Rotate uo vo so they point along LAT/LON lines
                                UoGLORYS = uGLORYS + 1i*vGLORYS;
                                Uo_ll_cmsGLORYS = UoGLORYS.*exp(1i*deg2rad(EASElon'));
                                gridmagGLORYS = sqrt(uGLORYS.^2+vGLORYS.^2);
                             

                                 % 3) CPOM/DOT
                                load('/Users/charlesbrunette/Dropbox/FREEDRIFT2020/Analysis/Freedrift/Data/CPOMDOT_EASE_2003_2014.mat');
   
                                ucpom = uocnEASE';
                                vcpom = vocnEASE';
   
                                %Rotate uo vo so they point along LAT/LON lines
                                Uocpom = ucpom + 1i*vcpom;
                                Uo_ll_cmscpom = Uocpom.*exp(1i*deg2rad(EASElon'));
                                gridmagcpom = sqrt(ucpom.^2+vcpom.^2);

                                % 3) PIOMAS
                                load('/Users/charlesbrunette/Dropbox/FREEDRIFT2020/Analysis/Freedrift/Data/PIOMAS_UO_EASE_2003_2014.mat');
   
                                ucpio = uocnEASE';
                                vcpio = vocnEASE';
   
                                %Rotate uo vo so they point along LAT/LON lines
                                Uocpio = ucpio + 1i*vcpio;
                                Uo_ll_cmscpio = Uocpio.*exp(1i*deg2rad(EASElon'));
                                gridmagcpio = sqrt(ucpio.^2+vcpio.^2);
                                     %iz = find(gridmagcpio==0); % Getting rid of zeroes on landmasses
                                     %gridmagcpio(iz)=NaN;
                                              
%%

                        %%%%%%% FIGURE Ocean current comparison %%%%%%%%%%%%%%

                                cmap = redblue; % Creating a red (+) to blue (-) colormap; 

                                fig23 = figure(23);   %Make figure handle
                                tlo = tiledlayout(3,4,'TileSpacing','compact','Padding','compact');
                                
                                
                               nexttile; %%%%%%%% BUOY DERIVED over 2003-2014

                                %- Colormap with unit scale vectors:
                                    ncpolarm('lat',66); hold on;
                                    f23hcol = pcolorm(EASElat',EASElon',gridmag);
                                    cb = colorbar('Westoutside');
                                    colormap(gcf,'parula')
                                    % Vector sequence:
                                    quivermc(EASElat',EASElon',real(Uo_ll_cms)./gridmag,-imag(Uo_ll_cms)./gridmag,'density',30,'reference', 1,'linewidth',1,'Fontsize',15);
                                    xlabel(cb,'U_w [cm/s]','Fontsize',15)
                                    cb.Ticks = [0:1:5]; % vector version was going to 10 cm/s
                                    caxis([0 5]);
                                    cb.TickLabels = {'0','1','2','3','4','5'};                                             
                                    title({'Clim (2003-2014)'},'FontSize', 15);
                                    textm(60,-135,'a)','Fontsize',25);
                                    box on;
                                    set(gca,'Fontsize',15);
                                
                                
                                
                                nexttile(5); %%%%%%%% ECCO

                                %- Colormap with unit scale vectors:
                                    ncpolarm('lat',66); hold on;
                                    f23acol = pcolorm(EASElat',EASElon',gridmagECCO);
                                    cb = colorbar('Westoutside');
                                    colormap(gcf,'parula')
                                    % Vector sequence:
                                    quivermc(EASElat',EASElon',real(Uo_ll_cmsECCO)./gridmagECCO,-imag(Uo_ll_cmsECCO)./gridmagECCO,'density',30,'reference', 1,'linewidth',1,'Fontsize',15);
                                    xlabel(cb,'U_w [cm/s]','Fontsize',15)
                                    cb.Ticks = [0:1:5]; % vector version was going to 10 cm/s
                                    caxis([0 5]);
                                    cb.TickLabels = {'0','1','2','3','4','5'};

                                    %gridm on;                                              
                                    title({'ECCO'},'FontSize', 15);
                                    textm(60,-135,'b)','Fontsize',25);
                                    box on;
                                    set(gca,'Fontsize',15);
                                    

                                nexttile(6); %%%%%%%% GLORYS

                                %- Colormap with unit scale vectors:
                                    ncpolarm('lat',66); hold on;
                                    f23dcol = pcolorm(EASElat',EASElon',gridmagGLORYS);
                                    %cb = colorbar('Southoutside');
                                    colormap(gcf,'parula')
                                    % Vector sequence:
                                    quivermc(EASElat',EASElon',real(Uo_ll_cmsGLORYS)./gridmagGLORYS,-imag(Uo_ll_cmsGLORYS)./gridmagGLORYS,'density',30,'reference', 1,'linewidth',1,'Fontsize',15);
                                    %xlabel(cb,'U_w [cm/s]','Fontsize',15)
                                    %cb.Ticks = [0:1:5]; % vector version was going to 10 cm/s
                                    caxis([0 5]);
                                    %cb.TickLabels = {'0','1','2','3','4','5'};

                                    %gridm on;                                              
                                    title({'GLORYS'},'FontSize', 15);
                                    textm(60,-135,'c)','Fontsize',25);
                                    box on;
                                    set(gca,'Fontsize',15);
                                    
                                nexttile(7); %%%%%%%% PIOMAS

                                %- Colormap with unit scale vectors:
                                    ncpolarm('lat',66); hold on;
                                    f23jcol = pcolorm(EASElat',EASElon',gridmagcpio);
                                    %cb = colorbar('Southoutside');
                                    colormap(gcf,'parula')
                                    % Vector sequence:
                                    quivermc(EASElat',EASElon',real(Uo_ll_cmscpio)./gridmagcpio,-imag(Uo_ll_cmscpio)./gridmagcpio,'density',30,'reference', 1,'linewidth',1,'Fontsize',15);
                                    %xlabel(cb,'U_w [cm/s]','Fontsize',15)
                                    %cb.Ticks = [0:1:5]; % vector version was going to 10 cm/s
                                    caxis([0 5]);
                                    %cb.TickLabels = {'0','1','2','3','4','5'};

                                    %gridm on;                                              
                                    title({'PIOMAS'},'FontSize', 15);
                                    textm(60,-135,'d)','Fontsize',25);
                                    box on;
                                    set(gca,'Fontsize',15);
                                    
                                    
                                 nexttile(8); %%%%%%%% CPOM/DOT

                                %- Colormap with unit scale vectors:
                                    ncpolarm('lat',66); hold on;
                                    f23gcol = pcolorm(EASElat',EASElon',gridmagcpom);
                                    %cb = colorbar('Southoutside');
                                    colormap(gcf,'parula')
                                    % Vector sequence:
                                    quivermc(EASElat',EASElon',real(Uo_ll_cmscpom)./gridmagcpom,-imag(Uo_ll_cmscpom)./gridmagcpom,'density',30,'reference', 1,'linewidth',1,'Fontsize',15);
                                    %xlabel(cb,'U_w [cm/s]','Fontsize',15)
                                    %cb.Ticks = [0:1:5]; % vector version was going to 10 cm/s
                                    caxis([0 5]);
                                    %cb.TickLabels = {'0','1','2','3','4','5'};

                                    %gridm on;                                              
                                    title({'CPOM/DOT'},'FontSize', 15);
                                    textm(60,-135,'e)','Fontsize',25);
                                    box on;
                                    set(gca,'Fontsize',15); 
                                    
                                    

                                nf=nexttile(9); %%%%%%%% DELTA Estimates - ECCO

                                ncpolarm('lat',66);
                                hold on;
                                gb3 = pcolorm(EASElat',EASElon',gridmag-gridmagECCO);
                                cb3 = colorbar('Westoutside');
                                set(nf,'colormap', cmap);
                                cb3.Ticks = [-5:1:5];
                                caxis([-5 5]);
                                cb3.TickLabels = {'-5','-4','-3','-2','-1','0','1','2','3','4','5'};
                                xlabel(cb3,'\Delta U_w [cm/s]','Fontsize',15)
                                title({'Clim - ECCO'});
                                textm(60,-135,'f)','Fontsize',25);
                                box on;
                                set(gca,'Fontsize',15);
  
                                    
                              ng=nexttile(10); %%%%%%%% DELTA Estimates - GLORYS

                                ncpolarm('lat',66);
                                hold on;
                                gb3 = pcolorm(EASElat',EASElon',gridmag-gridmagGLORYS);
                                %cb3 = colorbar('Southoutside');
                                set(ng,'colormap', cmap);
                                %cb3.Ticks = [-5:1:5];
                                caxis([-5 5]);
                                %cb3.TickLabels = {'-5','-4','-3','-2','-1','0','1','2','3','4','5'};
                                %xlabel(cb3,'\Delta U_w [cm/s]','Fontsize',15)
                                title({'Clim - GLORYS'});
                                textm(60,-135,'g)','Fontsize',25);
                                box on;
                                set(gca,'Fontsize',15);
                                
                                ni=nexttile(11);%%%%%%%% DELTA Estimates - PIOMAS

                                ncpolarm('lat',66);
                                hold on;
                                gb3 = pcolorm(EASElat',EASElon',gridmag-gridmagcpio);
                                %cb3 = colorbar('Southoutside');
                                set(ni,'colormap', cmap);
                                %cb3.Ticks = [-5:1:5];
                                caxis([-5 5]);
                                %cb3.TickLabels = {'-5','-4','-3','-2','-1','0','1','2','3','4','5'};
                                %xlabel(cb3,'\Delta U_w [cm/s]','Fontsize',15)
                                title({'Clim - PIOMAS'});
                                textm(60,-135,'h)','Fontsize',25);
                                box on;
                                set(gca,'Fontsize',15);
                                

                            nh=nexttile(12); %%%%%%%% DELTA Estimates - CPOM/DOT

                                ncpolarm('lat',66);
                                hold on;
                                gb3 = pcolorm(EASElat',EASElon',gridmag-gridmagcpom);
                                %cb3 = colorbar('Southoutside');
                                set(nh,'colormap', cmap);
                                %cb3.Ticks = [-5:1:5];
                                caxis([-5 5]);
                                %cb3.TickLabels = {'-5','-4','-3','-2','-1','0','1','2','3','4','5'};
                                %xlabel(cb3,'\Delta U_w [cm/s]','Fontsize',15)
                                title({'Clim - CPOM/DOT'});
                                textm(60,-135,'i)','Fontsize',25);
                                box on;
                                set(gca,'Fontsize',15);
                                
                                     
%% Seasonality of the ERA5 wind:

WindSiz = NaN(12,2019-1979+1);

for iy= 1979:2019
    for im=1:12
        id = find(DATA.year==iy & DATA.month==im);
        WindSiz(im,iy-1978)=nanmean(sqrt(DATA.u_PPFv4Wind(id).^2+DATA.v_PPFv4Wind(id).^2));
    end
end

WindSiz = NaN(12,1);

    for im=1:12
        id = find(DATA.month==im);
        WindSiz(im)=nanmean(sqrt(DATA.u_PPFv4Wind(id).^2+DATA.v_PPFv4Wind(id).^2));
    end
    
%% First figure of the paper update: Scatter plot PPF vs buoy & scatter plot Constant Alpha (no current) vs buoy:

Sb = sqrt(Y_input_mat(:,1).^2+Y_input_mat(:,2).^2);
Sf = sqrt(Ui_fit.^2+Vi_fit.^2);
Sppf = sqrt(DATA.u_PPFv4Wind(Ind_common).^2+DATA.v_PPFv4Wind(Ind_common).^2);

spacing = 10;
tlim=20;

fig01 = figure(1);   %Make figure handle
tlo = tiledlayout(1,3,'TileSpacing','compact','Padding','compact');

nexttile; % Scatter plot ppf vs iabp

hold on; box on; grid on;

scatter_kde(Sb(1:spacing:end),Sppf(1:spacing:end), '.');
 
xlim([0 tlim]); xlabel('|U_i|: IABP [cm/s]');
ylim([0 tlim]); ylabel('|U_i|: \alpha_{ppf} [cm/s]');
title('Polar Pathfinder');

set(gca,'Fontsize',15);
text(0.02,0.98,'a)','Units', 'Normalized', 'VerticalAlignment', 'Top','Fontsize',25);

nexttile; % Scatter plot constant alpha vs iabp

hold on; box on; grid on;

scatter_kde(Sb(1:spacing:end),Sf(1:spacing:end), '.');
 
xlim([0 tlim]); xlabel('|U_i|: IABP [cm/s]');
ylim([0 tlim]); ylabel('|U_i|: \alpha_o [cm/s]');
title('Constant \alpha_o');

set(gca,'Fontsize',15);
text(0.02,0.98,'b)','Units', 'Normalized', 'VerticalAlignment', 'Top','Fontsize',25);

nexttile; % error reduction histogram

hold on; box on; grid on;

hf = histogram(Sf-Sb,'Normalization','PDF','FaceColor','#0000FF','FaceAlpha',0.5); 
[maxcount, whichbin] = max(hf.Values); 
xf = xline((hf.BinEdges(whichbin)+hf.BinEdges(whichbin+1))/2,'-b','Linewidth',2);
xf2 = xline(mean(Sf-Sb),'--b','Linewidth',2);

hppf = histogram(Sppf-Sb,'Normalization','PDF','FaceColor','#FF0000','FaceAlpha',0.5); 
[maxcount, whichbin] = max(hppf.Values); 
xppf = xline((hppf.BinEdges(whichbin)+hppf.BinEdges(whichbin+1))/2,'-r','Linewidth',2);
xppf2 = xline(nanmean(Sppf-Sb),'--r','Linewidth',2);
 
xlim([-20 20]); xlabel('\Delta |U_i| [cm/s]');
ylim([0 0.12]); ylabel('PDF');
set(gca,'Fontsize',15);
title('Error distribution','Fontsize',15);
legend([xf xf2 xppf xppf2],{'\alpha_o: modal error', '\alpha_o: mean error', 'PPF: modal error', 'PPF: mean error'},'location','northeast','Fontsize',15); 

text(0.02,0.98,'c)','Units', 'Normalized', 'VerticalAlignment', 'Top','Fontsize',25);


