From f1553f548f9379710ce0d278fdfeb5bd298224de Mon Sep 17 00:00:00 2001 From: Nuwan Yapa Date: Fri, 21 Feb 2025 10:55:11 -0500 Subject: [PATCH] Reorganize dependency structure --- NuclearRMF.jl | 56 +++++++++++++++++++++++++++++++ common.jl | 2 ++ mesons.jl | 7 ++-- nucleons.jl | 5 ++- system.jl | 58 +-------------------------------- test/Linear_nucleon_spectrum.jl | 4 +-- test/Pb208.jl | 2 +- test/Pb208_meson_flds.jl | 4 +-- test/Pb208_nucleon_dens.jl | 2 +- test/Pb208_nucleon_spectrum.jl | 4 +-- test/Pb208_nucleon_wf.jl | 4 +-- test/Pb208_shooting.jl | 2 +- 12 files changed, 74 insertions(+), 76 deletions(-) create mode 100644 NuclearRMF.jl create mode 100644 common.jl diff --git a/NuclearRMF.jl b/NuclearRMF.jl new file mode 100644 index 0000000..cdf314a --- /dev/null +++ b/NuclearRMF.jl @@ -0,0 +1,56 @@ +using PolyLog, Plots +include("nucleons.jl") +include("mesons.jl") + +"Total binding energy of the system" +total_E(s::system, Es_p, occs_p, Es_n, occs_n) = sum(occs_p .* (M_p .- Es_p)) + sum(occs_n .* (M_n .- Es_n)) - meson_E(s) + +"Normalized Woods-Saxon form used for constructing an initial solution" +Woods_Saxon(r::Float64; R::Float64=7.0, a::Float64=0.5) = -1 / (8π * a^3 * reli3(-exp(R / a)) * (1 + exp((r - R) / a))) + +"Run the full program by self-consistent solution of nucleon and meson densities" +function solve_system(s::system; reinitialize_densities=true, monitor_print=true, monitor_plot=false) + if reinitialize_densities + dens_guess = Woods_Saxon.(rs(s)) + @. s.ρ_sp = s.Z * dens_guess + @. s.ρ_vp = s.Z * dens_guess + @. s.ρ_sn = s.N * dens_guess + @. s.ρ_vn = s.N * dens_guess + end + + if monitor_plot + p = plot(legends=false, size=(1024, 768), layout=(2, 4), title=["ρₛₚ" "ρᵥₚ" "ρₛₙ" "ρᵥₙ" "Φ₀" "W₀" "B₀" "A₀"]) + end + + previous_E_per_A = NaN + + while true + # mesons + @time "Meson fields" solveMesonFields!(s, isnan(previous_E_per_A) ? 50 : 15) + + # protons + @time "Proton spectrum" (κs_p, Es_p) = findAllOrbitals(true, s) + occs_p = fillNucleons(s.Z, κs_p, Es_p) + @time "Proton densities" (s.ρ_sp, s.ρ_vp) = calculateNucleonDensity(κs_p, Es_p, occs_p, true, s) + + # neutrons + @time "Neutron spectrum" (κs_n, Es_n) = findAllOrbitals(false, s) + occs_n = fillNucleons(s.N, κs_n, Es_n) + @time "Neutron densities" (s.ρ_sn, s.ρ_vn) = calculateNucleonDensity(κs_n, Es_n, occs_n, false, s) + + if monitor_plot + for s in p.series_list + s.plotattributes[:linecolor] = :gray + end + plot!(p, rs(s), hcat(s.ρ_sp, s.ρ_vp, s.ρ_sn, s.ρ_vn, s.Φ0, s.W0, s.B0, s.A0), linecolor=:red) + display(p) + end + + E_per_A = total_E(s, Es_p, occs_p, Es_n, occs_n) / A(s) + monitor_print && println("Total binding E per nucleon = $E_per_A") + + # check convergence + abs(previous_E_per_A - E_per_A) < 0.0001 && break + previous_E_per_A = E_per_A + end +end diff --git a/common.jl b/common.jl new file mode 100644 index 0000000..ccb93b5 --- /dev/null +++ b/common.jl @@ -0,0 +1,2 @@ +const ħc = 197.33 # MeVfm +const r_reg = 1E-8 # fm # regulator for R diff --git a/mesons.jl b/mesons.jl index 2eb9ca0..1481acd 100644 --- a/mesons.jl +++ b/mesons.jl @@ -1,6 +1,5 @@ -using DifferentialEquations - -const ħc = 197.33 # MeVfm +include("common.jl") +include("system.jl") # Values defined in C. J. Horowitz and J. Piekarewicz, Phys. Rev. Lett. 86, 5647 (2001) # Values taken from Hartree.f (FSUGarnet) @@ -17,8 +16,6 @@ const λ = -0.003551486718 # dimensionless # LambdaSS const ζ = 0.023499504053 # dimensionless # LambdaVV const Λv = 0.043376933644 # dimensionless # LambdaVR -const r_reg = 1E-8 # fm # regulator for Green's functions - "Green's function for Klein-Gordon equation in natural units" greensFunctionKG(m, r, rp) = -1 / (m * (r + r_reg) * (rp + r_reg)) * sinh(m * min(r, rp)) * exp(-m * max(r, rp)) diff --git a/nucleons.jl b/nucleons.jl index 35f61df..98f3c0b 100644 --- a/nucleons.jl +++ b/nucleons.jl @@ -1,12 +1,11 @@ using DifferentialEquations include("bisection.jl") +include("common.jl") +include("system.jl") -const ħc = 197.33 # MeVfm const M_n = 939.0 # MeV/c2 const M_p = 939.0 # MeV/c2 -const r_reg = 1E-8 # fm # regulator for the centrifugal term - "The spherical Dirac equation that returns du=[dg, df] in-place where u=[g, f] are the reduced radial components evaluated at r, κ is the generalized angular momentum, diff --git a/system.jl b/system.jl index 1c860ef..531e15a 100644 --- a/system.jl +++ b/system.jl @@ -1,4 +1,4 @@ -using Interpolations, PolyLog, Plots +using Interpolations "Defines a nuclear system containing relevant parameters and meson/nucleon densities" mutable struct system @@ -40,9 +40,6 @@ Z_or_N(s::system, p::Bool)::Int = p ? s.Z : s.N "Create an empty array for the size of the mesh" zero_array(s::system) = zeros(1 + s.divs) -"Normalized Woods-Saxon form used for constructing an initial solution" -Woods_Saxon(r::Float64; R::Float64=7.0, a::Float64=0.5) = -1 / (8π * a^3 * reli3(-exp(R / a)) * (1 + exp((r - R) / a))) - "Create linear interpolations for the meson fields" function fields_interp(s::system) S_interp = linear_interpolation(rs(s), s.Φ0) @@ -51,56 +48,3 @@ function fields_interp(s::system) A_interp = linear_interpolation(rs(s), s.A0) return (S_interp, V_interp, R_interp, A_interp) end - -include("nucleons.jl") -include("mesons.jl") - -"Run the full program by self-consistent solution of nucleon and meson densities" -function solve_system(s::system; reinitialize_densities=true, monitor_print=true, monitor_plot=false) - if reinitialize_densities - dens_guess = Woods_Saxon.(rs(s)) - @. s.ρ_sp = s.Z * dens_guess - @. s.ρ_vp = s.Z * dens_guess - @. s.ρ_sn = s.N * dens_guess - @. s.ρ_vn = s.N * dens_guess - end - - if monitor_plot - p = plot(legends=false, size=(1024, 768), layout=(2, 4), title=["ρₛₚ" "ρᵥₚ" "ρₛₙ" "ρᵥₙ" "Φ₀" "W₀" "B₀" "A₀"]) - end - - previous_E_per_A = NaN - - while true - # mesons - @time "Meson fields" solveMesonFields!(s, isnan(previous_E_per_A) ? 50 : 15) - - # protons - @time "Proton spectrum" (κs_p, Es_p) = findAllOrbitals(true, s) - occs_p = fillNucleons(s.Z, κs_p, Es_p) - @time "Proton densities" (s.ρ_sp, s.ρ_vp) = calculateNucleonDensity(κs_p, Es_p, occs_p, true, s) - - # neutrons - @time "Neutron spectrum" (κs_n, Es_n) = findAllOrbitals(false, s) - occs_n = fillNucleons(s.N, κs_n, Es_n) - @time "Neutron densities" (s.ρ_sn, s.ρ_vn) = calculateNucleonDensity(κs_n, Es_n, occs_n, false, s) - - if monitor_plot - for s in p.series_list - s.plotattributes[:linecolor] = :gray - end - plot!(p, rs(s), hcat(s.ρ_sp, s.ρ_vp, s.ρ_sn, s.ρ_vn, s.Φ0, s.W0, s.B0, s.A0), linecolor=:red) - display(p) - end - - E_per_A = total_E(s, Es_p, occs_p, Es_n, occs_n) / A(s) - monitor_print && println("Total binding E per nucleon = $E_per_A") - - # check convergence - abs(previous_E_per_A - E_per_A) < 0.0001 && break - previous_E_per_A = E_per_A - end -end - -"Total binding energy of the system" -total_E(s::system, Es_p, occs_p, Es_n, occs_n) = sum(occs_p .* (M_p .- Es_p)) + sum(occs_n .* (M_n .- Es_n)) - meson_E(s) diff --git a/test/Linear_nucleon_spectrum.jl b/test/Linear_nucleon_spectrum.jl index e8084a4..12cb119 100644 --- a/test/Linear_nucleon_spectrum.jl +++ b/test/Linear_nucleon_spectrum.jl @@ -1,5 +1,5 @@ -using DelimitedFiles, Interpolations, Plots -include("../system.jl") +using DelimitedFiles, Plots +include("../nucleons.jl") # test data generated from Hartree.f # format: x S(x) V(x) R(x) A(x) diff --git a/test/Pb208.jl b/test/Pb208.jl index e0b5d97..e493195 100644 --- a/test/Pb208.jl +++ b/test/Pb208.jl @@ -1,4 +1,4 @@ -include("../system.jl") +include("../NuclearRMF.jl") s = system(82, 126, 20.0, 400) diff --git a/test/Pb208_meson_flds.jl b/test/Pb208_meson_flds.jl index c438106..ccab61b 100644 --- a/test/Pb208_meson_flds.jl +++ b/test/Pb208_meson_flds.jl @@ -1,5 +1,5 @@ -using DelimitedFiles, Interpolations, Plots -include("../system.jl") +using DelimitedFiles, Plots +include("../mesons.jl") # test data generated from Hartree.f # format: x S(x) V(x) R(x) A(x) diff --git a/test/Pb208_nucleon_dens.jl b/test/Pb208_nucleon_dens.jl index 63575fe..d6887fc 100644 --- a/test/Pb208_nucleon_dens.jl +++ b/test/Pb208_nucleon_dens.jl @@ -1,5 +1,5 @@ using DelimitedFiles, Plots -include("../system.jl") +include("../nucleons.jl") # test data generated from Hartree.f # format: x S(x) V(x) R(x) A(x) diff --git a/test/Pb208_nucleon_spectrum.jl b/test/Pb208_nucleon_spectrum.jl index 0704593..538ec19 100644 --- a/test/Pb208_nucleon_spectrum.jl +++ b/test/Pb208_nucleon_spectrum.jl @@ -1,5 +1,5 @@ -using DelimitedFiles, Interpolations, Plots -include("../system.jl") +using DelimitedFiles, Plots +include("../nucleons.jl") # test data generated from Hartree.f # format: x S(x) V(x) R(x) A(x) diff --git a/test/Pb208_nucleon_wf.jl b/test/Pb208_nucleon_wf.jl index 0018384..70aa5ae 100644 --- a/test/Pb208_nucleon_wf.jl +++ b/test/Pb208_nucleon_wf.jl @@ -1,5 +1,5 @@ -using DelimitedFiles, Interpolations, Plots -include("../system.jl") +using DelimitedFiles, Plots +include("../nucleons.jl") # test data generated from Hartree.f # format: x S(x) V(x) R(x) A(x) diff --git a/test/Pb208_shooting.jl b/test/Pb208_shooting.jl index 05dc7ed..35d37cd 100644 --- a/test/Pb208_shooting.jl +++ b/test/Pb208_shooting.jl @@ -1,5 +1,5 @@ using DelimitedFiles, Plots -include("../system.jl") +include("../nucleons.jl") # test data generated from Hartree.f # format: x S(x) V(x) R(x) A(x)