classdef GliderDataHandler
    properties (Constant)
        dataFolder = '\\atlas.shore.mbari.org\ProjectLibrary\901805_Coastal_Biogeochemical_Sensing\Spray_Data\';
        outputFile = '\\atlas.shore.mbari.org\ProjectLibrary\901805_Coastal_Biogeochemical_Sensing\Locness\LocnessMapProduct.txt';
        vars = { ...
                    'Cruise', 'Platform', 'Layer', 'CastDirection', ...
                    'unixTimestamp', 'lat', 'lon', 'temperature',...
                    'salinity', 'pHin', 'pH25atm', 'rhodamine', 'MLD', ...
                };
        density_threshold = .03; % Adjust for MLD processing
        reference_pressure = 15; % Adjust for MLD processing
    end

    methods
        function processGliderData(obj, depID)
            % Construct full path to .mat file
            mFile = fullfile(obj.dataFolder, depID, [depID 'sat.mat']);

            if ~isfile(mFile)
                error("Data file not found: %s", mFile);
            end

            % Load data
            s = load(mFile);
            s = s.s;
            % Process to table
            T = obj.buildTable(s);

            % Read existing data and append new rows
            
            keyVars = ["unixTimestamp", "Platform", "Layer"];
            if isfile(obj.outputFile)
                T_old = readtable(obj.outputFile);
                T_old.Properties.VariableNames = obj.vars;
                iappend = ~ismember(T(:,keyVars), T_old(:,keyVars), 'rows');
                if any(iappend)
                    writetable(T(iappend,:), obj.outputFile, "WriteMode", "append");
                end
            else
                writetable(T, obj.outputFile);
            end
        end

        function T = buildTable(obj, s)
            nRows = length(s.sdn);
            % MLD calculation
            [s] = calcMLD(s, obj.density_threshold, obj.reference_pressure);
            % Check the timestamp and the gps point!!!!
            % Build table
            % Mean of the upcast and downcast with the gps location at the
            % diveend
%             if all(s.direction) == 1 || all(s.direction) ~=-1
            T1 = table( ...
                repmat(string(s.depID), nRows, 1), ...
                repmat("Glider", nRows, 1), ...
                repmat("MLD", nRows, 1), ...
                repmat("Mean", nRows, 1), ...
                posixtime(datetime(datevec(s.sdn(:)))), ...
                s.lat(:), ...
                s.lon(:), ...
                s.tc_mld', ...
                s.psal_mld', ...
                s.pHin_mld', ...
                s.pH25atm_mld', ...
                s.rhodamine_mld', ... 
                s.mld', ...
                'VariableNames', obj.vars);
            
            % Ascent data comes from the dive end location. Mixed layer
            % depth on the upcast is closest to dive end
            idx = s.direction == 1;
            T2 = table( ...
                repmat(string(s.depID), nRows, 1), ...
                repmat("Glider", nRows, 1), ...
                repmat("MLD", nRows, 1), ...
                repmat("Up", nRows, 1), ...
                posixtime(datetime(datevec(s.sdn_(2,:)))), ...
                s.lat_(2,:)', ... 
                s.lon_(2,:)', ...
                s.tc_ts(idx)', ...
                s.psal_ts(idx)', ...
                s.phin_ts(idx)', ...
                s.pH25atm_ts(idx)', ...
                s.rhodamine(idx)', ... 
                s.mld_ts(idx)', ...
                'VariableNames', obj.vars);

            % Descent data comes from the dive start location. Mixed layer
            % depth on the downcast is closest to dive start
            idx = s.direction == -1;
            T3 = table( ...
                repmat(string(s.depID), nRows, 1), ...
                repmat("Glider", nRows, 1), ...
                repmat("MLD", nRows, 1), ...
                repmat("Down", nRows, 1), ...
                posixtime(datetime(datevec(s.sdn_(1,:)))), ...
                s.lat_(1,:)', ...
                s.lon_(1,:)', ...
                s.tc_ts(idx)', ...
                s.psal_ts(idx)', ...
                s.phin_ts(idx)', ...
                s.pH25atm_ts(idx)', ...
                s.rhodamine(idx)', ... 
                s.mld_ts(idx)', ...
                'VariableNames', obj.vars);
            
            % Surface values
            depth = s.depth;     % n x m
            
            % Mask for valid depth range
            validMask = depth > 0 & depth < 10;
            % Set all other pHin values to NaN
            maskedPH = s.pHin; maskedPH(~validMask) = NaN;
            maskedpH25atm = s.pH25atm; maskedpH25atm(~validMask) = NaN;
            maskedtemperature = s.tc; maskedtemperature(~validMask) = NaN;
            maskedsalinity = s.psal; maskedsalinity(~validMask) = NaN;
            maskedrhodamine = s.rhodamine; maskedrhodamine(~validMask) = NaN;
            % Compute column-wise mean, omitting NaNs
            pHin_surface = mean(maskedPH, 'omitnan')';
            pH25atm_surface = mean(maskedpH25atm, 'omitnan')';
            temperature_surface = mean(maskedtemperature, 'omitnan')';
            salinity_surface = mean(maskedsalinity, 'omitnan')';
            rhodamine_surface = mean(maskedrhodamine, 'omitnan')';
            T4 = table( ...
                repmat(string(s.depID), nRows, 1), ...
                repmat("Glider", nRows, 1), ...
                repmat("Surface", nRows, 1), ...
                repmat("Mean", nRows, 1), ...
                posixtime(datetime(datevec(s.sdn_(2,:)))), ...
                s.lat_(2,:)', ...
                s.lon_(2,:)', ...
                temperature_surface,...
                salinity_surface,...
                pHin_surface, ...
                pH25atm_surface, ...
                rhodamine_surface, ...
                NaN(nRows,1), ...
                'VariableNames',obj.vars);
            
            T = vertcat(T1,T2,T3,T4); % Need to fix the upcast situation!!!!!!
%             if ~ismember(s.depID,'25706901')
%                 T = vertcat(T1,T2,T3,T4); % surface
%             else
%                 T = vertcat(T1,T2,T3,T4); % Just surface and upcast (might need to flip upcast)
%             end
%             T = vertcat(T1,T2,T3); % meancast, upcast, downcast
        end
    end
end
