using Plots
using ClassicalDFT
L = 10
dz = 1//100
geom = PlanarGeometry(L, dz; boundary=:walls)
functional = RosenfeldFMTFunctional(geom)
Vext(z) = 0
T, μ = 1.0, 3.0
ρ = solve(functional, Vext, T; μ)
plot(coordinates(geom), ρ, label="ρ(z)")

using Plots
using ClassicalDFT
R = 10
dr = 1//100
Rsolute = 2
geom = SphericalGeometry(R, dr, Rsolute)
functional = RosenfeldFMTFunctional(geom)
Vext(r) = 0
T, μ = 1.0, 3.0
ρ = solve(functional, Vext, T; μ)
plot(coordinates(geom), ρ, label="ρ(r)")

using Plots
using ClassicalDFT
L = 10
dz = 1//100
geom = PlanarGeometry(L, dz; boundary=:walls)
functional = CompositeFunctional(geom, RosenfeldFMTFunctional, LennardJonesMeanFieldFunctional)
Vext(z) = 0
T, μ = 1.5, 1.0
ρ = solve(functional, Vext, T; μ)
plot(coordinates(geom), ρ, label="ρ(z)")

using Plots
using ClassicalDFT
L = 50
dx = 1 // 100
geom = PlanarGeometry(L, dx; boundary=:periodic)
Vext(x) = 0
ρ_init(x) = x < L/4 || x > 3*L/4 ? 1 : 0
ρb = 0.3
T = 0.9
functional = CompositeFunctional(geom, RosenfeldFMTFunctional, LennardJonesMeanFieldFunctional)
solver = PicardIterationSolver(α=0.02; adaptive=true, αmin=0.000001)
callbacks = [PlotCallback()]
ρ = solve(functional, Vext, T; ρb, ρ_init)
plot(coordinates(geom), ρ, label="ρ(z)")
