From 50ff7ab08e691a88c94078161e5aafac094bedb2 Mon Sep 17 00:00:00 2001
From: Ryan Anderson <rymanderson@gmail.com>
Date: Thu, 30 Jan 2025 11:20:20 -0700
Subject: [PATCH] check vpm regularization

---
 scripts/vpm_regularization.jl | 43 +++++++++++++++++++++++++++++++++++
 1 file changed, 43 insertions(+)
 create mode 100644 scripts/vpm_regularization.jl

diff --git a/scripts/vpm_regularization.jl b/scripts/vpm_regularization.jl
new file mode 100644
index 0000000..50f76ba
--- /dev/null
+++ b/scripts/vpm_regularization.jl
@@ -0,0 +1,43 @@
+using Roots
+using SpecialFunctions: erf
+using PythonPlot
+using LaTeXStrings
+
+function upper_bound(σ, ω, ε)
+    return ω / (8 * pi * ε * σ) * (sqrt(2/pi) + sqrt(2/(pi*σ*σ) + 16 * pi * ε / ω))
+end
+
+function residual(ρ_σ, σ, ω, ε)
+    t1 = 4*pi*σ*σ*ε*ρ_σ*ρ_σ / ω
+    t2 = erf(ρ_σ / sqrt(2))
+    t3 = sqrt(2/pi) * ρ_σ * exp(-ρ_σ*ρ_σ*0.5)
+    return t1 + t2 - t3 - 1.0
+end
+
+function solve_ρ_over_σ(σ, ω, ε)
+    return Roots.find_zero((x) -> residual(x, σ, ω, ε), (0.0, upper_bound(σ, ω, ε)), Roots.Brent())
+end
+
+ω = 1.0
+σs = [10.0^(-i) for i in 0:5]
+εs = [10.0^(-i) for i in range(0,stop=9,length=50)]
+
+res = zeros(length(εs), length(σs))
+
+for (iσ,σ) in enumerate(σs)
+    for (iε,ε) in enumerate(εs)
+        res[iε,iσ] = solve_ρ_over_σ(σ, ω, ε)
+    end
+end
+
+fig = figure("regularization")
+fig.clear()
+fig.add_subplot(111, xlabel=L"\varepsilon_{tol}", ylabel=L"\rho / \sigma")
+ax = fig.get_axes()[0]
+ncolors = length(σs)
+colormap = get_cmap("RdBu",ncolors)
+for iσ in 1:length(σs)
+    ax.plot(εs, res[:,iσ], label="σ = $(round(σs[iσ], sigdigits=1))", color=colormap(iσ-1))
+end
+ax.legend()
+ax.set_xscale("log")