% PlotViscosity_vPP_vPE - Matlab script written by David Kazmer to read & plot viscosity data
% December 18, 2024
% Still need to fit constitutive model to include Arrhenius or WLF temperature dependence
warning off

% Define files and variables
sFile{1}='vHDPE L-D 10.xlsx';
sFile{2}='vHDPE L-D 30.xlsx';
sFile{3}='vPP profax 6523 FDM grade L-D 10.xlsx';
sFile{4}='vPP profax 6523 FDM grade L-D 30.xlsx';
sMaterial{1}='vPE';sColor{1}='k';sLD{1}='v';
sMaterial{2}='vPP';sColor{2}='b';sLD{2}='s';
i_LD=[1 2 1 2];;
i_LD_Actual=[10 30 10 30];
i_Mat=[1 1 2 2];
i_Temps=[230 230 230 230];   % Note: 210 and 250 also characterized

% Define global variables and options for fitting
global bFitLog iModel_G1G2 GDot Eta_Obs G1_Obs G2_Obs 
bFitLog=false;
iModel_G1G2=3;

% Close figures for replotting
close all

for i=1:4
    % Get Data
    Data_Table=readtable(sFile{i},'VariableNamingRule','preserve');
    Shear_Rates=table2array(Data_Table(:,27));
    Viscosities=table2array(Data_Table(:,29));
    Shear_Stresses=table2array(Data_Table(:,26));
    Temps=table2array(Data_Table(:,9));

    % Plot data for each trace
    figure(i);hold on;
    k=find(Temps==i_Temps(i));
    % Remove known outliers for 230 C and 10:1 L:D
    if i==1,k=k(3:end),end
    plot(Shear_Rates(k),Viscosities(k),'--','Marker',sLD{i_LD(i)},'Color',sColor{i_Mat(i)}, ...
        'DisplayName',[sMaterial{i_Mat(i)} ': ' num2str(i_LD_Actual(i)) ':1 L:D']);
    xlabel('Apparent Shear Rate, s^-^1');ylabel('Apparent Viscosity, Pa.s')
    set(gca,'xscale','log','yscale','log');grid on;
    legend('location','NorthEast');

    % Set constraints
    x0=[Eta0 n  kW];
    lb=[  10 0  0];
    ub=[ 1e6 1  1];
    options = optimset ('fmincon') ;
    options = optimset (options,...
        'LargeScale','off',...
        'Algorithm','interior-point',...
        'MaxFunEvals',100000,...
        'MaxIter',5000,...
        'Display','iter'); 

    % Fit the model and store statistics
    GDot=Shear_Rates(k);Eta_Obs=Viscosities(k);
    Coeffs=fmincon(@sumErrorViscosity,x0,[],[],[],[],lb,ub,[],options);
    Eta=Williamson(GDot,Coeffs);
    [B,BINT,R,RINT,STATS]=regress(log(Eta_Obs),[log(Eta) ones(size(Eta))]);
    R2=STATS(1);p=STATS(3);m=B(1);b=B(2);
       
    % Plot individual viscosity models
    figure(i)
    Shear_Rates_Model=10.^(1:0.25:4);
    Eta_Model=Williamson(Shear_Rates_Model,Coeffs);
    sWilliamson=sprintf('%s_0=%.0f, n=%.3f, k=%.4f; R^2=%.4f', ...
       [sMaterial{i_Mat(i)} ': ' eta_char],Coeffs(1),Coeffs(2),Coeffs(3),R2);

    plot(Shear_Rates_Model,Eta_Model,':','Color',sColor{i_Mat(i)},'DisplayName',sWilliamson);
    legend('location','SouthWest');
    set(gca,'XScale','log','YScale','log');grid on;
    ylim([1e1 1e3])

    % Aggregate data, fit, and plot combined L:D comparative viscosities
    figure(5);hold on;
    plot(Shear_Rates(k),Viscosities(k),sLD{i_LD(i)},'Color',sColor{i_Mat(i)}, ...
        'DisplayName',[sMaterial{i_Mat(i)} ': 230 C, ' num2str(i_LD_Actual(i))  ':1 L:D']);
    xlabel('Apparent Shear Rate, s^-^1');ylabel('Apparent Viscosity, Pa.s')
    set(gca,'xscale','log','yscale','log')
    legend('Location','southwest');
    grid on;
    if mod(i,2)
        GDot=Shear_Rates(k);Eta_Obs=Viscosities(k);
        Sigma10=Shear_Stresses(k(end-7:end));
    else
        Sigma30=Shear_Stresses(k(end-7:end));
        GDot=[GDot;Shear_Rates(k)];Eta_Obs=[Eta_Obs;Viscosities(k)];
        % Fit the model for combined data and store statistics
        Coeffs=fmincon(@sumErrorViscosity,x0,[],[],[],[],lb,ub,[],options);
        Eta=Williamson(GDot,Coeffs);
        [B,BINT,R,RINT,STATS]=regress(log(Eta_Obs),[log(Eta) ones(size(Eta))]);
        R2=STATS(1);p=STATS(3);m=B(1);b=B(2);
        % Plot viscosity data for each temperature
        Shear_Rates_Model=10.^(1:0.25:4);
        Eta_Model=Williamson(Shear_Rates_Model,Coeffs);
        eta_char = char(951);
        sWilliamson=sprintf('%s_0=%.0f, n=%.3f, k=%.4f; R^2=%.4f', ...
            [sMaterial{i_Mat(i)} ': ' eta_char],Coeffs(1),Coeffs(2),Coeffs(3),R2);
        plot(Shear_Rates_Model,Eta_Model,':','Color',sColor{i_Mat(i)},'DisplayName',sWilliamson);
        legend('location','SouthWest');
        set(gca,'XScale','log','YScale','log');grid on;
        ylim([1e1 1e3])
        % Make a plot of entrance/exit effects
        figure(6);hold on;
        Gamma=mean([Sigma10 Sigma30],2);
        CA=Gamma./(Sigma10-Sigma30);
        CA_mean(i/2)=mean(CA(1:4));
        CA_std(i/2)=std(CA(1:4));
        plot(Shear_Rates(k(end-7:end)),CA,'--','Marker',sLD{i_LD(i)},'Color',sColor{i_Mat(i)}, ...
        'DisplayName',sMaterial{i_Mat(i)})
        xlabel('Apparent Shear Rate, s^-^1');ylabel('Ratio of capillary to entrance/exit forces')
        legend('location','NorthWest');
        set(gca,'xscale','log','yscale','log');grid on;
    end
