%% Load data and Variables Definition

load('Colors') %JJ palette

%%%% OLD (not usefull) datasets %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

load('CaCO3_2') %: 55,56,58,59,63
load('CaCO3_61') %: 61 are in Ca_differentPS

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%% Active datasets %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

load ('CaCO3_activeData')
%it contains:
%load('CaCO3_3') : 71,65,66 in scan3
%load('StephJune'): 72,74,65C,66B in scanJune
%load('StephDecember')%: 80,82,65B E,66A C,85,86 in scanDec
%load('Ca_differentPS') %%61s
%load('WhitePaperComparsions') %from published by JJ

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%Graphs to show%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
close all;
%put 1 to the additive you want to see
PS=0;
ASP=0;
GLY=0;
VES=0;
unheated=1;
otherWhites=0;
comparisonAdditives=0;
heatedCalcite=0;

graph1=1; %put 1 here to see all T1 plots
graph2=0; %put 1 here to see all T2 plots
graph3=0;
graph4=0;
graph5=0;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%Variables Definition%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
wl=scan3(1).wl(:,1);
back=scan3(4).spec(:,1);
spectraUnheated=zeros(length(scan),2068);
spectra3=zeros(length(scan3),2068);
spectraJune=zeros(length(scanJune),2068);
spectraDec=zeros(length(scanDec),2068);
spectraRedoGly=zeros(length(scanRedoGly),2068);

for i=1:length(scan)
    spectraUnheated(i,:)=smooth((scan(i).spec(:,1)),20)';
end
for i=1:length(scan3)
    spectra3(i,:)=smooth((scan3(i).spec(:,1)),20)';
end
for i=1:length(scanJune)
    spectraJune(i,:)=smooth((scanJune(i).spec(:,1)),20)';
end
for i=1:length(scanDec)
    spectraDec(i,:)=smooth((scanDec(i).spec(:,1)),20)';
end
for i=1:length(scanRedoGly)
    spectraRedoGly(i,:)=smooth((scanDec(i).spec(:,1)),20)';
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%load CaCO3%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
b65t1=mean(spectra3(1:3,:));
c65t1 =mean(spectra3(5:7,:));
e65t1=mean(spectra3(8:11,:));
f65t1=mean(spectra3(12:13,:));
f65t2=mean(spectra3(17:19,:));
%a63 = smooth( ( (scanCalciteControl(53,:).spec(:,1) - scanCalciteControl(52,:).spec(:,1)) + (scanCalciteControl(54,:).spec(:,1) - scanCalciteControl(52,:).spec(:,1)) + (scanCalciteControl(55,:).spec(:,1) - scanCalciteControl(52,:).spec(:,1)) )/3, 20);
a63 = smooth( ( (scanCalciteControl(53,:).spec(:,1) - scanCalciteControl(52,:).spec(:,1)) + (scanCalciteControl(54,:).spec(:,1) - scanCalciteControl(52,:).spec(:,1)) )/2, 20);

