<p>I am trying to plot the path of a charged particle in a planet's magnetic field. For positive and negative charge (β=charge/mass) different solutions/paths are expected. But,I got the same solution (values). In the plots of the variables with time in 2D it too show the same graphs. However, plots in 3D clearly shows different solution/trajectory. Please somehow help me review the code.</p>
<pre><code>import numpy as np
import matplotlib.pyplot as plt
from math import sin, cos, pi
from scipy.integrate import odeint
scales = np.array([1e2, 0.1, 1, 1e-15, 10, 0.1])

β =  9.67e7   # charge/mass
def odes(p, t):
# assigning each ODE to a vector element
r,x,θ,y,ϕ,z = p*scales
    
# constants
R = 60268e3;     g_10 = 21141e-9;     Ω = 9.74e-3
B_θ = (R/r)**3*g_10*sin(θ)
B_r = 2*(R/r)**3*g_10*cos(θ)
β = + 9.67e5

# defining the ODEs only Lorentz Force
drdt = x
dxdt = r*(y**2 +(z+Ω)**2*sin(θ)**2-β*z*sin(θ)*B_θ)
dθdt = y
dydt = (-2*x*y +r*(z+Ω)**2*sin(θ)*cos(θ)+β*r*z*sin(θ)*B_r)/r 
dϕdt = z
dzdt = (-2*(z+Ω)*(x*sin(θ)+r*y*cos(θ))+β*(x*B_θ-r*y*B_r))/(r*sin(θ))

return np.array([drdt,dxdt,dθdt,dydt,dϕdt,dzdt])/scales

# initial conditions
r0 = 6.8e+07
x0 = 0.002
θ0 = 80.0*1.745e-2; y0 = 0.0*1.745e-2;  ϕ0 = 0.0*1.745e-2; z0 = 0.0202*1.745e-2

# time window
t = np.linspace(0,3600,4000)
p0 = np.array([r0,x0,θ0,y0,ϕ0,z0])
p = odeint(odes,p0,t, atol=1e-8, rtol=1e-8)

r,x,θ,y,ϕ,z = p.T*scales[:,None]
print (p.T)

#plot the results
fig,ax=plt.subplots(2,3,figsize=(8,4))
plt.xlabel('time in seconds')
plt.ylabel('parameters')
#plt.title('EM for Negative Charge ( β = - 9.67e5) ')
plt.title('EM for Positive Charge ( β = + 9.67e5) ')
for a,u in zip(ax.flatten(),[r,x,θ,y,ϕ,z]):
    a.plot(t,u); a.grid()
plt.tight_layout(); plt.show()


import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import axes3d
import matplotlib.pyplot as plt
from mpl_toolkits import mplot3d

# 3D Cartesian
fig = plt.figure()

# spherical to cartesian
x = (r * np.sin(θ)* np.cos(ϕ))/1000
y = (r * np.sin(θ)* np.sin(ϕ))/1000
z = (r * np.cos(θ))/1000
print (&quot;x =&quot;, x)
print(&quot;y =&quot;, y)
print (&quot;z =&quot;, z)

ax = plt.axes(projection=&quot;3d&quot;)
ax.plot3D(x,y,z, color=&quot;red&quot;,linewidth='0.5')
ax.set_xlabel('X-Axis')
ax.set_ylabel('Y-axis')
ax.set_zlabel('Z-axis')

