examples: cpr_sns: curves for different lattice spacing

examples/bulk_j_an.py

@@ -50,11 +50,11 @@ def j_anom(T, h, n=5, alpha_soc=0.01, perturbative=False):
 
     # Solver without SOC
     if perturbative:
-        sol0 = get_solver(soc_alpha=0, eta=0, h=h)
+        sol0 = get_solver(soc_alpha=0, eta=0, h=h, n=n)
 
     # Solver for evaluating current perturbatively in SOC
     eta = 1.0
-    sol = get_solver(soc_alpha=alpha_soc, eta=eta, h=h, D=1.0)
+    sol = get_solver(soc_alpha=alpha_soc, eta=eta, h=h, D=1.0, n=n)
 
     E_typical = 10 + 10 * T
     w, a = get_matsubara_sum(T, E_typical)
@@ -104,9 +104,9 @@ def main():
     T = 0.02
     alpha_soc = 0.05
     Da2 = alpha_soc**2
-    h = np.linspace(0, 1.75, 29)
+    h = np.linspace(0, 1.75, 79)
 
-    res = j_anom(T, h, alpha_soc=alpha_soc, perturbative=False, mem=mem)
+    res = j_anom(T, h, n=20, 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])
 

examples/cpr_sns.py

@@ -79,21 +79,49 @@ def j(T, h, phi, n=15, eta=0.1, alpha_soc=0.1, L=10):
 def main():
     T = 0.1
     alpha_soc = 0.9
-    h = np.r_[0.0, 0.25, 0.5, 1.0, 1.5]
+    h = np.r_[0.0, 1.0]
     phi = np.linspace(-pi, pi, 37)
 
-    res = j(T, h[:, None], phi[None, :], alpha_soc=alpha_soc, L=3, mem=mem)
+    res = j(T, h[:, None], phi[None, :], alpha_soc=alpha_soc, L=3, n=20, mem=mem)
+    res0 = j(T, h[:, None], phi[None, :], alpha_soc=alpha_soc, L=3, mem=mem)
+    res1 = j(T, h[:, None], phi[None, :], alpha_soc=alpha_soc, L=3, n=10, mem=mem)
 
     Jx_mean = np.asarray(
         [Res(*x).Jx[1:-1].sum(axis=1).mean(axis=0) for x in res.flat]
     ).reshape(res.shape)
 
+    Jx0_mean = np.asarray(
+        [Res(*x).Jx[1:-1].sum(axis=1).mean(axis=0) for x in res0.flat]
+    ).reshape(res.shape)
+
+    Jx1_mean = np.asarray(
+        [Res(*x).Jx[1:-1].sum(axis=1).mean(axis=0) for x in res1.flat]
+    ).reshape(res.shape)
+
     plt.plot(phi / pi, Jx_mean.T)
+    plt.plot(phi / pi, Jx0_mean.T, 'k:')
+    plt.plot(phi / pi, Jx1_mean.T, 'k:')
     plt.xlabel(r"$\varphi / \pi$")
     plt.ylabel(r"$I$")
     plt.legend(h, title="$h$")
     plt.savefig("cpr_sns.pdf")
 
+    Jx_sym = (Jx_mean + Jx_mean[:,::-1])/2
+    Jx0_sym = (Jx0_mean + Jx0_mean[:,::-1])/2
+    Jx1_sym = (Jx1_mean + Jx1_mean[:,::-1])/2
+    plt.clf()
+    plt.plot(phi / pi, Jx_sym.T)
+    plt.plot(phi / pi, Jx0_sym.T, 'k:')
+    plt.plot(phi / pi, Jx1_sym.T, 'k:')
+    plt.xlabel(r"$\varphi / \pi$")
+    plt.ylabel(r"$[I(\varphi)+I(-\varphi)]/2$")
+    plt.legend(h, title="$h$")
+    plt.savefig("cpr_sns_sym.pdf")
+
+    print(Jx_mean.mean(axis=1))
+    print(Jx0_mean.mean(axis=1))
+    print(Jx1_mean.mean(axis=1))
+
 
 if __name__ == "__main__":
     main()