%% **--**--**--**--**--**--** Section (1) **--**--**--**--**--**--**--**--
% example plots of daily maps:
% a. HFR-derived surface currents
% b. eddy locations
% c. retention times

clear; clc;
% use ncdisp to look into .nc file
ncdisp("RB_2023maps.nc");
% load in lat/lon and time
lat = ncread("RB_2023maps.nc","lat");
lon = ncread("RB_2023maps.nc","lon");
t = ncread("RB_2023maps.nc","t");

%% Plotting mean/total values
% load in each variable and take time mean/sum
% surface currents
u = mean(ncread("RB_2023maps.nc","u",[1 1 1],[length(lat) length(lon) length(t)]),3,'omitnan');
v = mean(ncread("RB_2023maps.nc","v",[1 1 1],[length(lat) length(lon) length(t)]),3,'omitnan');
% eddy position
eddypos = ncread("RB_2023maps.nc","eddypos",[1 1 1],[length(lat) length(lon) length(t)]);
ce = nan*eddypos; ce(eddypos > 0) = 1; ce = sum(ce,3,'omitnan')/length(t); ce(ce==0) = nan; % index seperately for plotting
ae = nan*eddypos; ae(eddypos < 0) = -1; ae = sum(ae,3,'omitnan')/length(t); ae(ae==0) = nan; 
% retention time
rt20 = mean(ncread("RB_2023maps.nc","rt20",[1 1 1],[length(lat) length(lon) length(t)]),3,'omitnan');

% plotting
% surface currents
figure('position',[1 1 800 500]); 
quiver(lon(1:2:end),lat(1:2:end),u(1:2:end,1:2:end),v(1:2:end,1:2:end),'k','linew',1);
axis equal; xlim([-120.75 -119]); ylim([33.75 34.75]);
title('surface currents'); set(gca,'linew',2,'fontsize',20);
% eddy positions
figure('position',[1 1 800 500]); 
ce(ce <= 0.1) = nan; ae(ae >= -0.1) = nan;
pcolor(lon,lat,ce*100); 
hold on; pcolor(lon,lat,ae*100); shading flat;
colorbar; colormap(redblue(11)); clim([-50 50]);
axis equal; xlim([-120.75 -119]); ylim([33.75 34.75]);
title('eddy positions (% time eddy present)'); set(gca,'linew',2,'fontsize',20);
% retention time
figure('position',[1 1 800 500]); 
pcolor(lon,lat,rt20/24); shading interp; colorbar; colormap(jet); clim([0 5]);
axis equal; xlim([-120.75 -119]); ylim([33.75 34.75]);
title('retention time (days)'); set(gca,'linew',2,'fontsize',20);

%% Plotting daily maps
time1 = datetime(0000,01,00) + t;
day = datetime(2019,04,23); % *** change date here ***
ind = find(time1 == day); 

% surface currents
u = ncread("RB_2023maps.nc","u",[1 1 ind],[length(lat) length(lon) 1]);
v = ncread("RB_2023maps.nc","v",[1 1 ind],[length(lat) length(lon) 1]);
% eddy position
eddypos = ncread("RB_2023maps.nc","eddypos",[1 1 ind],[length(lat) length(lon) 1]);
% retention time
rt20 = ncread("RB_2023maps.nc","rt20",[1 1 ind],[length(lat) length(lon) 1]);

% plotting
% surface currents
figure('position',[1 1 800 500]);
quiver(lon(1:2:end),lat(1:2:end),u(1:2:end,1:2:end),v(1:2:end,1:2:end),'k','linew',1);
axis equal; xlim([-120.75 -119]); ylim([33.75 34.75]);
title(['surface currents; ',datestr(day)]);  set(gca,'linew',2,'fontsize',20);
% eddy positions
figure('position',[1 1 800 500]);
pcolor(lon,lat,eddypos); shading flat; c = colorbar; colormap(redblue(5)); clim([-2.5 2.5]);
axis equal; xlim([-120.75 -119]); ylim([33.75 34.75]); c.Ticks = [-2 -1 1 2]; 
c.TickLabels = {'AE core','AE extent','CE extent','CE core'};
title(['eddy positions; ',datestr(day)]);  set(gca,'linew',2,'fontsize',20);
% retention time
figure('position',[1 1 800 500]);
pcolor(lon,lat,rt20/24); shading interp; colorbar; colormap(jet); clim([0 max(rt20/24,[],'all')]);
axis equal; xlim([-120.75 -119]); ylim([33.75 34.75]);
title(['retention time (hours); ',datestr(day)]);  set(gca,'linew',2,'fontsize',20);

%% **--**--**--**--**--**--** Section (2) **--**--**--**--**--**--**--**--
% example plots of monthly timeseries:
% wind stress; dSSP; eof amplitude & phase; vorticity;
% eddy areas; eddy presence; regional retention times

clear; clc;
% use ncdisp to look into .nc file
ncdisp("RB_2023timeseries.nc");
time = datetime(0000,01,00) + ncread("RB_2023timeseries.nc","tm"); % monthly time vector
% choose time series to load and plot
tau54 = ncread("RB_2023timeseries.nc","tau54"); % wind stress at buoy 46054
RT_E = ncread("RB_2023timeseries.nc","RT_E"); % retention time in the eastern SBC

%% plotting
figure('position',[1 1 1200 500]);
subplot(2,1,1); % wind stress at buoy 46054 
plot(time,tau54,'b-','linew',2);
ylabel('{\tau}_5_4 (Pa)');
set(gca,'linew',2,'fontsize',20);
subplot(2,1,2); % retention time in the eastern SBC 
plot(time,RT_E,'r-','linew',2);
ylabel('RT_E (days)');
set(gca,'linew',2,'fontsize',20);

