irfpy.util.leapflog

A module for particle tracing by leap flog.

Code author: Yoshifumi Futaana

Note

This module replaces irfpy.util.ptracing, though the interface is not identical. There is no reason using irfpy.util.ptracing from now on.

Leap flog method is implemented as a class LeapFlog. Indeed, this class can also be used as a simple forward method.

Use :class:`LeapFlog` as leap flog

  • pos(t) and vel(t+0.5) will give pos(t+1)

  • vel(t+0.5) and force(t+1) will give vel(t+1.5)

Use :class:`LeapFlog` as simple forward

  • pos(t) and vel(t) will give pos(t+1)

  • vel(t) and force(t) will give vel(t+1)

See details in LeapFlog.

Note

Development log: This is originally developed under irfpy.collision module. This module (irfpy.util) is the maintained version currently.

LeapFlog()

A leap flog (or simple forward) solver.

irfpy.util.keplernumeric

Kepler's trajectory using numeric solvers
class irfpy.util.leapflog.LeapFlog[source]

Bases: object

A leap flog (or simple forward) solver.

As a tutorial, we try to simulate the Moon motion around the Earth. Use MKSA for simplicity.

>>> G = 6.67e-11   # Grav. const.
>>> me = 5.97e+24   # Mass of earth
>>> mm = 7.35e+22   # Mass of moon

The equation of motion is

\[m_m \frac{d\mathbf{v}}{dt} = -G \frac{m_m m_e}{|\mathbf{r}|^3}\mathbf{r}\]

First task to be done is to implement the force term. The force term is the force acting to a unit mass of the trace object. In the case of the Moon, it is \(-G m_e \mathbf{r}/ r^3\).

>>> force = lambda r: -G * me * r / ((r * r).sum() ** 1.5)

Then, one need to specify the initial condition. Here, the velocity of the Moon 1.022 km/s at the distance of 3.84e5 km is assumed.

>>> pos0 = np.array([3.84e5, 0, 0]) * 1e3   # in m
>>> vel0 = np.array([0, 1.022, 0]) * 1e3    # in m/s

In addition, parameters related to the time steps are defined. In this tutorial, I would take just \(\Delta t=\) 1 hour (=3600 s).

>>> t = 0.
>>> dt = 3600.

Note that the time step should be chosen by the user’s responsibility affecting accuracy and calculation time and resource.

The initial velocity is, by intrinsic definition, given at the time \(t=0\). For the leap flog solver, we want the velocity at \(t=\Delta t/2\).

The simplest way is to assume \(v(t=0) = v(t=\Delta t/2)\). This looks too simplified, but not a bad idea for many cases in practice.

>>> vel0h = vel0    # velocity at the half step beyond.

More precise idea is to proceed the velocity by \(\Delta t/2\) following the equation of motion but in a forward method.

>>> vel0h = LeapFlog.progress_velocity(vel0, force(pos0), dt / 2.)

Now, all the needed information is defined. Let’s proceed the time.

>>> while t <= 86400 * 20:   # For 20 days.
...     # pos0 is at t, pos1 is at t+dt
...     pos1 = LeapFlog.progress_position(pos0, vel0h, dt)
...     # vel0h is at t+0.5dt, vel1h is at t+1.5dt.
...     vel1h = LeapFlog.progress_velocity(vel0h, force(pos1), dt)
...     # Update the present time, t.
...     pos0 = pos1
...     vel0h = vel1h
...     t = t + dt
>>> print(t)
1731600.0
>>> print(pos0)   
[-7.10934633e+07  -3.80723849e+08   0.00000000e+00]
>>> print(vel0h)  
[ 998.26945656 -174.1793034     0.        ]
classmethod progress_position(x, v, dt)[source]

Progress position according to the velocity.

The classmethod will progress 1 step of position.

The equation is

\[\frac{d\mathbf{x}}{dt} = \mathbf{v}\]

In the discrete form,

\[\mathbf{x}_{i+1} = \mathbf{x}_i + \mathbf{v}_{i+1/2} \times \Delta t\]
Parameters:

  • x (A np.array with shape (3,)) – Position vector at the time \(t\).

  • v (A np.array with shape (3,)) – Velocity vector at the time \(t + \Delta t / 2\). If you specify velocity at \(t\), this method can be used for a simple forward progress.

  • dt – \(\Delta t\)

Returns:

Updated position vector.

Return type:

A np.array with shape (3,)

A simple use is as follows.

>>> pos = np.array([1.5, 2.8, -2.7])
>>> vel = np.array([0.2, -1.4, 3.2])
>>> dt = 0.5
>>> newpos = LeapFlog.progress_position(pos, vel, dt)
>>> print(newpos)
[ 1.6  2.1 -1.1]
classmethod progress_velocity(v, force, dt)[source]

Progress velocity according to the force.

This progress 1 step (\(\Delta t\)) of velocity.

The equation to be solved is

\[\frac{d\mathbf{v}}{dt} = \frac{\mathbf{F}}{m} = \mathbf{f}\]

or

\[\mathbf{v}_{i+1} = \mathbf{v}_i + \mathbf{f}_{i+1/2} \times \Delta t\]

Note that the mass should be included in the force vector. This indicates that the force acting to the object is per unit mass.

Parameters:

  • v (np.array with shape (3,)) – Velocity vector

  • force – Force (per mass unit) vector. Mass of the particle should be considered in the force term. This means that the unit of this quantities is the acceleration, or [N/kg]. The force vector should be given at the time of \(t + \Delta t / 2\). This is ok if the force only depend on the position. However, if the force depends on the velocity, this calculation is difficult since we do not know the velocity at \(t + \Delta t / 2\). In such a case, one has to consider sacrifice the accuracy, or use other methods like Runge-Kutta.

>>> v = np.array([1.3, -1.9, 2.5])
>>> f = np.array([-0.8, 1.6, -2.4])
>>> dt = 0.125
>>> vnew = LeapFlog.progress_velocity(v, f, dt)
>>> print(vnew)
[ 1.2 -1.7  2.2]
irfpy.util.leapflog.doctests()[source]