''' Spacecraft position and attitude
The module :mod:`irfpy.cy1orb.pvat2` provides functionality of
- Spacecraft position
- Position in LSE frame
- Position in ME frame
- Spacecraft attitude
- Convert between LSE frame and spacecraft (SC) frame
For frame conversions, please refer also to :ref:`frame_conversion`.
*Setup*
See :ref:`get_start_positon` for the setup.
.. warning::
The returned values are in :mod:`irfpy.util.vector3d`, which should be
modified in some ways.
*Getting position*
To get the position of spacecraft, one can use either of the function.
- :func:`irfpy.cy1orb.pvat2.getlsepos`: To get LSE position.
- :func:`irfpy.cy1orb.pvat2.getmepos`: To get ME position.
>>> import datetime
>>> from irfpy.cy1orb import pvat2 as pvat
>>> t = datetime.datetime(2009, 5, 10, 3, 10)
>>> pos_lse = pvat.getlsepos(t, asarray=True)
>>> print(pos_lse)
[1396.34 1005.16 633.71]
>>> pos_me_lonlatalt = pvat.getmepos(t)
>>> print(pos_me_lonlatalt) # The returned value is longitude, latitude and altitude
(28.6, 21.2, 96.1)
>>> pos_me_xyz = pvat.getmepos(t, asarray=True)
>>> print(pos_me_xyz)
[1501.32901199 818.5511499 663.25562399]
*Conversion between LSE and SC*
The spacecraft frame and LSE frame conversion can be doen using
- :func:`lse2sc`, from LSE to SC
- :func:`sc2lse`, from SC to LSE
The example below will give the representation of the spacecraft x-axis in the LSE frame.
>>> vec = sc2lse(datetime.datetime(2009,1,25,14,0,0), [1, 0, 0])
>>> print(vec)
[-0.8204 0.5718 -0.0042]
The next example will give the representations of LSE x-axis in the SC frame.
Since LSE a-axis is toward the Sun, the returned vector is the Sun direction in the SC frame.
>>> vec = lse2sc(datetime.datetime(2009,1,25,14,0,0), [1, 0, 0])
>>> print(vec)
[-0.8204 0.0147 0.57163338]
To convert the LSE frame to ME (lunar fixed) frame, see the module :mod:`irfpy.cy1orb.lseme`.
'''
import datetime
import logging
import urllib.request
import urllib.parse
import urllib.error
import numpy as np
from irfpy.util import utc, ringcache, vector3d, julday
from irfpy.cy1orb.Cy1OrbitNr import Cy1OrbitNr
from irfpy.cy1orb.Cy1Orbit import Cy1Orbit, Cy1Attitude
import irfpy.util.irfpyrc
import irfpy.util.utc
import irfpy.util.julday
from irfpy.util.timeinterval import timeinterval
import irfpy.cy1orb.Cy1Orbit
logger = logging.getLogger(__name__)
#__cy1orbitnr = Cy1OrbitNr()
#__cy1orbitnr.setFromDefaultUri()
class __cy1orbitnr:
instance = None
@classmethod
def cache(cls):
if cls.instance is None:
cls.instance = Cy1OrbitNr()
cls.instance.setFromDefaultUri()
return cls.instance
# This is a ring buffer for an instance of Cy1Orbit.
# Even though Cy1Orbit has a cache support, but it is not ring cache.
# Therefore, in pvat2, an instance of Cy1Orbit is "1 instance 1 orbit".
# The Cy1Orbit instance will be cached to the ringcache with size 30.
__cy1orbitcache = ringcache.SimpleRingCache(30)
def __getCy1Orbit(orb):
''' Get the Cy1Orbit instance for the corresponding orbit number.
If the specified orbit number is in the cache, just return.
If not, make new instance of Cy1Orbit corresponding to the orbit,
register it to the cache, and return the Cy1Orbit instance.
'''
if __cy1orbitcache.hasKey(orb):
return __cy1orbitcache.get(orb)
else:
c = Cy1Orbit()
c.loadOrbit(orb)
__cy1orbitcache.add(orb, c)
return c
__cy1attcache = ringcache.SimpleRingCache(30)
def __getCy1Att(orb):
if __cy1attcache.hasKey(orb):
return __cy1attcache.get(orb)
else:
c = Cy1Attitude()
c.loadOrbit(orb)
__cy1attcache.add(orb, c)
return c
# Below cache is for more specific data.
#
#
__cy1lseorbitcache = ringcache.SimpleRingCache(30)
def __getCy1LseOrbitCache(orb):
if __cy1lseorbitcache.hasKey(orb):
return __cy1lseorbitcache.get(orb)
else:
c = __getCy1Orbit(orb)
lse = c.getPosOfOrb2(orb)
__cy1lseorbitcache.add(orb, lse)
return lse
__cy1meorbitcache = ringcache.SimpleRingCache(30)
def __getCy1MeOrbitCache(orb):
if __cy1meorbitcache.hasKey(orb):
return __cy1meorbitcache.get(orb)
else:
c = __getCy1Orbit(orb)
me = c.getSelenographicOfOrb2(orb)
__cy1meorbitcache.add(orb, me)
return me
__cy1sc2lsecache = ringcache.SimpleRingCache(30)
def __getSc2LseCache(orb):
if __cy1sc2lsecache.hasKey(orb):
return __cy1sc2lsecache.get(orb)
else:
c = __getCy1Att(orb)
cnv = c.getSc2LseMatrixOfOrb(orb)
__cy1sc2lsecache.add(orb, cnv)
return cnv
__cy1sc2lsematxcache = ringcache.SimpleRingCache(300) # This is a cache of Jd->matrix
def __getSc2LseMatxCache(t):
t0 = utc.convert(t, datetime.datetime)
if __cy1sc2lsematxcache.hasKey(t0):
return __cy1sc2lsematxcache.get(t0)
else:
o0 = __cy1orbitnr.cache().getOrbitNr(t0)
orbdata = __getSc2LseCache(o0)
matx = orbdata.getNeighbor(utc.convert(t0, julday.Julday)).getData()
__cy1sc2lsematxcache.add(t0, matx)
return matx
__cy1lse2scmatxcache = ringcache.SimpleRingCache(300)
def __getLse2ScMatxCache(t):
t0 = utc.convert(t, datetime.datetime)
if __cy1lse2scmatxcache.hasKey(t0):
return __cy1lse2scmatxcache.get(t0)
else:
o0 = __cy1orbitnr.cache().getOrbitNr(t0)
orbdata = __getSc2LseCache(o0)
matx = orbdata.getNeighbor(utc.convert(t0, julday.Julday)).getData().T
__cy1lse2scmatxcache.add(t0, matx)
return matx
__clear_cache_warned = False
[docs]def clear_cache():
''' Clear the cache.
Usually user do not need to use this function.
'''
global __clear_cache_warned
if not __clear_cache_warned:
logger.debug(('Because pvat2 module uses ring cache,'
'calling clear_cache() does not make '
'a significant memory release.'))
__clear_cache_warned = True
__cy1orbitcache.clearCache()
__cy1attcache.clearCache()
__cy1lseorbitcache.clearCache()
__cy1meorbitcache.clearCache()
__cy1sc2lsecache.clearCache()
__cy1sc2lsematxcache.clearCache()
[docs]def getlsepos(t, asarray=False):
''' Get LSE orbit at the time of t.
The LSE postiion of Chandrayaan-1.
The unit is km.
The returned position of the spacecraft is at the sampling time
(1 sec resolution) closest to the specified time.
:param t: Time of the data to inquire
:keyword asarray: If *True*, numpy array is returned. Otherwise, :class:`Vector3d` is used.
:returns: The position of spacecraft.
:rtype: :class:`Vector3d` (default) or :class:`np.narray` (*as_array=True*)
'''
t0 = utc.convert(t, datetime.datetime)
o0 = __cy1orbitnr.cache().getOrbitNr(t0)
logger.debug('T=%s @ %d' % (t0, o0))
orbdata = __getCy1LseOrbitCache(o0)
v = orbdata.getNeighbor(utc.convert(t0, julday.Julday)).getData()
if asarray:
return np.array([v.x, v.y, v.z])
else:
return vector3d.Vector3d(v.x, v.y, v.z)
[docs]def getsza(t):
''' Get the solar zenith angle at the specfied time.
'''
lsepos = getlsepos(t)
sundir = vector3d.Vector3d(1, 0, 0)
angle = lsepos.angle(sundir)
return angle * 180. / np.pi
[docs]def getmepos(t, as_vector=None, asarray=False, asvector=False):
''' Get s/c position in ME frame as an array of [lon,lat,height].
lon and lat is in degrees, height is in km above the surface.
>>> me=getmepos(datetime.datetime(2009,4,18,1,30,0))
:param: Time to inquire
:keyword as_vector: Deprecated. Same effect as ``asarray``
:keyword asarray: Return as numpy array, in km.
:keyword asvector: Return as :class:`Vector3d` class object.
:returns: The position of the spacecraft. Default is (longitude, latitude, and height).
Angles in degrees, and height in km.
If *asarray* is True, the returned is ``[x, y, z]`` in the ME frame in km.
If *asvector* is True, the returns is ``Vector3d(x, y, z)``.
'''
if as_vector is not None: # as_vector is explicitly stated
logger.warning('getmepos() takes keyword asvector instead of as_vector that is deprecated')
asarray = as_vector # For consistency
t0 = utc.convert(t, datetime.datetime)
o0 = __cy1orbitnr.cache().getOrbitNr(t0)
orbdata = __getCy1MeOrbitCache(o0)
v = orbdata.getNeighbor(utc.convert(t0, julday.Julday)).getData()
if asarray or asvector:
lon, lat, hei = v
r = hei + 1738.
lon = lon * np.pi / 180.
lat = lat * np.pi / 180.
x = r * np.cos(lat) * np.cos(lon)
y = r * np.cos(lat) * np.sin(lon)
z = r * np.sin(lat)
if asarray:
return np.array([x, y, z])
elif asvector:
return vector3d.Vector3d(x, y, z)
else:
return tuple(v)
[docs]def getlsevec(t, vec_in_sc):
''' Convert the SC -> LSE.
Get LSE-frame vector correspoinding to the SC-frame vector at the
time of t. The returned variable is vector3d.Vector3d instance.
>>> vec = getlsevec(datetime.datetime(2009,1,25,14,0,0),
... vector3d.Vector3d(1, 0, 0))
>>> print(vec.x, vec.y, vec.z)
-0.8204 0.5718 -0.0042
'''
t0 = utc.convert(t, datetime.datetime)
matx = __getSc2LseMatxCache(t0)
logger.debug('sc = %s' % vec_in_sc)
logger.debug('matx = %s' % matx)
if isinstance(vec_in_sc, vector3d.Vector3d):
vec_sc = [vec_in_sc.x, vec_in_sc.y, vec_in_sc.z]
else:
vec_sc = vec_in_sc[:]
vec_lse = np.dot(matx, vec_sc)
logger.debug('lse = %s' % vec_lse)
return vector3d.Vector3d(vec_lse[0], vec_lse[1], vec_lse[2])
[docs]def sc2lse(t, arr_sc):
""" Convert SC -> LSE
Get LSE-frame arrray correspoinding to the given SC-frame array at the
time of t.
:param t: Time
:param arr_sc: (..., 3) shaped array in SC frame
:returns: (..., 3) shaped array in LSE frame
>>> vec = sc2lse(datetime.datetime(2009,1,25,14,0,0),
... [1, 0, 0])
>>> print(vec)
[-0.8204 0.5718 -0.0042]
>>> arr_sc = [[1, 0, 0], [0, 1, 0]]
>>> arr_lse = sc2lse(datetime.datetime(2009,1,25,14,0,0), arr_sc)
>>> print(arr_lse)
[[-0.8204 0.5718 -0.0042]
[ 0.0147 0.0284 0.9995]]
"""
t0 = utc.convert(t, datetime.datetime)
matx = __getSc2LseMatxCache(t0)
arr_rev = np.transpose(arr_sc) # Now (3, ...) array (reversed)
arr_lse = np.dot(matx, arr_rev) # Now (3, ...) array
arr_lse = np.transpose(arr_lse)
return arr_lse
[docs]def get_sc2lse_matrix(t):
t0 = utc.convert(t, datetime.datetime)
matx = __getSc2LseMatxCache(t0)
return matx
[docs]def getscvec(t, vec_in_lse):
''' Convert LSE -> SC.
Inverse function of :func:`getlsevec`.
Conversion from LSE vector to SC vector at the time of t.
The returned is vector3d.Vector3d instance.
>>> vec = getscvec(datetime.datetime(2009, 1, 25, 14, 0, 0),
... vector3d.Vector3d(1, 0, 0))
>>> print('{:.4f} {:.4f} {:.4f}'.format(vec.x, vec.y, vec.z))
-0.8204 0.0147 0.5716
'''
t0 = utc.convert(t, datetime.datetime)
matx = __getLse2ScMatxCache(t0)
if isinstance(vec_in_lse, vector3d.Vector3d):
vec_lse = [vec_in_lse.x, vec_in_lse.y, vec_in_lse.z]
else:
vec_lse = vec_in_lse[:]
vec_sc = np.dot(matx, vec_lse)
return vector3d.Vector3d(vec_sc[0], vec_sc[1], vec_sc[2])
[docs]def lse2sc(t, arr_lse):
''' Convert LSE -> SC.
Inverse function of :func:`sc2lse`.
Conversion from LSE vector to SC vector at the time of t.
:param t: Time
:param arr_sc: (..., 3) shaped array in SC frame
:returns: (..., 3) shaped array in LSE frame
>>> vec_sc = lse2sc(datetime.datetime(2009, 1, 25, 14, 0, 0), [1, 0, 0])
>>> print('{v[0]:.4f} {v[1]:.4f} {v[2]:.4f}'.format(v=vec_sc))
-0.8204 0.0147 0.5716
'''
t0 = utc.convert(t, datetime.datetime)
matx = __getLse2ScMatxCache(t0)
arr_rev = np.transpose(arr_lse) # Now (3, ...) array (reversed order)
arr_sc = np.dot(matx, arr_rev) # Now (3, ...) array (reversed order)
arr_sc = np.transpose(arr_sc)
return arr_sc
[docs]def get_lse2sc_matrix(t):
t0 = utc.convert(t, datetime.datetime)
matx = __getLse2ScMatxCache(t0)
return matx
[docs]def getdatarange():
''' Return the data range.
:returns: The range of the data that can be accessed.
:rtype: ``util:irfpy.util.timeinterval.timeinterval`` instance.
>>> trange = getdatarange()
>>> print(datetime.datetime(2009, 3, 1, 0) in trange)
True
>>> print(datetime.datetime(2008, 12, 8, 0) in trange)
False
>>> print(datetime.datetime(2009, 9, 1, 0) in trange)
False
'''
rc = irfpy.util.irfpyrc.Rc()
path = rc.get('cy1orb', 'pvatbase')
logger.debug(path)
i = 0
while i < 10000: # At a certain point, you have to be stopped.
try:
uri = path + '/pvat_%04d.gz' % i
urllib.request.urlretrieve(uri)
except:
i += 1
continue
logger.debug('First orbit is %d' % i)
oi = irfpy.cy1orb.Cy1Orbit.Cy1Orbit()
t0 = oi.getPosOfOrb(i)[0, 1] # The second data as t0. This is because the first data
# tries to refer the data file before.
logger.debug(t0)
t0 = irfpy.util.utc.convert(irfpy.util.julday.Julday(t0), datetime.datetime)
logger.debug(t0)
break
while i < 10000:
try:
uri = path + '/pvat_%04d.gz' % i
urllib.request.urlretrieve(uri)
except:
logger.debug('Last orbit is %d' % (i - 1))
oi = irfpy.cy1orb.Cy1Orbit.Cy1Orbit()
t1 = oi.getPosOfOrb(i - 1)[0, -1]
logger.debug(t1)
t1 = irfpy.util.utc.convert(irfpy.util.julday.Julday(t1), datetime.datetime)
logger.debug(t1)
break
else:
i += 1
return timeinterval(t0, t1)
class __Doctest4Develop:
""" This is a purpose for doctesting. Not used for any purpose.
>>> t00 = datetime.datetime(2009, 1, 25, 14, 0, 0)
>>> __cy1sc2lsematxcache.clearCache()
>>> print(__cy1sc2lsematxcache.hasKey(t00))
False
>>> vec = getlsevec(t00, vector3d.Vector3d(1, 0, 0))
>>> print(vec)
Vector3d( -0.8204, 0.5718, -0.0042 )
>>> print(__cy1sc2lsematxcache.hasKey(t00))
True
>>> # This call should find from matx cache.
>>> vec = getlsevec(t00, vector3d.Vector3d(0, 1, 0))
>>> print(vec)
Vector3d( 0.0147, 0.0284, 0.9995 )
>>> vec = getlsevec(t00, vector3d.Vector3d(0, 0, 1))
>>> print(vec)
Vector3d( 0.571633, 0.819928, -0.0317048 )
>>> t11 = datetime.datetime(2009, 1, 25, 14, 10, 0)
>>> print(__cy1sc2lsematxcache.hasKey(t11))
False
>>> vec = getlsevec(t11, vector3d.Vector3d(0, 0, 1))
>>> print(__cy1sc2lsematxcache.hasKey(t11))
True
>>> print(__cy1sc2lsematxcache.hasKey(t00))
True
>>> clear_cache()
>>> print(__cy1sc2lsematxcache.hasKey(t00))
False
>>> vec = getlsevec(t00, vector3d.Vector3d(0, 0, 1))
>>> print(__cy1sc2lsematxcache.hasKey(t00))
True
>>> print(vec)
Vector3d( 0.571633, 0.819928, -0.0317048 )
"""