examples: make cpr_sns produce a plot

examples/cpr_sns.py

@@ -15,6 +15,8 @@ from usadelndsoc.matsubara import get_matsubara_sum
 from usadelndsoc.solver import *
 from usadelndsoc.util import vectorize_parallel
 
+import collections
+
 mem = joblib.Memory("cache")
 
 
@@ -49,9 +51,13 @@ def get_solver(soc_alpha, eta, phi, L=10, h=0.5, D=1.0, n=5):
     return sol
 
 
-# @usadelndsoc.with_log_level(logging.WARNING)
+Res = collections.namedtuple("Res", ["x", "y", "J", "Jx", "Jy"])
+
+
 @vectorize_parallel(returns_object=True, noarray=True)
+@usadelndsoc.with_log_level(logging.WARNING)
 def j(T, h, phi, n=15, eta=0.1, alpha_soc=0.1, L=10, perturbative=False):
+
     sol = get_solver(soc_alpha=alpha_soc, eta=eta, L=L, n=n, phi=phi, h=h)
 
     E_typical = 10 + 10 * T
@@ -67,29 +73,24 @@ def j(T, h, phi, n=15, eta=0.1, alpha_soc=0.1, L=10, perturbative=False):
     Jx = (J[:, :, 0] + J[:, :, 2]).real / 2
     Jy = (J[:, :, 1] + J[:, :, 3]).real / 2
 
-    return sol.x, sol.y, J, Jx, Jy
+    return Res(sol.x, sol.y, J, Jx, Jy)
 
 
 def main():
-    T = 0.02
-    alpha_soc = 0.05
-    Da2 = alpha_soc**2
-    h = np.linspace(0, 1.75, 29)
+    T = 0.1
+    alpha_soc = 0.1
+    h = 0.5
+    phi = np.linspace(-pi, pi, 37)
 
-    res = j_anom(T, h, alpha_soc=alpha_soc, perturbative=False, mem=mem)
-    j = np.asarray([x[0] for x in res])
-    j_an = np.asarray([x[1] for x in res])
+    res = j(T, h, phi, alpha_soc=alpha_soc, perturbative=False, mem=mem)
 
-    t3 = np.diag([1, 1, -1, -1])
-    j = tr(j @ t3)[:, 2, 2, UP]
+    Jx_mean = np.asarray([x.Jx[1:-1].sum(axis=1).mean(axis=0) for x in res])
 
-    plt.plot(h, -j.real, "-", label=r"numerics")
-    plt.plot(h, -j_an.real, "--", label=r"analytic")
-    plt.title(rf"$T = {T} \Delta$, $D\alpha^2 = {Da2:.4} \Delta$")
-    plt.xlabel(r"$h/\Delta$")
-    plt.ylabel(r"$j_{\mathrm{an}}$")
+    plt.plot(phi / pi, Jx_mean)
+    plt.xlabel(r"$\varphi / \pi$")
+    plt.ylabel(r"$I$")
     plt.legend()
-    plt.savefig("bulk_j_an.pdf")
+    plt.savefig("cpr_sns.pdf")
 
 
 if __name__ == "__main__":