4panelspringpot.jl

# using RHEOS; using PyPlot; using PyCall; @pyimport seaborn as sns;

# initialise seaborn package with correct parameters
pgf_with_latex = Dict("figure.figsize" => (10.5, 9.5), "mathtext.default" => "regular");
sns.set(rc=pgf_with_latex)
sns.set_style("ticks")

# initialise time data
step_size = 1e-3
t_len = 1.0
t_off = 0.5
loading_data = stepgen(t_len, 0.0; stepsize=step_size)
# initialise frequency data
step_size = 1e-2
ω₀ = step_size
ω₊ = 1e2
chirp = collect(ω₀:step_size:ω₊)
normbool = false

# figure
fig, ax = subplots(2,2)
fig[:subplots_adjust](left=0.09, right=0.96, bottom=0.11, top=0.95, wspace=0.3, hspace=0.235)

# Spring-pot relaxation response
# for α in vcat(collect(0.0:0.1:0.9), [0.99])
for α in collect(0.0:0.1:0.9)
    modelSP = SpringPot([1.0, α])
    result = modelsteppredict(loading_data, modelSP, :G; step_on = 0.0)
    if α==0.0
        ax[1, 1][:plot](result.t, result.σ/maximum(abs.(result.σ)), "-", color="red", zorder=11, linewidth=3)
    # elseif α==0.99
    #     ax[:plot](result.t, result.σ/maximum(abs.(result.σ)), "-", color="blue", label="")
    elseif α==0.3
        ax[1, 1][:plot](result.t[2:end], result.σ[2:end]/maximum(abs.(result.σ[2:end])), "-", color="#dc7633", linewidth=3)
    elseif α==0.7
        ax[1, 1][:plot](result.t[2:end], result.σ[2:end]/maximum(abs.(result.σ[2:end])), "-", color="#229954", linewidth=3)
    elseif α==0.5
        ax[1, 1][:plot](result.t[2:end], result.σ[2:end]/maximum(abs.(result.σ[2:end])), "-", color="#d4ac0d", linewidth=3)
    else
        ax[1, 1][:plot](result.t[2:end], result.σ[2:end]/maximum(abs.(result.σ[2:end])), "-", color="black", alpha=0.1, linewidth=2)
    end
end
# dash-pot relaxation response (special case)
dashpot_σ = [0.0 for i in loading_data.t]
dashpot_σ[1] = 1.0
ax[1, 1][:plot](loading_data.t, dashpot_σ, "-", color="blue", label="", zorder=10, linewidth=3)
#
ax[1, 1][:set_xlabel]("Time, t (s)", fontname="serif", fontsize=13)
ax[1, 1][:set_ylabel]("Normalized Stress, \$\\overline{\\sigma}(t)\$", fontname="serif", fontsize=13)
ax[1, 1][:locator_params](axis ="y", nbins=5)
ax[1, 1][:tick_params]("both", labelsize=13)
ax[1, 1][:grid]("on", alpha=0.4)
ax[1, 1][:annotate]("(a)", xy=(140, -52), fontweight="regular", xycoords="axes points")
ax[1, 1][:spines]["right"][:set_visible](false)
ax[1, 1][:spines]["top"][:set_visible](false)
ax[1, 1][:locator_params](axis="y", nbins=6)

# Spring-pot creep response
for α in collect(0.0:0.1:1.0)
    modelSP = SpringPot([1.0, α])
    loading = modelsteppredict(loading_data, modelSP, :J; step_on = 0.0)
    unloading = modelsteppredict(loading_data, modelSP, :J; step_on = t_off)
    result = loading - unloading
    if α==0.0
        ax[1, 2][:plot](result.t, result.ϵ/maximum(abs.(result.ϵ)), "-", color="red", linewidth=3, zorder=11)
    elseif α==1.0
        ax[1, 2][:plot](result.t, result.ϵ/maximum(abs.(result.ϵ)), "-", color="blue", linewidth=3, zorder=10)
    elseif α==0.5
        ax[1, 2][:plot](result.t, result.ϵ/maximum(abs.(result.ϵ)), "-", color="#d4ac0d", linewidth=3)
    elseif α==0.3
        ax[1, 2][:plot](result.t, result.ϵ/maximum(abs.(result.ϵ)), "-", color="#dc7633", linewidth=3)
    elseif α==0.7
        ax[1, 2][:plot](result.t, result.ϵ/maximum(abs.(result.ϵ)), "-", color="#229954", linewidth=3)
    else  
        ax[1, 2][:plot](result.t, result.ϵ/maximum(abs.(result.ϵ)), "-", color="black", alpha=0.1, linewidth=2)
    end
