clc
clear 
close all

%% Crystallographic and thermal data %%
omega = 1.09*10^-29;
b = 2.49*10^-10;
T_m = 1726;
k = 1.381*10^-23;
R = 8.314;
d = 1*10^-6;

%%  Modulus %%
mu_0 = 7.89*10^10;
mu_td = -0.64;

T = 0:5:1726;
 
sigma_s = 0:10^5:10^10;
      
mu = (mu_0*(1+(T-300)/T_m*mu_td)).'*ones(1,length(sigma_s));

%% Plasticity limited by discrete obstacles %%  
gamma_0 = 10*6;   
delta_F = 0.5*mu_0*b^3;         % Chapter 2 Ashby's handbook
tau = 6.3*10^-3*mu_0;
sigma_s_array = (sigma_s.'*ones(1,length(T))).';
gamma_2 = gamma_0.*exp(-delta_F./(k.*T.'*ones(1,length(sigma_s))).*(1-sigma_s_array./tau));     % Equation 2.9 Ashby's handbook

%% Power law creep %%  % section 2.4 Ashby's book 
A = 3*10^6;
n = 4.6;
D_ov = 1.9*10^-4;
Q_v = 284*10^3;
acD_oc = 3.1*10^-23;
Q_c = 170*10^3;
A_2 = sqrt(3)^(n+1)*A;
D_v = (D_ov*exp(-Q_v./(R.*T))).'*ones(1,length(sigma_s));
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);               % Equation 2.20 Ashby's Handbook

gamma_4 = A_2.*D_eff.*mu.*b./(k.*T.'*ones(1,length(sigma_s))).*(sigma_s_array./mu).^n;        % Equation 2.21 Ashby's Handbook
clearvars  D_eff 
%% Diffusional flow %%  Section 2.5 Ashby's Handbook
deltaD_ob = 3.5*10^-15;
Q_b = 115*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);       %Equation 2.30 Ashby's Handbook

gamma_7 = 42*sigma_s_array*omega./(k*T.'*ones(1,length(sigma_s)).*d^2).*D_eff1;      %Equation 2.29 Ashby's Handbook
clearvars  D_eff1 
%% Calculate for net strain rate %%
gamma_net = max(gamma_2,gamma_4)+gamma_7;
y = log10(sigma_s_array./mu);
x = (T./T_m).'*ones(1,length(sigma_s));

%% Find dominant mechanism %%

% plasticity limited by discrete obstacles dominant & boundary with
% diffusion flow

 test=((gamma_2>=gamma_4).*(gamma_2>=gamma_7)).*( abs(gamma_2-gamma_7) < 0.08*gamma_7);
    Y_1 = y(test==1);
    X_1 = x(test==1);
   
% plasticity limited by discrete obstacles dominant & boundary with
% power-law creep
test=((gamma_2>=gamma_4).*(gamma_2>=gamma_7)).* ( abs(gamma_2-gamma_4) < 0.08*gamma_4);
    Y_2 = y(test==1);
    X_2 = x(test==1);

% power law creep dominant & boundary with diffusional flow
test=((gamma_4>=gamma_2).*(gamma_4>=gamma_7)).*( abs(gamma_4-gamma_7) < 0.08*gamma_7);
    Y_3 = y(test==1);
    X_3 = x(test==1);

% between H.T. power-law creep and L.T power-law creep
test=( D_v > 10./b.^2.*(sigma_s_array./mu).^2.*acD_c) .*(abs(D_v-10./b.^2*(sigma_s_array./mu).^2.*acD_c) < 0.05*10/b.^2*(sigma_s_array./mu).^2.*acD_c) .* ((gamma_4>=gamma_2).*(gamma_4>=gamma_7));
    Y_4 = y(test==1);
    X_4 = 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);
    Y_5 = y(test==1);
    X_5 = x(test==1);  
clearvars test deltaD_b D_v sigma_s_array mu acD_c gamma_4 gamma_2 gamma_7 T
%%  Deformation mechanism map %%
v = [10^-10 10^-9 10^-8 10^-7 10^-6 10^-5 10^-4 10^-3 10^-2 10^-1 1];
cs = contour(x,y,gamma_net,v);
clabel(cs);
hold on
plot(X_1,Y_1,'r','LineWidth',2)
hold on
plot(X_2,Y_2,'m','LineWidth',2)
hold on
plot(X_3,Y_3,'g','LineWidth',2)
hold on
plot(X_4,Y_4,'b--','LineWidth',2)
hold on
% axis([-5 -2 -8 1])
% plot(X_5,Y_5,'k--','LineWidth',2)
% hold on
text(0.3,-1.5,'PLASTICITY')
text(0.6,-2.5,'POWER-LAW CREEP')
text(0.5,-2.2,'(L.T.CREEP)')
text(0.8,-3,'(H.T.CREEP)')
text(0.4,-4,'DIFFUSIONAL FLOW')
text(0.4,-4.5,'(coble creep)')
title('Pure nickel of grain size 1\mum')
xlabel('Homologous Temperature T/T_m')
ylabel('Normalised Shear Stress log(\sigma/\mu)')