clear all
close all
polyOrder = 5;
numValuesIndex = 1;

matFileName = 'Lionix03V2SN001_1129and3023_ChB.mat'
load(matFileName);
data = [dataTable1; dataTable2];

%Get rid of data from pressure holds
dP1 = diff(data.PressuredBar);
indexes = find(abs(dP1) > 0.5);

pressure = data.PressuredBar(indexes);
Vrs = data.VRSVOLTS(indexes);
tempDegC = data.ThermistorTemp(indexes);
timeStamp = data.LabviewDateTime(indexes);

%Separate into individual pressure Incr and pressure Decr ramps
dP2 = diff(pressure);
signDP2 = sign(dP2);
dP3 = diff(signDP2);

indexBreaks = find(dP3 ~= 0);
indexBreaks = indexBreaks + 1;
indexBreaks = [1; indexBreaks];

lastBreak = length(indexBreaks)-1;

for j=1:lastBreak
    ramp(j).pressure = pressure(indexBreaks(j):indexBreaks(j+1)-1);
    ramp(j).Vrs = Vrs(indexBreaks(j):indexBreaks(j+1)-1);
    ramp(j).tempDegC = tempDegC(indexBreaks(j):indexBreaks(j+1)-1);
    ramp(j).timeStamp = timeStamp(indexBreaks(j):indexBreaks(j+1)-1);
    numValues(numValuesIndex) = length(ramp(j).pressure);
    numValuesIndex = numValuesIndex + 1;
end
%indexBreaks ends at the pressure peak before the last ramp Decr to
%zero. This add the last ramp to zero to the ramp data
ramp(j+1).pressure = pressure(indexBreaks(end):length(pressure));
ramp(j+1).Vrs = Vrs(indexBreaks(end):length(pressure));
ramp(j+1).tempDegC = tempDegC(indexBreaks(end):length(pressure));
ramp(j+1).timeStamp = timeStamp(indexBreaks(end):length(pressure));
numValues(numValuesIndex) = length(ramp(j+1).pressure);

allVrsIncr = [];
allVrsDecr = [];
allPressuresIncr = [];
allPressuresDecr = [];
allTempIncr = [];
allTempDecr = [];
allResidualsIncr = [];
allResidualsDecr = [];
allResidualsIncr = [];
allResidualsDecr = [];

