%% Loading data
%
% Spray data from each mission is stored at:
% \\atlas.shore.mbari.org\ProjectLibrary\901805_Coastal_Biogeochemical_Sensing\Spray_Data
% Each mission folder has sat files and QC files. For science purposes you should use the file: mission_QC.mat
%

load('24403401_QC.mat','s') % load science data 's'

%
% This will load the science data into your workspace with the name s.
% Can also run load('24403401_QC.mat') to read in all of structs contained in the QC file such as calibration or engineering data.
%
%% Structs
%
% Working with a struct is similar to a table (dot indexing). Structs allow for many different data types to be stored within one variable in the workspace.
% The s struct contains parameters directly from the sensors and carbonate parameters that were calculated during the processing of the file.
% Example:
%

s.tc; % CTD temperature
s.pco2in; % Calculated pco2 insitu

%
%% Spray Data format
% Spray sensor data is formatted with measurements along each row and profiles across columns
% The final processed data are upcasts
% To link the measurements with depth or pressure readings use s.depth or s.pres
% Example 1: Plotting temperature vs depth for profile 1

temperature_profile_1 = s.tc(:,1); % Output will be all CTD temperature measurements for profile 1
depth_profile_1 = s.depth(:,1); % Output will be all depths for profile 1

figure(1)
clf
plot(temperature_profile_1,depth_profile_1)
title('Temperature')
subtitle('Profile 1')
xlabel('^o C')

%
% Example 2: Plotting temperature vs depth for profile 1. Only the first 10 measurements.

measurements = 1:10;
temperature_profile_1 = s.tc(measurements,1); % Output will be first 10 CTD temperature measurements for profile 1
depth_profile_1 = s.depth(measurements,1); % Output will be first 10 depths for profile 1

figure(2)
clf
plot(temperature_profile_1,depth_profile_1)
title('Temperature')
subtitle('Profile 1')
xlabel('^o C')

%
% Example 3: Plotting salinity (psal) vs depth for profiles 10 - 15 and measurements 1 - 100
%

profiles = 10:15;
measurements = 1:100;
temperature_profile_10_15 = s.psal(measurements,profiles); % Output will be first 10 CTD temperature measurements for profile 1
depth_profile_10_15 = s.depth(measurements,profiles); % Output will be first 10 depths for profile 1

figure(3)
clf
plot(temperature_profile_10_15,depth_profile_10_15)
title('Salinity')
subtitle('Profiles 10 - 15')
xlabel('^o C')

%
%% Data Operations
%

% Addition and subtraction
x = 10;
y = 10;
xy1 = x + y

%

% Multiplication and division
x = 10;
y = 15;
xy2 = x * y

x = [10 10 10 10 10];
y = 15;
xy3 = x * y

x = [10 10 10 10 10];
y = [15 16 17 18 19];
xy4 = x .* y % elementwise multiplication (x1 * y1, x2 * y2, ...)

%
% Example 1: Subtract Spray pHin from canyonB pHin for profile 120 and plot the residual vs depth then plot the location of Spray
%

profile = 120;
pHresidual = s.pHin_canb(:,profile) - s.pHin(:,profile);
depth = s.depth(:,profile);

figure(4)
clf
plot(pHresidual,depth)
title('pHin residuals (canB-pHin)')
subtitle('Profile 120')
xlabel('pH')
hold on; % Overlay other plots on top of the original
xline(0,'k--') % Dashed line at x = 0

%

% Now show the location for this profile
figure(5)
clf
geobasemap('satellite')
geoscatter(s.lat(profile),s.lon(profile),'or','filled')
geolimits([35.3 37.4],[-124.9 -121.6]) % set zoom limits
title('Spray location profile 120')

%
% Example 2: Calculate the pCO2 from profile 300 with CO2SYSv3 and plot the upper 5 surface values
%

profile = 300;
OutputMatrix = CO2SYSv3(s.ta_canb(:,profile), s.pHin(:,profile), 1, 3, s.psal(:,profile), s.tc(:,profile), s.tc(:,profile), 0, 0, s.sil_canb(:,profile), s.po4_canb(:,profile), 0, 0, 1, 10, 1, 2, 2);
% Read CO2SYSv3 for inputs and outputs. For example: TA = 1, pH = 3.
pCO2 = reshape(OutputMatrix(:,21), height(s.pres), 1); % 21 - pCO2 output (uatm)
pCO2(pCO2 == -999) = NaN; % CO2sysv3 outputs missing values as '-999'. Find all '-999' and replace with NaN.

%

measurements = 1:5;
figure(6)
clf
plot(pCO2(measurements),s.depth(measurements,profile))
title('pCO2')
subtitle(['Profile ', num2str(profile)])
xlabel('pCO2')
ylabel('depth (m)')