end
CA_mean
CA_std
figure(5);print('Fig-U-Viscosity.tiff','-r600','-dtiff')
figure(6);print('Fig-V-Ce.tiff','-r600','-dtiff')

function eta = Williamson(gdot,C)
    eta=C(1)./(1+(C(3)*abs(gdot)).^(1-C(2))); 
end

function fSumError=sumErrorViscosity(X)
global bFitLog GDot Eta_Obs
    Eta=Williamson(GDot,X);
    if bFitLog
        fSumError=sum((log(Eta)-log(Eta_Obs)).^2);
    else
        fSumError=sum((Eta-Eta_Obs).^2);
    end
end

function G1est = CalcG1(gdot,C)
global iModel_G1G2
    if iModel_G1G2==1
        G1est=C(1)*gdot.^2*C(2)^2./(1+C(2)^2*gdot.^2); 
    elseif iModel_G1G2==2
        G1est=C(1)*gdot.*((C(2)-C(3))*gdot.^2)./(1+C(2)^2*gdot.^2); 
    elseif iModel_G1G2==3
        G1est=C(2)*C(3)*gdot.^2*(C(4)^2-3)/C(4)^2./(1+C(3)^2*gdot.^2*(C(4)^2-3)^2/C(4)^4);
    end
end

function G2est = CalcG2(gdot,C)
global iModel_G1G2
    if iModel_G1G2==1
        G2est=C(3)*gdot+C(1).*gdot*C(2)./(1+C(2)^2*gdot.^2); 
    elseif iModel_G1G2==2
        G2est=C(1)*gdot.*(1+C(2)*C(3)*gdot.^2)./(1+C(2)^2*gdot.^2); 
    elseif iModel_G1G2==3    % FENE-P model from Purnode & Crochet, 1998
        G2est=C(1)*gdot+C(2)*gdot./(1+C(3)^2*gdot.^2*(C(4)^2-3)^2/C(4)^4); 
    end
end

function fSumError=sumErrorG1G2(X)
global bFitLog bFit_Alt_G1G2 GDot G1_Obs G2_Obs
    G1_Est=CalcG1(GDot,X);
    G2_Est=CalcG2(GDot,X);
    if bFitLog
        fSumError=sum((log(G1_Est)-log(G1_Obs)).^2+(log(G2_Est)-log(G2_Obs)).^2);
    else
        fSumError=sum((G1_Est-G1_Obs).^2+(G2_Est-G2_Obs).^2);
    end
end

function Etaest = CalcEta_FENEP(gdot,C)
    N=(length(C)-1)/3;
    Etaest=C(1);
    for i=1:N
        L=C(1+i*3);
        Etaest=Etaest+C(i*3-2)*(L^6./(4*C(i*3-1)^2*gdot.^2*(L^2-3)^2)).^(1/3); 
    end
end

function fSumError=sumErrorEta_FENEP(X)
global bFitLog GDot Eta_Obs
    Eta=CalcEta_FENEP(GDot,X);
    if bFitLog
        fSumError=sum((log(Eta)-log(Eta_Obs)).^2);
    else
        fSumError=sum((Eta-Eta_Obs).^2);
    end
end