# *** material constants ***
E = 200.0e9
nu = 0.3
G = E/(2*(1+nu))
SS = 0.5e6
s = 2
C = 2.0e9
sigma_f = 200.0e6
sigma_ultimate = 600.0e6
b = 2/15 * (4-5*nu)/(1-nu)
# b = 0.9999
Dc = 1
h = 0.2   #amount of damage activation in compression#
epsilon_pd = 0.1
# println("b = " , b)
# a = (1-beta)/(1+nu)
# k = (3*a*E)/(2*C)



# *** loading ***
delta_strain = 0.00304
strain_mean = +0.0000
num_substeps = 100
num_cycles = 10000000
Damage_0 = 0.0
acc_pl_strain_micro = 0 
Back_str_0 = 0


# *** Damage threshold ***
sigma_max_meso = (strain_mean + delta_strain / 2) * E
delta_sigma_meso = delta_strain * E
strain_f = sigma_f / E
strain_max = strain_mean + delta_strain / 2
strain_min = strain_mean - delta_strain / 2
strain_C = max(abs(strain_max), abs(strain_min))
sigma_max_micro = sigma_f + (strain_C - strain_f) * C
sigma_s = sigma_max_micro
# sigma_s = ((abs(delta_sigma_meso) - 2*sigma_f) * C / E + 2*sigma_f)/2
PD = epsilon_pd * (sigma_ultimate-sigma_f)/(sigma_s-sigma_f)

# *** NO Damage threshold ***
# PD = 0.0


function loadpath(strain_mean, delta_strain, num_substeps)
    # *** strain load cycle generation ***
    strain_max = strain_mean + delta_strain / 2  
    println("strain_max = ", strain_max)

    strain_min = strain_mean - delta_strain / 2
    println("strain_min = ", strain_min)

    println("strain_mean = ", strain_mean)

    num_path_steps = trunc(Int, (num_substeps/4))
    println("num_path_steps = ", num_path_steps)

    strain_icrement = (delta_strain / 2) / num_path_steps
    println("strain_icrement = ", strain_icrement)

    substeps = 4 * num_path_steps
    strain_cycle = zeros(Float64, substeps)

    strain_cycle[1] = strain_mean + strain_icrement
    for step in 2: num_path_steps
        strain_cycle[step] = strain_cycle[step-1] + strain_icrement
    end   
    count = num_path_steps
    # println(strain_cycle, count)

    for step in 1: num_path_steps
        strain_cycle[count + step] = strain_cycle[count + step-1] - strain_icrement  
    end
    count = 2 * num_path_steps   
    # println(strain_cycle, count)

    for step in 1: num_path_steps
        strain_cycle[count + step] = strain_cycle[count + step-1] - strain_icrement
    end  
    count = 3 * num_path_steps   
    # println(strain_cycle, count)

    for step in 1: num_path_steps
        strain_cycle[count + step] = strain_cycle[count + step-1] + strain_icrement
    end 
    count = 4 * num_path_steps   
    # println(strain_cycle, count) 

    return(strain_cycle, substeps) 

end

