
clc; clear all; close all;  fclose('all');
%vstupy I.
Sig_0=[0; 0; 0] ;
Eps_0=[0; 0; 0] ;

%LC CT1-1: 
d_Eps_m=[2.4e-3 2.4e-3 2.4e-3 2.4e-3 2.4e-3   2.4e-3 2.4e-3 2.4e-3 2.4e-3 2.4e-3  2.4e-3 2.4e-3 2.4e-3 2.4e-3 2.4e-3   2.4e-3 2.4e-3 2.4e-3 2.4e-3 2.4e-3 ; 0 0 0 0 0   0 0 0 0 0   0 0 0 0 0   0 0 0 0 0 ; 0 0 0 0 0   0 0 0 0 0   0 0 0 0 0   0 0 0 0 0] ; %konst. prirustek delta Epsilon
%LC CT1-2: d_Eps_m=[0.8e-3 0.8e-3 0.8e-3 0.8e-3 0.8e-3   0.8e-3 0.8e-3 0.8e-3 0.8e-3 0.8e-3  0 0 0 0 0   0 0 0 0 0 ; 0 0 0 0 0   0 0 0 0 0   0 0 0 0 0   0 0 0 0 0 ; 0 0 0 0 0   0 0 0 0 0   6e-3 6e-3 6e-3 6e-3 6e-3   6e-3 6e-3 6e-3 6e-3 6e-3] ; %konst. prirustek delta Epsilon
%LC CT1-3: d_Eps_m=[3e-3 0 3e-3 2e-3 2e-3   2e-3 2e-3 2e-3 2e-3 2e-3   2e-3 2e-3 2e-3 2e-3 2e-3   2e-3 2e-3 2e-3 2e-3 2e-3;0 0 0 0 0   0 0 0 0 0   0 0 0 0 0   0 0 0 0 0 ; 0 3e-3 3e-3 3e-3 3e-3    3e-3 3e-3 3e-3 3e-3 3e-3   3e-3 3e-3 3e-3 3e-3 3e-3    2e-3 2e-3 2e-3 2e-3 2e-3] ; %konst. prirustek delta Epsilon

%LC CT1-4: d_Eps_m=[5e-3 5e-3 0 0 0   0 0 0 0 0   0 0 5e-3 5e-3 0   0 0 0 0 0;0 0 0 0 0   0 0 0 0 0   0 0 0 0 0   0 0 0 0 0 ; 0 0 2e-3 2e-3 2e-3    2e-3 2e-3 -2e-3 -2e-3 -2e-3   2e-3 2e-3 0 0 2e-3    2e-3 2e-3 2e-3 2e-3 2e-3 ] ; 

%LC CT2-1  d_Eps_m=[0 0 0 0 0  0 0 0 0 0   0 0 0 0 0  0 0 0 0 0 ; 4e-4 4e-4 4e-4 4e-4 4e-4   4e-4 4e-4 4e-4 4e-4 4e-4  4e-4 4e-4 4e-4 4e-4 4e-4   4e-4 4e-4 4e-4 4e-4 4e-4 ; 0 0 0 0 0  0 0 0 0 0   0 0 0 0 0  0 0 0 0 0 ] ; 
%LC CT2-2: d_Eps_m=[0 0 0 0 0  0 0 0 0 0   0 0 0 0 0  0 0 0 0 0 ; 3e-4 3e-4 3e-4 3e-4 3e-4   3e-4 3e-4 3e-4 3e-4 3e-4  3e-4 3e-4 3e-4 3e-4 3e-4   3e-4 3e-4 3e-4 3e-4 3e-4 ; 5e-4 5e-4 5e-4 5e-4 5e-4   5e-4 5e-4 5e-4 5e-4 5e-4  1e-4 1e-4 1e-4 1e-4 1e-4   1e-4 1e-4 1e-4 1e-4 1e-4 ] ; 
%LC CT2-3: d_Eps_m=[0 0 0 0 0  0 0 0 0 0   0 0 0 0 0  0 0 0 0 0 ; 4e-4 4e-4 0 0 0   0 0 0 0 0  0 0 4e-4 4e-4 0  0 0 0 0 0 ; 0 0 6e-4 6e-4 6e-4  6e-4 6e-4 -6e-4 -6e-4 -6e-4  -6e-4 -6e-4 0 0 6e-4   6e-4 6e-4 6e-4 6e-4 6e-4 ] ; 
%LC CT2-4: d_Eps_m=[0 0 0 0 0  0 0 0 0 0   0 0 0 0 0  0 0 0 0 0 ; 0 0 0 0 0  0 0 0 0 0  0 0 0 0 0  0 0 0 0 0 ; 4e-3 4e-3 4e-3 3e-4 3e-4  3e-4 3e-4 3e-4 3e-4 3e-4  3e-4 3e-4 3e-4 3e-4 3e-4 3e-4 3e-4 3e-4 3e-4 3e-4 ] ; 

%LC UNL: d_Eps_m=[5e-3 2e-3 2e-3 2e-3 2e-3   -2e-3 -2e-3 -2e-3 -2e-3 2e-3  2e-3 2e-3 2e-3 0 0   0 0 2e-3 2e-3 2e-3 ; 0 0 0 0 0   0 0 0 0 0   0 0 0 0 0   0 0 0 0 0 ; 0 0 0 0 0   0 0 0 0 0   0 0 0 4e-3 4e-3   4e-3 4e-3 4e-3 4e-3 4e-3] ; %konst. prirustek delta Epsilon

Ex=14956; %spruce:14956 MPa;  
Ey=3088; %spruce: 3088 MPa;
Gxy=573;%spruce: 573 MPa;
ny_xy=0.023; %spruce: 0.023; 
ny_yx=Ex*ny_xy/Ey; %for symmetric matrix Ce and De: ny_xy/Ex=ny_yx/Ey
Ce=[1/Ex -ny_xy/Ey 0;-ny_yx/Ex 1/Ey 0;0 0 1/Gxy]; %Stiffness matrix: De, Compliance matrix: Ce
De_G=inv(Ce); 

%-INPUTS II. SCM
ft1=74.4   ;%spruce:74.4 ;  %ft1 = ftx; ft2 = fty; fs=tau_xy
ft2=2.2; %spruce:2.2;
fs=6.7; %spruce: 6.7;
alfa=((ft1*ft2)/(fs^2));
k=0;

%-INPUTS III. SCM
delta_crit_11=0.45; %spruce CT1: delta_crit_11=0.45;
delta_crit_12=0.15; %spruce CT2: delta_crit_12=0.15;
l_ch_1=10; %for CT1, spruce: 10
l_ch_2=30; %for CT2, spruce: 30

%Failure surface:
    syms X Y Z real
    f_Sig=(((X-ft1)+(Y+(k*X)-ft2))/2 + sqrt( (((X-ft1-(Y+(k*X)-ft2))/2)^2) + alfa*(Z^2)));
X=Sig_0(1); Y=Sig_0(2); Z=Sig_0(3); 
f_Sig_0=subs(f_Sig);


 if f_Sig_0>0
     f_Sig_0_negat=0 %Incorrect input. The initial value of the function of plasticity is f_Sig_0 > 0. Acceptable values are f_Sig_0 < 0
     return;
 end;
 if abs(f_Sig_0)<1E-8 
     f_Sig_0_negat=0 %Incorrect input. The initial value of the function of plasticity is f_Sig_0 = 0. Acceptable values are f_Sig_0 < 0
      return;
 end;
    
 if f_Sig_0<0 
     f_Sig_0_negat=1 %correct input. The initial value of the function of plasticity is f_Sig_0 < 0. 
 end;

 
