clc
clear
close all
%% LOADING EXPERIMENTAL RESULT %%
load('fitting_1')

%% ADJUST PARAMETERS FOR ALL LOADING CONDITIONS %%
eps_t1 = 0.0006;
c1 = 11800;

eps_t2 = 0.0007;
c2 = 7000;

eps_t3 = 0.0018;
c3 = 2800;

eps_t4 = 0.0032;
c4 = 2600;

%% STEADY-STATE STRAIN RATES FROM EXPERIMENTAL RESULTS %%
eps_ss1 = 4.1406*10^-8;
eps_ss2 = 6.0938*10^-8;
eps_ss3 = 9.01961*10^-8;
eps_ss4 = 1.46299*10^-7;

t = 0:19000;

%% POWER-LAW MODEL %%
eps1 = 6.3*10^-3+eps_t1*(1-exp(-c1*eps_ss1*t))+eps_ss1*t;
eps2 = 9.6*10^-3+eps_t2*(1-exp(-c2*eps_ss2*t))+eps_ss2*t;
eps3 = 14.6*10^-3+eps_t3*(1-exp(-c3*eps_ss3*t))+eps_ss3*t;
eps4 = 18.4*10^-3+eps_t4*(1-exp(-c4*eps_ss4*t))+eps_ss4*t;

t = [0 t];
eps1 = [0 eps1];
eps2 = [0 eps2];
eps3 = [0 eps3];
eps4 = [0 eps4];

%% 
sigma_cr = 47*10^6;
R = 8.314;
time = 5*3600;

%% MATERIAL PARAMETERS %%
A = 3*10^-10;

i = 0;
j = 0;
k = 0;
m = 0;

for time = [0 0.01 0.1 0.2 0.5 1 2 3 4 5]*3600;
    m = m+1;
for stress = [14.4 21.6 28.8 33.6]*10^6;
    i = i+1;
    for Temp = [523 523 523 523];
        j = j+1;
        
% power-law exponent n          
n = 3.1943019-0.002283*Temp;
        
% fit activition energy using bilinear interpolation
if Temp<=548
    Q = 123*10^3;
end

if Temp>548 && Temp<=558
    if stress<=0.45*sigma_cr
        Q = 123*10^3;
    end
    if stress>0.45*sigma_cr && stress<=0.6*sigma_cr
        a = 123-(stress/sigma_cr-0.45)/0.15*3;
        Q = (123-(Temp-548)/10*(123-a))*10^3;
    end
    if stress>0.6*sigma_cr
        a = 120-(stress/sigma_cr-0.6)/0.1*4;
        Q = (123-(Temp-548)/10*(123-a))*10^3;
    end
end

if Temp>558
    if stress<=0.3*sigma_cr
        Q = 123*10^3;
    end
    if stress>0.3*sigma_cr && stress<=0.45*sigma_cr
        a = 123-(stress/sigma_cr-0.3)/0.15*4;
        Q = (123-(Temp-558)/15*(123-a))*10^3;
    end
    if stress>0.45*sigma_cr && stress<=0.6*sigma_cr
        a = 119-(stress/sigma_cr-0.45)/0.15*5;
        b = 123-(stress/sigma_cr-0.45)/0.15*3;
        Q = (b-(Temp-558)/15*(b-a))*10^3;
    end
    if stress>0.6*sigma_cr
        a = 114-(stress/sigma_cr-0.6)/0.1*8;
        b = 120-(stress/sigma_cr-0.6)/0.1*4;
        Q = (b-(Temp-558)/15*(b-a))*10^3;
    end
end

x(i,j) = Temp;
y(i,j) = log10(stress);
eps_ss = A*stress^n*exp(-Q/(R*Temp));
Q1(i,j) = Q;

% elastic modulus
slope = 0.0068*Temp^2-7.1754*Temp+1861.8;
E_ave = -1.2148*10^7*Temp+8.5038*10^9;
E = slope.*(stress-0.51*sigma_cr)+E_ave;
eps_elastic = stress/E;

% transient constant
if Temp<=548
    C = 10^(6.6184)*(stress/10^6)^(-2.1575);
end

if Temp>548
    T1 = 548;
    a = -3.8947+(573-Temp)/(573-T1)*1.7372;
    b = 8.3528-(573-Temp)/(573-T1)*1.7344;
    C = 10^(b)*(stress/10^6)^(a);
end

% transient strain
eps_t = ((1.022*10^-9*exp(0.0227*Temp))*exp((0.00689*exp(0.00487*Temp))*(stress/10^6)));

% plastic strain calcualte by power-law model
eps_plastic = eps_t*(1-exp(-C*eps_ss*time))+eps_ss*time;

% total creep strain
eps(i,j) = eps_elastic+eps_plastic;

% total creep strain at loading conditions from experimental data
if i == j
    k = k+1;
    eps_1(k) = eps(i,j);
end
        
    end
    
    j = 0;
    
end

% total creep strain at loading conditions from experimental data
eps_total(:,m) = eps_1;
k = 0;
i = 0;

end

t_fit = [0 0.01 0.1 0.2 0.5 1 2 3 4 5]*3600;

%% PLOT FIGURE %%
figure
plot(time1(1:5:end),strain1(1:5:end),'ro','LineWidth',2)
hold on
plot(t_fit,eps_total(1,:),'k--','LineWidth',2)
hold on
% figure
plot(time2(1:5:end),strain2(1:5:end),'ro','LineWidth',2)
hold on

% figure
plot(time3(1:4:end),strain3(1:4:end),'ro','LineWidth',2)
hold on

axis([0 18000 0 0.027])
% xlabel('time (s)')
ylabel('strain')
legend('experiment','model') 

% figure
plot(time4(1:6:end),strain4(1:6:end),'ro','LineWidth',2)
hold on

plot(t_fit,eps_total(2,:),'k--','LineWidth',2)
hold on
plot(t_fit,eps_total(3,:),'k--','LineWidth',2)
hold on
plot(t_fit,eps_total(4,:),'k--','LineWidth',2)
hold on

text(9000,0.0055,'0.3\sigma_{cr}')
text(9000,0.013,'0.45\sigma_{cr}')
text(9000,0.019,'0.6\sigma_{cr}')
text(9000,0.024,'0.7\sigma_{cr}')
text(7000,0.026,'Temperature = 523K')

% save('fitting_1')