diff --git a/dirac.jl b/dirac.jl index 731a864..a2c0437 100644 --- a/dirac.jl +++ b/dirac.jl @@ -27,7 +27,7 @@ end divs is the number of mesh divisions so solution would be returned as a 2×(1+divs) matrix, shooting method divides the interval into two partitions at r_max/2, ensuring convergence at both r=0 and r=r_max, the other parameters are the same from dirac!(...)." -function solveWf(κ, p, E, Φ0, W0, B0, A0, r_max, divs; shooting=true, normalize=true) +function solveDiracWf(κ, p, E, Φ0, W0, B0, A0, r_max, divs; shooting=true, normalize=true) saveat = r_max / divs if shooting @@ -124,7 +124,7 @@ function calculateNucleonDensity(N, p, Φ0, W0, B0, A0, r_max, divs, E_min=0, E_ ρr2 = zeros(2, divs + 1) # ρ×r² values for (κ, E, occ) in zip(κs, Es, occs) - wf = solveWf(κ, p, E, Φ0, W0, B0, A0, r_max, divs; shooting=true, normalize=true) + wf = solveDiracWf(κ, p, E, Φ0, W0, B0, A0, r_max, divs; shooting=true, normalize=true) wf2 = wf .* wf ρr2 += (occ / (4 * pi)) * wf2 # 2j+1 factor is accounted in the occupancy number end diff --git a/mesons.jl b/mesons.jl index 675a34a..e2c2312 100644 --- a/mesons.jl +++ b/mesons.jl @@ -17,36 +17,47 @@ const λ = -0.003551486718 # LambdaSS const ζ = 0.023499504053 # LambdaVV const Λv = 0.043376933644 # LambdaVR -const r_reg = 1E-6 # regulator for the centrifugal term in fm +const r_reg = 1E-9 # regulator for Green's functions in fm -"Coupled radial meson equations that returns ddu=[ddΦ0, ddW0, ddB0, ddA0] in-place where - du=[dΦ0, dW0, dB0, dA0] are the first derivatives of meson fields evaluated at r, - u=[Φ0, W0, B0, A0] are the meson fields evaluated at r, - κ is the generalized angular momentum, - ρ_sp, ρ_vp are the scalar and vector proton densities as funtions of r in fm, - ρ_sn, ρ_vn are the scalar and vector neutron densities as funtions of r in fm, - Φ0, W0, B0, A0 are the mean-field potentials (couplings included) in MeV as functions of r in fm, +"Green's function for Klein-Gordon equation" +greensFunctionKG(m, r, rp) = -(1 / m) * rp / (r + r_reg) * sinh(m * min(r, rp)) * exp(-m * max(r, rp)) + +"Green's function for Poisson's equation" +greensFunctionP(r, rp) = -rp^2 / (max(r, rp) + r_reg) + +"Solve the Klein-Gordon equation (or Poisson's equation when m=0) for a given a source function (as an array) where + the other parameters are the same from mesons!(...)." +function solveKG(m, source, r_max) + greensFunction = m == 0 ? greensFunctionP : (r, rp) -> greensFunctionKG(m, r, rp) + mesh = range(0, r_max; length=length(source)) + greenMat = greensFunction.(transpose(mesh), mesh) + return greenMat * source +end + +"Iteratively solve meson equations and return the wave functions u(r)=[Φ0(r), W0(r), B0(r), A0(r)] where + divs is the number of mesh divisions so the arrays are of length (1+divs), + ρ_sp, ρ_vp are the scalar and vector proton densities as arrays, + ρ_sn, ρ_vn are the scalar and vector neutron densities as arrays, + Φ0, W0, B0, A0 are the mean-field potentials (couplings included) in MeV as as arrays, r is the radius in fm. Reference: P. Giuliani, K. Godbey, E. Bonilla, F. Viens, and J. Piekarewicz, Frontiers in Physics 10, (2023)" -function mesons!(ddu, du, u, (ρ_sp, ρ_vp, ρ_sn, ρ_vn), r) - (dΦ0, dW0, dB0, dA0) = du - (Φ0, W0, B0, A0) = u +function solveMesonWfs(ρ_sp, ρ_vp, ρ_sn, ρ_vn, r_max, divs, iterations=3) + (Φ0, W0, B0, A0) = (zeros(1 + divs) for _ in 1:4) + (src_Φ0, src_W0, src_B0, src_A0) = (zeros(1 + divs) for _ in 1:4) - ddΦ0 = -(2/(r + r_reg)) * dΦ0 + m_s^2 * Φ0 / ħc + g2_s * ((κ/2) * Φ0^2 + (λ/6) * Φ0^3 - (ρ_sp(r) + ρ_sn(r))) / ħc - ddW0 = -(2/(r + r_reg)) * dW0 + m_ω^2 * W0 / ħc + g2_v * ((ζ/6) * W0^3 + 2 * Λv * B0^2 * W0 - (ρ_vp(r) + ρ_vn(r))) / ħc - ddB0 = -(2/(r + r_reg)) * dB0 + m_ρ^2 * B0 / ħc + g2_ρ * (2 * Λv * W0^2 * B0 - (ρ_vp(r) - ρ_vn(r)) / 2) / ħc - ddA0 = -(2/(r + r_reg)) * dA0 - ρ_vp(r) / ħc # e goes here? + for _ in iterations + @. src_Φ0 = g2_s * ((κ/2) * Φ0^2 + (λ/6) * Φ0^3 - (ρ_sp + ρ_sn)) / ħc + Φ0 .= solveKG(m_s / sqrt(ħc), src_Φ0, r_max) - ddu .= [ddΦ0, ddW0, ddB0, ddA0] -end - -"Solve meson equations and return the wave functions u(r)=[Φ0(r), W0(r), B0(r), A0(r)] where - divs is the number of mesh divisions so solution would be returned as a 4×(1+divs) matrix, - the other parameters are the same from mesons!(...)." -function solveWfs(ρ_sp, ρ_vp, ρ_sn, ρ_vn, r_max, divs) - prob = SecondOrderODEProblem(mesons!, BigFloat.([0, 0, 0, 0]), BigFloat.([0, 0, 0, 0]), (r_max, 0)) - sol = solve(prob, Feagin14(), p=(ρ_sp, ρ_vp, ρ_sn, ρ_vn), saveat=r_max/divs) - wfs = hcat(sol.u...) - - return wfs + @. src_W0 = g2_v * ((ζ/6) * W0^3 + 2 * Λv * B0^2 * W0 - (ρ_vp + ρ_vn)) / ħc + W0 .= solveKG(m_ω / sqrt(ħc), src_W0, r_max) + + @. src_B0 = g2_ρ * (2 * Λv * W0^2 * B0 - (ρ_vp - ρ_vn) / 2) / ħc + W0 .= solveKG(m_ρ / sqrt(ħc), src_B0, r_max) + + @. src_A0 = ρ_vp / ħc + A0 .= solveKG(0, src_A0, r_max) # this doesn't need iterations + end + + return (Φ0, W0, B0, A0) end diff --git a/test/Pb208_wf.jl b/test/Pb208_dens.jl similarity index 90% rename from test/Pb208_wf.jl rename to test/Pb208_dens.jl index 1a8a5ee..46da8ab 100644 --- a/test/Pb208_wf.jl +++ b/test/Pb208_dens.jl @@ -25,7 +25,7 @@ groundE = findEs(κ, p, S_interp, V_interp, R_interp, A_interp, r_max, E_min, E_ println("ground state E = $groundE") divs = 400 -wf = solveWf(κ, p, groundE, S_interp, V_interp, R_interp, A_interp, r_max, divs) +wf = solveDiracWf(κ, p, groundE, S_interp, V_interp, R_interp, A_interp, r_max, divs) rs = range(0, r_max, length=divs+1) gs = wf[1, :] fs = wf[2, :] diff --git a/test/Pb208_flds.jl b/test/Pb208_flds.jl new file mode 100644 index 0000000..caa026b --- /dev/null +++ b/test/Pb208_flds.jl @@ -0,0 +1,19 @@ +using DelimitedFiles, Interpolations, Plots +include("../mesons.jl") + +# test data generated from Hartree.f +# format: x Rhos(n) Rhov(n) Rhot(n) Rhos(p) Rhov(p) Rhot(p) +test_data = readdlm("test/Pb208DensFSUGarnet.csv") +xs = test_data[:, 1] +ρ_sn = test_data[:, 2] +ρ_vn = test_data[:, 3] +ρ_sp = test_data[:, 5] +ρ_vp = test_data[:, 6] + +r_max = maximum(xs) +divs = length(xs) - 1 + +(Φ0, W0, B0, A0) = solveMesonWfs(ρ_sp, ρ_vp, ρ_sn, ρ_vn, r_max, divs, 3) + +plot(xs, hcat(Φ0, W0, B0, A0), layout=4, label=["Φ0" "W0" "B0" "A0"]) +xlabel!("r (fm)") diff --git a/test/Pb208_wfs.jl b/test/Pb208_wfs.jl deleted file mode 100644 index 58fbfb9..0000000 --- a/test/Pb208_wfs.jl +++ /dev/null @@ -1,26 +0,0 @@ -using DelimitedFiles, Interpolations, Plots -include("../mesons.jl") - -# test data generated from Hartree.f -# format: x Rhos(n) Rhov(n) Rhot(n) Rhos(p) Rhov(p) Rhot(p) -test_data = readdlm("test/Pb208DensFSUGarnet.csv") -xs = test_data[:, 1] -rho_sn = test_data[:, 2] -rho_vn = test_data[:, 3] -rho_sp = test_data[:, 5] -rho_vp = test_data[:, 6] - -ρ_sn = linear_interpolation(xs, rho_sn) -ρ_vn = linear_interpolation(xs, rho_vn) -ρ_sp = linear_interpolation(xs, rho_sp) -ρ_vp = linear_interpolation(xs, rho_vp) - -r_max = maximum(xs) -divs = 400 - -wfs = solveWfs(ρ_sp, ρ_vp, ρ_sn, ρ_vn, r_max, divs) - -rs = range(0, r_max, length=divs+1) - -plot(rs, transpose(wfs)) -xlabel!("r (fm)")