%%A]---ELASTIC CALCULATION 
f_Sig_n=f_Sig_0;
Sig_el_n=Sig_0;
StepN_el = 0;
StepN=StepN_el;
Eps_el_n=Eps_0;
%Write in Excel:
 soubor2 = fopen('SCM.csv', 'wt');
 soubor = fopen('SCM.csv', 'wt');
 fprintf(soubor, '%s\n', '1)Input - elastic calculation');
 fprintf(soubor, '%s;%s;%s;%s;%s;%s;%s;%s;%s\n', 'f_Sig_0_negat','Ex', 'Ey', 'Gxy', 'ny_xy', 'ny_yx', 'ft1', 'ft2', 'fs');
 fprintf(soubor, '%.6f;%.6f;%.6f;%.6f;%.6f;%.6f;%.6f;%.6f;%.6f\n', f_Sig_0_negat, Ex, Ey, Gxy, ny_xy, ny_yx, ft1, ft2, fs);
 fprintf(soubor, '%s;%s;%s;%s;\n', 'Ce(1,1)','Ce(1,2)','Ce(2,2)','Ce(3,3)');
 fprintf(soubor, '%.6f;%.6f;%.6f;%.6f;\n', Ce(1,1),Ce(1,2),Ce(2,2),Ce(3,3));
 fprintf(soubor, '%s;%s;%s;%s;\n', 'De(1,1)','De(1,2)','De(2,2)','De(3,3)');
 fprintf(soubor, '%.6f;%.6f;%.6f;%.6f;\n', De_G(1,1),De_G(1,2),De_G(2,2),De_G(3,3));

 fprintf(soubor, '%s\n', '2)Elastic calculation');
 fprintf(soubor, '%s;%s;%s;%s;%s;%s;%s;%s;%s;%s;%s\n', 'StepN_el', 'd_Eps(1)','d_Eps(2)', 'd_Eps(3)','Eps(1)','Eps(2)', 'Eps(3)', 'Sig_el_n(1)', 'Sig_el_n(2)', 'Sig_el_n(3)','f_Sig_n');
 fprintf(soubor, '%.6f;%s;%s;%s;%.6f;%.6f;%.6f;%.6f;%.6f;%.6f;%.6f\n', StepN_el, '-','-', '-',Eps_0(1), Eps_0(2), Eps_0(3), Sig_el_n(1), Sig_el_n(2), Sig_el_n(3),f_Sig_n);

 g_Sig_1=[];g_Sig_2=[];g_Sig_3=[];g_Eps_G_1=[];g_Eps_G_2=[];g_Eps_G_3=[];g_StepN=[];
 g_Eps_c_L_1=[];g_Eps_c_L_2=[];g_Eps_c_L_3=[];
 g_Sig_L_1=[];g_Sig_L_2=[];g_Sig_L_3=[];  g_Eps_L_1=[];g_Eps_L_2=[];g_Eps_L_3=[];

 g_Sig_1=[g_Sig_1 (0)];g_Sig_2=[g_Sig_2 (0)];g_Sig_3=[g_Sig_3 (0)];
 g_Eps_G_1=[g_Eps_G_1 (0)];g_Eps_G_2=[g_Eps_G_2 (0)];g_Eps_G_3=[g_Eps_G_3 (0)]; g_StepN=[g_StepN (StepN)];
 g_Eps_c_L_1=[g_Eps_c_L_1 (0)];g_Eps_c_L_2=[g_Eps_c_L_2 (0)];g_Eps_c_L_3=[g_Eps_c_L_3 (0)];
 g_Sig_L_1=[g_Sig_L_1 (0)];g_Sig_L_2=[g_Sig_L_2 (0)];g_Sig_L_3=[g_Sig_L_3 (0)];
 g_Eps_L_1=[g_Eps_L_1 (0)];g_Eps_L_2=[g_Eps_L_2 (0)];g_Eps_L_3=[g_Eps_L_3 (0)];

 
 
    while f_Sig_n<0
        StepN_el=StepN_el+1; %last step in elastic calculation is divided in el. and pl. part
        StepN=StepN_el;
        d_Eps=d_Eps_m(:,StepN_el);
        Sig_el_n=Sig_el_n+De_G*d_Eps; 
        X=Sig_el_n(1); Y=Sig_el_n(2); Z=Sig_el_n(3); 
        f_Sig_n=subs(f_Sig);
        Eps_el_n=Eps_el_n+d_Eps; 
        fprintf(soubor, '%f;%f;%f;%f;%f;%f;%f;%f;%f;%f;%f\n', StepN_el, d_Eps(1),d_Eps(2), d_Eps(3),Eps_el_n(1), Eps_el_n(2), Eps_el_n(3), Sig_el_n(1), Sig_el_n(2), Sig_el_n(3),f_Sig_n);
    end;
    Sig_el_n_above=Sig_el_n  
    Sig_el_n_under=Sig_el_n-De_G*d_Eps
    Eps_el_n_under=Ce*Sig_el_n_under;
    f_Sig_n_above=f_Sig_n

    X=Sig_el_n_under(1); Y=Sig_el_n_under(2); Z=Sig_el_n_under(3); 
    f_Sig_n_under=subs(f_Sig)


    Sig_el_n_above_x=Sig_el_n(1);
    Sig_el_n_above_y=Sig_el_n(2);
    Sig_el_n_above_z=Sig_el_n(3);
    Sig_el_n_under_x=Sig_el_n_under(1);
    Sig_el_n_under_y=Sig_el_n_under(2);
    Sig_el_n_under_z=Sig_el_n_under(3);

    if f_Sig_n>0
    syms t real positive
    Xp=Sig_el_n_above_x+(Sig_el_n_under_x-Sig_el_n_above_x)*t;  %Xp;Yp;Zp - coordinates of point PEAK (after crack initiation)
    Yp=Sig_el_n_above_y+(Sig_el_n_under_y-Sig_el_n_above_y)*t;
    Zp=Sig_el_n_above_z+(Sig_el_n_under_z-Sig_el_n_above_z)*t;
    f_Sig_p=(((Xp-ft1)+(Yp+(k*Xp)-ft2))/2 + sqrt( (((Xp-ft1-(Yp+(k*Xp)-ft2))/2)^2) + alfa*(Zp^2)));
    t_2reseni=solve(f_Sig_p,'t');size_t_2reseni=size(t_2reseni);  
    if size_t_2reseni(1)==2
    t_2reseni_1=eval(t_2reseni(1));t_2reseni_2=eval(t_2reseni(2));
       if (t_2reseni_1<1 && t_2reseni_1>0) || t_2reseni_1==1 || t_2reseni_1==0
         t=t_2reseni_1; %parameter t is in interval [0;1] 
       end; 
       if (t_2reseni_2<1 && t_2reseni_2>0) || t_2reseni_2==1 || t_2reseni_2==0
         t=t_2reseni_2; %parameter t is in interval [0;1]  
       end;
    end;
    if size_t_2reseni(1)==1
     t_2reseni=eval(t_2reseni(1))
       if (t_2reseni<1 && t_2reseni>0) || t_2reseni==1 || t_2reseni==0
        t=t_2reseni; %parameter t is in interval [0;1] 
       end;
    end;
    Sig_p=[subs(Xp) subs(Yp) subs(Zp)]'
    f_Sig_p=(((Sig_p(1)-ft1)+(Sig_p(2)+(k*Sig_p(1))-ft2))/2 + sqrt( (((Sig_p(1)-ft1-(Sig_p(2)+(k*Sig_p(1))-ft2))/2)^2) + alfa*(Sig_p(3)^2))); %potvrzeni ze f_Sig_p v bodu P je NULA
     end;
    
    
    Eps_p=Ce*Sig_p;
    d_Eps_step1_G=d_Eps-(Eps_p-Eps_el_n_under); % d_Eps_step1_G - 1st step of plastic calculation
    fprintf(soubor, '%.6f;%.6f;%.6f;%.6f;%.6f;%.6f;%.6f;%.6f;%.6f;%.6f;%.6f\n', StepN_el, d_Eps_step1_G(1),d_Eps_step1_G(2), d_Eps_step1_G(3),Eps_p(1), Eps_p(2), Eps_p(3), Sig_p(1), Sig_p(2), Sig_p(3),f_Sig_p);


    X=Sig_p(1);
    Y=Sig_p(2);
    Zp=Sig_p(3);
    theta=8; %[°]
    Z_f_CT=(tan(2*(theta/180*pi)))*((X-Y)/2); %(theta/180*pi) [rad]
    fi=0.5*atan(2*Sig_p(3)/(Sig_p(1)-Sig_p(2)));  %fi-deviation of principal stress and grain direction
    fi_deg=fi/pi*180 ;  
    if Zp>Z_f_CT
        MODE=2 %MODE=2 - it means CT2
        n=[0;1]; m=[-1;0];
        delta_crit_1=delta_crit_12;
        l_ch=l_ch_2;
        fi_deg=90;fi=fi_deg/180*pi;
    end;
    if Zp<=Z_f_CT
        MODE=1  %MODE=1 - it means CT1
       n=[cos(fi); sin(fi)]
       m=[-n(2); n(1)];
       fi_deg=fi/pi*180;
        delta_crit_1=delta_crit_11;
        l_ch=l_ch_1;
    end;
    ec_crit_1=delta_crit_1/l_ch;

    %Write in Excel:
     fprintf(soubor, '%s\n', '3)Output - elastic calculation');
     fprintf(soubor, '%s;%s;%s;%s;%s;%s;%s;%s;%s;%s ;%s;%s;%s;%s;%s;%s\n', 'MODE', 'f_Sig_p', 'Sig_p(1)','Sig_p(2)', 'Sig_p(3)', 'Eps_p_G(1)','Eps_p_G(2)','Eps_p_G(3)', 'd_Eps_step1_G(1)','d_Eps_step1_G(2)','d_Eps_step1_G(3)', 'fi_deg', 'n(1)', 'n(2)', 'm(1)', 'm(2)');
     fprintf(soubor, '%.6f;%.6f;%.6f;%.6f;%.6f;%.6f;%.6f;%.6f;%.6f;%.6f ;%.6f;%.6f;%.6f;%.6f;%.6f;%.6f\n', MODE, f_Sig_p, Sig_p(1), Sig_p(2), Sig_p(3), Eps_p(1),Eps_p(2), Eps_p(3),d_Eps_step1_G(1),d_Eps_step1_G(2),d_Eps_step1_G(3), fi_deg, n(1), n(2), m(1), m(2));
    %Write in vectors for graphs at the end:
     g_Sig_1=[g_Sig_1 (Sig_p(1))];g_Sig_2=[g_Sig_2 (Sig_p(2))];g_Sig_3=[g_Sig_3 (Sig_p(3))];g_Eps_G_1=[g_Eps_G_1 (Eps_p(1))];g_Eps_G_2=[g_Eps_G_2 (Eps_p(2))];g_Eps_G_3=[g_Eps_G_3 (Eps_p(3))];g_StepN=[g_StepN (StepN)]
  

