%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)');

%

%Plot the increasing pressure data

pressurePIncr = [];
temperaturePIncr = [];
VrsPIncr = [];

for i=1:length(cycle)
    if(cycle(i).cycleName == 'pIncr')
        pressurePIncr = [pressurePIncr cycle(i).pressureDBar];
        temperaturePIncr = [temperaturePIncr cycle(i).thermTemp];
        VrsPIncr = [VrsPIncr cycle(i).vrsVolts];
    end
end

figure
plot(pressurePIncr, temperaturePIncr, '+')
xlabel('Pressure (DBar)');
ylabel('Temperature (°C)');
title('Pressure Increasing vs. Temperature');

plot(pressurePIncr, VrsPIncr, '+')
xlabel('Pressure Increasing (DBar)');
ylabel('Vrs (Volts)');
title('Pressure Increasing vs. Vrs');

%

%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: pressure decreasing and
%temperature dependent variable: Vrs. All temperatures

x = [pressurePDecr; temperaturePDecr]';
[vrs2DFit, gof, output] = fit(x, VrsPDecr', 'poly32', 'Normalize', 'on')
plot(vrs2DFit, x, VrsPDecr')

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

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(pressurePDecr, output.residuals*1e6, ".")
xlabel('Pressure (dBar)')
ylabel('Residuals (uVolts)')
title('Residual vs P decr')
figure
plot(temperaturePDecr, output.residuals*1e6, ".")
xlabel('Temperature (degC)')
ylabel('Residuals (uVolts)')
title('Residual vs Temperatue (P decr)')

condMatrix = [pressurePDecr, temperaturePDecr, VrsPDecr];

%

plot(pressurePDecr)

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')
% hold off

%

%Calculate polyfit for decreasing pressure ramps

plotIndivFits = 1;

allData.pressurePDecr = [];
allData.temperaturePDecr = [];
allData.VrsPDecr = [];

clear ramp
polyOrder = 4;
rampIndex = 1;
for i=1:length(cycle)
    if(cycle(i).cycleName == 'pDecr')
        ramp(rampIndex).pressurePDecr = cycle(i).pressureDBar(:);
        ramp(rampIndex).temperaturePDecr = cycle(i).thermTemp(:);
        ramp(rampIndex).VrsPDecr = cycle(i).vrsVolts(:);

        [ramp(rampIndex).fitCoefPDecr, ramp(rampIndex).sPDecr, ramp(rampIndex).mu] = polyfit(ramp(rampIndex).pressurePDecr, ramp(rampIndex).VrsPDecr, polyOrder);
        [ramp(rampIndex).fitValuesPDecr, ramp(rampIndex).deltaPDecr] = polyval(ramp(rampIndex).fitCoefPDecr, ramp(rampIndex).pressurePDecr, ramp(rampIndex).sPDecr, ramp(rampIndex).mu);
        ramp(rampIndex).residuals = ramp(rampIndex).VrsPDecr - ramp(rampIndex).fitValuesPDecr;
        [S, M] = std(ramp(rampIndex).residuals);
        ramp(rampIndex).residualsStdDev = S;
        ramp(rampIndex).residualsMean = M;
        VrsRange = max(VrsPDecr) - min(VrsPDecr);
        
        if(plotIndivFits)
            figure
            plot(ramp(rampIndex).pressurePDecr, ramp(rampIndex).VrsPDecr, '+-')
            hold on
            plot(ramp(rampIndex).pressurePDecr, ramp(rampIndex).fitValuesPDecr, 'or-')
            xlabel('Pressure (dBar)')
            ylabel('Vrs (volts)')
            yyaxis right
            plot(ramp(rampIndex).pressurePDecr, ramp(rampIndex).residuals*1000, 'x-')
            ylabel('Residual Vrs (mV)')
            plot(ramp(rampIndex).pressurePDecr, ramp(rampIndex).deltaPDecr*1000, 'o-')
            hold off
            titleString = sprintf('f(p) (mV) = %6.3f Res StdDev = %6.3f mV   Res Mean = %6.3f mV', VrsRange, S*1000, M*1000);
            title(titleString)
            legend('Vrs (volts)', 'Vrs ', 'Resdiuals (uV)', 'Std Error')
        end

        allData.pressurePDecr = [allData.pressurePDecr; ramp(rampIndex).pressurePDecr(:)];
        allData.temperaturePDecr = [allData.temperaturePDecr; ramp(rampIndex).temperaturePDecr(:)];
        allData.VrsPDecr = [allData.VrsPDecr; ramp(rampIndex).VrsPDecr(:)];
        rampIndex = rampIndex + 1;
    end
end

[allData.fitCoefPDecr, allData.sPDecr, allData.mu] = polyfit(allData.pressurePDecr, allData.VrsPDecr, polyOrder);
[allData.fitValuesPDecr, allData.deltaPDecr] = polyval(allData.fitCoefPDecr, allData.pressurePDecr, allData.sPDecr, allData.mu);

figure
plot(allData.pressurePDecr, allData.VrsPDecr, '+')
hold on
plot(allData.pressurePDecr, allData.fitValuesPDecr, 'o')
hold off
title('All Vrs (volts) and All PDecr VrsFit (volts) vs. Pressure (dBar)')

%

figure
plot(dataTable.pressureDBar)
yyaxis right
plot(dataTable.flowRate)
legend('pressure', 'flow rate')

%


fprintf(1,'done')