%a66t1=mean(spectraRedoGly(1:3,:));
a66t1=mean(spectra3(21:22,:)); %crystal lost
%c66t1=mean(spectraRedoGly(4:6,:));
c66t1=mean(spectra3(26:28,:));
d66t1=mean(spectra3(29:32,:));
b66t2=mean(spectra3(33:35,:));
d66t2=mean(spectra3(36:37,:));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%load StephJune%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
step=4;

    initialIndex=4;
    b71=[mean(spectraJune(initialIndex:step:initialIndex+8,:))',mean(spectraJune(initialIndex+1:step:initialIndex+1+8,:))',mean(spectraJune(initialIndex+2:step:initialIndex+2+8,:))',mean(spectraJune(initialIndex+3:step:initialIndex+3+8,:))'];

    initialIndex=initialIndex+step*3;
    b66t1=[mean(spectraJune(initialIndex:step:initialIndex+8,:))',mean(spectraJune(initialIndex+1:step:initialIndex+1+8,:))',mean(spectraJune(initialIndex+2:step:initialIndex+2+8,:))',mean(spectraJune(initialIndex+3:step:initialIndex+3+8,:))'];

    initialIndex=initialIndex+step*3;
    c65t2=[mean(spectraJune(initialIndex:step:initialIndex+8,:))',mean(spectraJune(initialIndex+1:step:initialIndex+1+8,:))',mean(spectraJune(initialIndex+2:step:initialIndex+2+8,:))',mean(spectraJune(initialIndex+3:step:initialIndex+3+8,:))'];

    initialIndex=initialIndex+step*3;
    a72=[mean(spectraJune(initialIndex:step:initialIndex+8,:))',mean(spectraJune(initialIndex+1:step:initialIndex+1+8,:))',mean(spectraJune(initialIndex+2:step:initialIndex+2+8,:))',mean(spectraJune(initialIndex+3:step:initialIndex+3+8,:))'];
    
    initialIndex=initialIndex+step*3;
    b72=[mean(spectraJune(initialIndex:step:initialIndex+8,:))',mean(spectraJune(initialIndex+1:step:initialIndex+1+8,:))',mean(spectraJune(initialIndex+2:step:initialIndex+2+8,:))',mean(spectraJune(initialIndex+3:step:initialIndex+3+8,:))'];
  
    initialIndex=initialIndex+step*3;
    a74=[mean(spectraJune(initialIndex:step:initialIndex+8,:))',mean(spectraJune(initialIndex+1:step:initialIndex+1+8,:))',mean(spectraJune(initialIndex+2:step:initialIndex+2+8,:))',mean(spectraJune(initialIndex+3:step:initialIndex+3+8,:))'];

    initialIndex=initialIndex+step*3;
    b74=[mean(spectraJune(initialIndex:step:initialIndex+8,:))',mean(spectraJune(initialIndex+1:step:initialIndex+1+8,:))',mean(spectraJune(initialIndex+2:step:initialIndex+2+8,:))',mean(spectraJune(initialIndex+3:step:initialIndex+3+8,:))'];

    initialIndex=initialIndex+step*3;
    c74=[mean(spectraJune(initialIndex:step:initialIndex+8,:))',mean(spectraJune(initialIndex+1:step:initialIndex+1+8,:))',mean(spectraJune(initialIndex+2:step:initialIndex+2+8,:))',mean(spectraJune(initialIndex+3:step:initialIndex+3+8,:))'];

    initialIndex=initialIndex+step*3;
    d74=[mean(spectraJune(initialIndex:step:initialIndex+8,:))',mean(spectraJune(initialIndex+1:step:initialIndex+1+8,:))',mean(spectraJune(initialIndex+2:step:initialIndex+2+8,:))',mean(spectraJune(initialIndex+3:step:initialIndex+3+8,:))'];
         
    initialIndex=initialIndex+step*3;
    e74=[mean(spectraJune(initialIndex:step:initialIndex+8,:))',mean(spectraJune(initialIndex+1:step:initialIndex+1+8,:))',mean(spectraJune(initialIndex+2:step:initialIndex+2+8,:))',mean(spectraJune(initialIndex+3:step:initialIndex+3+8,:))'];

    initialIndex=initialIndex+step*3;
    f74=[mean(spectraJune(initialIndex:step:initialIndex+8,:))',mean(spectraJune(initialIndex+1:step:initialIndex+1+8,:))',mean(spectraJune(initialIndex+2:step:initialIndex+2+8,:))',mean(spectraJune(initialIndex+3:step:initialIndex+3+8,:))'];
       
    initialIndex=initialIndex+step*3+2;
    Skel=[mean(spectraJune(initialIndex:step:initialIndex+8,:))',mean(spectraJune(initialIndex+1:step:initialIndex+1+8,:))',mean(spectraJune(initialIndex+2:step:initialIndex+2+8,:))',mean(spectraJune(initialIndex+3:step:initialIndex+3+8,:))'];
 
    initialIndex=initialIndex+step*3;%150
    %Spline=[mean(spectraJune(initialIndex:step:initialIndex+8,:))',mean(spectraJune(initialIndex+1:step:initialIndex+1+8,:))',mean(spectraJune(initialIndex+2:step:initialIndex+2+8,:))',mean(spectraJune(initialIndex+3:step:initialIndex+3+8,:))'];
    Spline=mean(spectraJune(initialIndex:step:initialIndex+8,:))'; %crystal lost
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%load StephDec%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    initialIndex=1;
    step=1;
    a80t1=mean(spectraDec(initialIndex:step:initialIndex+2,:))';
    
    initialIndex=4;
    step=1;
    b80t1=mean(spectraDec(initialIndex:step:initialIndex+2,:))';
    
    initialIndex=7;
    step=1;
    c80t1=mean(spectraDec(initialIndex:step:initialIndex+2,:))';
    
    
    initialIndex=10;
    step=1;
    a82t1=mean(spectraDec(initialIndex:step:initialIndex+2,:))';

    initialIndex=13;
    step=1;
    b82t1=mean(spectraDec(initialIndex:step:initialIndex+2,:))';
    
    initialIndex=16;
    step=1;
    c82t1=mean(spectraDec(initialIndex:step:initialIndex+2,:))';    
    
    initialIndex=19;
    step=1;
    d82t1=mean(spectraDec(initialIndex:step:initialIndex+2,:))';
    
    initialIndex=22;
    step=1;
    b65t2=mean(spectraDec(initialIndex:step:initialIndex+2,:))';    
    
    initialIndex=25;
    step=1;
    e65t2=mean(spectraDec(initialIndex:step:initialIndex+2,:))';      
 
    initialIndex=28;
    step=1;
    a66t2=mean(spectraDec(initialIndex:step:initialIndex+2,:))';       
    
    initialIndex=31;
    step=1;
    c66t2=mean(spectraDec(initialIndex:step:initialIndex+2,:))';  
    
    initialIndex=34;
    step=1;
    a85t1=mean(spectraDec(initialIndex:step:initialIndex+2,:))';       
    
    initialIndex=37;
    step=1;
    b85t1=mean(spectraDec(initialIndex:step:initialIndex+2,:))'; 
    
    initialIndex=40;
    step=1;
    c85t1=mean(spectraDec(initialIndex:step:initialIndex+2,:))';       
 
    initialIndex=43;
    step=1;
    d85t1=mean(spectraDec(initialIndex:step:initialIndex+2,:))';       
    
    initialIndex=46;
    step=1;
    e85t1=mean(spectraDec(initialIndex:step:initialIndex+2,:))';     
    
    initialIndex=49;
    step=1;
    %a86t1=mean(spectraDec(initialIndex:step:initialIndex+2,:))'; 
    a86t1=mean(spectraDec(initialIndex:step:initialIndex+1,:))'; %crystal lost
 
    initialIndex=52;
    step=1;
    b86t1=mean(spectraDec(initialIndex:step:initialIndex+2,:))';   
    
    initialIndex=55;
    step=1;
    a85t2=mean(spectraDec(initialIndex:step:initialIndex+2,:))';         

    initialIndex=58;
    step=1;
    b85t2=mean(spectraDec(initialIndex:step:initialIndex+2,:))';      
  
    initialIndex=61;
    step=1;
    c85t2=mean(spectraDec(initialIndex:step:initialIndex+2,:))';    
    
    initialIndex=64;
    step=1;
    d85t2=mean(spectraDec(initialIndex:step:initialIndex+2,:))'; 
    
    initialIndex=67;
    step=1;
    %e85t2=mean(spectraDec(initialIndex:step:initialIndex+2,:))'; 
    e85t2=spectraDec(69,:)'; %crystal lost

    initialIndex=70;
    step=1;
    a86t2=mean(spectraDec(initialIndex:step:initialIndex+2,:))';       
    
    initialIndex=73;
    step=1;
    b86t2=mean(spectraDec(initialIndex:step:initialIndex+2,:))';    
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
%% PS
if PS
    
