% Copyright (c) 2012, Jan Philipp Kruse and Harro Schmeling, Universität Frankfurt
% All rights reserved.
% 
% Redistribution and use in source and binary forms, with or without modification, 
% are permitted provided that the following conditions are met:
% 
% Redistributions of source code must retain the above copyright notice, 
% this list of conditions and the following disclaimer.
% Redistributions in binary form must reproduce the above copyright notice, 
% this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution.
% 
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, 
% INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. 
% IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY,
% OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
% LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY,
% WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE,
% EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

clear all;
clc;

path([pwd,'/general'],path)
path([pwd,'/visc'],path)

eta_b0          = Inf;
eta_f           = 0;

eta_s0          = input('eta_s0 - shear viscosity (default is 10^18 [Pa s]) = ');   

if isempty(eta_s0)
    eta_s0      = 10^18;
end

missmax         = 10^-10*eta_s0;

alpha1          = 0;
while ~isempty(find(alpha1 == 0))
    alpha1      = str2num(input('alpha1 - aspect ratios of films (default is 0.01) = ','s'));
    if isempty(alpha1)
        alpha1  = .01;
    end
end
kappa           = str2num(input('kappa - tube shape, 0 for tapered, large for cylinders (default is 0) = ','s'));
if isempty(kappa)
    kappa       = 0;
end
alpha3          = str2num(input('alpha3 - aspect ratios of spheroids - (default is 1) = ','s'));
if isempty(alpha3)
    alpha3      = 1;
end
beta            = str2num(input('beta total [min max] - range of melt fractions (default is [10^-6 .16]) = ','s'));
if length(beta)<2 || length(beta)>2 || isempty(beta)
    beta        = [10^-6 .2];
end

betanum = input('number of values for beta (default is 10^3) = ');
if isempty(betanum)
    betanum  = 10^3;
end

betamin         = log10(beta(1));
betamax         = log10(beta(2));
beta            = logspace(betamin,betamax,betanum);

c_all           = cell(3,1);

c_all{1}        = 0;
c_all{2}        = 0;
c_all{3}        = 0;

while ~isempty(find(c_all{1}+c_all{2}+c_all{3} < 1)) || ~isempty(find(c_all{1}+c_all{2}+c_all{3} > 1))
    c_all{1}    = str2num(input('c1c - fraction of films of all melt (default is 0) = ','s'));
    c_all{2}    = str2num(input('c2c - fraction of tubes of all melt (default is 0) = ','s'));
    c_all{3}    = str2num(input('c3c - fraction of spheroids of all melt (default is 1) = ','s'));
    
    if isempty(c_all{1})
        c_all{1} = 0;
    end
    if isempty(c_all{2})
        c_all{2} = 0;
    end
    if isempty(c_all{3})
        c_all{3} = 1;
    end
end

nmax            = input('nmax - iterations, try 10 - 100 (default is 1000) = ');

if isempty(nmax)
    nmax        = 10^3;
end

nmax            = 2+nmax;

disp(['number of models = ',num2str(length(kappa)*length(alpha1)*length(alpha3)*length(c_all{1}))]);

models          = cell(1,length(kappa)*length(alpha1)*length(alpha3)*length(c_all{1}));
q               = 1;