%--- B] SMEARED CRACK MODEL (SCM):
%- INPUTS: n, m, Eps_p, Sig_p,d_Eps_step1_G 
 Sig_2ord_p=[Sig_p(1) Sig_p(3);Sig_p(3) Sig_p(2)];
 Eps_p_G=Eps_p;
  
 %-- G->L: 
       R_De=[(cos(fi))^2 (sin(fi))^2 -sin(fi)*cos(fi); (sin(fi))^2 (cos(fi))^2 sin(fi)*cos(fi); 2*cos(fi)*sin(fi) 2*(-sin(fi))*cos(fi) (cos(fi))^2-(sin(fi))^2];
       De_L=R_De'*De_G*R_De; 
       R_Eps_Sig=[cos(fi) -sin(fi); sin(fi) cos(fi)];   %rotation matrix countercloclwise for Epsilon or Sigma
       d_Eps_2ord_step1_L=R_Eps_Sig*[d_Eps_step1_G(1) d_Eps_step1_G(3); d_Eps_step1_G(3) d_Eps_step1_G(2)]*R_Eps_Sig';

       d_Eps_step1_L=[d_Eps_2ord_step1_L(1,1); d_Eps_2ord_step1_L(2,2); d_Eps_2ord_step1_L(1,2) ];  

       if MODE==1  
       Sig_p_hl_L_max=max( (Sig_p(1)+Sig_p(2))/2+sqrt( ((Sig_p(1)-Sig_p(2))/2)^2+Sig_p(3)^2), (Sig_p(1)+Sig_p(2))/2-sqrt( ((Sig_p(1)-Sig_p(2))/2)^2+Sig_p(3)^2) );
       Sig_nn_p_scalar=Sig_p_hl_L_max;
       Sig_nm_p_scalar=0;
     else
       tc_p_G=n'*Sig_2ord_p; % tc_p_G=tensorDot2D(n,Sig_2ord_p);  
       tc_p_L=[tc_p_G(2);tc_p_G(1)]
       Sig_nn_p_scalar=tc_p_L(1);
       Sig_nm_p_scalar=tc_p_L(2);
     end;  
       
       %-- Eps_2ord_p, Sig_2ord_p: G->L
       Eps_2ord_p_L=R_Eps_Sig*[Eps_p_G(1) Eps_p_G(3); Eps_p_G(3) Eps_p_G(2)]*R_Eps_Sig';
       Eps_p_L=[Eps_2ord_p_L(1,1); Eps_2ord_p_L(2,2); Eps_2ord_p_L(1,2) ];
       Sig_p_2ord_L=R_Eps_Sig*Sig_2ord_p*R_Eps_Sig';
       Sig_p_L=[Sig_p_2ord_L(1,1); Sig_p_2ord_L(2,2); Sig_p_2ord_L(1,2) ];
       
 %--transformation matrix N_L 
     N_L=[1 0 0 ; 0 0 1];
 %--Qe Acoustic matrix:
     Qe_2ord = N_L*De_L*N_L';  
 %-- Dc_2ord:
     c1=3;c2=6.93; %parameters of function t_n
     syms Eps_c_nn Gama_c_nm real
     t_n=Sig_nn_p_scalar*( (1+((c1*l_ch*Eps_c_nn/delta_crit_1)^3)) *exp(-c2*l_ch*Eps_c_nn/delta_crit_1) - exp(-c2)*(1+(c1^3))*l_ch*Eps_c_nn/delta_crit_1);  
     
     %-- ec_0 for defined Sig_nm_p_scalar=t_m
     t_m=fs*atan(abs(Gama_c_nm)*((1-((Eps_c_nn*l_ch)/delta_crit_1)^2)/(1-(1-((Eps_c_nn*l_ch)/delta_crit_1)^2))))/(pi/2);
     if Sig_nm_p_scalar==0; Sig_nm_p_scalar_q=(delta_crit_1/l_ch/100); else Sig_nm_p_scalar_q=Sig_nm_p_scalar; end;
     if (abs(Sig_nm_p_scalar-fs)<=1e-3);Sig_nm_p_scalar_q=fs-1e-5;end;
     q=solve(t_m-Sig_nm_p_scalar_q,Eps_c_nn,Gama_c_nm);
        if MODE==2; x=-ec_crit_1/1e9;else x=ec_crit_1/1e20;end;
        s_Eps_c_nn_obec=(eval(q.Eps_c_nn))
        Eps_c_nn=max(subs(s_Eps_c_nn_obec))
        Gama_c_nm=x
        ec_0=[Eps_c_nn;Gama_c_nm]
     %---
     syms Eps_c_nn Gama_c_nm real
     if MODE==2
          t_m=-fs*atan((-Gama_c_nm)*((1-((Eps_c_nn*l_ch)/delta_crit_1)^2)/(1-(1-((Eps_c_nn*l_ch)/delta_crit_1)^2))))/(pi/2);
     else t_m=fs*atan( (Gama_c_nm)*((1-((Eps_c_nn*l_ch)/delta_crit_1)^2)/(1-(1-((Eps_c_nn*l_ch)/delta_crit_1)^2))))/(pi/2);
     end;
    Dc_2ord=[diff(t_n,Eps_c_nn), diff(t_n,Gama_c_nm); diff(t_m,Eps_c_nn), diff(t_m,Gama_c_nm)];
    Dc_2ord_00=[0, 0; 0, 0];
    fc_ec=[t_n;t_m]; 
 %--derivative of s=[fc(ec)+Qe*ec]
 syms Eps_c_nn Gama_c_nm real
 Qe_times_ec=Qe_2ord*[Eps_c_nn; Gama_c_nm];
 diff_Qe_times_ec=[diff(Qe_times_ec(1),Eps_c_nn), diff(Qe_times_ec(1),Gama_c_nm);diff(Qe_times_ec(2),Eps_c_nn), diff(Qe_times_ec(2),Gama_c_nm)];
 s_slope=diff_Qe_times_ec+Dc_2ord;
 Eps_c_nn=0; Gama_c_nm=0;
 s_slope_eval=subs(s_slope);
     if s_slope_eval(1,1)<0;s_slope_11_positive=0
     else s_slope_11_positive=1; end;
     if s_slope_eval(1,2)<0;s_slope_12_positive=0
     else s_slope_12_positive=1; end;
     if s_slope_eval(2,2)<0;s_slope_22_positive=0
     else s_slope_22_positive=1;end; 
 
 %--pro vypocet Dec_L: De_L*N_L' a N_L*De_L
 De_L_times_N_L_T=De_L*N_L'; 
 N_L_times_De_L=N_L*De_L;
 
