# Conversion of 1D plot in 3D plot gives different physics

I am to track the path of a particle due to Lorentz force of the planet magnetic field, assuming there is no gravitational force on it. The equation of motion due to Lorentz force in rotating spherical polar coordinates, in three components are :

Radial components $$\frac{d^2 r}{dt^2} = r \left( \frac{d\theta}{dt}\right)^2 + r\left(\frac{d\phi}{dt} + \Omega\right)^2 \sin^2 \theta + \\ + \beta r\frac{d\phi}{dt} \sin\theta B_\theta$$ \

latitudinal/meridional components $$\frac{d^2 \phi}{dt^2} = \frac{1}{r} \left[-2 \frac{dr}{dt} \frac{d\theta}{dt} - r(t) \left(\frac{d\phi}{dt} + \Omega\right)^2 \sin \theta \cos \theta(t) + \\ + \beta r \frac{d\phi}{dt} \sin \theta) B_r \right]$$\

Azimuthal component $$\frac{d^2 \phi}{dt^2} = \frac{1}{r\sin\theta} \left[-2 \left(\frac{d\phi}{dt} + \Omega\right) (\frac{dr}{dt}\sin\theta + r\frac{d\phi}{dt}\cos\theta)+ \\ + q \left(r \frac{d\phi}{dt} \sin \theta B_\theta - r\frac{d\phi}{dt} B_r\right)\right] \\\\$$ $$\Omega = 0.0017 \ rad/s, \ \beta= q/m$$

The aim is to learn solving coupled ODES numerically through scipy and plots the path of the particle with time. The solutions for theta (co-latitude angle, in radian),is = 1.57079633 1.57079633 1.57079633 1.57079633 1.57079633 1.57079633...in time. While the radial (r) linearly increases and azimuthal angle linearly decreases This would mean that the particle is moving aways with the same co-latidude, meaning in the same plane. This too can be seen in the 1D graph as shown below:

However, the 3D plot shows that is is spirally moving upward/downward depending on the nature of the charge of particle. This would mean that the particle in spirally away with increasing latitude/co-latidude. I failed to understand,the difference. Is it somthing to do with the coding or there is something which I dont understand in the interpretation of the solutions ( 2D and 3D plots). Thank you. The code is as given below:

from scipy.integrate import solve_ivp
import numpy as np
from math import sin, cos,pi

# charge by mass ratio
r_ho = 1e3               # 1gm/cm^3 Headmen 2021= 10^3 kg/m^3,
b = 1e-9                 # 1nm = 1e-3 micro
V = 10                   # Volt
e_p = 8.85e-12           # Farad/metre, free space permittivity (ε)
q_m = (3*e_p*V)/(rho*b**2)                # q/m = -3.46 x 10^3 C/kg

def odes(t,p):
## assigning each ODE to a vector element
r,x,theta,y,phi,z = p

# constants
m_u = 379312077e8           # m^3/s^2, G= 6.67408e-11m^3/kg*s^2 M = 5.68e26 # m_u =GM
R = 60268e3                 # metre
j_2 = 1.629071e-2
g_10 = 21141e-9                      #
w = 17e-5                   # rad/s = 9.74e-7, angular speed
nu_0 = 4*pi*1e-7            # magnetic permeability
B_theta = nu_0*(R/r)**3*g_10*sin(theta)
B_r = nu_0*2*(R/r)**3*g_10*cos(theta)
q_m = (3*e_p*V)/(rho*b**2)

# defining the ODEs
drdt = x
dxdt = r*(y**2+(z+w)**2*sin(theta)**2-q_m*z*sin(theta)*B_theta)
d(theta)dt = y
dydt = (-2*x*y + r*(z+w)**2*sin(theta)*cos(theta)+ q_m*r*z*sin(theta)*B_r)/r
d(phi)dt = z
dzdt = (-2*(z+w)*(x*sin(theta)+r*y*cos(θ)) + q_m*(x*B_theta-r*y*B_r))/(r*sin(theta))

