clc
clear 
close all

%% TIME %%
t = 1*10^9;
tic;

%% Crystallographic and thermal data %%
omega = 1.21*10^-29;
b = 2.58*10^-10;
T_m = 1810;
k = 1.381*10^-23;
R = 8.314;
d = 100*10^-6;   %%grain size

%%  Modulus %%
mu_0 = 8.1*10^10;
mu_td = -0.85;

T = 0:5:1810;
sigma_s = 0:10^5:10^10;

mu = (mu_0*(1+(T-300)/T_m*mu_td)).'*ones(1,length(sigma_s));
%% Elastic Deformation %%
eps_1 = ((sigma_s.'*ones(1,length(T))).')./mu;

%% Low-temperature Plasticity %%
sigma_os = (mu_0*(7.5*10^-4-2.2*10^-7*(T-273.15))).'*ones(1,length(sigma_s));
m = (0.31+6*10^-5*(T-273.15)).'*ones(1,length(sigma_s));
K_s = (mu_0*(2.5*10^-3-5.7*10^-7*(T-273.15))).'*ones(1,length(sigma_s));

sigma_s_array = (sigma_s.'*ones(1,length(T))).';

eps_2 = ((sigma_s_array-sigma_os)./K_s).^(1./m);

Aa=(sigma_s_array-sigma_os)./K_s;
Bb=Aa(:,:)>=0;
eps_2=eps_2.*Bb;

 
% clearvars sigma_os m K_s Aa Bb
%% Power-law Creep %%
A = 1*10^10;
n = 7.9;
D_ov = 3.7*10^-5;
Q_v = 280*10^3;
A_2 = sqrt(3)^(n+1)*A;
D_v = (D_ov*exp(-Q_v./(R.*T))).'*ones(1,length(sigma_s));

acD_oc = 1.0*10^-24;
Q_c = 174*10^3;
acD_c = (acD_oc*exp(-Q_c./(R.*T))).'*ones(1,length(sigma_s));
D_eff = D_v.*(1+10./b^2.*(sigma_s_array./mu).^2.*acD_c./D_v); 

gamma_3 = (A_2.*D_eff.*mu.*b./(k.*T.'*ones(1,length(sigma_s))).*(sigma_s_array./mu).^n);

eps_t = 0.087;
C_s = 46.0;

eps_3 = (eps_t.*(1-exp(-C_s.*gamma_3.*t))+gamma_3.*t);

clearvars gamma_3 D_eff
%% Diffusional Flow %%
deltaD_ob = 2*10^-13;
Q_b = 167*10^3;
deltaD_b = (deltaD_ob.*exp(-Q_b./(R.*T))).'*ones(1,length(sigma_s));
D_eff1 = D_v.*(1+pi/d.*deltaD_b./D_v);

gamma_4 = 42*sigma_s_array*omega./(k*T.'*ones(1,length(sigma_s)).*d^2).*D_eff1;

eps_4 = sigma_s_array./mu.*(1-exp(-(gamma_4.*t)./(sigma_s_array./mu)))+gamma_4.*t;

clearvars gamma_4 D_eff1
%% Total strain %%
eps = eps_1+eps_2+eps_3+eps_4;
y = log10(sigma_s_array./mu);
x = (T./T_m).'*ones(1,length(sigma_s));


% elasticity deformation dominant & boundary with plasticity
test=((eps_1>=eps_2).* (eps_1>=eps_3) .* (eps_1>=eps_4)).*( abs(eps_1-eps_2) < 0.08.*eps_2);

    Y_1 = y(test==1);
    X_1 = x(test==1);
    

% plasticity dominent & boundary with power-law creep
test= ((eps_2>=eps_1).* (eps_2>=eps_3) .* (eps_2>=eps_4)).*  (abs(eps_2-eps_3) < 0.08.*eps_3);

    Y_2 = y(test==1);
    X_2 = x(test==1);
    
    
% diffusional flow dominent & boundary with plasticity (use this when time
% is large)
test= ((eps_4>=eps_1).* (eps_4>=eps_3) .* (eps_4>=eps_2)).* ( abs(eps_4-eps_2) < 0.08.*eps_2);

    Y_3 = y(test==1);
    X_3 = x(test==1);
    
% % % power-law creep dominent & boundary with elastic deformation (use this
% % when time is small)
% test= ((eps_3>=eps_1).* (eps_3>=eps_2) .* (eps_3>=eps_4)).* ( abs(eps_3-eps_1) < 0.01.*eps_1);
% 
%     Y_3 = y(test==1);
%     X_3 = x(test==1);


 
% power-law creep dominent & boundary with diffusional flow
test=((eps_3>=eps_1).* (eps_3>=eps_4) .* (eps_3>=eps_2)).* ( abs(eps_3-eps_4) < 0.08.*eps_4);

    Y_4 = y(test==1);
    X_4 = x(test==1);
[c,pos] = sort(X_4);
Y_4 = Y_4(pos);
X_4 = c;

% Diffusional flow dominent & boundary with elastic deformation
test=((eps_4>=eps_1).* (eps_4>=eps_3) .* (eps_4>=eps_2)).* (abs(eps_4-eps_1) < 0.15*eps_1);

    Y_5 = y(test==1);
    X_5 = x(test==1);

% between Nabarro-Herring creep and Coble creep
test= (D_v > pi/d*deltaD_b) .* (abs(D_v-pi/d*deltaD_b) < 0.05*pi/d*deltaD_b) .*((eps_4>eps_1).* (eps_4>eps_3) .* (eps_4>eps_2));

    Y_6 = y(test==1);
    X_6 = x(test==1);
    
   
% clearvars Aa Bb test D_v  deltaD_b d
%%  Transient Map %%
v = [10^-5 10^-4 10^-3 10^-2 10^-1 1 10];
cs = contour(x,y,eps,v);
clabel(cs);
hold on
plot(X_1,Y_1,'r','LineWidth',2)
hold on
plot(X_2,Y_2,'b','LineWidth',2)
hold on
plot(X_3,Y_3,'y','LineWidth',2)
hold on
plot(X_4,Y_4,'g','LineWidth',2)
hold on
plot(X_5,Y_5,'k','LineWidth',2)
hold on
plot(X_6,Y_6,'c--','LineWidth',2)
hold on
% axis([0 1 -6 -1])

% % for 3600s (1 hour)
% text(0.05,-3.7,'ELASTIC DEFORMATION')
% text(0.2,-2,'PLASTICITY')
% text(0.6,-2,'POWER-LAW CREEP')
% text(0.72,-3.7,'DIFFUSIONAL FLOW')
% text(0.75,-4,'(N-H)')
% text(0.87,-4,'(Coble)')
% title('steel 316 of grain size 100*10^-6, for a time of 3600 s')

% % for 1*10^8s (3 years)
% text(0.05,-3.7,'ELASTIC DEFORMATION')
% text(0.2,-2,'PLASTICITY')
% text(0.6,-2,'POWER-LAW CREEP')
% text(0.7,-3.7,'DIFFUSIONAL FLOW')
% text(0.75,-4,'(N-H)')
% text(0.87,-4,'(Coble)')
% title('steel 316 of grain size 100*10^-6, for a time of 3 years')

% for 1*10^9s (30 years)
text(0.03,-3.7,'ELASTIC DEFORMATION')
text(0.2,-2,'PLASTICITY')
text(0.6,-2,'POWER-LAW CREEP')
text(0.7,-3.7,'DIFFUSIONAL FLOW')
text(0.75,-4,'(N-H)')
text(0.87,-4,'(Coble)')
title('steel 316 of grain size 100*10^-6, for a time of 30 years')

% % for 10000s & 75 micro meter (3 hour)
% text(0.05,-3.7,'ELASTIC DEFORMATION')
% text(0.2,-2,'PLASTICITY')
% text(0.6,-2,'POWER-LAW CREEP')
% text(0.7,-3.7,'DIFFUSIONAL FLOW')
% text(0.75,-4,'(N-H)')
% text(0.9,-4,'(Coble)')
% title('steel 316 of grain size 75*10^-6, for a time of 1*10^4 s')

xlabel('Homologous Temperature T/T_m')
ylabel('Normalised Shear Stress log(sigma/mu)')
grid on
toc
