% 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,'/moduli'],path)

K0      = input('K0 - bulk modulus matrix (default is 0.66 [M Pa]) = ');
if isempty(K0)
    K0  = .66;
end

missmax1    = 10^-11*K0;

mu0     = input('mu0 - shear modulus matrix (default is 0.4 [MPa]) = ');
if isempty(mu0)
    mu0  = .4;
end

missmax2    = 10^-11*mu0;

Kf      = input('Kf - bulk modulus melt (default is 0.2 [MPa]) = ');            
if isempty(Kf)
    Kf  = .2;
end

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 .16];
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,10^3);


c_all   = cell(3,2);

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

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

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

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

nmax        = 2+nmax;

unit        = input('unit for plotting (default is MPa): ');
if isempty(unit)
    unit    = 'MPa';
end

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

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

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

                Kbar   = zeros(1,length(beta));
                Abar   = zeros(1,length(beta));
                nubar  = zeros(1,length(beta));

                Kiso   = zeros(1,length(beta));
                muisoo = zeros(1,length(beta));
                Aiso   = zeros(1,length(beta));
                nuiso  = zeros(1,length(beta));

                Kr     = zeros(1,length(beta));
                mur    = 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)
                    
                    %% UNRELAXED MODULI %%
                    %% ---------------- %%
                    
                    %% START INDEX %%
                    n                   = 2;
                    
                    %% SET INITIAL CONDITIONS %%
                    %% ---------------------- %%
                    
                    %% K %%
                    Ku_tmp(1)           = K0;
                    Ku_tmp(2)           = 0;
                    
                    %% MU %%
                    muu_tmp(1)          = mu0;
                    
                    %% INITIAL ERROR %%
                    miss1               = K0;
                    miss2               = mu0;
                    %% ---------------------- %%
                    
                    %% RECURSION FOR K_U & MU_U %%
                    %% ------------------------ %%
                    
                    while (((miss1>=missmax1) || (miss2>=missmax2)) && (n<=nmax))
                        
                        K_tmp           = 0;
                        mu_tmp          = 0;
                        
                        for j=1:3
                            nu_tmp      = poisson(Ku_tmp(n-1),muu_tmp(n-1));
                            R_tmp       = R(Ku_tmp(n-1),muu_tmp(n-1));
                            B_tmp       = B(j,Kf,Ku_tmp(n-1),nu_tmp,kappa(z),muu_tmp(n-1));
                            theta_tmp   = theta(j,Ku_tmp(n-1),muu_tmp(n-1),nu_tmp,alpha1(u),kappa(z),g_tmp,phi_tmp,R_tmp);
                            D_tmp       = D(j,Kf,Ku_tmp(1),theta_tmp);
                            A_tmp       = A(j,alpha1(u),alpha3(a),Ku_tmp(n-1),g_tmp,R_tmp,D_tmp,phi_tmp,B_tmp,muu_tmp(n-1),nu_tmp,kappa(z),Kf);
                            mu_tmp      = mu_tmp+muunrelaxed(A_tmp,c(j,1)*beta(i),c(j,2)*beta(i));
                            K_tmp       = K_tmp+KModulunrelaxed(Kf,K0,c(j,1)*beta(i),c(j,2)*beta(i),theta_tmp,Ku_tmp(n-1));
                        end
                        
                        Ku_tmp(n)       = 1/(1/K0+K_tmp);
                        muu_tmp(n)      = 1/(1/mu0+mu_tmp);
                        
                        %% ERROR FOR MY AND K %%
                        miss1           = abs(Ku_tmp(n)-Ku_tmp(n-1));
                        miss2           = abs(muu_tmp(n)-muu_tmp(n-1));
                        
                        n               = n+1;
                    end
                    