function my_plot(cycle_m_1, cycle_m_2, cycle_m_3, cycle_m_4, cycle, substeps, eps_meso,sigma_meso,eps_micro,sigma_micro)

    cycle_1_fs = 1
    cycle_1_ls = substeps + 1

    cycle_m_1_fs = cycle_m_1*substeps + 1
    cycle_m_1_ls = (cycle_m_1 + 1)*substeps + 1

    cycle_m_2_fs = cycle_m_2*substeps + 1
    cycle_m_2_ls = (cycle_m_2 + 1)*substeps + 1

    cycle_m_3_fs = cycle_m_3*substeps + 1
    cycle_m_3_ls = (cycle_m_3 + 1)*substeps + 1

    cycle_m_4_fs = cycle_m_4*substeps + 1
    cycle_m_4_ls = (cycle_m_4 + 1)*substeps + 1

    cycle_fs = (cycle-2)*substeps + 1
    cycle_ls = (cycle-1)*substeps + 1

    plot(eps_meso[cycle_1_fs:cycle_1_ls], sigma_meso[cycle_1_fs:cycle_1_ls], xlabel="strain", ylabel="stress", label="meso", title="stress_strain", grid=true)
    plot!(eps_micro[cycle_1_fs:cycle_1_ls], sigma_micro[cycle_1_fs:cycle_1_ls], xlabel="strain", ylabel="stress", label="micro_1", title="stress_strain", grid=true)
    plot!(eps_micro[cycle_m_1_fs:cycle_m_1_ls], sigma_micro[cycle_m_1_fs:cycle_m_1_ls], xlabel="strain", ylabel="stress", label="micro_m_1", title="stress_strain", grid=true)
    plot!(eps_micro[cycle_m_2_fs:cycle_m_2_ls], sigma_micro[cycle_m_2_fs:cycle_m_2_ls], xlabel="strain", ylabel="stress", label="micro_m_2", title="stress_strain", grid=true)
    plot!(eps_micro[cycle_m_3_fs:cycle_m_3_ls], sigma_micro[cycle_m_3_fs:cycle_m_3_ls], xlabel="strain", ylabel="stress", label="micro_m_3", title="stress_strain", grid=true)
    plot!(eps_micro[cycle_m_4_fs:cycle_m_4_ls], sigma_micro[cycle_m_4_fs:cycle_m_4_ls], xlabel="strain", ylabel="stress", label="micro_m_4", title="stress_strain", grid=true)
    plot!(eps_micro[cycle_fs:cycle_ls], sigma_micro[cycle_fs:cycle_ls], xlabel="strain", ylabel="stress", label="micro_last", title="stress_strain", seriestype=:scatter)

end


