begin cleaning up

green.pyx

@@ -37,8 +37,8 @@ cdef libcpp.map.map[double, Entry] cache
 
 cdef object sys
 cdef int n, m, nleads, norbs, index
-cdef double t, fermi_energy
-cdef double [:] pots, temps
+cdef double t
+cdef double [:] fermi_energies, temps
 
 
 cdef double real_f(double e, void* ignore):
@@ -71,9 +71,9 @@ cdef double complex complex_quad(double a, double b,
 
 cdef double fermi_of_lead(double E, int lead):
     if temps[lead] == 0:
-        return float(pots[lead] + fermi_energy - E > 0)
+        return float(fermi_energies[lead] - E > 0)
     else:
-        return 1. / (exp((E - (pots[lead] + fermi_energy) ) / temps[lead]) + 1)
+        return 1. / (exp((E - (fermi_energies[lead]) ) / temps[lead]) + 1)
 
 
 cdef double complex func(double E):
@@ -125,21 +125,19 @@ cdef double complex func(double E):
             g += (cached[0].psis[i * norbs + n] *
                   cached[0].psis[i * norbs + m].conjugate() * f)
     sincos(-E * t, &sin, &cos)
-    g *= (cos + 1j * sin)
+    g *= 1j * (cos + 1j * sin) / (2*np.pi)
     return g
 
 
-def green(sys_, norbs_, interval, pots_, temps_, fermi_energy_, times,
+def green(sys_, norbs_, interval, fermi_energies_, temps_, times,
           epsabs=1e-6, epsrel=1e-6, limit=1000):
-    global cache, sys, n, m, norbs, nleads, index, t, fermi_energy, \
-        pots, temps, workspace
+    global cache, sys, n, m, norbs, nleads, index, t, fermi_energies, temps, workspace
 
     cache.clear()
     sys = sys_
     norbs = norbs_
     nleads = len(sys.leads)
-    fermi_energy = fermi_energy_
-    pots = np.asarray(pots_)
+    fermi_energies = np.asarray(fermi_energies_)
     temps = np.asarray(temps_)
 
     try:
@@ -147,22 +145,13 @@ def green(sys_, norbs_, interval, pots_, temps_, fermi_energy_, times,
         if workspace == NULL:
             raise RuntimeError()
 
-        F = np.empty((norbs, norbs, len(times)), complex)
-        G = np.empty((norbs, norbs, len(times)), complex)
+        G = np.empty((len(times), 2, norbs, norbs), complex)
         for n in xrange(norbs):
             for m in xrange(norbs):
-                j = 0
-                for t in times:
-                    index = 0
-                    F[n, m, j] = complex_quad(interval[0], interval[1],
-                                              epsabs, epsrel, limit)
-                    index = 1
-                    G[n, m, j] = complex_quad(interval[0], interval[1],
-                                              epsabs, epsrel, limit)
-                    j += 1
-
-        F *= 1j / (2*np.pi)
-        G *= 1j / (2*np.pi)
+                for j, t in enumerate(times):
+                    for index in [0, 1]:
+                        G[j, index, n, m] = complex_quad(interval[0], interval[1],
+                                                         epsabs, epsrel, limit)
     finally:
         gsl_integration_workspace_free(workspace)
-    return F, G
+    return G

test.py

@@ -30,9 +30,8 @@ def make_system(a=1, t=1.0, W=10, L=30):
 def test():
     syst = make_system(W=2, L=1).finalized()
     # kwant.plotter.bands(syst.leads[0])
-    G_lesser, G_upper = green.green(syst, 2, [0., 8.], [2., 3.], [0., 0.], 0., [0.])
-    Gl = G_lesser[:, :, 0]
-    Gu = G_upper[:, :, 0]
+    G = green.green(syst, 2, [0., 8.], [2., 3.], [0., 0.], [0.])
+    Gl, Gu = G[0]
     np.testing.assert_almost_equal(-1j * (Gl - Gu), np.identity(len(Gl)))
 
 if __name__ == "__main__":