clear all
close all

%% Import data from text file
% Script for importing data from the following text file:
%
%    filename: C:\Users\magene\Documents\GDF\Data\AllDataFromPerformances2023Nov13\subsetAE\AE GDF 2000 Ti cvrplt_012623.csv
%
% Auto-generated by MATLAB on 15-Nov-2023 10:53:08

%% Set up the Import Options and import the data
opts = delimitedTextImportOptions("NumVariables", 36);

% 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", "IBAMPSSet_temp", "Set_pres", "StepNumH", "DataloggerSerialNum", "SN0777H", "ApplicationName", "APF_MAINC"];
opts.VariableTypes = ["datetime", "double", "double", "double", "double", "double", "double", "double", "string", "categorical", "string", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "string", "string", "string"];

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

% Specify variable properties
% opts = setvaropts(opts, ["CRC", "SN0777H", "ApplicationName", "APF_MAINC"], "WhitespaceRule", "preserve");
% opts = setvaropts(opts, ["CRC", "ID", "SN0777H", "ApplicationName", "APF_MAINC"], "EmptyFieldRule", "auto");
%opts = setvaropts(opts, "LabviewDateTime", "InputFormat", "dd/MM/yyyy hh:mm:ss aa");
%opts = setvaropts(opts, "DATETIME", "InputFormat", "dd/MM/yyyy HH:mm:ss");

path = 'C:\Users\magene\Documents\GDF\Data\AllDataFromPerformances2023Nov13\subsetFinalCalData\'
%path = 'C:\Users\magene\Documents\GDF\Data\GDFLionixV2\'
fullPath = strcat(path, '*.csv');

dataFiles = dir(strcat(path, '**\*.csv'));
% for i=1:length(dataFiles)
%     dataFiles(i).name = strrep(dataFiles(i).name, ' ', '_')
% end

fileDatenums = dataFiles.date;
[sortedTime, i] = sort([dataFiles.datenum]);
sortedDataFiles = dataFiles(i);

