yapanuwan
/
NuclearRMF
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

All units sorted out (works but unstable)

This commit is contained in:
Nuwan Yapa 2025-01-31 13:47:45 -05:00
parent 8a42579594
commit 20e6b5ba84
1 changed files with 10 additions and 12 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

22
mesons.jl
View File

@ -12,10 +12,10 @@ const g2_s = 110.349189097820 # dimensionless
const g2_v = 187.694676506801 # dimensionless const g2_v = 187.694676506801 # dimensionless
const g2_ρ = 192.927428365698 # dimensionless const g2_ρ = 192.927428365698 # dimensionless
const g2_γ = 0.091701236 # dimensionless # equal to 4πα const g2_γ = 0.091701236 # dimensionless # equal to 4πα
const κ = 3.260178893447 const κ = 3.260178893447 # MeV
const λ = -0.003551486718 # LambdaSS const λ = -0.003551486718 # dimensionless # LambdaSS
const ζ = 0.023499504053 # LambdaVV const ζ = 0.023499504053 # dimensionless # LambdaVV
const Λv = 0.043376933644 # LambdaVR const Λv = 0.043376933644 # dimensionless # LambdaVR
const r_reg = 1E-9 # fm # regulator for Green's functions const r_reg = 1E-9 # fm # regulator for Green's functions
@ -25,13 +25,13 @@ greensFunctionKG(m, r, rp) = -1 / (m * (r + r_reg) * (rp + r_reg)) * sinh(m * mi
"Green's function for Poisson's equation in natural units" "Green's function for Poisson's equation in natural units"
greensFunctionP(r, rp) = -1 / (max(r, rp) + r_reg) greensFunctionP(r, rp) = -1 / (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 "Solve the Klein-Gordon equation (or Poisson's equation when m=0) and return an array in MeV for a source array given in fm⁻³ where
m is the mass of the meson in MeV/c2, m is the mass of the meson in MeV/c2,
r_max is the r-cutoff in fm." r_max is the r-cutoff in fm."
function solveKG(m, source, r_max) function solveKG(m, source, r_max)
Δr = r_max / (length(source) - 1) Δr = r_max / (length(source) - 1)
rs = range(0, r_max; length=length(source)) rs = range(0, r_max; length=length(source))
int_measure = Δr .* rs .^ 2 int_measure = ħc .* Δr .* rs .^ 2
greensFunction = m == 0 ? greensFunctionP : (r, rp) -> greensFunctionKG(m / ħc, r, rp) greensFunction = m == 0 ? greensFunctionP : (r, rp) -> greensFunctionKG(m / ħc, r, rp)
greenMat = greensFunction.(rs, transpose(rs)) greenMat = greensFunction.(rs, transpose(rs))
@ -52,17 +52,15 @@ function solveMesonWfs(ρ_sp, ρ_vp, ρ_sn, ρ_vn, r_max, divs, iterations=3; in
(src_Φ0, src_W0, src_B0, src_A0) = (zeros(1 + divs) for _ in 1:4) (src_Φ0, src_W0, src_B0, src_A0) = (zeros(1 + divs) for _ in 1:4)
# A0 doesn't need iterations # A0 doesn't need iterations
@. src_A0 = -g2_γ * ρ_vp * ħc @. src_A0 = -g2_γ * ρ_vp
A0 .= solveKG(m_γ, src_A0, r_max) A0 .= solveKG(m_γ, src_A0, r_max)
for _ in 1:iterations for _ in 1:iterations
@. src_Φ0 = g2_s * ((κ/2) * (Φ0/ħc)^2 + (λ/6) * (Φ0/ħc)^3) - g2_s * ħc * (ρ_sp + ρ_sn) @. src_Φ0 = g2_s * ((κ/ħc)/2 * (Φ0/ħc)^2 + (λ/6) * (Φ0/ħc)^3) - g2_s * (ρ_sp + ρ_sn)
@. src_W0 = g2_v * ((ζ/6) * (W0/ħc)^3 + 2Λv * (B0/ħc)^2 * (W0/ħc)) - g2_v * (ρ_vp + ρ_vn)
@. src_B0 = 2Λv * g2_ρ * (W0/ħc)^2 * (B0/ħc) - g2_ρ/4 * (ρ_vp - ρ_vn)
Φ0 .= solveKG(m_s, src_Φ0, r_max) Φ0 .= solveKG(m_s, src_Φ0, r_max)
@. src_W0 = g2_v * ħc * ((ζ/6) * (W0/ħc)^3 + 2Λv * (B0/ħc)^2 * (W0/ħc)) - g2_v * (ρ_vp + ρ_vn) * ħc
W0 .= solveKG(m_ω, src_W0, r_max) W0 .= solveKG(m_ω, src_W0, r_max)
@. src_B0 = 2Λv * g2_ρ * W0^2 * B0 / ħc^2 - (g2_ρ / 4) * (ρ_vp - ρ_vn) * ħc
B0 .= solveKG(m_ρ, src_B0, r_max) B0 .= solveKG(m_ρ, src_B0, r_max)
end end