# Help with modified bessel functions

I'm trying to parse the following expression (Eq. 35 in The Magnetic Field in the Vicinity of Parallel and Twisted Three-Wire Cable Carrying Balanced Three-Phased Current) in Python and calculate its sum over indices (...-5, -2, 1, 4, 7...). Relations between Bessel functions explains that these are the modified Bessel functions of the first and second kind.

Question: For $$K_m$$ I plan to use scipy.special.kn(n, x) but for $$I_m$$ there does not seem to be a method specific for integer orders. Should I just use scipy.special.iv(v, z) for real orders with $$v=m$$?

Here is the expression inside the Sum over $$m = (...-5, -2, 1, 4, 7...)$$ from Equation 35 from the linked paper:

$$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 and likely to be arrays, $$a',p, q$$ as constant arguments and $$\eta_m=|mq|$$ as an abreviation.

• Are the arguments of the modified Bessel functions real or complex? Are there specific concerns about accuracy or performance? Dec 7 '20 at 20:50
• As far as i currently understand, the arguments to the functions are real. Only the currents are complex (but they are not in this specific equation, and I will multiply with them afterwards). Regarding accuracy and performance: no, no special concerns. That said, I do try to use vectorized computations when working with numpy/scypy. And I am wrecking my brain right now how to do this one and the other two components of the magnetic field all vectorized... But that is a secondary concern for an other question. The basic functionality is most important right now. Dec 7 '20 at 21:28
• Is there anything stopping you from scripting this in Python using kn() and iv(), setting up say a 100 x 100 meshgrid for (r, z) and summing over say -20 \le m \le 19 and just seeing how well it works? Is there any reason to suspect that iv() will give the wrong answer with integer orders? In other words is "it's okay to do this" the answer you need, or are there other issues you'd like an answer to address?
– uhoh
Dec 8 '20 at 9:00
• I am doing this, i feel uncertain and I will see what happens. Dec 8 '20 at 11:53
• It's always okay to post a helpful answer to your own question in Stack Exchange; looking forward to seeing some interesting fields! :-)
– uhoh
Dec 9 '20 at 1:10

This is the code I came up with. I still dont feel confident about the results, mostly because I didn't find the usual amount of errors yet.

import numpy as np
from scipy.special import kn, iv
import matplotlib.pyplot as plt

global size, mu_0, p, a, q, current, order
mu_0 = 4 * np.pi * 1e-7
size = 10
order = 5  # order of bessel functions considered
p = 1.8  # length of cable [m] after which we have 360 turn; "Ganghöhe der Leiterwindung"
a = 2.67e-2  # helix radius [m] = 2*a/sqrt(3)
q = 2 * np.pi * a / p
current = 50 + 0j  # peak current [A] of I * sin(omega * t)

def coordinate_transform_cy2ca(r, theta, z):
x = np.cos(theta) * r
y = np.sin(theta) * r
return x, y, z

def vector_transform_cy2ca(v_r, v_theta, v_z, R, T, Z):
x = v_r * np.cos(T) - R * np.sin(T) * v_theta
y = v_r * np.sin(T) + R * np.cos(T) * v_theta
return x, y, v_z

def calc_b_r(R, T, Z, m):
m_1 = m
m_0 = m - 1
m_2 = m + 1

eta_m = np.abs(m_1 * q)
eta_r_a = eta_m * R / a
exp_term = np.exp((0 + 1j) * m_1 * (T - 2 * np.pi * Z / p))
m_0_term = iv(m_0, eta_m) * kn(m_0, eta_r_a)
m_2_term = iv(m_2, eta_m) * kn(m_2, eta_r_a)
m_1_term = 2 * m_1 * iv(m_1, eta_m) * kn(m_1, eta_r_a)
exp_product = 2 * np.pi * R * m_1 * q / p * exp_term
m_02_sum = m_0_term + m_2_term
bessel_term = m_1_term + exp_product * m_02_sum
b_r = ((0 + 1j) * 3 * mu_0 * current / (4 * np.pi * R)) * bessel_term
return b_r

def calc_b_theta(R, T, Z, m):
m_1 = m
m_0 = m - 1
m_2 = m + 1

eta_m = np.abs(m_1 * q)
eta_r_a = eta_m * R / a
exp_term = np.exp((0 + 1j) * m_1 * (T - 2 * np.pi * Z / p))

