%% Check simulation in same coordinates as field data
define_paths;

% Load meta data
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

D = dir(fullfile(dirSimData, 'Agate_*'));

T_SIM = [];
I_SIM = [];
for i = 1:length(D)
   filename = fullfile(dirSimData,D(i).name);
   % Load time-evolving data
   times = ncread(filename,'Times');
   times = times(1:19,:)'; 
   % convert to datetime
   ts = datetime(times,'InputFormat','uuuu-MM-dd''T''HH:mm:ss','TimeZone','UTC');
   T_SIM = [T_SIM; ts];
   I_SIM = [I_SIM; i*ones(size(ts))];
end   

% Define reference location. 
latRef = meta.locSK(3,1);
lonRef = meta.locSK(3,2);

% INITIALIZE SIMULATION PROPERTIES - TIME-INVARIENT VARIABLES
filename = fullfile(dirSimData,D(1).name);

% Load constant grid parameters
x_n = ncread(filename,'x');
y_n = ncread(filename,'y');
%lat_n = ncread(filename,'lat');
%lon_n = ncread(filename,'lon');
h_n = ncread(filename,'h');
siglay_n = ncread(filename,'siglay'); % centers [1x20]
siglev_n = ncread(filename,'siglev'); % edges   [1x21]

x_c = ncread(filename,'xc');
y_c = ncread(filename,'yc');
%lat_c = ncread(filename,'latc');
%lon_c = ncread(filename,'lonc');
h_c = ncread(filename,'h_center');
siglay_c = ncread(filename,'siglay_center'); % [-1, 0]
siglev_c = ncread(filename,'siglev_center');  
nLayer = size(siglay_c,2);

% [x, y] = geodetic2ned(lat_n,lon_n,0,latRef,lonRef,0,referenceEllipsoid('wgs84'));
% Convert reference location to zonal offset. 
[xRef, yRef] = ll2utm(latRef, lonRef);

% Sim node coordinates: local
xn_l = x_n-xRef;
yn_l = y_n-yRef;
% Sim center coordinates: local
xc_l = x_c-xRef;
yc_l = y_c-yRef;


% PLOT!
f1 = figure(1); clf(f1); f1.Color = 'w';
ax1 = gca;
set(ax1,'NextPlot','add');
f1.Units = 'inches'; 
f1.Position = [1,1,6.5,7];

xOffset = 0; 
yOffset = 0;

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]);

p1 = plot(ax1, xc_l, yc_l, 'r.'); hold on;
p2 = plot(ax1, xn_l, yn_l, 'g.');
title(ax1, 'Model coordinates shifted to observation reference frame (SK3)');
legend(ax1,[p1,p2], {'Center Points','Node Points'},'Location','NorthWest');
xlabel(ax1, 'Easting [m]');
ylabel(ax1, 'Northing [m]');
%print(f1, fullfile(dirFig,'model_coord_2_local'),'-djpeg','-r600');