clear all;
clc;

KMod_0  = .66; % input('K_o (Bulk modulus matrix)= ');      
SMod_0  = .4; % input('My_o (Shear modulus matrix)= ');   
Kf      =  .2; % input('Kf (Bulk modulus melt) = ');            

missmax = 10^-10;

alpha1  = 0; % str2num(input('alpha1(aspect ratios of films) = ','s'));
alpha3  = [.01 .0316 .05 .1 .2 .3 1]; % str2num(input('alpha3(aspect ratios of spheroids) = ','s'));

beta = [0 .16]; % str2num(input('beta total [min max] (range of melt fractions)= ','s'));

betamin = beta(1);
betamax = beta(2);
beta = betamin:.0001:betamax;

c_all = cell(3,2);

c_all{1,1} = 0; % [0 0 0 0 0 0 0]; % str2num(input('c1i (fraction of isolated films of all melt)= ','s'));
c_all{1,2} = 0; % [1 .75 .5 .25 .1 .05 0]; % str2num(input('c1c (fraction of connected films of all melt)= ','s'));

c_all{2,1} = 0; % [0 0 0 0 0 0 0]; % str2num(input('c2i (fraction of isolated tubes of all melt)= ','s'));
c_all{2,2} = 0; % [0 0 0 0 0 0 0]; % str2num(input('c2c (fraction of connected tubes of all melt)= ','s'));

c_all{3,1} = 1; % [0 0 0 0 0 0 0]; % str2num(input('c3i (fraction of isolated spheroids of all melt)= ','s'));
c_all{3,2} = 0; % [0 .25 .5 .75 .9 .95 1]; % str2num(input('c3c (fraction of connected spheroids of all melt)= ','s'));

kappa = 0; % str2num(input('kappa (tube shape, 0 for tapered, large for cylinders)= ','s'));

nmax = 10^2; % input('nmax (iterations, try 10 - 100)= ');