polyOrder = 3;
numValuesIndex = 1;
for i = 1: length(dataFiles)
    disp("File Num: " + i + " Processing: " + sortedDataFiles(i).name);
    clearvars -except opts path cycle sortedDataFiles polyOrder numValuesIndex i
    data = readtable(strcat(sortedDataFiles(i).folder, "\", sortedDataFiles(i).name), opts);

    %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);
    timeStamp = data.LabviewDateTime(indexes);
    tempDegC = data.ThermistorTemp(indexes);

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

    indexBreaks = find(dP3 ~= 0);
    indexBreaks = indexBreaks + 1;
    indexBreaks = [1; indexBreaks];
    if( (indexBreaks(2) - indexBreaks(1)) < 10)
        disp("tossing first samples")
        pressure = pressure(indexBreaks(2)+1:end);
        Vrs = Vrs(indexBreaks(2)+1:end);
        tempDegC = TempDegC(indexBreaks(2)+1:end);
        timeStamp = timeStamp(indexBreaks(2)+1:end);
        indexBreaks = find(dP3 ~= 0);
        indexBreaks = indexBreaks + 1;
        indexBreaks = [1; indexBreaks];
    end

    if(i == 41)
        lastBreak = 15;
    else
        lastBreak = length(indexBreaks)-1;
    end

    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).timeStamp = timeStamp(indexBreaks(j):indexBreaks(j+1)-1);
        ramp(j).tempDegC = tempDegC(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 down 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).timeStamp = timeStamp(indexBreaks(end):length(pressure));
    ramp(j+1).tempDegC = tempDegC(indexBreaks(end):length(pressure));
    numValues(numValuesIndex) = length(ramp(j+1).pressure);

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

        cycle(i, cycleIndex).pressureDown = ramp(j+1).pressure;
        cycle(i, cycleIndex).VrsDown = ramp(j+1).Vrs;
        cycle(i, cycleIndex).timeStampDown = ramp(j+1).timeStamp;
        cycle(i, cycleIndex).tempDegCDown = ramp(j+1).tempDegC;

        lastwarn('')
        [cycle(i, cycleIndex).polyPUp, cycle(i, cycleIndex).polySUp] = polyfit(cycle(i, cycleIndex).pressureUp, cycle(i, cycleIndex).VrsUp, polyOrder);
        [warnMsg, warnId] = lastwarn;

        [cycle(i, cycleIndex).VrsFitUp, cycle(i, cycleIndex).polyDeltaUp] = polyval(cycle(i, cycleIndex).polyPUp, cycle(i, cycleIndex).pressureUp, cycle(i, cycleIndex).polySUp);
        cycle(i, cycleIndex).residualsUp = cycle(i, cycleIndex).VrsFitUp - cycle(i, cycleIndex).VrsUp;
        cycle(i, cycleIndex).maxResidualUp = max(abs(cycle(i, cycleIndex).residualsUp));
        cycle(i, cycleIndex).meanResidualUp = mean(abs(cycle(i, cycleIndex).residualsUp));

        lastwarn('')
        [cycle(i, cycleIndex).polyPDown, cycle(i, cycleIndex).polySDown] = polyfit(cycle(i, cycleIndex).pressureDown, cycle(i, cycleIndex).VrsDown, polyOrder);
        [warnMsg, warnId] = lastwarn;

        [cycle(i, cycleIndex).VrsFitDown, cycle(i, cycleIndex).polyDeltaDown] = polyval(cycle(i, cycleIndex).polyPDown, cycle(i, cycleIndex).pressureDown, cycle(i, cycleIndex).polySDown);
        cycle(i, cycleIndex).residualsDown = cycle(i, cycleIndex).VrsFitDown - cycle(i, cycleIndex).VrsDown;
        cycle(i, cycleIndex).maxResidualDown = max(abs(cycle(i, cycleIndex).residualsDown));
        cycle(i, cycleIndex).meanResidualDown = mean(abs(cycle(i, cycleIndex).residualsDown));

        pressureBoth = [cycle(i, cycleIndex).pressureUp; cycle(i, cycleIndex).pressureDown];
        VrsBoth = [cycle(i, cycleIndex).VrsUp; cycle(i, cycleIndex).VrsDown];
        [warnMsg, warnId] = lastwarn;
        [cycle(i, cycleIndex).polyPBoth, cycle(i, cycleIndex).polySBoth] = polyfit(pressureBoth, VrsBoth, polyOrder);
        [warnMsg, warnId] = lastwarn;

        [cycle(i, cycleIndex).VrsFitBoth, cycle(i, cycleIndex).polyDeltaBoth] = polyval(cycle(i, cycleIndex).polyPBoth, pressureBoth, cycle(i, cycleIndex).polySBoth);
        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));

        cycleIndex = cycleIndex + 1;
    end

    figure
    cmap = colormap(jet(500));

    hold on
    for k=1:width(cycle)
        x = [cycle(i,k).timeStampUp; cycle(i,k).timeStampDown]';
        y = [cycle(i,k).pressureUp; cycle(i,k).pressureDown]' ;
        z = [cycle(i,k).VrsUp; cycle(i,k).VrsDown]';

        %plot3(x, y, z);
        if(length(x) > 0)
            plot3(x, y, z, 'Color', cmap(   round(10 * cycle(i,k).tempDegCUp(1)),:));
        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;
    title(" File " + i + ": " + sortedDataFiles(i).name)
    view(3)
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
        size(cycle)
        if(length(cycle(i,j).maxResidualUp > 0))
            allMaxResidualUp(allIndex) = cycle(i,j).maxResidualUp;
            allMaxResidualDown(allIndex) = cycle(i,j).maxResidualDown;
            allMaxResidualBoth(allIndex) = cycle(i,j).maxResidualBoth;
            allMeanResidualUp(allIndex) = cycle(i,j).meanResidualUp;
            allMeanResidualDown(allIndex) = cycle(i,j).meanResidualDown;
            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
hold on
histogram(allMaxResidualUp, 100)
histogram(allMaxResidualDown, 100)
histogram(allMaxResidualBoth, 100)
help legend
legend('Max Vrs fit residual Pressure increasing', 'Max Vrs fit residual Pressure decreasing', 'Max Vrs fit residual Pressure both')

figure
hold on
histogram(allMeanResidualUp, 100)
histogram(allMeanResidualDown, 100)
histogram(allMeanResidualBoth, 100)
help legend
legend('Mean Vrs fit residual Pressure increasing', 'Mean Vrs fit residual Pressure decreasing', 'Mean Vrs fit residual Pressure both')




