% FitWilliamson - Matlab script written by David Kazmer to read & fit viscosity data
warning off

% Get files names and set styles
sFileCapillary{1}='PP1 Capillary.csv';
sFileCapillary{2}='PP2 Capillary.csv';
sFileParallel{1}='PP1 Parallel Plate.txt';
sFileParallel{2}='PP2 Parallel Plate.txt';
sStyle{1}='kv';sFilled{1}='w';sLine{1}='k--';sWeight{1}=2;
sStyle{2}='ko';sFilled{2}='w';sLine{2}='k-.';sWeight{2}=2;
sStyle{3}='rv';sFilled{3}='k';sLine{3}='r--';sWeight{3}=1;
sStyle{4}='ro';sFilled{4}='k';sLine{4}='r-.';sWeight{4}=1;
global bFitLog iModel_G1G2 GDot Eta_Obs G1_Obs G2_Obs 
bFitLog=true;
iModel_G1G2=3;

% Define temperature of interest
MyTemp=230;
figure(1);clf;hold on;

% Loop on each material
for iFile=1:2
    % Loop on capillary and (2&3) parallel plate data
    for iType=1:2
        % Define index iStyle for style and legend
        iStyle=(iFile-1)*2+iType;
        if iType==1
            sFileName=sFileCapillary{iFile};
            % Read the table, then extract columns and rows
            T=readtable(sFileName,'VariableNamingRule','preserve');
            d=table2array(T(:,[9 28 29]));
            i_MyTemp=find(d(:,1)==MyTemp);
            Eta_Obs=d(i_MyTemp,3);
            GDot=d(i_MyTemp,2);
            sXlabel='Shear rate, s^-^1';
        else
            sFileName=sFileParallel{iFile};
            T=readtable(sFileName,'VariableNamingRule','preserve');
            d=table2array(T(2:end,:));  % Remove first row w/units
            Eta_Obs=d(:,10);
            GDot=d(:,1);
            G1_Obs=1e6*d(:,7);
            G2_Obs=1e6*d(:,8);
            sXlabel='Angular frequency, rad/s';
        end
    
        if iType<3     % Fit Williamson
            % Set default Williamson coefficients (generic HDPE from CMold)
            Eta0=1000;
            n=0.5;
            k=0.1;
        
            % Set constraints
            x0=[Eta0 n  k];
            lb=[  10 0  0];
            ub=[ 1e6 1 1e3];
            options = optimset ('fmincon') ;
            options = optimset (options,...
                'LargeScale','off',...
                'Algorithm','interior-point',...
                'MaxFunEvals',100000,...
                'MaxIter',5000,...
                'Display','iter'); 
        
            % Fit the model 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
            figure(1);
            plot(GDot,Eta_Obs,sStyle{iStyle},'MarkerFaceColor',sFilled{iStyle},'DisplayName',sFileName(1:end-4));
            Shear_Rates=10.^(-1:0.25:4);
            Eta=Williamson(Shear_Rates,Coeffs);
            eta_char = char(951);
            sWilliamson=sprintf('%s_0=%.0f, n=%.3f, k=%.4f; R^2=%.4f', ...
                eta_char,Coeffs(1),Coeffs(2),Coeffs(3),R2);
            plot(Shear_Rates,Eta,sLine{iStyle},'DisplayName',sWilliamson);
            legend('location','SouthWest');
            set(gca,'XScale','log','YScale','log');grid on;
            xlabel('Apparent shear rate, 1/s || Oscillatory frequency, rad/s');
            ylabel('Shear viscosity, Pa.s || Dynamic viscosity, Pa.s');
            ylim([10 3000]);
        else
            % Set default Oldroyd coefficients (generic HDPE from CMold)
            if iModel_G1G2==1
                G=10000;
                tau=0.5;
                Eta0=100;
                x0=[  G tau Eta0];
                lb=[ 10   0    0];
                ub=[1e6  10  1e6];
            elseif iModel_G1G2==2
                Eta0=150;
                L1=0.02;
                L2=0.01;
                x0=[Eta0 L1  L2];
                lb=[  10 0  0];
                ub=[ 1e6 10 10];
            elseif iModel_G1G2==3      % FENE-P 
                EtaS=150;
                EtaG=1;
                Lambda=1;
                L=2;
                x0=[EtaS EtaG 2*Lambda L EtaG 0.2*Lambda L];
                lb=[  10 .1 0  1.74 .1 0  1.74];
                ub=[ 1e6 10000 10 2.50 10000 10 2.50];
            end

            % Set constraints
            options = optimset ('fmincon') ;
            options = optimset (options,...
                'LargeScale','off',...
                'Algorithm','interior-point',...
                'MaxFunEvals',1000000,...
                'MaxIter',50000,...
                'Display','iter'); 
        
            % Fit the model and store statistics
            Coeffs_G1G2=fmincon(@sumErrorG1G2,x0,[],[],[],[],lb,ub,[],options);
            G1_Est=CalcG1(GDot,Coeffs_G1G2);
            G2_Est=CalcG2(GDot,Coeffs_G1G2);
            [BL,BINTL,RL,RINTL,STATSL]=regress(log([G1_Obs;G2_Obs]),[log([G1_Est;G2_Est]) ones(size([GDot;GDot]))]);
            R2L=STATSL(1);pL=STATS(3);mL=B(1);bL=B(2);
            [B,BINT,R,RINT,STATS]=regress([G1_Obs;G2_Obs],[[G1_Est;G2_Est] ones(size([GDot;GDot]))]);
            R2=STATS(1);p=STATS(3);m=B(1);b=B(2);
            
            % Plot viscosity data for each temperature
            figure(1);clf;hold on;
            plot(GDot,G1_Obs,'ro');
            plot(GDot,G1_Est,'rx--');
            plot(GDot,G2_Obs,'bo');
            plot(GDot,G2_Est,'bx--');
            %Shear_Rates=10.^(-1:0.25:4);
            %Eta=Williamson(Shear_Rates,Coeffs);
            %plot(Shear_Rates,Eta,'k--');
            eta_char = char(951);
            %sWilliamson=sprintf('%s_0=%.0f, n=%.3f, k=%.4f; R^2=%.4f', ...
            %    eta_char,Coeffs(1),Coeffs(2),Coeffs(3),R2);
            legend('G1 Obs','G1 Est','G2 Obs','G2 Est','location','SouthWest');
            set(gca,'XScale','log','YScale','log');grid on;
            xlabel(sXlabel);ylabel('Storage & Loss Modulus, Pa');
            %ylim([10 3000]);

            % Fit the model and store statistics
            %EtaS=150;
            %EtaG=150;
            %Lambda=0.02;
            %L=2;
            %x0=[EtaS EtaG Lambda L Eta];
            %lb=[  10 10 0  1.74];
            %ub=[ 1e6 1e6 10 2.50];

            Coeffs_Eta=fmincon(@sumErrorEta_FENEP,x0,[],[],[],[],lb,ub,[],options);
            Eta_Est=CalcEta_FENEP(GDot,Coeffs_Eta);
            [B,BINT,R,RINT,STATS]=regress(log(Eta_Obs),[log(Eta_Est) ones(size(Eta_Est))]);
            R2=STATS(1);p=STATS(3);m=B(1);b=B(2);
            
            % Plot viscosity data for each temperature
            figure(2);clf;hold on;
            plot(GDot,Eta_Obs,'ko');
            Shear_Rates=10.^(-1:0.25:4);
            Eta=CalcEta_FENEP(Shear_Rates,Coeffs_Eta);
            plot(Shear_Rates,Eta,'k--');
            eta_char = char(951);
            sWilliamson=sprintf('%s_0=%.0f, n=%.3f, k=%.4f; R^2=%.4f', ...
                eta_char,Coeffs(1),Coeffs(2),Coeffs(3),R2);
            legend(sFileName(1:end-4),sWilliamson,'location','SouthWest');
            set(gca,'XScale','log','YScale','log');grid on;
            xlabel(sXlabel);ylabel('Dynamic viscosity, Pa.s');
            ylim([10 3000]);

        end
        Coeffs
        %pause
    end
%    Coeffs_G1G2
%    Coeffs_Eta
end    

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