function [f,R] = HysteresisCalcV3(T_many,FN_list)
% HyseresisCalcV3
%
% Function that grabs data from the ovm_TEST.m script and
% computes several histograms to evaluate the hysteresis of the "One" data
% set from the ovm_TEST.m scipt against the "Many"
%
% INPUT ARGUMENTS:
%
% T_Many - a table that contains the raw data output from the PTC labview
% program and uses only the Vrs, Pressure, and Datetime to determine the
% Vrs voltage signal against the calculated fit and parses everything by
% cycle 
%
% FN_list - the list of folder names of calibration runs that the "Many"
% data is compiled from
%
% OUTPUTS:
%
% f - struct of figures of the histogram
%
% R - struct of the hysteresis calculations that can be used to evaluate
% one or more calibration runs to create a histogram comparison
%
% Ben Davis
% 11/30/2023
%
% Derived from CalcHystersisV3.m by Gene Massion 2023
%
% Most recent update made on 11/30/2023 by BD

%% Parsing and Calculations

% sortedDataFiles = M.list_FN;
sortedDataFiles = FN_list;

polyOrder = 6;
numValuesIndex = 1;
for i = 1: length(T_many)
    disp("File Num: " + i + " Processing: " + char(sortedDataFiles(i)));

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


    pressure = T_many(i).g.pres(indexes);
    Vrs = T_many(i).g.Vrs(indexes);
    tempDegC = T_many(i).g.TC(indexes);
    timeStamp = T_many(i).g.DT(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)
        % Tossing out first samples
%         disp("tossing first samples");
        pressure = pressure(indexBreaks(2)+1:end);
        Vrs = Vrs(indexBreaks(2)+1:end);
        timeStamp = timeStamp(indexBreaks(2)+1:end);
        indexBreaks = find(dP3 ~= 0);
        indexBreaks = indexBreaks + 1;
        indexBreaks = [1; indexBreaks];
    end

    % Normalizing runs to merge the exceptions that have odd parsing from
    % the Labview program
    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).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 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).tempDegC = tempDegC(indexBreaks(end):length(pressure));
    ramp(j+1).timeStamp = timeStamp(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)

% DOWN IS DECREASING PRESSURE AND UP IN INCREASING PRESSURE %
        cycle(i, cycleIndex).pressureUp = ramp(j).pressure;
        cycle(i, cycleIndex).VrsUp = ramp(j).Vrs;
        cycle(i, cycleIndex).timeStampUp = ramp(j).timeStamp;

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

        lastwarn('')
        [cycle(i, cycleIndex).polyPUp, cycle(i, cycleIndex).polySUp, mu] = 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, mu);
        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, mu] = 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, mu);
        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, mu] = 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, 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).polyPDown, linPressure, cycle(i, cycleIndex).polySDown, mu);
        VrsFitIncr = polyval(cycle(i, cycleIndex).polyPUp, linPressure, cycle(i, cycleIndex).polySUp, mu);
        cycle(i, cycleIndex).hysteresis = abs(VrsFitDecr - VrsFitIncr);
        cycle(i, cycleIndex).hysteresisMean = mean(cycle(i, cycleIndex).hysteresis);
        cycle(i, cycleIndex).hysteresisMax = max(cycle(i, cycleIndex).hysteresis);

        cycleIndex = cycleIndex + 1;
    end

% Commenting out next few lines so plots do not pop up
    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]';
    end
 end

% Hysteresis Calculations for the number of tables in the T_many input
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;
            allMeanHysteresis(allIndex) = cycle(i,j).hysteresisMean;
            allMaxHysteresis(allIndex) = cycle(i,j).hysteresisMax;
            allIndex = allIndex + 1;
        else
%             disp("skipping resiudal statistics calc for zero length array i, j: " + i + j);
        end
    end
end

% Creating the edges for histogram for the "many" hysteresis data
% edges_mean = 0:0.00001:0.00032; edges_mean = [edges_mean 0.00038];
edges_mean = linspace(0,(max(allMeanResidualDown)+std(allMeanResidualDown)),16); 
edges_max = 0:0.0001:0.0023; edges_max = [edges_max 0.004];
R = struct();
R.allMaxResidualUp = allMaxResidualUp;
R.allMaxResidualDown = allMaxResidualDown;
R.allMaxResidualBoth = allMaxResidualBoth;