cycleIndex = 1;
i = 1;
for j = 1:2:length(ramp)-1
    %disp(char(9) + "Processing cycle: " + cycleIndex)
    cycle(i, cycleIndex).pressureIncr = ramp(j).pressure;
    cycle(i, cycleIndex).VrsIncr = ramp(j).Vrs;
    cycle(i, cycleIndex).tempDegCUp = ramp(j).tempDegC;
    cycle(i, cycleIndex).timeStampIncr = ramp(j).timeStamp;

    cycle(i, cycleIndex).pressureDecr = ramp(j+1).pressure;
    cycle(i, cycleIndex).VrsDecr = ramp(j+1).Vrs;
    cycle(i, cycleIndex).tempDegCDown = ramp(j+1).tempDegC;
    cycle(i, cycleIndex).timeStampDecr = ramp(j+1).timeStamp;

    allPressuresIncr = [allPressuresIncr; cycle(i, cycleIndex).pressureIncr];
    allPressuresDecr = [allPressuresDecr; cycle(i, cycleIndex).pressureDecr];
    allVrsIncr = [allVrsIncr; cycle(i, cycleIndex).VrsIncr];
    allVrsDecr = [allVrsDecr; cycle(i, cycleIndex).VrsDecr];
    allTempIncr = [allTempIncr; cycle(i, cycleIndex).tempDegCUp];
    allTempDecr = [allTempDecr; cycle(i, cycleIndex).tempDegCDown];
 
    lastwarn('')
    [cycle(i, cycleIndex).polyPIncr, cycle(i, cycleIndex).polySIncr, mu] = polyfit(cycle(i, cycleIndex).pressureIncr, cycle(i, cycleIndex).VrsIncr, polyOrder);
    [warnMsg, warnId] = lastwarn;

    [cycle(i, cycleIndex).VrsFitIncr, cycle(i, cycleIndex).polyDeltaIncr] = polyval(cycle(i, cycleIndex).polyPIncr, cycle(i, cycleIndex).pressureIncr, cycle(i, cycleIndex).polySIncr, mu);
    cycle(i, cycleIndex).residualsIncr = cycle(i, cycleIndex).VrsFitIncr - cycle(i, cycleIndex).VrsIncr;
    cycle(i, cycleIndex).maxResidualIncr = max(abs(cycle(i, cycleIndex).residualsIncr));
    cycle(i, cycleIndex).meanResidualIncr = mean(abs(cycle(i, cycleIndex).residualsIncr));
    allResidualsIncr = [allResidualsIncr; cycle(i, cycleIndex).residualsIncr]; 

    lastwarn('')
    [cycle(i, cycleIndex).polyPDecr, cycle(i, cycleIndex).polySDecr, mu] = polyfit(cycle(i, cycleIndex).pressureDecr, cycle(i, cycleIndex).VrsDecr, polyOrder);
    [warnMsg, warnId] = lastwarn;

    [cycle(i, cycleIndex).VrsFitDecr, cycle(i, cycleIndex).polyDeltaDecr] = polyval(cycle(i, cycleIndex).polyPDecr, cycle(i, cycleIndex).pressureDecr, cycle(i, cycleIndex).polySDecr, mu);
    cycle(i, cycleIndex).residualsDecr = cycle(i, cycleIndex).VrsFitDecr - cycle(i, cycleIndex).VrsDecr;
    cycle(i, cycleIndex).maxResidualDecr = max(abs(cycle(i, cycleIndex).residualsDecr));
    cycle(i, cycleIndex).meanResidualDecr = mean(abs(cycle(i, cycleIndex).residualsDecr));
    allResidualsDecr = [allResidualsDecr; cycle(i, cycleIndex).residualsDecr];

    cycle(i, cycleIndex).timeStampBoth = [cycle(i, cycleIndex).timeStampIncr; cycle(i, cycleIndex).timeStampDecr];
    cycle(i, cycleIndex).pressureBoth = [cycle(i, cycleIndex).pressureIncr; cycle(i, cycleIndex).pressureDecr];
    VrsBoth = [cycle(i, cycleIndex).VrsIncr; cycle(i, cycleIndex).VrsDecr];
    [warnMsg, warnId] = lastwarn;
    [cycle(i, cycleIndex).polyPBoth, cycle(i, cycleIndex).polySBoth, mu] = polyfit(cycle(i, cycleIndex).pressureBoth, VrsBoth, polyOrder);
    [warnMsg, warnId] = lastwarn;

    [cycle(i, cycleIndex).VrsFitBoth, cycle(i, cycleIndex).polyDeltaBoth] = polyval(cycle(i, cycleIndex).polyPBoth, cycle(i, cycleIndex).pressureBoth, cycle(i, cycleIndex).polySBoth, mu);
    cycle(i, cycleIndex).residualsBoth = cycle(i, cycleIndex).VrsFitBoth - VrsBoth;
    cycle(i, cycleIndex).maxResidualBoth = max(abs(cycle(i, cycleIndex).residualsBoth));
    cycle(i, cycleIndex).meanResidualBoth = mean(abs(cycle(i, cycleIndex).residualsBoth));

    linPressure = 0:5:2000;
    VrsFitDecr = polyval(cycle(i, cycleIndex).polyPDecr, linPressure, cycle(i, cycleIndex).polySDecr, mu);
    VrsFitIncr = polyval(cycle(i, cycleIndex).polyPIncr, linPressure, cycle(i, cycleIndex).polySIncr, mu);
    cycle(i, cycleIndex).hysteresis = VrsFitDecr - VrsFitIncr;
    cycle(i, cycleIndex).hysteresisMean = mean(cycle(i, cycleIndex).hysteresis);
    cycle(i, cycleIndex).hysteresisMax = max(cycle(i, cycleIndex).hysteresis);

    cycleIndex = cycleIndex + 1;
end

figure
cmap = colormap(jet(500));

hold on
for k=1:width(cycle)
    x = [cycle(i,k).timeStampIncr; cycle(i,k).timeStampDecr]';
    y = [cycle(i,k).pressureIncr; cycle(i,k).pressureDecr]' ;
    z = [cycle(i,k).VrsIncr; cycle(i,k).VrsDecr]';

    colorIndex = min(500, round(20 * cycle(i,k).tempDegCUp(1)));
    plot3(x, y, z, 'color', cmap(colorIndex,:));
end

set(gca, 'ydir', 'reverse');
zLimits = zlim;
zlabel('Vrs (volts)');
ylabel('Pressure (dBar)');
set(gca, 'ztick', [zLimits(1):0.004:zLimits(2)+.001], 'FontSize', 16);
view(3);
datetick('x', 'mmmdd HH:MM')
set(gca, 'xtick', [data.LabviewDateTime(1):1:data.LabviewDateTime(end)], 'FontSize', 16);
xlim([data.LabviewDateTime(1) data.LabviewDateTime(end)]);
grid on
colorbar('Ticks', [0, 1], 'TickLabels', {'0 degC', '50 degC'});
title(['SURpH Sensor V2 (pre-loaded at 2.0 in-lbf' newline 'Vrs (volts) vs Time and Pressure (dBar)'], 'FontSize', 18)
view(3)

