clc
clear 
close all

%% TIME %%
t = 30;
tic;

%% 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^4:10^9;

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

%% Low-temperature Plasticity %%
% yield stress
test=(T<400);
sigma_yield1 = (170-0.1*T(test==1))*10^6;

test=(T>=400 & T<600);
sigma_yield2 = (130-0.075*(T(test==1)-400))*10^6;

test=(T>=600 & T<800);
sigma_yield3 = (115-0.175*(T(test==1)-600))*10^6;

test=(T>=800 & T<1000);
sigma_yield4 = (80-0.225*(T(test==1)-800))*10^6;

test=(T>=1000 & T<1150);
sigma_yield5 = (35-0.167*(T(test==1)-1000))*10^6;

test=(T>=1150);
sigma_yield6 = (10-0*T(test==1))*10^6;

sigma_yield = [sigma_yield1 sigma_yield2 sigma_yield3 sigma_yield4 sigma_yield5 sigma_yield6];

% work-hardening parameters
sigma_os = (sigma_yield/sqrt(3)).'*ones(1,length(sigma_s));
m = 0.56.'*ones(length(T),length(sigma_s));
K_s = ((530*10^6)/(3^((0.56+1)/2))).'*ones(length(T),length(sigma_s));

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

% plastic deformation
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 = 11;
n = 3.8;
Q_v = 193*10^3;
D_ov = 2.0*10^-6;
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 = 122*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 = 8.484*10^-3*sqrt(3);
C_s = 277/sqrt(3);

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

clearvars gamma_3 D_eff
%% Diffusional Flow %%
deltaD_ob = 5.0*10^-14;   
Q_b = 104*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)/sqrt(3); %total transient strain
eps1 = (eps_2+eps_3+eps_4)/sqrt(3);   %total plastic strain

sigma_array = sigma_s_array*sqrt(3);   %uniaxial stress
y = log10(sigma_array./mu);    %normalized uniaxial stress

y_1 = sigma_array;   %uniaxial stress

x = (T./T_m).'*ones(1,length(sigma_s));   %homologous temperature
x1 = T.'*ones(1,length(sigma_s));   %temperature array

%% Calculate Boundaries between Mechanisms %%
% 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.001.*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.00005.*eps_3);

    Y_2 = y(test==1);
    X_2 = x(test==1);
% sort by stress
[c,pos] = sort(Y_2);
X_2 = X_2(pos);
Y_2 = c;   
X_2 = [max(X_1);X_2];
Y_2 = [min(Y_1);Y_2];
   
% % 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.1.*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.001.*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);
% sort by temperature
[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);
    
% temperature abd stress generate 0.2% plastic strain
test = (eps1 > 0.00199 & eps1<= 0.002);

X_7 = x(test==1);   %homologous temperature generate 0.2% plastic strain
X_71 = x1(test==1);   %temperature generate 0.2% plastic strain

Y_72 = y_1(test==1)/1.25;   %normalized stress with safety factor 1.25

% sort by temperature
[c,pos] = sort(X_7);
Y_72 = Y_72(pos);
X_7 = c;

mu1 = (mu_0+mu_1.*X_71);   %shear modulus corresponding with temperature generate 0.2% plastic strain
   
clearvars Aa Bb test D_v  deltaD_b d
%%  Transient Map %%
figure
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,'g','LineWidth',2)
hold on
plot(X_4,Y_4,'m','LineWidth',2)
hold on
plot(X_5,Y_5,'k','LineWidth',2)
hold on
plot(X_6,Y_6,'c--','LineWidth',2)
hold on
plot(X_7,log10(Y_72./mu1),'m','LineWidth',2)
hold on

text(0.03,-3.7,'ELASTIC DEFORMATION')
text(0.2,-2,'PLASTICITY')
text(0.7,-2.2,'POWER-LAW CREEP')
title('NARloy-Z of grain size 100 \mum, for a time of 30 s')

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

% plot design stress
figure
plot(((X_7*T_m-273.15)*9/5+32),Y_72/10^6,'b','LineWidth',2)
hold on
title('NARloy-Z of grain size 100 \mum, for a time of 30 s')
xlabel('Temperature T (F)')
ylabel('Design Stress \sigma_D (MPa)')
grid on

X30 = ((X_7*T_m-273.15)*9/5+32);
Y30 = Y_72/10^6;
save('30','X30','Y30')