end
ax[1, 2][:set_xlabel]("Time, t (s)", fontname="serif", fontsize=13)
ax[1, 2][:set_ylabel]("Normalized Strain, \$\\overline{\\epsilon}(t)\$", fontname="serif", fontsize=13)
ax[1, 2][:tick_params]("both", labelsize=13)
ax[1, 2][:grid]("on", alpha=0.4)
ax[1, 2][:annotate]("(b)", xy=(140, -52), fontweight="regular", xycoords="axes points")
ax[1, 2][:spines]["right"][:set_visible](false)
ax[1, 2][:spines]["top"][:set_visible](false)

# G' and G'' [0.0, 0.5, 1.0]
# for β in [0.0, 0.5]
for β in collect(0.0:0.1:0.9)
    model = SpringPot([1.0, β])
    result = model.Gp(chirp, [1.0, β])
    if !all(i -> i==0.0, result) && normbool
        result = result/maximum(result)
    end
    if β==0.0
        ax[2, 1][:loglog](chirp, result, "-", color="red", linewidth=3, zorder=11)
    elseif β==1.0
        ax[2, 1][:loglog](chirp, result, "-", color="blue", linewidth=3, zorder=10)
    elseif β==0.5
        ax[2, 1][:loglog](chirp, result, "-", color="#d4ac0d", linewidth=3)
    end
end
# for β in [0.0, 0.5, 1.0]
for β in collect(0.0:0.1:1.0)
    model = SpringPot([1.0, β])
    result = model.Gpp(chirp, [1.0, β])
    if !all(i -> i==0.0, result) && normbool
        result = result/maximum(result)
    end
    if β==0.0
        ax[2, 1][:loglog](chirp, result, "--", color="red", linewidth=3, zorder=11)
    elseif β==1.0
        ax[2, 1][:loglog](chirp, result, "--", color="blue", linewidth=3, zorder=10)
    elseif β==0.5
        ax[2, 1][:loglog](chirp, result, "--", color="#f7dc6f", linewidth=3)
    end
end
ax[2, 1][:tick_params]("both", labelsize=13)
ax[2, 1][:grid]("on", alpha=0.4, which="major")
ax[2, 1][:set_xlabel]("Frequency, \${\\omega}\$ (\$Hz\$)", fontname="serif", fontsize=13)
ax[2, 1][:set_ylabel]("Storage and Loss moduli, \$G'\$ & \$G''\$ (\$Pa\$)", fontname="serif", fontsize=13)
ax[2, 1][:set_ylim](bottom=1e-2, top=1e2)
ax[2, 1][:annotate]("(c)", xy=(140, -65), fontweight="regular", xycoords="axes points")
ax[2, 1][:spines]["right"][:set_visible](false)
ax[2, 1][:spines]["top"][:set_visible](false)

# G' and G'' [0.3, 0.5, 0.7]
for β in [0.3, 0.5, 0.7]
    model = SpringPot([1.0, β])
    result = model.Gp(chirp, [1.0, β])
    if !all(i -> i==0.0, result) && normbool
        result = result/maximum(result)
    end
    if β==0.3
        ax[2, 2][:loglog](chirp, result, "-", color="#dc7633", linewidth=3, zorder=11)
    elseif β==0.7
        ax[2, 2][:loglog](chirp, result, "-", color="#229954", linewidth=3, zorder=10)
    elseif β==0.5
        ax[2, 2][:loglog](chirp, result, "-", color="#d4ac0d", linewidth=3)
    end
end
for β in [0.3, 0.5, 0.7]
    model = SpringPot([1.0, β])
    result = model.Gpp(chirp, [1.0, β])
    if !all(i -> i==0.0, result) && normbool
        result = result/maximum(result)
    end
    if β==0.3
        ax[2, 2][:loglog](chirp, result, "--", color="#dc7633", linewidth=3, zorder=11)
    elseif β==0.7
        ax[2, 2][:loglog](chirp, result, "--", color="#229954", linewidth=3, zorder=10)
    elseif β==0.5
        ax[2, 2][:loglog](chirp, result, "--", color="#f7dc6f", linewidth=3)
    end
end
ax[2, 2][:set_xlabel]("Frequency, \${\\omega}\$ (\$Hz\$)", fontname="serif", fontsize=13)
ax[2, 2][:set_ylabel]("Storage and Loss moduli, \$G'\$ & \$G''\$ (\$Pa\$)", fontname="serif", fontsize=13)
ax[2, 2][:tick_params]("both", labelsize=13)
ax[2, 2][:grid]("on", alpha=0.4, which="major")
ax[2, 2][:set_ylim](bottom=1e-2, top=1e2)
ax[2, 2][:annotate]("(d)", xy=(140, -65), fontweight="regular", xycoords="axes points")
ax[2, 2][:spines]["right"][:set_visible](false)
ax[2, 2][:spines]["top"][:set_visible](false)