mission = 24403401;

%

 
filepath = ['\\atlas.shore.mbari.org\ProjectLibrary' ...
    '\901805_Coastal_Biogeochemical_Sensing\Spray_Data' ...
    '\',num2str(mission),'\sat_files\',num2str(mission),'_sat.mat'];
load(filepath,'s')
wmclose("all");
fig1 = figure;
% tl = tiledlayout(2,2);
%% MBARI
% mosslat = 36.80594819607352;
% mosslon = -121.786377002975;


%% CalCOFI Line 66
CalCOFI = readtable("CalCOFIStationOrder.csv");
CalCOFI67 = CalCOFI(CalCOFI.Line == 66.7,["Sta","Lat_dec_","Lon_dec_"]);
CalCOFI67.Properties.VariableNames = ["Station","Lat","Lon"];
CalCOFI67Description = strsplit(num2str((CalCOFI67.Station')));

%% M1 and C1
M1Lat = 36.75901;
M1Lon = -122.02019;
% need to add C1

%% Spray 034
SprayLatHist = s.lat(end-4:end-1);
SprayLonHist = s.lon(end-4:end-1);
SprayLat = s.lat(end);
SprayLon = s.lon(end);

%% Datetimes
date = datetime(datevec(s.sdn_(2,:))); 
date.TimeZone = "UTC"; % Set timezone for the datetime var type
date.TimeZone = "America/Los_Angeles"; % change timezone to local
SprayHistDescription = cellstr(date(end-4:end-1))';
SprayDescription = cellstr(date(end))';
% nexttile([1 2])
clf
%% Create a geographic axes
ax = geoaxes('Basemap','satellite');
%% Hold the current plot
hold(ax, 'on');
%% Plot CalCOFI Line 66
geoplot(ax, CalCOFI67.Lat, CalCOFI67.Lon, '.', 'Color', 'yellow', 'DisplayName', 'CalCOFI Line 67', 'MarkerSize', 10);
%% Plot M1
geoplot(ax, M1Lat, M1Lon, 'Color', [1 0.5 0], 'DisplayName', 'MBARI Mooring', 'MarkerSize', 10,'Marker','diamond');
%% Plot Spray 034 Location History
geoplot(ax, SprayLatHist, SprayLonHist, '.', 'Color', [0 0.75 0.75], 'DisplayName', 'Spray 034 Location History', 'MarkerSize', 10);
%% Plot Spray 034 Location
geoplot(ax, SprayLat, SprayLon, '*', 'Color', [1 0.5 0], 'DisplayName', 'Spray 034 Location', 'MarkerSize', 20);
%% Add a legend
legend(ax, 'Location', 'northwest');
%% Release the hold on the plot
hold(ax, 'off');

% wm = webmap("National Geographic Map");
% wmmarker(CalCOFI67.Lat,CalCOFI67.Lon,"Description",CalCOFI67Description,"Icon","CalCOFILogo.png","OverlayName",'CalCOFI Line 67')
% wmmarker(M1Lat,M1Lon,"Description",'MBARI Mooring',"Icon","M1_buoyriding_72.png","OverlayName",'M1','IconScale',2)
% wmmarker(SprayLatHist,SprayLonHist,"Icon","spray-glider.png","OverlayName","Spray 034 Location History","Alpha",.2,"IconScale",2,"Description",SprayHistDescription)
% wmmarker(SprayLat,SprayLon,"Icon","spray-glider.png","OverlayName","Spray 034 Location","Alpha",1,"IconScale",2,"Description",SprayDescription)

%

dielstart = 516; % start profile of diel cycle measurements
dielend = length(s.sdn);
dielmax = length(s.sdn);
% Get vars
s.o2satconc = reshape(calcO2sat(s.tc(:), s.psal(:)), size(s.pres));
s.o2anom = s.doxy - s.o2satconc;

idiel = dielstart:dielend;

date = datetime(datevec(s.sdn(idiel))); 
date.TimeZone = "UTC"; % Set timezone for the datetime var type
date.TimeZone = "America/Los_Angeles"; % change timezone to local

dielstart = 517;
dielend = 584;
% Get vars
s.o2satconc = reshape(calcO2sat(s.tc(:), s.psal(:)), size(s.pres));
s.o2anom = s.doxy - s.o2satconc;

idiel = dielstart:dielend;

date = datetime(datevec(s.sdn(idiel))); 
date.TimeZone = "UTC"; % Set timezone for the datetime var type
date.TimeZone = "America/Los_Angeles"; % change timezone to local
%
depthavg = mean(s.depth(:,idiel),2);
depthmax = max(depthavg);
depth1 = 1;
d1 = find(depthavg > depth1-1 & depthavg < depth1+1, 1);
depth2 = 3;
d2 = find(depthavg > depth2-1 & depthavg < depth2+1, 1);
depth3 = 20;
d3 = find(depthavg > depth3-1 & depthavg < depth3+1, 1);
%
fig2 = figure;
tl = tiledlayout(1,1,"TileSpacing","tight");
title(tl, [datestr(date(1)), ' to ', datestr(date(end))])
subtitle(tl,[num2str(length(s.o2anom(1,idiel))), ' profiles'])
ylabel(tl, 'O_2 Sat Anom [\mumol/kg]')
xlabel(tl, 'Time Local (PST)')

nexttile
for iupper = [d1, d2, d3]
    meano2anom = mean(s.o2anom(iupper,idiel),1);
    plot(date,meano2anom,'.')
    hold on
end
title('O2 Anom Diel')
% xticks
datetickstart = date(1);
datetickstart.Hour = 0;
datetickstart.Minute = 0;
datetickstart.Second = 0;
datetickend = date(end);
datetickend.Hour = 0;
datetickend.Minute = 0;
datetickend.Second = 0;

if numel(datetickstart:.25:datetickend+1) < 15
    xticks(datetickstart:.25:datetickend+1)
else
    xticks(datetickstart:.5:datetickend+1)
end
xtickformat('HH:mm')
xtickangle(45)
% xlabel(tl,'Local time (hours:minutes)')
legend(num2str(depthavg(d1)),num2str(depthavg(d2)),num2str(depthavg(d3)), Location='northeastoutside')
grid on

%

fig3 = figure;
tl = tiledlayout(1,1,"TileSpacing","tight");

%% Thermocline depth (max dt/dz) 
temp = s.tc(:,idiel);
doxy = s.doxy(:,idiel);
sal = s.psal(:,idiel);
pH = s.pH25atm(:,idiel);
depth = s.depth(:,idiel);
thermdepth = NaN(1,length(idiel));
for i = 1:length(idiel)
    [~, idx] = nanmax(abs(diff(temp(:,i))));
    thermdepth(i) = depth(idx);
end

nexttile 
plot(date,thermdepth,'.')
title('Max depth of thermocline [max(diff(temp)]')
ylabel('thermocline depth [m]')
set(gca,'YDIR','reverse')
if numel(datetickstart:.25:datetickend+1) < 15
    xticks(datetickstart:.25:datetickend+1)
else
    xticks(datetickstart:.5:datetickend+1)
end
xtickformat('HH:mm')
xtickangle(45)
grid on

%

fig4 = figure;
tl = tiledlayout(1,4,"TileSpacing","tight");

title(tl,['Diel Profile ', num2str(dielend - dielstart + 1)])
nexttile
plot(temp(:,(dielend - dielstart + 1)),depth(:,dielend - dielstart + 1))
title('temp')
set(gca,'YDIR','reverse')

nexttile
plot(doxy(:,(dielend - dielstart + 1)),depth(:,dielend - dielstart + 1))
title('doxy')
set(gca,'YDIR','reverse')

nexttile
plot(sal(:,(dielend - dielstart + 1)),depth(:,dielend - dielstart + 1))
title('sal')
set(gca,'YDIR','reverse')

nexttile
plot(pH(:,(dielend - dielstart + 1)),depth(:,dielend - dielstart + 1))
title('pH 25 atm')
set(gca,'YDIR','reverse')

%

%contour plots of all the vars
iupper = 1:10;
meano2anom = mean(s.o2anom(1:5,idiel),1);

% s.dic_canb(s.dic_canb_QC ~= 1) = [];
DIC = s.dic_canb(:,idiel);
fig5 = figure;
tl = tiledlayout(3,1,"TileSpacing","tight");

nexttile
plot(date,meano2anom,'.')
ylabel('meanO2Anom')
hold on
yyaxis right
ylabel('DIC')
plot(date,nanmean(DIC),'.')
xticklabels([])
title('mean of upper 10 m')

nexttile
pcolor(date,s.depth(iupper,idiel(1)),s.o2anom(iupper,idiel))
xticklabels([])
shading interp
cb1 = colorbar;
ylabel(cb1,'O_2 anom')
ylabel('depth')
if numel(date(1):.25:date(end)+1) < 15
    xticks(date(1):.25:date(end)+1)
else
    xticks(date(1):.5:date(end)+1)
end
set(gca,'YDIR','reverse')

nexttile
pcolor(date,s.depth(iupper,idiel(1)),s.dic_canb(iupper,idiel))
shading interp
cb2 = colorbar;
ylabel(cb2,'DIC')

xtickformat('HH:mm')
xtickangle(45)
xlabel(tl,'Local time (hours:minutes)')
ylabel('depth')
set(gca,'YDIR','reverse')