% ovm_analysis_V2.m
%
% Compares the k2 pressure corrected f(P) of a particular GDF run versus
% the same parameters of a subset of GDF runs
%
% USER MUST ENTER VARIABLES TO EXTRACT
%
% Ben Davis 
% 10/26/2023
% Updated: 12/5/2025
%
clear
close all

addpath(genpath(pwd));
% ================== USER INPUTS =======================================

% Enter variable names as displayed in .mat file
% e.g. Var1 = 'variable';
Var1 = 'g';
Var2 = 'iPdecrease';
Var3 = 'iPincrease';
Var4 = 'k2_1atm';
Var5 = 'k2_2000';

% ==================================================================

%% Selecting the dataset
% Select the "One" calibration run to compare

fprintf(['Select the calibration run to compare (One).',...
        '\nMust be in Performances directory in GDF Production\n'])

% Select .mat file
[fname_one, pathname_one] = uigetfile('*.mat',...
               'Select a .mat file from the Performances directory to compare (One)',...
               '\\atlas\Chem\DuraFET\GDF Production\Performances\');
fpath_one = [pathname_one,fname_one];

% Loading variables
myVars = {Var1,Var2,Var3,Var4,Var5};
O = load(fpath_one, myVars{:});

% Getting time elapsed for "one"
te = datenum(O.g.DT(1));
te = datenum(O.g.DT) - te;

% Writing specific table using these variables

% f(P) k2 pressure corrected and corrected data for decreasing pressure in
% table T: pressure, f(P) k2 corr, temperature, f(P) normal,
T_one = [O.g.pres(O.iPdecrease), O.g.fp_k2variable(O.iPdecrease), O.g.TC(O.iPdecrease), O.g.fP(O.iPdecrease), te(O.iPdecrease)];
% Now for increasing pressure data
T_one_inc = [O.g.pres(O.iPincrease), O.g.fp_k2variable(O.iPincrease), O.g.TC(O.iPincrease), O.g.fP(O.iPincrease), te(O.iPincrease)];

% Normalizing pressure from 1 to 2000 dBar
O.pres = linspace(1,2000)';

% Decreasing Pressure
% O.p1 = polyfit(T_one(:,1),T_one(:,2),7); % seventh degree, 6th order default
O.p1 = polyfit(T_one(:,1),T_one(:,2),6); % sixth degree, 5th order (same as HysteresisCalcV3.m)
O.fP1 = polyval(O.p1,O.pres);
% Increasing Pressure
% O.p2 = polyfit(T_one_inc(:,1),T_one_inc(:,2),7); % seventh degree, 6th order default
O.p2 = polyfit(T_one_inc(:,1),T_one_inc(:,2),6); % sixth degree, 5th order (same as HysteresisCalcV3.m)
O.fP2 = polyval(O.p2,O.pres);

% Residuals of fit vs measurements per calibration run and relevant
% diagnostics
O.res.dec = (O.g.fp_k2variable(O.iPdecrease)) - polyval(O.p1,O.g.pres(O.iPdecrease));
O.res.avg = mean(abs(O.res.dec),1,'omitnan');
O.res.med = median(abs(O.res.dec),1,'omitnan');
O.res.dev = std(O.res.dec,1,'omitnan');

% % title of plot
% title_one = fname_one(1:end-4);

%% Switch Cases

many = input(['\n',...
              'Enter a number to select what',...
              '\ndataset to plot against (1-5):',...
              '\n1 --- AE Sensors',...
              '\n2 --- GDF Sensors',...
              '\n3 --- All Deployed Sensors',...
              '\n4 --- Deployed AE Sensors',...
              '\n5 --- Deployed GDF Sensors',...
              '\n6 --- Ti Cover Plates',...
              '\n---OR---',...
              '\n7 --- Select your own dataset\n',...
              '\n']);
switch many
        case 1
            ind_many = 16;
            title_many = 'AE Sensors';
            disp(title_many)
        case 2
            ind_many = 9;
            title_many = 'GDF Sensors';
            disp(title_many)
        case 3
            ind_many = 10;
            title_many = 'All Deployed Sensors';
            disp(title_many)
        case 4
            ind_many = 17;
            title_many = 'Deployed AE Sensors';
            disp(title_many)
        case 5
            ind_many  = 11;
            title_many = 'Deployed GDF Sensors';
            disp(title_many)
        case 6
            ind_many  = 18;
            title_many = 'Ti Cover Plate';
            disp(title_many)
        case 7
            disp('Select your own .mat files')
            % Select .mat file
            % Prompt the user to select multiple files
            [fileName, pathName, filterIndex] = uigetfile('*.*', 'Select files from Data directory', ...
                '\\atlas\Chem\DuraFET\GDF Production\Data\mat Files\', ...
                'MultiSelect', 'on');
            
            % Check if the user clicked 'Cancel'
            if isequal(fileName, 0)
                disp('User canceled file selection.');
                return;
            end
            
            % If only one file is selected, uigetfile returns a char instead of a cell array
            if ischar(fileName)
                % Convert to a cell array for consistency
                fileName = {fileName};
            end
            
            % Display selected file names and path
            disp('Selected files:');
            disp(fileName);
            disp('Selected path:');
            disp(pathName);
    otherwise
            disp('Select a number between 1-6 ya filthy scallywag')
end

while many < 1 || many > 7 || round(many) ~= many
    many = input(['\n',...
              'Enter a number to select what',...
              '\ndataset to plot against (1-5):',...
              '\n1 --- AE Sensors',...
              '\n2 --- GDF Sensors',...
              '\n3 --- All Deployed Sensors',...
              '\n4 --- Deployed AE Sensors',...
              '\n5 --- Deployed GDF Sensors',...
              '\n---OR---',...
              '\n6 --- Select your own dataset\n',...
              '\n']);
    switch many
        case 1
            ind_many = 16;
            title_many = 'AE Sensors';
            disp(title_many)
        case 2
            ind_many = 9;
            title_many = 'GDF Sensors';
            disp(title_many)
        case 3
            ind_many = 10;
            title_many = 'All Deployed Sensors';
            disp(title_many)
        case 4
            ind_many = 17;
            title_many = 'Deployed AE Sensors';
            disp(title_many)
        case 5
            ind_many  = 11;
            title_many = 'Deployed GDF Sensors';
            disp(title_many)
        case 6
            ind_many  = 18;
            title_many = 'Ti Cover Plate';
            disp(title_many)
        case 7
            disp('Select your own .mat files')
                        % Select .mat file
            % Prompt the user to select multiple files
            [fileName, pathName, filterIndex] = uigetfile('*.*', 'Select files from Data directory', ...
                '\\atlas\Chem\DuraFET\GDF Production\Data\mat Files\', ...
                'MultiSelect', 'on');
            
            % Check if the user clicked 'Cancel'
            if isequal(fileName, 0)
                disp('User canceled file selection.');
                return;
            end
            
            % If only one file is selected, uigetfile returns a char instead of a cell array
            if ischar(fileName)
                % Convert to a cell array for consistency
                fileName = {fileName};
            end
            
            % Display selected file names and path
            disp('Selected files:');
            disp(fileName);
            disp('Selected path:');
            disp(pathName);
        otherwise
            disp('Select a number between 1-6 ya filthy scallywag')
    end
end

%% Importing the "Many" Data
% Either personally selected files or predetermined Testlog subsets

if many == 7
    % Pre-allocation
    M.fP = zeros(100,length(fileName));
    for i = 1:length(fileName)
        T_many(i) = load([pathName,char(fileName(i))], myVars{:});
        M.pM = polyfit(T_many(i).g.pres(T_many(i).iPdecrease),T_many(i).g.fp_k2variable(T_many(i).iPdecrease),7);
        M.fP(:,i) = polyval(M.pM,O.pres);
    end
    title_many = 'User Selected Dataset';
else
    % Paths to extract the test log data and .mat files from the
    % upload_pTC_cal.m script
    M.path = '\\atlas\Chem\DuraFET\GDF Production\Performances\';
    M.testlog_path = '\\atlas\Chem\DuraFET\GDF Production\Test Log\Test Log.xlsx';
    % Copying Files
    copyfile(M.testlog_path,pwd);
    % Reading the manually updated Test Log
    M.testlog = readtable('Test Log.xlsx','PreserveVariableNames',true); 
    
    % TEST CODE NEXT 6 LINES
    % Sorting through test log and taking out any that are not what was
    % selected. 'AE' is test case
    % ind = M.testlog.AE;
    % ind(isnan(ind)) = 0; ind = logical(ind);
    % M.list_SN = M.testlog.SN(ind); 
    % M.list_FN = M.testlog.NameOfFolder(ind);
    
    % To do this for an input number, make INPUT string tied to number that
    % corresponds to column. The use:
    ind = M.testlog(:,ind_many); ind = table2array(ind);
    ind(isnan(ind)) = 0; ind = logical(ind);
    M.list_SN = M.testlog.SN(ind); 
    M.list_FN = M.testlog.NameOfFolder(ind);
    
    % Pre-allocation
    M.fP = zeros(100,height(M.list_SN));
    M.fP_inc = zeros(100,height(M.list_SN));
    M.pfM = zeros(height(M.list_SN),8);
    % Getting the fP (with pressure corrected k2) for both increasing and
    % decreasing pressure of "many" data to find median and std of all fits. 
    for i = 1:height(M.list_SN)
        T_many(i) = load([M.path,char(M.list_SN(i)),'\',char(M.list_FN(i)),'\',char(M.list_FN(i)),'.mat'], myVars{:});
        % Decreasing Pressure
        M.pM = polyfit(T_many(i).g.pres(T_many(i).iPdecrease),T_many(i).g.fp_k2variable(T_many(i).iPdecrease),7);
        M.pfM(i,:) = M.pM;
        M.fP(:,i) = polyval(M.pM,O.pres);
        % Increasing Pressure
        M.pM_inc = polyfit(T_many(i).g.pres(T_many(i).iPincrease),T_many(i).g.fp_k2variable(T_many(i).iPincrease),7);
        M.fP_inc(:,i) = polyval(M.pM_inc,O.pres);
        % Residuals of fit vs measurements per calibration run and relevant
        % diagnostics
        M.res(i).dec = (T_many(i).g.fp_k2variable(T_many(i).iPdecrease)) - polyval(M.pfM(i,:),T_many(i).g.pres(T_many(i).iPdecrease));
        M.res(i).avg = mean(abs(M.res(i).dec),1,'omitnan');
        M.res(i).dev = std(M.res(i).dec,1,'omitnan');
    end
    M.avg_res = [];
    for i = 1:height(M.list_SN)
        M.avg_res = [M.avg_res;M.res(i).avg];
    end

end

clear Var1 Var2 Var3 Var4 Var5 myVars ind

%% Calculation for the "Many" Data

% For DECREASING Pressure
mX   = median(M.fP,2,'omitnan')'; 
sX   = std(M.fP,[],2,'omitnan')';

bot_many = mX - 1.*sX; % 1 sigma bottom boundary
top_many = mX + 1.*sX; % 1 sigma top boundary

% For INCREASING Pressure
mX_inc   = median(M.fP_inc,2,'omitnan')';
sX_inc   = std(M.fP_inc,[],2,'omitnan')';

bot_many_inc = mX_inc - 1.*sX_inc; % 1 sigma bottom boundary
top_many_inc = mX_inc + 1.*sX_inc; % 1 sigma top boundary


%% Plotting

% OvM Decreasing Pressure Fit against "Many"
figure('Name','f(P) Fits Decreasing Pressure','unit', 'inch', 'position', [0.5 0.5 10.5 6.8]);
   % Plotting the shaded region
    inBetweenRegionX = [1:20.2:2020, 2000:-20.19:1];
    inBetweenRegionY = [top_many, fliplr(bot_many)];
    % Display the area first so it will be in the background.
    fill(inBetweenRegionX, inBetweenRegionY, [0.15 0.65 0.75]); alpha(0.25);
    hold on; box on; grid on;
    plot(O.pres,top_many,'LineWidth',2,'Color',[0.25 0.65 0.55]); hold on;
    plot(O.pres,bot_many,'LineWidth',2,'Color',[0.25 0.65 0.55]); hold on;
%     plot(O.pres,mX,'LineWidth',3.5,'Color',[0.25 0.65 0.55]); hold on;  
    plot(O.pres,O.fP1, '<','MarkerEdgeColor',[0.3010 0.5450 0.9330],'MarkerSize',12,'MarkerFaceColor',[0.05 0.15 0.9]);
    ylabel('f(P) [mV]','FontSize',18);
    xlim([0 2100]);
    xlabel('Pressure [dBar]','FontSize',18);
    title('f(P) Comparisons (Decreasing Pressure Data)');
    subtitle('Comparing the fit of "One" run vs. "Many"')
    f.ax = gca; f.ax.Color = [0.9 0.9 0.9];
    f.ax.FontSize = 16; f.ax.FontWeight = 'Bold';
    O.l = legend('"Many" Decreasing Data','+/- 1 Std Dev Bounds','','"One" Calibration Run');
    set(O.l,'Interpreter','none');
    f.a = annotation('textbox', [0.15, 0.23, 0, 0], 'String', append('"One" data comes from ', ...
        erase(fname_one,'.mat')),'FitBoxToText','on','FontSize',12,'BackgroundColor','w',...
        'Color',[0 0 0]);
    f.b = annotation('textbox', [0.15, 0.18, 0, 0], 'String', append('"Many" data comes from ', ...
        title_many),'FitBoxToText','on','FontSize',12,'BackgroundColor','w',...
        'Color',[0 0 0]);
    set(f.a,'Interpreter','none'); set(f.b,'Interpreter','none');

% OvM Increasing Pressure Fit against "Many"
figure('Name','f(P) Fits Increasing Pressure','unit', 'inch', 'position', [0.5 0.5 10.5 6.8]);
   % Plotting the shaded region
    inBetweenRegionX = [1:20.2:2020, 2000:-20.19:1];
    inBetweenRegionY = [top_many_inc, fliplr(bot_many_inc)];
    % Display the area first so it will be in the background.
    fill(inBetweenRegionX, inBetweenRegionY, [0.9290 0.6940 0.1250]); alpha(0.25);
    hold on; box on; grid on;
    plot(O.pres,top_many_inc,'LineWidth',2,'Color',[0.9290 0.5940 0.1250]); hold on;
    plot(O.pres,bot_many_inc,'LineWidth',2,'Color',[0.9290 0.5940 0.1250]); hold on;
%     plot(O.pres,mX,'LineWidth',3.5,'Color',[0.25 0.65 0.55]); hold on;  
    plot(O.pres,O.fP2, '>','MarkerEdgeColor',[0.8500 0.5250 0.0980],'MarkerSize',12,'MarkerFaceColor',[0.815 0.4 0.25]);
    ylabel('f(P) [mV]','FontSize',18);
    xlim([0 2100]);
    xlabel('Pressure [dBar]','FontSize',18);
    title('f(P) Comparisons (Increasing Pressure Data)');
    subtitle('Comparing the fit of "One" run vs. "Many"')
    f.ax = gca; f.ax.Color = [0.9 0.9 0.9];
    f.ax.FontSize = 16; f.ax.FontWeight = 'Bold';
    O.l = legend('"Many" Increasing Data','+/-1 Std Dev Bounds','','"One" Calibration Run');
    set(O.l,'Interpreter','none');
    f.a = annotation('textbox', [0.15, 0.23, 0, 0], 'String', append('"One" data comes from ', ...
        erase(fname_one,'.mat')),'FitBoxToText','on','FontSize',12,'BackgroundColor','w',...
        'Color',[0 0 0]);
    f.b = annotation('textbox', [0.15, 0.18, 0, 0], 'String', append('"Many" data comes from ', ...
        title_many),'FitBoxToText','on','FontSize',12,'BackgroundColor','w',...
        'Color',[0 0 0]);
    set(f.a,'Interpreter','none'); set(f.b,'Interpreter','none');

% OvM Scatter of all "One" Data    
figure('Name','Scatter of One vs Many','unit', 'inch', 'position', [0.5 0.5 10.5 6.8]);
   % Plotting the shaded region
    inBetweenRegionX = [1:20.2:2020, 2000:-20.19:1];
    inBetweenRegionY = [top_many, fliplr(bot_many)];
    % Display the area first so it will be in the background.
    fill(inBetweenRegionX, inBetweenRegionY, [0.15 0.65 0.75]); alpha(0.25);
    hold on; box on; grid on;
    plot(O.pres,top_many,'LineWidth',2,'Color',[0.25 0.65 0.55]); hold on;
    plot(O.pres,bot_many,'LineWidth',2,'Color',[0.25 0.65 0.55]); hold on;
    plot(T_one(:,1),T_one(:,2), 'o','MarkerEdgeColor',[0.3010 0.5450 0.9330],'MarkerSize',12,'MarkerFaceColor',[0.05 0.15 0.9]);
    ylabel({'f(P) [mV]';'(Vrs-Vct-Vcp-k2(P)*T-k0)'},'FontSize',18);
    xlim([0 2100]);
    xlabel('Pressure [dBar]','FontSize',18);
    title('f(P) Comparisons');
    subtitle('All data of "One" run vs. "Many"')
    f.ax = gca; f.ax.Color = [0.9 0.9 0.9];
    f.ax.FontSize = 16; f.ax.FontWeight = 'Bold';
    O.l = legend('"Many" Decreasing Data','+/- 1 Std Dev Bounds','','"One" Calibration Run');
    set(O.l,'Interpreter','none');
    f.a = annotation('textbox', [0.15, 0.23, 0, 0], 'String', append('"One" data comes from ', ...
        erase(fname_one,'.mat')),'FitBoxToText','on','FontSize',12,'BackgroundColor','w',...
        'Color',[0 0 0]);
    f.b = annotation('textbox', [0.15, 0.18, 0, 0], 'String', append('"Many" data comes from ', ...
        title_many),'FitBoxToText','on','FontSize',12,'BackgroundColor','w',...
        'Color',[0 0 0]);
    set(f.a,'Interpreter','none'); set(f.b,'Interpreter','none');

% OvM Increasing and Decreasing Fits
figure('Name','Hysteresis Analysis','unit', 'inch', 'position', [0.5 0.5 10.5 6.8]);
    box on; grid on;
    plot(O.pres,O.fP1, '<','MarkerEdgeColor',[0.3010 0.5450 0.9330],'MarkerSize',12,'MarkerFaceColor',[0.05 0.15 0.9]);
    hold on; grid on;
    plot(O.pres,O.fP2, '>','MarkerEdgeColor',[0.8500 0.5250 0.0980],'MarkerSize',12,'MarkerFaceColor',[0.815 0.4 0.25]);
    ylabel('f(P) [mV]','FontSize',18);
    xlim([0 2100]);
    xlabel('Pressure [dBar]','FontSize',18);
    title('Hysteresis Analysis');
    subtitle('Comparing f(P) fits for "One" run')
    f.ax = gca; f.ax.Color = [0.9 0.9 0.9];
    f.ax.FontSize = 16; f.ax.FontWeight = 'Bold';
    O.l = legend('Decreasing Pressure','Increasing Pressure');
    set(O.l,'Interpreter','none');
    f.a = annotation('textbox', [0.25, 0.88, 0, 0], 'String', append('"One" data comes from ', ...
        erase(fname_one,'.mat')),'FitBoxToText','on','FontSize',12,'BackgroundColor','w',...
        'Color',[0 0 0]);
    set(f.a,'Interpreter','none'); 

% OvM by Time
figure('Name','All Decreasing Pressure Data by Time','unit', 'inch', 'position', [0.5 0.5 10.5 6.8]);
   % Plotting the shaded region
    inBetweenRegionX = [1:20.2:2020, 2000:-20.19:1];
    inBetweenRegionY = [top_many, fliplr(bot_many)];
    % Display the area first so it will be in the background.
    fill(inBetweenRegionX, inBetweenRegionY, [0.15 0.65 0.75]); alpha(0.25);
    hold on; box on; grid on;
    plot(O.pres,top_many,'LineWidth',2,'Color',[0.25 0.65 0.55]); hold on;
    plot(O.pres,bot_many,'LineWidth',2,'Color',[0.25 0.65 0.55]); hold on;
    scatter(T_one(:,1), T_one(:,2), 14, T_one(:,5),'SizeData',100,'LineWidth',2);
    colormap(flipud(cbrewer('div', 'RdYlBu', 1000, 'PCHIP')));   
    cbt = colorbar;
    ylabel(cbt,'Days Elapsed');
    ylabel({'f(P) [mV]';'(Vrs-Vct-Vcp-k2(P)*T-k0)'},'FontSize',18);
    xlim([0 2100]);
    xlabel('Pressure [dBar]','FontSize',18);
    title('All computed f(P) Data by Time (Decreasing Pressure)');
    subtitle('Comparing the fit of "One" run vs. "Many"')
    f.ax = gca; f.ax.Color = [0.9 0.9 0.9];
    f.ax.FontSize = 16; f.ax.FontWeight = 'Bold';
    O.l = legend('"Many" Decreasing Data','+/- 1 Std Dev Bounds','','All "One" Data by Time');
    set(O.l,'Interpreter','none');
    f.a = annotation('textbox', [0.15, 0.23, 0, 0], 'String', append('"One" data comes from ', ...
        erase(fname_one,'.mat')),'FitBoxToText','on','FontSize',12,'BackgroundColor','w',...
        'Color',[0 0 0]);
    f.b = annotation('textbox', [0.15, 0.18, 0, 0], 'String', append('"Many" data comes from ', ...
        title_many),'FitBoxToText','on','FontSize',12,'BackgroundColor','w',...
        'Color',[0 0 0]);
    set(f.a,'Interpreter','none'); set(f.b,'Interpreter','none');

%% Vrs  and Ik/Ib of "One"

% Vrs  and Ik/Ib
f5 = figure('unit', 'inch', 'position', [1 1 16 9]);
t = tiledlayout(2,1,'TileSpacing','compact');
% pres and temp
ax(1) = nexttile;
hold on; box on;
ax = gca; 
ax.FontSize = 18;    
ax.FontWeight = 'Bold';
yyaxis left
    plot(O.g.DT, O.g.Vrs, 'b','LineWidth',2);
    ylabel('Vrs');
yyaxis right
    plot(O.g.DT, O.g.pres, 'r','LineWidth',2);
    ylabel('pressure');
ax(2) = nexttile;
ax = gca; 
ax.FontSize = 18;    
ax.FontWeight = 'Bold';
yyaxis left
    plot(O.g.DT, O.g.Ik.*1E9, 'b','LineWidth',2);
    ylabel('Ik (nA)');
yyaxis right
    plot(O.g.DT, O.g.Ib.*1E9, 'r','LineWidth',2);
    ylabel('Ib (nA)'); 
linkaxes(ax, 'x');    
title(t,erase(fname_one,'.mat'),'FontSize',18,'FontWeight','bold');
xlabel(t,'Time','FontSize',18,'FontWeight','bold')

%% Residual Plotting

figure('Name','Decreasing Pressure Residuals for One Data','unit', 'inch', 'position', [0.5 0.5 10.5 6.8]);
    % Display the area first so it will be in the background.
    plot(O.g.pres(O.iPdecrease),abs(O.res.dec),'o','MarkerEdgeColor',[0.010 0.0450 0.0330],...
        'MarkerSize',12,'MarkerFaceColor',[0.25 0.55 0.9]);
    hold on; box on; grid on; 
    yline(mean(M.avg_res,1,'omitnan'),':',...
    ['Residual Mean of ',title_many],'LineWidth',5,'Color','r'); hold on;
    ylabel('Residuals of Fit vs. Measurements [mV]','FontSize',18);
    xlim([0 2100]);
    xlabel('Pressure [dBar]','FontSize',18);
    title('f(P) Residuals (Decreasing Pressure Data)');
    subtitle('Comparing the fit of "One" run vs. "Many"')
    f.ax = gca; f.ax.Color = [0.9 0.9 0.9];
    f.ax.FontSize = 16; f.ax.FontWeight = 'Bold';
    O.l = legend('"One" Residuals','Mean of "Many"');
    set(O.l,'Interpreter','none');

%% Hysteresis Calcualtions
% Using the function to grab histogram data of the residuals for each
% pressurization cycle to evaluate hysteresis for One against Many

% Function to return figures
[f_one,One] = HysteresisCalcV3(O,fname_one);
[f_many,Many] = HysteresisCalcV3(T_many,M.list_FN);

%% Hysteresis Visualization & Plotting

% f1 = figure('Name','Mean Vrs Fit Residuals Bi-Directional','unit', 'inch', 'position', [0.5 0.5 11.5 7]);
% yyaxis left
% histogram(Many.allMeanResidualBoth*1000, Many.edges_mean*1000); hold on;
% xline(1000*(Many.mean_mean_both + Many.std_mean_both),':',...
%     ['Std Dev ',title_many],'LineWidth',2.5,'Color','#0072BD'); hold on;
% xline(Many.mean_mean_both,'--','LineWidth',2.5,'Color','#0072BD');
% ylabel('Number of Cycles')
% yyaxis right
% histogram(One.allMeanResidualBoth*1000, One.edges_mean*1000);
% grid on;
% f.l = legend([title_many,' Mean Vrs Fit Residuals - Both'],'',['Mean Vrs Fit (Both) from ',...
%     title_many],[erase(fname_one,'.mat'),' Mean Vrs Fit Residuals - Both']);
% set(f.l,'Interpreter','none'); f.l.Location  = 'east';
% 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';
% title('Histogram of Mean Vrs Per Cycle vs. The Bi-Directional Fit')

f2 = figure('Name','Mean Vrs Fit Residuals Decreasing Pressure','unit', 'inch', 'position', [0.5 0.5 11.5 7]);
yyaxis left
histogram(Many.allMeanResidualDown*1000, Many.edges_mean*1000); hold on;
xline(1000*(Many.mean_mean_down + Many.std_mean_down),':',...
    ['Std Dev ',title_many],'LineWidth',2.5,'Color',[0 0 1]); hold on;
xline(Many.mean_mean_down*1000,'--','LineWidth',2.5,'Color','#0072BD');
ylabel('Number of Cycles')
yyaxis right
histogram(One.allMeanResidualDown*1000, One.edges_mean*1000);
grid on;
f.l = legend([title_many,' Mean Vrs Fit Residuals - Dec. Pres.'],['Std Dev of Mean Vrs Fit (Dec. Pres.) from ',...
    title_many],['Mean Vrs Fit (Dec. Pres.) from ',title_many],[erase(fname_one,'.mat'),...
    ' Mean Vrs Fit Residuals - Dec. Pres.']);
set(f.l,'Interpreter','none'); f.l.Location  = 'east';
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';
title('Histogram of Mean Vrs Per Cycle vs. The Decreasing Pressure Fit')

%% Mean OPS+ (Davis-Trout Number)
% 
% % OPS Adjusted score
% ops = 100.*((O.res.med)/(mean(M.avg_res,1,'omitnan')) + ...
%     One.mean_mean_down/Many.mean_mean_down ...
%     - 1);
% 
% if ops > 1
%     disp(['OPS+ score is ',num2str(ops)]);
% else
%     ops  = ops*-1;
%     disp(['OPS+ very low, recalculated score is ',num2str(ops)]);
% end

%% Median OPS+ (Davis-Trout Number)

% OPS Adjusted score
ops_med = 100.*((O.res.med)/(median(M.avg_res,1,'omitnan')) + ...
    One.med_mean_down/Many.med_mean_down ...
    - 1);

if ops_med > 1
    disp(['OPS+ (median) score is ',num2str(ops_med)]);
else
    ops_med  = ops_med*-1;
    disp(['OPS+ (median) very low, recalculated score is ',num2str(ops_med)]);
end

%% Residual Values

fprintf(['Median residual of decr. pressure from fit is ',num2str(O.res.med,'%1.3f'), ' mV ' ...
    '(',num2str(mean(M.avg_res,1,'omitnan'),'%1.3f'), ' mV for ',title_many,') \n']);

%% k2 value

% Displaying k2 value at atmospheric and at pressure
disp(['k2 (atmosphereic) is ',num2str(O.k2_1atm*1000,'%1.3f'), ' mV/C']);
disp(['k2 (2000 dBar) is ',num2str(O.k2_2000*1000,'%1.3f'), ' mV/C']);

%% Hysteresis Values

% Displaying mean and max hysteresis for each cycle
fprintf(['Mean hysteresis for each cycle is ',num2str(One.hysteresisMean*1000,'%1.3f'), ' mV '...
    '(',num2str(Many.hysteresisMean*1000,'%1.3f'),' mV for ',title_many,') \n']);
disp(['Median of max hysteresis value per cycle is ',num2str(One.hysteresisMax*1000,'%1.3f'), ' mV']);

%% f(P) magnitude

disp(['f(P) decr. pressure magnitude is ',num2str(range(O.fP2),'%1.3f'), ' mV']);