%                     disp(['K_u calculated for beta = ',num2str(beta(i)),' ! -> error @ iteration ',num2str(n-1),' = ',num2str(miss1)]);
                    
                    %% SET VALUES FOR UNRELAXED TERMS %%
                    %% ------------------------------ %%
                    Ku(i)               = Ku_tmp(n-1);
                    muu(i)              = muu_tmp(n-1);
                    Au(i)               = A_tmp;
                    nuu(i)              = nu_tmp;
                    %% ------------------------------ %%
                    
         
                    %% KMODULI' and RELAXED SHEARMODULI %%
                    %% -------------------------------- %%
                    
                    %% START INDEX %%
                    n                   = 2;
                    
                    %% SET INITIAL CONDITIONS %%
                    %% ---------------------- %%
                    
                    %% K'                               %%
                    Kbar_tmp(1)         = K0;
                    Kbar_tmp(2)         = 0;
                    
                    %% MY_R                             %%
                    mubar_tmp(1)        = mu0;
                    
                    %% INITIAL ERROR                    %%
                    miss1               = K0;
                    miss2               = mu0;
                    %% -------------------------------- %%
                    
                    %% RECURSION FOR K' & MY_R          %%
                    %% -------------------------------- %%
                    
                    while (((miss1>=missmax1) || (miss2>=missmax2)) && (n<=nmax))
                        K_tmp           = 0;
                        mu_tmp          = 0;
                        
                        for j=1:3
                            nu_tmp      = poisson(Kbar_tmp(n-1),mubar_tmp(n-1));
                            R_tmp       = R(Kbar_tmp(n-1),mubar_tmp(n-1));
                            B1_tmp      = B(j,Kf,Kbar_tmp(n-1),nu_tmp,kappa(z),mubar_tmp(n-1));
                            B2_tmp      = B(j,0,Kbar_tmp(n-1),nu_tmp,kappa(z),mubar_tmp(n-1));
                            theta_tmp   = theta(j,Kbar_tmp(n-1),mubar_tmp(n-1),nu_tmp,alpha1(u),kappa(z),g_tmp,phi_tmp,R_tmp);
                            D1_tmp      = D(j,Kf,Kbar_tmp(1),theta_tmp);
                            D2_tmp      = D(j,0,Kbar_tmp(1),theta_tmp);
                            A1_tmp      = A(j,alpha1(u),alpha3(a),Kbar_tmp(n-1),g_tmp,R_tmp,D1_tmp,phi_tmp,B1_tmp,mubar_tmp(n-1),nu_tmp,kappa(z),Kf);
                            A2_tmp      = A(j,alpha1(u),alpha3(a),Kbar_tmp(n-1),g_tmp,R_tmp,D2_tmp,phi_tmp,B2_tmp,mubar_tmp(n-1),nu_tmp,kappa(z),0);
                            mu_tmp      = mu_tmp+murelaxed(A1_tmp,A2_tmp,c(j,1)*beta(i),c(j,2)*beta(i));
                            K_tmp       = K_tmp+KModulbar(Kf,K0,theta_tmp,Kbar_tmp(n-1),c(j,1)*beta(i),c(j,2)*beta(i));
                        end
                        
                        Kbar_tmp(n)     = 1/(1/K0+K_tmp);
                        mubar_tmp(n)    = 1/(1/mu0+mu_tmp);
                        
                        %% ERROR FOR MY_R AND K'            %%
                        miss1           = abs(Kbar_tmp(n)-Kbar_tmp(n-1));
                        miss2           = abs(mubar_tmp(n)-mubar_tmp(n-1));
                        %% -------------------------------- %%
                        
                        n               = n+1;
                    end
                    
