''' Module for getting moon position w.r.t Earth. Database is from Par-Ola's work.
Default database is from .irfpyrc. Specify [cy1orb] moongsebaseuri entry.
..
[cy1orb]
moongsebaseuri=file:///Volumes/scidata/data/moon/support/peje/
'''
import urllib.request
import urllib.parse
import urllib.error
import numpy
import logging
import datetime
from irfpy.util.julday import Julday
from irfpy.util.utc import convert
from irfpy.util.irfpyrc import Rc
[docs]class MoonGse:
''' Class of Moon position in GSE coordinate system.
'''
def __init__(self, baseurl=None):
''' Constructor
@param[in] baseurl Set default URL base which contains the data. Default is in .irfpyrc [cy1orb] moongsebaseuri.
'''
if baseurl is None:
rc = Rc()
baseurl = rc.get('cy1orb', 'moongsebaseuri')
self._baseurl=baseurl
[docs] def setFromDefaultUri(self):
''' Data will be loaded from the default URI.
'''
datafiles=['2008/2008-12.dat', '2009/2009-01.dat', '2009/2009-02.dat', '2009/2009-03.dat', '2009/2009-04.dat', '2009/2009-05.dat', '2009/2009-06.dat', '2009/2009-07.dat', '2009/2009-08.dat']
self._darr=None
for di in datafiles:
url=self._baseurl+'/'+di;
logging.debug('Load data from: %s'%url)
f=urllib.request.urlretrieve(url)
logging.debug('Load data from local file: %s'%f[0])
d=numpy.loadtxt(f[0]) # Nx4
logging.debug('Data shape %s'%str(d.shape))
if self._darr is None:
self._darr=d
else:
self._darr=numpy.concatenate([self._darr, d])
logging.debug('Merged data shpe %s'%str(self._darr.shape))
self._darr=self._darr.transpose()
logging.debug('Merged data shpe %s'%str(self._darr.shape))
self._datetime=[]
for i in range(len(self._darr[0,:])):
jd=self._darr[0,i]+0.5*Julday.SECOND
jd=Julday(jd)
dt=jd.getDatetime()
dt=datetime.datetime(dt.year, dt.month, dt.day, dt.hour, dt.minute, dt.second, 0) # Round millisecs
self._datetime.append(dt)
### self._darr is 4xN elements float: jd,x,y,z
### self._datetime is N elements datetime instance. Index corrensponding to _darr
[docs] def setFromUri(self, baseurl=None):
if baseurl is None:
rc = Rc()
baseurl = rc.get('cy1orb', 'moongsebaseuri')
self._baseurl=baseurl
setFromDefaultUri()
[docs] def getPosition(self, start=datetime.datetime.min, stop=datetime.datetime.max):
''' Returns array of the position (Nx4 array) of the Moon in GSE system.
@param[in] start Start time of the data retrieval. Instance of datetime.datetime
@param[in] stop Stop time of the data retrieval. Instance of datetime.datetime
@retval Nx4 array of the position of the Moon in the GSE coordinates. N is the number of data, 4 is [jd,x,y,z]
'''
val=[]
for i in range(len(self._datetime)):
if self._datetime[i] >= start and self._datetime[i] <= stop:
val.append(self._darr[:,i])
return numpy.array(val).transpose().copy()
[docs] def getGse(self, t):
''' Get a position of the Moon in the GSE coordinates at the specfied t.
@param t Time of the inquiry. Time resolution is only 1 hr.
'''
t_dt = convert(t, datetime.datetime) # Convert to datetime obj
logging.debug('getGse: t=%s'%t_dt)
onehr=datetime.timedelta(0, 4000)
t0 = t_dt-onehr # Data inquire for the datebase before 1 hr
t1 = t_dt+onehr # Data inquire for the datebase after 1 hr
pos = self.getPosition(start=t0, stop=t1) # Position list
logging.debug(pos.shape)
if(len(pos)==0): # Postion not obtained
raise RuntimeError('Failed to load lunar GSE data for time %s.'%t_dt)
ndat=len(pos[0,:]) # Number of data available
logging.debug('Ndata = %d'%ndat)
delta=[abs(convert(Julday(juld), datetime.datetime)-t_dt) for juld in pos[0,:]] # Time difference between the database and the specified.
f_delta=numpy.array([delta[i].days*86400+delta[i].seconds+delta[i].microseconds*1e-6 for i in range(len(delta))]) # Time difference in seconds between the database and the specified.
pos=pos[1:4, f_delta.argmin()] # Get the nearest point to the specified.
logging.debug(pos)
return (pos[0], pos[1], pos[2])
if __name__=='__main__':
logging.basicConfig(level=logging.DEBUG)
m=MoonGse()
# m=MoonGse()
m.setFromDefaultUri();
dat=m.getPosition()
print(dat.shape)
dat=m.getPosition(start=datetime.datetime(2009,3,10), stop=datetime.datetime(2009,3,11))
print(dat.shape) # Only 25 data is fetched.
print(m.getGse(datetime.datetime(2009,1,25,9,59)))
print(m.getGse(datetime.datetime(2009,1,25,10)))
print(m.getGse(datetime.datetime(2009,1,25,10,1)))