''' A function to get a terminator in ME frame.
The terminator, the day-night boundary of the planet, is calculated for Moon case.
In theory, the Moon has a different polar radius than the equator radius,
but for the first trial, the Moon is assumed to be a sphere.
The terminator is simply described in LSE frame: x=0, y=Rm cos(t), z=Rm sin(t).
Here Rm is the radius of the Moon, and t is the angle ranging from 0 to 2pi.
x,y, and z are functions of the parameter of t.
This coordinates are converted to in ME frame: X,Y,Z by lse2me() method.
'''
import datetime
from irfpy.util.utc import convert
from irfpy.util.julday import Julday
from irfpy.cy1orb import lseme
from irfpy.cy1orb.lseme import lse2me, me2lse
import logging
import math
import numpy as np
[docs]def terminator(jd, nSample=361):
''' Get list of terminator coordinates (x,y,z) in ME frame.
:param jd: Julian day or the time in ``datetime.datetime`` instance.
:keyword nSample: Number of samples.
:returns: The terminator position vector in the ME frame.
:rtype: List of ``numpy.array``.
>>> t0 = datetime.datetime(2009, 3, 10, 0, 0, 0)
>>> tpos = terminator(t0)
>>> print(len(tpos))
361
>>> print(tpos[0].shape)
(3,)
>>> print(tpos[150]) # May change if DB changes.
[ 450.13942892 -1456.29471099 834.99925307]
'''
Rm=1738.
term=[]
for i in range(nSample):
t = 2. / (nSample-1) * math.pi * i
### First, calcuate the terminator coordinates in LSE frame.
x=0.
y=Rm * math.cos(t)
z=Rm * math.sin(t)
v=[x,y,z]
jd2 = convert(jd, Julday)
### Then, convert it to ME frame.
vme = lse2me(jd2, v)
term.append(vme)
return term
[docs]def terminator_asdegrees(jd, nSample=361):
''' Return the array of the coordinate system, (longitude, latitude, height).
:param jd: Julian day or the time in ``datetime.datetime`` instance.
:keyword nSample: Number of samples.
:returns: ``np.array``. Shape with [nSample, 3]. Longitude ranges -180 to 180.
Latitude ranges -90 to 90. Height should be almost zero.
>>> t0 = datetime.datetime(2009, 3, 10, 0, 0, 0)
>>> tpos = terminator_asdegrees(t0)
>>> print(tpos.shape)
(361, 3)
>>> print(tpos[150])
[-7.28236664e+01 2.87139399e+01 -1.31371162e-04]
'''
Rm = 1738.
pos = np.array(terminator(jd, nSample=nSample)) # (nSample, 3)
heis = np.sqrt((pos * pos).sum(1))
lons = np.rad2deg(np.arctan2(pos[:, 1], pos[:, 0]))
lats = np.rad2deg(np.arcsin(pos[:, 2] / heis))
heis -= Rm
return np.array([lons, lats, heis]).T