% - STEP N (2D):
 d_Eps_StepN_L=d_Eps_step1_L; 
 d_Eps_StepN_G=d_Eps_step1_G;
 d_Eps_iter_L=d_Eps_step1_L; 
 s_previous=N_L*De_L*Eps_p_L;
 d_s_iter=N_L*De_L*d_Eps_iter_L; 

 
  ec=ec_0; Eps_c_nn=ec(1); Gama_c_nm=ec(2);
       Dec_StepN_L=De_L-De_L*N_L'*inv(Qe_2ord+subs(Dc_2ord))*N_L*De_L; % Dec(ec_0)
       Dec_StepN_G=inv(R_De')*Dec_StepN_L*R_De';
 delta=[0;0];
 fc_ec_StepN=subs(fc_ec);
 fc_ec_plus_Qe_ec=fc_ec_StepN+Qe_2ord*ec;
 s_StepN=N_L*De_L*(Eps_p_L+d_Eps_StepN_L) ;
 d_s_StepN=N_L*De_L*d_Eps_StepN_L ;
 d_ec_StepN=[0;0];
 Sig_StepN=Sig_p;   
 d_Sig_StepN_L= [0; 0; 0]; 
 Sig_StepN_L=[0; 0; 0]; 
 Eps_StepN_G=Eps_p_G;
 Eps_StepN_L=Eps_p_L;
 d_Eps_StepN_L=d_Eps_step1_L;
 d_Eps_c_L=[0; 0; 0]; 
 Eps_c_L=[0; 0 ;0];
 iKT=inv(Qe_2ord+subs(Dc_2ord)); 
 StepN_pl=0;
 StepN=StepN_el;
    fprintf(soubor, '%s\n', '4)Input - FCM calculation');
    fprintf(soubor, '%s;%s;%s;%s;%s;%s;%s\n', 'delta_crit_1', 'Sig_nn_p_scalar', 'Sig_nm_p_scalar',  'l_ch','s_slope_11_pos','s_slope_12_pos','s_slope_22_pos' );
 	fprintf(soubor, '%.6f;%.6f;%.6f;%f;%f;%f;%f\n', delta_crit_1, Sig_nn_p_scalar, Sig_nm_p_scalar, l_ch,s_slope_11_positive,s_slope_12_positive,s_slope_22_positive );
    
    fprintf(soubor, '%s\n', '5)Output - FCM calculation');
    fprintf(soubor, '%s;%s;%s;%s;%s; %s;%s;%s;%s;%s; %s;%s;%s;%s;%s; %s;%s;%s;%s;%s; %s;%s;%s;%s;%s; %s;%s;%s;%s;%s; %s;%s;%s;%s;%s; %s;%s;%s;%s;%s; %s\n',  'StepN_pl','d_Eps_L(1)','d_Eps_L(2)', 'd_Eps_L(3)','d_Eps_G(1)','d_Eps_G(2)', 'd_Eps_G(3)','Eps_G(1)','Eps_G(2)', 'Eps_G(3)','d_ec(1)', 'd_ec(2)', 'ec(1)', 'ec(2)', 'delta(1)', 'delta(2)','fc_ec(1)' , 'fc_ec(2)','d_Sig(1)', 'd_Sig(2)', 'd_Sig(3)', 'Sig(1)', 'Sig(2)',  'SigG(3)', 'Dec_G(1,1)', 'Dec_G(1,2)', 'Dec_G(2,2)', 'Dec_G(3,3)','Dec_L(1,1)','Dec_L(1,2)','Dec_L(2,2)','Dec_L(3,3)', 's_StepN(1)','s_StepN(2)','d_Eps_c_L(1)','d_Eps_c_L(2)','d_Eps_c_L(3)', 'iKT(1,1)', 'iKT(1,2)', 'iKT(2,2)','StepN' );
 	fprintf(soubor, '%f;%.6s;%.6s;%.6s;%.6s;  %.6s;%.6s;%.6f;%.6f;%.6f;  %.6f;%.6f;%.6f;%.6f;%.6f;  %.6f;%.6f;%.6f;%s;%s;  %s;%.6f;%.6f ;%.6f;%s;  %s;%s;%s;%s;%s;  %s;%s;%.6f;%.6f;%.6f;  %.6f;%.6f;%s;%s;%s;  %f\n', StepN_pl, '-','-', '-','-','-', '-',Eps_StepN_G(1),Eps_StepN_G(2), Eps_StepN_G(3), d_ec_StepN(1), d_ec_StepN(2),ec(1), ec(2),delta(1), delta(2),fc_ec_StepN(1) ,fc_ec_StepN(2),'-', '-', '-', Sig_StepN(1), Sig_StepN(2), Sig_StepN(3), '-', '-', '-', '-', '-', '-', '-', '-',s_previous(1), s_previous(2), d_Eps_c_L(1),d_Eps_c_L(2),d_Eps_c_L(3),'-', '-', '-', StepN );
 
    fprintf(soubor2, '%s\n', '6)Plastic iterations within a StepN_pl');
    fprintf(soubor2, '%s;%s;%s;%s;%s;%s;%s;%s;%s\n', 'StepN_pl', 'No_pl_iter', 's_StepN(1)', 's_StepN(2)','s(1)','s(2)','ec(1)','ec(2)', 'StepN'); 
    fprintf(soubor2, '%f;%s;%s;%s;%.6f;%.6f;%.6f;%.6f;%f\n',  StepN_pl, '-', '-', '-',s_previous(1),s_previous(2),ec(1),ec(2),StepN);
    
    g_fc_1=[];g_fc_2=[];g_ec_1=[];g_ec_2=[];g_s_1=[];g_s_2=[];g_delta_1=[];g_delta_2=[];g_Dec_G_11=[];g_Dec_G_22=[];g_Dec_G_33=[];g_Dec_G_12=[];
    
    g_fc_1=[g_fc_1 (fc_ec_StepN(1))];g_fc_2=[g_fc_2 (fc_ec_StepN(2))];g_ec_1=[g_ec_1 (ec(1))];g_ec_2=[g_ec_2 (ec(2))];  g_s_1=[g_s_1 (s_previous(1))];g_s_2=[g_s_2 (s_previous(2))];g_delta_1=[g_delta_1 (ec(1)*l_ch)];g_delta_2=[g_delta_2 (ec(2)*l_ch)];
    g_Eps_c_L_1=[g_Eps_c_L_1 (Eps_c_L(1))];g_Eps_c_L_2=[g_Eps_c_L_2 (Eps_c_L(2))];g_Eps_c_L_3=[g_Eps_c_L_3 (Eps_c_L(3))];
    g_Sig_L_1=[g_Sig_L_1 (Sig_p_L(1))];g_Sig_L_2=[g_Sig_L_2 (Sig_p_L(2))];g_Sig_L_3=[g_Sig_L_3 (Sig_p_L(3))];
    g_Eps_L_1=[g_Eps_L_1 (Eps_p_L(1))];g_Eps_L_2=[g_Eps_L_2 (Eps_p_L(2))];g_Eps_L_3=[g_Eps_L_3 (Eps_p_L(3))];
     
  UNL=0; Kapa=[0; 0];  %Kapa-internal variable, sum of positive d_ec
  fc_ec_S_1=0;
     
  StepN_pl=1;
  Cond=1;
 
  while StepN<=20    
      
    if ( (abs(ec(1))<ec_crit_1) && (abs(ec(2))<ec_crit_1) && (UNL==0) )  % && (fc_ec_StepN(2)<fs))    
          No_pl_iter=1;
          while Cond > 1e-4
            d_ec_iter=(inv(Qe_2ord+subs(Dc_2ord)))*d_s_iter;
            ec=ec+d_ec_iter;
               Eps_c_nn=ec(1); Gama_c_nm=ec(2);
            d_ec_StepN=d_ec_StepN+d_ec_iter
            s=subs(fc_ec)+Qe_2ord*ec;
            d_s_iter=(s_StepN-s);
            Cond=( sqrt(d_s_iter(1)^2+d_s_iter(2)^2)/sqrt(s_previous(1)^2+s_previous(2)^2) ) ;

            fprintf(soubor2, '%f;%f;%.6f;%.6f;%.6f;%.6f;%.6f;%.6f;%f\n',  StepN_pl, No_pl_iter, s_StepN(1), s_StepN(2),s(1),s(2),ec(1),ec(2),StepN);
            No_pl_iter=No_pl_iter + 1;
          end;
    end;
     
    if ( (abs(ec(1))<ec_crit_1) && (abs(ec(2))<ec_crit_1) && UNL==1)  % && (fc_ec_StepN(2)<fs))    
          No_pl_iter=1;
          
          while Cond > 1e-4
            d_ec_iter=(inv(Qe_2ord+subs(Dc_2ord_unl)))*d_s_iter;
            ec=ec+d_ec_iter;
               Eps_c_nn=ec(1); Gama_c_nm=ec(2);
            d_ec_StepN=d_ec_StepN+d_ec_iter;
            s=subs(fc_ec_unl)+Qe_2ord*ec;
            d_s_iter=(s_StepN-s);
            Cond=( sqrt(d_s_iter(1)^2+d_s_iter(2)^2)/sqrt(s_previous(1)^2+s_previous(2)^2) ) ;

            fprintf(soubor2, '%f;%f;%.6f;%.6f;%.6f;%.6f;%.6f;%.6f;%f\n',  StepN_pl, No_pl_iter, s_StepN(1), s_StepN(2),s(1),s(2),ec(1),ec(2),StepN);
            No_pl_iter=No_pl_iter + 1;
          end;
    end;
    
    %--- unloading UNL:
   
    if ( (d_ec_StepN(1)>-1e-8 && d_ec_StepN(2)>-1e-8 && UNL==0) || (((d_ec_StepN(1))>-1e-8 && d_ec_StepN(2)>-1e-8 && ec(1)>ec_S(1) && ec(2)>ec_S(2))) ); %nacitani Kapa
        Kapa=Kapa+d_ec_StepN ;
    end;

    disp(StepN);disp(UNL); disp(ec); disp(d_ec_StepN);  disp(Kapa); 
    
    if ( UNL==1 && ((ec(1)-Kapa(1))<1e-8) && ((ec(2)-Kapa(2))<1e-8) &&   (Kapa(1)-ec(1)>d_ec_StepN(1))  && abs(Kapa(2)-ec(2)-d_ec_StepN(2))<1e-5 )  %krok 2->(N-1) 
        disp('Step_UNL 2->(N-1)'); UNL==1; 
    end; 
    
    if ( UNL==1 && (d_ec_StepN(1)>-1e-8) && (d_ec_StepN(2)>-1e-8) &&  (abs(Kapa(1)-ec(1))<1e-6)  && (abs(Kapa(2)-ec(2))<1e-6) ) 
          disp('Step_UNL N'); UNL=0 
    end;   
    
    if ( (UNL==0) && (ec(1)-Kapa(1))<1e-8 && (abs(ec(2)-Kapa(2)))<1e-8 && (abs(Kapa(1)-(ec(1)-d_ec_StepN(1)))<1e-8) && (abs(Kapa(2)-(ec(2)-d_ec_StepN(2))<1e-8)) ) %1.krok odtezovani
        disp('Step_UNL 1') 
        ec_S=ec-d_ec_StepN
          Eps_c_nn=ec_S(1); Gama_c_nm=ec_S(2);
        fc_ec_S=(subs(fc_ec)) ; 
        fc_ec_S_1=fc_ec_S(1);ec_S_1=ec_S(1);
       %--data pro iterace:
        syms Eps_c_nn Gama_c_nm real
        t_n_unl=fc_ec_S(1)/ec_S(1)*Eps_c_nn
        t_m=fs*atan( (Gama_c_nm)*((1-((Eps_c_nn*l_ch)/delta_crit_1)^2)/(1-(1-((Eps_c_nn*l_ch)/delta_crit_1)^2))))/(pi/2);
        fc_ec_unl=[t_n_unl;t_m];
        Dc_2ord_unl=[(fc_ec_S(1)/ec_S(1)) 0; Dc_2ord(2,1), Dc_2ord(2,2)];
        
        Eps_c_nn=ec_S(1); Gama_c_nm=ec_S(2);
        %--iter.:
          while Cond > 1e-8
            d_ec_iter=(inv(Qe_2ord+subs(Dc_2ord_unl)))*d_s_iter;
            ec=ec+d_ec_iter
               Eps_c_nn=ec(1); Gama_c_nm=ec(2);
            d_ec_StepN=d_ec_StepN+d_ec_iter;
            s=subs(fc_ec_unl)+Qe_2ord*ec;
            d_s_iter=(s_StepN-s);
            Cond=( sqrt(d_s_iter(1)^2+d_s_iter(2)^2)/sqrt(s_previous(1)^2+s_previous(2)^2) ) ;

            fprintf(soubor2, '%f;%f;%.6f;%.6f;%.6f;%.6f;%.6f;%.6f;%f\n',  StepN_pl, No_pl_iter, s_StepN(1), s_StepN(2),s(1),s(2),ec(1),ec(2),StepN);
            No_pl_iter=No_pl_iter + 1;
          end;
           disp(ec); disp(d_ec_StepN);  
        UNL=1
    end;
    

    if ( (abs(ec(1))<ec_crit_1) && (abs(ec(2))<ec_crit_1)  )  
      d_Eps_c_L=N_L'*d_ec_StepN;
      Eps_c_L=Eps_c_L+d_Eps_c_L;
      
      d_Sig_StepN_L=De_L*(d_Eps_StepN_L-d_Eps_c_L);
      Sig_StepN_L=Sig_StepN_L+d_Sig_StepN_L;
      d_Sig_2ord_StepN_L=[d_Sig_StepN_L(1) d_Sig_StepN_L(3); d_Sig_StepN_L(3) d_Sig_StepN_L(2)];
      d_Sig_2ord_StepN_G=inv(R_Eps_Sig)*d_Sig_2ord_StepN_L*inv(R_Eps_Sig'); %rotace po smeru hodin. rucicek L->G
      d_Sig_StepN=[d_Sig_2ord_StepN_G(1,1); d_Sig_2ord_StepN_G(2,2); d_Sig_2ord_StepN_G(1,2)];
      Sig_StepN=Sig_StepN+d_Sig_StepN
      if UNL==1  ; fc_ec_StepN=subs(fc_ec_unl); 
      else fc_ec_StepN=subs(fc_ec);
      end; 

    else
      fc_ec_StepN=[0;0];

      d_Eps_c_L=N_L'*d_ec_StepN;
      Eps_c_L=Eps_c_L+d_Eps_c_L;
      d_Sig_StepN_L=[0;0;0];
      Sig_StepN_L=[0;0;0];
      d_Sig_StepN=[0;0;0];
      Sig_StepN=[0;0;0];
    end;

      delta=ec*l_ch;  
      Eps_StepN_G=Eps_StepN_G+d_Eps_StepN_G;
      Eps_StepN_L=Eps_StepN_L+d_Eps_StepN_L;
      
      fprintf(soubor, '%f;%.6f;%.6f;%.6f;%.6f;%.6f;%.6f;%.6f;%.6f;%.6f; %.6f;%.6f;%.6f;%.6f;%.6f;%.6f;%.6f;%.6f;%.6f;%.6f;  %.6f;%.6f;%.6f;%.6f;%.6f;%.6f;%.6f;%.6f;%.6f;%.6f;  %.6f;%.6f;%.6f;%.6f;%.6f;%.6f;%.6f;%.6f;%.6f;%.6f; %f \n', StepN_pl, d_Eps_StepN_L(1),d_Eps_StepN_L(2), d_Eps_StepN_L(3),d_Eps_StepN_G(1),d_Eps_StepN_G(2), d_Eps_StepN_G(3),Eps_StepN_G(1),Eps_StepN_G(2), Eps_StepN_G(3),d_ec_StepN(1), d_ec_StepN(2),ec(1),  ec(2),delta(1), delta(2),fc_ec_StepN(1) ,fc_ec_StepN(2), d_Sig_StepN(1), d_Sig_StepN(2), d_Sig_StepN(3), Sig_StepN(1), Sig_StepN(2), Sig_StepN(3), Dec_StepN_G(1,1), Dec_StepN_G(1,2), Dec_StepN_G(2,2), Dec_StepN_G(3,3), Dec_StepN_L(1,1),Dec_StepN_L(1,2),Dec_StepN_L(2,2),Dec_StepN_L(3,3), s_StepN(1),s_StepN(2) , d_Eps_c_L(1),d_Eps_c_L(2),d_Eps_c_L(3),iKT(1,1), iKT(1,2), iKT(2,2),StepN );
  
      g_Sig_1=[g_Sig_1 (Sig_StepN(1))];g_Sig_2=[g_Sig_2 (Sig_StepN(2))];g_Sig_3=[g_Sig_3 (Sig_StepN(3))];
      g_Eps_G_1=[g_Eps_G_1 (Eps_StepN_G(1))];g_Eps_G_2=[g_Eps_G_2 (Eps_StepN_G(2))];g_Eps_G_3=[g_Eps_G_3 (Eps_StepN_G(3))]; 
      g_StepN=[g_StepN (StepN)];
      g_fc_1=[g_fc_1 (fc_ec_StepN(1))];g_fc_2=[g_fc_2 (fc_ec_StepN(2))];g_ec_1=[g_ec_1 (ec(1))];g_ec_2=[g_ec_2 (ec(2))];  g_s_1=[g_s_1 (s_StepN(1))];g_s_2=[g_s_2 (s_StepN(2))];g_delta_1=[g_delta_1 (ec(1)*l_ch)];g_delta_2=[g_delta_2 (ec(2)*l_ch)];
      g_Dec_G_11=[g_Dec_G_11 (Dec_StepN_G(1,1))];g_Dec_G_22=[g_Dec_G_22 (Dec_StepN_G(2,2))];g_Dec_G_33=[g_Dec_G_33 (Dec_StepN_G(3,3))];g_Dec_G_12=[g_Dec_G_12 (Dec_StepN_G(1,2))];
      g_Eps_c_L_1=[g_Eps_c_L_1 (Eps_c_L(1))];g_Eps_c_L_2=[g_Eps_c_L_2 (Eps_c_L(2))];g_Eps_c_L_3=[g_Eps_c_L_3 (Eps_c_L(3))];
      g_Sig_L_1=[g_Sig_L_1 (Sig_StepN_L(1))];g_Sig_L_2=[g_Sig_L_2 (Sig_StepN_L(2))];g_Sig_L_3=[g_Sig_L_3 (Sig_StepN_L(3))];
      g_Eps_L_1=[g_Eps_L_1 (Eps_StepN_L(1))];g_Eps_L_2=[g_Eps_L_2 (Eps_StepN_L(2))];g_Eps_L_3=[g_Eps_L_3 (Eps_StepN_L(3))];
            
      if StepN==20
          break;
      end;    
        
      
      %---For the next STEP:
            if ( (abs(ec(1))>=ec_crit_1) || (abs(ec(2))>=ec_crit_1) )
               No_pl_iter=0;
               Dc_2ord=Dc_2ord_00; fc_ec=[0;0];
               iKT=inv(Qe_2ord+Dc_2ord_00);
               Dec_StepN_L=[0 0 0;0 0 0;0 0 0]; 
               Dec_StepN_G=[0 0 0;0 0 0;0 0 0];
            else
               iKT=inv(Qe_2ord+subs(Dc_2ord));
               Dec_StepN_L=De_L-De_L*N_L'*inv(Qe_2ord+subs(Dc_2ord))*N_L*De_L;
               Dec_StepN_G=inv(R_De')*Dec_StepN_L*R_De';
            end;
         StepN_pl=StepN_pl + 1;
         StepN=StepN+1; 
         s_previous=s_previous+d_s_StepN;
         d_Eps_StepN_G=[d_Eps_m(1,StepN); d_Eps_m(2,StepN); d_Eps_m(3,StepN)];
         d_Eps_2ord_L=R_Eps_Sig*[d_Eps_StepN_G(1) d_Eps_StepN_G(3); d_Eps_StepN_G(3) d_Eps_StepN_G(2)]*R_Eps_Sig';
         d_Eps_StepN_L=[d_Eps_2ord_L(1,1); d_Eps_2ord_L(2,2); d_Eps_2ord_L(1,2) ];        
         d_Eps_iter_L=d_Eps_StepN_L;       
         d_s_StepN=N_L*De_L*d_Eps_StepN_L ;
         s_StepN=s_StepN+d_s_StepN;
         d_s_iter=N_L*De_L*d_Eps_iter_L; 
         d_ec_StepN=[0;0];
         Cond=1; 
      
     if ( (abs(ec(1))>=ec_crit_1) || (abs(ec(2))>=ec_crit_1) )
           No_pl_iter=1;
           Dc_2ord=Dc_2ord_00; fc_ec=[0;0];
           iKT=inv(Qe_2ord+Dc_2ord_00);
                d_ec_StepN=iKT*d_s_StepN;
                ec=ec+d_ec_StepN;
                Eps_c_nn=ec(1); Gama_c_nm=ec(2);
                s=subs(fc_ec)+Qe_2ord*ec;

            fprintf(soubor2, '%f;%f;%.6f;%.6f;%.6f;%.6f;%.6f;%.6f;%f\n',  StepN_pl, No_pl_iter, s_StepN(1), s_StepN(2),s_StepN(1),s_StepN(2),ec(1),ec(2),StepN);         
     end;
  end;



figure %1
     subplot(2,2,1); hold on
      plot(1:length(g_Eps_G_1),g_Eps_G_1,'-xb'), plot(1:length(g_Eps_G_2),g_Eps_G_2,':ob'), plot(1:length(g_Eps_G_3),g_Eps_G_3,'--.b');
      xlabel('StepN [-]'); ylabel('\epsilon_x ; \epsilon_y ; \gamma_x_y [-]'); legend('\epsilon_x','\epsilon_y','\gamma_x_y');
     subplot(2,2,2); hold on
      plot(g_Eps_G_1,g_Sig_1,'-xb');
      xlabel('\epsilon_x [-]'); ylabel('\sigma_x [MPa]'); 
     subplot(2,2,3); hold on
      plot(g_Eps_G_1,g_Sig_2,'-xb');
      xlabel('\epsilon_x [-]'); ylabel('\sigma_y [MPa]'); 
     subplot(2,2,4); hold on
      plot(g_Eps_G_1,g_Sig_3,'-xb');
      xlabel('\epsilon_x [-]'); ylabel('\tau_x_y [MPa]'); 
   figure%2   
     subplot(2,2,1); hold on
      plot(g_Eps_G_2,g_Sig_1,'-xb');
      xlabel('\epsilon_y [-]'); ylabel('\sigma_x [MPa]');  
     subplot(2,2,2); hold on
      plot(g_Eps_G_2,g_Sig_2,'-xb');
      xlabel('\epsilon_y [-]'); ylabel('\sigma_y [MPa]');  
     subplot(2,2,3); hold on
      plot(g_Eps_G_2,g_Sig_3,'-xb');
      xlabel('\epsilon_y [-]'); ylabel('\tau_x_y [MPa]');  
   figure%3  
     subplot(2,2,1); hold on;
      plot(g_Eps_G_3,g_Sig_1,'-xb');
      xlabel('\gamma_x_y [-]'); ylabel('\sigma_x [MPa]');  
     subplot(2,2,2); hold on;
      plot(g_Eps_G_3,g_Sig_2,'-xb');
      xlabel('\gamma_x_y [-]'); ylabel('\sigma_y [MPa]');  
     subplot(2,2,3); hold on;
      plot(g_Eps_G_3,g_Sig_3,'-xb');
      xlabel('\gamma_x_y [-]'); ylabel('\tau_x_y [MPa]');     
  
  figure%4 
  if MODE==1
     subplot(2,2,1); hold on
      X=g_ec_1; Y=g_ec_2; Z=g_fc_1;
      plot3(X, Y, Z,'-db','MarkerSize',7,'MarkerFaceColor','c','MarkerEdgeColor','k');
      [X,Y]=meshgrid(0:(ec_crit_1/15):ec_crit_1);
      Z=Sig_nn_p_scalar.*( (1+((c1.*l_ch.*X./delta_crit_1).^3)) .*exp(-c2.*l_ch.*X./delta_crit_1) - exp(-c2).*(1+(c1.^3)).*l_ch.*X./delta_crit_1);    
      mesh(X,Y,Z);
      [X,Y]=meshgrid(ec_crit_1:(ec_crit_1/15):1.2*ec_crit_1, 0:(ec_crit_1/15):ec_crit_1);Z=(X-X)+(Y-Y) ; mesh(X,Y,Z);
      [X,Y]=meshgrid(0:(ec_crit_1/15):1.2*ec_crit_1, ec_crit_1:(ec_crit_1/15):1.2*ec_crit_1 );Z=(X-X)+(Y-Y) ; mesh(X,Y,Z);

      if fc_ec_S_1>0
          [X,Y]=meshgrid(0:(ec_crit_1/10):ec_crit_1/4,0:(ec_crit_1/10):ec_crit_1 );

          Z=fc_ec_S_1/ec_S_1.*X;  %t_n_unl
          mesh(X,Y,Z);
      end;
      xlabel('\epsilon^c_n_n [-]');ylabel('\gamma^c_n_m [-]'); zlabel('t_n [MPa] '); legend('t_n_,_i_t_e_r','t_n_,_p_r_e_c_i_s_e');
     subplot(2,2,2); hold on
      X=g_ec_1; Y=g_ec_2; Z=g_fc_2;
      plot3(X, Y, Z,'-db','MarkerSize',7,'MarkerFaceColor','c','MarkerEdgeColor','k');
      [X,Y]=meshgrid(0:(ec_crit_1/15):ec_crit_1); 
      Z=fs.*atan(Y.*((1-((X.*l_ch)./delta_crit_1).^2)./(1-(1-((X.*l_ch)./delta_crit_1).^2))))./(pi./2);
          mesh(X,Y,Z);
      [X,Y]=meshgrid(ec_crit_1:(ec_crit_1/15):1.2*ec_crit_1, 0:(ec_crit_1/15):ec_crit_1);Z=(X-X)+(Y-Y) ; mesh(X,Y,Z);
      [X,Y]=meshgrid(0:(ec_crit_1/15):1.2*ec_crit_1, ec_crit_1:(ec_crit_1/15):1.2*ec_crit_1 );Z=(X-X)+(Y-Y) ; mesh(X,Y,Z);
          xlabel('\epsilon^c_n_n [-]');ylabel('\gamma^c_n_m [-]'); zlabel('t_m [MPa] '); legend('t_m_,_i_t_e_r','t_m_,_p_r_e_c_i_s_e');      
     subplot(2,2,3); hold on
      X=g_ec_1; Y=g_ec_2; Z=g_s_1;
      plot3(X,Y,Z,'-db','MarkerSize',7,'MarkerFaceColor','c','MarkerEdgeColor','k');  
      [X,Y]=meshgrid(0:(ec_crit_1/15):ec_crit_1); 
      t_n=Sig_nn_p_scalar.*( (1+((c1.*l_ch.*X./delta_crit_1).^3)) .*exp(-c2.*l_ch.*X./delta_crit_1) - exp(-c2).*(1+(c1.^3)).*l_ch.*X./delta_crit_1);  %.*exp(-c2.*l_ch.*Eps_c_nn./delta_crit_1)    
      Qe_times_ec_1=Qe_2ord(1,1).*X+Qe_2ord(1,2).*Y;
      Z=t_n+Qe_times_ec_1;
      mesh(X,Y,Z);
      xlabel('\epsilon^c_n_n [-]');ylabel('\gamma^c_n_m [-]'); zlabel('s(1) [-]'); legend('s_i_t_e_r(1)','s_p_r_e_c_i_s_e(1)');      
     subplot(2,2,4); hold on
      X=g_ec_1; Y=g_ec_2; Z=g_s_2;
      plot3(X,Y,Z,'-db','MarkerSize',7,'MarkerFaceColor','c','MarkerEdgeColor','k');  
      [X,Y]=meshgrid(0:(ec_crit_1/15):ec_crit_1); 
      t_m=fs.*atan(Y.*((1-((X.*l_ch)./delta_crit_1).^2)./(1-(1-((X.*l_ch)./delta_crit_1).^2))))./(pi./2);
      Qe_times_ec_2=Qe_2ord(2,1).*X+Qe_2ord(2,2).*Y;
      Z=t_m+Qe_times_ec_2;
      mesh(X,Y,Z);
      xlabel('\epsilon^c_n_n [-]');ylabel('\gamma^c_n_m [-]'); zlabel('s(2) [-]'); legend('s_i_t_e_r(2)','s_p_r_e_c_i_s_e(2)');      
  else %Figure 4
  subplot(2,2,1); hold on
      X=g_ec_1; Y=g_ec_2; Z=g_fc_1;
      plot3(X, Y, Z,'-db','MarkerSize',7,'MarkerFaceColor','c','MarkerEdgeColor','k');
      [X,Y]=meshgrid(0:(ec_crit_1/15):ec_crit_1, -ec_crit_1:(ec_crit_1/15):0);
      Z=Sig_nn_p_scalar.*( (1+((c1.*l_ch.*X./delta_crit_1).^3)) .*exp(-c2.*l_ch.*X./delta_crit_1) - exp(-c2).*(1+(c1.^3)).*l_ch.*X./delta_crit_1);  %.*exp(-c2.*l_ch.*Eps_c_nn./delta_crit_1)    
      mesh(X,Y,Z);
      [X,Y]=meshgrid(ec_crit_1:(ec_crit_1/15):1.2*ec_crit_1, -ec_crit_1:(ec_crit_1/15):0);Z=(X-X)+(Y-Y) ; mesh(X,Y,Z);
      [X,Y]=meshgrid(0:(ec_crit_1/15):1.2*ec_crit_1, -ec_crit_1:(-ec_crit_1/15):-1.2*ec_crit_1 );Z=(X-X)+(Y-Y) ; mesh(X,Y,Z);
 
      xlabel('\epsilon^c_n_n [-]');ylabel('\gamma^c_n_m [-]'); zlabel('t_n [MPa] '); legend('t_n_,_i_t_e_r','t_n_,_p_r_e_c_i_s_e');
    subplot(2,2,2); hold on
      X=g_ec_1; Y=g_ec_2; Z=g_fc_2;
      plot3(X, Y, Z,'-db','MarkerSize',7,'MarkerFaceColor','c','MarkerEdgeColor','k'); 
      [X,Y]=meshgrid(0.00001:(ec_crit_1/30):ec_crit_1+0.00001,-0.00001:(-ec_crit_1/15):-ec_crit_1-0.00001);
      Z=-fs.*atan((-Y).*((1-((X.*l_ch)./delta_crit_1).^2)./(1-(1-((X.*l_ch)./delta_crit_1).^2))))./(pi./2);
      mesh(X,Y,Z);
           [X,Y]=meshgrid(ec_crit_1+0.00001:(ec_crit_1/15):1.2*ec_crit_1+0.00001, -ec_crit_1-0.00001:(ec_crit_1/15):-0.00001);Z=(X-X)+(Y-Y) ; mesh(X,Y,Z);
      [X,Y]=meshgrid(0.00001:(ec_crit_1/15):1.2*ec_crit_1+0.00001, -ec_crit_1-0.00001:(-ec_crit_1/15):-1.2*ec_crit_1-0.00001 );Z=(X-X)+(Y-Y) ; mesh(X,Y,Z);

      xlabel('\epsilon^c_n_n [-]');ylabel('\gamma^c_n_m [-]'); zlabel('t_m [MPa] '); legend('t_m_,_i_t_e_r','t_m_,_p_r_e_c_i_s_e');      
     subplot(2,2,3); hold on
      X=g_ec_1; Y=g_ec_2; Z=g_s_1;
      plot3(X,Y,Z,'-db','MarkerSize',7,'MarkerFaceColor','c','MarkerEdgeColor','k');  
      [X,Y]=meshgrid(0:(ec_crit_1/15):ec_crit_1, -ec_crit_1:(ec_crit_1/15):0); 
      t_n=Sig_nn_p_scalar.*( (1+((c1.*l_ch.*X./delta_crit_1).^3)) .*exp(-c2.*l_ch.*X./delta_crit_1) - exp(-c2).*(1+(c1.^3)).*l_ch.*X./delta_crit_1);  %.*exp(-c2.*l_ch.*Eps_c_nn./delta_crit_1)    
      Qe_times_ec_1=Qe_2ord(1,1).*X+Qe_2ord(1,2).*Y;
      Z=t_n+Qe_times_ec_1;
      mesh(X,Y,Z);
      [X,Y]=meshgrid(ec_crit_1+0.00001:(ec_crit_1/30):1.2*ec_crit_1+0.00001,-0.00001:(-ec_crit_1/15):-ec_crit_1-0.00001); 
      t_m= (X-X)+(Y-Y);Qe_times_ec_2=Qe_2ord(1,1).*X+Qe_2ord(1,2).*Y;
      Z=t_m+Qe_times_ec_2;
      mesh(X,Y,Z);
      [X,Y]=meshgrid(0.00001:(ec_crit_1/30):1.2*ec_crit_1+0.00001,-ec_crit_1-0.00001:(-ec_crit_1/15):-1.2*ec_crit_1-0.00001); 
      t_m= (X-X)+(Y-Y);Qe_times_ec_2=Qe_2ord(1,1).*X+Qe_2ord(1,2).*Y;
      Z=t_m+Qe_times_ec_2;
      mesh(X,Y,Z);
      xlabel('\epsilon^c_n_n [-]');ylabel('\gamma^c_n_m [-]'); zlabel('s(1) [-]'); legend('s_i_t_e_r(1)','s_p_r_e_c_i_s_e(1)');      
     subplot(2,2,4); hold on
      X=g_ec_1; Y=g_ec_2; Z=g_s_2;
      plot3(X,Y,Z,'-db','MarkerSize',7,'MarkerFaceColor','c','MarkerEdgeColor','k');  
      [X,Y]=meshgrid(0.00001:(ec_crit_1/30):ec_crit_1+0.00001,-0.00001:(-ec_crit_1/15):-ec_crit_1-0.00001); 
      t_m=-fs.*atan((-Y).*((1-((X.*l_ch)./delta_crit_1).^2)./(1-(1-((X.*l_ch)./delta_crit_1).^2))))./(pi./2);
      Qe_times_ec_2=Qe_2ord(2,1).*X+Qe_2ord(2,2).*Y;
      Z=t_m+Qe_times_ec_2;
      mesh(X,Y,Z);
      [X,Y]=meshgrid(ec_crit_1+0.00001:(ec_crit_1/30):1.2*ec_crit_1+0.00001,-0.00001:(-ec_crit_1/15):-ec_crit_1-0.00001); 
      t_m= (X-X)+(Y-Y);Qe_times_ec_2=Qe_2ord(2,1).*X+Qe_2ord(2,2).*Y;
      Z=t_m+Qe_times_ec_2;
      mesh(X,Y,Z);
      [X,Y]=meshgrid(0.00001:(ec_crit_1/30):1.2*ec_crit_1+0.00001,-ec_crit_1-0.00001:(-ec_crit_1/15):-1.2*ec_crit_1-0.00001); 
      t_m= (X-X)+(Y-Y);Qe_times_ec_2=Qe_2ord(2,1).*X+Qe_2ord(2,2).*Y;
      Z=t_m+Qe_times_ec_2;
      mesh(X,Y,Z);
      xlabel('\epsilon^c_n_n [-]');ylabel('\gamma^c_n_m [-]'); zlabel('s(2) [-]'); legend('s_i_t_e_r(2)','s_p_r_e_c_i_s_e(2)');      
  end;   
          
  figure%5
     subplot(2,2,1); hold on
      plot(g_Sig_3,g_Sig_1,'-xb')
      xlabel('\tau_x_y [MPa]'); ylabel('\sigma_x [MPa]'); 
     subplot(2,2,2); hold on
      plot(g_Sig_3,g_Sig_2,'-xb');
     xlabel('\tau_x_y [MPa]'); ylabel('\sigma_y [MPa] '); 
     subplot(2,2,3); hold on
      plot(g_delta_1,g_delta_2,'-xb'), 
      xlabel('\delta_n [mm]'); ylabel('\delta_m [mm]'); 
     subplot(2,2,4); hold on
      plot(1:length(g_Dec_G_11),g_Dec_G_11,'-b'), 
      plot(1:length(g_Dec_G_22),g_Dec_G_22,'--xb'),
      plot(1:length(g_Dec_G_33),g_Dec_G_33,':b'),
      plot(1:length(g_Dec_G_12),g_Dec_G_12,'-.b'),
      xlabel('Step_s_c'); ylabel('D_e_c(i,j) [MPa]'); legend('D_e_c(1,1)','D_e_c(2,2)','D_e_c(3,3)','D_e_c(1,2)');
