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:co-latitude vs time

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)

# print('radial =',sol.y[0,:],'theta=', sol.y[2,:],'phi=',sol.y[4,:])

import matplotlib.pyplot as plt
plt.suptitle('Equatorial launch micro sized charged grain',fontsize=13)

plt.plot(sol.t, sol.y[0,:]);plt.grid()
plt.xlabel('time in sec',)
plt.ylabel('r in m')
plt.plot(sol.t, sol.y[1,:]);plt.grid()
plt.xlabel('time in sec', fontsize=10)
plt.ylabel('velcoity in m',fontsize=10)
plt.plot(sol.t, sol.y[2,:]);plt.grid()
plt.xlabel('time in sec',fontsize=10)
plt.ylabel('Co-laitude in radian',fontsize=10)
plt.plot(sol.t, sol.y[3,:]);plt.grid()
plt.xlabel('time in sec',fontsize=10)
plt.ylabel('Co-laitude speed in rad/s',fontsize=10)
plt.plot(sol.t, sol.y[4,:]);plt.grid()
plt.xlabel('time in sec',fontsize=10)
plt.ylabel('Azimuthal in radian',fontsize=10)

plt.plot(sol.y[0,:], sol.y[2,:]);plt.grid()

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)

3D plot of r, theta, phi in Cartesian coordinates

  • 1
    $\begingroup$ There are neither 2D nor 3D plots in your question. $\endgroup$ Commented Apr 14, 2023 at 10:10
  • 1
    $\begingroup$ Now there is one 3-D plot and six 1-D plots. $\endgroup$ Commented Apr 14, 2023 at 16:10
  • 3
    $\begingroup$ 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. $\endgroup$
    – uhoh
    Commented Apr 15, 2023 at 0:54
  • 1
    $\begingroup$ @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. $\endgroup$
    – uhoh
    Commented Apr 15, 2023 at 12:27
  • 1
    $\begingroup$ 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. $\endgroup$
    – James K
    Commented Apr 16, 2023 at 9:42


You must log in to answer this question.