r''' Frame definition. Draft version.
This is a draft version of the frame definition.
Thus, I added "0" in the trail of the module name.
**Convert between the JDC0 angles and vectors**
The converion of the angles to JDC0 frame is implemented.
See :func:`angles2jdc` and :func:`jdc2angles`.
The definition follows the definition of myself.
Thus, the coordinate system may be chagned in the future.
The definition is consistent to the one defined in :mod:`fov0 <irfpy.jdc.fov0>`.
The elevation angle, :math:`\theta`, 0 is the zenith (z axis) and 180 for the nadir looking direction.
The azimutha angle, :math:`\phi`, 0 is x axis and 90 for y axis.
Angle :math:`\theta` and :math:`\phi` is the polar coordinate system
of the JDC0 frame.
**Convert between the JDC0 vectors and NSC vectors**
The frame conversion between JDC0 and NSC vectors is also implemented.
See :func:`jdc2nsc` and :func:`nsc2jdc`.
NSC means the *nadir looking spacecraft* coordinate system defined
at :mod:`irfpy.pep.pep_attitude`.
The definition is
- z-axis: Nadir direction, i.e. spacecraft to the body center.
- y-axis: y-z plane contains the velocity vector. The y-component
of the velocity should have positive.
- x-axis: Completes the right hand system.
The JDC0 frame has fixed relation with NSC, as the rigid mount of the sensor.
I define as follows for draft frame conversion.
.. math::
x_{NSC} = y_{JDC} \\
y_{NSC} = x_{JDC} \\
z_{NSC} = -z_{JDC}
This means that the positive velocity (:math:`y_{NSC}`) is in +x in JDC,
thus the velocity vector (ram directioN) is in between the CH-0 and CH-15,
i.e. :math:`\phi=0`.
**Convert from NSC vectors to physical frame**
The conversion is implemented in :mod:`irfpy.pep.pep_attitude` module.
'''
import numpy as np
[docs]def angles2jdc(theta, phi):
r''' Return the unit vector as ``np.array`` object corresponding to the angles.
:param theta: Theta in degrees. Float or (N,) shape np.array.
:type theta: float
:param phi: Phi in degrees. Float or (N,) shape np.array.
:type phi: float
:returns: The corresponding vector with shape of (3,) or (3, N)
:rtype: ``np.array``
>>> print(angles2jdc(0, 0))
[0. 0. 1.]
The following returns the corresponding vector for
:math:`(\theta, \phi) = (0, 0), (1, 1), (2, 2), ...`
at once.
>>> print(angles2jdc(np.arange(181), np.arange(181)).shape)
(3, 181)
'''
t = np.deg2rad(theta); p = np.deg2rad(phi)
x = np.sin(t) * np.cos(p)
y = np.sin(t) * np.sin(p)
z = np.cos(t)
return np.array([x, y, z])
[docs]def jdc2angles(vec):
''' Return the JDC angle for the corresponding JDC vector.
:param vec: np.array, with (3,) or (3, N) shape.
:type vec: np.array
:returns: Theta and Phi
:rtype: np.array with (2,) or (2, N) shape.
>>> print(jdc2angles([1, 0, 0]))
[90. 0.]
For multiple vectors, you can convert to angles.
Note that the shape should be (3, N), thus, the array should look like
``np.array([[x0, x1, x2, ...], [y0, y1, y2, ...], [z0, z1, z2, ...]])``
>>> print(jdc2angles([[1, 10, 0, 1], [0, 0, 1, 0], [0, 0, 0, 1]]))
[[90. 90. 90. 45.]
[ 0. 0. 90. 0.]]
'''
# vec is (3,) or (3, N)
veclen = np.sqrt((np.array(vec) ** 2).sum(0)) # Now, scalar or (N,)
vu = vec / veclen # Now, (3,) or (3, N)
t = np.arccos(vu[2])
p = np.arctan2(vu[1], vu[0])
return np.array([np.rad2deg(t), np.rad2deg(p)])
[docs]def nsc2jdc(vecnsc):
''' Convert from NSC to JDC.
vec should be (3,) or (3,N) array.
>>> nsc = [1, 2, -3]
>>> print(nsc2jdc(nsc))
[2 1 3]
>>> nscs = [[1, 1, 1, 1, 1], [2, 1, 0, -1, 2], [-3, -1, 0, 1, 3]]
>>> print(nsc2jdc(nscs))
[[ 2 1 0 -1 2]
[ 1 1 1 1 1]
[ 3 1 0 -1 -3]]
'''
mat = np.array([[0, 1, 0], [1, 0, 0], [0, 0, -1]])
vecjdc = mat.dot(vecnsc)
return vecjdc
[docs]def jdc2nsc(vecjdc):
''' Convert from NSC to JDC.
vec should be (3,) or (3,N) array.
>>> jdc = [1, 2, -3]
>>> print(jdc2nsc(jdc))
[2 1 3]
>>> jdcs = [ [ 2, 1, 0, -1, 2],
... [ 1, 1, 1, 1, 1],
... [ 3, 1, 0, -1, -3]]
>>> print(jdc2nsc(jdcs))
[[ 1 1 1 1 1]
[ 2 1 0 -1 2]
[-3 -1 0 1 3]]
'''
mat = np.array([[0, 1, 0], [1, 0, 0], [0, 0, -1]])
vecnsc = mat.dot(vecjdc)
return vecnsc
import unittest
import doctest
[docs]def doctests():
return unittest.TestSuite((
doctest.DocTestSuite(),
))
if __name__ == '__main__':
unittest.main(defaultTest='doctests')