clear all
close all
%% Import data from text file
% Script for importing data from the following text file:
%
%    filename: \\atlas\projectlibrary\ChemSensors\GeneM\GDF pH Sensor\Data\Performances\GDF2001\AE_GDF2001_122022_Ch F\AE_GDF2001_122022_Ch F.csv
%
% Auto-generated by MATLAB on 03-Jan-2023 13:19:14

%% Set up the Import Options and import the data
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\Performances\CSVFiles\';
%path = '.\';
path = '\\atlas\ProjectLibrary\ChemSensors\GeneM\GDF pH Sensor\Data\Cal Data\Calcs\DataFiles\AE GDF 2003 Ti cvrplt_012623\'
%path = 'C:\Users\gene\Documents\GDF\Data\Cal Data\Calcs\DataFiles\AE GDF 2003 Ti cvrplt_012623\';
fullPath = strcat(path, '*.csv');

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


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

    %% Plot Environmental Conditions
    figure('Units', 'Normalized', 'Position', [.01 + i*figIncr, .01 + i*figIncr, .7 .7], 'Name', sortedDataFiles(i).name, 'NumberTitle', 'off')
    subplot(3,2,1)
    yyaxis right
    plot(data{i}.dtLV, data{i}.pressureDBar)
    ylabel("Pressure (dBar)")
    yyaxis left
    plot(data{i}.dtLV, data{i}.thermTemp, data{i}.dtLV, data{i}.bathTemp, data{i}.dtLV, data{i}.housingTemp)
    ylabel("Temperature (degC)")
    legend(["Thermistor Temp" "Bath Temp" " Housing Temp" "Pressure"]);
    title("Pressure, Temperature vs Time");

    subplot(3,2,2)
    pRate{i} = (data{i}.pressureDBar(2:end) - data{i}.pressureDBar(1:end-1))./seconds(data{i}.dtLV(2:end) - data{i}.dtLV(1:end-1));
    plot(data{i}.dtLV(2:end), pRate{i})
    title("Pressurization Rate vs. Time");
    ylabel("Pressurization Rate (dBar/sec)")

    %% Plot sensor data vs time
    subplot(3,2,3)
    title("Vrs and Vrs Std Dev vs Time");
    yyaxis left
    plot(data{i}.dtLV, data{i}.vrsVolts)
    ylabel("Vrs (volts)");
    yyaxis right
    plot(data{i}.dtLV, data{i}.vrsStdDev)
    ylim([0 0.001]);
    ylabel("Vrs Std Dev (volts)");

    subplot(3,2,4)
    title("Vk and Vk Std Dev vs Time");
    yyaxis left
    plot(data{i}.dtLV, data{i}.vkVolts)
    ylabel("Vk (volts)");
    yyaxis right
    plot(data{i}.dtLV, data{i}.vkStdDev)
    ylim([0 0.001]);
    ylabel("Vk Std Dev (volts)");

    subplot(3,2,5)
    title("Ik and Ib vs Time");
    yyaxis left
    plot(data{i}.dtLV, data{i}.ikAmps)
    ylabel("Ik (amps)");
    yyaxis right
    plot(data{i}.dtLV, data{i}.ibAmps)
    ylabel("Ib (amps)");
end

tempSteps = diff(data{1}.SetTemp);
stepIndexes = find((tempSteps>1) | (tempSteps<-1));

startIndex = 1;
for i=1:length(stepIndexes)
    cycle{i}.dt = data{1}.dtLV(startIndex:stepIndexes(i)-1);
    cycle{i}.pressureDBar = data{1}.pressureDBar(startIndex:stepIndexes(i)-1);
    cycle{i}.vrsVolts = data{1}.vrsVolts (startIndex:stepIndexes(i)-1);
    cycle{i}.vrsStdDev = data{1}.vrsStdDev(startIndex:stepIndexes(i)-1);
    cycle{i}.thermTemp = data{1}.thermTemp(startIndex:stepIndexes(i)-1);
    cycle{i}.vkVolts = data{1}.vkVolts(startIndex:stepIndexes(i)-1);
    cycle{i}.vkStdDev = data{1}.vkStdDev(startIndex:stepIndexes(i)-1);
    cycle{i}.ikAmps = data{1}.ikAmps(startIndex:stepIndexes(i)-1);
    cycle{i}.ibAmps = data{1}.ibAmps(startIndex:stepIndexes(i)-1);

    startIndex = stepIndexes(i);
end

