<p>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 :</p>
<p><strong>Radial components</strong>
<span class="math-container">$$\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 $$</span>  \</p>
<p><strong>latitudinal/meridional components</strong>
<span class="math-container">$$
\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] $$</span>\</p>
<p><strong>Azimuthal component</strong>
<span class="math-container">$$
\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] \\\\$$</span>
<span class="math-container">$\Omega = 0.0017 \ rad/s, \ \beta= q/m$</span></p>
<p>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:<a href="https://i.sstatic.net/9vKuO.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/9vKuO.png" alt="co-latitude vs time " /></a></p>
<p>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:</p>
<pre><code>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= &quot;DOP853&quot;,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
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.ylabel('Co-laitude in radian',fontsize=10)
plt.subplot(2,3,4)
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.subplot(2,3,5)
plt.plot(sol.t, sol.y[4,:]);plt.grid()
plt.xlabel('time in sec',fontsize=10)
plt.ylabel('Azimuthal in radian',fontsize=10)

plt.subplot(2,3,6)
plt.plot(sol.y[0,:], sol.y[2,:]);plt.grid()
plt.xlabel('radial',fontsize=10)
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()
</code></pre>
<p><a href="https://i.sstatic.net/sLtVR.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/sLtVR.png" alt="3D plot of r, theta, phi in Cartesian coordinates" /></a></p>

CC BY-SA 4.0
2023-04-14T10:10:40.597Z
There are neither 2D nor 3D plots in your question.
AtmosphericPrisonEscape
CC BY-SA 4.0
2023-04-14T15:34:28.380Z
@AtmosphericPrisonEscape, please recheck.
Lunthang Peter
CC BY-SA 4.0
2023-04-14T16:10:56.823Z
Now there is one 3-D plot and six 1-D plots.
AtmosphericPrisonEscape
CC BY-SA 4.0
2023-04-14T22:41:09.300Z
@uhoh, You the man. Thank you so much. One thing though, the question was related to theta (co-latitude), which remain constant (at about 1,575 rad ) with time . But the 3D plot show the path spiraling vertically means, there is an increase of altitude (co-latitude) vertically. The nature of trajectory should remain the same irrespective of the co-ordinate. I shall surely follow your suggestion to restart from simpler problem. However, I fell the problem is something to do with the Transformation of Co-ordinates to x, y,z in matplotlib rather the physics.
Lunthang Peter
CC BY-SA 4.0
2023-04-15T00:54:07.670Z
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](https://math.meta.stackexchange.com/q/5020/284619) 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
CC BY-SA 4.0
2023-04-15T10:42:45.530Z
@uhoh, Thank you once I again. I finally noticed that z is in a factor of e-17 ! Yes the code was copied and pasted from my spyder ide but not someone's else.  The word "scale" in your first comment now strike me. How do I mentioned you in case I want to  acknowledge in a paper ? Can get your 'linked' or 'Github' or place of work.
Lunthang Peter
CC BY-SA 4.0
2023-04-15T12:27:54.557Z
@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
CC BY-SA 4.0
2023-04-15T12:28:52.603Z
I am now learning numerical/ mathematical scaling and normalization. Any advice or comment  would be highly solicited.
Lunthang Peter
CC BY-SA 4.0
2023-04-16T09:42:20.923Z
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.
James K
CC BY-SA 4.0
2023-04-16T10:24:10.427Z
@JamesK, the helix is not constrained in a plane. It is vertically rising albeit at the rate e-7 which is not shown in the top right plot.
Lunthang Peter
CC BY-SA 4.0
2023-04-16T10:29:55.740Z
No, it is constrained to the surface of a cone, that is what constant latitude means.
James K
CC BY-SA 4.0
2023-04-16T14:35:01.600Z
@Jamesk, Thank you. I now understand. I was thinking in terms of Z-axis. Related to the plot, how do I draw a sphere, representing the spherical magnetic field.
Lunthang Peter
CC BY-SA 4.0
2023-04-16T14:45:05.633Z
Don't know, but with such uneven axis scales, a sphere wouldn't look like much (it would be visibly indistinguishable from a cylinder at that scale).  Also, can magnetic fields be spherical? Doesn't that break Gauss's law,  the zero divergence theorem?
James K
CC BY-SA 4.0
2023-04-16T20:00:51.297Z
@JamesK, I meant the source of the magnetic field. such as planet with radius, say, 60,268km. There is a question related to this: https://stackoverflow.com/questions/16158339/plotting-a-sphere-in-python-for-an-orbital-trajectory.
Lunthang Peter
CC BY-SA 4.0
2023-04-16T20:49:15.920Z
@JamesK,. You are right. The sphere actually looks like a cylinder. By scaling with : scales = np.array([1e2, 0.1, 1e-10, 1e-3, 1, 10]) I the sphere appears but the trajectories (red helix) dissapears.
Lunthang Peter
CC BY-SA 4.0
2023-04-16T21:01:46.957Z
are you surprised by that.  I'm not.   Fundamentally there is no problem in the plotting.  You need to go back and solve the underlying maths.  Uhoh says "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."  I agree.  I don't think at the moment there is an answerable question here.
James K
<astrophysics><magnetic-field><planetary-science><modeling><spherical-astronomy>
Conversion of 1D plot in 3D plot gives different physics
Lunthang Peter
# Conversion of 1D plot in 3D plot gives different physics

(Asked by: username_0 on 2023-04-14)

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)

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

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.ylabel('Co-laitude in radian',fontsize=10)
plt.subplot(2,3,4)
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.subplot(2,3,5)
plt.plot(sol.t, sol.y[4,:]);plt.grid()
plt.xlabel('time in sec',fontsize=10)
plt.ylabel('Azimuthal in radian',fontsize=10)

