# Build your own nuclear DFT -- step 1: the grid and a guess density.
#
# Everything lives on a uniform radial mesh r_i = (i+1)h that EXCLUDES the
# origin: the 1/r terms coming later would blow up at r = 0, and the boundary
# condition u(0) = 0 handles the origin for us. The only other tool we need is
# a volume integral: for a spherically symmetric f(r),
#     integral f d^3r  =  4 pi  integral f(r) r^2 dr      (trapezoid rule).
#
# With those two pieces we can already write down a nucleus-shaped guess: the
# Woods-Saxon profile 1/(1 + exp((r-R0)/a)) -- flat interior, soft surface --
# normalized so it holds exactly A = 16 nucleons.

pi = 3.141592653589793
h = 0.1              # grid spacing [fm]
nr = 140             # 14 fm box
amass = 16.0

r = zeros(nr)
for i to nr { r[i] = (i + 1.0) * h }

def integrate_r2(f, n, hh) {
    s = 0.5 * f[0] * r[0] * r[0] + 0.5 * f[n - 1.0] * r[n - 1.0] * r[n - 1.0]
    for i to n {
        if i > 0 {
            if i < n - 1.0 { s += f[i] * r[i] * r[i] }
        }
    }
    return 4.0 * pi * s * hh
}

# Woods-Saxon shape: R0 = 1.2 A^(1/3) fm is the nuclear radius, a = 0.6 fm the
# surface thickness. Normalize it to A particles.
r0ws = 1.2 * pow(amass, 1.0 / 3.0)
rho = zeros(nr)
for i to nr { rho[i] = 1.0 / (1.0 + exp((r[i] - r0ws) / 0.6)) }
norm = integrate_r2(rho, nr, h)
for i to nr { rho[i] = rho[i] * amass / norm }

# check: integrating the density must count the nucleons
a_check = integrate_r2(rho, nr, h)
show a_check         # 16 -- the normalization worked

# the central density, in nucleons/fm^3 -- compare nuclear saturation ~0.16
rho_central = rho[0]
show rho_central

show rho             # the guess density: flat middle, soft edge