function plot_yield_surface(sigma_max_micro, sigma_min_micro, D_cycle, cycle, sigma_f)
    tmin = 0
    tmax = 2π
    tmin_t = -3π/4
    tmax_t = 3π/4
    tmin_c = 3π/4
    tmax_c = 5π/4
    θ = π/4
    # θ = 0
    tvec = range(tmin, tmax, length = 100)
    tvec_t = range(tmin_t, tmax_t, length = 100)
    tvec_c = range(tmin_c, tmax_c, length = 100)

    cycle_y_1 = trunc(Int, (cycle * 0.6))
    cycle_y_2 = trunc(Int, (cycle * 0.7))
    cycle_y_3 = trunc(Int, (cycle * 0.8))
    cycle_y_4 = trunc(Int, (cycle * 0.9))


    a = sqrt(2) * sigma_f
    b = sqrt(6)/3 * sigma_f
    a_t_1 = abs(sqrt(2) * sigma_max_micro[cycle_y_1])
    b_t_1 = abs(sqrt(6)/3 * sigma_max_micro[cycle_y_1])
    a_c_1 = abs(sqrt(2) * sigma_min_micro[cycle_y_1])
    b_c_1 = abs(sqrt(6)/3 * sigma_min_micro[cycle_y_1])

    a_t_2 = abs(sqrt(2) * sigma_max_micro[cycle_y_2])
    b_t_2 = abs(sqrt(6)/3 * sigma_max_micro[cycle_y_2])
    a_c_2 = abs(sqrt(2) * sigma_min_micro[cycle_y_2])
    b_c_2 = abs(sqrt(6)/3 * sigma_min_micro[cycle_y_2])

    a_t_3 = abs(sqrt(2) * sigma_max_micro[cycle_y_3])
    b_t_3 = abs(sqrt(6)/3 * sigma_max_micro[cycle_y_3])
    a_c_3 = abs(sqrt(2) * sigma_min_micro[cycle_y_3])
    b_c_3 = abs(sqrt(6)/3 * sigma_min_micro[cycle_y_3])

    a_t_4 = abs(sqrt(2) * sigma_max_micro[cycle_y_4])
    b_t_4 = abs(sqrt(6)/3 * sigma_max_micro[cycle_y_4])
    a_c_4 = abs(sqrt(2) * sigma_min_micro[cycle_y_4])
    b_c_4 = abs(sqrt(6)/3 * sigma_min_micro[cycle_y_4])

    println(D_cycle[1], "   ", D_cycle[cycle_y_1], "   ", D_cycle[cycle_y_2], "   ", D_cycle[cycle_y_3], "   ", D_cycle[cycle_y_4], "   ", D_cycle[cycle])


    plot(((sqrt(a)*cos(θ)*cos.(tvec))-(sqrt(b)*sin(θ)*sin.(tvec))),((sqrt(a)*sin(θ)*cos.(tvec))+(sqrt(b)*cos(θ)*sin.(tvec))), aspect_ratio=:equal, label=D_cycle[1])
    plot!(((sqrt(a_t_1)*cos(θ)*cos.(tvec_t))-(sqrt(b_t_1)*sin(θ)*sin.(tvec_t))),((sqrt(a_t_1)*sin(θ)*cos.(tvec_t))+(sqrt(b_t_1)*cos(θ)*sin.(tvec_t))), label=D_cycle[cycle_y_1])
    plot!(((sqrt(a_c_1)*cos(θ)*cos.(tvec_c))-(sqrt(b_c_1)*sin(θ)*sin.(tvec_c))),((sqrt(a_c_1)*sin(θ)*cos.(tvec_c))+(sqrt(b_c_1)*cos(θ)*sin.(tvec_c))), label=D_cycle[cycle_y_1])
    plot!((((sqrt(a)*cos(θ)*cos.(tvec_t))-(sqrt(b)*sin(θ)*sin.(tvec_t)))*(1- D_cycle[cycle_y_1])^2),(((sqrt(a)*sin(θ)*cos.(tvec_t))+(sqrt(b)*cos(θ)*sin.(tvec_t)))*(1- D_cycle[cycle_y_1])^2), label=D_cycle[cycle_y_1])

    # plot!(((sqrt(a_t_2)*cos(θ)*cos.(tvec_t))-(sqrt(b_t_2)*sin(θ)*sin.(tvec_t))),((sqrt(a_t_2)*sin(θ)*cos.(tvec_t))+(sqrt(b_t_2)*cos(θ)*sin.(tvec_t))), label=D_cycle[cycle_y_2])
    # plot!(((sqrt(a_c_2)*cos(θ)*cos.(tvec_c))-(sqrt(b_c_2)*sin(θ)*sin.(tvec_c))),((sqrt(a_c_2)*sin(θ)*cos.(tvec_c))+(sqrt(b_c_2)*cos(θ)*sin.(tvec_c))), label=D_cycle[cycle_y_2])
    # plot!(((sqrt(a_t_3)*cos(θ)*cos.(tvec_t))-(sqrt(b_t_3)*sin(θ)*sin.(tvec_t))),((sqrt(a_t_3)*sin(θ)*cos.(tvec_t))+(sqrt(b_t_3)*cos(θ)*sin.(tvec_t))), label=D_cycle[cycle_y_3])
    # plot!(((sqrt(a_c_3)*cos(θ)*cos.(tvec_c))-(sqrt(b_c_3)*sin(θ)*sin.(tvec_c))),((sqrt(a_c_3)*sin(θ)*cos.(tvec_c))+(sqrt(b_c_3)*cos(θ)*sin.(tvec_c))), label=D_cycle[cycle_y_3])
    # plot!(((sqrt(a_t_4)*cos(θ)*cos.(tvec_t))-(sqrt(b_t_4)*sin(θ)*sin.(tvec_t))),((sqrt(a_t_4)*sin(θ)*cos.(tvec_t))+(sqrt(b_t_4)*cos(θ)*sin.(tvec_t))), label=D_cycle[cycle_y_4])
    # plot!(((sqrt(a_c_4)*cos(θ)*cos.(tvec_c))-(sqrt(b_c_4)*sin(θ)*sin.(tvec_c))),((sqrt(a_c_4)*sin(θ)*cos.(tvec_c))+(sqrt(b_c_4)*cos(θ)*sin.(tvec_c))), label=D_cycle[cycle_y_4])



