clc
clear 
close all

%% Crystallographic and thermal data %%
omega = 1.18*10^-29;
b = 2.56*10^-10;
T_m = 1350;
k = 1.381*10^-23;
R = 8.314;
d = 100*10^-6;

%%  Modulus %%
mu_0 = 54000*10^6;
mu_1 = -17*10^6;

T = 0:5:1350;
 
sigma_s = 0:10^5:10^9;
      
mu = (mu_0+mu_1*T).'*ones(1,length(sigma_s));

%% Plasticity limited by discrete obstacles %%  
gamma_0 = 10*6;   
delta_F = 0.5*mu_0*b^3;
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 
n = 3.8;
Q_v = 193*10^3;
Q_c = 122*10^3;
% constants subject to change
A = 11;
D_ov = 2.0*10^-6;
acD_oc = 1.0*10^-24;

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 = 5.0*10^-14;   %change this
Q_b = 104*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));

%% Experimental Results %%
a = log10([62.1 72.9 84.0]*10^6/(mu_0+mu_1*773));
b = [773 773 773]/T_m;

c = log10([17.8 27.6 37.4]*10^6/(mu_0+mu_1*923));
d = [923 923 923]/T_m;

e = log10([6.2 10.4 14.6]*10^6/(mu_0+mu_1*1073));
f = [1073 1073 1073]/T_m;

v_1 = [7.65*10^-5 2.71*10^-4 2.43*10^-4 1.09*10^-3]/3600;
v_2 = [1.42*10^-3 1.93*10^-3 1.62*10^-3 7.06*10^-4]/3600;
v_3 = [5.46*10^-3 8.66*10^-3 9.32*10^-3 7.43*10^-3]/3600;

w_1 = [3.72*10^-5 1.58*10^-4 2.16*10^-4 1.53*10^-4]/3600;
w_2 = [6.10*10^-4 8.43*10^-4 1.54*10^-3 7.78*10^-4]/3600;
w_3 = [4.43*10^-2 1.02*10^-2 3.86*10^-3 2.59*10^-2]/3600;

u_1 = [2.61*10^-4 7.00*10^-4 1.68*10^-3 4.42*10^-4]/3600;
u_2 = [3.27*10^-2 2.36*10^-2 2.39*10^-2 1.79*10^-2]/3600;
u_3 = [1.33*10^-1 2.98*10^-2 1.17*10^-1 1.94*10^-2]/3600;
%% Plot Deformation Mechanism Maps and Compare with Experimental Data
figure
contour(x,y,gamma_net,v_1,'b','LineWidth',2);
hold on
plot(b(1),a(1),'ob','LineWidth',2)
hold on
% xlabel('Homologous Temperature T/T_m')
ylabel('Normalized Stress log(\sigma/\mu)')
text(0.5,-1.5,'Loading Condition 773K, 62.1MPa')

figure
contour(x,y,gamma_net,v_2,'r','LineWidth',2);
hold on
plot(b(2),a(2),'or','LineWidth',2)
hold on
% xlabel('Homologous Temperature T/T_m')
% ylabel('Normalised Stress log(\sigma/\mu)')
text(0.5,-1.5,'Loading Condition 773K, 72.9MPa')

figure
contour(x,y,gamma_net,v_3,'g','LineWidth',2);
hold on
plot(b(3),a(3),'og','LineWidth',2)
hold on
% xlabel('Homologous Temperature T/T_m')
% ylabel('Normalised Stress log(\sigma/\mu)')
text(0.5,-1.4,'Loading Condition 773K, 84.0MPa')

figure
contour(x,y,gamma_net,w_1,'b','LineWidth',2);
hold on
plot(d(1),c(1),'ob','LineWidth',2)
hold on
% xlabel('Homologous Temperature T/T_m')
ylabel('Normalized Stress log(\sigma/\mu)')
text(0.5,-1.5,'Loading Condition 923K, 17.8MPa')

figure
contour(x,y,gamma_net,w_2,'r','LineWidth',2);
hold on
plot(d(2),c(2),'or','LineWidth',2)
hold on
% xlabel('Homologous Temperature T/T_m')
% ylabel('Normalised Stress log(\sigma/\mu)')
text(0.5,-1.5,'Loading Condition 923K, 27.6MPa')

figure
contour(x,y,gamma_net,w_3,'g','LineWidth',2);
hold on
plot(d(3),c(3),'og','LineWidth',2)
hold on
% xlabel('Homologous Temperature T/T_m')
% ylabel('Normalised Stress log(\sigma/\mu)')
text(0.5,-1.45,'Loading Condition 923K, 37.4MPa')

figure
contour(x,y,gamma_net,u_1,'b','LineWidth',2);
hold on
plot(f(1),e(1),'ob','LineWidth',2)
hold on
xlabel('Homologous Temperature T/T_m')
ylabel('Normalized Stress log(\sigma/\mu)')
text(0.5,-1.5,'Loading Condition 1073K, 6.2MPa')

figure
contour(x,y,gamma_net,u_2,'r','LineWidth',2);
hold on
plot(f(2),e(2),'or','LineWidth',2)
hold on
xlabel('Homologous Temperature T/T_m')
% ylabel('Normalised Stress log(\sigma/\mu)')
text(0.5,-1.4,'Loading Condition 1073K, 10.4MPa')

figure
contour(x,y,gamma_net,u_3,'g','LineWidth',2);
hold on
plot(f(3),e(3),'og','LineWidth',2)
hold on
% title('NARloy Z of grain size 100*10^-6')
xlabel('Homologous Temperature T/T_m')
% ylabel('Normalised Stress log(\sigma/\mu)')
text(0.5,-1.4,'Loading Condition 1073K, 14.6MPa')