for t=1:length(c_all{1})
    
    c(1) = c_all{1}(t);
    c(2) = c_all{2}(t);
    c(3) = c_all{3}(t);
    
    for z=1:length(kappa)
        for u=1:length(alpha1)
            for a=1:length(alpha3)
                
                disp(['Model No.',num2str(q)])
                
                Ar          = zeros(1,length(beta));
                eta_b       = zeros(1,length(beta));
                eta_s       = zeros(1,length(beta));
                nur         = zeros(1,length(beta));

                phi_tmp     = phi(alpha3(a));
                g_tmp       = g(alpha3(a),phi_tmp);

                for i=1:length(beta)
                    
                    
                    %% ETA' and RELAXED SHEARVISCOSITY %%
                    %% -------------------------------- %%
                    
                    %% START INDEX %%
                    n = 2;
                    
                    %% SET INITIAL CONDITIONS %%
                    %% ---------------------- %%
                    %% ETA_BR %%
                    eta_btmp1(1) = eta_s0; % c1*eta_s0*(c2(u)-beta(1))/beta(1);
                    eta_btmp1(2) = 0;
                    %% ETA_SR %%
                    eta_stmp1(1) = eta_s0;
                    %% INITIAL ERROR %%
                    miss1 = eta_s0;
                    miss2 = eta_s0;
                    %% ---------------------- %%
                    
                    %% RECURSION FOR ETA_BR & ETA_SR %%
                    %% ----------------------- %%
                    while (((miss1>=missmax) || (miss2>=missmax)) && (n<=nmax))
                        eta_btmp2   = 0;
                        eta_stmp2   = 0;
                        
                        for j=1:3
                            nu_tmp      = poisson(eta_btmp1(n-1),eta_stmp1(n-1));
                            R_tmp       = R(eta_btmp1(n-1),eta_stmp1(n-1));
                            B_tmp       = B(j,eta_f,eta_btmp1(n-1),nu_tmp,kappa(z),eta_stmp1(n-1));
                            theta_tmp   = theta(j,eta_btmp1(n-1),eta_stmp1(n-1),nu_tmp,alpha1(u),kappa(z),g_tmp,phi_tmp,R_tmp);
                            A_tmp       = A(j,alpha1(u),alpha3(a),eta_btmp1(n-1),g_tmp,R_tmp,1,phi_tmp,B_tmp,eta_stmp1(n-1),nu_tmp,kappa(z),eta_f);
                            eta_stmp2   = eta_stmp2+eta_srelaxed(A_tmp,c(j)*beta(i));
                            eta_btmp2   = eta_btmp2+eta_brelaxed(theta_tmp,c(j)*beta(i));
                        end

                        eta_btmp1(n)    = 1/eta_btmp2;
                        eta_stmp1(n)    = 1/(1/eta_stmp1(1)+eta_stmp2);

                        %% ERROR FOR ETA_SR AND ETA_BR %%
                        miss1           = abs(eta_btmp1(n)-eta_btmp1(n-1));
                        miss2           = abs(eta_stmp1(n)-eta_stmp1(n-1));
                        
                        n = n + 1;
                    end
                    %% ----------------------- %%
                    % disp(['ETA_BRE calculated for beta = ',num2str(beta(i)),' ! -> error @ iteration ',num2str(n-1),' = ',num2str(miss1)]);

                    %% SET VALUES FOR BAR AND MY_RELAXED TERMS %%
                    %% --------------------------------------- %%
                    
                    eta_s(i)      = eta_stmp1(n-1);
                    eta_b(i)      = eta_btmp1(n-1);
                    Ar(i)         = A_tmp;
                    nur(i)        = poisson(eta_b(i),eta_s(i));
                    %% ---------------------------- %%
                end
                models{q} = struct('No',q,'kappa',kappa(z),'alpha1',alpha1(u),'alpha3',alpha3(a), ...
                                   'c1',c(1),'c2',c(2),'c3',c(3), ...
                                   'eta_b',eta_b,'eta_s',eta_s,'nu_r',nur,'A_r',Ar);
                q = q + 1;
            end
        end
    end
end

eta_b       = zeros(length(models),length(beta));
eta_s       = zeros(length(models),length(beta));
Ar          = zeros(length(models),length(beta));
nur         = zeros(length(models),length(beta));

for i=1:length(models)
    eta_b(i,:)      = models{i}.eta_b;
    eta_s(i,:)      = models{i}.eta_s;
    nur(i,:)        = models{i}.nu_r;
    Ar(i,:)         = models{i}.A_r;
end

logplot     = input('0 for normal plot, 1 for semilogy plot, 2 for loglog plot (default is 2) = ');
if isempty(logplot)
    logplot = 2;
