%% INPUT THE DATA


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\ProcessThese\';
fullPath = strcat(path, '*.csv');

[file, location] = uigetfile(fullPath);

currentFile = strcat(location, file);
disp("Reading: " + currentFile);
dataTable = readtable(currentFile, opts);
disp('Read data done')

%
%% Separate out cycles

[cycle, numCycles] = getCyclesFromPressure3(dataTable, i);
disp('Separate out cycle Done')

%
%% Plot pressure and temperature vs time

figure
plot(dataTable.dtLV, dataTable.pressureDBar)
ylabel('Pressure (dBar)')
yyaxis right
plot(dataTable.dtLV, dataTable.thermTemp)
ylabel("Temperature (degC)")
title('Pressure (dBar) and Temperature (degC) vs Time')

ax = gca;
ax.XAxis.MinorTick = 'on';
ax.XAxis.TickValues = dataTable.dtLV(1):hours(12):dataTable.dtLV(end);
ax.XAxis.MinorTickValues = dataTable.dtLV(1):hours(2):dataTable.dtLV(end);
grid on

%
%% Plot Vrs vs Time and Waterfall,

 
titleString = "Vrs vs Time " + file;
figure('name', titleString, 'NumberTitle', 'off');
hold on
yyaxis right
plot(dataTable.dtLV, dataTable.pressureDBar)
yyaxis left
scatter(dataTable.dtLV, dataTable.vrsVolts, 10, dataTable.thermTemp, 'filled')
colormap(turbo);                                   
colorbar;                                    
clim([min(dataTable.thermTemp) max(dataTable.thermTemp)]);   
ax = gca;
ax.XAxis.MinorTick = 'on';
ax.XAxis.TickValues = dataTable.dtLV(1):hours(12):dataTable.dtLV(end);
ax.XAxis.MinorTickValues = dataTable.dtLV(1):hours(2):dataTable.dtLV(end);
grid on
grid minor
title(titleString)
hold off

% crossCorr = xcorr(dataTable.pressureDBar, dataTable.vrsVolts, 'normalized');
% figure
% plot(crossCorr)
% title('Cross Correlation of pressure and Vrs')

% Waterfall plot
plotPressureWaterfall(cycle, numCycles)
title(titleString)
disp('Vrs vs Time and Waterfall plots done')
xlabel('Time');
ylabel('Pressure (dBar)');
zlabel('Vrs (volts)');

%
%% Separate cycle struct into lowPHold, pIncr, highPHold and pDecr structs

lowPHold = [];
pIncr= [];
highPHold = [];
pDecr = [];

lowPHoldIndex = 1;
pIncrIndex = 1;
highPHoldIndex = 1;
pDecrIndex = 1;

for i=1:length(cycle)
    if(cycle(i).cycleName == 'lowPHold')
        lowPHold(lowPHoldIndex).pressure = cycle(i).pressureDBar;
        lowPHold(lowPHoldIndex).temperature = cycle(i).thermTemp;
        lowPHold(lowPHoldIndex).Vrs = cycle(i).vrsVolts;
        lowPHold(lowPHoldIndex).cycleNum = i;
        lowPhold(lowPHoldIndex).datetime = cycle(i).dt;
        lowPHoldIndex = lowPHoldIndex + 1;
    end
end

for i=1:length(cycle)
    if(cycle(i).cycleName == 'pIncr')
        pIncr(pIncrIndex).pressure = cycle(i).pressureDBar;
        pIncr(pIncrIndex).temperature = cycle(i).thermTemp;
        pIncr(pIncrIndex).Vrs = cycle(i).vrsVolts;
        pIncr(pIncrIndex).cycleNum = i;
        pIncr(pIncrIndex).datetime = cycle(i).dt;
        pIncrIndex = pIncrIndex + 1;
    end
end

for i=1:length(cycle)
    if(cycle(i).cycleName == 'highPHold')
        highPHold(highPHoldIndex).pressure = cycle(i).pressureDBar;
        highPHold(highPHoldIndex).temperature = cycle(i).thermTemp;
        highPHold(highPHoldIndex).Vrs = cycle(i).vrsVolts;
        highPHold(highPHoldIndex).cycleNum = i;
        highPHold(highPHoldIndex).datetime = cycle(i).dt;
        highPHoldIndex = highPHoldIndex + 1;
    end
end