m_0_term = iv(m_0, eta_m) * kn(m_0, eta_r_a)
m_2_term = iv(m_2, eta_m) * kn(m_2, eta_r_a)
first_term = (2 * np.pi * m_1 * q / p) * (m_0_term - m_2_term)

mixed_term = (2 * eta_m / a) * iv(m_1, eta_m) * (kn(m_0, eta_r_a) + m_1 * kn(m_1, eta_r_a) / eta_r_a)
second_term = mixed_term * exp_term

bessel_term = second_term - first_term
b_theta = (3 * mu_0 * current / (4 * np.pi)) * bessel_term
return b_theta

def calc_b_z(R, T, Z, m):
m_1 = m
m_0 = m - 1
m__1 = m - 2  # the subscript "_1" is supposed to indicate a "-1"  -  its the K_ m-2  in the paper
m_2 = m + 1

eta_m = np.abs(m_1 * q)
eta_r_a = eta_m * R / a
exp_term = np.exp((0 + 1j) * m_1 * (T - 2 * np.pi * Z / p))

m__1_term = (eta_m / a) * iv(m_0, eta_m) * (m_1 * kn(m_0, eta_r_a) / eta_r_a + kn(m__1, eta_r_a))
# same term as before, just two orders up
mirror_term = (eta_m / a) * iv(m_2, eta_m) * (m_1 * kn(m_2, eta_r_a) / eta_r_a + kn(m_1, eta_r_a))

m_0_term = iv(m_0, eta_m) * kn(m_0, eta_r_a) / R
# same as before, just two orders up
m_2_term = iv(m_2, eta_m) * kn(m_2, eta_r_a) / R

m_02_term = (iv(m_0, eta_m) * kn(m_0, eta_r_a) - iv(m_2, eta_m) * kn(m_2, eta_r_a)) * m_1 / R

bessel_term = m_0_term + m__1_term - m__1_term - mirror_term - m_02_term
b_z = 3 * mu_0 * q * current / (4 * np.pi) * bessel_term * exp_term
return b_z

def main():
r = np.arange(1e-1, 2, (2-1e-1)/size)
theta = np.arange(0, 2 * np.pi/2, 2 * np.pi / (size*2))
z = np.arange(0, p, p / size)

R, T, Z = np.meshgrid(r, theta, z)

m = np.arange(int(-order / 3) * 3 - 2, order, 3)
m_num = len(m)

b_r = np.zeros([m_num, R.shape, R.shape, R.shape]).astype(np.complex_)
b_theta = np.zeros(b_r.shape).astype(np.complex_)
b_z = np.zeros(b_r.shape).astype(np.complex_)

for i, m in enumerate(m):
b_r[i] = calc_b_r(R, T, Z, m)
b_theta[i] = calc_b_theta(R, T, Z, m)
b_z[i] = calc_b_z(R, T, Z, m)

b_r = np.sum(b_r, axis=0)
b_theta = np.sum(b_theta, axis=0)
b_z = np.sum(b_z, axis=0)

b_x, b_y, b_z = vector_transform_cy2ca(b_r, b_theta, b_z, R, T, Z)
X, Y, Z = coordinate_transform_cy2ca(R, T, Z)

fig = plt.figure()
ax = fig.gca(projection='3d')
ax.quiver(X, Y, Z, b_x, b_y, b_z, length=0.1, normalize=True)

plt.show()

if __name__ == '__main__':
main()


Here is a plot of a cross-section of the field in 3D: To answer my direct question: yes, it is possible to use the modified bessel function for orders in real numbers even with natural, whole numbers.

There were two major challanges besides the bessel functions:

• Matplotlib does not permit plots in cylinder coordinates, so I had to do the coordinate and vector transformations myself. Aparently most people dont need to do vector transformations.
• The vectorisation was not as flexible as I would have liked. I ended up looping over the orders of the bessel equations in order to be able to use the vectorisation for the spatial dimensions.

Still missing:

• Color coding the strength of the magnetic field into the arrows of the quiver field.
• plot the helix of the three conductors
• Congratulations and a happy ending, but in fact it's only the beginning. Now that you have working code that already does what you want, you are ready for Code Review SE. Here are some examples of help I've gotten there: (1, 2, 3). They sometimes complain that my PEP-8 is poor, but for interesting problems and nicely formatted script like yours I think you can get a lot of improvement.
– uhoh
Dec 23 '20 at 13:06
• I think you can almost post this answer as-is as a question there, since you've already described areas you thing might need improvement.
– uhoh
Dec 23 '20 at 13:08