''' A module for particle tracing by leap flog.
.. codeauthor:: Yoshifumi Futaana
.. note::
This module replaces :mod:`irfpy.util.ptracing`,
though the interface is not identical.
There is no reason using :mod:`irfpy.util.ptracing`
from now on.
Leap flog method is implemented as a class :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 :class:`LeapFlog`.
.. note::
Development log: This is originally developed under
``irfpy.collision`` module.
This module (irfpy.util) is the maintained version currently.
.. autosummary::
LeapFlog
irfpy.util.keplernumeric
'''
import numpy as np
[docs]class LeapFlog:
r''' 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
.. math::
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 :math:`-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 :math:`\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 :math:`t=0`.
For the leap flog solver, we want the velocity at :math:`t=\Delta t/2`.
The simplest way is to assume :math:`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 :math:`\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) # doctest: +NORMALIZE_WHITESPACE
[-7.10934633e+07 -3.80723849e+08 0.00000000e+00]
>>> print(vel0h) # doctest: +NORMALIZE_WHITESPACE
[ 998.26945656 -174.1793034 0. ]
'''
[docs] @classmethod
def progress_position(self, x, v, dt):
r''' Progress position according to the velocity.
The classmethod will progress 1 step of position.
The equation is
.. math::
\frac{d\mathbf{x}}{dt} = \mathbf{v}
In the discrete form,
.. math::
\mathbf{x}_{i+1} = \mathbf{x}_i + \mathbf{v}_{i+1/2} \times \Delta t
:param x: Position vector at the time :math:`t`.
:type x: A np.array with shape (3,)
:param v: Velocity vector at the time :math:`t + \Delta t / 2`.
If you specify velocity at :math:`t`, this method
can be used for a simple forward progress.
:type v: A np.array with shape (3,)
:param dt: :math:`\Delta t`
:returns: Updated position vector.
:rtype: 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]
'''
### xnew = xold + vel * dt.
xnew = x + v * dt
return xnew
[docs] @classmethod
def progress_velocity(self, v, force, dt):
r''' Progress velocity according to the force.
This progress 1 step (:math:`\Delta t`) of velocity.
The equation to be solved is
.. math::
\frac{d\mathbf{v}}{dt} = \frac{\mathbf{F}}{m} = \mathbf{f}
or
.. math::
\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.
:param v: Velocity vector
:type v: np.array with shape (3,)
:param 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
:math:`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 :math:`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]
'''
vnew = v + force * dt
return vnew
import unittest
import doctest
[docs]def doctests():
return unittest.TestSuite((
doctest.DocTestSuite(),
))
if __name__ == '__main__':
unittest.main(defaultTest='doctests')