I want to use vectorized computation (with python and numpy) to do some calculation involving bessel functions.
As is often the case, I need to calculate the sums of bessel functions of (infinit) differnt order. I need to do that for a whole vector field, and would like to use vectorized computation. Actually, both, the bessel orders and the dimensions of the field would be handy to do vectorized, but they intersect. Here is a minimal, dummed down demo of the issue:
import numpy as np
from scipy.special import kn, iv
order = 20
size = 200
r = np.arange(-size, size, 1.0)
theta = np.arange(0, 2 * np.pi, 2 * np.pi / size)
z = np.arange(0, 1, 1 / size)
R, T, Z = np.meshgrid(r, theta, z)
m = np.arange(order)
eta = m * R
bessel_term = iv(m, eta) * kn(m, eta)
print(sum(bessel_term))
What would need to happen here is, that of each value of R the bessel functions with the order-vector would need to be calculated. It would need to create an additional dimension in the vector, so to speak.
Instead I get the following error message which I understand, but want to avoid:
Traceback (most recent call last):
eta = m * R
ValueError: operands could not be broadcast together with shapes (20,) (200,400,200)
How can I do that? One way would be to iterate through the space dimensions manually (not vectorized) and use the vectorisation for the inner loop where the sums over the bessel functions are calculated. But to me it would be preferable if there was a way to indicate to numpy that at this point an additional dimention in the involved vectors is needed and collapse this dimension again when the infinite sum is calculated.
Background, the math:
In reality I want to calulate a magnetic field with three vector components in every point in space of my meshgrid. In every point I need to calculate an infinit sum of besselfunctions of infinit orders:
$$ B_r = j\frac{3 \mu_0 I_1}{4 \pi r} \sum_{m= ... -5, -2, 1, 4, 7...}^{\infty} 2 m I_m(\eta_m) K_m\left( \eta_m \frac{r}{a’} \right) + \frac{2 \pi r m q}{p} \left[ I_{m-1}(\eta_m) K_{m-1}\left(\eta_m \frac{r}{a’} \right) + I_{m+1}(\eta_m) K_{m+1}\left( \eta_m \frac{r}{a'} \right) \right] \exp[jm(\theta - 2 \pi z/p)]$$
with $r$, $z$ and $\theta$ as input parameters (space coordinates), $a',p, q$ as constant arguments and $\eta_m=|mq|$ as an abreviation.
