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;   %fi4x 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.'*ones(1,length(sigma_s));

%% Deformation Mechanism Map with Temperature Contours %%
v = [773 923 1073];
cs = contour(y,log10(gamma_net),x,v);
hold on
axis([-4.5 -2 -10 -2])
clabel(cs)
grid on

text(-2.7,-2.5,'1073 K')
text(-2.7,-4,'923 K')
text(-2.7,-5.5,'773 K')

%% Experimental Data %%
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;

a = log10([62.1 62.0 62.0 62.2 72.9 72.8 73.0 72.8 84.0 84.0 84.1 83.9]*10^6/(mu_0+mu_1*773));
b = log10([7.65*10^-5 2.71*10^-4 2.43*10^-4 1.09*10^-3 1.42*10^-3 1.93*10^-3 1.62*10^-3 7.06*10^-4 5.46*10^-3 8.66*10^-3 9.32*10^-3 7.43*10^-3]/3600);

scatter(a,b,'x','b','LineWidth',2)
hold on

c = log10([17.9 17.7 18.0 17.8 27.7 27.6 27.7 27.7 37.3 37.4 37.6 37.4]*10^6/(mu_0+mu_1*923));
d = log10([3.72*10^-5 1.58*10^-4 2.16*10^-4 1.53*10^-4 6.10*10^-4 8.43*10^-4 1.54*10^-3 7.78*10^-4 4.43*10^-2 1.02*10^-2 3.86*10^-3 2.59*10^-2]/3600);

scatter(c,d,'x','r','LineWidth',2)
hold on

e = log10([6.2 6.2 5.7 6.2 10.3 10.3 10.3 10.4 14.7 14.6 14.6 14.6]*10^6/(mu_0+mu_1*1073));
f = log10([2.61*10^-4 7.00*10^-4 1.68*10^-3 4.42*10^-4 3.27*10^-2 2.36*10^-2 2.39*10^-2 1.79*10^-2 1.33*10^-1 2.98*10^-2 1.17*10^-1 1.94*10^-2]/3600);

scatter(e,f,'x','g','LineWidth',2)
hold on

plot([log10(62.1*10^6/(mu_0+mu_1*773)),log10(72.9*10^6/(mu_0+mu_1*773)),log10(84.0*10^6/(mu_0+mu_1*773))],[mean(log10(v_1)),mean(log10(v_2)),mean(log10(v_3))],'LineWidth',2)
hold on

plot([log10(17.8*10^6/(mu_0+mu_1*923)),log10(27.6*10^6/(mu_0+mu_1*923)),log10(37.4*10^6/(mu_0+mu_1*923))],[mean(log10(w_1)),mean(log10(w_2)),mean(log10(w_3))],'LineWidth',2)
hold on

plot([log10(6.2*10^6/(mu_0+mu_1*1073)),log10(10.4*10^6/(mu_0+mu_1*1073)),log10(14.6*10^6/(mu_0+mu_1*1073))],[mean(log10(u_1)),mean(log10(u_2)),mean(log10(u_3))],'LineWidth',2)
hold on

title('NARloy Z of grain size 100\mum')
xlabel('Normalised Stress log(\sigma/\mu)')
ylabel('Strain Rate log(epsdot)')

%% 95% Confidence Interval %%
v_1 = log10(v_1);
v_2 = log10(v_2);
v_3 = log10(v_3);
w_1 = log10(w_1);
w_2 = log10(w_2);
w_3 = log10(w_3);
u_1 = log10(u_1);
u_2 = log10(u_2);
u_3 = log10(u_3);

[H_1,P_1,C_1] = ttest(v_1,mean(v_1));
[H_2,P_2,C_2] = ttest(v_2,mean(v_2));
[H_3,P_3,C_3] = ttest(v_3,mean(v_3));
[H_4,P_4,C_4] = ttest(w_1,mean(w_1));
[H_5,P_5,C_5] = ttest(w_2,mean(w_2));
[H_6,P_6,C_6] = ttest(w_3,mean(w_3));
[H_7,P_7,C_7] = ttest(u_1,mean(u_1));
[H_8,P_8,C_8] = ttest(u_2,mean(u_2));
[H_9,P_9,C_9] = ttest(u_3,mean(u_3));
L = [mean(v_1)-C_1(1) mean(v_2)-C_2(1) mean(v_3)-C_3(1) mean(w_1)-C_4(1) mean(w_2)-C_5(1) mean(w_3)-C_6(1) mean(u_1)-C_7(1) mean(u_2)-C_8(1) mean(u_3)-C_9(1)];
U = [C_1(2)-mean(v_1) C_2(2)-mean(v_2) C_3(2)-mean(v_3) C_4(2)-mean(w_1) C_5(2)-mean(w_2) C_6(2)-mean(w_3) C_7(2)-mean(u_1) C_8(2)-mean(u_2) C_9(2)-mean(u_3)];

x_1 = log10([62.1 72.9 84.0]*10^6/(mu_0+mu_1*773));
x_2 = log10([17.8 27.6 37.4]*10^6/(mu_0+mu_1*923));
x_3 = log10([6.2 10.4 14.6]*10^6/(mu_0+mu_1*1073));

x = [x_1 x_2 x_3];
y_data = [mean(v_1) mean(v_2) mean(v_3) mean(w_1) mean(w_2) mean(w_3) mean(u_1) mean(u_2) mean(u_3)];
errorbar(x,y_data,L,U,'.')
hold on

