% Matlab script DoFTIR - Authored February 8, 2024 by David Kazmer
% Performs analysis of Sixtus Nzeh's recycled PP films by two FTIR methods
% Note that there are two excellent master curves fom cited literature
% The set of FTIR data and master curve are switched by setting bUseATR

% Set Boolean values for execution options
bUseATR=true;
bPause=false;
bMarkAllLabels=false;

% Prepare figures for plotting
for iPlot=1:3
    figure(iPlot);clf;hold on;
end

% Define string styles for plotting
sStyle{1}='k-';
sStyle{2}='k--';
sStyle{3}='k-.';
sStyle{4}='k-';
sStyle{5}='r-';
sStyle{6}='r--';
sStyle{7}='r-.';
sStyle{8}='r:';

% Define the wavenumbers of interest and the range
wavenumbers_of_interest = [720, 1168, 719, 1376];
range = 0.5;

% Define string styles for critical peaks
sPeakStyles{1}='ko';sPeakLabels{1}='PP @ 720';
sPeakStyles{2}='rd';sPeakLabels{2}='PE @ 1168';
sPeakStyles{3}='ko';sPeakLabels{3}='PP @ 719';
sPeakStyles{4}='rd';sPeakLabels{4}='PE @ 1376';

for iMat=1:2
    for iTest=1:4
        % Load spectal data
        if bUseATR
            sFile = ['PP' num2str(iMat) ' ' num2str(iTest) ' ATR'];
        else
            sFile = ['PP' num2str(iMat) ' ' num2str(iTest)];
        end
        if iTest<4
            data = readmatrix([sFile '.CSV']);
            wavenumber = data(:, 1);
            absorbance = data(:, 2)-min(data(:,2));
            all_waves(:,iTest)=wavenumber;
            all_absorbances(:,iTest)=absorbance;
            absorbance = smooth(absorbance);
        else
            wavenumber=(mean(all_waves'))';
            absorbance=(mean(all_absorbances'))';
        end
        
        range_of_PP = [2950, 2915, 2850, 1450, 1375, 1160, 998, 840, 810];
        % 2950 cm⁻¹: Asymmetric CH₃ stretching
        % 2915 cm⁻¹: Asymmetric CH₂ stretching
        % 2850 cm⁻¹: Symmetric CH₂ stretching
        % 1450 cm⁻¹: CH₂ bending
        % 1375 cm⁻¹: CH₃ bending
        % 1160 cm⁻¹, 998 cm⁻¹: CH rocking and twisting
        % 840 cm⁻¹, 810 cm⁻¹: CH out-of-plane bending
        
        range_of_PE = [2914, 2847, 1470, 718];
        % 2914 cm⁻¹ and 2847 cm⁻¹: C-H stretching
        % 1470 cm⁻¹: C-C stretching
        % 718 cm⁻¹: C-C bending
        
        [peaks, locations] = findpeaks(absorbance, wavenumber, 'MinPeakProminence', 0.05);
        identified_polymers = cell(size(peaks));
        
        for i = 1:length(peaks)
            if any(abs(locations(i) - range_of_PP) <= 10)  % Allow for slight variations
                identified_polymers{i} = 'PP';
            elseif any(abs(locations(i) - range_of_PE) <= 10)  % Allow for slight variations
                identified_polymers{i} = 'PE';
            elseif bMarkAllLabels
                identified_polymers{i} = 'Unknown';
            end
        end

        % Calculate curve number for style and make the plot
        iCurve=4*(iMat-1)+iTest;

        % Loop over specific wavenumbers of interest
        for i = 1:length(wavenumbers_of_interest)
            % Find the indices of the wavenumbers within the range of interest
            indices = find(wavenumber >= wavenumbers_of_interest(i) - range & wavenumber <= wavenumbers_of_interest(i) + range);
        
            % Calculate the average absorbance within this range
            average_absorbance(iCurve,i) = mean(absorbance(indices));

            % Print result
            fprintf('Curve %d: average absorbance near a wavenumber of %d cm^-1 (±%d cm^-1) is %.2f.\n', iCurve,wavenumbers_of_interest(i), range, average_absorbance(iCurve,i));
        end

        % Plot three sets of figures: (1) single curves, (2) all curves for
        % one material, and (3) final master curves for PP1 and PP2
        for iPlot=1:3
            if (iPlot<3)||(iTest==4)
                figure(iPlot);
                % Clear plot for each figure(1) and first figure(2)
                if (iPlot==1)||((iPlot==2)&&(iTest==1)),clf;hold on;end

                if iTest==4
                    plot(wavenumber, absorbance,sStyle{iCurve},'linewidth',2,'displayname',['PP' num2str(iMat) ' average, unfiltered']);
                elseif (iPlot==2)
                    plot(wavenumber, absorbance,sStyle{iCurve},'displayname',sFile);
                else
                    plot(wavenumber, absorbance,sStyle{iCurve});
                end
                if (iPlot==1)||(iTest==4)
                    if bMarkAllLabels
                        scatter(locations, peaks, 'r','displayname','Identified peaks');
                        text(locations, peaks, identified_polymers);
                    else
                        if iMat==1
                            scatter(wavenumbers_of_interest,average_absorbance(iCurve,:),'ko',...
                           'MarkerFaceColor','y','Displayname','Idenitified peaks for PP1');
                            text(wavenumbers_of_interest,average_absorbance(iCurve,:),sPeakLabels);
                        else
                            scatter(wavenumbers_of_interest,average_absorbance(iCurve,:),'kh',...
                           'MarkerFaceColor','g','Displayname','Idenitified peaks for PP2');
                            text(wavenumbers_of_interest,average_absorbance(iCurve,:),sPeakLabels);
                        end
                    end
                end
                xlabel('Wavenumber, cm^-^1');
                ylabel('Absorbance');
                if iPlot==1
                    title(['Spectral Data with Identified Polymers for ' sFile]);
                else
                    legend('location','northwest')
                end
                grid on;xlim([400 2000]);
               
            end
        end
        if bPause;
            pause
        end
    end

    % Plot master curve and points on figure(4)
    figure(4);
    if iMat==1
        x=0:100;
        if bUseATR
            % Plot master curve from Figure 2 of Xiao & Ernst, ATR Infrared Spectrometer
            % Field Report, PerkinElmer, 2012
            subplot(1,2,1);hold on;
            A_Master=x;                         % Impressively simple!
            %plot(x,A_Master,'k-','linewidth',2,'displayname',{'Master curve';'Xiao & Ernst w/ATR'});grid on;
            plot(x,0.01*A_Master,'k-','linewidth',2,'displayname','Xiao & Ernst ATR');grid on;
            xlabel('PP content, % wt');
            ylabel('1-A(719)/A(1376)')
        else
            % Plot master curve from Fig. 1 of Camacho and Karlsson, PES 2001 (very nice)
            subplot(1,2,2);hold on
            A_Master=4e-5*x.^2+0.0058*x;
            plot(x,A_Master,'k-','linewidth',2,'displayname','Camacho & Karlsson');grid on;
            xlabel('PP content, % wt');
            ylabel('A(1168)/[A(1168)+A(720)]')
        end
    end

    if bUseATR
        % Calculate attenuation fraction  per Xiao & Ernst
        A_PP(iMat)=1-average_absorbance(iCurve,3)/(average_absorbance(iCurve,4));
    
        % Calculate percent hdpe based on inverse of master curve
        c_PP(iMat)=100*A_PP(iMat);
    else
        % Calculate attenuation fraction  per Camacho and Karlsson
        A_PP(iMat)=average_absorbance(iCurve,2)/(average_absorbance(iCurve,2)+average_absorbance(iCurve,1));
    
        % Calculate percent hdpe based on inverse of master curve
        c_PP(iMat)=5/2*(sqrt(841 + 4000 * A_PP(iMat)) - 29);
    end
    
    % Plot and label the result
    if bUseATR
        sPPLabel=sprintf('PP%d - E_X_E[%.1f%% PP]',iMat,c_PP(iMat));
    else
        sPPLabel=sprintf('PP%d - E_C_K[%.1f%% PP]',iMat,c_PP(iMat));
    end

    if iMat==1
        plot(c_PP(iMat),A_PP(iMat),'k^','MarkerFaceColor','w','MarkerSize',6,'DisplayName',sPPLabel)
    else
        plot(c_PP(iMat),A_PP(iMat),'r^','MarkerFaceColor','k','MarkerSize',6,'DisplayName',sPPLabel)
    end

end

figure(4)
% Add larsen
if ~bUseATR
    z=0.01*(0:0.1:100);
    plot(100-(120.9*z.^2-217.97*z+100),z,'k-.','linewidth',2,'displayname','Larsen et al')
    plot(98.0,0.862,'kv','MarkerFaceColor','w','MarkerSize',6,'DisplayName','PP1 - E_L[98.0% PP]')
    plot(69.3,0.398,'rv','MarkerFaceColor','k','MarkerSize',6,'DisplayName','PP2 - E_L[69.3% PP]')
end

legend('location','northwest')