end

strain_cycle, substeps = loadpath(strain_mean, delta_strain, num_substeps)
# println(strain_cycle, substeps)

# *** varibles difinition ***
eps_meso = Float64[] 
sigma_meso = Float64[] 
eps_micro = Float64[] 
eps_el_micro = Float64[] 
eps_pl_micro = Float64[] 
sigma_micro = Float64[] 
sigma_bar_micro = Float64[] 
beta_micro = Float64[] 
D_micro = Float64[]
D_cycle = Float64[]
sigma_max_micro = Float64[] 
sigma_min_micro = Float64[] 

# *** initail values of variables ***
step = 1
cycle = 0
Damage = 0


append!(eps_meso, strain_mean)
append!(sigma_meso, E*strain_mean)
append!(eps_micro, strain_mean)
append!(eps_el_micro, strain_mean)
append!(eps_pl_micro, 0)
append!(sigma_micro, E*strain_mean)
append!(sigma_bar_micro, E*strain_mean)
append!(beta_micro, Back_str_0)
append!(D_micro, Damage_0)


for cyc in 1:num_cycles
    global cycle = cyc
    global sigma_max_mic = 0
    global sigma_min_mic = 0
    global delta_sigma_meso = 0
    for strain in strain_cycle

        global step = step + 1
       
        # println("cycle: ", cycle, "         step: ", step)
        
        # println("delta_p_micro = ", delta_p_micro)
        # *** Elastic solution in meso scale ***
        eps_meso_current = strain
        # println("eps_meso_current: ", eps_meso_current)
        delta_eps_meso = eps_meso_current - eps_meso[step-1]
        delta_sigma_meso = E * delta_eps_meso
        sigma_meso_current = sigma_meso[step-1] + delta_sigma_meso 
        # println("sigma_meso_current: ",sigma_meso_current)

        # *** Transition from meso scale to micro scale ***
        # *** trial elastic micro step ***
        delta_p_micro = 0
        delta_eps_micro = delta_eps_meso
        eps_pl_micro_trail = eps_pl_micro[step-1]
        delta_eps_pl_micro = 0.0
        delta_eps_el_micro = delta_eps_micro
        D_micro_trail = D_micro[step-1]
        beta_micro_trail = beta_micro[step-1]
        eps_el_micro_trail = eps_el_micro[step-1] + delta_eps_el_micro
        sigma_bar_micro_trail = sigma_bar_micro[step-1] + E * delta_eps_el_micro
        phi_micro_trail = abs(sigma_bar_micro_trail - beta_micro_trail) - sigma_f
        sign_current = sign(sigma_bar_micro_trail - beta_micro_trail)
        # println("sign_current = ", sign_current)

        if phi_micro_trail <= 0
            # println("elastic")
            eps_pl_micro_current = eps_pl_micro_trail
            # println("eps_pl_micro_current: ", eps_pl_micro_current)
            eps_el_micro_current = eps_el_micro_trail
            # println("eps_el_micro_current: ", eps_el_micro_current)
            eps_micro_current = eps_el_micro_current + eps_pl_micro_current
            # println("eps_micro_current: ", eps_micro_current)
            D_micro_current = D_micro_trail
            # println("D_micro_current: ", D_micro_current)
            beta_micro_current = beta_micro_trail
            # println("beta_micro_current: ", beta_micro_current)
            sigma_bar_micro_current = sigma_bar_micro_trail
            # println("sigma_bar_micro_current: ", sigma_bar_micro_current)
            if sigma_bar_micro_current >= 0
                sigma_micro_current = (1-D_micro_current) * sigma_bar_micro_current
            else 
                sigma_micro_current = (1-h*D_micro_current) * sigma_bar_micro_current
            end
            # sigma_micro_current = (1-D_micro_current) * sigma_bar_micro_current
            # println("sigma_micro_current = ", sigma_micro_current)
            if  sigma_max_mic < sigma_micro_current     
                global sigma_max_mic = sigma_micro_current
            end
            if  sigma_min_mic > sigma_micro_current
                global sigma_min_mic = sigma_micro_current
            end
            # println("max-min:    ",  sigma_max_mic,"     ",  sigma_min_mic)
        else
            # *** plastic micro corrector ***
            # println("plastic")
            Qs = -(sigma_bar_micro[step-1] + delta_sigma_meso) + beta_micro[step-1]
            # println("Qs = ", Qs)
            q = (2/3*C*(1-D_micro[step-1]) + E*(1-b))*sign_current
            # println("q = ", q)
            S_micro_current = sign_current * sigma_f
            # println("S_micro_current = ", S_micro_current)
            delta_p_micro = (-1/q)*(Qs + S_micro_current)
            # println("delta_p_micro = ", delta_p_micro)
            delta_eps_pl_micro = delta_p_micro * sign_current
            eps_pl_micro_current = eps_pl_micro[step-1] + delta_eps_pl_micro
            # println("eps_pl_micro_current = ", eps_pl_micro_current)
            beta_micro_current = beta_micro[step-1] + 2/3* C*(1-D_micro[step-1])*delta_p_micro *sign_current
            # println("beta_micro_current = ",  beta_micro_current)
            sigma_bar_micro_current = S_micro_current + beta_micro_current
            # println("sigma_bar_micro_current = ", sigma_bar_micro_current)

            # eps_el_micro_current = eps_el_micro[step-1] + delta_eps_meso + (b-1) * delta_p_micro * sign_current
            # eps_el_micro_current = sigma_bar_micro_current / E
            # eps_el_micro_current = eps_meso_current + (b-1) * eps_pl_micro_current
            delta_eps_el_micro = delta_eps_meso + (b-1) * delta_eps_pl_micro
            eps_el_micro_current = eps_el_micro[step-1] + delta_eps_el_micro
            # println("eps_el_micro_current = ", eps_el_micro_current)
            eps_micro_current = eps_el_micro_current + eps_pl_micro_current
            # println("eps_micro_current = ", eps_micro_current)
            # *** damage deactivation effect ***
            if sigma_bar_micro_current >= 0
                Y_micro_cuurent = 1/(2*E) * (sigma_bar_micro_current)^2
            else 
                Y_micro_cuurent = 1/(2*E) * (sigma_bar_micro_current)^2 * h*(1-D_micro_trail)^2/(1-h*D_micro_trail)^2
            end
            # println("Y_micro_cuurent = ", Y_micro_cuurent)
            
            if (acc_pl_strain_micro + delta_p_micro) >= PD
                delta_D_micro = (Y_micro_cuurent / SS)^s * delta_p_micro
            else 
                delta_D_micro = 0.0
            end        
            # println("delta_D_micro = ", delta_D_micro)
            D_micro_current = D_micro[step-1] + delta_D_micro
            # println("D_micro_current = ", D_micro_current)
            if sigma_bar_micro_current >= 0
                sigma_micro_current = (1-D_micro_current) * sigma_bar_micro_current
            else 
                sigma_micro_current = (1-h*D_micro_current) * sigma_bar_micro_current
            end
            # sigma_micro_current = (1-D_micro_current) * sigma_bar_micro_current
            # println("sigma_micro_current = ", sigma_micro_current)
            if  sigma_max_mic < sigma_micro_current     
                global sigma_max_mic = sigma_micro_current
            end
            if  sigma_min_mic > sigma_micro_current
                global sigma_min_mic = sigma_micro_current
            end
            # println("max-min:    ",  sigma_max_mic,"     ",  sigma_min_mic)
        
        end
        
        # *** update values of variables ***
        append!(eps_meso, eps_meso_current)
        append!(sigma_meso, sigma_meso_current)
        append!(eps_micro, eps_micro_current)
        append!(eps_el_micro, eps_el_micro_current)
        append!(eps_pl_micro, eps_pl_micro_current)
        append!(sigma_micro, sigma_micro_current)
        append!(sigma_bar_micro, sigma_bar_micro_current)
        append!(beta_micro, beta_micro_current)
        append!(D_micro, D_micro_current)
        global Damage = D_micro_current
        global Back_str = beta_micro_current
        global acc_pl_strain_micro = acc_pl_strain_micro + delta_p_micro  

    end
    append!(sigma_max_micro, sigma_max_mic)
    # println(cycle, "   max   ",sigma_max_mic)
    append!(sigma_min_micro, sigma_min_mic)
    # println(cycle, "   min   ",sigma_min_mic)
    append!(D_cycle, Damage)
    if  Damage >= Dc
        break
    end