return np.array([drdt,dxdt,d(theta)dt,dydt,d(phi)dt,dzdt])

## time window
t_span = (0,86400)
t = np.linspace(t_span[0],t_span[1],400)

#initial conditions
p0 = np.array([1.11*R,0.0,90.0*(pi/180),0.0*(pi/180), 0.0*(pi/180),0.0206*(pi/180)])

#  Solve IVP
sol = solve_ivp(odes,t_span ,p0,t_eval= t,method= "DOP853",dense_output=False, vectorised= False, atol=1e-6,rtol=1e-8)

import matplotlib.pyplot as plt
fig,ax=plt.subplots(2,3,figsize=(12,6),layout='constrained')
plt.suptitle('Equatorial launch micro sized charged grain',fontsize=13)

plt.subplot(2,3,1)
plt.plot(sol.t, sol.y[0,:]);plt.grid()
plt.xlabel('time in sec',)
plt.ylabel('r in m')
plt.subplot(2,3,2)
plt.plot(sol.t, sol.y[1,:]);plt.grid()
plt.xlabel('time in sec', fontsize=10)
plt.ylabel('velcoity in m',fontsize=10)
plt.subplot(2,3,3)
plt.plot(sol.t, sol.y[2,:]);plt.grid()
plt.xlabel('time in sec',fontsize=10)
plt.subplot(2,3,4)
plt.plot(sol.t, sol.y[3,:]);plt.grid()
plt.xlabel('time in sec',fontsize=10)
plt.subplot(2,3,5)
plt.plot(sol.t, sol.y[4,:]);plt.grid()
plt.xlabel('time in sec',fontsize=10)

plt.subplot(2,3,6)
plt.plot(sol.y[0,:], sol.y[2,:]);plt.grid()
plt.ylabel('theta',fontsize=10)
plt.tight_layout()
plt.show()

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

# 3D Cartesian
fig = plt.figure()

# spherical to cartesian
x1= sol.y[0,:] * np.sin(sol.y[2,:])* np.cos(sol.y[4,:])
y1= sol.y[0,:] * np.sin(sol.y[2,:])* np.sin(sol.y[4,:])
z1= sol.y[0,:] * np.cos(sol.y[2,:])

plt.suptitle('Equatorial launch micro sized charged grain',fontsize=12)
plt.title('1e-9nm,10V, 1.11R,0.0,90.0,0, 0,0.0206,0',fontsize=10)
plt.show()


• There are neither 2D nor 3D plots in your question. Commented Apr 14, 2023 at 10:10
• Now there is one 3-D plot and six 1-D plots. Commented Apr 14, 2023 at 16:10
• 1) You're comparing 1 million kilometer radius and 100 nanometer change in z? That's not meaningful numerically. 2) You should first present the differential equations that you want to solve in MathJax and state your assumptions. Don't just paste code and expect others to work backwards to understand your physics. Otherwise it's easy to wonder if you just got someone else's code you don't understand and are asking how it works. 3) As posted your script has several problems, better fix it.
– uhoh
Commented Apr 15, 2023 at 0:54
• @LunthangPeter nah don't worry about it, and here in Stack Exchange we really don't do personal information or outreach, unless someone says so explicitly in their profile. Just keep asking good SE questions and feel free to ping me here (or in The Observatory chat room for this site) if you have more to discuss.
– uhoh
Commented Apr 15, 2023 at 12:27
• I don't see any inconsistency between the graphs. The top graphs show that the co-latitude is constant, That means that the particle is constrained to the surface of a cone. The "r" is increasing, so the particle moves away from the centre, and the azimuth is decreasing, so it is moving clockwise around the pole. So the motion should be a helical motion outwards from the centre on a cone, which is what is shown in the 3d plot. The physics is the same. I don't know if the physics is right, but the same things seem to be shown in the two plots. Commented Apr 16, 2023 at 9:42