% function [P, A]=ReadMsgFile_mc(path_data,filename);

% Purpose: Read MSG file and extract LR profile data as well as BGC data
% from Anchor Points (i.e. when CPF is parked at depth). 

% clear all
% close all
% path_data='./';
% filename='CPF003.001r1.msg.txt'; % test file from Gene, no modifications needed from original .msg file
% 
% filename=msg_file;
% fid = fopen([path_data,filename]);
addpath(genpath('/Users/Mau/Documents/MATLABhelp_new'));

clear all
% msg_file ='/Volumes/CPF/CPF003/Data/ParagonTests/2022-3-14T16-30-13/CPF003.005.msg'; % MC, profiles that have Anchor data (~6h), 001,002,,003,005,013

% test anchor data from previous deployment
msg_file ='/Volumes/CPF/CPF003/Data/ParagonTests/2022-2-23T16-59-6/CPF003.001.msg'; % works
msg_file ='/Volumes/CPF/CPF003/Data/ParagonTests/2022-2-23T16-59-6/CPF003.002.msg'; % no file exists
msg_file ='/Volumes/CPF/CPF003/Data/ParagonTests/2022-2-23T16-59-6/CPF003.003.msg'; % works
msg_file ='/Volumes/CPF/CPF003/Data/ParagonTests/2022-2-23T16-59-6/CPF003.004.msg'; % number of columns is 19 or 18?
msg_file ='/Volumes/CPF/CPF003/Data/ParagonTests/2022-2-23T16-59-6/CPF003.006.msg'; % AirCal messing it up?

fid = fopen(msg_file);

% extract date stamp from msg_file name
ddate=regexp(msg_file,['(?<year>\d{4})-(?<month>\d+)-(?<day>\d+)'],'names');
dirs.date=[ddate.year,'-',ddate.month,'-',ddate.day]
cast=regexp(msg_file,'\d+(?=\.msg)','once','match');

rowIndex = 1;
getHeaderLine = true;
while ~feof(fid)
    inLine = fgetl(fid);
    if(getHeaderLine)
        if strncmp(inLine, "Comment", 7)
          dataLines{rowIndex} = vertcat(strsplit(inLine, '\t'))
          rowIndex = rowIndex + 1;
          getHeaderLine = false;
        end
    end
    
    if(strncmp(inLine, "AirCal", 6) || ...
        strncmp(inLine, "AnchorPt", 8) || ... 
        strncmp(inLine, "AscendPt", 8))
          dataLines{rowIndex} = vertcat(strsplit(inLine, '\t'))
          rowIndex = rowIndex + 1;
    end    
end

fclose(fid);

dataLines = dataLines';

% MC, modified the following lines to account for different columns for
% AscendPt/AnchorPt and AirCal data. Gene will eventually add NaNs to the
% two columns that have missing data (blanks)

for f=1:length(dataLines)
    tsize(f,:)=size(dataLines{f});
end

dataLines = vertcat(dataLines{find(tsize(:,2)==17)});
dataTable = array2table(dataLines(tsize(:,2)==17,:), 'VariableNames', dataLines(1,:));

dataTable.dt = datetime(dataTable.DateTime, 'InputFormat', 'yyyy-MMM-dd''T''HH:mm:ss''Z')

% plot(dataTable.dt, str2double(dataTable.pres_dBar))

%% Create index for AnchorPt lines
C=dataTable(:,1);
indA=zeros(size(C));

%% Will do this here, but eventually have to add to calc_BGC.m or create a function calc_PHYS.m
Pres=dataTable(:,3);
Temp=dataTable(:,4);
Sal=dataTable(:,5);
dtime=datenum(dataTable.dt);

nM=size(Pres); nM=nM(1);

for m=2:nM-1
    if strcmp(cell2mat(C{m,1}),'AnchorPt:')
    P(m-1,1)=str2num(cell2mat(Pres{m,1}));
    T(m-1,1)=str2num(cell2mat(Temp{m,1}));
    S(m-1,1)=str2num(cell2mat(Sal{m,1}));
    dt(m-1,1)=dtime(m,1);
    elseif strcmp(cell2mat(C{m,1}),'AscendPt:')
        break
    end
end

Lat=NaN;
Lon=NaN;

% Mean sampling rate:
dt_minutes=nanmean(diff(dt))*1440; % sampling rate in minutes

%%
p=1;

    if isnan(Lat(p))
        Lon(p)=-121.936; % aprox for the MB in case no GPS fix
        Lat(p)=36.7329;
    end

%     % Convert pressure to depth (only if there is a GPS fix):
%     if ~isnan(Lat(p))
%         DEPTH(:,p)=-1*gsw_z_from_p(P(:,p),Lat(p)); % NO GPS Fix in the tank?
%     else
%         DEPTH(:,p)=P(:,p); % aprox with pressure
%     end

    % 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

    %% Close all

    rangeRHO=range(RHO(1:end));
    deltaRHO=rangeRHO(2)-rangeRHO(1);

    rangeT=range(T(1:end));
    deltaT=rangeT(2)-rangeT(1);

    rangeP=range(P(1:end));
    deltaP=rangeP(2)-rangeP(1);

    %% Plot Timeseries of Anchor data
    
    close all
    subplot(311)
    plot(dt,RHO); 
    datetick
    ylabel(['\rho [kg/m^3]'])

    title([dirs.date,' - ', cast ,', Sampling rate = ',num2str(dt_minutes), ' minutes, T = ',num2str(nansum(diff(dt))*1440), ' minutes'])
   
    subplot(312)
    plot(dt,T);
    datetick
    ylabel(['T [^{\circ} C]'])

    subplot(313)
    plot(dt,S);
    datetick
    ylabel(['S [psu]'])

    xlabel(['Time (GMT)'])

    % legend({'T','S','\rho'})
    print_setup(60,30)

    path_plots='/Users/mau/Documents/MBARIprojects/CPF/CPFPROCESSING/PLOTS/ParagonTests/AnchorData/';
    print([path_plots,'TSRHO_vs_Time_',dirs.date,'_',cast,'.jpg'],'-djpeg')

%% autocorrelation analysis
    close all

    maxlags=60;
    [r, lags]= xcorr(detrend(RHO),'normalized');
    stem(lags*dt_minutes,r); hold on

    [r, lags]= xcorr(detrend(T),'normalized');
    stem(lags*dt_minutes,r)

    [r, lags]= xcorr(detrend(S),'normalized');
    stem(lags*dt_minutes,r)

    ylabel(['Autocorrelation, r'])
    xlabel(['Lag (minutes)'])
    
    xlim([0 360])
    legend({'\rho','T','S'})

    title([dirs.date,' - ', cast ,', Sampling rate = ',num2str(dt_minutes), ' minutes, T = ',num2str(nansum(diff(dt))*1440), ' minutes'])
    print_setup(60,30)

    path_plots='/Users/mau/Documents/MBARIprojects/CPF/CPFPROCESSING/PLOTS/ParagonTests/AnchorData/';
    print([path_plots,'DecorrelationTimeScale_TSRHO_',dirs.date,'_',cast,'.jpg'],'-djpeg')


    % Note ~zero crossing at ~ lag=10; for a sampling rate at the bottom of
    %% 
    alpha=6.9*10^(-5);
    gamma=2.4*10^(-6);

%    %% save data to append
%    anchor_hdr= dataLines{1,3:5};
%    anchor_d=[P T S];
% 
%    %truncate NaN
% 
%    anchor_d=anchor_d(~isnan(anchor_d),:)
