%% Plot Processed Data
% Trevor Harrison
% 

define_paths;
clear dirSimData model_flag % Not relevent to this code. Removed for clarity. 

testID = 9;

%% Load Data
% Using full-resolution, non-time-shifted data sets.
% For time-averaged values, change the data loaded below. 

model_data = 1; %[ 1 = original run; 2 = refined run]

if ~exist('fs','var')
    if model_data == 1
    load('sim_estimates_model_v1_0min_sensor_avg.mat');
    load('assembled_comparison_statistics_v1.mat');
    elseif model_data == 2
    load('sim_estimates_model_v2_0min_sensor_avg.mat');
    load('assembled_comparison_statistics_v2.mat');
    end
    clear meta;
end

if ~exist('meta','var')   
    libFolder = fullfile(dirObservationalData, 'lib');
    load(fullfile(libFolder, 'AP2_metadata.mat'));  
    load(fullfile(libFolder, 'AgatePass_SiteBoundary.mat'));
    load(fullfile(libFolder, 'AgatePassBathy.mat'));
    tTests = meta.tTests;
    nTest = length(tTests);
    nLocSK = size(meta.locSK,1);
end

dirSensor = fullfile(dirObservationalData, 'sensors');

% Load station keeping data
if ~exist('SK','var')
    SK = importStationKeepingData(fullfile(dirSensor,'ADCP_StationKeeping','compiled'), meta);
end

% Load SWIFT data
if ~exist('SWIFT','var')
    SWIFT = importSWIFTdata(fullfile(dirSensor,'ADCP_SWIFT','compiled'), meta);
end

% Load FLOAT SWARM data
if ~exist('FS','var')
    FILES = dir(fullfile(dirSensor,'FloatSwarm','compiled','Float*'));
    for i = 1:length(FILES)
    FS(i) = load(fullfile(FILES(i).folder, FILES(i).name)); 
    end
    clear FILES i;
end

