clear all
close all

addpath(genpath('/Users/Mau/Documents/MBARIprojects/w_ctd/cpf'));
addpath(genpath('/Users/Mau/Documents/MBARIprojects/CPF/CPFPROCESSING')); % parse_CPFmsg_v4

addpath(genpath('/Users/Mau/Documents/MATLABhelp_new')); % print_setup.m for saving plots
% cbdate, maximize

addpath(genpath('\\atlas\Chem\ARGO_PROCESSING\MFILES\FLOATS\'));
addpath(genpath('C:\Users\magdalena\Documents\ARGO_PROCESSING-master\'))
addpath(genpath('\\atlas\CarbonLab\magdalena\MATLABhelp_new\'));

%%  PURPOSE: Compute float's rising velocities from CTD data
% msg files on atlas for the CPF TankTests data are under
% /Volumes/CPF/CPF003/Data/TankTests/

%% Select path to data

% mission=['TankTests/2021-8-25T21-57-38/'];
% mission=['TankTests/2021-8-31T16-23-33/'];  % latest tank test
% mission=['TankTests/2021-9-3T18-48-26/'];
% mission='TankTests/2021-9-9T17-33-2/';
% mission='TankTests/2021-9-9T21-6-13/';
mission='ParagonTests/2021-9-16T16-33-30/';

mission='ParagonTests\MICRORIDER\OvernightMission\2021-10-21T18-34-19\';

maxp=100; % max depth for plotting

% Define root directory to data and plots
data_root=['/Volumes/CPF/CPF003/Data/']; % from local folder on macosx
dir_plots=['/Volumes/CPF/CPF003/Data/data_processing/plots/'];

data_root=['\\atlas\ProjectLibrary\CPF\CPF003\Data\'];% from local folder on windows
dir_plots=[data_root,'data_processing\plots\'];

path_data=[data_root,mission];

%% Create directory to save plots (if it doesn't exist already)
% % subdir=path_data(64:end); if path_data='/Users/Mau/Documents/..'
% subdir=path_data(26:end); % if /Volumes/CPF/CPF003/Data/
% % path_plots=['/Users/Mau/Documents/MBARIprojects/w_ctd/cpf/plots/',subdir];

path_plots=[dir_plots,mission];
if isfolder(path_plots)==0 & numel(dir(path_plots))<=2 % will create directory only if it is empty!
    mkdir(path_plots)
end

%%  List files for a given mission/folder

files=dir([path_data,'*.msg']);
np=length(files); % number of profiles

% Preallocate memory
P=nan(50,np); T=P; S=P; N=P; DEPTH=P;
W_ctd=P; W_ctd1=P;
d_axis=P;
SA=P; CT=P; RHO=P;
O2_air=nan(5,4,np);

Lat=nan(np,1);
Lon=nan(np,1);
Date=nan(np,1);

%% Process each profile from a given mission/folder

for p=1:np
    
filename=files(p).name;

% msg_file=['./test_files/',filename];
msg_file=[path_data,filename];

display(msg_file)

data=parse_CPFmsg_v4(msg_file); % for latest tank test

if isempty(data.hr_d)==1
    display('No hr data, skiping')
    continue
end

% See data.hr_hdr for header info
nz=length(data.hr_d);
P(1:nz,p)=data.hr_d(:,1);
T(1:nz,p)=data.hr_d(:,2);
S(1:nz,p)=data.hr_d(:,3);
N(1:nz,p)=data.hr_d(:,4);

if ~isempty(data.aircal)
O2_air(:,:,p)=data.aircal;
end

Lon(p)=data.gps(1);
Lat(p)=data.gps(2);
Date(p)=data.sdn;

% Calculate vertical velocity
[nz,np]=size(N);
w_ctd=nan(nz-1,np);
w_ctd2=w_ctd; wp_ctd=w_ctd; wp_ctd2=w_ctd; dT=w_ctd; dS=w_ctd; dz=dS;
p_axis=w_ctd;
% d_axis=w_ctd;
    
% Convert pressure to depth:
DEPTH(:,p)=-1*gsw_z_from_p(P(:,p),Lat(p)); % NO GPS Fix in the tank?

% Compute SA, CT and RHO:
[SA(:,p) in_ocean]= gsw_SA_from_SP(S(:,p),P(:,p),Lon(p),Lat(p)); % Absolute Salinity from Pactical salinity 
CT(:,p)=gsw_CT_from_t(SA(:,p),T(:,p),P(:,p)); % Conservative Temp
RHO(:,p)=gsw_rho(SA(:,p),CT(:,p),0); % Potential density, ref to the surface

% Bin size in db;
dp=diff(P(:,p));
p_axis(:,p)=P(2:end,p)-dp/2;

% Bin size in dz;
dz(:,p)=diff(DEPTH(:,p));
d_axis(1:nz-1,p)=DEPTH(2:end,p)-dz(:,p)/2;

% Compute Vertical velocity from the CTD
W_ctd1(1:nz-1,p)=dz(:,p)./N(1:end-1,p);
W_ctd(1:nz-1,p)=dz(:,p)./N(2:end,p);
% discrding the last N (from the last CTD pressure towards the surface,
% assuming data comes from ascending profiles). The last Nbin (closer) to the
% surface would have 1-2 samples per bin.. and, thus give unrealistacally
% large surfacing velocities?

% Compute Vertical velocity from the CTD
wp_ctd1(:,p)=dp./N(1:end-1,p);
wp_ctd(:,p)=dp./N(2:end,p);

clear data 
end
 
% NOTE: Floats' sample always at 1Hz *SBE 41cp CTD).
% They are either set to bin data in 1-m or 2-m bins
% Thus, larger N per bin, implies float was slowed down by the flow and
% thus were able to collect more samples per bin.

%% ------------ Plot ALL profile data (from a mission) --------------------

% Plot T,S and Nbin
figure(1)

subplot(131)
plot(T,DEPTH,'-o'); axis ij
xlim(range(T))
ylim([0 maxp])
xlabel('T (C)');

subplot(132)
plot(S,DEPTH,'-o'); axis ij
xlim(range(S))
ylim([0 maxp])
xlabel('S (psu)');

title(mission)

subplot(133)
plot(N,DEPTH,'-o'); axis ij
xlim(range(N))
ylim([0 maxp])
xlabel('N_{bin}');

print_setup(40,30)
% print([path_plots,'TS_Nbin_HR_allprofiles.eps'],'-depsc2');
print([path_plots,'TS_Nbin_HR_allprofiles.jpg'],'-djpeg');

%% Plot DENSITY, and W_ctd

figure(2)

subplot(131)
plot(N,DEPTH,'-o'); axis ij
xlim([0 60]) %xlim(range(N))
ylim([0 maxp])
xlabel('N_{bin}');
ylabel('Depth [m]');

subplot(132)
plot(RHO-1000,DEPTH,'-o'); axis ij
xlim(range(RHO-1000))
ylim([0 maxp])
xlabel('\sigma_o [kg/m^3]');

title(mission)

subplot(133)
plot(W_ctd,d_axis,'-o'); axis ij
xlim(range(W_ctd(1:end,:)))
ylim([0 maxp])
xlabel('w_{ctd} [m/s]');

print_setup(40,30)
% print([path_plots,'RHO_Wctd_HR_allprofiles.eps'],'-depsc2');
print([path_plots,'RHO_Wctd_HR_allprofiles.jpg'],'-djpeg');

%% Plot LON,LAT

figure(3)
scatter(Lon,Lat,100,Date,'filled');
% hc=colormap(m_colmap('jet','step',10));
hc=colormap(cbrewer('seq','YlOrRd',21));
h=colorbar;
% datetick(h,'keeplimits','keepticks')
cbdate('dd-mmm-yy HH:MM')
set(get(h,'Title'),'String','GMT','fontsize',20);


ylabel('Latitude [\circ]');
xlabel('Longitude [\circ]');

print_setup(60,40)
% maximize;
% print([path_plots,'RHO_Wctd_HR_allprofiles.eps'],'-depsc2');
print([path_plots,'Lon_vs_Lat_allprofiles.jpg'],'-djpeg');

disp(['Plots saved in ',path_plots])


%% Plot AirCal

figure(4)

subplot(311)
plot(squeeze(O2_air(:,1,:)),squeeze(O2_air(:,2,:)),'o'); 
ylabel('Topt [C]');
datetick('x','mmm-dd-yy HH:MM','keepticks','keeplimits')

title(['Air O_2 ', mission]);

subplot(312)
plot(squeeze(O2_air(:,1,:)),squeeze(O2_air(:,3,:)),'o'); axis ij
ylabel('Tphase');
datetick('x','mmm-dd-yy HH:MM','keepticks','keeplimits')

subplot(313)
plot(squeeze(O2_air(:,1,:)),squeeze(O2_air(:,4,:)),'o'); axis ij
ylabel('Rphase');
datetick('x','mmm-dd-yy HH:MM','keepticks','keeplimits')

print_setup(60,30)
print([path_plots,'O2_air_vs_time_allprofiles.jpg'],'-djpeg');

%% Plot DENSITY, and W_ctd

figure(2)

subplot(131)
plot(N,DEPTH,'-o'); axis ij
xlim(range(N))
ylim([0 maxp])
xlabel('N_{bin}');

subplot(132)
plot(RHO-1000,DEPTH,'-o'); axis ij
xlim(range(RHO-1000))
ylim([0 maxp])
xlabel('\sigma_o (kg/m^3)');

title(mission)

subplot(133)
plot(W_ctd,d_axis,'-o'); axis ij
xlim(range(W_ctd(1:end,:)))
ylim([0 maxp])
xlabel('w_{ctd} (m/s)');

print_setup(40,30)
% print([path_plots,'RHO_Wctd_HR_allprofiles.eps'],'-depsc2');
print([path_plots,'RHO_Wctd_HR_allprofiles.jpg'],'-djpeg');
