Added the short-range correction equations to the FMM pdf documentation.

theory/Multipoles/mesh_summary.tex

 \subsection{Coupling the FMM to a mesh for periodic long-range forces}
 \label{ssec:mesh_summary}
 
+\begin{equation}
+  S(x) = \frac{e^x}{1 + e^x}
+\end{equation}
+
+\begin{align}
+  \varphi_s(r) &= \frac{1}{r}\left[2 - 2S\left(\frac{2r}{r_s}\right)\right] \nonumber\\
+  &= \frac{1}{r}\left[2 - \frac{2e^{\frac{2r}{r_s}}}{1+e^{\frac{2r}{r_s}}}\right] 
+\end{align}
+\begin{align}
+  |\mathbf{f}_s(r)| &= \frac{1}{r^2}\left[\frac{4r}{r_s}S'\left(\frac{2r}{r_s}\right) - 2S\left(\frac{2r}{r_s}\right) + 2\right] \nonumber \\
+  &= \frac{1}{r^2}\left[\frac{4r}{r_s}\frac{e^{\frac{2r}{r_s}}}{(1+e^{\frac{2r}{r_s}})^2} - \frac{2e^{\frac{2r}{r_s}}}{1+e^{\frac{2r}{r_s}}} + 2\right]
+\end{align}
+
+\begin{equation}
+  \tilde\varphi_l(k) = \frac{1}{k^2}\left[\frac{\upi}{2}kr_s\textrm{csch}\left(\frac{\upi}{2}kr_s\right) \right]
+\end{equation}
+
 \begin{figure}
 \includegraphics[width=\columnwidth]{potential_short.pdf}
 \caption{aa}
@@ -9,7 +26,14 @@
 
 
 \begin{figure}
-\includegraphics[width=\columnwidth]{potential_long.pdf}
+\includegraphics[width=\columnwidth]{force_short.pdf}
 \caption{bb}
+\label{fig:fmm:force_short}
+\end{figure}
+
+
+\begin{figure}
+\includegraphics[width=\columnwidth]{potential_long.pdf}
+\caption{cc}
 \label{fig:fmm:potential_long}
 \end{figure}

theory/Multipoles/plot_mesh.py

@@ -67,56 +67,42 @@ phi_newton = 1. / r
 phit_newton = 1. / k**2
 force_newton = 1. / r**2
 
-def smoothstep(x): #S_2(x)
-    ret = 6*x**5 - 15*x**4 + 10*x**3
-    #ret = 3*x**2 - 2*x**3
-    ret[x < 0] = 0.
-    ret[x > 1] = 1.
-    return ret
-
 def my_exp(x):
-    #return 1. + x + (x**2 / 2.) + (x**3 / 6.) + (x**4 / 24.) + (x**5 / 120.)# + (x**6 / 720.)
-    return exp(x)
+    return 1. + x + (x**2 / 2.) + (x**3 / 6.) + (x**4 / 24.) + (x**5 / 120.)# + (x**6 / 720.)
+    #return exp(x)
+
+def csch(x): # hyperbolic cosecant
+    return 1. / sinh(x)
 
 def sigmoid(x):
-    return 1. / (1. + 1./my_exp(x))
-    #return x / sqrt(1. + x**2)
+    return my_exp(x) / (my_exp(x) + 1.)
 
 def d_sigmoid(x):
     return my_exp(x) / ((my_exp(x) + 1)**2)
 
 def swift_corr(x):
-    #return 2. * smoothstep(x/4. + 1./2.) - 1.
-    #return sigmoid(4. * x)
     return 2 * sigmoid( 4 * x ) - 1
 
-def d_swift_corr(x):
-    return 2 * d_sigmoid( 4 * x )
-
-def csch(x): # hyperbolic cosecant
-    return 1. / sinh(x)
-
-figure()
-x = linspace(-4, 4, 100)
-plot(x, special.erf(x), '-', color=colors[0])
-plot(x, swift_corr(x), '-', color=colors[1])
-plot(x, x, '-', color=colors[2])
-ylim(-1.1, 1.1)
-xlim(-4.1, 4.1)
-#plot(x, exp(x), '-', color=colors[0])
-#plot(x, my_exp(x), '-', color=colors[1])
-savefig("temp.pdf")
+#figure()
+#x = linspace(-4, 4, 100)
+#plot(x, special.erf(x), '-', color=colors[0])
+#plot(x, swift_corr(x), '-', color=colors[1])
+#plot(x, x, '-', color=colors[2])
+#ylim(-1.1, 1.1)
+#xlim(-4.1, 4.1)
+#savefig("temp.pdf")
 
 # Correction in real space
 corr_short_gadget2 = special.erf(r / (2.*r_s))
-corr_long_gadget2 = exp(-k**2*r_s**2)
 corr_short_swift = swift_corr(r / (2.*r_s)) 
-#corr_long_swift = k * r_s * math.pi / (2. * sinh(0.5 * math.pi * r_s * k))
-corr_long_swift = math.pi * k * r_s * csch(0.5 * math.pi * r_s * k) / 2.
 eta_short_gadget2 = special.erfc(r / 2.*r_s) + (r / (r_s * math.sqrt(math.pi))) * exp(-r**2 / (4.*r_s**2))
-eta_short_swift = (2. * r * d_swift_corr(r / (2.*r_s)) / r_s) - 2*sigmoid(2*r / (r_s)) + 2.
+eta_short_swift = 4. * (r / r_s) * d_sigmoid(2. * r / r_s) - 2. * sigmoid(2 * r / r_s) + 2.
+
+# Corection in Fourier space
+corr_long_gadget2 = exp(-k**2*r_s**2)
+corr_long_swift = math.pi * k * r_s * csch(0.5 * math.pi * r_s * k) / 2.
 
-# Shortrange term2
+# Shortrange term
 phi_short_gadget2 = (1.  / r ) * (1. - corr_short_gadget2)
 phi_short_swift = (1.  / r ) * (1. - corr_short_swift)
 force_short_gadget2 = (1. / r**2) * eta_short_gadget2
@@ -149,18 +135,17 @@ legend(loc="upper right", frameon=True, handletextpad=0.1, handlelength=3.2, fon
 
 # Correction
 subplot(312, xscale="log", yscale="log")
-#plot(r_rs, np.ones(np.size(r)), '--', lw=1.4, color=colors[0])
+plot(r_rs, np.ones(np.size(r)), '--', lw=1.4, color=colors[0])
 plot(r_rs, 1. - corr_short_gadget2, '-', lw=1.4, color=colors[2])
 plot(r_rs, 1. - corr_short_swift, '-', lw=1.4, color=colors[3])
-#plot(r_rs, np.zeros(np.size(r)), 'k--', alpha=0.5, lw=0.5)
-plot(r_rs, np.ones(np.size(r)), 'k--', alpha=0.5, lw=0.5)
-plot(r_rs, np.ones(np.size(r))*0.01, 'k--', alpha=0.5, lw=0.5)
+plot(r_rs, np.ones(np.size(r))*0.01, 'k:', alpha=0.5, lw=0.5)
 plot([1., 1.], [-1e5, 1e5], 'k-', alpha=0.5, lw=0.5)
 
 yticks([1e-2, 1e-1, 1], ["$0.01$", "$0.1$", "$1$"])
 xlim(1.1*r_min/r_s, 0.9*r_max/r_s)
 ylim(3e-3, 1.5)
-ylabel("$\\chi_s(r)$", labelpad=-3)
+#ylabel("$\\chi_s(r)$", labelpad=-3)
+ylabel("$\\varphi_s(r) \\times r$", labelpad=-2)
 
 # 1 - Correction
 subplot(313, xscale="log", yscale="log")
@@ -168,12 +153,13 @@ plot(r_rs, corr_short_gadget2, '-', lw=1.4, color=colors[2])
 plot(r_rs, corr_short_swift, '-', lw=1.4, color=colors[3])
 
 plot([1., 1.], [1e-5, 1e5], 'k-', alpha=0.5, lw=0.5)
-plot(r_rs, np.ones(np.size(r)), 'k--', alpha=0.5, lw=0.5)
-plot(r_rs, np.ones(np.size(r))*0.01, 'k--', alpha=0.5, lw=0.5)
+plot(r_rs, np.ones(np.size(r)), 'k:', alpha=0.5, lw=0.5)
+plot(r_rs, np.ones(np.size(r))*0.01, 'k:', alpha=0.5, lw=0.5)
 
 xlim(1.1*r_min/r_s, 0.9*r_max/r_s)
 ylim(3e-3, 1.5)
-ylabel("$1 - \\chi_s(r)$", labelpad=-2)
+#ylabel("$1 - \\chi_s(r)$", labelpad=-2)
+ylabel("$1 - \\varphi_s(r) \\times r$", labelpad=-2)
 yticks([1e-2, 1e-1, 1], ["$0.01$", "$0.1$", "$1$"])
 xlabel("$r / r_s$", labelpad=-3)
 
@@ -200,16 +186,17 @@ legend(loc="upper right", frameon=True, handletextpad=0.1, handlelength=3.2, fon
 
 # Correction
 subplot(312, xscale="log", yscale="log")
+plot(r_rs, np.ones(np.size(r)), '--', lw=1.4, color=colors[0])
 plot(r_rs, eta_short_gadget2, '-', lw=1.4, color=colors[2])
 plot(r_rs, eta_short_swift, '-', lw=1.4, color=colors[3])
-plot(r_rs, np.ones(np.size(r)), 'k--', alpha=0.5, lw=0.5)
-plot(r_rs, np.ones(np.size(r))*0.01, 'k--', alpha=0.5, lw=0.5)
+plot(r_rs, np.ones(np.size(r))*0.01, 'k:', alpha=0.5, lw=0.5)
 plot([1., 1.], [-1e5, 1e5], 'k-', alpha=0.5, lw=0.5)
 
 yticks([1e-2, 1e-1, 1], ["$0.01$", "$0.1$", "$1$"])
 xlim(1.1*r_min/r_s, 0.9*r_max/r_s)
 ylim(3e-3, 1.5)
-ylabel("$\\eta_s(r)$", labelpad=-3)
+#ylabel("$\\eta_s(r)$", labelpad=-3)
+ylabel("$|\\mathbf{f}_s(r)|\\times r^2$", labelpad=-2)
 
 # 1 - Correction
 subplot(313, xscale="log", yscale="log")
@@ -217,12 +204,13 @@ plot(r_rs, 1. - eta_short_gadget2, '-', lw=1.4, color=colors[2])
 plot(r_rs, 1. - eta_short_swift, '-', lw=1.4, color=colors[3])
 
 plot([1., 1.], [1e-5, 1e5], 'k-', alpha=0.5, lw=0.5)
-plot(r_rs, np.ones(np.size(r)), 'k--', alpha=0.5, lw=0.5)
-plot(r_rs, np.ones(np.size(r))*0.01, 'k--', alpha=0.5, lw=0.5)
+plot(r_rs, np.ones(np.size(r)), 'k:', alpha=0.5, lw=0.5)
+plot(r_rs, np.ones(np.size(r))*0.01, 'k:', alpha=0.5, lw=0.5)
 
 xlim(1.1*r_min/r_s, 0.9*r_max/r_s)
 ylim(3e-3, 1.5)
-ylabel("$1 - \\eta_s(r)$", labelpad=-2)
+#ylabel("$1 - \\eta_s(r)$", labelpad=-2)
+ylabel("$1 - |\\mathbf{f}_s(r)|\\times r^2$", labelpad=-3)
 yticks([1e-2, 1e-1, 1], ["$0.01$", "$0.1$", "$1$"])
 xlabel("$r / r_s$", labelpad=-3)
 
@@ -231,7 +219,7 @@ savefig("force_short.pdf")
 ##################################################################################################
 
 figure()
-subplot(211, xscale="log", yscale="log")
+subplot(311, xscale="log", yscale="log")
 
 # Potential
 plot(k_rs, phit_newton, '--', lw=1.4, label="${\\rm Newtonian}$", color=colors[0])
@@ -245,22 +233,35 @@ legend(loc="lower left", frameon=True, handletextpad=0.1, handlelength=3.2, font
 xlim(1.1*r_min/ r_s, 0.9*r_max / r_s)
 ylim(1.1/r_max**2, 0.9/r_min**2)
 ylabel("$\\tilde{\\varphi_l}(k)$", labelpad=-3)
+yticks([1e-4, 1e-2, 1e0, 1e2], ["$10^{-4}$", "$10^{-2}$", "$10^{0}$", "$10^{2}$"])
 
-
-subplot(212, xscale="log", yscale="log")
+subplot(312, xscale="log", yscale="log")
 
 # Potential normalized
 plot(k_rs, phit_newton * k**2, '--', lw=1.4, label="${\\rm Newtonian}$", color=colors[0])
 plot(k_rs, phit_long_gadget2 * k**2, '-', lw=1.4, label="${\\rm Gadget}$", color=colors[2])
 plot(k_rs, phit_long_swift * k**2, '-', lw=1.4, label="${\\rm SWIFT}$", color=colors[3])
-plot(k_rs, -phit_long_swift * k**2, ':', lw=1.4, label="${\\rm SWIFT}$", color=colors[3])
 plot([1., 1.], [1e-5, 1e5], 'k-', alpha=0.5, lw=0.5)
-plot(r_rs, np.ones(np.size(r))*0.01, 'k--', alpha=0.5, lw=0.5)
+plot(r_rs, np.ones(np.size(r))*0.01, 'k:', alpha=0.5, lw=0.5)
 
 xlim(1.1*r_min/ r_s, 0.9*r_max / r_s)
 ylim(3e-3, 1.5)
 ylabel("$k^2 \\times \\tilde{\\varphi_l}(k)$", labelpad=-3)
 yticks([1e-2, 1e-1, 1], ["$0.01$", "$0.1$", "$1$"])
+
+subplot(313, xscale="log", yscale="log")
+
+plot(k_rs, 1. - phit_long_gadget2 * k**2, '-', lw=1.4, label="${\\rm Gadget}$", color=colors[2])
+plot(k_rs, 1. - phit_long_swift * k**2, '-', lw=1.4, label="${\\rm SWIFT}$", color=colors[3])
+plot([1., 1.], [1e-5, 1e5], 'k-', alpha=0.5, lw=0.5)
+plot(r_rs, np.ones(np.size(r)), 'k:', alpha=0.5, lw=0.5)
+plot(r_rs, np.ones(np.size(r))*0.01, 'k:', alpha=0.5, lw=0.5)
+
+xlim(1.1*r_min/ r_s, 0.9*r_max / r_s)
+ylim(3e-3, 1.5)
+ylabel("$1 - k^2 \\times \\tilde{\\varphi_l}(k)$", labelpad=-3)
+yticks([1e-2, 1e-1, 1], ["$0.01$", "$0.1$", "$1$"])
+
 xlabel("$k \\times r_s$", labelpad=0)
 
 savefig("potential_long.pdf")