% Load CTD data
if ~exist('CTD','var')
    FILES = dir(fullfile(dirSensor,'CTD','compiled','CC*'));
    for i = 1:length(FILES)
        CTD(i) = load(fullfile(FILES(i).folder, FILES(i).name));
    end
    CTD = CTD(~[CTD.IsInvalid]);
    CTD = CTD([3:7,9:12]); % Grab relevant profiles (skips doubles, and 1st)   
    castTimes = datetime([CTD.CastTimeUtc]', 'TimeZone','UTC');
    for i = 1:length(CTD)
        meanSOS(i) = mean(CTD(i).Sound_velocity);
    end
    %sos =  interp1(castTimes, meanSOS, meta.tTests(testID)+minutes(10));
end  

loc_ADCP_NOAA =  [47.71102, -122.56715];
% https://tidesandcurrents.noaa.gov/cdata/StationInfo?id=PUG1501#

%% Convert Meta Data into Local Coordinates
latRef = meta.locSK(3,1);
lonRef = meta.locSK(3,2);

% Define station keeping locations (local coordinates)
[locSK(:,2), locSK(:,1)] = geodetic2ned(meta.locSK(:,1),meta.locSK(:,2),0,latRef,lonRef,0,referenceEllipsoid('wgs84'));

% Need coastlines for plotting
[yCLwest, xCLwest] =  geodetic2ned(CLwest(:,1), CLwest(:,2), 0,...
    latRef, lonRef, 0, referenceEllipsoid('wgs84'));
[yCLeast, xCLeast] =  geodetic2ned(CLeast(:,1), CLeast(:,2), 0,...
    latRef, lonRef, 0, referenceEllipsoid('wgs84'));

% NOAA ADCP to local
[locNOAA(:,2), locNOAA(:,1)] = geodetic2ned(loc_ADCP_NOAA(:,1),loc_ADCP_NOAA(:,2),0,latRef,lonRef,0,referenceEllipsoid('wgs84'));

% Bathymetry (natively in UTM, convert to lat lon, then to local 
xB = bathyData.xB;
yB = bathyData.yB;
zB = -bathyData.zB; % Our bathymetry 
[latBs, lonBs] = utm2ll(xB, yB, 10);
[yBs, xBs] = geodetic2ned(latBs, lonBs, 0, latRef, lonRef, 0, referenceEllipsoid('wgs84'));
% Crop section we care about
xCrop = [-1000, 1500];
yCrop = [-1000, 2500];
bInside = xBs>xCrop(1) & xBs<xCrop(2) & yBs>yCrop(1) & yBs<yCrop(2);
xBs = xBs(bInside);
yBs = yBs(bInside);
zB = zB(bInside);
G_B = scatteredInterpolant(xBs, yBs, zB,'linear','none');
kBounds = boundary(xBs, yBs,1);

% Compiled Observational data is all in [m] reference to SK3. Do not
% require any adjustment here. 

%% Figure 17/18: PLOTTING DIFFERENCES INTO DISTRIBUTIONS
% Plots horizontal velocity magnitude differences clouds and distribtuions. 
bit_print = 1;

iSurvey = testID;
% Get subset of survey data
ux_f = FS(iSurvey).UX;
uy_f = FS(iSurvey).UY;
ux_fq = FS(iSurvey).UXQ;
uy_fq = FS(iSurvey).UYQ;

bSub = SWIFT.testID == iSurvey;
ux_s = SWIFT.UX(bSub);
uy_s = SWIFT.UY(bSub);
ux_sq = SWIFT.UXQ(bSub);
uy_sq = SWIFT.UYQ(bSub);

bSub = SK.testID == iSurvey;
ux_sk = SK.UX(bSub);
uy_sk = SK.UY(bSub);
ux_skq = SK.UXQ(bSub);
uy_skq = SK.UYQ(bSub);

% Assemble for bounds on the plots
allData = [ux_sk; ux_s; ux_f(:); uy_sk; uy_s; uy_f(:);...
           ux_skq; ux_sq; ux_fq(:); uy_skq; uy_sq; uy_fq(:)];
          
xL = [floor(min(allData)*10)/10, ceil(max(allData)*10)/10];

% Set correct bounds
if iSurvey <= 3
    xL = [-0.25, 2.25];
    x3 = [-0.5, 0.5];
elseif iSurvey == 4
    xL = [-0.75 0.75];
    x3 = [-0.5, 0.5];
elseif iSurvey == 5
    xL = [-2, 0.5];
    x3 = [-0.5, 0.5];
elseif iSurvey == 6
    xL = [-2.5, 0.5];
elseif iSurvey == 7
    xL = [-2.5, 1];
elseif iSurvey == 8
    xL = [-2.5, 1];
end
    
    
f2 = figure(2); clf(f2); f2.Color = 'w';
ax = buildAxes(f2, 3, 4);
f2.Units = 'inches'; 
f2.Position = [1.0000    7    6.5000    5.3417];


% Rows: SK, SWIFTS, uFloat
% Column 1: Ux v Uy
mSize = 3;
set(ax(:,1), 'NextPlot','add');
plot(ax(1,1), ux_f(:), uy_f(:), '.','MarkerSize',mSize);
plot(ax(1,1), ux_fq(:), uy_fq(:), '.','MarkerSize',mSize);

plot(ax(2,1), ux_s, uy_s, '.','MarkerSize',mSize);
plot(ax(2,1), ux_sq, uy_sq, '.','MarkerSize',mSize);

plot(ax(3,1), ux_sk, uy_sk, '.','MarkerSize',mSize);
plot(ax(3,1), ux_skq, uy_skq, '.','MarkerSize',mSize);

set(ax(:,1), 'XLim', xL, 'YLim', xL);

% Column 2: Uh_m v Uh_obs
set(ax(:,2),'XLim',[0, max(abs(xL))], 'YLim',[0, max(abs(xL))]); 
for i = 1:3
    plot(ax(i,2), [0, 3], [0,3], '-', 'Color', 0.2*[1,1,1]);
end

set(ax, 'NextPlot','add','ColorOrderIndex',1,'DataAspectRatio',[1,1,1]);
mSize = 3;
plot(ax(1,2), sqrt(ux_f(:).^2+uy_f(:).^2), sqrt(ux_fq(:).^2+uy_fq(:).^2), '.',...
    'MarkerSize', mSize);
plot(ax(2,2), sqrt(ux_s.^2+uy_s.^2), sqrt(ux_sq.^2+uy_sq.^2), '.',...
    'MarkerSize', mSize);
plot(ax(3,2), sqrt(ux_sk.^2+uy_sk.^2), sqrt(ux_skq.^2+uy_skq.^2), '.',...
    'MarkerSize', mSize);
set(ax(:,2),'XLim',[0, max(abs(xL))], 'YLim',[0, max(abs(xL))]); 

%
% Column 3: Uh_m - Uh_obs: Time shifting
for k = 1:3
    bRow = dev_type == k-1 & survey_id == iSurvey & dev_avg == 0 & model_avg == 0; 
    cBox = 0.5*[1,1,1];
    ax_temp = ax(k,3); 
    ax_temp.NextPlot = 'add';
        for i = 1:length(t_shift)
            % Plot interdecile range
            % FOR HORIZONTAL BOX PLOTS
            bSub = bRow & time_shift == t_shift(i);
            temp = stat_UH(bSub, :);
            n_temp = mode(n_samp(bSub));
            % Plot interdecile
            x = temp(1, [2,6]);
            y = t_shift(i)*[1, 1];
            plot(ax_temp, x, y, '-', 'Color', cBox, 'LineWidth', 1);

            % Plot interquartile
            x = temp(1, [3,5]);
            plot(ax_temp, x, y, '-', 'Color', cBox, 'LineWidth', 4);

            % Plot median
            x = temp(1, 4) * [1, 1];
            y = y + 0.3*[-1, 1];
            plot(ax_temp, x, y, '-', 'Color', [1,1,1], 'LineWidth', 1);
            
            % Plot heavy 0 line
            plot(ax_temp, [0, 0], [-16, 16], ':', 'Color', 0.2*[1,1,1], 'LineWidth', 1);
        end
end


% Column 4: Uh_m - Uh_obs: Averaging
c_model = [20, 20, 20;...
           44,17,87;...
           17,23,87;...
           17,87,74]/256 + 0.2*[1,1,1];

model_opt = [0, 60, 120, 300];

for k = 1:3   % Each device
    if k == 1    
        % Floats
        dev_avg_opt = [0, 30, 60, 120];      
    elseif k == 2
        % SWIFTS
        dev_avg_opt = [0, 60, 120];
    elseif k == 3
        % Station-keeping
        dev_avg_opt = [0, 30, 60, 120, 300];
    end

    n_data_rows = length(model_opt)*(length(dev_avg_opt)+1);
    ax_temp = ax(k,4);
    ax_temp.NextPlot = 'add';

    for m = 1:4  % Each model averaging
        % Model averaging   
        for n = 1:length(dev_avg_opt)  % Each device averaging
        bRow = dev_type == k-1 & survey_id == iSurvey & dev_avg == dev_avg_opt(n) & model_avg == model_opt(m) & time_shift == 0;
        if m == 1
        bRow = dev_type == k-1 & survey_id == iSurvey & dev_avg == dev_avg_opt(n) & model_avg == model_opt(m) & time_shift == 100;
        end
        temp = stat_UH(bRow, :);
        %n_temp = mode(n_samp(bRow));
        cBox = c_model(m,:);

        % Plot interdecile
        x = temp(1, [2,6]);
        y = (n_data_rows-((m-1)*(length(dev_avg_opt)+1)+n))*[1, 1];
        plot(ax_temp, x, y, '-', 'Color', cBox, 'LineWidth', 1);

        % Plot interquartile
        x = temp(1, [3,5]);
        plot(ax_temp, x, y, '-', 'Color', cBox, 'LineWidth', 4);

        % Plot median
        x = temp(1, 4) * [1, 1];
        y = y + 0.3*[-1, 1];
        plot(ax_temp, x, y, '-', 'Color', [1,1,1], 'LineWidth', 1);
        end
    end
end

set(ax(:,3), 'DataAspectRatio', [1,30,1],...
    'XLim', [-0.5, 0.5], 'YLim', [-16 16],...
    'XTick', [-5:0.25:0.5],  'XGrid', 'on',...
    'YTick', [-15 -10 -5 0 5 10 15],'YGrid','on');
set(ax(:,4), 'DataAspectRatio', [1,20,1], ...
     'XLim', [-0.5, 0.5],...
    'XTick', [-5:0.25:0.5],  'XGrid', 'on');
for i = 1:3
    ax(i,3).Position([2,4]) = ax(i,1).Position([2,4]);
    ax(i,4).Position([2,4]) = ax(i,1).Position([2,4]);
end

% ROW LABELS
ar1 = text(ax(1,1), ax(1,1).XLim(1)-0.6*diff(ax(1,1).XLim), ax(1,1).YLim(1)+0.5*diff(ax(1,1).YLim), '\muFloats',...
    'FontSize',12,'FontWeight','bold',...
    'HorizontalAlignment', 'center',...
    'Rotation', 90);

ar2 = text(ax(2,1), ax(2,1).XLim(1)-0.6*diff(ax(2,1).XLim), ax(2,1).YLim(1)+0.5*diff(ax(2,1).YLim), 'SWIFTS',...
    'FontSize',12,'FontWeight','bold',...
    'HorizontalAlignment', 'center',...
    'Rotation', 90);

ar3 = text(ax(3,1), ax(3,1).XLim(1)-0.6*diff(ax(3,1).XLim), ax(3,1).YLim(1)+0.5*diff(ax(3,1).YLim), 'Station-Keeping',...
    'FontSize',12,'FontWeight','bold',...
    'HorizontalAlignment', 'center',...
    'Rotation', 90);

for i = 1:3
xlabel(ax(i,1), 'U_x [m/s]');
ylabel(ax(i,1), 'U_y [m/s]');

xlabel(ax(i,2), 'Model [m/s]');
ylabel(ax(i,2), 'Obsv. [m/s]');

xlabel(ax(i,3), 'Diff. [m/s]');
ylabel(ax(i,3), 'Shift [min.]');

xlabel(ax(i,4), 'Diff. [m/s]');
yticks(ax(i,4),'');
end

for i = 1:3
    ax(i,1).Position(1) = 0.12;
    ax(i,2).Position(1) = 0.35;
    ax(i,3).Position(1) = 0.6;
    ax(i,4).Position(1) = 0.8;
end

% COLUMN LABELS
ac1 = annotation(f2, 'textbox', 'Position',...
    [ax(1,1).Position(1), 0.9, ax(1,1).Position(3), 0.1],...
    'string', 'U_x and U_y',...
    'FontSize',12,'FontWeight','bold',...
    'HorizontalAlignment', 'center',...
    'EdgeColor','none');

ac2 = annotation(f2, 'textbox', 'Position',...
    [ax(1,2).Position(1), 0.9, ax(1,2).Position(3), 0.1],...
    'string', 'U_h',...
    'FontSize',12,'FontWeight','bold',...
    'HorizontalAlignment', 'center',...
    'EdgeColor','none');

ac3 = annotation(f2, 'textbox', 'Position',...
    [ax(1,3).Position(1)-0.1, 0.9, ax(1,3).Position(3)+0.2, 0.1],...
    'string', 'Time Shifted',...
    'FontSize',12,'FontWeight','bold',...
    'HorizontalAlignment', 'center',...
    'EdgeColor','none');

ac4 = annotation(f2, 'textbox', 'Position',...
    [ax(1,4).Position(1), 0.9, ax(1,4).Position(3), 0.1],...
    'string', 'Averaged',...
    'FontSize',12,'FontWeight','bold',...
    'HorizontalAlignment', 'center',...
    'EdgeColor','none');

ac5 = annotation(f2, 'textbox', 'Position',...
    [0, 0.95, 0.05, 0.05],...
    'string', meta.testNames{iSurvey},...
    'FontSize',12,'FontWeight','bold',...
    'Color', 0.3*[1,1,1],...
    'HorizontalAlignment', 'center',...
    'VerticalAlignment','middle',...
    'EdgeColor','none');

% Axes management

if iSurvey == 5
    set(ax(3,[3,4]), 'XLim',  [-0.5 0.75], 'XTick', [-0.5:0.5:1]);
    set(ax([1,2],[3,4]), 'XLim',  [-0.75 0.75], 'XTick', [-0.5:0.5:0.5]);
    set(ax([1,2],[3,4]), 'DataAspectRatioMode','auto');
end

if iSurvey == 6
    set(ax(3,[3,4]), 'XLim',  [-0.25 0.75], 'XTick', [-0.25:0.25:0.75]);
    set(ax([1,2],[3,4]), 'XLim',  [-0.5 1], 'XTick', [0:0.5:1]);
    set(ax([1,2],[3,4]), 'DataAspectRatioMode','auto');
end

if iSurvey == 7
    set(ax(:,2), 'XLim',  [0, 2.75], 'YLim', [0, 2.75]);
    set(ax(3,[3,4]), 'XLim',  [-0.25 0.75], 'XTick', [-0.25:0.25:0.75]);
    set(ax([1,2],[3,4]), 'XLim',  [-0.5 1], 'XTick', [0:0.5:1]);
    set(ax([1,2],[3,4]), 'DataAspectRatioMode','auto');
end

if iSurvey == 8
    set(ax(:,2), 'XLim',  [0, 2.75], 'YLim', [0, 2.75]);
    set(ax(3,[3,4]), 'XLim',  [-1 0.5], 'XTick', [-1:0.5:1]);
    set(ax([1,2],[3,4]), 'XLim',  [-1 0.5], 'XTick', [-1:0.5:1]);
    set(ax(:,[3,4]), 'DataAspectRatioMode','auto');
end

if iSurvey == 9
    set(ax(:,2), 'XLim',  [0, 2.75], 'YLim', [0, 2.75]);
    set(ax(3,[3,4]), 'XLim',  [-1 0.5], 'XTick', [-1:0.5:1]);
    set(ax([1,2],[3,4]), 'XLim',  [-1 0.5], 'XTick', [-1:0.5:1]);
    set(ax(:,[3,4]), 'DataAspectRatioMode','auto');
end

grid(ax, 'on');
cm = lines(7);
p1 = plot(ax(1,1), nan, nan, '.','MarkerSize', 16, 'Color',cm(1,:));
p2 = plot(ax(1,1), nan, nan, '.','MarkerSize', 16, 'Color',cm(2,:));
lg = legend(ax(1,1), [p1, p2], {'Obsv.', 'Model'}, 'Box','off',...
    'Orientation','vertical', 'Position',[0.12 0.882 0.08 0.0571],...
    'FontSize', 8);

if bit_print == 1
print(f2, fullfile(dirFig,sprintf('%s_Compiled_Stat', meta.testNames{iSurvey})), '-djpeg','-r600');
end

%% FIGURE 17/18: Examining spatial differences
% Plots, in 3D space, the horizontal velocity magnitude differences

f1 = figure(1); clf(f1); f1.Color = 'w';
ax1 = gca;
set(ax1,'NextPlot','add');
f1.Units = 'inches'; 
%f1.Position = [1,1,6.5,7];
f1.Position = [1.0000    2.6583    6.5000    5.3417];
xOffset = 0; 
yOffset = 0;

% Plot the site
ax = ax1;
pLandE = patch(ax, xCLeast-xOffset, yCLeast-yOffset, 0.9*[1,1,1],'EdgeColor','none');
pLandW = patch(ax, xCLwest-xOffset, yCLwest-yOffset, 0.9*[1,1,1],'EdgeColor','none');
  
set(ax1, 'XLim', [-1000, 1500], 'YLim', [-1250, 2000]);
%set(ax1(3:4), 'XLim', [-500, 1250], 'YLim', [-1250, 1000]);
set(ax1,'XTick', -1000:500:2000, 'YTick', -1000:500:2000,...
    'XGrid','on', 'YGrid', 'on',...
    'DataAspectRatio',[1,1,1]);

xlabel(ax1, 'Easting [m]');
ylabel(ax1, 'Northing [m]');
%print(f1, fullfile(dirFig,'model_coord_2_local'),'-djpeg','-r600');

plot(ax1, locNOAA(1), locNOAA(2), 'k.', 'MarkerSize', 20);

bSub = SK.testID == testID;
x = SK.X(bSub);
y = SK.Y(bSub);
z = SK.Z(bSub);
ux = SK.UX(bSub);
uy = SK.UY(bSub);
uxq = SK.UXQ(bSub);
uyq = SK.UYQ(bSub);   
X = x(:); Y = y(:); Z = z(:); UX = ux(:); UY = uy(:); UXQ = uxq(:); UYQ = uyq(:);
devID = 3+0*x(:);

bSub = SWIFT.testID == testID;
x = SWIFT.X(bSub);
y = SWIFT.Y(bSub);
z = SWIFT.Z(bSub);
ux = SWIFT.UX(bSub);
uy = SWIFT.UY(bSub);
uxq = SWIFT.UXQ(bSub);
uyq = SWIFT.UYQ(bSub); 
X = [X; x(:)]; Y = [Y; y(:)]; Z = [Z; z(:)];
UX = [UX; ux(:)]; UY = [UY; uy(:)]; UXQ = [UXQ; uxq(:)]; UYQ = [UYQ; uyq(:)]; 
devID = [devID; 2+0*x(:)];

x = FS(testID).X;
y = FS(testID).Y;
z = FS(testID).Z;
ux = FS(testID).UX;
uy = FS(testID).UY;
uxq = FS(testID).UXQ;
uyq = FS(testID).UYQ; 
X = [X; x(:)]; Y = [Y; y(:)]; Z = [Z; z(:)];
UX = [UX; ux(:)]; UY = [UY; uy(:)]; UXQ = [UXQ; uxq(:)]; UYQ = [UYQ; uyq(:)]; 
devID = [devID; 1+0*x(:)];

x = X; y = Y; z = Z; ux = UX; uy = UY; uxq = UXQ; uyq = UYQ; 

% Difference in magnitude
c_map_choice = 'diff_mag'; % ['diff_mag' | 'diff_dir']
switch c_map_choice
    case 'diff_mag'
        cdata = (uxq .^2 + uyq .^2) .^0.5 - (ux .^2 + uy .^2) .^0.5;
        data_mapping = 'U_{H_{model}} - U_{H_{obs.}}';
        cb_label = 'Difference [m/s]';
        c_lim = [-0.5, 0.5];
    case 'diff_dir'
        cdata = angdiff(atan2(uyq,uxq),atan2(uy, ux)).*180/pi;
        data_mapping = 'Difference in horizontal direction';
        cb_label = 'Difference [^o]';
        c_lim =[-30, 30];
end

dev_ind = 4;
dev_name = {'Float', 'SWIFT', 'Station-keeping',  'All'};
bSub = devID == 3;
if dev_ind == 4
    bSub = devID < 4;
end
sc = scatter3(ax1, x(bSub), y(bSub), z(bSub), 12, cdata(bSub), 'filled');
caxis(ax1, c_lim);
cm = colormap(cmocean('balance','pivot',0));
cb = colorbar(ax1, 'south');
cb.Position = [0.1471    0.08    0.7408    0.0414];
ylabel(cb, cb_label);

ax1.ZLabel.String = 'Depth [m]';

ax1.DataAspectRatio = [100,100,1];
ax1.ZDir = 'rev';
title(ax1, sprintf('Survey %s  |  Platform(s): %s   |   %s',...
   meta.testNames{testID}, dev_name{dev_ind}, data_mapping));

% Can we find locations where SWIFT and float data are in the same
% locations (space and time)?

pLandE.EdgeColor = 0.2*[1,1,1];
pLandE.FaceColor = 'none';
pLandE.LineWidth = 01;

pLandW.EdgeColor = 0.2*[1,1,1];
pLandE.EdgeColor = 0.2*[1,1,1];

% For Ebb
ax1.ZLim = [-0.5 10];
ax1.View = [32.7569 25.9990];
ax1.CameraTarget = [250 375 4.7500];
ax1.CameraPosition = [1.0542e+04 -1.5622e+04 -88.0195];
ax1.CameraViewAngle = 10.4240;
ax1.Position = [0.00    0.050    1   1];

if bit_print == 1
filename = sprintf('spatial_diff_survey_%s', meta.testNames{testID});
print(f1, fullfile(dirFig,filename), '-djpeg','-r600');
end
%% Supplemental - Plot Comparisons
if 0 == 1

bit_print = 1;
for iSurvey = 2

% Get subset of survey data
ux_f = FS(iSurvey).UX;
uy_f = FS(iSurvey).UY;
ux_fq = FS(iSurvey).UXQ;
uy_fq = FS(iSurvey).UYQ;

bSub = SWIFT.testID == iSurvey;
ux_s = SWIFT.UX(bSub);
uy_s = SWIFT.UY(bSub);
ux_sq = SWIFT.UXQ(bSub);
uy_sq = SWIFT.UYQ(bSub);

bSub = SK.testID == iSurvey;
ux_sk = SK.UX(bSub);
uy_sk = SK.UY(bSub);
ux_skq = SK.UXQ(bSub);
uy_skq = SK.UYQ(bSub);

allData = [ux_sk; ux_s; ux_f(:); uy_sk; uy_s; uy_f(:);...
           ux_skq; ux_sq; ux_fq(:); uy_skq; uy_sq; uy_fq(:)];
          
xL = [floor(min(allData)*10)/10, ceil(max(allData)*10)/10];

f1 = figure(1); clf(f1);
ax = buildAxes(f1, 3, 3);


set(ax, 'XLim', xL, 'YLim', xL);
set(ax(3,1:3),'XLim',[0, max(abs(xL))], 'YLim',[0, max(abs(xL))]); 
for i = 1:numel(ax)
    plot(ax(i), ax(i).XLim, ax(i).YLim, '-', 'Color', 0.2*[1,1,1]);
end
mSize = 3;
set(ax, 'NextPlot','add','ColorOrderIndex',1,'DataAspectRatio',[1,1,1]);
plot(ax(1,1), ux_sk, ux_skq, '.', 'MarkerSize', mSize);
plot(ax(2,1), uy_sk, uy_skq, '.', 'MarkerSize', mSize);
plot(ax(3,1), sqrt(ux_sk.^2+uy_sk.^2), sqrt(ux_skq.^2+uy_skq.^2), '.',...
    'MarkerSize', mSize);

plot(ax(1,2), ux_s, ux_sq, '.', 'MarkerSize', mSize);
plot(ax(2,2), uy_s, uy_sq, '.', 'MarkerSize', mSize);
plot(ax(3,2), sqrt(ux_s.^2+uy_s.^2), sqrt(ux_sq.^2+uy_sq.^2), '.',...
    'MarkerSize', mSize);

plot(ax(1,3), ux_f(:), ux_fq(:), '.', 'MarkerSize', mSize);
plot(ax(2,3), uy_f(:), uy_fq(:), '.', 'MarkerSize', mSize);
plot(ax(3,3), sqrt(ux_f(:).^2+uy_f(:).^2), sqrt(ux_fq(:).^2+uy_fq(:).^2), '.',...
    'MarkerSize', mSize);

ylabel(ax(1,1), 'UX');
ylabel(ax(2,1), 'SIMULATION\newline         UY')
ylabel(ax(3,1), '|U|');
xlabel(ax(3,2), 'OBSERVATION');
title(ax(1,1), 'Station Keeping');
title(ax(1,2), 'SWIFT');
title(ax(1,3), '\muFloat');
set(ax, 'Color', 0.96*[1,1,1]);
f1.Color = 'w';

set(ax(:,2:3),'YTickLabels','');
set(ax(1:2,:),'XTickLabels','');

set(ax(1:2,1:3), 'XLim', xL, 'YLim', xL);
grid(ax, 'on');
annotation(f1, 'textbox', 'Position', [0.01, 0.9, 0.1, 0.1],...
    'String',meta.testNames{iSurvey},...
    'EdgeColor', 'none', 'HorizontalAlignment','left',...
    'FontSize',14,'FontWeight','bold');

for i = 1:numel(ax)
    ax(i).Position([3,4]) = 0.28*[1, f1.Position(4)/f1.Position(3)];
    ax(i).Position(1) = ax(i).Position(1)-0.04;
    ax(i).Position(2) = ax(i).Position(2)-0.04;
end
set(ax(3,1:3),'XLim',[0, max(abs(xL))], 'YLim',[0, max(abs(xL))]); 

if bit_print ==1
print(f1, fullfile(dirFig,sprintf('%s_obs_v_sim_corr', meta.testNames{iSurvey})), '-djpeg','-r600');
end



%{
annotation(f1,'textbox',...
    [0.1 0.96 0.8 0.04],...
    'String','Observed vs. Simulated Velocity Components',...
    'FitBoxToText','off','HorizontalAlignment','center',...
    'FontWeight','bold','EdgeColor','none');
%}
    
%
f2 = figure(2); clf(f2);
%ax = buildAxes(f2, 1, 3); 
ax(1) = polaraxes('Position', [0.1, 0.1, 0.2, 0.8]);
ax(2) = polaraxes('Position', [0.4, 0.1, 0.2, 0.8]);
ax(3) = polaraxes('Position', [0.7, 0.1, 0.2, 0.8]);
set(ax, 'NextPlot','add');

mSize = 1;
theta = atan2(uy_sk, ux_sk);
rho = sqrt(ux_sk.^2+uy_sk.^2);
theta_q = atan2(uy_skq, ux_skq);
rho_q = sqrt(ux_skq.^2+uy_skq.^2);
polarscatter(ax(1), theta, rho, mSize, 'filled');
polarscatter(ax(1), theta_q, rho_q, mSize, 'filled');

theta = atan2(uy_s, ux_s);
rho = sqrt(uy_s.^2+ux_s.^2);
theta_q = atan2(uy_sq, ux_sq);
rho_q = sqrt(ux_sq.^2+uy_sq.^2);
polarscatter(ax(2), theta, rho, mSize, 'filled');
polarscatter(ax(2), theta_q, rho_q, mSize, 'filled');

theta = atan2(uy_f(:), ux_f(:));
rho = sqrt(uy_f(:).^2+ux_f(:).^2);
theta_q = atan2(uy_fq(:), ux_fq(:));
rho_q = sqrt(ux_fq(:).^2+uy_fq(:).^2);
polarscatter(ax(3), theta, rho, mSize, 'filled');
polarscatter(ax(3), theta_q, rho_q, mSize, 'filled');

title(ax(1), 'Station Keeping');
title(ax(2), 'SWIFT');
title(ax(3), '\muFloat');

%
f3 = figure(3); clf(f3); f3.Color = 'w';
mSize = 3;
ax = buildAxes(f3, 1, 3); set(ax, 'NextPlot','add');
plot(ax(1), ux_sk, uy_sk, '.','MarkerSize',mSize);
plot(ax(1), ux_skq, uy_skq, '.','MarkerSize',mSize);

plot(ax(2), ux_s, uy_s, '.','MarkerSize',mSize);
plot(ax(2), ux_sq, uy_sq, '.','MarkerSize',mSize);

plot(ax(3), ux_f(:), uy_f(:), '.','MarkerSize',mSize);
plot(ax(3), ux_fq(:), uy_fq(:), '.','MarkerSize',mSize);

title(ax(1), 'Station Keeping');
title(ax(2), 'SWIFT');
title(ax(3), '\muFloat');
grid(ax, 'on');
xlabel(ax(2), 'U_x [m/s]');
ylabel(ax(1), 'U_y [m/s]');

cLines = lines(7);
pl(1) = plot(ax(2),nan, nan, '.','MarkerSize',12,'Color',cLines(1,:));
pl(2) = plot(ax(2),nan, nan, '.','MarkerSize',12,'Color',cLines(2,:));

lg = legend(ax(2),pl, {'Obs','Sim'},'Box','Off');
lg.Position([1,2]) = [0.39, 0.75];
linkaxes(ax,'xy');
set(ax, 'Color', 0.96*[1,1,1]);

annotation(f3, 'textbox', 'Position', [0.01, 0.9, 0.1, 0.1],...
    'String',meta.testNames{iSurvey},...
    'EdgeColor', 'none', 'HorizontalAlignment','left',...
    'FontSize',14,'FontWeight','bold');

if bit_print == 1
print(f3, fullfile(dirFig,sprintf('%s_obs_v_sim_vel_cloud',meta.testNames{iSurvey})),'-djpeg','-r600');
end
end



end
















