clear all
close all

addpath(genpath('/Users/Mau/Documents/MATLABhelp_new/'));
% PURPOSE: MC clenaing up Gene's code that reads data .eng files
% original version (processCPFEngrLog3.m) of this code copied from
% //atlas.shore.mbari.org/ProjectLibrary/CPF/Data/AnalysisScripts/processCPFEngrLog3.m
% on Sep 13th, 2021

% Code Modifications:
% 09/13/2021, MC:
%    - modified fileList and FileName (using dir instead of ls to list files).
%    - included OCR data header in varNames (i.e. varNames and tableData
%      columns didn't match)

%     if(fileNum == 5)
plotEachProfile = true;
%     else
%         plotEachProfile = false;
%     end

plotWaterfalls = true;
plotWaterfalls = false;

% % startingDir = '//atlas/ProjectLibrary/CPF/CPF003/Data/TankTests/2021-8-25T21-57-38/';
% startingDir = '/Volumes/CPF/CPF003/Data/TankTests/2021-9-9T21-6-13/'; %
% Testing, code runs, thorughs a billion figures!
% fileList = ls([startingDir,'*.eng']); MC, modifying below to use dir and
% define mission/folder

%% MC, Path to data (CPF mission)

mission='/ParagonTests/MICRORIDER/2021-11-16T19-16-49'; %OVERNIGHT mission, 19 profiles
mission='/ParagonTests/MICRORIDER/2021-10-19T16-3-9/'; % Half-day mission (6h, 8 profiles)
% mission='/ParagonTests/MICRORIDER/2021-10-15T16-10-6/';
% % NOTE mission needs '/' at the beginning, so that code works both in macOS
% % and windows

% mission='TankTests/MICRORIDERTESTS/2021-10-5T22-38-16/'; % TankTest CPF+MR, SH1 & T1 probes
% mission='ParagonTests/2021Sep27/2021-9-27AllProfiles/';
% mission='ParagonTests/2021Sep27/2021-9-27T18-12-43/';
% mission='BenchTests/2021Sep22/2021-9-22T21-58-0/';
% mission='BenchTests/2021Sep22/2021-9-22T21-30-59/';
% mission='ParagonTests/2021-9-16T16-33-30/';
% mission='Paragon Tests/2021-9-16T16-33-30/';
% mission='TankTests/2021-9-9T21-6-13/';
% mission='Paragon Tests/2021Apr02/2021-4-2T17-48-47/';

mission='/ParagonTests/2022-2-23T16-59-6/';


maxp=100; % MC,for plotting purposes

data_root='/Volumes/CPF/CPF003/Data';
plots_root=['/Volumes/CPF/CPF003/Data/data_processing/plots'];
path_data=[data_root,mission];

data_local='/Users/Mau/Documents/MBARIprojects/CPF/CPF003/Data/';
% path_mat=[data_local,'matfiles',mission]; % NOTE mission needs / at the beginning
% path_mat=[data_local,'matfiles',mission(1:35),'/']; % works if inside /MICRORIDER directory
path_mat=[data_local,'matfiles',mission,'/']; % not truncating mission name anymore (to avoid futture confusion!)

matfiles=dir([path_mat,'*.mat']); % to append ENG structure to .mat file
%NOTE need to run process_CPFprofiles.m first (which saves '*.mat' files)

startingDir = path_data;

fileList = dir([startingDir,'/*.eng']); % MC, fileList = ls([startingDir,'*.eng']);
numFiles = size(fileList, 1);

%%  MC, Create directory to save plots (if it doesn't exist already)
% path_plots=[plots_root,mission];
path_plots=['/Users/Mau/Documents/MBARIprojects/w_ctd/cpf/plots',mission]; %(1:35)];

if isfolder(path_plots)==0 && numel(dir(path_plots))<=2 % will create directory only if it is empty!
    mkdir(path_plots)
end

% % path to engineering plots (create subdirectory in the above path_plots)
% path_plots_eng=[plots_root,mission,'/eng/'];
% if isfolder(path_plots_eng)==0 && numel(dir(path_plots_eng))<=2 % will create directory only if it is empty!
%     cd(path_plots)
%     mkdir('./eng')
% end

%% Variable names
fileFormat = 3;  %1 = pre 2021Feb26, 2 = 1 to 2021Mar15, 3 = post 2021Mar15

if(fileFormat == 3)
    varNames = {'stx', 'date', 'time', 	'state', 	'comment', 	...
        'pressure', 'bellowsPos', ...
        'pumpVolts', 	'pumpMilliamps', 'pumpCoulombs', 'pumpJoules', 'pumpMotor'...
        'hotelVolts', 	'hotelMilliamps', 'hotelCoulombs', 'hotelJoules', ...
        'CTDTemp', 	'CTDSal', 	...
        'RawPress', 	'LPFPress', 'LPFVel', 	'diffVel', 	'diff2Vel', 'slopeOfLPFPress', 	'CTDSamplePeriod', ...
        'LPS331Press', 	'LPS311Temp', 	'SHT21Humid', 	'SHT21Temp', ...
        'pidSP', 	'pidPV', 	'pidY1',	'pidY2', 'pidI', 	'pid_v', 'pid_u', 'pid_uCPS', 'pid_h',  ...
        'FLBB', 'Topt', 'Tphase', 'Rphase', 'NO3', 'pH_Vrs', 'OCR', 'GPS', 'IridModem', 'MSC'};
    varTypes = {'string', 'string', 'string', 'string', 'string', ...
        'double', 'double', ...
        'double', 'double','double', 'double', 'string'...
        'double', 'double','double',  'double', ...
        'string', 'string', ...
        'double', 'double','double', 'double', 'double','double', 'double', ...
        'double', 'double','double', 'double', ...
        'double', 'double','double', 'double', 'double', 'double', 'double','double', 'double', ...
        'string', 'string', 'string', 'string', 'string', 'string', 'string', 'string', 'string', 'string'};
elseif(fileFormat == 2)
    varNames = {'stx', 'date', 'time', 	'state', 	'comment', 	...
        'pressure', 'bellowsPos', ...
        'pumpVolts', 	'pumpMilliAmps', 	'pumpJoules', ...
        'hotelVolts', 	'hotelMilliAmps', 	'hotelJoules', ...
        'CTDTemp', 	'CTDSal', 	...
        'RawPress', 	'LPFPress', 'LPFVel', 	'diffVel', 	'diff2Vel', 'slopeOfLPFPress', 	'CTDSamplePeriod', ...
        'LPS331Press', 	'LPS311Temp', 	'SHT21Humid', 	'SHT21Temp', ...
        'pidSP', 	'pidPV', 	'pidY1',	'pidY2', 'pidI', 	'pid_v', 'pid_u', 'pid_uCPS', 'pid_h',  ...
        'FLBB', 'Topt', 'Tphase', 'Rphase', 'NO3', 'pH_Vrs', 'OCR', 'GPS', 'IridModem', 'MSC'};
    varTypes = {'string', 'string', 'string', 'string', 'string', ...
        'double', 'double', ...
        'double', 'double','double',  ...
        'double', 'double','double',  ...
        'string', 'string', ...
        'double', 'double','double', 'double', 'double','double', 'double', ...
        'double', 'double','double', 'double', ...
        'double', 'double','double', 'double', 'double', 'double', 'double','double', 'double', ...
        'string', 'string', 'string', 'string', 'string', 'string', 'string', 'string', 'string', 'string'};
else
    varNames = {'stx', 'date', 'time', 	'state', 	'comment', 	'pressure', 'bellowsPos', ...
        'pumpVolts', 	'pumpMAmps', 	'pumpJoules', ...
        'hotelVolts', 	'hotelMAmps', 	'hotelJoules', ...
        'CTDTemp', 	'CTDSal', 	'RawPress', 	'LPFPress', 'LPFVel', 	'diffVel', 	'diff2Vel', 'slopeOfLPFPress', 	'CTDSamplePeriod', ...
        'LPS331Press', 	'LPS311Temp', 	'SHT21Humid', 	'SHT21Temp', ...
        'PID_SP', 	'PID_PV', 	'PID_Y1',	'PID_Y2', 'PID_I', 	'PID_v', 'PID_u', 'PID_uCPS', 'PID_h',  ...
        'FLBB', 'Topt', 'Tphase', 'Rphase', 'NO3', 'pH_Vrs', 'OCR', 'GPS', 'IridModem', 'MSC'};
    varTypes = {'string', 'string', 'string', 'string', 'string', ...
        'double', 'double', ...
        'double', 'double','double',  ...
        'double', 'double','double',  ...
        'string', 'string', ...
        'double', 'double','double', 'double', 'double','double', 'double', ...
        'double', 'double','double', 'double', ...
        'double', 'double','double', 'double', 'double', 'double', 'double','double', 'double', ...
        'string', 'string', 'string', 'string', 'string', 'string', 'string', 'string', 'string', 'string'};

    dateIndex = 2;
    timeIndex = 3;
    stateIndex = 4;
    commentIndex = 5;
    pressureIndex = 6;
    bellowsPosIndex = 7;
    pumpVoltsIndex = 8;
    pumpAmpsIndex = 9;
    pumpJoulesIndex = 10;
    hotelVoltsIndex = 11;
    hotelAmpsIndex = 12;
    hotelJoulesIndex = 13;
    ctdIndex = 15;
    rawPressureIndex = 16;
    InternalPressureIndex = 23;
    IPTemperatureIndex = 20;
    HumidityIndex = 21;
    HumidityTemperatureIndex = 22;
    pumpSpeedIndex = 30;
    gpsIndex = 39;
    iridiumIndex = 40;
end

%% Process files
figX = 0.05;
figXIncr = 0.01;
figY = 0.3;
figYIncr = -0.02;

for fileNum = 1:numFiles-1
    %% Open and read file
    %     if(fileNum == 5)
    %       plotEachProfile = true;
    %     else
    % plotEachProfile = false;
    %     end

    % fullFileName = [startingDir fileList(fileNum,:)];
    fullFileName = [startingDir,'/',fileList(fileNum).name]; % MC
    fprintf('%s\n', [datestr(now, 'HH:MM:SS ') 'Importing Filename: ' fullFileName ', file number: ' int2str(fileNum) ' of ' int2str(length(fileList))]);

    opts = delimitedTextImportOptions('VariableNames', varNames, 'VariableTypes', varTypes, 'Delimiter', '\t');
    tableData = readtable(fullFileName, opts);

    profile(fileNum).fileName = fullFileName;

    %% Process date/time
    fprintf('%s\n', [datestr(now, 'HH:MM:SS ') 'Processing date/time ...']);

    dt = tableData.date + tableData.time;
    DT = datetime(dt, 'InputFormat', 'dd-MMM-yyyyHH:mm:ss.SSS');
    DT.Format='yyyyMMMdd HH:mm:ss';

    %% Process pressure
    fprintf('%s\n', [datestr(now, 'HH:MM:SS ') 'Processing pressure ...']);
    notNaN = ~isnan(tableData.RawPress);
    profile(fileNum).pressure = tableData.RawPress(notNaN);
    profile(fileNum).pressureDT = DT(notNaN);

    %% Process bellows position
    fprintf('%s\n', [datestr(now, 'HH:MM:SS ') 'Processing bellows position ...']);
    notNaN = ~isnan(tableData.bellowsPos);
    profile(fileNum).bellowsPosition = tableData.bellowsPos(notNaN);
    profile(fileNum).bellowsPositionDT = DT(notNaN);

    %% Find glide on
    glideOnFind = strfind(tableData.comment, 'glide on');
    glideOnMsgIndexes = find(isempty(glideOnFind) == 0);


    %% Find glide off
    glideOffFind = strfind(tableData.comment, 'glide off');
    glideOffMsgIndexes = find(isempty(glideOffFind) == 0);


    %% Plot pressure and bellows position
    %     if plotEachProfile
    %         fig = figure('Name',fullFileName,'NumberTitle','off');
    %         fig.Units = 'normalized';
    %         fig.Position = [figX figY 0.6 0.6];
    %         figX = figX + figXIncr;
    %         figY = figY + figYIncr;
    %
    %         hAx = subplot(3,1,1);
    %         yyaxis left
    %         plot(profile(fileNum).bellowsPositionDT, profile(fileNum).bellowsPosition);
    %         ylabel('Bellows Position (%IB Full)');
    %         yyaxis right
    %         plot(profile(fileNum).pressureDT, profile(fileNum).pressure);
    %         ylabel('Pressure (dBar)');
    %         hold on;
    %
    %         plotStateTransitions(tableData, DT);
    %
    %         %         if length(glideOnMsgIndexes) > 0
    %         %             plot(glideOntime, glideOnDepth, 'gO')
    %         %             text(glideOntime, glideOnDepth, 'glide on');
    %         %         end
    %         %         glideOnMsgIndexes = [];
    %         %         if length(glideOffMsgIndexes) > 0
    %         %             plot(glideOfftime, glideOffDepth, 'rX')
    %         %             text(glideOfftime, glideOffDepth, 'glide off');
    %         %         end
    %         %         glideOffMsgIndexes = [];
    %         hold off;
    % end

    %% Process pump speed
    fprintf('%s\n', [datestr(now, 'HH:MM:SS ') 'Processing pump speed ...']);

    if(~ismissing(tableData.pid_uCPS(1)))
        notNaN = ~isnan(tableData.pid_uCPS);
        profile(fileNum).pumpSpeed = tableData.pid_uCPS(notNaN);
        profile(fileNum).pumpSpeedDT = DT(notNaN);

        %         if plotEachProfile
        %             hAx = subplot(3,1,2);
        %             yyaxis left
        %             plot(profile(fileNum).pumpSpeedDT, profile(fileNum).pumpSpeed);
        %             ylabel('Pump Speed (cps)');
        %             yyaxis right
        %             plot(profile(fileNum).pressureDT, profile(fileNum).pressure);
        %             ylabel('Pressure (dBar)');
        %             hold on;
        %
        %             plotStateTransitions(tableData, DT);
        %         end
    end


    %% Process pressure and humidity
    fprintf('%s\n', [datestr(now, 'HH:MM:SS ') 'Processing internal pressure and humidity ...']);
    notNaN = ~isnan(tableData.LPS331Press);
    profile(fileNum).intPressure = tableData.LPS331Press(notNaN);
    profile(fileNum).intPressureDT = DT(notNaN);

    notNaN = ~isnan(tableData.SHT21Humid);
    profile(fileNum).humidity = tableData.SHT21Humid(notNaN);
    profile(fileNum).humidityDT = DT(notNaN);

    %     if plotEachProfile
    %         hAx = subplot(3,1,3);
    %         yyaxis left
    %         plot(profile(fileNum).intPressureDT, profile(fileNum).intPressure);
    %         ylabel('Internal Pressure (mBar)');
    %         yyaxis right
    %         plot(profile(fileNum).humidityDT, profile(fileNum).humidity);
    %         ylabel('Humidity (%RH)');
    %         hold on;
    %
    %         plotStateTransitions(tableData, DT);
    %     end
    %
    %     % MC, print figure
    %     print_setup(80,30)
    %     print([path_plots_eng,'fig0_',fileList(fileNum).name,'.jpg'],'-djpeg');

    %% Process pump power
    fprintf('%s\n', [datestr(now, 'HH:MM:SS ') 'Calculating Pump Power ...']);

    notNaN = ~isnan(tableData.pumpVolts);
    profile(fileNum).pumpVolts = tableData.pumpVolts(notNaN);
    profile(fileNum).pumpVoltsDT = DT(notNaN);

    notNaN = ~isnan(tableData.pumpMilliamps);
    profile(fileNum).pumpAmps = tableData.pumpMilliamps(notNaN) / 1000.0;
    profile(fileNum).pumpAmpsDT = DT(notNaN);

    notNaN = ~isnan(tableData.pumpCoulombs);
    profile(fileNum).pumpCoulombs = tableData.pumpCoulombs(notNaN);
    profile(fileNum).pumpCoulombsDT = DT(notNaN);

    notNaN = ~isnan(tableData.pumpJoules);
    profile(fileNum).pumpJoules = tableData.pumpJoules(notNaN);
    profile(fileNum).pumpJoulesDT = DT(notNaN);

    vecLength = min(length(profile(fileNum).pumpAmpsDT), length(profile(fileNum).pumpVoltsDT));
    profile(fileNum).pumpPowerCalc = profile(fileNum).pumpVolts(1:vecLength) .* profile(fileNum).pumpAmps(1:vecLength);
    profile(fileNum).pumpPowerCalcDT = profile(fileNum).pumpVoltsDT(1:vecLength);


    %     if plotEachProfile
    %         fig = figure('Name',fullFileName,'NumberTitle','off');
    %         fig.Units = 'normalized';
    %         fig.Position = [figX figY 0.6 0.6];
    %         figX = figX + figXIncr;
    %         figY = figY + figYIncr;
    %
    %         hAx = subplot(4,1,1);  % Auto-fitted to the figure.
    %         yyaxis left
    %         plot(profile(fileNum).pumpVoltsDT, profile(fileNum).pumpVolts);
    %         ylabel('Pump Volts');
    %         yyaxis right
    %         plot(profile(fileNum).pumpAmpsDT, profile(fileNum).pumpAmps);
    %         ylabel('Pump Amps');
    %         hold on;
    %
    %         plotStateTransitions(tableData, DT)
    %     end

    fprintf('%s\t', [datestr(now, 'HH:MM:SS ') 'Calculating Energy ...']);
    profile(fileNum).pumpEnergyCalc = profile(fileNum).pumpPowerCalc(1:end-1) .* seconds( profile(fileNum).pumpPowerCalcDT(2:end) - profile(fileNum).pumpPowerCalcDT(1:end-1));
    profile(fileNum).pumpEnergyCalcDT = profile(fileNum).pumpAmpsDT(1:end-1);

    profile(fileNum).pumpEnergyCalcSum = cumsum(profile(fileNum).pumpEnergyCalc);
    profile(fileNum).pumpEnergyCalcSumDT = profile(fileNum).pumpEnergyCalcDT;

    profile(fileNum).totalPumpEnergy = sum(profile(fileNum).pumpEnergyCalc);
    fprintf('Total Energy this profile = %.0f Joules\n\n', profile(fileNum).totalPumpEnergy);

    %     if plotEachProfile
    %         hAx = subplot(4,1,2);  % Auto-fitted to the figure.
    %         yyaxis left
    %         plot(profile(fileNum).pumpPowerCalcDT, profile(fileNum).pumpPowerCalc);
    %         ylabel('Pump Power (Watts)');
    %         yyaxis right
    %         plot(profile(fileNum).pressureDT, profile(fileNum).pressure);
    %         ylabel('Pressure (dBar)');
    %         hold on;
    %
    %         xlimits = xlim;
    %         text(xlimits(1), 0.95*max(profile(fileNum).pumpPowerCalc), sprintf('Energy = %.0f Joules', profile(fileNum).totalPumpEnergy));
    %
    %         plotStateTransitions(tableData, DT)
    %     end

    %% Process hotel power
    fprintf('%s\n', [datestr(now, 'HH:MM:SS ') 'Calculating Hotel Power ...']);

    notNaN = ~isnan(tableData.hotelVolts);
    profile(fileNum).hotelVolts = tableData.hotelVolts(notNaN);
    profile(fileNum).hotelVoltsDT = DT(notNaN);

    notNaN = ~isnan(tableData.hotelMilliamps);
    profile(fileNum).hotelAmps = tableData.hotelMilliamps(notNaN) / 1000.0;
    profile(fileNum).hotelAmpsDT = DT(notNaN);

    notNaN = ~isnan(tableData.hotelCoulombs);
    profile(fileNum).hotelCoulombs = tableData.hotelCoulombs(notNaN);
    profile(fileNum).hotelCoulombsDT = DT(notNaN);

    notNaN = ~isnan(tableData.hotelJoules);
    profile(fileNum).hotelJoules = tableData.hotelJoules(notNaN);
    profile(fileNum).hotelJoulesDT = DT(notNaN);

    vecLength = min(length(profile(fileNum).hotelAmpsDT), length(profile(fileNum).hotelVoltsDT));
    profile(fileNum).hotelPowerCalc = profile(fileNum).hotelVolts(1:vecLength) .* profile(fileNum).hotelAmps(1:vecLength);
    profile(fileNum).hotelPowerCalcDT = profile(fileNum).hotelVoltsDT(1:vecLength);

    %     if plotEachProfile
    %         hAx = subplot(4,1,3);  % Auto-fitted to the figure.
    %         yyaxis left
    %         plot(profile(fileNum).hotelVoltsDT, profile(fileNum).hotelVolts);
    %         ylabel('Hotel Volts');
    %         yyaxis right
    %         plot(profile(fileNum).hotelAmpsDT, profile(fileNum).hotelAmps);
    %         ylabel('Hotel milliAmps');
    %         hold on;
    %
    %         plotStateTransitions(tableData, DT)
    %     end
    %
    %     fprintf('%s\t', [datestr(now, 'HH:MM:SS ') 'Calculating Energy ...']);
    %     profile(fileNum).hotelEnergyCalc = profile(fileNum).hotelPowerCalc(1:end-1) .* seconds( profile(fileNum).hotelPowerCalcDT(2:end) - profile(fileNum).hotelPowerCalcDT(1:end-1));
    %     profile(fileNum).hotelEnergyCalcDT = profile(fileNum).hotelPowerCalcDT(2:end);
    %
    %     profile(fileNum).hotelEnergyCalcSum = cumsum(profile(fileNum).hotelEnergyCalc);
    %     profile(fileNum).hotelEnergyCalcSumDT = profile(fileNum).hotelEnergyCalcDT;
    %
    %     profile(fileNum).totalHotelEnergy = sum(profile(fileNum).hotelEnergyCalc);
    %     fprintf('Total Energy this profile = %.0f Joules\n\n', profile(fileNum).totalHotelEnergy);
    %
    %     if plotEachProfile
    %         hAx = subplot(4,1,4);  % Auto-fitted to the figure.
    %         yyaxis left
    %         plot(profile(fileNum).hotelPowerCalcDT, profile(fileNum).hotelPowerCalc);
    %         ylabel('Hotel Power (Watts)');
    %         yyaxis right
    %         plot(profile(fileNum).pressureDT, profile(fileNum).pressure);
    %         ylabel('Pressure (dBar)');
    %         hold on;
    %
    %         xlimits = xlim;
    %         text(xlimits(1), 0.95*max(profile(fileNum).hotelPowerCalc), sprintf('Energy = %.0f Joules', profile(fileNum).totalHotelEnergy));
    %
    %         plotStateTransitions(tableData, DT)
    %     end
    %
    % %     % MC, print figure
    % %     print_setup(80,30)
    % %     print([path_plots_eng,'fig1_',fileList(fileNum).name,'.jpg'],'-djpeg');
    %
    %     figure
    %     hold on
    %     yyaxis left
    %     plot(profile(fileNum).pumpJoulesDT, profile(fileNum).pumpJoules)
    %     plot(profile(fileNum).hotelJoulesDT, profile(fileNum).hotelJoules)
    %     %plot(profile(fileNum).pumpEnergyCalcSumDT, profile(fileNum).pumpEnergyCalcSum)
    %     plot(profile(fileNum).hotelEnergyCalcSumDT, profile(fileNum).hotelEnergyCalcSum)
    %     yyaxis right
    %     plot(profile(fileNum).pumpCoulombsDT, profile(fileNum).pumpCoulombs)
    %     plot(profile(fileNum).hotelCoulombsDT, profile(fileNum).hotelCoulombs)
    %     legend
    %
    %     hold off
    %
    %
    %     %     if(plotEachProfile)
    %     %         figure
    %     %         plot(profile(fileNum).pumpEnergyCalcSumDT, profile(fileNum).pumpEnergyCalcSum, ...
    %     %             profile(fileNum).hotelEnergyCalcSumDT, profile(fileNum).hotelEnergyCalcSum, ...
    %     %             profile(fileNum).hotelEnergyCalcSumDT, (profile(fileNum).pumpEnergyCalcSum + profile(fileNum).hotelEnergyCalcSum));
    %     %         plotStateTransitions(tableData, DT);
    %     %     end
    %
    % %    % MC, print figure
    % %     print_setup(40,30)
    % %     print([path_plots_eng,'fig2_',fileList(fileNum).name,'.jpg'],'-djpeg');
    %
    %     figure
    %     hold on
    %
    %     %% Process UBX messages
    %     fprintf('%s\n', [datestr(now, 'HH:MM:SS ') 'Process GPS Messages ...']);
    %
    %     ubxFind = strfind(tableData.GPS, 'UBX');
    %     ubxMsgIndexes = find(~cellfun('isempty', ubxFind));
    %     profile(fileNum).ubxMsgs = tableData.GPS(ubxMsgIndexes);
    %     profile(fileNum).ubxMsgsDT = DT(ubxMsgIndexes);
    %     for i=1:length(profile(fileNum).ubxMsgs)
    %         fprintf('%s %s\n', profile(fileNum).ubxMsgsDT(i),  profile(fileNum).ubxMsgs(i));
    %     end
    %     fprintf('\n');
    %
    %     %% Process SBDIX messages
    %     fprintf('%s\n', [datestr(now, 'HH:MM:SS ') 'Process Iridium Messages ...']);
    %
    %     %     sbdixFind = strfind(tableData.IridModem, '+SBDIX:');
    %     %     sbdixMsgIndexes = find(~cellfun('isempty', sbdixFind));
    %     %     profile(fileNum).sbdixMsgs = tableData.IridModem(sbdixMsgIndexes);
    %     %     profile(fileNum).sbdixMsgsDT = DT(sbdixMsgIndexes);
    %
    %     %% Find error messages
    %     fprintf('%s\n', [datestr(now, 'HH:MM:SS ') 'Find Error Messages ...']);
    %
    %     errorFind = strfind(tableData.comment, '[ERROR]');
    %     errorMsgIndexes = find(~cellfun(@isempty, errorFind));
    %     profile(fileNum).errorMsg = tableData.comment(errorMsgIndexes);
    %     profile(fileNum).errorMsgDT = DT(errorMsgIndexes);
    %
    %     fprintf(1, 'Errors in %d files processed: \r', i );
    %
    %     for i = 1:length(errorMsgIndexes)
    %         fprintf(1, 'Error %d at time %s : %s \r', i, datestr(profile(fileNum).errorMsgDT(i)), profile(fileNum).errorMsg(i));
    %     end
    %
    %% Other stuff
    %         figure(numSatsFig);
    %         clf;
    %         try
    %             plot(numSatstime, numSats, '+-');
    %             title('Number of GPS Satellites vs time')
    %         catch me
    %             title('Number of GPS Satellites not processed');
    %         end
    %
    %         figure(numIridiumTrysFig);
    %         clf;
    %         try
    %             plot(numSBDIXtime, numSBDIX, '+-');
    %             title('Number of SBDIX Attempts vs time')
    %         catch me
    %             title('Number of SBDIX Attempts not processed');
    %         end

    %        fprintf(1, 'Enter 2 points on descend plot and 2 points on ascend plots to calculate platform velocity');
    %        [x y] = ginput(4)
    %
    %        dt = 24 * 3600 * datenum(pressureDT(round(x(1)):round(x(2))));
    %        dt = dt - dt(1);
    %        [p, s] = polyfit(dt, pressure(round(x(1)):round(x(2))), 1);
    %        fprintf(1, 'platform velocity (slope) between first pair of points = %f    ', p);
    %        dt = 24 * 3600 * datenum(pressureDT(round(x(3)):round(x(4))));
    %        dt = dt - dt(1);
    %        [p, s] = polyfit(dt, pressure(round(x(3)):round(x(4))), 1);
    %        fprintf(1,'\r\n');
    %        fprintf(1, 'platform velocity (slope) between second pair of points = %f  ', p);

    close all

    TTif(fileNum,:)=[DT(3) DT(end)];
    TT(fileNum)=(datenum(DT(end))-datenum(DT(3)))*1440; % Total Time [minutes]


    %% ------------------------ MC plots --------------------------------------
    % Volume of oil and Bellows Positions vs time

    close all

    % Bellows Position

    subplot(211)
    x=profile(fileNum).bellowsPositionDT;
    y=profile(fileNum).bellowsPosition;

    plot(x,y,'o')
    datetick

    hold on
    % Volume of oil displaced
    % integrate JV over time, between instances when pump is on/off

    subplot(212)
    x=profile(fileNum).pumpJoulesDT;
    y=profile(fileNum).pumpJoules;
%% 

    plot(x,y,'o')


    %% MC, Integrate Jog Velocity

    % Lines to identify
    %     % commtxt=char(string(table2array(tableData(7221,5))));
    %     commtxt=char(string(table2array(tableData(7221,5))));
    %     strcmp(commtxt(1:11),"BE::sendJV:") % pump on
    %
    %     commtxt=char(string(table2array(tableData(10086,5))));
    %     strcmp(commtxt(1:11),"Normal State Exit") % pump off
    %
    %     % Note: JV stands for Jog Velocity (units revolutions/s or counts/s)
    %     % Gene will pass on the [ml of oil]/count conversion.
    %     % For now, integrating between pump_on_pump_off lines
    %
    %% Loop trhough the table

    nl=height(tableData); % number of lines in table

    pump_on=zeros(nl,3); pump_on(:,2:3)=NaN;
    pump_off=zeros(nl,2); pump_off(:,2)=NaN;

    flag_up=logical(zeros(nl,1));

    n=1;
    for l=1:nl

        state=table2array(tableData(l,4));

        if state=="Ascend" && ~ismissing(table2array(tableData(l,5)))
            commtxt=char(string(table2array(tableData(l,5)))); % comment column

            flag_up(l)=true; % flag for lines with "Ascend" data

            if length(commtxt)>=11 && strcmp(commtxt(1:11),"BE::sendJV:")
                pump_on(l,1)=n; % using a different number
                pump_on(l,2)=str2num(commtxt(12:end)); % JV values in counts/sec
                pump_on(l,3)=datenum(DT(l)); % matlab time
            elseif length(commtxt)==17 && strcmp(commtxt,"Normal State Exit")
                pump_off(l,1)=n;
                pump_off(l,2)=datenum(DT(l));
                n=n+1;
            end
        end

        clear commtxt
    end


    %% Plot pump_on/pump_off switch
    close all

    plot(pump_on(~isnan(pump_on(:,3)),1),pump_on(~isnan(pump_on(:,3)),3),'g');
    hold on
    plot(pump_off(~isnan(pump_off(:,2)),1),pump_off(~isnan(pump_off(:,2)),2),'ro');
    datetick('y','mmm-dd HH:MM','keepticks','keeplimits')

    legend('ON','OFF')
    title('Ascend Mode')

    %     print_setup(20,30)
    %     print([path_plots_eng,'PumpONOFF_upcast_',fileList(fileNum).name,'.jpg'],'-djpeg');
    %
    %% Integrate JV from 1st pump_on to

    dt=diff(pump_on(~isnan(pump_on(:,3)),3)).*86400; % in fraction of day
    JV=pump_on(~isnan(pump_on(:,3)),2);

    Vb=sum(JV(2:end).*dt); % total counts (from beginning to end of profile)

    Vb=cumsum(JV(2:end).*dt); % total counts (from beginning to end of profile)

    cumsum(Vb);

    %% Total "profile duration", from descend to ascend

    % TT(fileNum,:)=[DT(2) DT(end)];

    %% MC, Save profile data (upcast) in structure array

    iP   = find(strcmp('pressure', varNames) == 1); % CTD P
    iB   = find(strcmp('bellowsPos', varNames) == 1); % Bellows Position

    profile_up.pressure=table2array(tableData(flag_up,iP));
    profile_up.pressure(profile_up.pressure<=1)=NaN; % confirmed with Gene his is sensible to do

    profile_up.bellowsPosition=table2array(tableData(flag_up,iB));
    profile_up.bellowsPosition(profile_up.bellowsPosition<=1)=NaN; % confirmed with Gene this is sensible to do

    profile_up.JV=pump_on(flag_up,2); % JV
    profile_up.DT=datenum(DT(flag_up)); % DateTime


    % Platform Velocity (as computed by Gene)
    iPlowpass   = find(strcmp('LPFPress', varNames) == 1); % LowPassFiltered Pres
    iWlowpass   = find(strcmp('LPFVel', varNames) == 1); % LowPassFiltered Vel
    iWdiff   = find(strcmp('diffVel', varNames) == 1); % W from dP/dt
    iWdiff2   = find(strcmp('diff2Vel', varNames) == 1); % W from dP/dt (filtered?)
    iW = find(strcmp('slopeOfLPFPress', varNames) == 1); % W from slope of P vs t
    % {'slopeOfLPFPress'}

    profile_up.pressureLP=table2array(tableData(flag_up,iPlowpass));
    profile_up.WLP=table2array(tableData(flag_up,iWlowpass));
    profile_up.Wdiff=table2array(tableData(flag_up,iWdiff));
    profile_up.Wdiff2=table2array(tableData(flag_up,iWdiff2));
    profile_up.W=table2array(tableData(flag_up,iW));

    %% MC, Plot Bellows Position and JV for the Ascent MOde only

    % close all
    figure(2)
    subplot(121)
    plot(profile_up.bellowsPosition,profile_up.pressure,'o'); axis ij
    ylim([0 maxp])
    ylabel(' Pressure [dbar]')
    xlabel(' Bellows Position')
    title('Ascend Mode')

    subplot(122)
    plot(profile_up.JV(1:end-1),profile_up.pressure(2:end),'o'); axis ij
    ylim([0 maxp])
    ylabel(' Pressure [dbar]')
    xlabel(' Jog Velocity [counts/s]')
    title('Ascend Mode')

    %     print_setup(40,30)
    %     print([path_plots_eng,'BellowsPos&JogVelocity_upcast_',fileList(fileNum).name,'.jpg'],'-djpeg');

    %%
    close all
    subplot(121)
    plot(profile_up.bellowsPosition,profile_up.DT,'o');
    ylabel(' Time')
    xlabel(' Bellows Position [IB % Full] ')
    ylim(range(profile_up.DT))
    datetick('y','keeplimits')
    title('Ascend Mode')

    subplot(122)
    plot(profile_up.JV(1:end),profile_up.DT(1:end),'o');
    ylabel(' Time')
    xlabel(' Jog Velocity [counts/s]')
    ylim(range(profile_up.DT))
    datetick('y','keeplimits')
    title('Ascend Mode')

    %     print_setup(40,30)
    %     print([path_plots_eng,'BellowsPos&JogVelocity_vs_time_upcast_',fileList(fileNum).name,'.jpg'],'-djpeg');

    %%
    figure(4)
    yyaxis left
    % plot(profile_up.bellowsPosition,profile_up.pressure,'o');
    plot(profile_up.DT,profile_up.bellowsPosition,'o');
    xlabel(' Time')
    ylabel(' Bellows Position [IB % Full] ')
    %
    datetick('x')
    title('Ascend Mode')

    yyaxis right
    % plot(profile_up.JV(1:end),profile_up.pressure(1:end),'o');
    plot(profile_up.DT,profile_up.pressure,'o');
    xlabel(' Time')
    ylabel(' Pressure [dbar]')
    datetick('x')

    xlim(range(profile_up.DT))
    % title('Ascend Mode')

    %     print_setup(40,30)
    %     print([path_plots_eng,'BellowsPos&Pressure_vs_time_upcast_',fileList(fileNum).name,'.jpg'],'-djpeg');
    %

    %% MC, SAVE .ENG data for upward profile in .mat file

    ENG.fileName=profile.fileName; % Data Source
    ENG.P1hz=profile_up.pressure; % Pressure (1Hz?) in dbar
    ENG.BellPos=profile_up.bellowsPosition; % Bellows Position (%)
    ENG.JV=profile_up.JV; % Jogg Velocity
    ENG.DT=profile_up.DT; % DateTime (matlab time)

    % Other parameters (i.e. platform velocity, computed by Gene)
    ENG.Plp=profile_up.pressureLP; % ask Gene the details
    ENG.Wlp=profile_up.WLP; %
    ENG.Wdiff=profile_up.Wdiff; %
    ENG.Wdiff2=profile_up.Wdiff; %
    ENG.W=profile_up.W; %

    save([path_mat,matfiles(fileNum).name],'ENG','-append')
    display(['Appending ENG structure to files in ',path_mat])
end


%% MC, estimate ranges for BBVolts, CanPressure and CanHumidity
%% 

%% 
CanHumidity=profile.humidity; CanHumidity(CanHumidity==0)=NaN;
rangeH=range(CanHumidity);
display(['Can Humidity range = ',num2str(rangeH)])

CanPressure=profile.intPressure;
rangeP=range(CanPressure);
display(['Can Pressure range = ',num2str(rangeP)])
%% 

%%
%% Waterfall plots
if(plotWaterfalls)
    fprintf('%s\n', [datestr(now, 'HH:MM:SS ') 'Plotting waterfalls ...']);

    figX = 0.05;
    figXIncr = 0.1;
    figY = 0.3;
    figYIncr = -0.05;

    %Depth
    fig = figure('Name','Depth','NumberTitle','off');
    fig.Units = 'normalized';
    fig.Position = [figX figY 0.6 0.6];
    figX = figX + figXIncr;
    figY = figY + figYIncr;
    title('Pressure')
    hold on;

    for profileNum = 1:length(profile)
        for i=1:length(profile(profileNum).pressureDT)
            x(i) = profile(profileNum).pressureDT(1);
        end
        y = profile(profileNum).pressureDT - profile(profileNum).pressureDT(1);
        plot3(x, y, profile(profileNum).pressure)
        clear x;
    end
    xlabel('Profile Start time');
    ylabel('Time from Profile Start')
    zlabel('Depth (dBar)')
    grid
    view(3)

    % Bellows Position
    fig = figure('Name','Bellows Pos','NumberTitle','off');
    fig.Units = 'normalized';
    fig.Position = [figX figY 0.6 0.6];
    figX = figX + figXIncr;
    figY = figY + figYIncr;
    title('Bellows Position')
    hold on;

    for profileNum = 1:length(profile)
        for i=1:length(profile(profileNum).bellowsPositionDT)
            x(i) = profile(profileNum).bellowsPositionDT(1);
        end
        y = profile(profileNum).bellowsPositionDT - profile(profileNum).bellowsPositionDT(1);
        plot3(x, y, profile(profileNum).bellowsPosition)
        clear x;
    end
    xlabel('Profile Start time');
    ylabel('Time from Profile Start')
    zlabel('Bellows Position (IB % Full)')
    grid
    view(3)

    %Pump Volts
    fig = figure('Name','Pump Volts','NumberTitle','off');
    fig.Units = 'normalized';
    fig.Position = [figX figY 0.6 0.6];
    figX = figX + figXIncr;
    figY = figY + figYIncr;
    title('Pump Bus Volts')
    hold on;

    for profileNum = 1:length(profile)
        for i=1:length(profile(profileNum).pumpVoltsDT)
            x(i) = profile(profileNum).pumpVoltsDT(1);
        end
        y = profile(profileNum).pumpVoltsDT - profile(profileNum).pumpVoltsDT(1);
        plot3(x, y, profile(profileNum).pumpVolts)
        clear x;
    end

    xlabel('Profile Start time');
    ylabel('Time from Profile Start')
    zlabel('Pump Bus Volts')
    grid
    view(3)

    %Pump power
    fig = figure('Name','Pump Power','NumberTitle','off');
    fig.Units = 'normalized';
    fig.Position = [figX figY 0.6 0.6];
    figX = figX + figXIncr;
    figY = figY + figYIncr;
    title('Pump Bus Power')
    hold on;

    for profileNum = 1:length(profile)
        for i=1:length(profile(profileNum).pumpPowerDT)
            x(i) = profile(profileNum).pumpPowerDT(1);
        end
        y = profile(profileNum).pumpPowerDT - profile(profileNum).pumpPowerDT(1);
        plot3(x, y, profile(profileNum).pumpPower)
        clear x;
    end
    xlabel('Profile Start time');
    ylabel('Time from Profile Start')
    zlabel('Pump Bus Power (Watts)')
    grid
    view(3)

    %Internal Pressure
    fig = figure('Name','Internal Pressure','NumberTitle','off');
    fig.Units = 'normalized';
    fig.Position = [figX figY 0.6 0.6];
    figX = figX + figXIncr;
    figY = figY + figYIncr;
    title('Internal Pressure (mBar)')
    hold on;

    for profileNum = 1:length(profile)
        for i=1:length(profile(profileNum).intPressureDT)
            x(i) = profile(profileNum).intPressureDT(1);
        end
        y = profile(profileNum).intPressureDT - profile(profileNum).intPressureDT(1);
        plot3(x, y, profile(profileNum).intPressure)
        clear x;

        y = profile(profileNum).pumpVoltsDT - profile(profileNum).pumpVoltsDT(1);
        plot3(x, y, profile(profileNum).pumpVolts)
        clear x;
    end
    xlabel('Profile Start time');
    ylabel('Time from Profile Start')
    zlabel('Pump Bus Volts')
    grid
    view(3)

    %Pump power
    fig = figure('Name','Pump Power (calculated)','NumberTitle','off');
    fig.Units = 'normalized';
    fig.Position = [figX figY 0.6 0.6];
    figX = figX + figXIncr;
    figY = figY + figYIncr;
    title('Pump Bus Power (calculated)')
    hold on;

    for profileNum = 1:length(profile)
        for i=1:length(profile(profileNum).pumpPowerCalcDT)
            x(i) = profile(profileNum).pumpPowerCalcDT(1);
        end
        y = profile(profileNum).pumpPowerCalcDT - profile(profileNum).pumpPowerCalcDT(1);
        plot3(x, y, profile(profileNum).pumpPowerCalc)
        clear x;
    end
    xlabel('Profile Start time');
    ylabel('Time from Profile Start')
    zlabel('Pump Bus Power (calculated Watts)')
    grid
    view(3)

    %Internal Pressure
    fig = figure('Name','Internal Pressure','NumberTitle','off');
    fig.Units = 'normalized';
    fig.Position = [figX figY 0.6 0.6];
    figX = figX + figXIncr;
    figY = figY + figYIncr;
    title('Internal Pressure (mBar)')
    hold on;

    for profileNum = 1:length(profile)
        for i=1:length(profile(profileNum).intPressureDT)
            x(i) = profile(profileNum).intPressureDT(1);

        end
        xlabel('Profile Start time');
        ylabel('Time from Profile Start')
        zlabel('Internal Pressure (mBar)')
        grid
        view(3)

        % Humidity
        fig = figure('Name','Humidity','NumberTitle','off');
        fig.Units = 'normalized';
        fig.Position = [figX figY 0.6 0.6];
        figX = figX + figXIncr;
        figY = figY + figYIncr;
        title('Internal Humidity (%RH)')
        hold on;

        for profileNum = 1:length(profile)
            for i=1:length(profile(profileNum).humidityDT)
                x(i) = profile(profileNum).humidityDT(1);
            end
            y = profile(profileNum).humidityDT - profile(profileNum).humidityDT(1);
            plot3(x, y, profile(profileNum).humidity)
            clear x;
        end

        xlabel('Profile Start time');
        ylabel('Time from Profile Start')
        zlabel('Humidity (%RH)')
        grid
        view(3)

        %Pump Energy sum
        fig = figure('Name','Pump Cumul. Energy','NumberTitle','off');
        fig.Units = 'normalized';
        fig.Position = [figX figY 0.6 0.6];
        figX = figX + figXIncr;
        figY = figY + figYIncr;
        title('Pump Cumulative Energy (Joules)')
        hold on;

        for profileNum = 1:length(profile)
            for i=1:length(profile(profileNum).pumpEnergySumDT)
                x(i) = profile(profileNum).pumpEnergySumDT(1);
            end
            y = profile(profileNum).pumpEnergySumDT - profile(profileNum).pumpEnergySumDT(1);
            plot3(x, y, profile(profileNum).pumpEnergySum)
            clear x;
        end
        xlabel('Profile Start time');
        ylabel('Time from Profile Start')
        zlabel('Pump Cumulative Energy (Joules')
        grid
        view(3)

        %Hotel Energy sum
        fig = figure('Name','Hotel Cumul. Energy','NumberTitle','off');
        fig.Units = 'normalized';
        fig.Position = [figX figY 0.6 0.6];
        figX = figX + figXIncr;
        figY = figY + figYIncr;
        title('Hotel Cumulative Energy (Joules)')
        hold on;

        for profileNum = 1:length(profile)
            for i=1:length(profile(profileNum).hotelEnergySumDT)
                x(i) = profile(profileNum).hotelEnergySumDT(1);
            end
            y = profile(profileNum).hotelEnergySumDT - profile(profileNum).hotelEnergySumDT(1);
            plot3(x, y, profile(profileNum).hotelEnergySum)
            clear x;
        end
        xlabel('Profile Start time');
        ylabel('Time from Profile Start')
        zlabel('Hotel Cumulative Energy (Joules')
        grid
        view(3)

        y = profile(profileNum).humidityDT - profile(profileNum).humidityDT(1);
        plot3(x, y, profile(profileNum).humidity)
        clear x;
    end
    xlabel('Profile Start time');
    ylabel('Time from Profile Start')
    zlabel('Humidity (%RH)')
    grid
    view(3)

    %Pump Energy sum
    fig = figure('Name','Pump Cumul. Energy','NumberTitle','off');
    fig.Units = 'normalized';
    fig.Position = [figX figY 0.6 0.6];
    figX = figX + figXIncr;
    figY = figY + figYIncr;
    title('Pump Cumulative Energy (Joules)')
    hold on;

    for profileNum = 1:length(profile)
        for i=1:length(profile(profileNum).pumpEnergyCalcSumDT)
            x(i) = profile(profileNum).pumpEnergyCalcSumDT(1);
        end
        y = profile(profileNum).pumpEnergyCalcSumDT - profile(profileNum).pumpEnergyCalcSumDT(1);
        plot3(x, y, profile(profileNum).pumpEnergyCalcSum)
        clear x;
    end
    xlabel('Profile Start time');
    ylabel('Time from Profile Start')
    zlabel('Pump Cumulative Energy (Joules')
    grid
    view(3)

    %Hotel Energy sum
    fig = figure('Name','Hotel Cumul. Energy','NumberTitle','off');
    fig.Units = 'normalized';
    fig.Position = [figX figY 0.6 0.6];
    figX = figX + figXIncr;
    figY = figY + figYIncr;
    title('Hotel Cumulative Energy (Joules)')
    hold on;

    for profileNum = 1:length(profile)
        for i=1:length(profile(profileNum).hotelEnergyCalcSumDT)
            x(i) = profile(profileNum).hotelEnergyCalcSumDT(1);
        end
        y = profile(profileNum).hotelEnergyCalcSumDT - profile(profileNum).hotelEnergyCalcSumDT(1);
        plot3(x, y, (profileNum.hotelEnergyCalcSum))
        clear x;

    end
end

hold off;

function plotStateTransitions(tableData, DT)

%figure(pressureFig);
ylimits = ylim;

textY1 = ylimits(1) + 0.1 * (ylimits(2) - ylimits(1));
textY2 = ylimits(1) + 0.9 * (ylimits(2) - ylimits(1));

indexes = find(strcmp(tableData.state, 'Surface Ops') == 1);
if(~isempty(indexes))
    index = indexes(1);
    plot( [DT(index) DT(index)], [ylimits(1) ylimits(2)], 'g-.');
    text( DT(index), textY1, 'Surface Ops')
end

indexes = find(strcmp(tableData.state, 'AB Dec Buoycy Fast') == 1);
if(~isempty(indexes))
    index = indexes(1);
    plot( [DT(index) DT(index)], [ylimits(1) ylimits(2)], 'g-.');
    text( DT(index), textY2, 'Dec Buoycy')
end

indexes = find(strcmp(tableData.state, 'Descend')== 1);
if(~isempty(indexes))
    index = indexes(1);
    plot( [DT(index) DT(index)], [ylimits(1) ylimits(2)], 'g-.');
    text( DT(index), textY1, 'Descend')
end

indexes = find(strcmp(tableData.state, 'Anchor')== 1);
if(~isempty(indexes))
    index = indexes(1);
    plot( [DT(index) DT(index)], [ylimits(1) ylimits(2)], 'g-.');
    text( DT(index), textY2, 'Anchor')
end

indexes = find(strcmp(tableData.state, 'Park')== 1);
if(~isempty(indexes))
    index = indexes(1);
    plot( [DT(index) DT(index)], [ylimits(1) ylimits(2)], 'g-.');
    text( DT(index), textY1, 'Park')
end

indexes = find(strcmp(tableData.state, 'Start CP Mode')== 1);
if(~isempty(indexes))
    index = indexes(1);
    plot( [DT(index) DT(index)], [ylimits(1) ylimits(2)], 'g-.');
    text( DT(index), textY2, 'Start CP')
end

indexes = find(strcmp(tableData.state, 'Deanchor')== 1);
if(~isempty(indexes))
    index = indexes(1);
    plot( [DT(index) DT(index)], [ylimits(1) ylimits(2)], 'g-.');
    text( DT(index), textY1, 'Deanchor')
end

indexes = find(strcmp(tableData.state, 'Ascend')== 1);
if(~isempty(indexes))
    index = indexes(1);
    plot( [DT(index) DT(index)], [ylimits(1) ylimits(2)], 'g-.');
    text( DT(index), textY2, 'Ascend')
end

indexes = find(strcmp(tableData.state, 'Init Profile') == 1);
if(~isempty(indexes))
    index = indexes(1);
    plot( [DT(index) DT(index)], [ylimits(1) ylimits(2)], 'g-.');
    text( DT(index), textY1, 'Init Profile')
end

indexes = find(strcmp(tableData.state, 'Surface Ops Set Bellows') == 1);
if(~isempty(indexes))
    index = indexes(1);
    plot( [DT(index) DT(index)], [ylimits(1) ylimits(2)], 'g-.');
    text( DT(index), textY2, 'SO Set Bel')
end

indexes = find(strcmp(tableData.state, 'Dump CP Data') == 1);
if(~isempty(indexes))
    index = indexes(1);
    plot( [DT(index) DT(index)], [ylimits(1) ylimits(2)], 'g-.');
    text( DT(index), textY1, 'Dump CP')
end

indexes = find(strcmp(tableData.state, 'RecoveryMode') == 1);
if(~isempty(indexes))
    index = indexes(1);
    plot( [DT(index) DT(index)], [ylimits(1) ylimits(2)], 'g-.');
    text( DT(index), textY1, 'Recovery')
end

indexes = find(strcmp(tableData.state, 'Init Mission') == 1);
if(~isempty(indexes))
    index = indexes(1);
    plot( [DT(index) DT(index)], [ylimits(1) ylimits(2)], 'g-.');
    text( DT(index), textY1, 'Init Mission')
end

end

function createFigure(fileName, figNum)
figX = 0.05;
figXIncr = 0.05;
figY = 0.6;
figYIncr = -0.05;

pressureFig = figure;
pressureFig.Units = 'normalized';
pressureFig.Position = [figX figY 0.6 0.4];
title('Pressure & Bellows Position vs time');
hold on;
grid on;

pumpVoltsFig = figure;
pumpVoltsFig.Units = 'normalized';
figX = figX + figXIncr;
figY = figY + figYIncr;
pumpVoltsFig.Position = [figX figY 0.6 0.4];
title('Pump Voltage & Current vs time');
hold on;
grid on;

pumpSpeedFig = figure;
pumpSpeedFig.Units = 'normalized';
figX = figX + figXIncr;
figY = figY + figYIncr;
pumpSpeedFig.Position = [figX figY 0.6 0.4];
title('Pressure & Pump Speed vs time');
hold on;
grid on;

pressureOnlyFig = figure;
pressureOnlyFig.Units = 'normalized';
figX = figX + figXIncr;
figY = figY + figYIncr;
pressureOnlyFig.Position = [figX figY 0.6 0.4];
title('Pressure & Pump Speed vs time');
hold on;
grid on;
end