plt.subplot(2,3,6)
plt.plot(sol.y[0,:], sol.y[2,:]);plt.grid()
plt.xlabel('radial',fontsize=10)
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()

**Tags:** astrophysics, magnetic-field, planetary-science, modeling, spherical-astronomy
> **Comments:**
> **username_1**: There are neither 2D nor 3D plots in your question.  
> **username_0**: @AtmosphericPrisonEscape, please recheck.  
> **username_1**: Now there is one 3-D plot and six 1-D plots.  
> **username_0**: @uhoh, You the man. Thank you so much. One thing though, the question was related to theta (co-latitude), which remain constant (at about 1,575 rad ) with time . But the 3D plot show the path spiraling vertically means, there is an increase of altitude (co-latitude) vertically. The nature of trajectory should remain the same irrespective of the co-ordinate. I shall surely follow your suggestion to restart from simpler problem. However, I fell the problem is something to do with the Transformation of Co-ordinates to x, y,z in matplotlib rather the physics.  
> **username_2**: 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](https://math.meta.stackexchange.com/q/5020/284619) 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.  
> **username_0**: @uhoh, Thank you once I again. I finally noticed that z is in a factor of e-17 ! Yes the code was copied and pasted from my spyder ide but not someone's else.  The word "scale" in your first comment now strike me. How do I mentioned you in case I want to  acknowledge in a paper ? Can get your 'linked' or 'Github' or place of work.  
> **username_2**: @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.  
> **username_0**: I am now learning numerical/ mathematical scaling and normalization. Any advice or comment  would be highly solicited.  
> **username_3**: 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.  
> **username_0**: @JamesK, the helix is not constrained in a plane. It is vertically rising albeit at the rate e-7 which is not shown in the top right plot.  
> **username_3**: No, it is constrained to the surface of a cone, that is what constant latitude means.  
> **username_0**: @Jamesk, Thank you. I now understand. I was thinking in terms of Z-axis. Related to the plot, how do I draw a sphere, representing the spherical magnetic field.  
> **username_3**: Don't know, but with such uneven axis scales, a sphere wouldn't look like much (it would be visibly indistinguishable from a cylinder at that scale).  Also, can magnetic fields be spherical? Doesn't that break Gauss's law,  the zero divergence theorem?  
> **username_0**: @JamesK, I meant the source of the magnetic field. such as planet with radius, say, 60,268km. There is a question related to this: https://stackoverflow.com/questions/16158339/plotting-a-sphere-in-python-for-an-orbital-trajectory.  
> **username_0**: @JamesK,. You are right. The sphere actually looks like a cylinder. By scaling with : scales = np.array([1e2, 0.1, 1e-10, 1e-3, 1, 10]) I the sphere appears but the trajectories (red helix) dissapears.  
> **username_3**: are you surprised by that.  I'm not.   Fundamentally there is no problem in the plotting.  You need to go back and solve the underlying maths.  Uhoh says "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."  I agree.  I don't think at the moment there is an answerable question here.  

---

## 0 answers

