Plot Weizsäcker-Bethe Mass Formula

WeizsackerBethePlot.jl

"""
Weizsäcker-Betheの質量公式
$$
    B(Z,A)
    = a_1 A
    - a_2 A^{2/3}
    - a_3 \frac{(A-2Z)^2}{A}
    - a_4 \frac{Z^2}{A^{1/3}}
    + \delta (A)
$$
をプロットする.
係数$a_1$-$a_4$は,
$$
    a_1 = 15.7   MeV,
    a_2 = 17.8   MeV,
    a_3 = 23.6   MeV,
    a_4 =  0.712 MeV.
$$
$Z$は
$$
    Z_0
    = \frac{
        A/2
    }{
        1 + \frac{a_4}{4a_3} A^{2/3}
    }
    = \frac{A/2}{1 + 0.00754A^{2/3}}
$$
を使う.
"""

using Plots
using Plots.PlotMeasures

# Plotの範囲は(質量数)A = 1~250

A = 1:1:250

a₁ = 15.7  # MeV
a₂ = 17.8  # MeV
a₃ = 23.6  # MeV
a₄ = 0.712 # MeV

Z₀(A) = (A/2) / (1 + a₄ / (4*a₃) * A^(2/3))

Z = Z₀.(A)
volume = (a₁.*A) ./ A
surface = a₂ .* A.^(2/3) ./ A
sym = a₃ .* (A - 2 .* Z).^2 ./ A ./ A
coulomb = a₄ .* Z.^2 ./ (A.^(1/3)) ./ A

labelfont = Plots.font(20, "Helvetica")
tickfont = Plots.font(18, "Helvetica")
legendfont = Plots.font(12, "Helvetica")

plt=plot(
    A,
    volume,
    label="Volume term",
    linewidth=2,
    xlims=(0, 250),
    ylims=(0,16),
    xticks=0:50:250,
    yticks=0:2:16,
    legend = :bottomright,
    framestyle = :zerolines,
    size=(800,600),
    titlefont=labelfont,
    guidefont=labelfont,
    tickfont=tickfont,
    legendfont=legendfont,
    top_margin=5mm,
    margin=10mm,
    grid=true,
    gridstyle=:dash,
    gridalpha=0.8,
)
plot!(
    A,
    volume .- surface,
    fillrange=volume,
    fillalpha=0.5,
    label="Surface term",
    linewidth=2,
)
plot!(
    A,
    volume .- surface .- sym,
    fillrange=volume .- surface,
    fillalpha=0.5,
    label="Symmetry term",
    linewidth=2,
)
plot!(
    A,
    volume .- surface .- sym .- coulomb,
    fillrange=volume .- surface .- sym,
    fillalpha=0.5,
    label="Coulomb term",
    linewidth=2,
)
title!("Weizsäcker-Bethe Mass Formula")
xlabel!("A")
ylabel!("B/A (MeV)")
savefig(plt,"Weiszacker-Bethe.pdf")