%% Graph1: Including all samples for PS
if graph1
    
    % Create two columns of data
    A = wl;
    B = smooth(f74(:,4)*(2*f74(1200,3)/f74(1200,4)),30);
    % Create a table with the data and variable names
    T = table(A, B,'VariableNames',{'A', 'B'});
    % Write data to text file
    writetable(T,"74F_T1.txt")

    plot(wl,smooth(a80t1,30),'LineWidth',7,'Color',"#F8E7B9"); hold on;
    plot(wl,smooth(b80t1,30),'LineWidth',7,'Color',"#F5CE8D"); hold on;
    plot(wl,smooth(c80t1,30),'LineWidth',7,'Color',"#F2B153"); hold on;
    k=4;
    plot(wl,smooth(a74(:,k)*(2*a74(1200,3)/a74(1200,4)),30),'LineWidth',7,'Color',"#FA9B2F"); hold on;
    plot(wl,smooth(b74(:,k)*(2*b74(1200,3)/b74(1200,4)),30),'LineWidth',7,'Color',"#FD8209"); hold on;
    plot(wl,smooth(c74(:,k)*(2*c74(1200,3)/c74(1200,4)),30),'LineWidth',7,'Color',"#FF4B01"); hold on;
    plot(wl,smooth(d74(:,k)*(2*d74(1200,3)/d74(1200,4)),30),'LineWidth',7,'Color',"#F91304"); hold on;
    plot(wl,smooth(e74(:,k)*(2*e74(1200,3)/e74(1200,4)),30),'LineWidth',7,'Color',"#BC0002"); hold on;
    plot(wl,smooth(f74(:,k)*(2*f74(1200,3)/f74(1200,4)),30),'LineWidth',7,'Color',"#630601"); hold on;
    xlim([450,750]);
    %legend({'80 A, t=7µm', '80 B, t=4µm', '80 C, t=13µm','74A, t=1.8µm','74B, t=2µm','74C, t=3.5µm','74D, t=2.5µm','74E, t=5.5µm','74F, t=7µm'},'FontSize',25,'Location','northwest')
    legend({'0.005 wt.% PS', '0.01 wt.% PS', '0.1 wt.% PS', '0.5 wt.% PS', '1 wt.% PS', '2 wt.% PS', '3 wt.% PS', '4 wt.% PS', '5 wt.% PS'}, 'FontSize',25,'Location','northwest')

end

if graph2

    k=4;
    plot(wl,smooth(a74(:,k)*(2*a74(1200,3)/a74(1200,4)),30),'LineWidth',7,'Color',b); hold on;
    plot(wl,smooth(b74(:,k)*(2*b74(1200,3)/b74(1200,4)),30),'LineWidth',7,'Color',c); hold on;
    plot(wl,smooth(c74(:,k)*(2*c74(1200,3)/c74(1200,4)),30),'LineWidth',7,'Color',g); hold on;
    plot(wl,smooth(d74(:,k)*(2*d74(1200,3)/d74(1200,4)),30),'LineWidth',7,'Color',r); hold on;
    plot(wl,smooth(e74(:,k)*(2*e74(1200,3)/e74(1200,4)),30),'LineWidth',7,'Color',o); hold on;
    plot(wl,smooth(f74(:,k)*(2*f74(1200,3)/f74(1200,4)),30),'LineWidth',7,'Color',p); hold on;
    xlim([450,750]);
    legend({'74A, t=1.8µm','74B, t=2µm','74C, t=3.5µm','74D, t=2.5µm','74E, t=5.5µm','74F, t=7µm'},'FontSize',25,'Location','northwest')


end

if graph3

    plot(wl,smooth(a80t1,30),'--','LineWidth',7,'Color',b); hold on;
    plot(wl,smooth(b80t1,30),'--','LineWidth',7,'Color',c); hold on;
    plot(wl,smooth(c80t1,30),'--','LineWidth',7,'Color',r); hold on;
    legend({'80 A, t=7µm', '80 B, t=4µm', '80 C, t=13µm'},'FontSize',25,'Location','northwest')


end

if graph5

    thicknesses = [7.039, 3.866, 2.776, 1.785, 2.012, 3.541, 2.457, 5.388, 6.996];
    ze = [2.5, 2.6, 2, 1.8, 2, 1.76, 2.1, 2, 2.2];

    plot(wl,smooth(a80t1,30),'--','LineWidth',7,'Color',b); hold on;
    plot(wl,smooth(b80t1,30),'--','LineWidth',7,'Color',c); hold on;
    plot(wl,smooth(c80t1,30),'--','LineWidth',7,'Color',r); hold on;
    k=4;
    plot(wl,smooth(a74(:,k)*(2*a74(1200,3)/a74(1200,4)),30),'LineWidth',7,'Color',b); hold on;
    plot(wl,smooth(b74(:,k)*(2*b74(1200,3)/b74(1200,4)),30),'LineWidth',7,'Color',c); hold on;
    plot(wl,smooth(c74(:,k)*(2*c74(1200,3)/c74(1200,4)),30),'LineWidth',7,'Color',g); hold on;
    plot(wl,smooth(d74(:,k)*(2*d74(1200,3)/d74(1200,4)),30),'LineWidth',7,'Color',r); hold on;
    plot(wl,smooth(e74(:,k)*(2*e74(1200,3)/e74(1200,4)),30),'LineWidth',7,'Color',o); hold on;
    plot(wl,smooth(f74(:,k)*(2*f74(1200,3)/f74(1200,4)),30),'LineWidth',7,'Color',p); hold on;
    xlim([450,750]);
    legend({'80 A, t=7µm', '80 B, t=4µm', '80 C, t=13µm','74A, t=1.8µm','74B, t=2µm','74C, t=3.5µm','74D, t=2.5µm','74E, t=5.5µm','74F, t=7µm'},'FontSize',25,'Location','northwest')

    % Preparing the data for scatter plot
    % Colors for each dataset
    colors = {'b', 'c', 'r', 'b', 'c', 'g', 'r', 'o', 'p'};
    mfp_results = zeros(1, length(thicknesses)); % Store mfp results
    for i = 1:length(thicknesses)
        if i == 1
            X = smooth(a80t1,30);
        elseif i == 2
            X = smooth(b80t1,30);
        elseif i == 3
            X = smooth(c80t1,30);
        elseif i == 4
            X = smooth(a74(:,k)*(2*a74(1200,3)/a74(1200,4)),30);
        elseif i == 5
            X = smooth(b74(:,k)*(2*b74(1200,3)/b74(1200,4)),30);
        elseif i == 6
            X = smooth(c74(:,k)*(2*c74(1200,3)/c74(1200,4)),30);
        elseif i == 7
            X = smooth(d74(:,k)*(2*d74(1200,3)/d74(1200,4)),30);
        elseif i == 8
            X = smooth(e74(:,k)*(2*e74(1200,3)/e74(1200,4)),30);
        else % i == 9
            X = smooth(f74(:,k)*(2*f74(1200,3)/f74(1200,4)),30);
        end
        
        % Apply mfp function
        indices = wl >= 400 & wl <= 800;
        meanX = mean(X(indices));
        mfp_results(i) = mfp(meanX, thicknesses(i), ze(i));
    end
    
    % Define the labels for each dataset, excluding the 't=' part
    datasetNames = {'80 A', '80 B', '80 C', '74A', '74B', '74C', '74D', '74E', '74F'};
    
    % Scatter plot for each dataset's mfp result
    figure;
    hold on;
    
    % Calculate positions for each dataset label on the x-axis
    xPositions = 1:length(mfp_results);
    
    for i = 1:length(mfp_results)
        scatter(xPositions(i), 1/mfp_results(i), 200, colors{i}, 'filled'); % Using 100 as arbitrary size for visibility
    end
    
    hold off;
    
    % Adjusting the scatter plot appearance
    xlim([0, length(mfp_results)+1]);
    set(gca, 'xtick', xPositions, 'xticklabel', datasetNames); % Set custom x-axis tick labels
    xlabel('Dataset');
    ylabel('Scattering efficiency');
    %title('Scatter Plot of mfp Results for Each Dataset');
    legend(datasetNames, 'FontSize', 12, 'Location', 'northwest');
        legend({'80 A, t=7µm', '80 B, t=4µm', '80 C, t=13µm','74A, t=1.8µm','74B, t=2µm','74C, t=3.5µm','74D, t=2.5µm','74E, t=5.5µm','74F, t=7µm'}, 'FontSize', 12, 'Location', 'northwest');
