nanograv / PINT

PINT is not TEMPO3 -- Software for high-precision pulsar timing
Other
121 stars 101 forks source link

Consider the jplephem package to replace the dependency on PySPICE #28

Closed paulray closed 8 years ago

paulray commented 9 years ago

A possible lighter-weight alternative to PySPICE is the python-jplephem package. This needs to be evaluated to see if it meets PINT's needs.

paulray commented 9 years ago

An alternative was suggested in #45, called SpiceyPy.

Both of these should be evaluated and presented at a PINT telecon at some point.

scottransom commented 9 years ago

I've started investigating the python version of NOVAS from USNO. It looks very high-quality and has high-precision routines (using two doubles for time, just like Astropy time objects) to interpolate the JPL ephemerides. It can also be installed via "pip install novas"! So I think huge wins all around. I'll investigate making a full PySPICE replacement with it over the next 2 weeks while I'm at Caltech.

http://aa.usno.navy.mil/software/novas/novas_info.php

paulray commented 8 years ago

NOVAS may be the best choice for PINT, but I thought I'd point out this development using jplepehem done for astropy: https://github.com/astropy/astropy/pull/4746

paulray commented 8 years ago

astropy v1.2 has now been released and includes code that requires jplephem, so maybe we can use astropy now and remove the PySPICE dependency? See: http://docs.astropy.org/en/stable/coordinates/solarsystem.html#astropy-coordinates-solarsystem

luojing1211 commented 8 years ago

@paulray I got this error when trying the example. Do you have the same problem. I am in version '1.2.dev15837'

>>>with solar_system_ephemeris.set('builtin'):
...     jup = get_body('jupiter', t, loc)
Traceback (most recent call last):
  File "<stdin>", line 2, in <module>
  File "astropy/coordinates/solar_system.py", line 346, in get_body
    obsgeovel=obsgeovel))
  File "astropy/coordinates/baseframe.py", line 849, in transform_to
    return trans(self, new_frame)
  File "astropy/coordinates/transformations.py", line 915, in __call__
    curr_coord = t(curr_coord, curr_toframe)
  File "astropy/coordinates/transformations.py", line 706, in __call__
    res = self.func(fromcoord, toframe)
  File "astropy/coordinates/builtin_frames/icrs_cirs_transforms.py", line 151, in icrs_to_gcrs
    gcrs_ra, gcrs_dec = atciqz(i_ra, i_dec, astrom)
  File "astropy/coordinates/builtin_frames/utils.py", line 240, in atciqz
    pco = erfa.s2c(rc, dc)
  File "astropy/_erfa/core.py", line 24440, in s2c
    stat_ok = _core._s2c(it)
AttributeError: 'module' object has no attribute '_s2c'
paulray commented 8 years ago

@luojing1211 it works OK for me. I'm using version 1.2

In [1]: import astropy

In [2]: astropy.__version__
Out[2]: u'1.2'

In [3]: from astropy.time import Time

In [4]: from astropy.coordinates import solar_system_ephemeris, EarthLocation

In [5]: from astropy.coordinates import get_body_barycentric, get_body, get_moon

In [6]: t = Time("2014-09-22 23:22")

In [7]: loc = EarthLocation.of_site('greenwich') 

In [8]: with solar_system_ephemeris.set('builtin'):
   ...:     jup = get_body('jupiter', t, loc) 
   ...:     

In [9]: jup
Out[9]: 
<SkyCoord (GCRS: obstime=2014-09-22 23:22:00.000, obsgeoloc=[ 3949481.6898847   -550931.9118969   4961151.73733442] m, obsgeovel=[  40.1745933   288.00078051   -0.        ] m / s): (ra, dec, distance) in (deg, deg, AU)
    (136.91116201, 17.02935408, 5.94386022)>
demorest commented 8 years ago

I just tried this with a fresh install of astropy 1.2.1 and it worked fine for me. Now that astropy is more stable I think it would be good to try to stick with the actual release versions (rather than dev versions) for testing with PINT.

luojing1211 commented 8 years ago

I changed my version of astropy to 1.2.1. It is working. But the I realize some problems:

>>> from astropy.time import Time
>>> from astropy.coordinates import solar_system_ephemeris, EarthLocation
>>> from astropy.coordinates import get_body_barycentric, get_body, get_moon
>>> t = Time("2014-09-22 23:22")
>>> get_body_barycentric('jupiter', t, ephemeris='de430')
<CartesianRepresentation (x, y, z) in km
    (-470381070.15935105, 579865017.4584614, 259985311.71687874)>
>>> get_body_barycentric('jupiter', t, ephemeris='de421')
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/Users/jingluo/anaconda/lib/python2.7/site-packages/astropy/coordinates/solar_system.py", line 207, in get_body_barycentric
    kernel = _get_kernel(ephemeris)
  File "/Users/jingluo/anaconda/lib/python2.7/site-packages/astropy/coordinates/solar_system.py", line 174, in _get_kernel
    return SPK.open(download_file(value, cache=True))
  File "/Users/jingluo/anaconda/lib/python2.7/site-packages/astropy/utils/data.py", line 1030, in download_file
    remote_url, timeout=timeout)) as remote:
  File "/Users/jingluo/anaconda/lib/python2.7/urllib2.py", line 154, in urlopen
    return opener.open(url, data, timeout)
  File "/Users/jingluo/anaconda/lib/python2.7/urllib2.py", line 421, in open
    protocol = req.get_type()
  File "/Users/jingluo/anaconda/lib/python2.7/urllib2.py", line 283, in get_type
    raise ValueError, "unknown url type: %s" % self.__original
ValueError: unknown url type: de421

It seems astropy can not find where to download de421. The old version of ephemeris is stored at http://naif.jpl.nasa.gov/pub/naif/generic_kernels/spk/planets/a_old_versions/ Which is one subdirectory of the current planets ephemeris. Maybe three line of astropy code change will make a difference.

We found out that the old version ephemeris can be accessed from url of those files. So these statement is not really true now.

luojing1211 commented 8 years ago

Hi all, I got reply from astropy solar system developer Marten. The plan for astropy 1.3 would have the velocity calculation provided. Actually the current development version already has a routine to do so.

I did a test of comparing the result with astropy and spice. The position difference is at around 1e-2 m/cm level, the velocity is even less. In addition, Astropy solar system posvel API takes an array of time, which is not the case for current PINT.

In conclusion, I think we should really consider to switch to astropy solar system module for PINT. By doing so, we can reduce our dependency which leads to an easy and clean installation. Secondly, it would potentially increase our computing speed.

demorest commented 8 years ago

Closing since this change was implemented in PR #184