clc
clear
close all

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

eps_dot_1 = [v_1 v_2 v_3];
eps_dot_2 = [w_1 w_2 w_3];
eps_dot_3 = [u_1 u_2 u_3];

T = [773 923 1073];

Q = 450000;
Tm = 1350;
k = 8.314;

v = exp(-Q./(k.*T));

Z_1 = eps_dot_1./v(1);
Z_2 = eps_dot_2./v(2);
Z_3 = eps_dot_3./v(3);

Z = [Z_1 Z_2 Z_3];

sigma = [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];
S_magnitude = sqrt(1/3).*sigma*10^6 ;

%% CHANGE THESE PARAMETERS %%
A = 4*10^16;
n = 3.7;
C = 4.5*10^6;

%% PLOT FIGURES %%
i = 0;
for sigma_1 = 2*10^6:10^6:55*10^6
    i = i+1;
    Z_1 = A*(sinh(sqrt(1/3)*sigma_1/C))^n;
    sigma1(i) = sigma_1;
    Z1(i) = Z_1;
end

figure
loglog(S_magnitude,Z,'*')
hold on
loglog(sigma1,Z1)
xlabel('|S| (MPa)')
ylabel('Z (1/s)')
legend('experiment','model')
grid on