#plt.title('EM for Negative Charge ( β = - 9.67e5)')
plt.title('EM for Positive Charge ( β = + 9.67e5) ')
plt.show()

 **Solutions for Positive beta ( +q/m):**
 [[ 6.80000000e+07  6.80025533e+07  6.80102123e+07 ...  2.35930802e+09
   2.35989788e+09  2.36048774e+09]
   [ 2.00000000e-03  5.67230607e+06  1.13433370e+07 ...  6.55236454e+08
   6.55236597e+08  6.55236740e+08]
 [ 1.39600000e+00  1.39600663e+00  1.39602652e+00 ...  1.56554064e+00
   1.56554189e+00  1.56554314e+00]
 [ 0.00000000e+00  1.47303014e+10  2.94505438e+10 ...  1.38803448e+09
   1.38734094e+09  1.38664793e+09]
 [ 0.00000000e+00  3.15050355e-06  6.16484138e-06 ... -3.35040490e+00
  -3.35128101e+00 -3.35215712e+00]
 [ 3.52490000e-04  3.44925730e-04  3.22240365e-04 ... -9.73211386e-02
  -9.73211780e-02 -9.73212174e-02]]
x = [ 6.69638146e+06  6.69664074e+06  6.69741848e+06 ... -1.16679833e+08
 -1.18501460e+08 -1.20314886e+08]
y = [ 0.00000000e+00  2.10977904e+02  4.12885226e+02 ... -2.05055169e+08
 -2.04076079e+08 -2.03080807e+08]
z = [1182571.48669918 1182571.48610586 1182571.4845813  ... 1239972.95449587
 1239988.16352907 1240003.37257228]
</code></pre>
<p><a href="https://i.sstatic.net/zF1Rm.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/zF1Rm.png" alt="2D plots for Positive charge" /></a>
<a href="https://i.sstatic.net/a1RQq.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/a1RQq.png" alt="3D plots for Positive" /></a></p>
<pre><code>****Solution for Negative charge****

[[ 6.80000000e+07  6.80025533e+07  6.80102124e+07 ...  2.35332758e+09
   2.35391578e+09  2.35450398e+09]
 [ 2.00000000e-03  5.67236481e+06  1.13434496e+07 ...  6.53392994e+08
   6.53393123e+08  6.53393253e+08]
 [ 1.39600000e+00  1.39600663e+00  1.39602652e+00 ...  1.56601504e+00
   1.56601630e+00  1.56601755e+00]
 [ 0.00000000e+00  1.47299932e+10  2.94499349e+10 ...  1.39926205e+09
   1.39856314e+09  1.39786474e+09]
 [ 0.00000000e+00  3.15047326e-06  6.16459918e-06 ... -3.35041308e+00
  -3.35128919e+00 -3.35216529e+00]
 [ 3.52490000e-04  3.44915637e-04  3.22200020e-04 ... -9.73209613e-02
  -9.73210008e-02 -9.73210403e-02]]
x = [ 6.69638146e+06  6.69664074e+06  6.69741849e+06 ... -1.16401081e+08
 -1.18218003e+08 -1.20026744e+08]
y = [ 0.00000000e+00  2.10975876e+02  4.12869006e+02 ... -2.04526354e+08
 -2.03549583e+08 -2.02556671e+08]
z = [1182571.48669918 1182571.48749931 1182571.49009402 ... 1125190.09460013
 1125174.89571226 1125159.69681527]}
</code></pre>
<p><a href="https://i.sstatic.net/pr0dM.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/pr0dM.png" alt="P2D plots for Negativelots for Negataive charge" /></a></p>

