clc
clear 
close all

%% Crystallographic and thermal data %%
omega = 1.7*10^-29;
b = 3.0*10^-10;
T_m = 1933;
k = 1.381*10^-23;
R = 8.314;
d = 5*10^-6;

%%  Modulus %%
mu_0 = 3.4*10^10;
mu_td = -1.02;

T = 0:5:1900;
sigma_s = 0:2*10^5:2*10^9;
      
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;
tau = 1.5*10^-2*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 
n = 4.85;
Q_v = 169*10^3;
Q_c = 101*10^3;
% constants subject to change
A = 7;
D_ov = 5*10^-8;
acD_oc = 1*10^-31;

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
% D_eff = D_v; 

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 = 10^-8*6.0*10^-10;   %change this
Q_b = Q_c*10^3;   %fix this

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)/sqrt(3);
sigma_array = sqrt(3).*sigma_s_array;
y = log10(sigma_array./mu);
% x = (T./T_m).'*ones(1,length(sigma_s));
x = T.'*ones(1,length(sigma_s));

%%  Deformation mechanism map %%
v = [617 658 714 783];
cs = contour(y,log10(gamma_net),x,v);
hold on
axis([-3 -1.5 -11 -5])
clabel(cs)
grid on

text(-2.5,-7.5,'783 K')
text(-2.5,-8.5,'714 K')
text(-2.5,-9.5,'658 K')
text(-2.5,-10.5,'617 K')

%% Experimental results %%
v_1 = [4.24*10^-9 6.06*10^-10];
v_2 = [7.04*10^-10 1.58*10^-9 1.85*10^-9 3.85*10^-8];
v_3 = [2.66*10^-8 3.82*10^-9 3.68*10^-9 3.08*10^-8 1.52*10^-10 3.19*10^-9 1.81*10^-8];
v_4 = [9.29*10^-9 1.83*10^-7];

a = log10([475.83 317.2]*10^6/(mu_0*(1+(617-300)/T_m*mu_td)));
b = log10([4.24*10^-9 6.06*10^-10]);

plot(a,b,'xb','LineWidth',2)
hold on

c = log10([165.5 317.2 317.2 475.8]*10^6/(mu_0*(1+(658-300)/T_m*mu_td)));
d = log10([7.04*10^-10 1.58*10^-9 1.85*10^-9 3.85*10^-8]);

plot(c,d,'xr','LineWidth',2)
hold on

e = log10([317.2 165.5 165.5 317.2 48.3 165.5 317.2]*10^6/(mu_0*(1+(714-300)/T_m*mu_td)));
f = log10([2.66*10^-8 3.82*10^-9 3.68*10^-9 3.08*10^-8 1.52*10^-10 3.19*10^-9 1.81*10^-8]);

plot(e,f,'xg','LineWidth',2)
hold on

g = log10([48.3 165.5]*10^6/(mu_0*(1+(783-300)/T_m*mu_td)));
h = log10([9.29*10^-9 1.83*10^-7]);

plot(g,h,'xm','LineWidth',2)
hold on

title('Ti-6Al-4V of grain size 5\mum')
xlabel('Normalised Stress log(\sigma/\mu)')
ylabel('Strain Rate log(epsdot)')
