clc
clear 
close all

%% TIME %%
t = 1800;
tic;

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

%% Loading conditions from experimental results %%
T = [617 617 658 658 658 714 714 714 714 658 714 714 714 783 783];
sigma_s = [475.83 317.2 165.5 317.2 317.2 317.2 165.5 165.5 317.2 475.8 48.3 165.5 317.2 48.3 165.5]/sqrt(3)*10^6;

mu = mu_0*(1+(T-300)/T_m*mu_td);
%% Elastic Deformation %%
eps_1 = sigma_s./mu;

%% Low-temperature Plasticity %%
test=(T<366.4833);
sigma_yield1 = (1420.1-1.3163*T(test==1))*10^6;

test=(T>=366.4833 & T<477.5944);
sigma_yield2 = (1256.1-0.8688*T(test==1))*10^6;

test=(T>=477.5944 & T<588.7056);
sigma_yield3 = (1196.8-0.7447*T(test==1))*10^6;

test=(T>=588.7056 & T<699.8167);
sigma_yield4 = (1269.9-0.8688*T(test==1))*10^6;

test=(T>=699.8167 & T<810.9278);
sigma_yield5 = (1530.5-1.2411*T(test==1))*10^6;

test=(T>=810.9278 & T<1088.706);
sigma_yield6 = (1878.7-1.6705*T(test==1))*10^6;

test=(T>=1088.706 & T<1123.15);
sigma_yield7 = (495.86-0.4004*T(test==1))*10^6;

test=(T>=1123.15 & T<1218.15);
sigma_yield8 = (331.51-0.254*T(test==1))*10^6;

test=(T>=1218.15 & T<1366.483);
sigma_yield9 = (167.59-0.1195*T(test==1))*10^6;

test=(T>=1366.483 & T<1588.706);
sigma_yield10 = (20.879-0.0121*T(test==1))*10^6;

test=(T>=1588.706);
sigma_yield11 = (1.6548-0*T(test==1))*10^6;

sigma_yield = [sigma_yield1 sigma_yield2 sigma_yield3 sigma_yield4 sigma_yield5 sigma_yield6 sigma_yield7 sigma_yield8 sigma_yield9 sigma_yield10 sigma_yield11];

sigma_os = sigma_yield/sqrt(3);
m = ((0.5).*ones(length(T),1))';
K_s = ((600*10^6/(3^((0.5+1)/2))).*ones(length(T),1))';

sigma_s_array = sigma_s;

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 = 7;
n = 4.85;
D_ov = 5*10^-8;
Q_v = 169*10^3;
A_2 = sqrt(3)^(n+1)*A;
D_v = (D_ov*exp(-Q_v./(R.*T)));

acD_oc = 1*10^-31;
Q_c = 101*10^3;
acD_c = (acD_oc*exp(-Q_c./(R.*T)));
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).*(sigma_s_array./mu).^n);

eps_t = 0.04*sqrt(3);
C_s = 44.35/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 = 10^-8*6.0*10^-10;   
Q_b = 101*10^3;   
deltaD_b = (deltaD_ob.*exp(-Q_b./(R.*T)));
D_eff1 = D_v.*(1+pi/d.*deltaD_b./D_v);

gamma_4 = 42*sigma_s_array*omega./(k*T.*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);
   
clearvars Aa Bb test D_v  deltaD_b d
%%  Transient Map %%
x = 1:15;
% total strain from experimental data
y = [0.571 0.394 0.231 0.4135 0.403 0.4695 0.227 0.2275 0.5115 0.64 0.073 0.242 0.496 0.094 0.344];
semilogy(x,y,'*b','LineWidth',2)
hold on
y_data = eps*100;
semilogy(x,y_data,'*r','LineWidth',2)
hold on
ylabel('strain %')
set(gca,'XTick',[])
title('compare plastic strain at 30mins')
legend('experimental data','results from transient map')