###############################################################################
# NonInertialFrameForce: Class that implements the fictitious forces
# present when integrating orbits in a non-intertial
# frame
###############################################################################
import numpy
import numpy.linalg
import hashlib
from ..util import coords, conversion
from .DissipativeForce import DissipativeForce
[docs]class NonInertialFrameForce(DissipativeForce):
"""Class that implements the fictitious forces present when integrating
orbits in a non-intertial frame. Coordinates in the inertial frame
:math:`\mathbf{x}` and in the non-inertial frame :math:`\mathbf{r}` are
related through rotation and linear motion as
.. math::
\mathbf{x} = \mathbf{R}\,\left(\mathbf{r} + \mathbf{x}_0\\right)
where :math:`\mathbf{R}` is a rotation matrix and :math:`\mathbf{x}_0`
is the motion of the origin. The rotation matrix has angular frequencies
:math:`\\boldsymbol{\Omega}` with time derivative :math:`\dot{\\boldsymbol{\Omega}}`;
:math:`\\boldsymbol{\Omega}` can be any function of time. The motion of the
origin can also be any function of time.
This leads to the fictitious force
.. math::
\mathbf{F} = -\mathbf{a}_0 - \\boldsymbol{\Omega} \\times ( \\boldsymbol{\Omega} \\times \left[\mathbf{r} + \mathbf{x}_0\\right]) - \dot{\\boldsymbol{\Omega}} \\times \left[\mathbf{r}+\mathbf{x}_0\\right] -2\\boldsymbol{\Omega}\\times \left[\dot{\mathbf{r}}+\mathbf{v}_0\\right]
where :math:`\mathbf{a}_0`, :math:`\mathbf{v}_0`, and :math:`\mathbf{x}_0` are
the acceleration, velocity, and position of the origin of the non-inertial frame,
respectively, as a function of time. Note that if the non-inertial frame is not
rotating, it is not necessary to specify :math:`\mathbf{v}_0` and :math:`\mathbf{x}_0`.
In that case, the fictitious force is simply
.. math::
\mathbf{F} = -\mathbf{a}_0\quad (\\boldsymbol{\Omega} = 0)
If the non-inertial frame only rotates without any motion of the origin, the
fictitious force is the familiar combination of the centrifugal force
and the Coriolis force (plus an additional term if :math:`\dot{\\boldsymbol{\Omega}}`
is not constant)
.. math::
\mathbf{F} = - \\boldsymbol{\Omega} \\times ( \\boldsymbol{\Omega} \\times \mathbf{r}) - \dot{\\boldsymbol{\Omega}} \\times \mathbf{r} -2\\boldsymbol{\Omega}\\times \dot{\mathbf{r}}\quad (\mathbf{a}_0=\mathbf{v}_0=\mathbf{x}_0=0)
The functions of time are passed to the C code for fast orbit integration
by attempting to build fast ``numba`` versions of them. Significant
speed-ups can therefore be obtained by making sure that the provided
functions can be turned into ``nopython=True`` ``numba`` functions (try
running ``numba.njit`` on them and then evaluate them to check).
"""
[docs] def __init__(self,amp=1.,Omega=None,Omegadot=None,
x0=None,v0=None,a0=None,
ro=None,vo=None):
"""
NAME:
__init__
PURPOSE:
initialize a NonInertialFrameForce
INPUT:
amp= (1.) amplitude to be applied to the potential (default: 1)
Omega= (1.) Angular frequency of the rotation of the non-inertial frame in an inertial one; can either be a function of time or a number (when the frequency is assumed to be Omega + Omegadot x t) and in each case can be a list [Omega_x,Omega_y,Omega_z] or a single value Omega_z (when not a function, can be a Quantity; when a function, need to take input time in internal units and output the frequency in internal units; see galpy.util.conversion.time_in_Gyr and galpy.util.conversion.freq_in_XXX conversion functions)
Omegadot= (None) Time derivative of the angular frequency of the non-intertial frame's rotation. format should match Omega input ([list of] function[s] when Omega is one, number/list if Omega is a number/list; when a function, need to take input time in internal units and output the frequency derivative in internal units; see galpy.util.conversion.time_in_Gyr and galpy.util.conversion.freq_in_XXX conversion functions)
x0= (None) Position vector x_0 (cartesian) of the center of mass of the non-intertial frame (see definition in the class documentation); list of functions [x_0x,x_0y,x_0z]; only necessary when considering both rotation and center-of-mass acceleration of the inertial frame (functions need to take input time in internal units and output the position in internal units; see galpy.util.conversion.time_in_Gyr and divided physical positions by the `ro` parameter in kpc)
v0= (None) Velocity vector v_0 (cartesian) of the center of mass of the non-intertial frame (see definition in the class documentation); list of functions [v_0x,v_0y,v_0z]; only necessary when considering both rotation and center-of-mass acceleration of the inertial frame (functions need to take input time in internal units and output the velocity in internal units; see galpy.util.conversion.time_in_Gyr and divided physical positions by the `vo` parameter in km/s)
a0= (None) Acceleration vector a_0 (cartesian) of the center of mass of the non-intertial frame (see definition in the class documentation); constant or a list of functions [a_0x,a_0y, a_0z] (functions need to take input time in internal units and output the acceleration in internal units; see galpy.util.conversion.time_in_Gyr and galpy.util.conversion.force_in_XXX conversion functions [force is actually acceleration in galpy])
OUTPUT:
(none)
HISTORY:
2022-03-02 - Started - Bovy (UofT)
2022-03-26 - Generalized Omega to any function of time - Bovy (UofT)
"""
DissipativeForce.__init__(self,amp=amp,ro=ro,vo=vo)
self._rot_acc= not Omega is None
self._omegaz_only= len(numpy.atleast_1d(Omega)) == 1
self._const_freq= Omegadot is None
if (self._omegaz_only and callable(Omega)) \
or (not self._omegaz_only and callable(Omega[0])):
self._Omega_as_func= True
self._Omega= Omega
self._Omegadot= Omegadot
# Convenient access in Python
if not self._omegaz_only:
self._Omega_py= lambda t: numpy.array(\
[self._Omega[0](t),
self._Omega[1](t),
self._Omega[2](t)])
self._Omegadot_py= lambda t: numpy.array(\
[self._Omegadot[0](t),
self._Omegadot[1](t),
self._Omegadot[2](t)])
else:
self._Omega_py= self._Omega
self._Omegadot_py= self._Omegadot
else:
self._Omega_as_func= False
self._Omega= conversion.parse_frequency(Omega,ro=self._ro,vo=self._vo)
self._Omegadot= conversion.parse_frequency(Omegadot,ro=self._ro,vo=self._vo)
self._lin_acc= not (a0 is None)
if self._lin_acc:
if not callable(a0[0]):
self._a0= [lambda t, copy=a0[0]: copy,
lambda t, copy=a0[1]: copy,
lambda t, copy=a0[2]: copy]
else:
self._a0= a0
# Convenient access in Python
self._a0_py= lambda t: [self._a0[0](t),
self._a0[1](t),
self._a0[2](t)]
if self._lin_acc and self._rot_acc:
self._x0= x0
self._v0= v0
# Convenient access in Python
self._x0_py= lambda t: numpy.array(\
[self._x0[0](t),
self._x0[1](t),
self._x0[2](t)])
self._v0_py= lambda t: numpy.array(\
[self._v0[0](t),
self._v0[1](t),
self._v0[2](t)])
# Useful derived quantities
self._Omega2= numpy.linalg.norm(self._Omega)**2. \
if self._rot_acc and not self._Omega_as_func\
else 0.
if not self._omegaz_only and not self._Omega_as_func:
self._Omega_for_cross= numpy.array([[0.,-self._Omega[2],self._Omega[1]],
[self._Omega[2],0.,-self._Omega[0]],
[-self._Omega[1],self._Omega[0],0.]])
if not self._const_freq:
self._Omegadot_for_cross= \
numpy.array([[0.,-self._Omegadot[2],self._Omegadot[1]],
[self._Omegadot[2],0.,-self._Omegadot[0]],
[-self._Omegadot[1],self._Omegadot[0],0.]])
self._force_hash= None
self.hasC= True
return None
def _force(self,R,z,phi,t,v):
"""Internal function that computes the fictitious forces in rectangular
coordinates"""
new_hash= hashlib.md5(numpy.array([R,phi,z,v[0],v[1],v[2],t])).hexdigest()
if new_hash == self._force_hash:
return self._cached_force
x,y,z= coords.cyl_to_rect(R,phi,z)
vx,vy,vz= coords.cyl_to_rect_vec(v[0],v[1],v[2],phi)
force= numpy.zeros(3)
if self._rot_acc:
if self._const_freq:
tOmega= self._Omega
tOmega2= self._Omega2
elif self._Omega_as_func:
tOmega= self._Omega_py(t)
tOmega2= numpy.linalg.norm(tOmega)**2.
else:
tOmega= self._Omega+self._Omegadot*t
tOmega2= numpy.linalg.norm(tOmega)**2.
if self._omegaz_only:
force+= -2.*tOmega*numpy.array([-vy,vx,0.])\
+tOmega2*numpy.array([x,y,0.])
if self._lin_acc:
force+= -2.*tOmega*numpy.array([-self._v0[1](t),self._v0[0](t),0.])\
+tOmega2*numpy.array([self._x0[0](t),self._x0[1](t),0.])
if not self._const_freq:
if self._Omega_as_func:
force-= self._Omegadot_py(t)*numpy.array([-y,x,0.])
if self._lin_acc:
force-= self._Omegadot_py(t)\
*numpy.array([-self._x0[1](t),self._x0[0](t),0.])
else:
force-= self._Omegadot*numpy.array([-y,x,0.])
if self._lin_acc:
force-= self._Omegadot\
*numpy.array([-self._x0[1](t),self._x0[0](t),0.])
else:
if self._Omega_as_func:
self._Omega_for_cross= \
numpy.array([[0.,-self._Omega[2](t),self._Omega[1](t)],
[self._Omega[2](t),0.,-self._Omega[0](t)],
[-self._Omega[1](t),self._Omega[0](t),0.]])
if not self._const_freq:
self._Omegadot_for_cross= \
numpy.array([[0.,-self._Omegadot[2](t),self._Omegadot[1](t)],
[self._Omegadot[2](t),0.,-self._Omegadot[0](t)],
[-self._Omegadot[1](t),self._Omegadot[0](t),0.]])
force+= -2.*numpy.dot(self._Omega_for_cross,numpy.array([vx,vy,vz]))\
+tOmega2*numpy.array([x,y,z])\
-tOmega*numpy.dot(tOmega,numpy.array([x,y,z]))
if self._lin_acc:
force+= -2.*numpy.dot(self._Omega_for_cross,self._v0_py(t))\
+tOmega2*self._x0_py(t)\
-tOmega*numpy.dot(tOmega,self._x0_py(t))
if not self._const_freq:
if not self._Omega_as_func: # Already included above when Omega=func
force-= 2.*t*numpy.dot(self._Omegadot_for_cross,numpy.array([vx,vy,vz]))
force-= numpy.dot(self._Omegadot_for_cross,numpy.array([x,y,z]))
if self._lin_acc:
if not self._Omega_as_func:
force-= 2.*t*numpy.dot(self._Omegadot_for_cross,self._v0_py(t))
force-= numpy.dot(self._Omegadot_for_cross,self._x0_py(t))
if self._lin_acc:
force-= self._a0_py(t)
self._force_hash= new_hash
self._cached_force= force
return force
def _Rforce(self,R,z,phi=0.,t=0.,v=None):
"""
NAME:
_Rforce
PURPOSE:
evaluate the radial force for this Force
INPUT:
R - Galactocentric cylindrical radius
z - vertical height
phi - azimuth
t - time
v= current velocity in cylindrical coordinates
OUTPUT:
the radial force
HISTORY:
2022-03-02 - Written - Bovy (UofT)
"""
force= self._force(R,z,phi,t,v)
return numpy.cos(phi)*force[0]+numpy.sin(phi)*force[1]
def _phiforce(self,R,z,phi=0.,t=0.,v=None):
"""
NAME:
_phiforce
PURPOSE:
evaluate the azimuthal force for this Force
INPUT:
R - Galactocentric cylindrical radius
z - vertical height
phi - azimuth
t - time
v= current velocity in cylindrical coordinates
OUTPUT:
the azimuthal force
HISTORY:
2022-03-02 - Written - Bovy (UofT)
"""
force= self._force(R,z,phi,t,v)
return R*(-numpy.sin(phi)*force[0]+numpy.cos(phi)*force[1])
def _zforce(self,R,z,phi=0.,t=0.,v=None):
"""
NAME:
_zforce
PURPOSE:
evaluate the vertical force for this Force
INPUT:
R - Galactocentric cylindrical radius
z - vertical height
phi - azimuth
t - time
v= current velocity in cylindrical coordinates
OUTPUT:
the vertical force
HISTORY:
2022-03-02 - Written - Bovy (UofT)
"""
return self._force(R,z,phi,t,v)[2]