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;

for i = 0:8;
    i = i+1;
    for j = 0:3;
        j = j+1;
        
% information from experimental results (temperature, stress and time to 1%
% strain)
Temp = [773 773 773 773;773 773 773 773;773 773 773 773;923 923 923 923;923 923 923 923;923 923 923 923;1073 1073 1073 1073;1073 1073 1073 1073;1073 1073 1073 1073];
t1 = [78.63 32.68 33.17 7.67;6.46 3.93 5.24 12.43;1.54 1.03 0.62 1.29;96.19 29.13 2.39 38.50;8.59 7.96 1.87 9.56;0.34 0.56 2.38 0.42;18.24 3.35 4.10 6.54;0.31 0.44 0.38 0.70;0.07 0.27 0.08 0.44]*3600;
sigma_s1 = [62.1 62.0 62.0 62.2;72.9 72.8 73.0 72.8;84.0 84.0 84.1 83.9;17.9 17.7 18.0 17.8;27.7 27.6 27.7 27.7;37.3 37.4 37.6 37.4;6.2 6.2 5.7 6.2;10.3 10.3 10.3 10.4;14.7 14.6 14.6 14.6]/sqrt(3)*10^6;
        
T = Temp(i,j);
sigma_s = sigma_s1(i,j);
t = t1(i,j);
mu = mu_0+mu_1*T;

%% Elastic Deformation %%
eps_1(i,j) = sigma_s/mu;

%% Low-temperature Plasticity %%
if T<400;
    sigma_yield = (170-0.1*T)*10^6;
end
if T>=400 && T<600;
    sigma_yield = (130-0.075*(T-400))*10^6;
end
if T>=600 && T<800;
    sigma_yield = (115-0.175*(T-600))*10^6;
end
if T>=800 && T<1000;
    sigma_yield = (80-0.225*(T-800))*10^6;
end
if T>=1000 && T<1150;
    sigma_yield = (35-0.167*(T-1000))*10^6;
end
if T>=1150;
    sigma_yield = 10*10^6;
end
sigma_os = sigma_yield/sqrt(3);
m = 0.56;
K_s = (530*10^6)/(3^((m+1)/2));

if (sigma_s-sigma_os)/K_s >= 0;
    eps_2(i,j) = ((sigma_s-sigma_os)/K_s)^(1/m);
end

if (sigma_s-sigma_os)/K_s < 0;
    eps_2(i,j) = 0;
end

%% Power-law Creep %%
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));
acD_c = acD_oc*exp(-Q_c/(R*T));
D_eff = D_v*(1+10/b^2*(sigma_s/mu)^2*acD_c/D_v);

gamma_3 = A_2*D_eff*mu*b/(k*T)*(sigma_s/mu)^n;

eps_t = 8.484*10^-3*sqrt(3);
C_s = 277/sqrt(3);

eps_3(i,j) = eps_t*(1-exp(-C_s*gamma_3*t))+gamma_3*t;

%% Diffusional Flow %%
deltaD_ob = 5.0*10^-14;   %change this
Q_b = 104*10^3;   %fix this 
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*omega/(k*T*d^2)*D_eff1;

eps_4(i,j) = sigma_s/mu*(1-exp(-(gamma_4*t)/(sigma_s/mu)))+gamma_4*t;

%% Total strain %%
eps(i,j) = ((eps_1(i,j)+eps_2(i,j)+eps_3(i,j)+eps_4(i,j))/sqrt(3))*100
ave = mean(eps,2)

    end
        
    j = 0;
    
end

%% 95% Confidence Interval %%
figure
[H_1,P_1,C_1] = ttest(eps(1,:),mean(eps(1,:)));
[H_2,P_2,C_2] = ttest(eps(2,:),mean(eps(2,:)));
[H_3,P_3,C_3] = ttest(eps(3,:),mean(eps(3,:)));
[H_4,P_4,C_4] = ttest(eps(4,:),mean(eps(4,:)));
[H_5,P_5,C_5] = ttest(eps(5,:),mean(eps(5,:)));
[H_6,P_6,C_6] = ttest(eps(6,:),mean(eps(6,:)));
[H_7,P_7,C_7] = ttest(eps(7,:),mean(eps(7,:)));
[H_8,P_8,C_8] = ttest(eps(8,:),mean(eps(8,:)));
[H_9,P_9,C_9] = ttest(eps(9,:),mean(eps(9,:)));
L = [C_1(1) C_2(1) C_3(1) C_4(1) C_5(1) C_6(1) C_7(1) C_8(1) C_9(1)];
U = [C_1(2) C_2(2) C_3(2) C_4(2) C_5(2) C_6(2) C_7(2) C_8(2) C_9(2)];
L = [mean(eps(1,:))-C_1(1) mean(eps(2,:))-C_2(1) mean(eps(3,:))-C_3(1) mean(eps(4,:))-C_4(1) mean(eps(5,:))-C_5(1) mean(eps(6,:))-C_6(1) mean(eps(7,:))-C_7(1) mean(eps(8,:))-C_8(1) mean(eps(9,:))-C_9(1)];
U = [C_1(2)-mean(eps(1,:)) C_2(2)-mean(eps(2,:)) C_3(2)-mean(eps(3,:)) C_4(2)-mean(eps(4,:)) C_5(2)-mean(eps(5,:)) C_6(2)-mean(eps(6,:)) C_7(2)-mean(eps(7,:)) C_8(2)-mean(eps(8,:)) C_9(2)-mean(eps(9,:))];

x = 1:9;
y = [1 1 1 1 1 1 1 1 1];
scatter(x,y)
hold on
y_data = [mean(eps(1,:)) mean(eps(2,:)) mean(eps(3,:)) mean(eps(4,:)) mean(eps(5,:)) mean(eps(6,:)) mean(eps(7,:)) mean(eps(8,:)) mean(eps(9,:))];
errorbar(x,y_data,L,U,'.')
hold on
ylabel('strain %')
set(gca,'XTick',[])
legend('1% strain (from experiment)','average strain calcualted from transient map')

