From 45cc39289d574316502c4264ce5e6ae7f66db7e1 Mon Sep 17 00:00:00 2001
From: Nuwan Yapa <yapanuwan@gmail.com>
Date: Fri, 14 Jun 2024 17:22:11 -0400
Subject: [PATCH] Barebones Dirac equation

---
 Project.toml  |  2 ++
 dirac.jl      | 25 +++++++++++++++++++++++++
 test/dirac.jl |  3 +++
 3 files changed, 30 insertions(+)
 create mode 100644 Project.toml
 create mode 100644 dirac.jl
 create mode 100644 test/dirac.jl

diff --git a/Project.toml b/Project.toml
new file mode 100644
index 0000000..03f6800
--- /dev/null
+++ b/Project.toml
@@ -0,0 +1,2 @@
+[deps]
+DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa"
diff --git a/dirac.jl b/dirac.jl
new file mode 100644
index 0000000..13a7f46
--- /dev/null
+++ b/dirac.jl
@@ -0,0 +1,25 @@
+using DifferentialEquations
+
+ħc = 197.327 # ħc in MeVfm
+M_n = 939.5654133 # Neutron mass in MeV/c2
+M_p = 938.2720813 # Proton mass in MeV/c2
+
+"The spherical Dirac equation that returns du=[dg, df] in-place where
+(g, f) are the reduced radial components evaluated at r,
+κ is the generalized angular momentum,
+M is the mass in MeV/c2,
+E in the energy in MeV,
+S(r) & V(r) are functions corresponding to scalar and vector potentials in MeV,
+r is the radius in fm.
+Reference: P. Giuliani, K. Godbey, E. Bonilla, F. Viens, and J. Piekarewicz, Frontiers in Physics 10, (2023)."
+function dirac!(du, (g, f), (κ, M, E, S, V), r)
+    du[1] = -(κ/r) * g + (E + M - S(r) - V(r)) * f / ħc
+    du[2] =  (κ/r) * f - (E - M + S(r) - V(r)) * g / ħc
+end
+
+"Solve the Dirac equation and return g(r=r_max)"
+function boundaryValue(κ, M, E, S, V, r_max, r_min=r_max/1000)
+    prob = ODEProblem(dirac!, [0, 1], (r_min, r_max))
+    sol = solve(prob, RK4(), p=(κ, M, E, S, V))
+    return sol(r_max)[1]
+end
\ No newline at end of file
diff --git a/test/dirac.jl b/test/dirac.jl
new file mode 100644
index 0000000..8b35112
--- /dev/null
+++ b/test/dirac.jl
@@ -0,0 +1,3 @@
+include("../dirac.jl")
+
+boundaryValue(-1, M_n, 940, r -> 0, r -> 0, 20)
\ No newline at end of file