maroba / multipoles Goto Github PK

A Python package for multipole expansions of electrostatic or gravitational potentials

License: MIT License

Python 100.00%

physics scientific-computing mathematics electrostatic-potential

multipoles's Issues

Multipoles in cartesian coordinates

Hello! I was wondering if there is any way to extract the multipoles in cartesian coordinates rather than spherical coordinates.
Thanks!

Phi(x,y,z) gives different answers if entire system is shifted

charge_dist = {'discrete': True, 'charges': [{'q': 1, 'xyz': (0, 0, 0.5)}, {'q': -1, 'xyz': (0, 0, -1)}]}
l_max = 3
Phi = MultipoleExpansion(charge_dist, l_max)
Phi(0,0,0)

and

charge_dist = {'discrete': True, 'charges': [{'q': 1, 'xyz': (5, 5, 5.5)}, {'q': -1, 'xyz': (5, 5, 4)}]}
l_max = 3
Phi = MultipoleExpansion(charge_dist, l_max)
Phi(5,5,5)

give different answers despite being the same arrangement of charges around the evaluation point, is this correct?

Document error bounds

So Beatson/Greengard list the L1 error as bounded here:

$$L1 < frac{\sum_i |q_i|}{r-a} (\frac{a}{r})^{lmax+1}$$

But this is the exterior expansion. I don't have as many resources on the interior expansion.

I think it would be nice to show

This project keeps these bound
Add documentation on the bounds

How would one adapt this to perform interior instead of exterior expansions?

As stated in the title. I'm interested in expanding distributions from a central point r onto a spherical interior. Any help would be appreciated.

Inaccurate results

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.

Feature request: evaluate multipole expansion on single axis

Discussed in #10

^{Originally posted by nordic-node July 14, 2023}
Hi,

instead of evaluating the multipole expansion point by point, like

phi = np.array([mpe(r, 0, 0, l_max=l_max) for r in rr])

where mpe is a multipole expansion object.

I would rather like to write something like this:

phi = mpe(rr, 0, 0, l_max=l_max)

As far as I can see, the call function just allows to pass a point. and the getitem function wants slices or masks on the same grid
that was used for the charge distribution. Is there a possibility to do the same for outside of that grid?

Cheers, Sara

Weird Harmonic Artifacts

Discussed in #15

^{Originally posted by levi2234 December 5, 2023}

Astrometric Halo Merger artifacts

While applying the Multipole Expansion on a set of 2 Gaussian distributed clusters I encountered a weird artifact dependent on the value of L_min. When the Gaussian cluster has a separation of multiple cluster standard deviations a deep potential well occurs exactly between the two clusters. As the clusters move inward this well resolves itself. I would like to know if this is a problem that should be adressed or if it is expected behaviour. If it is the latter I think I must look to some other solution.

To sketch the problem I have included a video mock collision for multiple modes for both haloes and single points. As a reference I have also included an discrete example of what I would expect the potential field to look like.

HaloCollisionMultipoleEvolution_l_max_4.mp4