# Poincare map for Arnold-Beltrami-Childress Magnetic Field in Python

I want to plot the Poincare map for Arnold-Beltrami-Childress magnetic field for parameters $A=1, B=0.816, C=0.5773$ in Python for the Poincare section $z=0$.

Also, I am not able to understand what initial conditions should I vary for $x$, $y$, and $z$ so as to plot the entire map? I am able to get only a subset of the plot for $x=0$, $y=0$, and $z=0$ but how should I vary $x$, $y$, and $z$?

ODEs that I have used for the same:

$$\frac{d\mathbf{r}}{dt}=\frac{\bf{B}}{|\mathbf{B}|}$$

where $\mathbf{r}=(x,y,z)$ and $\mathbf{B}=\left( \begin{array}{c} A\sin(z)+C\cos(y)\\ B\sin(x)+A\cos(z)\\ C\sin(y)+B\cos(x) \end{array} \right)$.

Can you please provide me with the Python code?

• Would Computational Science be a better home for this question? Jul 7 '18 at 16:39
• Yes, it belongs there
– Nakul Aggarwal
Jul 7 '18 at 22:20

One first observation is that as $$x,y,z$$ are angles, one should take their values mod $$2\pi$$. Without that, the trajectory of the initial point $$[0,0,0]$$ exhibits this interesting pattern

def ABC_ode(u,t):
A, B, C = 1, 0.816, 0.5773
x, y, z = u
B = np.array([A*np.sin(z)+C*np.cos(y), B*np.sin(x)+A*np.cos(z), C*np.sin(y)+B*np.cos(x)])
return B/(1e-12+sum(B**2))**0.5

def mysolver(u0, tspan): return odeint(ABC_ode, u0, tspan, atol=1e-10, rtol=1e-11) So one has to reduce the coordinates by removing multiples of $$2\pi$$, but then also exclude the so-introduced "long jumps" by about $$2\pi$$ from the root computation.

To fill the cross section at $$z=0$$ maximally in a reasonable time, that is without computing similar trajectories, divide the square into smaller squares and remember for each small square in the grid if it was already visited. Only start trajectories from squares without visit. Here linear interpolation is used to find the intersection points, one could also use cubic interpolation with the derivative values or an additional integration step with smaller step size to get a more accurate intersection location. But at the level of the plot that is not necessary.

u0 = [0,0,0]
t = np.arange(0, 1500, 0.2)

N=54
grid = np.zeros([N,N], dtype=int)
for i in range(N):
for j in range(N):
if grid[i,j]>0: continue;
x0, y0 = (2*i+1)*pi/N-pi, (2*j+1)*pi/N-pi
u0 = [x0,y0,0]
px = []
py = []

u = mysolver(u0, t); u0 = u[-1]
u = np.mod(u+pi,2*pi)-pi
x,y,z = u.T

for k in range(len(z)-1):
if z[k]<=0 and z[k+1]>=0 and z[k+1]-z[k]<pi:
s = -z[k]/(z[k+1]-z[k])  # 0 = z[k] + s*(z[k+1]-z[k])
rx, ry = (1-s)*x[k]+s*x[k+1], (1-s)*y[k]+s*y[k+1]
px.append(rx);
py.append(ry);
m, n = int((rx+pi)*N/(2*pi)), int((ry+pi)*N/(2*pi))
grid[m,n]=1

#plt.plot(px, py, '-k', lw=0.2)
plt.plot(px, py, '.', ms=3)
plt.show() 