%%Note: Not final format. I am extrapolating over pressure parsing errors instead of NaNing out for now.
t = load("phTest_2504020901.mat");
t.ctd.t = cellfun(@double, t.ctd.t, 'UniformOutput', false); % Convert unit 16 to double
t.ph.t = cellfun(@double, t.ph.t, 'UniformOutput', false); % Convert unit 16 to double

% now fill in the matrix variables that exist. each parameter is on a
% different pres grid; interpolate them onto the 'pH' pressure axis

pHvars = {'pres'    'P'
    'vrs'           'Vrse'
    'vrs_std'       'Vrse_std'
    'vk'            'Vk'
    'vk_std'        'Vk_std'
    'Ik'            'Ik'
    'Ib'            'Ib'};

% Extract the downcast and create an ODV style table then repeat for the
% upcast
Cast_Direction = 1; % 1: down, 2: up
c = max(length(t.ctd.t)); % number of profiles
% Checks that time and the variable to interpolate are both non empty
% values before attempting to index into them or interpolate.

for i = 1:c
    % Extract variables in 'pH' from JSON
    l = length(t.ph.P{i,Cast_Direction});
    for j = 1:size(pHvars,1)
        s.([pHvars{j,1}])(1:l,i) = t.ph.(pHvars{j,2}){i,Cast_Direction};
    end

    % CTD temp
    if ~isempty(t.ctd.t{i}) && ~isempty(t.ctd.T{i})
        [xt, idx] = unique(t.ctd.t{i,1}, 'stable');
        T = t.ctd.T{i,1};
        yt = T(idx);
        s.tc(1:l,i) = interp1(xt, yt, t.ph.t{i,Cast_Direction}, 'linear', 'extrap');
    else
        s.tc(:,i) = NaN;
    end

    % CTD salinity
    if ~isempty(t.ctd.t{i}) && ~isempty(t.ctd.S{i})
        [xt, idx] = unique(t.ctd.t{i,1}, 'stable');
        yt = t.ctd.S{i,1}(idx);
        s.psal(1:l,i) = interp1(xt, yt, t.ph.t{i,Cast_Direction}, 'linear', 'extrap');
    else
        s.psal(:,i) = NaN;
    end
end

sz = size(s.vrs);
dnum = 1:sz(2);
profs = repmat(dnum,sz(1),1);
s.Cast_Direction = 'D';

ODV_down = table;
ODV_down.Station = reshape(profs,[],1);
ODV_down.Cast_Direction = repmat(s.Cast_Direction,height(ODV_down),1);
ODV_down.pres = reshape(s.pres,[],1);
ODV_down.tc = reshape(s.tc,[],1);
ODV_down.psal = reshape(s.psal,[],1);
ODV_down.vrs = reshape(s.vrs,[],1);
ODV_down.vrs_std = reshape(s.vrs_std,[],1);
ODV_down.vk = reshape(s.vk,[],1);
ODV_down.vk_std = reshape(s.vk_std,[],1);
ODV_down.Ik = reshape(s.Ik,[],1);
ODV_down.Ib = reshape(s.Ib,[],1);
%%
clear s;

Cast_Direction = 2; % 1: down, 2: up
c = max(length(t.ctd.t)); % number of profiles
% Checks that time and the variable to interpolate are both non empty
% values before attempting to index into them or interpolate.

for i = 1:c
    % Extract variables in 'pH' from JSON
    l = length(t.ph.P{i,Cast_Direction});
    for j = 1:size(pHvars,1)
        s.([pHvars{j,1}])(1:l,i) = t.ph.(pHvars{j,2}){i,Cast_Direction};
    end

    % CTD temp
    if ~isempty(t.ctd.t{i}) && ~isempty(t.ctd.T{i})
        [xt, idx] = unique(t.ctd.t{i,1}, 'stable');
        T = t.ctd.T{i,1};
        yt = T(idx);
        s.tc(1:l,i) = interp1(xt, yt, t.ph.t{i,Cast_Direction}, 'linear', 'extrap');
    else
        s.tc(:,i) = NaN;
    end

    % CTD salinity
    if ~isempty(t.ctd.t{i}) && ~isempty(t.ctd.S{i})
        [xt, idx] = unique(t.ctd.t{i,1}, 'stable');
        yt = t.ctd.S{i,1}(idx);
        s.psal(1:l,i) = interp1(xt, yt, t.ph.t{i,Cast_Direction}, 'linear', 'extrap');
    else
        s.psal(:,i) = NaN;
    end