end

end 

%% ASP
if ASP
    
%% Graph1: Including all T1 samples for ASP
if graph1
    
    % Create two columns of data
    A = wl;
    B = smooth(c85t1,30);
    % Create a table with the data and variable names
    T = table(A, B,'VariableNames',{'A', 'B'});
    % Write data to text file
    writetable(T,"85C_T1.txt")

    plot(wl,smooth(d85t1,30),'LineWidth',7,'Color','#DDDB64'); hold on;
    plot(wl,smooth(b65t1,30),'LineWidth',7,'Color','#A5CC59'); hold on;
    plot(wl,smooth(c65t1,30),'LineWidth',7,'Color','#55DB84'); hold on;
    plot(wl,smooth(e85t1,30),'LineWidth',7,'Color','#F0D469'); hold on;
    plot(wl,smooth(e65t1,30),'LineWidth',7,'Color','#EECF65'); hold on;
    plot(wl,smooth(f65t1,30),'LineWidth',7,'Color','#EF8328'); hold on;
    plot(wl,smooth(a85t1,30),'LineWidth',7,'Color','#FEC0C1'); hold on;
    plot(wl,smooth(b85t1,30),'LineWidth',7,'Color','#FA8274'); hold on;
    plot(wl,smooth(c85t1,30),'LineWidth',7,'Color','#F95367'); hold on;

    %legend({'85D T1, t=9µm','65B T1, t=9µm','65C T1, t=9µm','85E T1, t=13µm','65E T1, t=3µm','65F T1, t=4µm','85A T1, t=10µm','85B T1, t=9µm','85C T1, t=18µm'},'FontSize',25,'Location','northeast')
    legend({'[Ca^2^+] = 1.25 mM, [Asp] = 10 mM','[Ca^2^+] = 1.25 mM, [Asp] = 20 mM','[Ca^2^+] = 1.25 mM, [Asp] = 50 mM','[Ca^2^+] = 5 mM, [Asp] = 10 mM','[Ca^2^+] = 5 mM, [Asp] = 20 mM','[Ca^2^+] = 5 mM, [Asp] = 50 mM','[Ca^2^+] = 10 mM, [Asp] = 10 mM','[Ca^2^+] = 10 mM, [Asp] = 20 mM','[Ca^2^+] = 10 mM, [Asp] = 50 mM'},'FontSize',12,'Location','northeast')

end

%% Graph2: Including all T2 samples for ASP
if graph2
    
    % Create two columns of data
    A = wl;
    B = smooth(c85t2,30);
    % Create a table with the data and variable names
    T = table(A, B,'VariableNames',{'A', 'B'});
    % Write data to text file
    writetable(T,"85C_T2.txt")

    plot(wl,smooth(d85t2,30),'LineWidth',7,'Color','#DDDB64'); hold on;
    plot(wl,smooth(b65t2,30),'LineWidth',7,'Color','#A5CC59'); hold on;
    %plot(wl,smooth(c65t2(:,4),30),'--','LineWidth',7,'Color',g); hold on;
    plot(wl,smooth(e85t2,30),'LineWidth',7,'Color','#F0D469'); hold on;
    plot(wl,smooth(e65t2,30),'LineWidth',7,'Color','#EECF65'); hold on;
    plot(wl,smooth(f65t2,30),'LineWidth',7,'Color','#EF8328'); hold on;
    plot(wl,smooth(a85t2,30),'LineWidth',7,'Color','#FEC0C1'); hold on;
    plot(wl,smooth(b85t2,30),'LineWidth',7,'Color','#FA8274'); hold on;
    plot(wl,smooth(c85t2,30),'LineWidth',7,'Color','#F95367'); hold on;

    %legend({'85D t2, t=5.4µm','65B t2, t=11µm','85E t2, t=8µm','65E t2, t=4.5µm','65F t2, t=4µm','85A t2, t=10µm','85B t2, t=8µm','85C t2, t=12µm'},'FontSize',25,'Location','northeast')
    legend({'[Ca^2^+] = 1.25 mM, [Asp] = 10 mM','[Ca^2^+] = 1.25 mM, [Asp] = 20 mM','[Ca^2^+] = 5 mM, [Asp] = 10 mM','[Ca^2^+] = 5 mM, [Asp] = 20 mM','[Ca^2^+] = 5 mM, [Asp] = 50 mM','[Ca^2^+] = 10 mM, [Asp] = 10 mM','[Ca^2^+] = 10 mM, [Asp] = 20 mM','[Ca^2^+] = 10 mM, [Asp] = 50 mM'},'FontSize',12,'Location','northeast')

end