%                     disp(['K'' calculated for beta = ',num2str(beta(i)),' ! -> error @ iteration ',num2str(n-1),' = ',num2str(miss1)]);
                    
                    %% SET VALUES FOR BAR AND MY_RELAXED TERMS %%
                    %% --------------------------------------- %%
                    mur(i)              = mubar_tmp(n-1);
                    Kbar(i)             = Kbar_tmp(n-1);
                    Abar(i)             = A_tmp;
                    nubar(i)            = nu_tmp;
                    %% --------------------------------------- %%
                    
                    
                    %% --------------------------------------- %%
                    %% KMODULI ISOLATED                        %%
                    %% --------------------------------------- %%
                    
                    %% START INDEX                             %%
                    n                   = 2;
                    
                    %% SET INITIAL CONDITIONS                  %%
                    %% --------------------------------------- %%
                    
                    %% KISO                                    %%
                    Kiso_tmp(1)         = K0;
                    Kiso_tmp(2)         = 0;
                    
                    %% MYISO                                   %%
                    muiso_tmp(1)        = mu0;
                    
                    %% INITIAL ERROR                           %%
                    miss1               = K0;
                    miss2               = mu0;
                    
                    %% --------------------------------------- %%
                    %% RECURSION FOR KISO & MUISO              %%
                    %% --------------------------------------- %%
                    
                    while (((miss1>=missmax1) || (miss2>=missmax2)) && (n<=nmax))
                        K_tmp           = 0;
                        mu_tmp          = 0;
                        
                        for j=1:3
                            betamod     = (c(j,1)*beta(i))/(1-c(j,2)*beta(i));
                            nu_tmp      = poisson(Kiso_tmp(n-1),muiso_tmp(n-1));
                            R_tmp       = R(Kiso_tmp(n-1),muiso_tmp(n-1));
                            B_tmp       = B(j,Kf,Kiso_tmp(n-1),nu_tmp,kappa(z),muiso_tmp(n-1));
                            theta_tmp   = theta(j,Kiso_tmp(n-1),muiso_tmp(n-1),nu_tmp,alpha1(u),kappa(z),g_tmp,phi_tmp,R_tmp);
                            D_tmp       = D(j,Kf,Kiso_tmp(1),theta_tmp);
                            A_tmp       = A(j,alpha1(u),alpha3(a),Kiso_tmp(n-1),g_tmp,R_tmp,D_tmp,phi_tmp,B_tmp,muiso_tmp(n-1),nu_tmp,kappa(z),Kf);
                            mu_tmp      = mu_tmp+muiso(A_tmp,betamod);
                            K_tmp       = K_tmp+KModuliso(Kf,K0,betamod,theta_tmp,Kiso_tmp(n-1));
                        end
                        
                        Kiso_tmp(n)     = 1/(1/K0+K_tmp);
                        muiso_tmp(n)    = 1/(1/mu0+mu_tmp);
                        
                        %% ERROR FOR MY_ISO AND K_ISO              %%
                        miss1           = abs(Kiso_tmp(n)-Kiso_tmp(n-1));
                        miss2           = abs(muiso_tmp(n)-muiso_tmp(n-1));
                        %% --------------------------------------- %%
                        
                        n               = n+1;
                    end
                    
%                     disp(['K_iso calculated for betamod = ',num2str(beta(i)),' ! -> error @ iteration ',num2str(n-1),' = ',num2str(miss1)]);
                    
                    %% SET VALUES FOR ISOLATED TERMS           %%
                    %% --------------------------------------- %%
                    Kiso(i)             = Kiso_tmp(n-1);
                    muisoo(i)           = muiso_tmp(n-1);
                    Aiso(i)             = A_tmp;
                    nuiso(i)            = nu_tmp;
                    %% --------------------------------------- %%
                    
                    %% SET VALUES FOR RELAXED TERMS            %%
                    %% --------------------------------------- %%
                    Kr(i)               = KModulrelaxed(Kf,Kiso(i),Kbar(i),(c(1,2)+c(2,2)+c(3,2))*beta(i));
                    nur(i)              = poisson(Kr(i),mur(i));
                    %% --------------------------------------- %%
                    
                end
                
                models{q}               = struct('No',q,'betamin',10^betamin,'betamax',10^betamax,'beta',beta,'K_0',K0,'mu_0',mu0,'K_f',Kf,...
                                                 'alpha1',alpha1(u),'alpha3',alpha3(a),'kappa',kappa(z),...
                                                 'c1i',c(1,1),'c1c',c(1,2),'c2i',c(2,1),'c2c',c(2,2),'c3i',c(3,1),'c3c',c(3,2),'nmax',nmax,...
                                                 'K_u',Ku,'K_bar',Kbar,'K_iso',Kiso,'K_r',Kr,...
                                                 'mu_u',muu,'mu_iso',muisoo,'mu_r',mur,...
                                                 'nu_u',nuu,'nu_iso',nuiso,'nu_bar',nubar,'nu_r',nur,...
                                                 'A_u',Au,'A_bar',Abar,'A_iso',Aiso,...
                                                 'unit',unit,'betanum',betanum);
                q                       = q+1;
            end
        end
    end
end

Ku     = zeros(length(models),length(beta));
muu    = zeros(length(models),length(beta));
nuu    = zeros(length(models),length(beta));

Kbar   = zeros(length(models),length(beta));
nubar  = zeros(length(models),length(beta));