end

sz = size(s.vrs);
dnum = 1:sz(2);
profs = repmat(dnum,sz(1),1);
s.Cast_Direction = 'U';

ODV_up = table;
ODV_up.Station = reshape(profs,[],1);
ODV_up.Cast_Direction = repmat(s.Cast_Direction,height(ODV_up),1);
ODV_up.pres = reshape(s.pres,[],1);
ODV_up.tc = reshape(s.tc,[],1);
ODV_up.psal = reshape(s.psal,[],1);
ODV_up.vrs = reshape(s.vrs,[],1);
ODV_up.vrs_std = reshape(s.vrs_std,[],1);
ODV_up.vk = reshape(s.vk,[],1);
ODV_up.vk_std = reshape(s.vk_std,[],1);
ODV_up.Ik = reshape(s.Ik,[],1);
ODV_up.Ib = reshape(s.Ib,[],1);
%%
% Combine
ODV = [ODV_down; ODV_up];
% Sort by Station and Cast_Direction ('U' before 'D')
[~, idx] = sortrows(ODV, {'Station', 'Cast_Direction'});
ODV = ODV(idx, :);

% Trim the file by removing empty rows. This is an artifact of having data
% stored in a matrix
idx_trim = ODV.pres == 0;
ODV(idx_trim,:) = [];
%% Calc pHin
% Pull cal data from excel sheet and join tables with most recent info
MissionID = '25420901'; % Set MissionID. This links SN and GliderFETID in the CalibrationLog
cd('\\atlas.shore.mbari.org\ProjectLibrary\901805_Coastal_Biogeochemical_Sensing\nanoFET\'); % Go into folder of cal file
fname_CalibrationLog = 'CalibrationLog.xlsx';
% Read each worksheet as a table
Deployments = readtable(fname_CalibrationLog, 'Sheet', 'Deployments');
Deployments = Deployments(contains(Deployments.MissionID,MissionID),:); % Get only data from active sensor
SN = char(Deployments.SN);
SNID = ['0', SN];
GliderFETID = Deployments.GliderFETID;
PTC = readtable(fname_CalibrationLog, 'Sheet', 'PTC');
PTC = PTC(ismember(PTC.GliderFETID,GliderFETID),:); % Get only active sensor data
k0 = readtable(fname_CalibrationLog, 'Sheet', 'k0');
k0 = k0(ismember(k0.GliderFETID,GliderFETID),:); % Get only data from active sensor
k0 = sortrows(k0, 'CalibrationDate', 'descend'); % Sort by CalibrationDate (most recent first)
[~, firstIdx] = unique(k0.GliderFETID, 'stable'); % Keep only the most recent entry
k0 = k0(firstIdx, :);
% Join tables on the common key column 'GliderFETID'
Cal = join(PTC, k0, 'Keys', 'GliderFETID'); % Join PTC and k0 on 'GliderFETID'
Cal = join(Deployments, Cal, 'Keys','GliderFETID'); % Join Deployments and Cal on 'GliderFETID'
% Calc pHin
k0 = Cal.k0_seawater;
k2 = Cal.k2F_P_C0;
fp = [Cal.f_P_K1, Cal.f_P_K2, Cal.f_P_K3, Cal.f_P_K4, Cal.f_P_K5, Cal.f_P_K6]';
fp(isnan(fp)) = 0; % NaNs do not work in phcalc_jp. Need to make 0's.
ODV.pHin = phcalc_jp(ODV.vrs, ODV.pres, ODV.tc, ODV.psal, k0, k2, fp); % pHin accounting for pressure affects on vrs
%% Write ODV file
folder = 'C:\Users\bwerb\Documents\dep_25420901';
savefilename = fullfile(folder,"2504020901ODV.txt");
writetable(ODV,'2504020901ODV.txt');