end

etaplot     = input('0 for eta_b and eta_s in different plots, 1 for both in one plot (default is 0) = ');
if isempty(etaplot)
    etaplot = 0;
end

plotuntil   = input('plot values until eta_b resp. eta_s > eta_s0 * (default is 10^-3) ');
if isempty(plotuntil)
    plotuntil   = 10^-3;
end

for i=1:length(models)
    col(i,:)= uisetcolor(['model',num2str(i),' color']);
end

if etaplot==1
    leg     = cell(1,length(models)*2);
    pind    = 1;
    for i=1:length(models)
        leg{pind} = ['model',num2str(i),', \eta_b'];
        if models{i}.c1 ~= 0
            leg{pind} = [leg{pind},', \alpha_1 = ',num2str(models{i}.alpha1),', c_1 = ',num2str(models{i}.c1*100),'%'];
        end
        if models{i}.c2 ~= 0
            leg{pind} = [leg{pind},', \kappa = ',num2str(models{i}.kappa),', c_2 = ',num2str(models{i}.c2*100),'%'];
        end
        if models{i}.c3 ~= 0
            leg{pind} = [leg{pind},', \alpha_3 = ',num2str(models{i}.alpha3),', c_3 = ',num2str(models{i}.c3*100),'%'];
        end
        pind    = pind+1;
        leg{pind} = ['model',num2str(i),', \eta_s'];
        if models{i}.c1 ~= 0
            leg{pind} = [leg{pind},', \alpha_1 = ',num2str(models{i}.alpha1),', c_1 = ',num2str(models{i}.c1*100),'%'];
        end
        if models{i}.c2 ~= 0
            leg{pind} = [leg{pind},', \kappa = ',num2str(models{i}.kappa),', c_2 = ',num2str(models{i}.c2*100),'%'];
        end
        if models{i}.c3 ~= 0
            leg{pind} = [leg{pind},', \alpha_3 = ',num2str(models{i}.alpha3),', c_3 = ',num2str(models{i}.c3*100),'%'];
        end
        pind    = pind+1;
    end
else
    leg = cell(1,length(models));
    for i=1:length(models)
        leg{i} = ['model',num2str(i)];
        if models{i}.c1 ~= 0
            leg{i} = [leg{i},', \alpha_1 = ',num2str(models{i}.alpha1),', c_1 = ',num2str(models{i}.c1*100),'%'];
        end
        if models{i}.c2 ~= 0
            leg{i} = [leg{i},', \kappa = ',num2str(models{i}.kappa),', c_2 = ',num2str(models{i}.c2*100),'%'];
        end
        if models{i}.c3 ~= 0
            leg{i} = [leg{i},', \alpha_3 = ',num2str(models{i}.alpha3),', c_3 = ',num2str(models{i}.c3*100),'%'];
        end
    end
end
if logplot==0
    if etaplot==0
        figbulk         = figure('Name','bulk-viscosity','NumberTitle','Off');
        for i=1:length(models)
            plot(beta(find(eta_b(i,:) >= eta_s0*plotuntil)),eta_b(i,find(eta_b(i,:) >= eta_s0*plotuntil)),'.','Color',col(i,:))
        end
        xlabel('total melt fraction \beta');
        ylabel('\eta_b [Pa s]');
        legend(leg);
        
        figshear        = figure('Name','shear-viscosity','NumberTitle','Off');
        for i=1:length(models)
            plot(beta(find(eta_b(i,:) >= eta_s0*plotuntil)),eta_s(i,find(eta_s(i,:) >= eta_s0*plotuntil)),'.','Color',col(i,:))
        end
        xlabel('total melt fraction \beta');
        ylabel('\eta_s [Pa s]');
        legend(leg);
    elseif etaplot==1
        figbs           = figure('Name','bulk- and shear-viscosity','NumberTitle','Off');
        for i=1:length(models)
            plot(beta(find(eta_b(i,:) >= eta_s0*plotuntil)),eta_b(i,find(eta_b(i,:) >= eta_s0*plotuntil)),'-','Color',col(i,:))
            hold on;
            plot(beta(find(eta_b(i,:) >= eta_s0*plotuntil)),eta_s(i,find(eta_s(i,:) >= eta_s0*plotuntil)),'--','Color',col(i,:))
        end
        xlabel('total melt fraction \beta');
        ylabel('\eta [Pa s]');
        title('bulk (\eta_b)- and shear (\eta_s) viscosities')
        legend(leg);
    end
