Stricter testing

README.md

@@ -5,7 +5,7 @@
 
 ## Background
 
-For a given function $\rho(x,y,z)$, the solution $\Ph   i(x,y,z)$ of the Poisson equation $\nabla^2\Phi=4\pi \rho$ with vanishing Dirichlet boundary conditions at infinity is
+For a given function $\rho(x,y,z)$, the solution $\Ph   i(x,y,z)$ of the Poisson equation $\nabla^2\Phi=-4\pi \rho$ with vanishing Dirichlet boundary conditions at infinity is
 
 $$\Phi(x,y,z)=-\int d^3r'\frac{\rho(r')}{|r-r'|}$$
 

multipoles/expansion.py

@@ -291,16 +291,15 @@ def _multipole_contribs(self, xyz):
                 else:
                     phi_l += np.sqrt(4 * np.pi / (2 * l + 1)) * q_lm * Y_lm * r ** l
             mp_contribs.append(phi_l.real)
-
         return mp_contribs
 
     def _calc_multipole_moments(self):
         moments = {}
         for l in range(0, self.l_max + 1):
-            for m in range(0, l + 1):
+            for m in range(-l, l + 1):
                 moments[(l, m)] = self._calc_multipole_coef(l, m)
-                if m != 0:
-                    moments[(l, -m)] = (-1) ** m * np.conj(moments[(l, m)])
+                #if m != 0:
+                #    moments[(l, -m)] = (-1) ** m * np.conj(moments[(l, m)])
         return moments
 
     def _calc_multipole_coef(self, l, m):
@@ -375,15 +374,15 @@ def cartesian_to_spherical(*coords):
     np.seterr(all='ignore')
 
     if hasattr(R_xy, '__len__'):
-        Phi = np.arccos(X / R_xy)
+        Phi = np.arctan2(Y, X)
         Phi[R_xy == 0] = 0
         Theta = np.arccos(Z / R)
         Theta[R == 0] = np.pi / 2
     else:
         if R_xy == 0:
             Phi = 0
         else:
-            Phi = np.arccos(X / R_xy)
+            Phi = np.arctan2(Y, X)
         if R == 0:
             Theta = np.pi / 2
         else:

tests/test_expansion.py

@@ -8,40 +8,79 @@
 from scipy.special import erf
 
 from multipoles.expansion import (
-    MultipoleExpansion, InvalidChargeDistributionException, InvalidExpansionException
+    MultipoleExpansion, InvalidChargeDistributionException, InvalidExpansionException, cartesian_to_spherical
 )
 
 
 class TestMultipoleExpansion(unittest.TestCase):
 
     def test_gaussian_monopole_at_center(self):
-        x, y, z = [np.linspace(-10, 10, 51)] * 3
+        x, y, z = [np.linspace(-7, 7, 51)] * 3
         X, Y, Z = np.meshgrid(x, y, z, indexing='ij')
         sigma = 1.5
         rho = gaussian((X, Y, Z), (0, 0, 0), sigma)
-        mpe = MultipoleExpansion(charge_dist={'discrete': False, 'rho': rho, 'xyz': (X, Y, Z)}, l_max=8)
+        mpe = MultipoleExpansion(charge_dist={'discrete': False, 'rho': rho, 'xyz': (X, Y, Z)}, l_max=4)
 
         self.assertAlmostEqual(1, mpe.total_charge, places=4)
         np.testing.assert_array_almost_equal(mpe.center_of_charge, (0, 0, 0))
 
         self.assertAlmostEqual(0.1, mpe._multipole_contribs((10, 0, 0))[0], places=4)
 
-        #R = np.sqrt(X**2 + Y**2 + Z**2)
-        #b = np.ones_like(R, dtype=bool)
-        #b[1:-1, 1:-1, 1:-1] = False
-        #npt.assert_array_almost_equal(
-        #    erf(R[b] / sigma) / R[b],
-        #    mpe[b]
-        #)
+        R = np.sqrt(X ** 2 + Y ** 2 + Z ** 2)
+        b = np.ones_like(R, dtype=bool)
+        b[1:-1, 1:-1, 1:-1] = False
+
+        r = 20
+        self.assertAlmostEqual(
+            erf(r / sigma) / r,
+            mpe(r, 0, 0), places=5
+        )
+        self.assertAlmostEqual(
+            erf(r / sigma) / r,
+            mpe(0, r, 0), places=5
+        )
+        self.assertAlmostEqual(
+            erf(r / sigma) / r,
+            mpe(0, 0, r), places=5
+        )
 
     def test_gaussian_monopole_at_off_center(self):
         x, y, z = [np.linspace(-5, 5, 51)] * 3
         X, Y, Z = np.meshgrid(x, y, z, indexing='ij')
         sigma = 1.5
         rho = gaussian((X, Y, Z), (1, 0, 0), sigma)
-        mpe = MultipoleExpansion(charge_dist={'discrete': False, 'rho': rho, 'xyz': (X, Y, Z)}, l_max=2)
+        mpe = MultipoleExpansion(charge_dist={'discrete': False, 'rho': rho, 'xyz': (X, Y, Z)}, l_max=8)
 
-        self.assertAlmostEqual(0.1, mpe._multipole_contribs((11, 0, 0))[0], places=4)
+        def dist(xyz1, xyz2):
+            return np.sqrt(sum((np.array(xyz1) - np.array(xyz2)) ** 2))
+
+        P = (20, 0, 0)
+        r = dist(P, (1, 0, 0))
+        npt.assert_almost_equal(
+            erf(r / sigma) / r,
+            mpe(*P), decimal=5
+        )
+
+        P = (0, 20, 0)
+        r = dist(P, (1, 0, 0))
+        npt.assert_almost_equal(
+            erf(r / sigma) / r,
+            mpe(*P), decimal=5
+        )
+
+        P = (0, 0, 20)
+        r = dist(P, (1, 0, 0))
+        npt.assert_almost_equal(
+            erf(r / sigma) / r,
+            mpe(*P), decimal=5
+        )
+
+        P = (20, 20, 20)
+        r = dist(P, (1, 0, 0))
+        npt.assert_almost_equal(
+            erf(r / sigma) / r,
+            mpe(*P), decimal=5
+        )
 
     def test_gaussian_dipole_at_center(self):
         x, y, z = [np.linspace(-5, 5, 51)] * 3
@@ -247,6 +286,36 @@ def test_issue_7(self):
 
         self.assertAlmostEqual(phi_1, phi_2)
 
+    def test_internal_coordinates_spherical(self):
+        x, y, z = [np.linspace(-5, 5, 51)] * 3
+        X, Y, Z = np.meshgrid(x, y, z, indexing='ij')
+        sigma = 1.5
+        rho = gaussian((X, Y, Z), (0, 0, 0), sigma)
+        mpe = MultipoleExpansion(charge_dist={'discrete': False, 'rho': rho, 'xyz': (X, Y, Z)}, l_max=3)
+
+        R, Phi, Theta = mpe.internal_coords_spherical
+        self.assertAlmostEqual(np.max(R), np.sqrt(3) * 5)
+        self.assertAlmostEqual(np.min(Phi), -np.pi)
+        self.assertAlmostEqual(np.max(Theta), np.pi)
+
+
+class TestCartesianToSpherical(unittest.TestCase):
+
+    def test_works(self):
+        pi = np.pi
+
+        actual = cartesian_to_spherical(1, 0, 0)
+        expected = (1, 0, pi / 2)
+        npt.assert_array_equal(expected, actual)
+
+        actual = cartesian_to_spherical(0, 1, 0)
+        expected = (1, pi / 2, pi / 2)
+        npt.assert_array_equal(expected, actual)
+
+        actual = cartesian_to_spherical(0, -1, 0)
+        expected = (1, -pi / 2, pi / 2)
+        npt.assert_array_equal(expected, actual)
+
 
 def gaussian(XYZ, xyz0, sigma):
     g = np.ones_like(XYZ[0])