Reverse engineering units (almost working)

This commit is contained in:
Nuwan Yapa 2025-01-30 18:10:41 -05:00
parent 220dd0eae3
commit 174dc40f20
1 changed files with 19 additions and 19 deletions

View File

@ -3,19 +3,19 @@ using DifferentialEquations
const ħc = 197.327 # MeVfm const ħc = 197.327 # MeVfm
# Values defined in C. J. Horowitz and J. Piekarewicz, Phys. Rev. Lett. 86, 5647 (2001) # Values defined in C. J. Horowitz and J. Piekarewicz, Phys. Rev. Lett. 86, 5647 (2001)
# Values taken from Hartree.f (FSUGarnet) and # Values taken from Hartree.f (FSUGarnet)
const m_ρ = 763.0 # MeV/c2
const m_ω = 782.5 # MeV/c2
const m_s = 496.939473213388 # MeV/c2 const m_s = 496.939473213388 # MeV/c2
const m_γ = 0.000001000 # MeV/c2 # not defined in paper const m_ω = 782.5 # MeV/c2
const g2_s = 110.349189097820 # units? const m_ρ = 763.0 # MeV/c2
const g2_v = 187.694676506801 # units? const m_γ = 0.000001000 # MeV/c2 # should be 0?
const g2_ρ = 192.927428365698 # units? const g2_s = 110.349189097820 # dimensionless
const g2_v = 187.694676506801 # 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 # units? const κ = 3.260178893447
const λ = -0.003551486718 # units? # LambdaSS const λ = -0.003551486718 # LambdaSS
const ζ = 0.023499504053 # units? # LambdaVV const ζ = 0.023499504053 # LambdaVV
const Λv = 0.043376933644 # units? # LambdaVR const Λv = 0.043376933644 # LambdaVR
const r_reg = 1E-9 # fm # regulator for Green's functions const r_reg = 1E-9 # fm # regulator for Green's functions
@ -47,21 +47,21 @@ end
r is the radius in fm. r is the radius in fm.
Reference: P. Giuliani, K. Godbey, E. Bonilla, F. Viens, and J. Piekarewicz, Frontiers in Physics 10, (2023)" Reference: P. Giuliani, K. Godbey, E. Bonilla, F. Viens, and J. Piekarewicz, Frontiers in Physics 10, (2023)"
function solveMesonWfs(ρ_sp, ρ_vp, ρ_sn, ρ_vn, r_max, divs, iterations=3) function solveMesonWfs(ρ_sp, ρ_vp, ρ_sn, ρ_vn, r_max, divs, iterations=3)
# A0 doesn't need iterations (Φ0, W0, B0, A0) = (zeros(1 + divs) for _ in 1:4)
src_A0 = -g2_γ .* ρ_vp .* ħc (src_Φ0, src_W0, src_B0, src_A0) = (zeros(1 + divs) for _ in 1:4)
A0 = solveKG(0, src_A0, r_max)
(Φ0, W0, B0) = (zeros(1 + divs) for _ in 1:3) # A0 doesn't need iterations
(src_Φ0, src_W0, src_B0) = (zeros(1 + divs) for _ in 1:3) @. src_A0 = -g2_γ * ρ_vp * ħc
A0 .= solveKG(m_γ, src_A0, r_max)
for _ in 1:iterations for _ in 1:iterations
@. src_Φ0 = g2_s * ((κ/2) * Φ0^2 + (λ/6) * Φ0^3 - (ρ_sp + ρ_sn)) * ħc @. src_Φ0 = g2_s * ((κ/2) * (Φ0/ħc)^2 + (λ/6) * (Φ0/ħc)^3) - g2_s * (ρ_sp + ρ_sn) * ħc
Φ0 .= solveKG(m_s, src_Φ0, r_max) Φ0 .= solveKG(m_s, src_Φ0, r_max)
@. src_W0 = g2_v * ((ζ/6) * W0^3 + 2 * Λv * B0^2 * W0 - (ρ_vp + ρ_vn)) * ħc @. src_W0 = g2_v * ((ζ/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 = g2_ρ * (2 * Λv * W0^2 * B0 - (ρ_vp - ρ_vn) / 2) * ħc @. src_B0 = g2_ρ * (2 * Λv * (W0/ħc)^2 * (B0/ħc)) - (g2_ρ / 2) * (ρ_vp - ρ_vn) * ħc
B0 .= solveKG(m_ρ, src_B0, r_max) B0 .= solveKG(m_ρ, src_B0, r_max)
end end