for i=1:length(cycle)
    if(cycle(i).cycleName == 'pDecr')
        pDecr(pDecrIndex).pressure = cycle(i).pressureDBar;
        pDecr(pDecrIndex).temperature = cycle(i).thermTemp;
        pDecr(pDecrIndex).Vrs = cycle(i).vrsVolts;
        pDecr(pDecrIndex).cycleNum = i;
        pDecr(pDecrIndex).datetime = cycle(i).dt;
        pDecrIndex = pDecrIndex + 1;
    end
end

%
%% Plot the decreasing pressure data


pressurePDecr = [];
temperaturePDecr = [];
VrsPDecr = [];

for i=1:length(cycle)
    if(cycle(i).cycleName == 'pDecr')
        pressurePDecr = [pressurePDecr cycle(i).pressureDBar];
        temperaturePDecr = [temperaturePDecr cycle(i).thermTemp];
        VrsPDecr = [VrsPDecr cycle(i).vrsVolts];
    end
end

figure
plot(pressurePDecr, temperaturePDecr, '+')
xlabel('Pressure (DBar)');
ylabel('Temperature (°C)');
title('Pressure Decreasing vs. Temperature');

plot(pressurePDecr, VrsPDecr, '+')
xlabel('Pressure (DBar)');
ylabel('Vrs (Volts)');
title('Pressure Decreasing vs. Vrs');

%
%% Calculate a polynomial fit for independent variables: all pressure decreasing and
%% all temperature dependent variable: all Vrs


pressure = [];
temperature = [];
Vrs = [];

for i=1:length(pDecr)
    pressure = [pressure pDecr(i).pressure];
    temperature = [temperature pDecr(i).temperature];
    Vrs = [Vrs pDecr(i).Vrs];
end