elseif logplot==1
    if etaplot==0
        figbulk         = figure('Name','bulk-viscosity','NumberTitle','Off');
        for i=1:length(models)
            semilogy(beta(find(eta_b(i,:) >= eta_s0*plotuntil)),eta_b(i,find(eta_b(i,:) >= eta_s0*plotuntil)),'.','Color',col(i,:))
        end
        xlabel('total melt fraction \beta');
        ylabel('\eta_b [Pa s]');
        legend(leg);
        
        figshear        = figure('Name','shear-viscosity','NumberTitle','Off');
        for i=1:length(models)
            semilogy(beta(find(eta_b(i,:) >= eta_s0*plotuntil)),eta_s(i,find(eta_s(i,:) >= eta_s0*plotuntil)),'.','Color',col(i,:))
        end
        xlabel('total melt fraction \beta');
        ylabel('\eta_s [Pa s]');
        legend(leg);
    elseif etaplot==1
        figbs           = figure('Name','bulk- and shear-viscosity','NumberTitle','Off');
        for i=1:length(models)
            semilogy(beta(find(eta_b(i,:) >= eta_s0*plotuntil)),eta_b(i,find(eta_b(i,:) >= eta_s0*plotuntil)),'-','Color',col(i,:))
            hold on;
            semilogy(beta(find(eta_b(i,:) >= eta_s0*plotuntil)),eta_s(i,find(eta_s(i,:) >= eta_s0*plotuntil)),'--','Color',col(i,:))
        end
        xlabel('total melt fraction \beta');
        ylabel('\eta [Pa s]');
        title('bulk (\eta_b)- and shear (\eta_s) viscosities')
        legend(leg);
    end
elseif logplot==2
    if etaplot==0
        figbulk         = figure('Name','bulk-viscosity','NumberTitle','Off');
        for i=1:length(models)
            loglog(beta(find(eta_b(i,:) >= eta_s0*plotuntil)),eta_b(i,find(eta_b(i,:) >= eta_s0*plotuntil)),'-','Color',col(i,:))
        end
        xlabel('total melt fraction \beta');
        ylabel('\eta_b [Pa s]');
        legend(leg);
        
        figshear        = figure('Name','shear-viscosity','NumberTitle','Off');
        for i=1:length(models)
            loglog(beta(find(eta_s(i,:) >= eta_s0*plotuntil)),eta_s(i,find(eta_s(i,:) >= eta_s0*plotuntil)),'-','Color',col(i,:))
        end
        xlabel('total melt fraction \beta');
        ylabel('\eta_s [Pa s]');
        legend(leg);
    elseif etaplot==1
        figbs           = figure('Name','bulk- and shear-viscosity','NumberTitle','Off');
        for i=1:length(models)
            loglog(beta(find(eta_b(i,:) >= eta_s0*plotuntil)),eta_b(i,find(eta_b(i,:) >= eta_s0*plotuntil)),'-','Color',col(i,:))
            hold on;
            loglog(beta(find(eta_s(i,:) >= eta_s0*plotuntil)),eta_s(i,find(eta_s(i,:) >= eta_s0*plotuntil)),'--','Color',col(i,:))
        end
        xlabel('total melt fraction \beta');
        ylabel('\eta [Pa s]');
        title('bulk (\eta_b)- and shear (\eta_s) viscosities')
        legend(leg);
    end
end