%% Graph3: Including all T1 1.25mM Ca samples samples for ASP
if graph3

    plot(wl,smooth(d85t1,30),'--','LineWidth',7,'Color',b); hold on;
    plot(wl,smooth(b65t1,30),'--','LineWidth',7,'Color',c); hold on;
    plot(wl,smooth(c65t1,30),'--','LineWidth',7,'Color',r); hold on;
    legend({'85D T1, t=9µm','65B T1, t=9µm','65C T1, t=9µm','85E T1, t=13µm','65E T1, t=3µm','65F T1, t=4µm','85A T1, t=10µm','85B T1, t=9µm','85C T1, t=18µm'},'FontSize',25,'Location','northeast')


end

%% Graph4: Including all T2 1.25mM Ca samples samples for ASP
if graph4

    plot(wl,smooth(d85t2,30),'--','LineWidth',7,'Color',b); hold on;
    plot(wl,smooth(b65t2,30),'--','LineWidth',7,'Color',c); hold on;
    plot(wl,smooth(c65t2(:,4),30),'--','LineWidth',7,'Color',r); hold on;
    legend({'85D t2, t=5.4µm','65B t2, t=11µm','65C t2, t=6.4µm','85E t2, t=8µm','65E t2, t=4.5µm','65F t2, t=4µm','85A t2, t=10µm','85B t2, t=8µm','85C t2, t=12µm'},'FontSize',25,'Location','northeast')


end

