Inaccurate results

[leon-vv]

Describe the bug

The packages give inaccurate answers for some charge distributions.

To Reproduce

Consider the following code.

import numpy as np
from multipoles import MultipoleExpansion

# Prepare the charge distribution dict for the MultipoleExpansion object:
charge_dist = {
    'discrete': True,     # point charges are discrete charge distributions
    'charges': [
        {'q': 1, 'xyz': (0, 0, 1)},
        {'q': -1.53, 'xyz': (0, 0, -1)},
        {'q': -1, 'xyz': (0, 1.3, -1)}
    ]
}
def true_potential(x, y, z):
    sum_ = 0.0

    for c in charge_dist['charges']:
        p = np.array(c['xyz'])
        r = np.linalg.norm(np.array([x,y,z]) - p)
        sum_ += c['q']/r

    return sum_

x, y, z = 30.5, 30.6, 30.7
Phi = MultipoleExpansion(charge_dist, 15)
print(true_potential(x, y, z))
print(Phi(x, y, z))
print('Accuracy: ', Phi(x, y, z)/true_potential(x, y, z) - 1)

The accuracy is only 5e-3, and does not improve when increasing l_max.

Expected behavior

I expect an accuracy down to floating point limit (1e-15) for large value of l_max.

Desktop (please complete the following information):

Linux, Debian, Python 3.9

Additional context

I believe there is a sign error related to the spherical harmonics functions used by this package. To be specific, I believe for some values of l, m the sign of the imaginary part of the spherical harmonics functions used is incorrect. I'm still investigating.. I will let you know when I figure it out.

[leon-vv]

Here is a multipole expansion that does give accurate results, within the theoretical error bounds:

import math

import numpy as np
from scipy.special import lpmv

# Prepare the charge distribution dict for the MultipoleExpansion object:
charge_dist = {
    'discrete': True,     # point charges are discrete charge distributions
    'charges': [
        {'q': 1, 'xyz': (0, 0, 1)},
        {'q': -1.53, 'xyz': (0, 0, -1)},
        {'q': -1, 'xyz': (0, 3.2, -1)}
    ]
}

charges = np.array([c['q'] for c in charge_dist['charges']])
positions = np.array([c['xyz'] for c in charge_dist['charges']])

def true_potential(x, y, z):
    sum_ = 0.0

    for c in charge_dist['charges']:
        p = np.array(c['xyz'])
        r = np.linalg.norm(np.array([x,y,z]) - p)
        sum_ += c['q']/r

    return sum_

l_max = 8
lm_indices = np.array([(l, m) for l in range(0, l_max+1) for m in range(-l, l+1)])

def harmonic(l, m, phi, theta):
    f = math.factorial
    return math.sqrt( f(l - abs(m)) / f(l + abs(m)) ) * lpmv(abs(m), l, np.cos(theta)) * np.exp(1j * m * phi)

def cartesian_to_spherical(x, y, z):
    r = np.linalg.norm([x,y,z])

    if(np.isclose(r, 0.)):
        return 0., 0., np.pi/2

    theta = math.acos(z/r)
    phi = math.atan2(y, x)
    return r, phi, theta

def multipole_coeff(charges, positions, l, m):
    N = len(charges)
    assert charges.shape == (N,)
    assert positions.shape == (N, 3)

    coeff = 0.

    for q, p in zip(charges, positions):
        r, phi, theta = cartesian_to_spherical(*p)
        coeff += q * r**l * harmonic(l, -m, phi, theta)

    return coeff

def multipole_expansion(charges, positions):
    exp = np.zeros( len(lm_indices), dtype=np.cdouble )

    for i, (l, m) in enumerate(lm_indices):
        exp[i] = multipole_coeff(charges, positions, l, m)

    return exp

def expansion_to_voltage(expansion, target):
    voltage = 0.

    for i, (l, m) in enumerate(lm_indices):
        r, phi, theta = cartesian_to_spherical(*target)
        voltage += harmonic(l, m, phi, theta) * expansion[i] / (r**(l+1))

    return np.real(voltage)

expansion = multipole_expansion(charges, positions)

x, y, z = 30.5, 30.6, 30.7

correct = true_potential(x, y, z)
print('Accuracy: ', expansion_to_voltage(expansion, (x, y, z))/true_potential(x, y, z) - 1)

max_x = np.max(np.linalg.norm(positions, axis=1))
y = np.linalg.norm( (x, y, z) )
print('Error bound: ', np.sum(np.abs(charges))/(y - max_x) * (max_x/y)**l_max)