end

println("cycle: ", cycle, "         step: ", step)
println("acc_pl_strain_micro: ", acc_pl_strain_micro)
println("Damage_micro: ", Damage)
println("Back_stress_micro: ", Back_str)


# cycle_m_1 = trunc(Int, (cycle * 0.6))
# cycle_m_2 = trunc(Int, (cycle * 0.7))
# cycle_m_3 = trunc(Int, (cycle * 0.8))
# cycle_m_4 = trunc(Int, (cycle * 0.9))

# cycle_m_1 = 1
# cycle_m_2 = 2
# cycle_m_3 = 3
# cycle_m_4 = 4

using Plots
# my_plt = my_plot(cycle_m_1, cycle_m_2, cycle_m_3, cycle_m_4, cycle, substeps, eps_meso,sigma_meso,eps_micro,sigma_micro)
# my_plt_2 = plot_yield_surface(sigma_max_micro, sigma_min_micro, D_cycle, cycle, sigma_f)

cycle_1_fs = 1
cycle_1_ls = substeps + 1
cycle_fs = (cycle-2)*substeps + 1
cycle_ls = (cycle-1)*substeps + 1

# for i in cycle_fs:cycle_ls
#     println(i, "   ", eps_micro[i], "   ", sigma_micro[i], "   ", sigma_bar_micro[i]  )
# end