CC BY-SA 4.0
2022-12-15T23:50:31.867Z
At the moment it's tough to say what's wrong. But a general advice to find bugs in your code is always to first analyze the simplest system possible. You have not chosen the simplest system possible. Your magnetic field is a dipole. Try a constant magnetic field instead, aligned with the z-axis (or any other), with v_z = 0 along that axis. That way you disentangle dimensions, and hence disentangle equations.
AtmosphericPrisonEscape
CC BY-SA 4.0
2022-12-16T10:11:05.313Z
Can you mention why you have BOTH cartesian AND spherical coordinates in your integration loop? There are only three degrees of freedom; why are you integrating six? Actually can you either add the equations you are integrating [using MathJax](https://math.meta.stackexchange.com/q/5020/284619) or just add a link to them? *Thanks!*
uhoh
CC BY-SA 4.0
2022-12-17T07:21:16.523Z
@AtmosphericPrisonEscape, Thank you all for your valuable suggestions. uhoh, the coordinates transformation are done to make the plot easier.
Lunthang Peter
<astrophysics><gravity><planetary-ring><dust><planetary-systems>
The path of a particle in planet's magnetic field doesn't not seem to change with the charge of the particle. Can someone please check what is wrong?
Lunthang Peter
# The path of a particle in planet's magnetic field doesn't not seem to change with the charge of the particle. Can someone please check what is wrong?

(Asked by: username_0 on 2022-12-15)

I am trying to plot the path of a charged particle in a planet's magnetic field. For positive and negative charge (β=charge/mass) different solutions/paths are expected. But,I got the same solution (values). In the plots of the variables with time in 2D it too show the same graphs. However, plots in 3D clearly shows different solution/trajectory. Please somehow help me review the code.
import numpy as np
import matplotlib.pyplot as plt
from math import sin, cos, pi
from scipy.integrate import odeint
scales = np.array([1e2, 0.1, 1, 1e-15, 10, 0.1])

β =  9.67e7   # charge/mass
def odes(p, t):
# assigning each ODE to a vector element
r,x,θ,y,ϕ,z = p*scales
    
# constants
R = 60268e3;     g_10 = 21141e-9;     Ω = 9.74e-3
B_θ = (R/r)**3*g_10*sin(θ)
B_r = 2*(R/r)**3*g_10*cos(θ)
β = + 9.67e5

# defining the ODEs only Lorentz Force
drdt = x
dxdt = r*(y**2 +(z+Ω)**2*sin(θ)**2-β*z*sin(θ)*B_θ)
dθdt = y
dydt = (-2*x*y +r*(z+Ω)**2*sin(θ)*cos(θ)+β*r*z*sin(θ)*B_r)/r 
dϕdt = z
dzdt = (-2*(z+Ω)*(x*sin(θ)+r*y*cos(θ))+β*(x*B_θ-r*y*B_r))/(r*sin(θ))

return np.array([drdt,dxdt,dθdt,dydt,dϕdt,dzdt])/scales

# initial conditions
r0 = 6.8e+07
x0 = 0.002
θ0 = 80.0*1.745e-2; y0 = 0.0*1.745e-2;  ϕ0 = 0.0*1.745e-2; z0 = 0.0202*1.745e-2

# time window
t = np.linspace(0,3600,4000)
p0 = np.array([r0,x0,θ0,y0,ϕ0,z0])
p = odeint(odes,p0,t, atol=1e-8, rtol=1e-8)

r,x,θ,y,ϕ,z = p.T*scales[:,None]
print (p.T)

#plot the results
fig,ax=plt.subplots(2,3,figsize=(8,4))
plt.xlabel('time in seconds')
plt.ylabel('parameters')
#plt.title('EM for Negative Charge ( β = - 9.67e5) ')
plt.title('EM for Positive Charge ( β = + 9.67e5) ')
for a,u in zip(ax.flatten(),[r,x,θ,y,ϕ,z]):
    a.plot(t,u); a.grid()
plt.tight_layout(); plt.show()


import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import axes3d
import matplotlib.pyplot as plt
from mpl_toolkits import mplot3d

# 3D Cartesian
fig = plt.figure()

# spherical to cartesian
x = (r * np.sin(θ)* np.cos(ϕ))/1000
y = (r * np.sin(θ)* np.sin(ϕ))/1000
z = (r * np.cos(θ))/1000
print ("x =", x)
print("y =", y)
print ("z =", z)

ax = plt.axes(projection="3d")
ax.plot3D(x,y,z, color="red",linewidth='0.5')
ax.set_xlabel('X-Axis')
ax.set_ylabel('Y-axis')
ax.set_zlabel('Z-axis')

#plt.title('EM for Negative Charge ( β = - 9.67e5)')
plt.title('EM for Positive Charge ( β = + 9.67e5) ')
plt.show()

 **Solutions for Positive beta ( +q/m):**
 [[ 6.80000000e+07  6.80025533e+07  6.80102123e+07 ...  2.35930802e+09
   2.35989788e+09  2.36048774e+09]
   [ 2.00000000e-03  5.67230607e+06  1.13433370e+07 ...  6.55236454e+08
   6.55236597e+08  6.55236740e+08]
 [ 1.39600000e+00  1.39600663e+00  1.39602652e+00 ...  1.56554064e+00
   1.56554189e+00  1.56554314e+00]
 [ 0.00000000e+00  1.47303014e+10  2.94505438e+10 ...  1.38803448e+09
   1.38734094e+09  1.38664793e+09]
 [ 0.00000000e+00  3.15050355e-06  6.16484138e-06 ... -3.35040490e+00
  -3.35128101e+00 -3.35215712e+00]
 [ 3.52490000e-04  3.44925730e-04  3.22240365e-04 ... -9.73211386e-02
  -9.73211780e-02 -9.73212174e-02]]
x = [ 6.69638146e+06  6.69664074e+06  6.69741848e+06 ... -1.16679833e+08
 -1.18501460e+08 -1.20314886e+08]
y = [ 0.00000000e+00  2.10977904e+02  4.12885226e+02 ... -2.05055169e+08
 -2.04076079e+08 -2.03080807e+08]
z = [1182571.48669918 1182571.48610586 1182571.4845813  ... 1239972.95449587
 1239988.16352907 1240003.37257228]



****Solution for Negative charge****

[[ 6.80000000e+07  6.80025533e+07  6.80102124e+07 ...  2.35332758e+09
   2.35391578e+09  2.35450398e+09]
 [ 2.00000000e-03  5.67236481e+06  1.13434496e+07 ...  6.53392994e+08
   6.53393123e+08  6.53393253e+08]
 [ 1.39600000e+00  1.39600663e+00  1.39602652e+00 ...  1.56601504e+00
   1.56601630e+00  1.56601755e+00]
 [ 0.00000000e+00  1.47299932e+10  2.94499349e+10 ...  1.39926205e+09
   1.39856314e+09  1.39786474e+09]
 [ 0.00000000e+00  3.15047326e-06  6.16459918e-06 ... -3.35041308e+00
  -3.35128919e+00 -3.35216529e+00]
 [ 3.52490000e-04  3.44915637e-04  3.22200020e-04 ... -9.73209613e-02
  -9.73210008e-02 -9.73210403e-02]]
x = [ 6.69638146e+06  6.69664074e+06  6.69741849e+06 ... -1.16401081e+08
 -1.18218003e+08 -1.20026744e+08]
y = [ 0.00000000e+00  2.10975876e+02  4.12869006e+02 ... -2.04526354e+08
 -2.03549583e+08 -2.02556671e+08]
z = [1182571.48669918 1182571.48749931 1182571.49009402 ... 1125190.09460013
 1125174.89571226 1125159.69681527]}

**Tags:** astrophysics, gravity, planetary-ring, dust, planetary-systems

> **Comments:**
> **username_1**: At the moment it's tough to say what's wrong. But a general advice to find bugs in your code is always to first analyze the simplest system possible. You have not chosen the simplest system possible. Your magnetic field is a dipole. Try a constant magnetic field instead, aligned with the z-axis (or any other), with v_z = 0 along that axis. That way you disentangle dimensions, and hence disentangle equations.  
> **username_2**: Can you mention why you have BOTH cartesian AND spherical coordinates in your integration loop? There are only three degrees of freedom; why are you integrating six? Actually can you either add the equations you are integrating [using MathJax](https://math.meta.stackexchange.com/q/5020/284619) or just add a link to them? *Thanks!*  
> **username_0**: @username_1, Thank you all for your valuable suggestions. username_2, the coordinates transformation are done to make the plot easier.  

---

## 0 answers