nmax = 2 + nmax;

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)])
                
                Kun = zeros(1,length(beta));
                MYun = zeros(1,length(beta));
                Aun = zeros(1,length(beta));
                NUun = zeros(1,length(beta));

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

                Kiso = zeros(1,length(beta));
                MYiso = zeros(1,length(beta));
                Aiso = zeros(1,length(beta));
                NUiso = zeros(1,length(beta));

                Kre = zeros(1,length(beta));
                MYre = zeros(1,length(beta));
                NUre = zeros(1,length(beta));

                PHI_temp = phi(alpha3(a));
                G_temp = g(alpha3(a),PHI_temp);

                for i=1:length(beta)
                    
                    %% UNRELAXED MODULI %%
                    %% ---------------- %%
                    
                    %% START INDEX %%
                    n = 2;
                    
                    %% SET INITIAL CONDITIONS %%
                    %% ---------------------- %%
                    %% K %%
                    K1(1) = KMod_0;
                    K1(2) = 0;
                    %% MY %%
                    MY1(1) = SMod_0;
                    %% INITIAL ERROR %%
                    miss1 = 1;
                    miss2 = 1;
                    %% ---------------------- %%
                    
                    %% RECURSION FOR K_U & MY_U %%
                    %% ------------------------ %%
                    
                    while ((miss1>=missmax) && (miss2>=missmax) && (n<=nmax))
                        
                        K_temp = 0;
                        MY_temp = 0;
                        
                        for j=1:3
                            NU_temp     = poisson(K1(n-1),MY1(n-1));
                            R_temp      = R(MY1(n-1),K1(n-1));
                            B_temp      = B(j,Kf,K1(n-1),NU_temp,kappa(z),MY1(n-1));
                            theta_temp  = theta(j,K1(n-1),NU_temp,alpha1(u),kappa(z),G_temp,PHI_temp,R_temp);
                            D_temp      = D(j,Kf,K1(1),theta_temp);
                            A_temp      = A(j,alpha1(u),alpha3(a),K1(n-1),G_temp,R_temp,D_temp,PHI_temp,B_temp,MY1(n-1),NU_temp,kappa(z),Kf);
                            MY_temp     = MY_temp+muunrelaxed(A_temp,c(j,1)*beta(i),c(j,2)*beta(i));
                            K_temp      = K_temp+KModulunrelaxed(Kf,K1(1),c(j,1)*beta(i),c(j,2)*beta(i),theta_temp,K1(n-1));
                        end
                        
                        K1(n) = 1/(1/K1(1)+K_temp);
                        MY1(n) = 1/(1/MY1(1)+MY_temp);
                        
                        %% ERROR FOR MY AND K %%
                        miss1 = abs(K1(n)-K1(n-1));
                        miss2 = abs(MY1(n)-MY1(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 %%
                    %% ------------------------------ %%
                    Kun(i) = K1(n-1);
                    MYun(i) = MY1(n-1);
                    Aun(i) = A_temp;
                    NUun(i) = NU_temp;
                    %% ------------------------------ %%
                    
                    
                    
                    %% KMODULI' and RELAXED SHEARMODULI %%
                    %% -------------------------------- %%
                    
                    %% START INDEX %%
                    n = 2;
                    
                    %% SET INITIAL CONDITIONS %%
                    %% ---------------------- %%
                    %% K' %%
                    K2(1) = KMod_0;
                    K2(2) = 0;
                    %% MY_R %%
                    MY2(1) = SMod_0;
                    %% INITIAL ERROR %%
                    miss1 = 1;
                    miss2 = 1;
                    %% ---------------------- %%
                    
                    %% RECURSION FOR K' & MY_R %%
                    %% ----------------------- %%
                    while (miss1>=missmax && miss2>=missmax && n<=nmax)
                        K_temp  = 0;
                        MY_temp = 0;
                        for j=1:3
                            NU_temp     = poisson(K2(n-1),MY2(n-1));
                            R_temp      = R(MY2(n-1),K2(n-1));
                            B1_temp     = B(j,Kf,K2(n-1),NU_temp,kappa(z),MY2(n-1));
                            B2_temp     = B(j,0,K2(n-1),NU_temp,kappa(z),MY2(n-1));
                            theta_temp  = theta(j,K2(n-1),NU_temp,alpha1(u),kappa(z),G_temp,PHI_temp,R_temp);
                            D1_temp     = D(j,Kf,K2(1),theta_temp);
                            D2_temp     = D(j,0,K2(1),theta_temp);
                            A1_temp     = A(j,alpha1(u),alpha3(a),K2(n-1),G_temp,R_temp,D1_temp,PHI_temp,B1_temp,MY2(n-1),NU_temp,kappa(z),Kf);
                            A2_temp     = A(j,alpha1(u),alpha3(a),K2(n-1),G_temp,R_temp,D2_temp,PHI_temp,B2_temp,MY2(n-1),NU_temp,kappa(z),0);
                            MY_temp     = MY_temp+murelaxed(A1_temp,A2_temp,c(j,1)*beta(i),c(j,2)*beta(i));
                            K_temp      = K_temp+KModulbar(Kf,K2(1),theta_temp,K2(n-1),c(j,1)*beta(i),c(j,2)*beta(i));
                        end

                        K2(n)   = 1/(1/K2(1)+K_temp);
                        MY2(n)  = 1/(1/MY2(1)+MY_temp);

                        %% ERROR FOR MY_R AND K' %%
                        miss1 = abs(K2(n)-K2(n-1));
                        miss2 = abs(MY2(n)-MY2(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 %%
                    %% --------------------------------------- %%
                    MYre(i)     = MY2(n-1);
                    Kbar(i)     = K2(n-1);
                    Abar(i)     = A1_temp;
                    NUbar(i)    = NU_temp;
                    %% --------------------------------------- %%
                    
                    %% KMODULI ISOLATED %%
                    %% ---------------- %%
                    
                    %% START INDEX %%
                    n = 2;
                    
                    %% SET INITIAL CONDITIONS %%
                    %% ---------------------- %%
                    %% K' %%
                    K3(1)   = KMod_0;
                    K3(2)   = 0;
                    %% MY_R %%
                    MY3(1)  = SMod_0;
                    %% INITIAL ERROR %%
                    miss1   = 1;
                    miss2   = 1;
                    %% ---------------------- %%

                    %% RECURSION FOR K_ISO & MY_ISO %%
                    %% ---------------------------- %%
                    while (miss1 >= missmax && miss2 >= missmax && n <= nmax)

                        K_temp  = 0;
                        MY_temp = 0;

                        for j=1:3
                            betamod     = (c(j,1)*beta(i))/(1-c(j,2)*beta(i));
                            NU_temp     = poisson(K3(n-1),MY3(n-1));
                            R_temp      = R(MY3(n-1),K3(n-1));
                            B_temp      = B(j,Kf,K3(n-1),NU_temp,kappa(z),MY3(n-1));
                            theta_temp  = theta(j,K3(n-1),NU_temp,alpha1(u),kappa(z),G_temp,PHI_temp,R_temp);
                            D_temp      = D(j,Kf,K3(1),theta_temp);
                            A_temp      = A(j,alpha1(u),alpha3(a),K3(n-1),G_temp,R_temp,D_temp,PHI_temp,B_temp,MY3(n-1),NU_temp,kappa(z),Kf);
                            MY_temp     = MY_temp+muiso(A_temp,betamod);
                            K_temp      = K_temp+KModuliso(Kf,K3(1),betamod,theta_temp,K3(n-1));
                        end

                        K3(n)   = 1/(1/K3(1)+K_temp);
                        MY3(n)  = 1/(1/MY3(1)+MY_temp);

                        %% ERROR FOR MY_ISO AND K_ISO %%
                        miss1 = abs(K3(n)-K3(n-1));
                        miss2 = abs(MY3(n)-MY3(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)     = K3(n-1);           
                    MYiso(i)    = MY3(n-1);
                    Aiso(i)     = A_temp;
                    NUiso(i)    = NU_temp;
                    %% ----------------------------- %%
                    
                    %% SET VALUES FOR RELAXED TERMS %%
                    %% ---------------------------- %%
                    Kre(i)      = KModulrelaxed(Kf,Kiso(i),Kbar(i),(c(1,2)+c(2,2)+c(3,2))*beta(i));
                    NUre(i)     = poisson(Kre(i),MYre(i));
                    %% ---------------------------- %%
                end
%                 plot(beta,Kun)
%                 pause
                models{q} = struct('No',q,'kappa',kappa(z),'alpha1',alpha1(u),'alpha3',alpha3(a), ...
                                   'c1i',c(1,1),'c1c',c(1,2),'c2i',c(2,1),'c2c',c(2,2),'c3i',c(3,1),'c3c',c(3,2), ...
                                   'Kun',Kun,'Kbar',Kbar,'Kiso',Kiso,'Kre',Kre, ...
                                   'MYun',MYun,'MYiso',MYiso,'MYre',MYre, ...
                                   'NUun',NUun,'NUiso',NUiso,'NUbar',NUbar,'NUre',NUre, ...
                                   'Aun',NUun,'Aiso',Aiso,'Are',Abar);
                q = q + 1;
            end
        end
    end
end

Kun = zeros(length(models),length(beta));
MYun = zeros(length(models),length(beta));
Aun = zeros(length(models),length(beta));
NUun = zeros(length(models),length(beta));

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

Kiso = zeros(length(models),length(beta));
MYiso = zeros(length(models),length(beta));
Aiso = zeros(length(models),length(beta));
NUiso = zeros(length(models),length(beta));

Kre = zeros(length(models),length(beta));
MYre = zeros(length(models),length(beta));
Are = zeros(length(models),length(beta));
NUre = zeros(length(models),length(beta));

for i=1:length(models)
    Kun(i,:) = models{i}.Kun;
    MYun(i,:) = models{i}.MYun;
    Aun(i,:) = models{i}.Aun;
    NUun(i,:) = models{i}.NUun;
    Kbar(i,:) = models{i}.Kbar;
    NUbar(i,:) = models{i}.NUbar;
    Kiso(i,:) = models{i}.Kiso;
    MYiso(i,:) = models{i}.MYiso;
    Aiso(i,:) = models{i}.Aiso;
    NUiso(i,:) = models{i}.NUiso;
    Kre(i,:) = models{i}.Kre;
    MYre(i,:) = models{i}.MYre;
    NUre(i,:) = models{i}.NUre;
    Are(i,:) = models{i}.Are;
end

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

for i=1:length(models)
    leg{i} = ['\kappa = ',num2str(models{i}.kappa), ...
            ', \alpha_1 = ',num2str(models{i}.alpha1), ...
            ', \alpha_3 = ',num2str(models{i}.alpha3)];
    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
    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
    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

fig1 = figure(1);

% subplot(2,1,1)
plot(beta,Kun);
hold on;
plot(beta,MYun)
text(.001,Kun(end,1)+.01,'K_u','FontSize',14);
text(.001,MYun(end,1)+.01,'\mu_u','FontSize',14);
hold off;
xlabel('TOTAL MELT FRACTION \beta');
ylabel('MODULI [MBAR]');
title('UNRELAXED K- AND \mu-MODULI')

legend(leg,'FontSize',8);

fig2 = figure(2);

% subplot(2,1,2)
plot(beta,Kre);
hold on;
plot(beta,MYre)
text(.001,Kre(end,1)+.01,'K_r','FontSize',14);
text(.001,MYre(end,1)+.01,'\mu_r','FontSize',14);
hold off;
xlabel('TOTAL MELT FRACTION \beta');
ylabel('MODULI [MBAR]');
title('RELAXED K- AND \mu-MODULI')
set(gca,'XLim',[0 beta(end)])

legend(leg,'FontSize',8);

fig3 = figure(3);

subplot(2,1,1)
plot(beta,Kiso);
hold on;
plot(beta,MYiso)
text(.001,Kiso(end,1)+.02,'K_{iso}','FontSize',7);
text(.001,MYiso(end,1)+.03,'\mu_{iso}','FontSize',7);
hold off;
xlabel('MODIFED TOTAL MELT FRACTION \beta_m^{''}');
ylabel('MODULI [MBAR]');
title('K_{iso} - AND \mu_{iso}-MODULI')
set(gca,'XLim',[0 beta(end)])

subplot(2,1,2)
plot(beta,Kbar);
xlabel('TOTAL MELT FRACTION \beta_m^{''}');
ylabel('MODULI [MBAR]');
title('K^{''}-MODULI')
set(gca,'XLim',[0 beta(end)])

legend(leg,'FontSize',8);

fig4 = figure(4);

subplot(3,1,1)
plot(beta,Aun);
xlabel('TOTAL MELT FRACTION \beta');
ylabel('A_u');
title('A')
set(gca,'XLim',[0 beta(end)])

subplot(3,1,2)
plot(beta,Aiso);
xlabel('MODIFIED TOTAL MELT FRACTION \beta');
ylabel('A_{iso}');
set(gca,'XLim',[0 beta(end)])

subplot(3,1,3)
plot(beta,Are);
xlabel('TOTAL MELT FRACTION \beta');
ylabel('A_r');
set(gca,'XLim',[0 beta(end)])

legend(leg,'FontSize',8);

fig5 = figure(5);

subplot(4,1,1)
plot(beta,NUun);
% xlabel('TOTAL MELT FRACTION \beta');
ylabel('\nu_u');
title('POISSON-RATIO \nu')
set(gca,'XLim',[0 beta(end)])
set(gca,'XTickLabel',[])

subplot(4,1,2)
plot(beta,NUbar);
xlabel('TOTAL MELT FRACTION \beta');
ylabel('\nu^{''}');
set(gca,'XLim',[0 beta(end)])

subplot(4,1,3)
plot(beta,NUiso);
xlabel('MODIFIED TOTAL MELT FRACTION \beta');
ylabel('\nu_{iso}');
set(gca,'XLim',[0 beta(end)])

subplot(4,1,4)
plot(beta,NUre);
xlabel('TOTAL MELT FRACTION \beta');
ylabel('\nu_r');
set(gca,'XLim',[0 beta(end)])

legend(leg,'FontSize',8);