plot(eps_meso[cycle_1_fs:cycle_1_ls], sigma_meso[cycle_1_fs:cycle_1_ls], xlabel="strain", ylabel="stress", label="meso", title="stress_strain", grid=true)
plot!(eps_micro[cycle_1_fs:cycle_1_ls], sigma_micro[cycle_1_fs:cycle_1_ls], xlabel="strain", ylabel="stress", label="micro_1", title="stress_strain", grid=true)
plot!(eps_micro[cycle_fs:cycle_ls], sigma_bar_micro[cycle_fs:cycle_ls], xlabel="strain", ylabel="stress", label="micro_bar_last", title="stress_strain", markershape=:star5)
plot!(eps_micro[cycle_fs:cycle_ls], sigma_micro[cycle_fs:cycle_ls], xlabel="strain", ylabel="stress", label="micro_last", title="stress_strain", markershape=:circle)

# plot(1:cycle,  sigma_max_micro[1:cycle], xlabel="cycle", ylabel="sigma", label="sigma_max", title="max-min", grid=true)
# plot!(1:cycle,  sigma_min_micro[1:cycle], xlabel="cycle", ylabel="sigma", label="sigma_min", title="max-min", grid=true)

# plot(1:cycle,  D_cycle[1:cycle], xlabel="cycle", ylabel="Damage", label="Damage", title="Damage_Cycle", grid=true)



کد بالا رو به پایتون تبدیل کن .\دقت کن اشتباه نکنی تا ارور نگیرم