Added functional form of the gravity long-range force correction to the plotting script.

.gitignore

@@ -104,6 +104,7 @@ theory/Multipoles/fmm_standalone.pdf
 theory/Multipoles/potential.pdf
 theory/Multipoles/potential_long.pdf
 theory/Multipoles/potential_short.pdf
+theory/Multipoles/force_short.pdf
 
 m4/libtool.m4
 m4/ltoptions.m4

theory/Multipoles/plot_mesh.py

@@ -65,6 +65,7 @@ k_rs = k * r_s
 # Newtonian solution
 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
@@ -73,15 +74,28 @@ def smoothstep(x): #S_2(x)
     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)
+
 def sigmoid(x):
-    return 1. / (1. + exp(-x))
+    return 1. / (1. + 1./my_exp(x))
     #return x / sqrt(1. + x**2)
-    
+
+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])
@@ -89,18 +103,24 @@ 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")
 
 # 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 = (-15. / 1024.) * (12. * r_s * k * cos(4. * r_s * k) + (16. * r_s**2 * k**2 - 3.) * sin(4. * r_s * k)) / (k**5 * r_s**5)
-corr_long_swift = k * r_s * math.pi / (2. * sinh(0.5 * math.pi * r_s * k))
+#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.
 
-# Shortrange term
+# Shortrange term2
 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
+force_short_swift = (1. / r**2) * eta_short_swift
 
 # Long-range term
 phi_long_gadget2 = (1.  / r ) * corr_short_gadget2
@@ -108,6 +128,9 @@ phi_long_swift = (1.  / r ) * corr_short_swift
 phit_long_gadget2 = corr_long_gadget2 / k**2
 phit_long_swift = corr_long_swift / k**2
 
+
+
+
 figure()
 
 # Potential
@@ -139,15 +162,11 @@ xlim(1.1*r_min/r_s, 0.9*r_max/r_s)
 ylim(3e-3, 1.5)
 ylabel("$\\chi_s(r)$", labelpad=-3)
 
+# 1 - Correction
 subplot(313, xscale="log", yscale="log")
-
-#print corr_short_gadget2
-
-#plot(r_rs, np.abs(1. - np.ones(np.size(r))), '--', lw=1.4, color=colors[0])
 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)
@@ -158,12 +177,58 @@ ylabel("$1 - \\chi_s(r)$", labelpad=-2)
 yticks([1e-2, 1e-1, 1], ["$0.01$", "$0.1$", "$1$"])
 xlabel("$r / r_s$", labelpad=-3)
 
-#ylim(0, 0.95)
-
 savefig("potential_short.pdf")
 
+##################################################################################################
+
+
+# Force
+figure()
+subplot(311, xscale="log", yscale="log")
+
+plot(r_rs, force_newton, '--', lw=1.4, label="${\\rm Newtonian}$", color=colors[0])
+plot(r_rs, force_short_gadget2, '-', lw=1.4, label="${\\rm Gadget}$", color=colors[2])
+plot(r_rs, force_short_swift, '-', lw=1.4, label="${\\rm SWIFT}$", color=colors[3])
+plot([1., 1.], [1e-5, 1e5], 'k-', alpha=0.5, lw=0.5)
+
+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("$|\\mathbf{f}_s(r)|$", labelpad=-3)
+yticks([1e-4, 1e-2, 1e0, 1e2], ["$10^{-4}$", "$10^{-2}$", "$10^{0}$", "$10^{2}$"])
+
+legend(loc="upper right", frameon=True, handletextpad=0.1, handlelength=3.2, fontsize=8)
+
+# Correction
+subplot(312, xscale="log", yscale="log")
+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([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)
+
+# 1 - Correction
+subplot(313, xscale="log", yscale="log")
+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)
+
+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)
+yticks([1e-2, 1e-1, 1], ["$0.01$", "$0.1$", "$1$"])
+xlabel("$r / r_s$", labelpad=-3)
 
+savefig("force_short.pdf")
 
+##################################################################################################
 
 figure()
 subplot(211, xscale="log", yscale="log")