clc
clear
close all

%% INPUT MATERIAL PARAMETERS %%
A = 8*10^18;
C = 14*10^6;
n = 4;
Q = 450000;
Tm = 1350;
mu_0 = 54000*10^6;
mu_1 = -17*10^6;
k = 8.314;
delta = 0.035;
f = 0.6;
h = 50*10^6;
l = 3;
D_0 = C/100;
poisson = 0.34;
alpha = 16.5*10^-6;
delta_T = 0;

T_0 = 538+273;
rho = 9.131*10^3;
c = 3860;

T = 538+273;
if  T>0.5*Tm
    v = exp(-Q/(k*T));
end 

if T<0.5*Tm
    v = exp(-Q/(k*0.5*Tm)*(log(0.5*Tm/T)+1));
end

mu = mu_0+mu_1*T;
E = 2*mu*(1+poisson);

epsdot_f = 0.02;
epsdot_s = 10^-9;

n_hat = log2(epsdot_f/epsdot_s)/log2(asinh((epsdot_f/(A*v))^(1/n))/asinh((epsdot_s/(A*v))^(1/n)));
eps_a = exp(0.5*(log(epsdot_f)+log(epsdot_s)));
v_hat = eps_a/(asinh((eps_a/(A*v))^(1/n)))^n_hat;

%% CALCULATE CYCLIC STRESS %%
i = 0;
j = 0;
m = 0;

% initial guess value
sigma = 0; 
beta = 0;
epsdot_p_mag = 0;
epsdot_p = sign(sigma-beta)*epsdot_p_mag/sqrt(3);
D = D_0;
L = 0;
xi = 0;
lambda = l*xi;
r = 0;
H = 0.1*mu;

delta_t = 0.01;

    for t = 0:0.01:300;
        i = i+1
         
if  T>0.5*Tm
    v = exp(-Q/(k*T));
end 

if T<0.5*Tm
    v = exp(-Q/(k*0.5*Tm)*(log(0.5*Tm/T)+1));
end

mu = mu_0+mu_1*T;
E = 2*mu*(1+poisson);

n_hat = log2(epsdot_f/epsdot_s)/log2(asinh((epsdot_f/(A*v))^(1/n))/asinh((epsdot_s/(A*v))^(1/n)));
eps_a = exp(0.5*(log(epsdot_f)+log(epsdot_s)));
v_hat = eps_a/(asinh((eps_a/(A*v))^(1/n)))^n_hat;

iter = 1;
fun = 1;

%% INPUT STRAIN %%
eps_dot = 0.01;

a = 0.013; % max strain
b = a/eps_dot;

if mod(fix(t/b),4) == 3  
    eps = -(fix(t/b)+1)*a+eps_dot*t;
end
if mod(fix(t/b),4) == 0
    eps = -fix(t/b)*a+eps_dot*t;
end
if mod(fix(t/b),4) == 1
    eps = (fix(t/b)+1)*a-eps_dot*t;
end
if mod(fix(t/b),4) == 2
    eps = fix(t/b)*a-eps_dot*t;
end

%% VALUES FROM LAST TIME STEP %%
L = f*(D-D_0)*(C-D)/(delta*C);

if L == 0
    xi = 0;
end
if L ~= 0
xi = 0.5*abs(sign(sigma-beta)-beta/(sqrt(3)*L));
end

H = (0.1+0.9*xi)*mu;

lambda = l*xi;
r = v_hat*sign(((D-D_0)/(delta*C)))*(abs((D-D_0)/(delta*C)))^n_hat;

if L == 0 
    beta_dot = 0;
end
if L ~= 0
beta_dot = 3*H*(epsdot_p-beta/(3*L)*epsdot_p_mag);
end

D_dot = h*((1-lambda)*epsdot_p_mag-r);
    
if i == 2
    S_magnitude = sqrt(1/3)*4.3108*10^6;   %initial magnitude of deviatoric stress
    beta = f*(C-D_0-delta*S_magnitude)*S_magnitude/C/sqrt(1/3);
    D = D_0+delta*S_magnitude;
end