x = [pressure; temperature]';
[vrs2DFit, gof, output] = fit(x, Vrs', 'poly32', 'Normalize', 'on')
plot(vrs2DFit, x, Vrs')

residualMean = mean(output.residuals);
residualStdDev = std(output.residuals);
VrsRange = max(Vrs) - min(Vrs);

xlabel('Pressure (dBar)')
ylabel('Temperature (degC)')
zlabel('Vrs (volts)')
titleString = sprintf('Vrs vs Pressure and Temperature\n f(p) (mV): %6.2f  residual mean (mV): %5.3f   residual stdDev (mV): %5.3f', VrsRange * 1000, residualMean*1000, residualStdDev*1000)
title(titleString)

figure
plot(pressure, output.residuals*1e6, ".")
xlabel('Pressure (dBar)')
ylabel('Residuals (uVolts)')
title('Residual vs P decr')
figure
plot(temperature, output.residuals*1e6, ".")
xlabel('Temperature (degC)')
ylabel('Residuals (uVolts)')
title('Residual vs Temperatue (P decr)')

condMatrix = [pressurePDecr, temperaturePDecr, VrsPDecr];

%

Hd = PressureDerivativeFilter();
pressureDeriv = filtfilt(Hd.Numerator, 1, pressurePDecr(348:612));
figure
hold on
plot(pressureDeriv)
plot(diff(pressurePDecr(348:612)))
grid
yyaxis right
plot(pressurePDecr(348:612))
legend('Pdecr Deriv Filter', 'Pdecr Diff', 'Pdecr')
title 'Various calculated pressure derivatives'
% hold off

%
%% Calculate polyfit for decreasing pressure ramps


plotIndivFits = 0;

allData.pressure = [];
allData.temperature = [];
allData.Vrs = [];

clear ramp
polyOrder = 4;
for i=1:length(pDecr)
    [pDecr(i).fitCoef, pDecr(i).sPDecr, pDecr(i).mu] = polyfit(pDecr(i).pressure, pDecr(i).Vrs, polyOrder);
    [pDecr(i).vrsFitValues, pDecr(i).delta] = polyval(pDecr(i).fitCoef, pDecr(i).pressure, pDecr(i).sPDecr, pDecr(i).mu);
    pDecr(i).vrsResiduals = pDecr(i).Vrs - pDecr(i).vrsFitValues;
    [S, M] = std(pDecr(i).vrsResiduals);
    pDecr(i).vrsResidualsStdDev = S;
    pDecr(i).vrsResidualsMean = M;
    pDecr(i).VrsRange = max(pDecr(i).Vrs) - min(pDecr(i).Vrs);
    pDecr(i).meanTemp = mean(pDecr(i).temperature);

    if(plotIndivFits)
        figure
        plot(pDecr(i).pressure, pDecr(i).Vrs, '+-')
        hold on
        plot(pDecr(i).pressure, pDecr(i).vrsFitValues, 'or-')
        xlabel('Pressure (dBar)')
        ylabel('Vrs (volts)')
        yyaxis right
        plot(pDecr(i).pressure, pDecr(i).vrsResiduals*1000, 'x-')
        ylabel('Residual Vrs (mV)')
        plot(pDecr(i).pressure, pDecr(i).delta*1000, 'o-')
        hold off
        titleString = sprintf('f(p) (mV) = %6.3f Res StdDev = %6.3f mV   Res Mean = %6.3f mV', pDecr(i).VrsRange, S*1000, M*1000);
        title(titleString)
        legend('Vrs (volts)', 'Vrs ', 'Resdiuals (uV)', 'Std Error')
    end
end

%
%% Plot Vrs and Polyfit Grouped by Temperature

    fig1 = figure();
    title('Vrs and Vrs Fit    Temp: <8')
    ax1 = gca;
    fig1ColorIndex = 1;
    hold on
    fig2 = figure();
    title('Vrs and Vrs Fit    Temp: 8-15')
    ax2 = gca;
    fig2ColorIndex = 1;
    hold on
    fig3 = figure();
    title('Vrs and Vrs Fit    Temp: >15')
    ax3 = gca;
    fig3ColorIndex = 1;
    hold on

    for i=1:length(pDecr)
        if( pDecr(i).meanTemp < 8)
            figure(fig1)
            ax1.ColorOrderIndex = fig1ColorIndex;
            plot(pDecr(i).pressure, pDecr(i).Vrs)
            ax1.ColorOrderIndex = fig1ColorIndex;
            plot(pDecr(i).pressure, pDecr(i).vrsFitValues, '--')
            fig1ColorIndex = fig1ColorIndex + 1;
        elseif( (pDecr(i).meanTemp >= 5) && (pDecr(i).meanTemp < 15))
            figure(fig2)
            ax2.ColorOrderIndex = fig2ColorIndex;
            plot(pDecr(i).pressure, pDecr(i).Vrs)
            hold on
            ax2.ColorOrderIndex = fig2ColorIndex;
            plot(pDecr(i).pressure, pDecr(i).vrsFitValues, '--')
            fig2ColorIndex = fig2ColorIndex + 1;
        else
            figure(fig3)
            ax3.ColorOrderIndex = fig3ColorIndex;
            plot(pDecr(i).pressure, pDecr(i).Vrs)
            hold on
            ax3.ColorOrderIndex = fig3ColorIndex;
            plot(pDecr(i).pressure, pDecr(i).vrsFitValues, '--')
            fig3ColorIndex = fig3ColorIndex + 1;
        end
    end

    figure
    histogram([pDecr.meanTemp])

%% Calculate polyfit for increasing pressure ramps


plotIndivFits = 0;

allData.pressure = [];
allData.temperature = [];
allData.Vrs = [];

clear ramp
polyOrder = 4;
for i=1:length(pIncr)
    [pIncr(i).fitCoef, pIncr(i).sPIncr, pIncr(i).mu] = polyfit(pIncr(i).pressure, pIncr(i).Vrs, polyOrder);
    [pIncr(i).vrsFitValues, pIncr(i).delta] = polyval(pIncr(i).fitCoef, pIncr(i).pressure, pIncr(i).sPIncr, pIncr(i).mu);
    pIncr(i).vrsResiduals = pIncr(i).Vrs - pIncr(i).vrsFitValues;
    [S, M] = std(pIncr(i).vrsResiduals);
    pIncr(i).vrsResidualsStdDev = S;
    pIncr(i).vrsResidualsMean = M;
    pIncr(i).VrsRange = max(pIncr(i).Vrs) - min(pIncr(i).Vrs);
    pIncr(i).meanTemp = mean(pIncr(i).temperature);

    if(plotIndivFits)
        figure
        plot(pIncr(i).pressure, pIncr(i).Vrs, '+-')
        hold on
        plot(pIncr(i).pressure, pIncr(i).vrsFitValues, 'or-')
        xlabel('Pressure (dBar)')
        ylabel('Vrs (volts)')
        yyaxis right
        plot(pIncr(i).pressure, pIncr(i).vrsResiduals*1000, 'x-')
        ylabel('Residual Vrs (mV)')
        plot(pIncr(i).pressure, pIncr(i).delta*1000, 'o-')
        hold off
        titleString = sprintf('f(p) (mV) = %6.3f Res StdDev = %6.3f mV   Res Mean = %6.3f mV', pIncr(i).VrsRange, S*1000, M*1000);
        title(titleString)
        legend('Vrs (volts)', 'Vrs ', 'Resdiuals (uV)', 'Std Error')
    end
end

%
%% Calculate hysteresis
%

pressures = 0:10:500;
pressures = [pressures 550:50:2000];

for i = 1:length(pIncr)
    [vrsPdecr, vrsPdecrDelta] = polyval(pDecr(i).fitCoef, pressures, pDecr(i).sPDecr, pDecr(i).mu);
    [vrsPincr, vrsPincrDelta] = polyval(pIncr(i).fitCoef, pressures, pIncr(i).sPIncr, pIncr(i).mu);
    hysteresis(i).hyst = vrsPincr - vrsPdecr;
    hysteresis(i).mean = mean(hysteresis(i).hyst);
    hysteresis(i).stdDev = std(hysteresis(i).hyst);
    hysteresis(i).max = max(hysteresis(i).hyst);
    hysteresis(i).min = min(hysteresis(i).hyst);

    figure
    hold on
    plot(pressures, vrsPdecr)
    plot(pressures, vrsPincr)
    yyaxis right
    plot(pressures, hysteresis(i).hyst)
    ax = gca;
    ax.YAxis(2).Exponent = -3;
    hold off
    figure
    histogram(hysteresis(i).hyst)
end

%
%% Plot all the increasing and decreasing data


lengthCMap = 64;
cMap = turbo(lengthCMap);
minTemp = 1;
maxTemp = 32.5;

figure
hold on
for i=1:length(pIncr)
    pDecr(i).meanTemp;
    colorIndex = round((2*pDecr(i).meanTemp)) - 1;
    if(colorIndex > 64)
        colorIndex = 64;
    elseif(colorIndex < 1)
            colorIndex = 1;
    end
    colorIndex = max(1, colorIndex);
    plot(pDecr(i).pressure, pDecr(i).Vrs, 'color', cMap(colorIndex,:))
    plot(pDecr(i).pressure, pDecr(i).vrsFitValues, '-k')
    plot(pIncr(i).pressure, pIncr(i).Vrs, 'color', cMap(colorIndex,:), LineWidth=3)
    plot(pIncr(i).pressure, pIncr(i).vrsFitValues, '-k')
    xlabel('Pressure (dBar)')
    ylabel('Vrs (volts)')
end
colorbar
clim([minTemp maxTemp]);
title('All Decr and Incr Data and PolyFits')

%
%% Calculate pressurization rate


dP = diff(dataTable.pressureDBar);
dT = diff(dataTable.dateTime);
dTSec = seconds(dT);
pRate = dP./dTSec;
figure
plot(pRate)
ylabel('Press. Rate (dBar/sec)')
hold on
yyaxis right
plot(dataTable.flowRate)
ylabel('Flow Rate')
legend('Press. Rate', 'Flow Rate')
title('Press. Rate and Flow Rate vs. Time')

%
%% Plot mean hysteresis data
%

figure
plot([hysteresis.mean])
ax = gca;
ax.YAxis.Exponent = -3;
figTitle = sprintf('Hysteresis mean vs Pressure cycle   File:\r\n%s\r\n', string(currentFile));
title(figTitle, 'Interpreter', 'none')
ylabel('Hysteresis mean (volts)')
xlabel('Cycle num')
grid

%
%% Plot residual data
%

figure
ax = gca;
ax.YAxis.Exponent = -3;
hold on
for i=1:length(pDecr)
    plot([pDecr(i).vrsResiduals])
end

figTitle = sprintf('Residuals Pressure Decr vs Cycle Num   File:\r\n%s\r\n', string(currentFile));
title(figTitle, 'Interpreter', 'none')
ylabel('Residuals mean (volts)')
xlabel('Cycle num')

figure
ax = gca;
ax.YAxis.Exponent = -3;
hold on
for i=1:length(pIncr)
    plot([pIncr(i).vrsResiduals])
end

figTitle = sprintf('Residuals Pressure Incr vs Cycle Num   File:\r\n%s\r\n', string(currentFile));
title(figTitle, 'Interpreter', 'none')
ylabel('Residuals mean (volts)')
xlabel('Cycle num')
grid

%

fprintf(1,'done')