%This model uses the micromechanical model and texture evolution for
%olivine aggregates described by Hansen et al., (2016, JGR) for modelling
%shear stresses and related deformation paths in the asthenosphere (Kiraly
%et al.). The example here is a simple shear case where the shear stress
%is constant (like model1 and model2 in Kiraly et al., ()).

%Created by A. Kiraly and L.N. Hansen

clearvars

%%%Variables
    yrtosec=60*60*24*365;

    Fnspull=4.1e12; %N/m  net slab pull which transfers from slab bouyancy to pulling the plate, lower bound after Schellart 2004 GRL
    Platelength=6000000; %m
    asthenosphere=200000; %m
    T = 1600;%Temperature (K) 
    d_o = 1000;%Grain size in microns
    %isotropic portion of viscosity
    n = 3.5; %stress exponent
    A = 1.1e5*exp(-530000/8.314/T); 
    %Slip system activities
    tau = [.3 .27 1.29]; %from Hansen et al., 2016a
    A_ss = (1./tau).^n;
    %Input for texture evolution model
    f = [1.35  1.45 .135 .068]; %from Hansen et al., 2016a
    %Number of grains
    ngrain=1000;
    %model time parameters
    max_strain=12; %model will stop when the accumulated strain reaches 21
    max_timestep=500;
    %strain rate amplitude input for the micromechanical model as described
    %in eq 6
    eps0=5e-15;
    E=[eps0 -eps0/2 -eps0/2 0 0 0;-eps0/2 eps0 -eps0/2 0 0 0;-eps0/2 -eps0/2 eps0 0 0 0;0 0 0 sqrt(2)*eps0 0 0;0 0 0 0 sqrt(2)*eps0 0;0 0 0 0 0 sqrt(2)*eps0];
    
    %Definition of the deformation tensor (D=DT*2*epsdot) as described in
    %eq(9)
    DT=[0 1 1; 
        0 0 1; 
        0 0 0]; %

%%%Shear stress
    given_stress=Fnspull/Platelength*1e-6;%From pull forca and plate area in MPa
    stress_g=[0 0 0 0 given_stress*sqrt(2) 0]; %In Kelvin notation #Representing F1 on Fig1 and Fig4
    snm=[0 0 given_stress;0 0 0;given_stress 0 0];
    Snew=eig(snm.*10^6);%Eigen values if the imposed stress
    SigmaEnew=sqrt(3*0.5*(Snew(1)^2+Snew(2)^2+Snew(3)^2));
    stress_n=(stress_g*10^6*SigmaEnew^(n-1)).*3^(-(n-1)/2);%stress component in eq(5)

%Initial values (strain and time start from 0)
    strain_store(1)=0; 
    t(1)=0*yrtosec;
    ts=1*1e5*yrtosec; %initial timestep. ts evolves according to the strain rate, so the strain step is around stst
    stst=0.075;



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Model calculation
%%
%generate aggregate with X crystals or load it previously
temp = randq(ngrain); %This requires MTEX
ori_quat = [temp(:).a temp(:).b temp(:).c temp(:).d];



for j = 1:max_timestep
    
    t(j+1)=t(j)+ts;
        
%convert orientations in quarternions to Euler angles requires MTEX
        euler=Euler(orientation('Quaternion',ori_quat(:,1,j),ori_quat(:,2,j),ori_quat(:,3,j),ori_quat(:,4,j),crystalSymmetry('mmm')),'Bunge')*180./pi;

%Micro-mechanical model going through E
        stress_coln=zeros(6,6);
        nstress_col=zeros(6,6);
for col=2:6
            nrate_ij=[E(1,col) E(6,col)/sqrt(2) E(5,col)/sqrt(2);E(6,col)/sqrt(2) E(2,col) E(4,col)/sqrt(2);E(5,col)/sqrt(2) E(4,col)/sqrt(2) E(3,col)];
            nstress_ij= Rate_Aggregate_Plasticity_CRATE(nrate_ij,euler,A_ss,A,d_o);
            nstress_ij(1,1)=-nstress_ij(2,2)-nstress_ij(3,3); %get rid of tiny numerical errors
            nstress_ij_Pa=nstress_ij.*10^6;
            nS=eig(nstress_ij_Pa); %eigen values of the deviatoric stress tensor
            nSigmaE=sqrt(3*0.5*(nS(1)^2+nS(2)^2+nS(3)^2)); %square root of the second invariant of the deviatoric stress
            nstress_col(:,col)=[nstress_ij_Pa(1,1) nstress_ij_Pa(2,2) nstress_ij_Pa(3,3) nstress_ij_Pa(2,3) nstress_ij_Pa(1,3) nstress_ij_Pa(1,2)]';
            stress_coln(:,col)=nstress_col(:,col)*(nSigmaE^(n-1)).*3^(-(n-1)/2);
end
  
        nA=E(2:6,2:6)/stress_coln(2:6,2:6); %A'=E'*S'^-1, eq(8)
        A_save(:,:,j)=nA;
        rate_new(2:6)=nA*stress_n(2:6)';rate_new(1)=-rate_new(2)-rate_new(3);
        rate_newsave(:,j)=rate_new;
        rate_ij=[rate_new(1) rate_new(6)/sqrt(2) rate_new(5)/sqrt(2); rate_new(6)/sqrt(2) rate_new(2) rate_new(4)/sqrt(2); rate_new(5)/sqrt(2) rate_new(4)/sqrt(2) rate_new(3)];
        Anisviscosity(j,:)=given_stress*1e6./rate_new(4:6)'; %eta2313, eta1313 eta 1213


    if j == max_timestep   
        break
    end
    
% % Texture calculation for time j+1
    rate_def=rate_ij.*DT;
    D = rate_def(:,:).*2./norm(rate_def(:,:).*2); %normalized velocity gradient tensor
    temp = linspace(t(j),t(j+1),3).*norm(rate_def(:,:)).*2;
    t_dir = temp - temp(1); %normalized time for director to run, effectively the shear strain for this step, broken into three steps
    ts=stst/(norm(rate_newsave(4:6,j)).*2);%update timestep based on the strain step and the last ts-s strain rate
    strain_store = [strain_store t_dir(end)+strain_store(end)];
    ori_quat_temp = director_quat(D,t_dir,ori_quat(:,:,j),f); %calculate next step of orientations
    ori_quat(:,:,j+1) = ori_quat_temp(:,:,end); %store orientations

% % Terminate process when reached <max_strain>
    if strain_store(end)>= max_strain
        max_timestep=j+1;
    end

end