%% VALUES FOR CURRENT TIME STEP %%
beta = beta+delta_t*beta_dot;
D = D+D_dot*delta_t;
Y = (1-f)*(D-D_0)*(C-D)/(delta*C);

%% ITERATION %%
while abs(fun)>=10^-4 & iter<400
   
if (abs(sigma-beta)-sqrt(3)*Y)/(sqrt(3)*D) >= 0
Z = ((abs(sigma-beta)-sqrt(3)*Y)/(sqrt(3)*D))^n_hat;

epsdot_p_mag = v_hat*Z;
epsdot_p = epsdot_p_mag/sqrt(3)*sign(sigma-beta);

if i==1
    eps_p(i) = 0;
end
if i>1
eps_p(i) = delta_t*epsdot_p+eps_p(i-1);
end

    fun = sigma-E*(eps-eps_p(i)-alpha*delta_T);

if sigma > beta
    dfun = 1+E*delta_t/sqrt(3)*n_hat*v_hat*sign((sigma-beta-sqrt(3)*Y)/(sqrt(3)*D))*(abs((sigma-beta-sqrt(3)*Y)/(sqrt(3)*D)))^(n_hat-1)*(1/(sqrt(3)*D));
end

if sigma < beta
    dfun = 1-E*delta_t/sqrt(3)*n_hat*v_hat*sign((-sigma+beta-sqrt(3)*Y)/(sqrt(3)*D))*(abs((-sigma+beta-sqrt(3)*Y)/(sqrt(3)*D)))^(n_hat-1)*(-1/(sqrt(3)*D));
end 

if sigma == beta
    dfun = 1;
end

end

if (abs(sigma-beta)-sqrt(3)*Y)/(sqrt(3)*D) < 0
Z = 0;

epsdot_p_mag = 0;
epsdot_p = 0;

if i==1
    eps_p(i) = 0;
end
if i>1
eps_p(i) = eps_p(i-1);
end

    fun = sigma-E*(eps-eps_p(i)-alpha*delta_T);
    dfun = 1;
end

    sigma_1 = sigma;
    sigma = sigma_1-fun/dfun;

    iter = iter+1;

end

trial_3(i) = epsdot_p_mag; % don't change
trial_4(i) = beta_dot;

stress(i) = sigma_1/10^6;
strain(i) = eps;
time(i) = t;
% i = 0;

%% DISSIPATION RATE %%
if mod(fix(t/b),4) == 3 || mod(fix(t/b),4) == 0
    W_dot(i) = sigma_1*eps_dot;
end

if mod(fix(t/b),4) == 1 || mod(fix(t/b),4) == 2
    W_dot(i) = sigma_1*(-eps_dot);
end

%% ENERGY DISSIPATION %%
W(1) = 0;
if i>=2
t_1(i) = t;   % time vector
W(i) = trapz(t_1,W_dot);
end

%% TEMPERATURE INCREASE %%
Tem(i) = T_0+W(i)/(rho*c);
T = Tem(i);
delta_T = T-T_0;
trial(i) = delta_T;

%% FIND STRESS RANGE %%

if mod(t,(b*4)) == 0
   m = m+1;
   j = 0;
end

j = j+1;
stress_1(j,m) = sigma_1/10^6;
strain_1(j,m) = eps;
cycle(m) = m;

    end
    
%%    
figure
plot(strain,stress,'LineWidth',2)
xlabel('strain \epsilon')
ylabel('stress \sigma (MPa)')
axis([-0.015 0.015 -170 170])
grid on

%% 
max = max(stress_1);
min = min(stress_1);
stress_range = (max-min);
stress_range(m) = [];
cycle(m) = [];

%%
figure
semilogx(cycle,stress_range*10^6*pi*(6.35*10^-3)^2/4/10^3,'LineWidth',2);
axis([1 1000 0 15])
xlabel('number of cycles')
ylabel('force range, kN (F = \sigmaA)')
set(gca,'YMinorGrid','on')
grid on

figure
plot(time,W,'LineWidth',2)
xlabel('time')
ylabel('energy dissipation')
grid on

figure
plot(time,Tem,'LineWidth',2)
xlabel('time')
ylabel('temperature (K)')
grid on