figure
plot(data{1}.dtLV, data{1}.vrsVolts)
title('time vs Vrs')

figure
plot(cycle{1}.dt, cycle{1}.thermTemp)
hold on
plot(cycle{2}.dt, cycle{2}.thermTemp)
title('time vs thermTemp cycles 1 and 2')

%calculate fit for each temp cycle
for i=1:length(stepIndexes)
    cycle{i}.fitCoef = polyfit(cycle{i}.thermTemp', cycle{i}.vrsVolts', 3);
    cycle{i}.yEstimate = polyval(cycle{i}.fitCoef, cycle{i}.thermTemp);
end

varIndex = 1;
for i=1:2:length(stepIndexes)-1
    x = linspace(min(cycle{i}.thermTemp), max(cycle{i}.thermTemp), 200);
    y1Estimate = polyval(cycle{i}.fitCoef, x);
    y2Estimate = polyval(cycle{i+1}.fitCoef, x);
    variance = y1Estimate - y2Estimate;
    maxVariance(varIndex) = max(variance);
    varIndex = varIndex + 1;
end
figure
plot(maxVariance)
title("Max variance between Temp decreasing fit and Temp increasing fit")

figure
hold on
plot(cycle{1}.thermTemp, cycle{1}.vrsVolts)
plot(cycle{2}.thermTemp, cycle{2}.vrsVolts)
plot(cycle{1}.thermTemp, cycle{1}.yEstimate)
plot(cycle{2}.thermTemp, cycle{2}.yEstimate)
title('Temperature vs Vrs and 3rd order fit cycles 1 and 2')
legend('cycle 1 Vrs', 'cycle 2 Vrs', 'cycle 1 poly fit', 'cycle 2 poly fit')
text(0.1, 0.1, ['Max variance between fits = ' num2str(maxVariance(1))], "Units", "normalized")

figure
hold on
plot(cycle{23}.thermTemp, cycle{23}.vrsVolts)
plot(cycle{24}.thermTemp, cycle{24}.vrsVolts)
plot(cycle{23}.thermTemp, cycle{23}.yEstimate)
plot(cycle{24}.thermTemp, cycle{24}.yEstimate)
title('Temperature vs Vrs and 3rd order fit cycles 23 and 24')
legend('cycle 23 Vrs', 'cycle 24 Vrs', 'cycle 23 poly fit', 'cycle 24 poly fit')
text(0.1, 0.1, ['Max variance between fits = ' num2str(maxVariance(23))], "Units", "normalized")

figure
hold on
[Vct, Vcp, Vc] = vrsHClCorr(cycle{1}.thermTemp, cycle{1}.pressureDBar);
Vc = Vc - Vc(1);
yyaxis left
plot(cycle{1}.dt, cycle{1}.vrsVolts)
plot(cycle{1}.dt, cycle{1}.vrsVolts + Vc)
yLimits =ylim;
ticks = yLimits(1):0.002:yLimits(2);
yticks(ticks);
yyaxis right
plot(cycle{1}.dt, Vc)
yLimits =ylim;
ticks = yLimits(1):0.002:yLimits(2);
yticks(ticks);
title('Vrs corrected for VrsHCl cycle 1')
legend('vrs', 'vrs+Vc', 'Vc');


% Waterfall plot
figure
hold on;
colormap(jet);
colorValue = -5;

for i=1:length(stepIndexes)
    x = datenum(cycle{i}.dt)';
    y = cycle{i}.thermTemp' ;
    z = cycle{i}.vrsVolts';

    lineColor = colorValue * ones(1, length(x));
    plot3(x, y, z);
    patch([x nan], [y nan], [z nan], [lineColor nan],'EdgeColor','interp','FaceColor','none')
    if(mod(i,2) == 0)
        colorValue = colorValue + 5;
    end
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}.dt(1));
xEnd = datenum(cycle{length(cycle)}.dt(end));
xlim([xStart xEnd]);

x = 1:length(stepIndexes);
for i=1:length(stepIndexes)
    c0(i) = cycle{i}.fitCoef(1);
    c1(i) = cycle{i}.fitCoef(2);
    c2(i) = cycle{i}.fitCoef(3);
    c3(i) = cycle{i}.fitCoef(4);
end

figure
hold on
yyaxis left
plot(x, c0)
plot(x, c1)
plot(x, c2)
ylim([-0.005 0.005]);
yyaxis right
plot(x, c3)
title('3rd order fit coefficients')
legend('coef0', 'coef1', 'coef2', 'coef3')



