Tracing particles (RKF45)¶
Introduction¶
(To be written)
You can use irfpy.util.rungekutta module.
Moon motion around the Earth¶
Moon motion around the Earth can be simulated here.
The equation to solve is
under certain initial conditions. Assume the Earth at the center.
The right-hand side is the force acting on the unit mass.
M is the mass of the earth that is 5.972e24 kg.
For the case of gravity, you can use irfpy.util.rungekutta.GravityForce class
to represent the force field.
>>> from irfpy.util import rungekutta
>>> gravity = rungekutta.GravityForce(5.972e24) # Earth gravity field
As the initial position of the Moon, let’s use 384000 km. The initial velocity is 1.022 km/s, which is a rough average speed of the Moon.
>>> import numpy as np
>>> x0 = np.array([3.84e8, 0, 0]) # Moon at 3.84e5 km = 3.84e8 m
>>> v0 = np.array([0, 1.022e3, 0]) # Moon with velocity 1.022 km/s (=1022 m/s)
The tolerance (error allowed) and time resoulution are difficult to determine. In this first tutorial, let’s just give the torrelance for the position as 1000 m and the velocity as 0.01 m/s. Then, try the time resolutin of 3600 s.
>>> rk = rungekutta.RKF45_EqM(x0, v0, gravity, 3600, tolx=1000, tolv=0.01)
Let’s see the time is propagated one step.
>>> rk.progress()
INFO:rkf45:T=0 : DT=3600 : MAXRELERR=1.4407569548939354
INFO:rkf45:T=0 : DT=1800.0 : MAXRELERR=0.7203994515470695
The INFO says that the time resolution 3600s does not fit the given tolerance. (MAXRELERR is the maximum relative error; if it exceed 1, the calculated error is more than the tolerance.) Therefore, the time resolution was automatically reduced to 1800s, in which case the calculation error is less than the given tolerance.
You can see the history of position and velocity of the Moon.
>>> rk.pos
[array([ 3.84000000e+08, 0.00000000e+00, 0.00000000e+00]),
array([ 3.83995623e+08, 1.83959301e+06, 0.00000000e+00])]
>>> rk.vel
[array([ 0., 1022., 0.]),
array([ -4.86372383, 1021.98834984, 0. ])]
>>> rk.t
[0, 1800.0]
Let’s proceed one step more.
>>> rk.progress()
INFO:rkf45:T=1800.0 : DT=1818.0 : MAXRELERR=0.727579533783719
The time resolution is 1818 sec now. To adjust the time resolution optimal, the time resolution is increased by 1% every step automatically. You can see the position, velocity and time.
>>> rk.pos
[array([ 3.84000000e+08, 0.00000000e+00, 0.00000000e+00]),
array([ 3.83995621e+08, 1.83959301e+06, 0.00000000e+00]),
array([ 3.83982309e+08, 3.69753922e+06, 0.00000000e+00])]
>>> rk.vel
[array([ 0., 1022., 0.]),
array([ -4.86535323, 1021.98834594, 0. ]),
array([ -9.77924484, 1021.95291671, 0. ])]
>>> rk.t
[0, 1800.0, 3618.0]
Now, let’s proceed 100 times!
>>> for i in range(100):
... rk.progress()
INFO:rkf45:T=3618.0 : DT=1836.18 : MAXRELERR=0.7348136115297675
INFO:rkf45:T=5454.18 : DT=1854.5418000000004 : MAXRELERR=0.7421013404829427
INFO:rkf45:T=7308.721800000001 : DT=1873.0872180000008 : MAXRELERR=0.7494423421751708
INFO:...
INFO:rkf45:T=74988.49608562831 : DT=2549.884960856282 : MAXRELERR=0.9991239715367556
INFO:rkf45:T=77538.3810464846 : DT=2575.3838104648516 : MAXRELERR=1.0076583436131477
INFO:rkf45:T=77538.3810464846 : DT=1287.6919052324258 : MAXRELERR=0.5041685480773449
INFO:rkf45:T=78826.07295171703 : DT=1300.568824284747 : MAXRELERR=0.5088352167010307
INFO:...
INFO:rkf45:T=192202.43517117843 : DT=2434.3324464793154 : MAXRELERR=0.8476346109211444
INFO:rkf45:T=194636.76761765775 : DT=2458.675770944119 : MAXRELERR=0.8529364400207996
You may realize that the time resolution is increased every step. At a certain point, the error exceeds the tolerance (here T=77538.38). As soon as it happens, the time resolution is automatically adjusted by half to keep the error below tolerance. You can see the latest position, velocity and time as follows.
>>> print(rk)
t=197095.44338860188 x=[332705561.53518194 192383094.1471222 0.0] v=[-508.3868366850197 885.5967600160176 0.0]
ExB drift¶
ExB drift is the famous motion of charged particle. When the electric field is perpendicular to the magnetic field, the equation of motion can be written as
Here as example, we trace the proton (m=1.67e-27 kg, q=1.60e-19) under the uniform magnetic field (0, 0, 2e-9) Tesla (2 nano-tesla) and electric field (0, 5e-4, 0) V/m. Initially proton is at the origin at rest. In this case, the gyrofrequency for the proton is \(\omega=qB/m=0.192 s^{-1}\). Corresponding gyro period is \(2\pi/\omega=32.7\) s. The E-B drift velocity is \(\vec{E}\times\vec{B}/|\vec{B}|^2 = 2.5e5\) m/s. Considering this, the time range to simulate here is 0 to 60 sec, and the spatial domain is <2e8 [m].
>>> import numpy as np
>>> mp = 1.67e-27
>>> qe = 1.60e-19
>>> x0 = np.array([0, 0, 0])
>>> v0 = np.array([0, 0, 0])
>>> b = np.array([0, 0, 2e-9])
>>> e = np.array([0, 5e-4, 0])
To represent the field, you are expected to use irfpy.util.fields.Vector3dField.
This is a powerful way of representing vector field from gridded data, function, or other way.
In this tutorial, using irfpy.util.fields.UniformVector3dField is a good idea
to represent the uniform vector field.
Todo
Create a tutorial of irfpy.util.fields.
>>> from irfpy.util import fields
>>> efield = fields.UniformVector3dField(e[0], e[1], e[2])
>>> bfield = fields.UniformVector3dField(b[0], b[1], b[2])
Then, the uniform electromagnetic force field can be easily prepared using
irfpy.util.rungekutta.ElectroMagneticForce.
>>> from irfpy.util import rungekutta
>>> emforce = rungekutta.ElectroMagneticForce(mp, qe, efield, bfield)
Ok, force field is ready. It is now time to instance a particle.
>>> rk = rungekutta.RKF45_EqM(x0, v0, emforce, 1e-5, tolx=10, tolv=0.1, dtmax=6)
Here, we put the original time step as 1e-5 s. This could be very short, considering the particle motion is tracked over 60 sec. But let’s try to look at what happens.
>>> rk.progress()
INFO:rkf45:T=0 : DT=1e-05 : MAXRELERR=0.001875977849181698
INFO:rkf45:T=0 : DT=1.5000000000000002e-05 : MAXRELERR=0.002813966773766996
INFO:rkf45:T=0 : DT=2.2500000000000005e-05 : MAXRELERR=0.00422095016063162
INFO:...
INFO:rkf45:T=0 : DT=0.0008649755859375004 : MAXRELERR=0.1622675020312414
INFO:rkf45:T=0 : DT=0.0012974633789062506 : MAXRELERR=0.24340124948608377
>>> print(rk)
t=0.0012974633789062506 x=[3.341497167312416e-06 0.040321226609146936 0.0] v=[0.007726223062830551 62.1539336785754 0.0]
INFO says that the time resolution 1e-5 sec results in the error with 1.88% of the given tolerance. This is over qualified, so that a longer time solution is examined automatically, and the programme reached up to use 0.0013 s here (the error is 24.3% of the tolerance). You can save 130 steps!!
Let’s see the next step. You can progress again.
>>> rk.progress()
INFO:rkf45:T=0.0012974633789062506 : DT=0.001310438012695313 : MAXRELERR=0.2458352400867625
>>> print(rk)
t=0.002607901391601564 x=[2.713496498343062e-05 0.16290178507273637 0.0] v=[0.03121471330735068 124.92940278139895 0.0]
Particle moves further by 0.00131 sec. The time step was increased slightly, because the error relative to tolerance is 24% so that longer time step yet can satisfy the tolerance.
Ok, now we try to proceed the particle up to 60 seconds.
>>> while rk.t[-1] < 60.0:
... rk.progress()
INFO:rkf45:...
>>> print(rk)
t=60.0019835998448 x=[16144331.743594196 677113.5717036109 0.0] v=[129746.31314081445 -219178.1257260302 0.0]
Here rk.t is the list (history) of the calculated time. Index -1 specifies the “present time”.
Let’s plot the calculated trajectory. First, time series of positions are plotted.
>>> import matplotlib.pyplot as plt
>>> plt.ion()
>>> pos = np.array(rk.pos)
>>> plt.plot(rk.t, pos[:, 0], 'b')
>>> plt.plot(rk.t, pos[:, 1], 'g')
>>> plt.plot(rk.t, pos[:, 2], 'r')
Of course the particle moves to x-direction. Motion in y-direction is oscillation. One cycle (gyromotion) takes ~33s. You may look at the trajectory as well.
>>> plt.figure()
>>> plt.plot(pos[:, 0], pos[:, 1])
>>> plt.gca().set_aspect(1)
The velocity can also be plotted.
>>> import matplotlib.pyplot as plt
>>> plt.ion()
>>> vel = np.array(rk.vel)
>>> plt.plot(vel[:, 0], vel[:, 1])
>>> plt.gca().set_aspect(1)
ExB drift velocity in this case is 250 km/s, so the velocity plot should be the circle centered at 250 km/s.
Error analysis¶
The error is also recorded in rk.
>>> err = np.array(rk.err)
The error, with (N, 6) shaped array now where N is the number of steps, is just a difference between the 5th order and 4th order Runge-Kutta solver. Considering the accumulation, the total error at the end could be estimated as
Let’s see how much the error could be.
>>> print(np.sqrt((err ** 2).sum(0))) # sum(0) is to sum over the time steps.
[ 64.92358004 36.46324778 0. 6.98696963 6.74559555 0. ]
The 1-sigma is indeed 65 m for x-direction. Extremely small compared to 1.6e7 m (6 orders of magnitude smaller). However, remember that you have specified the tolerance to be 10 m (for a single step) and 0.1 m/s.
It may be interesting to plot the error (but invisible due to very high accuracy).
>>> plt.figure()
>>> err_x = err[:, 0]
>>> err_x2 = err_x ** 2
>>> err_x2_cumsum = err_x2.cumsum()
>>> err_x_array = np.sqrt(err_x2_cumsum)
>>> plt.plot(rk.t, pos[:, 0], 'b')
>>> plt.plot(rk.t, pos[:, 0] - err_x_array, 'b:')
>>> plt.plot(rk.t, pos[:, 0] + err_x_array, 'b:')
Relative tolerance¶
For comparison, we may also try relative tolerance, allowing 1% of the error. The same setting (ExB drift) is used.
>>> import numpy as np
>>> mp = 1.67e-27
>>> qe = 1.60e-19
>>> x0 = np.array([0, 0, 0])
>>> v0 = np.array([0, 0, 0])
>>> b = np.array([0, 0, 2e-9])
>>> e = np.array([0, 5e-4, 0])
>>> from irfpy.util import fields
>>> efield = fields.UniformVector3dField(e[0], e[1], e[2])
>>> bfield = fields.UniformVector3dField(b[0], b[1], b[2])
>>> from irfpy.util import rungekutta
>>> emforce = rungekutta.ElectroMagneticForce(mp, qe, efield, bfield)
>>> rk = rungekutta.RKF45_EqM(x0, v0, emforce, 1e-5, reltolx=0.01, reltolv=0.01, dtmax=np.inf)
Try progressing.
>>> rk.progress()
INFO:rkf45:T=0 : DT=1e-05 : MAXRELERR@4=0.039161037601691914
INFO:rkf45:T=0 : DT=1.5000000000000002e-05 : MAXRELERR@4=0.03916103760164462
INFO:rkf45:T=0 : DT=2.2500000000000005e-05 : MAXRELERR@4=0.0391610376015369
...
INFO:rkf45:T=0 : DT=2.876265888493262 : MAXRELERR@0=0.13561430468470423
INFO:rkf45:T=0 : DT=4.314398832739894 : MAXRELERR@0=0.44068725383770563
In this case, the time resolution of 4.3 s is used for the first step.
>>> while rk.t[-1] < 60.0:
... rk.progress()
INFO:rkf45:T=4.314398832739894 : DT=4.357542821067293 : MAXRELERR@4=0.18243816700337456
INFO:rkf45:T=4.314398832739894 : DT=6.536314231600939 : MAXRELERR@4=2.763637148973773
INFO:rkf45:T=4.314398832739894 : DT=3.2681571158004696 : MAXRELERR@4=0.05802640196626537
INFO:rkf45:T=4.314398832739894 : DT=4.902235673700704 : MAXRELERR@4=0.34614063850775195
INFO:rkf45:T=9.216634506440599 : DT=4.951258030437712 : MAXRELERR@3=0.16890430264235867
...
INFO:rkf45:T=54.688410054512175 : DT=4.294045522172377 : MAXRELERR@4=0.3286771762904559
INFO:rkf45:T=58.98245557668455 : DT=4.336985977394099 : MAXRELERR@3=0.1128833003912918
INFO:rkf45:T=58.98245557668455 : DT=6.5054789660911485 : MAXRELERR@1=0.5988161923772433
>>> print(rk)
t=65.4879345427757 x=[16418898.678595228 -16291.124124205708 0.0] v=[-3121.6525267938778 -8989.708819410907 0.0]
Force field of your own¶
You may want to define the force field of your own.
The force field can be a function (or precisely speaking callable) that eats 7 argument
(t [s], x [m], y, z, vx [m/s], vy, vz) returning the force vector (Fx, Fy, Fz) in N/kg.
Here let’s implement the following magnetic field.
Only Bz. (Bx=0, By=0)
Bz increases monotonically as Bz = t * 3e-8
The corresponding gyro period is \(T(t) = \frac{2\pi m}{qB} = 2.19 / t\).
Let’s implement the increasing magnetic field as follows.
import numpy as np
def increasing_magforce(t, x, y, z, vx, vy, vz):
v = np.array([vx, vy, vz]) # velocity vector
b = np.array([0, 0, 3e-8 * t]) # Magnetic field.
m = 1.67e-27 # Mass
q = 1.60e-19 # Charge
f = (q / m) * np.cross(v, b)
return f
You just use this function as the third argument of the irfpy.util.rungekutta.RKF45_EqM.
>>> x0 = np.array([0, 0, 0])
>>> v0 = np.array([100, 0, 0]) # 100 m/s in x direction
>>> from irfpy.util import rungekutta
>>> rk = rungekutta.RKF45_EqM(x0, v0, increasing_magforce, 0.1, tolx=10, tolv=0.1, dtmax=6)
Let’s step 100 steps.
>>> for i in range(100): rk.progress()
>>> print(rk)
t=7.407898341700332 x=[50.97878817476706 56.844392817730196 0.0] v=[-97.41650683152639 -30.68566545705031 0.0]
After 100 steps, 7.4 s was ellapsed. Let’s see the trajectory.
>>> import matplotlib.pyplot as plt
>>> pos = np.array(rk.pos)
>>> vel = np.array(rk.vel)
>>> plt.ion()
>>> plt.figure()
>>> plt.plot(rk.t, pos[:, 0])
>>> plt.plot(rk.t, pos[:, 1])
>>> plt.figure()
>>> plt.plot(pos[:, 0], pos[:, 1])
>>> plt.gca().set_aspect(1)
As an exercise, you can try different tolerance (e.g. tolx=0.1, tolv=0.001) and see how the results differe.