R.allMeanResidualUp = allMeanResidualUp;
R.allMeanResidualDown = allMeanResidualDown;
R.allMeanResidualBoth = allMeanResidualBoth;

R.edges_mean = edges_mean;
R.edges_max = edges_max;

R.std_max_up = std(allMaxResidualUp);
R.std_max_down = std(allMaxResidualDown);
R.std_max_both = std(allMaxResidualBoth);
R.mean_max_both = mean(allMaxResidualBoth);
R.mean_max_down = mean(allMaxResidualDown);
R.med_max_both = median(allMaxResidualBoth);
R.med_max_down = median(allMaxResidualDown);

R.std_mean_up = std(allMeanResidualUp);
R.std_mean_down = std(allMeanResidualDown);
R.std_mean_both = std(allMeanResidualBoth);
R.mean_mean_both = mean(allMeanResidualBoth);
R.mean_mean_down = mean(allMeanResidualDown);
R.med_mean_both = median(allMeanResidualBoth);
R.med_mean_down = median(allMeanResidualDown);

R.hysteresisMean = median(abs(allMeanHysteresis));
R.hysteresisMax = median(abs(allMaxHysteresis));

%% Plotting Hysteresis

% 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).residualsUp > 0))
%         
%         figure(hysteresis)
%         plot(linPressure, cycle(i,j).hysteresis)
%     else
%         disp("skipping resiudal statistics calc for zero length array i, j: " + i + j)
%     end
% end

%% Plotting histograms

f1 = figure('Name','Max Vrs Fit Residuals','unit', 'inch', 'position', [0.5 0.5 10.5 6.8]);
histogram(allMaxResidualUp, edges_max); hold on;
histogram(allMaxResidualDown, edges_max); hold on;
histogram(allMaxResidualBoth, edges_max);
xline(R.mean_max_both + R.std_max_both,'-.','Std Dev of Max Vrs Fit Both','LineWidth',2.5);
grid on;
legend('Max Vrs fit residual Pressure increasing', 'Max Vrs fit residual Pressure decreasing', 'Max Vrs fit residual Pressure both')
xlabel('Residual from Fit (mV)')
ylabel('Number of Cycles')
f.ax = gca; f.ax.Color = [0.9 0.9 0.9];
f.ax.FontSize = 16; f.ax.FontWeight = 'Bold';
f.F1  = ancestor(f1, 'figure');
close(f1);

f2 = figure('Name','Mean Vrs Fit Residuals','unit', 'inch', 'position', [0.5 0.5 10.5 6.8]);
histogram(allMeanResidualUp, edges_mean); hold on;
histogram(allMeanResidualDown, edges_mean); hold on;
histogram(allMeanResidualBoth, edges_mean);
xline(R.mean_mean_both + R.std_mean_both,'-.','Std Dev of Mean Vrs Fit Both','LineWidth',2.5);
grid on;
legend('Mean Vrs fit residual Pressure increasing', 'Mean Vrs fit residual Pressure decreasing', 'Mean Vrs fit residual Pressure both')
xlabel('Residual from Fit (mV)')
ylabel('Number of Cycles')
f.ax = gca; f.ax.Color = [0.9 0.9 0.9];
f.ax.FontSize = 16; f.ax.FontWeight = 'Bold';
f.F2  = ancestor(f2, 'figure');
close(f2);

f3 = figure('Name','Mean Vrs Fit Residuals','unit', 'inch', 'position', [0.5 0.5 10.5 6.8]);
yyaxis left
histogram(allMeanResidualBoth, edges_mean); hold on;
xline(R.mean_mean_both + R.std_mean_both,'-.','Std Dev of Mean Vrs Fit Both','LineWidth',2.5);
ylabel('Number of Cycles')
yyaxis right
histogram(allMeanResidualDown, edges_mean);
grid on;
legend('Mean Vrs fit residual Pressure both','','Mean Vrs fit residual Pressure decreasing')
xlabel('Residual from Fit (mV)')
ylabel('Number of Cycles')
f.ax = gca; f.ax.Color = [0.9 0.9 0.9];
f.ax.FontSize = 16; f.ax.FontWeight = 'Bold';
f.F3  = ancestor(f3, 'figure');
close(f3);

return;
end