if graph5

    thicknesses = [8.8,9.02,13.4,13.1,3.1,3.8,9.6,9.1,18.5];
    ze = [2.5,2.3,2.6,2.4,2.6,2.3,2.4,2.4,1.9];
    

    plot(wl,smooth(d85t1,30),'--','LineWidth',7,'Color',b); hold on;
    plot(wl,smooth(b65t1,30),'--','LineWidth',7,'Color',c); hold on;
    plot(wl,smooth(c65t1,30),'--','LineWidth',7,'Color',r); hold on;
    plot(wl,smooth(e85t1,30),'LineWidth',7,'Color',b); hold on;
    plot(wl,smooth(e65t1,30),'LineWidth',7,'Color',c); hold on;
    plot(wl,smooth(f65t1,30),'LineWidth',7,'Color',g); hold on;
    plot(wl,smooth(a85t1,30),'LineWidth',7,'Color',r); hold on;
    plot(wl,smooth(b85t1,30),'LineWidth',7,'Color',o); hold on;
    plot(wl,smooth(c85t1,30),'LineWidth',7,'Color',p); hold on;

    legend({'85D T1, t=9µm','65B T1, t=9µm','65C T1, t=9µm','85E T1, t=13µm','65E T1, t=3µm','65F T1, t=4µm','85A T1, t=10µm','85B T1, t=9µm','85C T1, t=18µm'},'FontSize',25,'Location','northeast')

    % Colors for each dataset
    colors = {'b', 'c', 'r', 'b', 'c', 'g', 'r', 'o', 'p'};
    
    % Assuming you have a similar way to calculate mfp_results for these new datasets
    mfp_results = zeros(1, length(thicknesses)); % Initialize mfp_results for new datasets
    
    % You would need a similar loop or method to fill mfp_results based on the new data
    % Loop through each dataset to calculate mfp results
    for i = 1:length(thicknesses)
        if i == 1
            X = smooth(d85t1,30); % Assuming d85t1 is your dataset
        elseif i == 2
            X = smooth(b65t1,30); % Assuming b65t1 is your dataset
        elseif i == 3
            X = smooth(c65t1,30); % Assuming c65t1 is your dataset
        elseif i == 4
            X = smooth(e85t1,30); % Assuming e85t1 is your dataset
        elseif i == 5
            X = smooth(e65t1,30); % Assuming e65t1 is your dataset
        elseif i == 6
            X = smooth(f65t1,30); % Assuming f65t1 is your dataset
        elseif i == 7
            X = smooth(a85t1,30); % Assuming a85t1 is your dataset
        elseif i == 8
            X = smooth(b85t1,30); % Assuming b85t1 is your dataset
        else % i == 9
            X = smooth(c85t1,30); % Assuming c85t1 is your dataset
        end
        
        % Apply mfp function (assuming it's defined elsewhere)
        indices = wl >= 400 & wl <= 800;
        meanX = mean(X(indices));
        mfp_results(i) = mfp(meanX, thicknesses(i), ze(i)); % Here you'd use your actual mfp calculation function
    end    

    % Define the labels for each dataset, corresponding to the new legend entries
    datasetNames = {'85D T1', '65B T1', '65C T1', '85E T1', '65E T1', '65F T1', '85A T1', '85B T1', '85C T1'};
    
    % Scatter plot for each dataset's mfp result
    figure;
    hold on;
    
    % Calculate positions for each dataset label on the x-axis
    xPositions = 1:length(mfp_results);
    
    for i = 1:length(mfp_results)
        scatter(xPositions(i), 1/mfp_results(i), 200, colors{i}, 'filled'); % Adjust size if necessary
    end
    
    hold off;
    
    % Adjusting the scatter plot appearance
    xlim([0, length(mfp_results)+1]);
    set(gca, 'xtick', xPositions, 'xticklabel', datasetNames); % Set custom x-axis tick labels
    xlabel('Dataset');
    ylabel('Scattering efficiency');
    legend(datasetNames, 'FontSize', 12, 'Location', 'northeast'); 


end

end 

%% GLY
if GLY
    
%% Graph1: Including all T1 samples for GLY
if graph1
    
    % Create two columns of data
    A = wl;
    B = smooth(b86t1,30);
    % Create a table with the data and variable names
    T = table(A, B,'VariableNames',{'A', 'B'});
    % Write data to text file
    writetable(T,"86B_T1.txt")

    plot(wl,smooth(a66t1,30),'LineWidth',7,'Color','#B0D3E6'); hold on;
    %plot(wl,smooth(b66t1(:, 4),30),'--','LineWidth',7,'Color',c); hold on;
    plot(wl,smooth(c66t1,30),'LineWidth',7,'Color','#BFE79E'); hold on;
    plot(wl,smooth(d66t1,30),'LineWidth',7,'Color','#65D161'); hold on;
    plot(wl,smooth(a86t1,30),'LineWidth',7,'Color','#FEA3A2'); hold on;
    plot(wl,smooth(b86t1,30),'LineWidth',7,'Color','#E4635D'); hold on;

    %legend({'66A T1, t=11µm','66C T1, t=3µm','66D T1, t=5µm','86A T1, t=14µm','86B T1, t=15µm'},'FontSize',25,'Location','northeast')
    legend({'[Ca^2^+] = 1.25 mM, [Gly] = 100 mM','[Ca^2^+] = 5 mM, [Gly] = 100 mM','[Ca^2^+] = 5 mM, [Gly] = 200 mM','[Ca^2^+] = 10 mM, [Gly] = 100 mM','[Ca^2^+] = 10 mM, [Gly] = 200 mM'},'FontSize',15,'Location','northeast')

end

%% Graph2: Including all T2 samples for GLY
if graph2
    
    % Create two columns of data
    A = wl;
    B = smooth(b86t2,30);
    % Create a table with the data and variable names
    T = table(A, B,'VariableNames',{'A', 'B'});
    % Write data to text file
    writetable(T,"86B_T2.txt")

    %plot(wl,smooth(a66t2,30),'--','LineWidth',7,'Color',b); hold on;
    plot(wl,smooth(b66t2,30),'LineWidth',7,'Color','#59A7D8'); hold on;
    %plot(wl,smooth(c66t2,30),'LineWidth',7,'Color',b); hold on;
    plot(wl,smooth(d66t2,30),'LineWidth',7,'Color','#65D161'); hold on;
    plot(wl,smooth(a86t2,30),'LineWidth',7,'Color','#FEA3A2'); hold on;
    plot(wl,smooth(b86t2,30),'LineWidth',7,'Color','#E4635D'); hold on;

    %legend({'66A t2, t=25µm','66B t2, t=8.5µm','66D t2, t=4µm','86A t2, t=4µm','86B t2, t=10µm'},'FontSize',25,'Location','northeast')
    legend({'[Ca^2^+] = 1.25 mM, [Gly] = 200 mM','[Ca^2^+] = 5 mM, [Gly] = 200 mM','[Ca^2^+] = 10 mM, [Gly] = 100 mM','[Ca^2^+] = 10 mM, [Gly] = 200 mM'},'FontSize',15,'Location','northeast')

end

%% Graph3: Including all T1 1.25mM Ca samples samples for ASP
if graph3

    plot(wl,smooth(a66t1,30),'--','LineWidth',7,'Color',b); hold on;
    plot(wl,smooth(b66t1(:, 4),30),'--','LineWidth',7,'Color',c); hold on;
    legend({'66A T1, t=11µm','66B T1, t=4µm','66C T1, t=3µm','66D T1, t=5µm','86A T1, t=14µm','86B T1, t=15µm'},'FontSize',25,'Location','northeast')


end

%% Graph4: Including all T2 1.25mM Ca samples samples for ASP
if graph4

    plot(wl,smooth(a66t2,30),'--','LineWidth',7,'Color',b); hold on;
    plot(wl,smooth(b66t2,30),'--','LineWidth',7,'Color',c); hold on;
    legend({'66A t2, t=?µm','66B t2, t=8.5µm','66C t2, t=?µm','66D t2, t=4µm','86A t2, t=4µm','86B t2, t=10µm'},'FontSize',25,'Location','northeast')


end

if graph5
    thickness = [11.3, 4.4, 2.9, 4.96,13.6,   15.3];
    ze = [2.6, 2.6, 2.4, 2.5, 2.3,  2.4];

    plot(wl,smooth(a66t1,30),'--','LineWidth',7,'Color',b); hold on;
    plot(wl,smooth(b66t1(:, 4),30),'--','LineWidth',7,'Color',c); hold on;
    plot(wl,smooth(c66t1,30),'LineWidth',7,'Color',b); hold on;
    plot(wl,smooth(d66t1,30),'LineWidth',7,'Color',c); hold on;
    plot(wl,smooth(a86t1,30),'LineWidth',7,'Color',g); hold on;
    plot(wl,smooth(b86t1,30),'LineWidth',7,'Color',r); hold on;
    legend({'66A T1, t=11µm','66B T1, t=4µm','66C T1, t=3µm','66D T1, t=5µm','86A T1, t=14µm','86B T1, t=15µm'},'FontSize',25,'Location','northeast')

    % Initialize mfp_results array for the new datasets
    mfp_results = zeros(1, length(thickness));
    
    % Assuming the existence of a function mfp(meanX, thickness, ze) for calculating mfp
    % Loop through each dataset
    for i = 1:length(thickness)
        switch i
            case 1
                X = smooth(a66t1, 30);
            case 2
                X = smooth(b66t1(:, 4), 30);
            case 3
                X = smooth(c66t1, 30);
            case 4
                X = smooth(d66t1, 30);
            case 5
                X = smooth(a86t1, 30);
            case 6
                X = smooth(b86t1, 30);
        end
    
        % Selecting the subset of values where wl is between 400 and 800
        indices = wl >= 400 & wl <= 800;
    
        % Calculate mean only for the subset of X corresponding to those indices
        meanX = mean(X(indices));
    
        % Store the calculated mfp result
        mfp_results(i) = mfp(meanX, thickness(i), ze(i)); % Use your actual mfp calculation function
    end
    
    % Plotting the scatter plot for mfp results
    figure; hold on;
    colors = {'b', 'c', 'b', 'c', 'g', 'r'}; % Colors for each dataset
    datasetNames = {'66A T1', '66B T1', '66C T1', '66D T1', '86A T1', '86B T1'}; % Dataset names
    for i = 1:length(mfp_results)
        scatter(i, 1/mfp_results(i), 100, colors{i}, 'filled'); % Plot each mfp result
    end
    set(gca, 'xtick', 1:length(datasetNames), 'xticklabel', datasetNames); % Custom x-axis tick labels
    xlabel('Dataset');
    ylabel('Inverse mfp Result');
    title('Scatter Plot of Inverse mfp Results for Each Dataset');
    legend(datasetNames, 'FontSize', 12, 'Location', 'northeast');
    hold off;
    

end

end 

%% VES
if VES
    
%% Graph1: Including all samples for VES
if graph1

   ze_values(1) = 2.47;
ze_values(6) = 1.55;
tt = linspace(0.1,35,100);
% Calculate R for l = 1
R_l1 = tt ./ (2 * ze_values(6) + tt);
% Calculate R for l = 1.5
R_l1_5 = tt ./ (3 * ze_values(1) + tt);
% Plot the curves
plot(tt, R_l1, 'k', 'DisplayName','l = 1');
plot(tt, R_l1_5, 'k', 'DisplayName','l = 1.5');
% Shade the area between the curves
fill([tt, fliplr(tt)], [R_l1, fliplr(R_l1_5)], 'r', 'FaceAlpha', 0.2, 'EdgeColor', 'none');
 
    plot(wl,smooth(a82t1,30),'LineWidth',7,'Color','#F9D49D'); hold on;
    plot(wl,smooth(b82t1,30),'LineWidth',7,'Color','#FD9C33'); hold on;
    plot(wl,smooth(c82t1,30),'LineWidth',7,'Color','#FE3C00'); hold on;
    %%plot(wl,smooth(d82t1,30),'--','LineWidth',7,'Color',r); hold on;

    %legend({'82 A, t=13µm','82 B, t=3µm','82 C, t=7µm','82 D, t=7.5µm' },'FontSize',25,'Location','northwest')
    legend({'0.05 wt.% Ves','0.1 wt.% Ves','0.25 wt.% Ves'},'FontSize',25,'Location','northwest')

end

%% Graph2: Including all Ves samples excluding 0.5% Ves 82D
if graph2

    plot(wl,smooth(a82t1,30),'LineWidth',7,'Color',b); hold on;
    plot(wl,smooth(b82t1,30),'LineWidth',7,'Color',c); hold on;
    plot(wl,smooth(c82t1,30),'LineWidth',7,'Color',g); hold on;


    legend({'82 A, t=13µm','82 B, t=3µm','82 C, t=7µm'},'FontSize',25,'Location','northwest')


end

if graph5
    thickness = [13,3.5,7.1]
    ze = [1.89,2.3,1.9]

    plot(wl,smooth(a82t1,30),'LineWidth',7,'Color',b); hold on;
    plot(wl,smooth(b82t1,30),'LineWidth',7,'Color',c); hold on;
    plot(wl,smooth(c82t1,30),'LineWidth',7,'Color',g); hold on;

    legend({'82 A, t=13µm','82 B, t=3µm','82 C, t=7µm'},'FontSize',25,'Location','northwest')

    % Initialize mfp_results array for the new datasets
    mfp_results = zeros(1, length(thickness));
    
    % Assuming the existence of a function mfp(meanX, thickness, ze) for calculating mfp
    % Loop through each dataset
    for i = 1:length(thickness)
        switch i
            case 1
                X = smooth(a82t1, 30);
            case 2
                X = smooth(b82t1, 30);
            case 3
                X = smooth(c82t1, 30);
        end
    
        % Selecting the subset of values where wl is between 400 and 800
        indices = wl >= 400 & wl <= 800;
    
        % Calculate mean only for the subset of X corresponding to those indices
        meanX = mean(X(indices));
    
        % Store the calculated mfp result
        mfp_results(i) = mfp(meanX, thickness(i), ze(i)); % Use your actual mfp calculation function
    end
    
    % Plotting the scatter plot for mfp results
    figure; hold on;
    colors = {'b', 'c', 'g'}; % Colors for each dataset
    datasetNames = {'82 A', '82 B', '82 C'}; % Dataset names
    for i = 1:length(mfp_results)
        scatter(i, 1/mfp_results(i), 100, colors{i}, 'filled'); % Plot each mfp result
    end
    set(gca, 'xtick', 1:length(datasetNames), 'xticklabel', datasetNames); % Custom x-axis tick labels
    xlabel('Dataset');
    ylabel('Inverse mfp Result');
    title('Scatter Plot of Inverse mfp Results for Each Dataset');
    legend(datasetNames, 'FontSize', 12, 'Location', 'northwest');
    hold off;


end
end 


%% Unheated

if unheated
    
%% Graph1: Including all T1 samples for unheated
if graph1
    
    
    
    buh61=mean(spectraUnheated(11:15,:));
    %duh61=mean(spectraUnheated(26:30,:));
    %fuh61=mean(spectraUnheated(41:45,:));

    % Create two columns of data
    A = wl;
    B = smooth(buh61,30);
    % Create a table with the data and variable names
    T = table(A, B,'VariableNames',{'A', 'B'});
    % Write data to text file
    writetable(T,'unheatPS.txt')

    %xal1=(spectraUnheated(11,:));
    %xal2=(spectraUnheated(12,:));
    %xal3=(spectraUnheated(13,:));
    %xal4=(spectraUnheated(14,:));
    %xal5=(spectraUnheated(15,:));

    %%plot(wl,smooth(buh61,30),'LineWidth',7,'Color',g); hold on;
    %plot(wl,smooth(duh61,30),'--','LineWidth',7,'Color',g); hold on;
    %plot(wl,smooth(fuh61,30),'--','LineWidth',7,'Color',g); hold on;

    %plot(wl,smooth(xal1,30),'--','LineWidth',7,'Color',g); hold on;
    %plot(wl,smooth(xal2,30),'--','LineWidth',7,'Color',g); hold on;
    %plot(wl,smooth(xal3,30),'--','LineWidth',7,'Color',g); hold on;
    %plot(wl,smooth(xal4,30),'--','LineWidth',7,'Color',g); hold on;
    %plot(wl,smooth(xal5,30),'--','LineWidth',7,'Color',g); hold on;

    %legend({'61_buh, t=6µm','61_duh, t=3µm','61_fuh, t=4µm'},'FontSize',25,'Location','northeast')
end
end

%% Other whites
if otherWhites
    plot(wl2,smooth(1-AdvFunctMaterials,20),'LineWidth',7,'Color','k');xlim([400 800]);hold on;
    plot(wl,smooth(AdvMater,20),'--','LineWidth',7,'Color',gray);
    legend({'AdvFunctMater, t=15µm','AdvMater, t=10µm'},'FontSize',25,'Location','southeast')
end
%% Heated calcite

if heatedCalcite
    
%% Graph1: 
if graph1
    
    % Create two columns of data
    A = wl;
    B = smooth(a63,30);
    % Create a table with the data and variable names
    T = table(A, B,'VariableNames',{'A', 'B'});
    % Write data to text file
    writetable(T,'63A_T1.txt')

    %plot(wl,smooth(a63,30),'LineWidth',7,'Color',b); hold on;
        
    %legend({'63A Heated Calcite'},'FontSize',25,'Location','northeast')
end
end

%% Comparison Different Additives
if comparisonAdditives
    %%% old data selected for best
    if 0
        plot(wl,smooth( B,30),'LineWidth',7,'Color',r); hold on;
        plot(wl, smooth(d85t1,30),'LineWidth',7,'Color',g); hold on;
    
        plot(wl,smooth(a66t2,30),'LineWidth',7,'Color',c); hold on;
        plot(wl, smooth(a82t1,30),'LineWidth',7,'Color',b); hold on;
    
        if otherWhites==0
        legend({'[PS]=1%, [Ca2+]=1.25mM , 61B T1, t=7µm', '[Asp]=10, [Ca2+]=1.25mM, 85 D, t=9µm','[Gly]= 100, [Ca2+]=1.25mM, 66 A T2, t=25µm', '[Ves]=0.05%, [Ca2+]=10mM, 82 A T1, t=13µm'},'FontSize',25,'Location','northeast')
        end
        legend({'AdvFunctMater, t=15µm','AdvMater, t=10µm','[PS]=1%, [Ca2+]=1.25mM , 61B T1, t=7µm', '[Asp]=10, [Ca2+]=1.25mM, 85 D, t=9µm','[Gly]= 100, [Ca2+]=1.25mM, 66 A T2, t=25µm', '[Ves]=0.05%, [Ca2+]=10mM, 82 A T1, t=13µm'},'FontSize',25,'Location','northeast')
    end

    if 1
        plot(wl,smooth(d74(:,4)*(2*d74(1200,3)/d74(1200,4)),30),'LineWidth',7,'Color',r); hold on;
        plot(wl, smooth(d85t1,30),'LineWidth',7,'Color',g); hold on;
    
        plot(wl,smooth(b66t1(:, 4),30),'LineWidth',7,'Color',c); hold on;
        plot(wl, smooth(a82t1,30),'LineWidth',7,'Color',b); hold on;
    
        if otherWhites==0
        legend({'[PS]=3%, [Ca2+]=1.25mM , 74D, t=2.5µm', '[Asp]=10, [Ca2+]=1.25mM, 85D, t=9µm','[Gly]= 200, [Ca2+]=1.25mM, 66B, t=4.5µm', '[Ves]=0.05%, [Ca2+]=10mM, 82A, t=13µm'},'FontSize',25,'Location','northeast')
        end
        legend({'AdvFunctMater, t=15µm','AdvMater, t=10µm','[PS]=3%, [Ca2+]=1.25mM , 74D, t=2.5µm', '[Asp]=10, [Ca2+]=1.25mM, 85D, t=9µm','[Gly]= 200, [Ca2+]=1.25mM, 66B, t=4.5µm', '[Ves]=0.05%, [Ca2+]=10mM, 82A, t=13µm'},'FontSize',25,'Location','northeast')

            if 1
            thickness = [15,10,2.5,9,4.5,13];
            ze = [1.5,1.5,2.1,2.5,2.6,1.9];
    
            % Initialize mfp_results array
            mfp_results = zeros(1, length(thickness));
            
            % Loop through each dataset
            datasets = {'advF','advM','d74_modified', 'd85t1', 'b66t1_double_smoothed', 'a82t1'}; % Example dataset identifiers
            for i = 1:length(datasets)
                switch datasets{i}
                    case 'advF'
                        X = smooth(1-AdvFunctMaterials,20);
                        %X = 1 - 0.16;
                     case 'advM'
                        X = smooth(AdvMater,20);                   
                    case 'd74_modified'
                        X = smooth(d74(:,4) * (2 * d74(1200,3) / d74(1200,4)), 30);
                    case 'd85t1'
                        X = smooth(d85t1, 30);
                    case 'b66t1_double_smoothed'
                        X = smooth(b66t1(:, 4), 30);
                    case 'a82t1'
                        X = smooth(a82t1, 30);
                end
            
                % Selecting the subset of values where wl is between 400 and 800
                indices = wl >= 400 & wl <= 800;
            
                % Calculate mean only for the subset of X corresponding to those indices
                meanX = mean(X(indices))
            
                % Store the calculated mfp result
                mfp_results(i) = mfp(meanX, thickness(i), ze(i)) % Use your actual mfp calculation function
            end
            
            graph5=1
            mfp_results(1) = 1
            % Plotting the scatter plot for mfp results
            figure; hold on;
            colors = {'k',gray,'r', 'g', 'c', 'b'}; % Colors corresponding to the plot colors
            datasetNames = {'advFun','advMat','PS', 'Asp', 'Gly', 'Ves'}; % Custom dataset names
            for i = 1:length(mfp_results)
                scatter(i, mfp_results(i), 250, colors{i}, 'filled'); % Plot each mfp result
            end
            set(gca, 'xtick', 1:length(datasetNames), 'xticklabel', datasetNames); % Custom x-axis tick labels
            xlabel('Dataset');
            %ylabel('Scattering Efficiency [1/lt]');
            ylabel('Transport mean free path');
            title('Comparison among best samples');
            legend(datasetNames, 'FontSize', 25, 'Location', 'northwest');
            hold off;
            xlim([0, length(datasetNames) + 1]);
            %ylim ([0 3.5]);
            end


    end


end



%% Figure properties
ax = gca;
set(gca,'XMinorTick','off','YMinorTick','off')
ax.Box = 'on';
ax.LineWidth = 7;
set(gca,'FontSize',50);
set(findall(gcf,'type','text'),'FontSize',50);
%set(gca,'color','none');
%set(gcf,'color','none');

pbaspect([1 1 1]);
frame_h = get(handle(gcf),'JavaFrame');
set(frame_h,'Maximized',1);
if graph5 == 0
    ylim([0 1]);
    xlim([400 800]);
    xlabel('Wavelength (nm)');
    ylabel('Reflectance');
end
legend boxoff;



% legend('Location','southeast')
set(legend,'color','none');
% title ('Model, ff=0.55, aligned particles','FontSize', 25)

%% MFP calculator
function l = mfp(R,t,ze)

l = t*(1-R)./(2*ze*(R));

end