Kiso   = zeros(length(models),length(beta));
muisoo = zeros(length(models),length(beta));
nuiso  = zeros(length(models),length(beta));

Kr     = zeros(length(models),length(beta));
mur    = zeros(length(models),length(beta));
nur    = zeros(length(models),length(beta));

for i=1:length(models)
    Ku(i,:)     = models{i}.K_u;
    muu(i,:)    = models{i}.mu_u;
    nuu(i,:)    = models{i}.nu_u;
    Kbar(i,:)   = models{i}.K_bar;
    nubar(i,:)  = models{i}.nu_bar;
    Kiso(i,:)   = models{i}.K_iso;
    muisoo(i,:) = models{i}.mu_iso;
    nuiso(i,:)  = models{i}.nu_iso;
    Kr(i,:)     = models{i}.K_r;
    mur(i,:)    = models{i}.mu_r;
    nur(i,:)    = models{i}.nu_r;
end

leg     = cell(1,length(models));

for i=1:length(models)
    leg{i}  = ['model',num2str(i)];
    if models{i}.c1i~=0 || models{i}.c1c~=0
        leg{i} = [leg{i},', \alpha_1=',num2str(models{i}.alpha1)];
        if models{i}.c1i~=0
            leg{i} = [leg{i},', c_{1i}=',num2str(models{i}.c1i*100),'%'];
        end
        if models{i}.c1c~=0
            leg{i} = [leg{i},', c_{1c}=',num2str(models{i}.c1c*100),'%'];
        end
    end
    if models{i}.c2i~=0 || models{i}.c2c~=0
        leg{i} = [leg{i},', \kappa=',num2str(models{i}.kappa)];
        if models{i}.c2i~=0
            leg{i} = [leg{i},', c_{2i}=',num2str(models{i}.c2i*100),'%'];
        end
        if models{i}.c2c~=0
            leg{i} = [leg{i},', c_{2c}=',num2str(models{i}.c2c*100),'%'];
        end
    end
    if models{i}.c3i~=0 || models{i}.c3c~=0
        leg{i} = [leg{i},', \alpha_3=',num2str(models{i}.alpha3)];
        if models{i}.c3i~=0
            leg{i} = [leg{i},', c_{3i}=',num2str(models{i}.c3i*100),'%'];
        end
        if models{i}.c3c~=0
            leg{i} = [leg{i},', c_{3c}=',num2str(models{i}.c3c*100),'%'];
        end
    end
end

plunbulk        = input('plot unrelaxed bulk-modulus - y/n (default is y): ','s');
plrbulk         = input('plot relaxed bulk-modulus - y/n (default is y): ','s');
plunshear       = input('plot unrelaxed shear-modulus - y/n (default is y): ','s');
plrshear        = input('plot relaxed shear-modulus - y/n (default is y): ','s');

if isempty(plunbulk)
    plunbulk    = 'y';
end

if strcmp(plunbulk,'y')==1
    figKu       = figure('Name','unrelaxed bulk-modulus','NumberTitle','Off');
    
    plot(beta,Ku,'-')
    xlabel('total melt fraction \beta');
    ylabel(['unrelaxed bulk-modulus K_u [',unit,']'])
    legend(leg)
end

if isempty(plrbulk)
    plrbulk    = 'y';
end

if strcmp(plrbulk,'y')==1
    figKr       = figure('Name','relaxed bulk-modulus','NumberTitle','Off');
    
    plot(beta,Kr,'-')
    xlabel('total melt fraction \beta');
    ylabel(['relaxed bulk-modulus K_r [',unit,']'])
    legend(leg)
end

if isempty(plunshear)
    plunshear    = 'y';
end

if strcmp(plunshear,'y')==1
    figmuu       = figure('Name','unrelaxed shear-modulus','NumberTitle','Off');
    
    plot(beta,muu,'-')
    xlabel('total melt fraction \beta');
    ylabel(['unrelaxed shear-modulus \mu_u [',unit,']'])
    legend(leg)
end

if isempty(plrshear)
    plrshear    = 'y';
end

if strcmp(plrshear,'y')==1
    figmur       = figure('Name','relaxed shear-modulus','NumberTitle','Off');
    
    plot(beta,mur,'-')
    xlabel('total melt fraction \beta');
    ylabel(['relaxed shear-modulus \mu_r [',unit,']'])
    legend(leg)
end