clear all
close all

opts = delimitedTextImportOptions("NumVariables", 30);

% Specify range and delimiter
opts.DataLines = [25, Inf];
opts.Delimiter = ",";

% Specify column names and types
%opts.VariableNames = ["LabviewDateTime", "ThermistorTemp", "BathTemp", "BathSetPoint", "PressuredBar", "MaxFlow", "FlowRate", "PumpVolume", "CRC", "ID", "DATETIME", "UTCSECS", "CTDDEPTH", "CTDTEMP", "CTDSALINITY", "SAMPLECOUNTER", "RESETCOUNTER", "ERRORCOUNTER", "HOUSINGTEMPERATURE", "HOUSINGHUMIDITY", "INPUTVOLTAGE", "INPUTCURRENT", "FOOBARPHVALUE", "BACKUPBATTERY", "VRSVOLTS", "VRSSTDDEV", "VKVOLTS", "VKSTDDEV", "IKAMPS", "IBAMPS"];
%opts.VariableNames = ["dtLV",   "thermTemp", "bathTemp", "bathSP", "pressureDBar", "maxFlow", "flowRate", "cycleNum", "pumpVolume", "CRC", "ID", "dateTime", "utcSecs", "ctdDepth", "ctdTemp", "ctdSalinity", "sampleCounter", "resetCounter", "errorCounter", "housingTemp", "housingHumidity", "inputVoltage", "inputCurrent", "foobarPHValue", "backupBattery", "vrsVolts", "vrsStdDev", "vkVolts", "vkStdDev", "ikAmps", "ibAmps"];
%opts.VariableTypes = ["datetime", "double", "double",   "double", "double",         "double", "double",     "double",   "double", "string", "categorical", "datetime", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double"];
%For the new VI with SetTemp, SetPressure and StepNum
opts.VariableNames = ["dtLV",   "thermTemp", "bathTemp", "bathSP", "pressureDBar", "maxFlow", "flowRate", "pumpVolume", "CRC", "ID", "dateTime", "utcSecs", "ctdDepth", "ctdTemp", "ctdSalinity", "sampleCounter", "resetCounter", "errorCounter", "housingTemp", "housingHumidity", "inputVoltage", "inputCurrent", "foobarPHValue", "backupBattery", "vrsVolts", "vrsStdDev", "vkVolts", "vkStdDev", "ikAmps", "ibAmps", "SetTemp", "SetPressure", "StepNum"];
opts.VariableTypes = ["datetime", "double", "double",   "double", "double",         "double", "double",    "double", "string", "categorical", "datetime", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double"];

% Specify file level properties
opts.ExtraColumnsRule = "ignore";
opts.EmptyLineRule = "read";

% Specify variable properties
opts = setvaropts(opts, "CRC", "WhitespaceRule", "preserve");
opts = setvaropts(opts, ["CRC", "ID"], "EmptyFieldRule", "auto");
%opts = setvaropts(opts, "dtLV", "InputFormat", "MM/dd/uuuu hh:mm:ss aa");
%opts = setvaropts(opts, "dateTime", "InputFormat", "MM/dd/uuuu hh:mm:ss aa");

%% Process files
path = '\\atlas\ProjectLibrary\ChemSensors\GeneM\GDF pH Sensor\Data\Cal Data\Calcs\DataFiles\'
fullPath = strcat(path, '*.csv');

dataFiles = dir(strcat(path, '*.csv'));
fileDatenums = dataFiles.date;
[sortedTime, i] = sort([dataFiles.datenum]);
sortedDataFileList = dataFiles(i);

for i=1:length(dataFiles)
    % Import the data
    disp("Processing: " + dataFiles(1).name);
    dataTable{i} = readtable(strcat(path, sortedDataFileList(i).name), opts);
end

% Create overview plot
createOverviewPlot2(dataTable, sortedDataFileList, 1);

% Separate out cycles
%[numCycles] = getCycles2(dataTable);
cycleNum = 0;
diffNum = 6;

dataTable{1}.cycleNum = zeros(length(dataTable{1}.pressureDBar),1);
pressureDBar = dataTable{1}.pressureDBar;

for i=1+diffNum:length(pressureDBar)
    diff = pressureDBar(i-diffNum)-pressureDBar(i);
    if( abs(diff) > 1.0)
        if(dataTable{1}.cycleNum(i-1)) == 0
            cycleNum = cycleNum + 1;
        end
        dataTable{1}.cycleNum(i) = cycleNum;
    else
        dataTable{1}.cycleNum(i-6:i) = 0;
    end
end

% Waterfall plot
%plotPressureWaterfall(cycle, numCycles)

figure('name', "Vrs vs Time", 'NumberTitle', 'off');
plot(dataTable{1}.dtLV, dataTable{1}.vrsVolts)
title('time vs Vrs')

%Plot individual cycles
numBins = 100;{i}
xIncr = 0.2;
yIncr = 0.3;
left = 0.05;
bottom = 0.65;
width = 0.25;
height = 0.25;

polyOrder = 3;

for i=2:2:max(dataTable{1}.cycleNum)
    cycleIndex = i/2;
    hold on

    indexes = find(dataTable{1}.cycleNum == i);
    cycle{cycleIndex}.pIncr = dataTable{1}.pressureDBar(indexes);
    cycle{cycleIndex}.VrsIncr = dataTable{1}.vrsVolts(indexes);
    cycle{cycleIndex}.dtIncr = dataTable{1}.dtLV(indexes);

    lastwarn('')
    [cycle{cycleIndex}.polyFitCoefIncr, cycle{cycleIndex}.polyFitSIncr, mu] = polyfit(cycle{cycleIndex}.pIncr, cycle{cycleIndex}.VrsIncr, polyOrder);
    [warnMsg, warnId] = lastwarn;
    if isempty(warnMsg)

    cycle{cycleIndex}.VrsFitIncr = polyval(cycle{cycleIndex}.polyFitCoefIncr, cycle{cycleIndex}.pIncr);
    residualIncr = cycle{cycleIndex}.VrsFitIncr - cycle{cycleIndex}.VrsIncr;
    cycle{cycleIndex}.maxResidualIncr = max(abs(residualIncr));
    cycle{cycleIndex}.meanResidualIncr = mean(abs(residualIncr));
    cycle{cycleIndex}.RMSResidualIncr = rms(residualIncr);
    cycle{cycleIndex}.STDResidualIncr = std(residualIncr);

    indexes = find(dataTable{1}.cycleNum == i+1);
    cycle{cycleIndex}.pDecr = dataTable{1}.pressureDBar(indexes);
    cycle{cycleIndex}.VrsDecr = dataTable{1}.vrsVolts(indexes);
    cycle{cycleIndex}.dtDecr = dataTable{1}.dtLV(indexes);

    [cycle{cycleIndex}.polyFitCoefDecr, cycle{cycleIndex}.polyFitSDecr, mu] = polyfit(cycle{cycleIndex}.pDecr, cycle{cycleIndex}.VrsDecr, polyOrder);
    cycle{cycleIndex}.VrsFitDecr = polyval(cycle{cycleIndex}.polyFitCoefDecr, cycle{cycleIndex}.pDecr);
    residualDecr = cycle{cycleIndex}.VrsFitDecr - cycle{cycleIndex}.VrsDecr;
    cycle{cycleIndex}.maxResidualDecr = max(abs(residualDecr));
    cycle{cycleIndex}.meanResidualDecr = mean(abs(residualDecr));
    cycle{cycleIndex}.RMSResidualDecr = rms(residualDecr);
    cycle{cycleIndex}.STDResidualDecr = std(residualDecr);

    cycle{cycleIndex}.pBoth = [cycle{cycleIndex}.pIncr; cycle{cycleIndex}.pDecr];
    cycle{cycleIndex}.VrsBoth = [cycle{cycleIndex}.VrsIncr; cycle{cycleIndex}.VrsDecr];
    cycle{cycleIndex}.dtBoth = [cycle{cycleIndex}.dtIncr; cycle{cycleIndex}.dtDecr];

    cycle{cycleIndex}.polyFitCoefBoth = polyfit(cycle{cycleIndex}.pBoth, cycle{cycleIndex}.VrsBoth, polyOrder);
    [cycle{cycleIndex}.polyFitCoefIB, cycle{cycleIndex}.polyFitSIncr, mu] = polyval(cycle{cycleIndex}.polyFitCoefBoth, cycle{cycleIndex}.pBoth);
    residualBoth = cycle{cycleIndex}.VrsFitBoth - cycle{cycleIndex}.VrsBoth;
    cycle{cycleIndex}.maxResidualBoth = max(abs(cycle{cycleIndex}.VrsFitBoth - cycle{cycleIndex}.VrsBoth));
    cycle{cycleIndex}.meanResidualBoth = mean(abs(cycle{cycleIndex}.VrsFitBoth - cycle{cycleIndex}.VrsBoth));
    cycle{cycleIndex}.RMSResidualBoth = rms(residualBoth);
    cycle{cycleIndex}.STDResidualBoth = std(residualBoth);

    left = left + xIncr;
    if(mod(i, 10) == 0)
        bottom = bottom - yIncr;
        left = 0.05;
    end
    else
        disp 'skipping this cycle'
        cycle{cycleIndex}.VrsFitIncr = NaN;
        cycle{cycleIndex}.maxResidualIncr = NaN;
        cycle{cycleIndex}.meanResidualIncr = NaN;

        cycle{cycleIndex}.pDecr = NaN;
        cycle{cycleIndex}.VrsDecr = NaN;
        cycle{cycleIndex}.polyFitCoefDecr = NaN;
        cycle{cycleIndex}.VrsFitDecr = NaN;
        cycle{cycleIndex}.maxResidualDecr = NaN;
        cycle{cycleIndex}.meanResidualDecr = NaN;

        cycle{cycleIndex}.pBoth = NaN;
        cycle{cycleIndex}.VrsBoth = NaN;
        cycle{cycleIndex}.polyFitCoefBoth = NaN;
        cycle{cycleIndex}.VrsFitBoth = NaN;
        cycle{cycleIndex}.maxResidualBoth = NaN;
        cycle{cycleIndex}.meanResidualBoth = NaN;
    end
end

for i=1:26
    x1(i) = cycle{1, i}.maxResidualIncr;
    x2(i) = cycle{1, i}.meanResidualIncr;
    x3(i) = cycle{1, i}.RMSResidualIncr;
    x4(i) = cycle{1, i}.STDResidualIncr;

    x5(i) = cycle{1, i}.maxResidualDecr;
    x6(i) = cycle{1, i}.meanResidualDecr;
    x7(i) = cycle{1, i}.RMSResidualDecr;
    x8(i) = cycle{1, i}.STDResidualDecr;

    x9(i) = cycle{1, i}.maxResidualBoth;
    x10(i) = cycle{1, i}.meanResidualBoth;
    x11(i) = cycle{1, i}.RMSResidualBoth;
    x12(i) = cycle{1, i}.STDResidualBoth;
end

figure
hold on
plot(x1)
plot(x2)
plot(x3)
plot(x4)
plot(x5)
plot(x6)
plot(x7)
plot(x8)
plot(x9)
plot(x10)
plot(x11)
plot(x12)

legend('max Incr P', 'mean Incr P', 'RMS Incr P', 'STD Incr P', ...
    'max Decr P', 'mean Decr P', 'RMS Decr P', 'STD Decr P', ...
    'max Both P', 'mean Both P', 'RMS Both P', 'STD Both P');


figure
hold on
colormap(jet);
colorValue = -5;
for i=1:length(cycle)
    x = datenum(cycle{i}.dtIncr)';
    y = cycle{i}.pIncr' ;
    z = cycle{i}.VrsIncr';

    %    lineColor = colorValue * ones(1, length(x));
    plot3(x, y, z, 'g');
    z = cycle{i}.VrsFitIncr';
    plot3(x, y, z, 'b');

 %   patch([x nan], [y nan], [z nan], [lineColor nan],'EdgeColor','interp','FaceColor','none')
end
set(gca, 'ydir', 'reverse');
zLimits = zlim;
set(gca, 'ztick', [zLimits(1):0.001:zLimits(2)+.001]);
view(3);
datetick('x', 'mmmdd HH:MM')
grid on
colorbar;
xStart = datenum(cycle{1}.dtIncr(1));
xEnd = datenum(cycle{length(cycle)}.dtIncr(end));
xlim([xStart xEnd]);


figure
hold on
colormap(jet);
colorValue = -5;
for i=1:length(cycle)   
    x = datenum(cycle{i}.dtDecr)';
    y = cycle{i}.pDecr' ;
    z = cycle{i}.VrsDecr';

%    lineColor = colorValue * ones(1, length(x));
    plot3(x, y, z, 'g');
    z = cycle{i}.VrsFitDecr';
    plot3(x, y, z, 'b');
 
 %   patch([x nan], [y nan], [z nan], [lineColor nan],'EdgeColor','interp','FaceColor','none')
end
set(gca, 'ydir', 'reverse');
zLimits = zlim;
set(gca, 'ztick', [zLimits(1):0.001:zLimits(2)+.001]);
view(3);
datetick('x', 'mmmdd HH:MM')
grid on
colorbar;
xStart = datenum(cycle{1}.dtDecr(1));
% xEnd = datenum(cycle{length(cycle)}.dtDecr(end));
% xlim([xStart xEnd]);


figure
hold on
colormap(jet);
colorValue = -5;
for i=1:length(cycle)   
    x = datenum(cycle{i}.dtBoth)';
    y = cycle{i}.pBoth' ;
    z = cycle{i}.VrsBoth';

%    lineColor = colorValue * ones(1, length(x));
    plot3(x, y, z, 'g');
    z = cycle{i}.VrsFitBoth';
    plot3(x, y, z, 'b');
 
 %   patch([x nan], [y nan], [z nan], [lineColor nan],'EdgeColor','interp','FaceColor','none')
end
set(gca, 'ydir', 'reverse');
zLimits = zlim;
set(gca, 'ztick', [zLimits(1):0.001:zLimits(2)+.001]);
view(3);
datetick('x', 'mmmdd HH:MM')
grid on
colorbar;
xStart = datenum(cycle{1}.dtBoth(1));
xEnd = datenum(cycle{length(cycle)}.dtBoth(end));
xlim([xStart xEnd]);