lineNumber = 1;
colors = ['r', 'g', 'r', 'g', 'r', 'g', 'r', 'g', 'r', 'g', 'r', 'g', 'b', 'k', 'b', 'k', 'b', 'k', 'b', 'k', 'b', 'k', 'b', 'k', 'b', 'k'];
figure
hold on
for k=1:width(cycle)
    if( (cycle(i,k).tempDegCDown(1) > 14) && (cycle(i,k).tempDegCDown(1) < 16))
        x = cycle(i,k).timeStampIncr;
        y = cycle(i,k).pressureIncr;
        z = cycle(i,k).VrsIncr;
        plot3(x, y, z, '+-', 'color', colors(lineNumber));
        z = cycle(i,k).VrsFitIncr;
        plot3(x, y, z, 'c');
        lineNumber = lineNumber + 1;

        x = cycle(i,k).timeStampDecr;
        y = cycle(i,k).pressureDecr;
        z = cycle(i,k).VrsDecr;
        plot3(x, y, z, '+-', 'color', colors(lineNumber));
        z = cycle(i,k).VrsFitDecr;
        plot3(x, y, z, 'c');
        lineNumber = lineNumber + 1;
   
        x = cycle(i,k).timeStampBoth;
        y = cycle(i,k).pressureBoth;
        z = cycle(i,k).VrsFitBoth;
        plot3(x, y, z, 'm')

    end
end
set(gca, 'ydir', 'reverse');
zLimits = zlim;
zlabel('Vrs (volts)');
ylabel('Pressure (dBar)');
set(gca, 'ztick', [zLimits(1):0.004:zLimits(2)+.001], 'FontSize', 16);
view(3);
datetick('x', 'mmmdd HH:MM')
set(gca, 'xtick', [data.LabviewDateTime(1):1:data.LabviewDateTime(end)], 'FontSize', 16);
xlim([data.LabviewDateTime(1) data.LabviewDateTime(end)]);
grid on
legend('first set P incr', 'first set P incr fit', 'first set P decr', 'first set P decr fit', ...
    'second set P incr', 'second set P incr fit', 'second set P decr', 'second set P decr fit', 'first set both fit', 'second set both fit');

title(['SURpH Sensor V2 (pre-loaded at 2.0 in-lbf' newline 'Vrs (volts) vs Time and Pressure (dBar) 15 degC'], 'FontSize', 18)
view(3)

sizeCycle = size(cycle);
numFiles = sizeCycle(1);
numCycles = sizeCycle(2);

fitDecr = figure;
figure(fitDecr)
hold on
title('Vrs Fit (volts) Decr vs. Pressure (dBar)   fit order: ' + string(polyOrder))

residualsDecr = figure;
figure(residualsDecr)
hold on
title('Residuals (volts) Decr vs. Pressure (dBar)   fit order: ' + string(polyOrder))

residualsBoth = figure;
figure(residualsBoth)
hold on
title('Residuals (volts) Both vs. Pressure (dBar)   fit order: ' + string(polyOrder))

hysteresis = figure;
figure(hysteresis)
hold on
title('Hysteresis (volts) Both vs. Pressure (dBar)   fit order: ' + string(polyOrder))




for j = 1:length(cycle)
    if(length(cycle(i,j).residualsIncr > 0))
        figure(fitDecr)
        plot(cycle(i,j).pressureDecr, cycle(i,j).VrsFitDecr)
        figure(residualsDecr)
        plot(cycle(i,j).pressureDecr, cycle(i,j).residualsDecr)
        figure(residualsBoth)
        plot(cycle(i,j).pressureBoth, cycle(i,j).residualsBoth)
    else
        disp("skipping resiudal statistics calc for zero length array i, j: " + i + j)
    end
end

disp("Starting hysteresis statistics calc")
allIndex = 1;
sizeCycle = size(cycle);
numFiles = sizeCycle(1);
numCycles = sizeCycle(2);
for i = 1:numFiles
    for j = 1:numCycles
        if(length(cycle(i,j).maxResidualIncr > 0))
            allMaxResidualIncr(allIndex) = cycle(i,j).maxResidualIncr;
            allMaxResidualDecr(allIndex) = cycle(i,j).maxResidualDecr;
            allMaxResidualBoth(allIndex) = cycle(i,j).maxResidualBoth;
            allMeanResidualIncr(allIndex) = cycle(i,j).meanResidualIncr;
            allMeanResidualDecr(allIndex) = cycle(i,j).meanResidualDecr;
            allMeanResidualBoth(allIndex) = cycle(i,j).meanResidualBoth;
            allIndex = allIndex + 1;
        else
            disp("skipping resiudal statistics calc for zero length array i, j: " + i + j)
        end
    end
end

figure
histogram(allMaxResidualIncr, 100)
title('MaxResidualsIncr    fit order: ' + string(polyOrder))
figure
histogram(allMaxResidualDecr, 100)
title('MaxResidualsDecr    fit order: ' + string(polyOrder))
figure
histogram(allMaxResidualBoth, 100)
title('MaxResidualsBoth    fit order: ' + string(polyOrder))

figure
histogram(allMeanResidualIncr, 100)
title('MeanResidualsIncr    fit order: ' + string(polyOrder))
figure
histogram(allMeanResidualDecr, 100)
title('MeanResidualsDecr    fit order: ' + string(polyOrder))
figure
histogram(allMeanResidualBoth, 100)
title('MeanResidualsBoth    fit order: ' + string(polyOrder))