If you're not doing it already (didn't see it anywhere, but didn't look that hard), incorporate the SOFA libraries into skyfield:

http://www.iausofa.org/

Description (from site):

The International Astronomical Union's SOFA service has the task of establishing and maintaining an accessible and authoritative set of algorithms and procedures that implement standard models used in fundamental astronomy. The service is managed by an international panel, the SOFA Board, appointed through IAU Division A — Fundamental Astronomy. SOFA also works closely with the International Earth Rotation and Reference Systems Service (IERS).

Hi, what is the usage of the hipparcos.py functions?

I see from the tests that I can do parsing of hipparcos lines.

line = 'H|       11767| |02 31 47.08|+89 15 50.9| 1.97|1|H|037.94614689|+89.26413805| |   7.56|   44.22|  -11.74|  0.39|  0.45|  0.48|  0.47|  0.55|-0.16| 0.05| 0.27|-0.01| 0.08| 0.05| 0.04|-0.12|-0.09|-0.36|  1| 1.22| 11767| 2.756|0.003| 2.067|0.003| | 0.636|0.003|T|0.70|0.00|L| | 2.1077|0.0021|0.014|102| | 2.09| 2.13|   3.97|P|1|A|02319+8915|I| 1| 1| | | |  |   |       |     |     |    |S| |P|  8890|B+88    8 |          |          |0.68|F7:Ib-IIv SB|G\n'
star = hipparcos.parse(line)
print(star.ra.hours, star.dec.degrees)

When I use get, I can download the hipparcos file. I think there needs to be another check for whether the file was downloaded. By default, when I use load, I get an error with TypeError: coercing to Unicode: need string or buffer, file found. If I instead use the following for the already downloaded file, I can open it without error.

gzip_file_handle = gzip.open('hip_main.dat.gz')
data = gzip_file_handle.readlines()

The above uses the shorthand of gzip.open(f). If I substitute gzip.GzipFile(fileobj=f), then I am able to read the lines of the file.

The New Cross Index (IV/27A) could be used to populate the Star name. The particular instance I wanted to use this was to look up both the Bayer designation and visual magnitude using the Hipparcos number.

from skyfield.data import hipparcos
s = hipparcos.get('25930')

s.names is currently [('HIP', 25930)].

from skyfield.units import Angle

def cross_index(catalog, name):
    '''Uses the New Cross Index (`IV/27A <http://cdsarc.u-strasbg.fr/viz-bin/Cat?IV/27A>`_) to look up designations between catalogs. It also includes the visual magnitude, right ascension, magnitude and constellation abbreviation.

    Designations:

    Henry Draper Catalog Number ``HD`` (`III/135 <http://vizier.u-strasbg.fr/viz-bin/VizieR?-source=III/135>`_), 
    Durchmusterung Identification from HD Catalog ``DM`` (see `IV/27A note 1 <http://cdsarc.u-strasbg.fr/viz-bin/Cat?IV/27A#sRM3.1>`_), 
    General Catalogue of 33342 stars ``GC`` (`I/113 <http://vizier.u-strasbg.fr/viz-bin/VizieR?-source=I/113>`_), 
    Harvard Revised Number - BSC5 ``HR`` (`V/50 <http://vizier.u-strasbg.fr/viz-bin/VizieR?-source=V/50>`_), 
    Hipparcos Catalog ``HIP`` (`I/196 <http://vizier.u-strasbg.fr/viz-bin/VizieR?-source=I/196>`_)

    You can look solely for the Bayer ``Bayer`` *or* Flamsteed number ``Fl``, or return either with ``BFD``.

    `IV/27A note 2 <http://cdsarc.u-strasbg.fr/viz-bin/Cat?IV/27A#sRM3.2>`_ on HIP and CSI (`IV/9 <http://vizier.u-strasbg.fr/viz-bin/VizieR?-source=IV/9>`_) right ascensions and visual magnitude:

        Right ascensions, declinations and visual magnitudes for all stars were taken from the Hipparcos catalog and from the CSI for the stars that has no number in catalog Hipparcos.

    Here is an example of looking up a star in Orion's belt by using the Hipparcos number:

    ::

        sn_dict = cross_index('HIP','25930')
        print(sn_dict['Bayer'],sn_dict['Vmag'],sn_dict['Cst'])

    The greek alphabet dictionary was appreciatively taken from Pystaratlas' `catalogues.py <https://code.google.com/p/pystaratlas/source/browse/catalogues.py>`_.

    Excluding DE to Dec as a key abbreviation and the addition of epoch, the labels are the same as the New Cross Index.
    '''

    alphabet={'alf  ':'α','alf01':'α¹','alf02':'α²',\
          'bet  ':'β','bet01':'β¹','bet02':'β²','bet03':'β³',\
          'gam  ':'γ','gam01':'γ¹','gam02':'γ²','gam03':'γ³',\
          'del  ':'δ','del01':'δ¹','del02':'δ²','del03':'δ³',\
          'eps  ':'ε','eps01':'ε¹','eps02':'ε²','eps03':'ε³',\
          'zet  ':'ζ','zet01':'ζ¹','zet02':'ζ²','zet03':'ζ³',\
          'eta  ':'η','eta01':'η¹','eta02':'η²','eta03':'η³',\
          'the  ':'θ','the01':'θ¹','the02':'θ²','the03':'θ³',\
          'iot  ':'ι','iot01':'ι¹','iot02':'ι²',\
          'kap  ':'κ','kap01':'κ¹','kap02':'κ²',\
          'lam  ':'λ','lam01':'λ¹','lam02':'λ²','lam03':'λ³',\
          'mu.  ':'μ','mu.01':'μ¹','mu.02':'μ²','mu.03':'μ³',\
          'nu.  ':'ν','nu.01':'ν¹','nu.02':'ν²','nu.03':'ν³',\
          'ksi  ':'ξ','ksi01':'ξ¹','ksi02':'ξ²','ksi03':'ξ³',\
          'omi  ':'ο','omi01':'ο¹','omi02':'ο²','omi03':'ο³',\
          'pi.  ':'π','pi.01':'π¹','pi.02':'π²','pi.03':'π³','pi.04':'π⁴','pi.05':'π⁵','pi.06':'π⁶',\
          'rho  ':'ρ','rho01':'ρ¹','rho02':'ρ²','rho03':'ρ³',\
          'sig  ':'σ','sig01':'σ¹','sig02':'σ²','sig03':'σ³',\
          'tau  ':'τ','tau01':'τ¹','tau02':'τ²','tau03':'τ³','tau04':'τ⁴','tau05':'τ⁵','tau06':'τ⁶','tau07':'τ⁷','tau08':'τ⁸','tau09':'τ⁹',\
          'ups  ':'υ','ups01':'υ¹','ups02':'υ²','ups03':'υ³',\
          'phi  ':'φ','phi01':'φ¹','phi02':'φ²','phi03':'φ³','phi04':'φ⁴',\
          'chi  ':'χ','chi01':'χ¹','chi02':'χ²','chi03':'χ³',\
          'psi  ':'ψ','psi01':'ψ¹','psi02':'ψ²','psi03':'ψ³','psi04':'ψ⁴','psi05':'ψ⁵','psi06':'ψ⁶','psi07':'ψ⁷','psi08':'ψ⁸','psi09':'ψ⁹',\
          'ome  ':'ω','ome01':'ω¹','ome02':'ω²','ome03':'ω³'}

    cross_index_url = 'ftp://cdsarc.u-strasbg.fr/pub/cats/IV/27A/catalog.dat'
    star_dict = {}

    #TODO: if used in skyfield load/cache can be used instead of wget. https://github.com/brandon-rhodes/python-skyfield/blob/master/skyfield/data/hipparcos.py and https://github.com/brandon-rhodes/python-skyfield/blob/master/skyfield/io.py

    try:
        data = open('catalog.dat','r')
    except:
        import wget
        wget.download(cross_index_url)
        data = open('catalog.dat','r')

    for l in data.readlines():
        s = {}
        s['HD']   = l[0:6]
        s['DM']   = l[7:19]
        s['GC']   = l[20:25]
        s['HR']   = l[26:30]
        if(s['HR'] == '    '): s['HR'] = None
        s['HIP']  = l[31:37]
        if(s['HIP'] == '      '): s['HIP'] = None
        ra = float(l[38:40])+float(l[40:42])/60.+float(l[42:47])/3600.
        s['RA']   = Angle(degrees=float(ra))
        if(l[48]=='+'): 
            sign=1
        else:
            sign=-1
        dec = sign*float(l[49:51])+float(l[51:53])/60.+float(l[53:57])/3600.
        s['Dec']  = Angle(degrees=float(dec))
        s['Vmag'] = l[58:63]
        s['Fl']   = l[64:67]
        s['Bayer']= l[68:73]
        if(s['Bayer'] == '     '):
            s['Bayer'] = None
        else:
            try:
                s['Bayer'] = alphabet[s['Bayer']]
            except:
                pass
        if(s['Bayer']):
            s['BFD'] = s['Bayer']
        elif(s['Fl']):
            s['BFD'] = s['Fl']
        else:
            s['BFD'] = None
        s['Cst']  = l[74:77]
        s['epoch'] = 2000 # The New Cross Index (IV/27A) uses epoch 2000.
        # This function could be modified to look at a range of magnitudes, positions, or by constellation.
        if(s[catalog] is not None and s[catalog].lower().strip() == name.lower().strip()):
            return s
    return False

Dictionary of names returned that matches the Hipparcos number 25930

n = cross_index('HIP','25930')
print(n)
print(n['Bayer'])

And the Bayer designation '\xce\xb4' = δ.

{'HIP': ' 25930', 'DM': 'BD-00   983 ', 'Vmag': ' 2.25', 'HR': '1852', 'BFD': '\xce\xb4', 'epoch': 2000, 'GC': ' 6847', 'RA': <Angle 05deg 32' 00.4">, 'Cst': 'Ori', 'Bayer': '\xce\xb4', 'Dec': <Angle 00deg 17' 56.7">, 'Fl': ' 34', 'HD': ' 36486'}
δ

s.names = n
print(s.names)
print(s.names['Bayer'])

{'HIP': ' 25930', 'DM': 'BD-00   983 ', 'Vmag': ' 2.25', 'HR': '1852', 'BFD': '\xce\xb4', 'epoch': 2000, 'GC': ' 6847', 'RA': <Angle 05deg 32' 00.4">, 'Cst': 'Ori', 'Bayer': '\xce\xb4', 'Dec': <Angle 00deg 17' 56.7">, 'Fl': ' 34', 'HD': ' 36486'}
δ

My apologies for all the issues this week. I've just been integrating Skyfield into a project and it has been useful enough that I've read through the code a few times :).

What do you think of instead of setting names=[('HIP','25930')], the format could be a dictionary with the catalog abbreviation as the key? Or of course left as is. The function could possibly be used as a standalone method, or in a Star class method. If used internally with the Hipparcos functions, since position is already populated, the values of RA and Dec should be popped out so that there arent duplicates like in the above example.

If this fits in with what you want for Skyfield, let me know and I can make a pull request for it.

Is there a plan to support conversion to/from galactic coordinates?
e.g. in PyEphem I could monitor the alt/az of the galactic center at my observing location (see code below).
In the short-term, I can use PyEphem to generate RA/Dec of various locations in galactic coodinates, and then use skyfield the rest of the way. But it would be great to have a Galactic <--> equatorial utility within skyfield.

eq00 = ephem.Equatorial(ephem.Galactic(0,0))
print 'RA:', eq00.ra, 'dec:', eq00.dec

g00 = ephem.FixedBody()
g00._ra = eq00.ra
g00._dec = eq00.dec
g00._epoch = eq00.epoch

haystack = ephem.Observer()
haystack.lat = '42.623220'
haystack.lon = '-71.488210'
haystack.date = '2014/10/10 13:00:00'  #UTC

g00.compute(haystack)
rToD = 180./math.pi
print "%20s  %7.2f %7.2f" % (ephem.localtime(haystack.date).ctime(), g00.alt*rToD, g00.az*rToD)

setup.py uses the install_requires param to define the dependencies, but that argument is ignored by distutils.core.setup. The script should import and use setuptools.setup.

Brandon, not sure if you want to add this to skyfield, but it might be
useful to know the geocentric position of a given observing site
(latitude/longitude/elevation) to high precision for a given
time. This would be an alternate way of compensating for
precession. Thoughts:

• The position of the North Pole at the J2000 epoch is, of course,
  (0,0,epr) by definition, where epr is the Earth's polar radius
  (6356.752 km per HORIZONS)
• I used HORIZONS to find the position of the North Pole from 9999
  BC to 9999 AD. HORIZONS won't let you compute the North Pole from
  the geocenter, but you can do the opposite: compute the geocenter
  from the north Pole (so the following values should be negated):

https://github.com/barrycarter/bcapps/blob/master/ASTRO/north-pole-from-geocenter.txt.bz2

• As expected, all 3 coordinates (x,y,z) follow a sinusoidal pattern
  over time with a period matching the Earth's precession. I
  considered retrofitting the data to sine waves, but I'm sure you can
  use a more accurate precession model.
• The position of non-polar locations of course changes as the Earth
  rotates, and repeats every sidereal day. However, it should follow a
  reasonably simple formula, even allowing for the Earth's ellipticity
  and precession tilt.
• You may or may not find other helpful files at:

https://github.com/barrycarter/bcapps/tree/master/ASTRO/

It would be useful at some point to have a method to calculate apparent angular separation.

I'm calculating solar eclipse visibility on the earth and on planes, and to track something like "path of totality" I'm minimizing the apparent separation between moon and sun.

mpos = place.at(jd).observe(moon).apparent().position.km
spos = place.at(jd).observe(sun).apparent().position.km
mlen = np.sqrt((mpos**2).sum())
slen = np.sqrt((spos**2).sum())
sepa = ((3600.*180./np.pi) *
        np.arccos(np.dot(mpos, spos)/(mlen*slen)))  # arcsec

I've asked the question here:
http://stackoverflow.com/questions/36111098/better-way-to-calculate-apparent-angular-separation-of-two-objects-in-skyfield
and discussed minimization here:
http://stackoverflow.com/questions/36110998/why-does-scipy-optimize-minimize-fail-when-calling-a-skyfield-based-function

Does skyfield provide position data for non-Terran satellites, in particular the 4 large Jovian satellites, Titan, and other non-Terran satellites "easily" visible from Earth?

On http://rhodesmill.org/skyfield/positions.html the following code doesn't work:

from skyfield.api import earth, now
sat = earth.satellite(tle_text)
position = sat(now())

Error is TypeError: 'EarthSatellite' object is not callable.

How do I calculate the GCRS (xyz) position of a satellite?

Good morning. I've found this project after digging through stackoverflow.

I'm just wondering how mature the code you recommend using is for converting ECI to ECEF?

The reason for me asking is that I've been working on a little project outside of work to track one of our cubesats, UKube-1. I want to calculate a dozen or so positions of the satellite based off its TLEs and draw a line on a map to project the satellite's position.

Many thanks for all your work,
Colin

Currently the docs contain an API reference section, but it's almost empty.

It would be great to have API docs available so that one can explore what is available in skyfield and what the inputs / outputs to each function / method are.

I wanted to make a PR for #81, but running py.test as per HACKING fails on current master with:

$ py.test skyfield/tests/test_timelib.py
========= test session starts ============================
platform darwin -- Python 3.5.0, pytest-2.9.0, py-1.4.31, pluggy-0.3.1
[...]
collected 22 items 

skyfield/tests/test_timelib.py EEEEEEEEEEEEEEEEEEEEEE

==================== ERRORS =======================
_____ ERROR at setup of test_JulianDate_init ________________
file /Users/amici/devel/python-skyfield/skyfield/tests/test_timelib.py, line 17
  def test_JulianDate_init(time_parameter, time_value):
        fixture 'time_parameter' not found
        available fixtures: capfd, cache, monkeypatch, spawnu, recwarn, tmpdir_factory, pytestconfig, tmpdir, capsys, record_xml_property
        use 'py.test --fixtures [testpath]' for help on them.

/Users/amici/devel/python-skyfield/skyfield/tests/test_timelib.py:17
[...]

Tests don't look right for use with py.test, all fixtures miss the @pytest.fixture decorator, I'd like to help, but I'm pretty confused.

I am trying to gauge the error in underlying sgp4 model by comparing derived satellite positions to those predicted by sgp4. I have written a little script to do this comparisons and it seems to work beautifully using the PyEphem package but not with Skyfield. Below is the output of the script which is available here

--------------------------------------------------------------------------------
Time diff between position measurement and TLE epoch (s): -2.984833
--------------------------------------------------------------------------------
----------------------------------- pyephem ------------------------------------
Total error in position sqrt(xd^2+yd^2+zd^2)(m): 444.325703545
Difference between observered and computed radius (m): -9.31322574615e-10
--------------------------------------------------------------------------------
----------------------------------- skyfield -----------------------------------
Total error in position in GCRS sqrt(xd^2+yd^2+zd^2)(m): 1891045.97619
Difference between observered and computed radius (m): 71.4918435588
Total error in position in ITRF sqrt(xd^2+yd^2+zd^2)(m): 34642.4471312
Difference between observered and computed radius (m): 71.4918435551
Total error in position in ITRF + polar sqrt(xd^2+yd^2+zd^2)(m): 34642.6084994
Difference between observered and computed radius (m): 71.4918435561
----------------------------------- skyfield -----------------------------------

For PyEphem the error is on the order of 0.5 km at epoch which from what I understand is normal but from Skyfield the error is 34 km? I hope I'm just doing something wrong that I can see or just don't understand some aspect of the code use?

But if not I think there are two options going forward to try and track this issue down:

1. See if there is a way to get the prediction out of PyEphem (libastro) in TEME system and compare directly to SPG4 output bypass skyfield.
2. Implement the transform from TEME to ITRF myself and see if they match then.

Any help or suggestion of which or an alternative approach would be great.

Thanks

The 0xc2 character in README.rst prevents building of the package under python 3.4.

The patch to fix it can be find at
https://github.com/jochym/conda-rcp/blob/master/skyfield/fix_unicode.patch

P.

I would have expected NumPy arrays to display using the same formatting and the same number of significant digits regardless of the version of Python that was being used to run NumPy, but it looks like Travis is failing us only for Python 3.4 because of a number of textual mismatches in our documentation:

https://travis-ci.org/brandon-rhodes/python-skyfield/jobs/27500929

Why is this — does anyone know enough about NumPy to understand the cause behind this? Or should we just turn off doctesting for Python 3.4 for now?

Is there anywhere in Skyfield that there will be less precision than a float or np.float64? I didn't see any, but I thought I'd ask anyway. I was just curious what my worst-case rounding error would be. https://en.wikipedia.org/wiki/Machine_epsilon

import skyfield
import numpy as np
#--------------------------------
print('\ndefault type is float which has same precision as numpy.float64')
from skyfield.constants import ASEC360
eps_a0 = np.finfo(ASEC360).eps
print(type(ASEC360))
print(eps_a0)
#--------------------------------
print('\nRadians Epsilon')
from skyfield.data import hipparcos
star = hipparcos.get('11767')
r_type = type(star.ra.radians)
print(r_type)
eps_f64 = np.finfo(r_type).eps
print(eps_f64)
#--------------------------------
#from skyfield.units import Angle
#eps_a = np.finfo(skyfield.units.Angle).eps

output:

default type is float which has same precision as numpy.float64
<type 'float'>
2.22044604925e-16

Radians Epsilon
'11767'
<type 'numpy.float64'>
2.22044604925e-16

After perusing the tests, I see that the distance precision (excluding Astropy or NOVAS tests) is expected to be less than 1e-15, and the time epsilon is one_second * 42.0e-6 with a note that 20.1e-6 is a theoretical best precision. Are there any other important ones I missed?

I'm evaluating the use of skyfield for an application to earth observation satellites, I understand this is not a target use case at the moment, but I easily managed to extract the satellite orbital mechanics out of EarthSatellite with a modicum of code duplication, so I'd like to know if you are interested in supporting some basic form of this use case.

As an example the EarthSatellite class compute the ITRF position of the satellite for a given time object inside _compute_GCRS, but doesn't expose it in the API:

https://github.com/skyfielders/python-skyfield/blob/master/skyfield/sgp4lib.py#L67

Could you please add a EarthSatellite.itrf method that returns the rITRF and vITRF quantities already computed?

Comparing a JulianData object to anything else than another JulianDate object crashes with:

AttributeError: <type> object has no attribute 'tt'

I would have expected to simply return False.

I didn't test it, but it seems to me that even comparing two JulianDate objects with an array of times would crash due to the fact that numpy.ndarray's do not have a defined truth value.

With the new way of loading ephemeris data, it loads the files to the "current folder", which may work in most cases but definitely not in all. I think it would be good to have another optional argument for load() where you can specify the path to a download folder. This is especially important for more production-like environments where you cannot control what the "current directory" is from which a script is run, and then the permissions might be set to prevent any creation of files in that one.

Loaded, built, and installed latest version of skyfield on Windows 7 using the Anaconda version of python. Attempting to follow the example, I issued the following import command in ipython

In [5]: from skyfield.api import load, now

ImportError Traceback (most recent call last)
in ()
----> 1 from skyfield.api import load, now

c:\Users\ast0528\Documents\Astronomy\skyfield\build\lib\skyfield\api.py in ()
10 from .constants import tau
11 from .errors import DeprecationError
---> 12 from .iokit import load
13 from .starlib import Star
14 from .timelib import JulianDate, T0, now, utc

c:\Users\ast0528\Documents\Astronomy\skyfield\build\lib\skyfield\iokit.py in ()
6 from time import time
7
----> 8 from .jpllib import SpiceKernel
9
10 try:

c:\Users\ast0528\Documents\Astronomy\skyfield\build\lib\skyfield\jpllib.py in ()
4 from numpy import max, min
5
----> 6 from jplephem.spk import SPK
7 from jplephem.names import target_names as _names
8

ImportError: No module named jplephem.spk

For my first go with Skyfield i tried running the example at the website. It results in a Permission error when loading the planets.

PermissionError: [WinError 32] The process cannot access the file because it is being 
used by another process: 'de421.bsp.download' -> 'de421.bsp'

It seems that it's trying to rename the downloaded file, while the handle is still open. This doesn't seem correct to me, maybe it a Windows thing, since I'm running Windows 7 x64. I installed Skyfield 0.6 with pip inside a Conda environment with Python 3.4 and Numpy 1.10.1.

Here is a notebook replicating the issue. I also tested if the handle really is the cause, moving the rename statement outside the with block (so closing) makes it work for me. This could be a possible solution, but i haven't inspected the entire code.

http://nbviewer.ipython.org/gist/RutgerK/811320c5222cbcfa49ee

A little side note, should Pandas be added to the requirements in setup.py?

An email correspondent tells me that Earth satellites do not accept date arrays. This should be fixed.

from skyfield.api import now, earth
jd = now()
tle = ["1 25544U 98067A   14231.33679667  .00006701  00000-0  12490-3 0  9227",
       "2 25544  51.6463 146.5620 0004567  38.0995  53.3525 15.50090279901121"]
satellite = earth.satellite("\n".join(tle))
apparent = satellite.gcrs(jd)
print apparent.radec()
# (Angle(22h 52m 14.47s), Angle(-22deg 27' 57.9"), Distance(4.54602e-05 AU))
print apparent.radec()[2].km
#6800.7482096125714

import ephem
date = ephem.Date(jd.utc_datetime()[0])
tledata = ephem.readtle('foo', tle[0], tle[1])
tledata.compute(date)
print tledata.a_ra, tledata.a_dec
#22:52:14.27 -22:27:52.6
print tledata.elevation/1000+6378.16
#6800.774875

By the way, have you though about using named tuples? Then it would be possible to do apparent.radec().ra for instance.

Created an array of dates:
when = JulianDate(utc=(2015,12,range(1,31))

Compute location of moon on each of those dates and convert to ra/dec
ra,dec,dist = moon.at(when).radec()

Print dec as dms string:

print dec.dstr()
31 values from +21deg 39' 13.5" to +23deg 24' 59.8"

Print ra as hms string:

print ra.hstr()
C:\Users\ast0528\AppData\Local\Continuum\Anaconda\lib\site-packages\skyfield\units.pyc in hstr(self, places, warn)
    207         if getattr(hours, 'shape', None):
    208             return "{0} values from {1} to {2}".format(
--> 209                 len(degrees),
    210                 _hstr(min(degrees),places,signed),
    211                 _hstr(max(degrees),places,signed)

NameError: global name 'degrees' is not defined

It appears the problem is in units.py in the hstr function. Instead of referencing 'degrees', it should have been 'hours'. Note the signed argument needs to be removed as well.

        if getattr(hours, 'shape', None):
            return "{0} values from {1} to {2}".format(
                len(hours),
                _hstr(min(hours),places),
                _hstr(max(hours),places)
            )

What is the best way to represent arbitrary positions in conjunction with Skyfield?

Here is an example I am thinking about. It takes the positions of a few stars, and averages the ra/dec/dist parameters. This, without any physics basis, finds the mean position; a visual barycenter I am thinking.

import skyfield
from skyfield.api import Star, earth, now
from skyfield.units import Angle, Distance
from numpy import array

# Gather star positions and put coords in an array
algol = Star(ra_hours=( 3,  8, 10.1315),  dec_degrees=(40, 57, 20.332))  #approximately Algol
mizar = Star(ra_hours=(13, 23, 55.5),     dec_degrees=(54, 55, 31))      #approximately Mizar
vega  = Star(ra_hours=(18, 36, 56.33635), dec_degrees=(38, 47, 01.2802)) #approximately Vega
stars = [algol, mizar, vega]
positions = []
for star in stars:
    star_pos = earth(now()).observe(star)
    positions.append(star_pos.radec())
position_array = array(positions)
print('Array of Posititions')
print(position_array)

# Average the coordinates to get an unweighted center
mean_ra = skyfield.units.Angle(hours=0.0).radians
mean_dec = skyfield.units.Angle(degrees=0.0).radians
mean_dist = skyfield.units.Distance(AU=0.0).AU
for coord in positions:
    mean_ra = mean_ra + coord[0].radians
    mean_dec = mean_dec + coord[1].radians
    mean_dist = mean_dist + coord[2].AU
mean_ra = mean_ra / len(position_array)
mean_dec = mean_dec / len(position_array)
mean_dist = mean_dist / len(position_array)

center = (Angle(radians=mean_ra, preference='hours'), Angle(radians=mean_dec, signed=True), Distance(AU=mean_dist))
print('Center of Array of Positions')
print(center) # What is the best way to express a position with Skyfield when it is for a made-up object?
#print(center.position.AU)

# Make a fake star to show getting position.AU attribute
# A downside to this is possibly can't specify distance?
fake_star = Star(center[0],center[1])
test_center = earth(now()).observe(fake_star)
print('Fake Star\'s Position')
print(test_center.radec())
print(test_center.position.AU)

Output:

Array of Positions
[[<Angle 03h 08m 10.13s> <Angle +40deg 57' 20.3"> <Distance 2.06265e+14 AU>]
 [<Angle 13h 23m 55.50s> <Angle +54deg 55' 31.0"> <Distance 2.06265e+14 AU>]
 [<Angle 18h 36m 56.34s> <Angle +38deg 47' 01.3"> <Distance 2.06265e+14 AU>]]
Center of Array of Positions
(<Angle 11h 43m 00.66s>, <Angle +44deg 53' 17.5">, <Distance 2.06265e+14 AU>)
Fake Star's Position
(<Angle 11h 43m 00.66s>, <Angle +44deg 53' 17.5">, <Distance 2.06265e+14 AU>)
[ -1.45734220e+14   1.08229331e+13   1.45566382e+14]

I see radec() returns a tuple. How do I express in Astronometic (and other?) position types in addition to this without actually creating a Star?

print(type(center))
print(type(test_center))
print(type(test_center.radec()))

Output:

<type 'tuple'>
<class 'skyfield.positionlib.Astrometric'>
<type 'tuple'>

I'm using sgp4 to get TEME coordinates. I could use skyfield to convert these to GCRS. My problem is that this adds a massive dependency to my own library, caused by the de421 dependency of skyfield which is 26MB.

The main use case of skyfield probably requires de421 or similar to do its work, and the satellite ephemeris calculations are a side use case I guess. So removing de421 as dependency and letting the user install such package himself may not be what you want (loosing out-of-the box experience). On the other hand, letting the user have more freedom in choosing the package he wants might be good as well. After all, it's just pip install skyfield de421 or pip install skyfield de423 or similar.

What do you think?

I wanted to address the question "How does the Earth move in the sky as seen from the Moon?" (http://astronomy.stackexchange.com/q/8100/7982) but body.topos() failed for bodies other than the earth (e.g. moon, mars). Is there a way to access at least the librations of the moon within Skyfield?

IDEAL: implement .topos() for "popular" solid bodies (e.g. moon, mars) Mars may not be possible because Mars rotation information is not part of the ephemerides, but lunar nutation is in most of them.

HELPFUL: at least a method to access the lunar libration information in the ephemeris when it is explicitly present.

Edit: at a bare minimum, at least give more info when the method fails on bodies that are not earth.

Question: is the majority of the effect (motion of Earth's apparent position in the lunar sky from lunar surface) just the result of the constant lunar rotation vs "non-constant" motion of the moon in it's elliptical orbit? ln that case I could approximate the motion of earth's position in the lunar sky for illustration just by using simple geometry and constant lunar rotation (e.g. Wikipedia's 27.321661 days).

import numpy as np
import matplotlib.pyplot as plt
from skyfield.api import load, now, JulianDate
import json

data = load('de405.bsp')  # de421 also tested

jd = JulianDate(utc=(2016,1,1))

earth = data['earth']
moon  = data['moon']

earth_0_0 = earth.topos(0.0, 0.0)
earth_pos = earth_0_0.at(jd).position.km
print "earth_pos: ", earth_pos

moon_0_0  = moon.topos(0.0, 0.0)    # also fails with ("0.0 N", "0.0 E")
moon_pos  = moon_0_0.at(jd).position.km
print "moon_pos: ", moon_pos

gives:

earth_pos:  [ -2.43910085e+07   1.33218189e+08   5.77225486e+07]
Traceback (most recent call last):
  File "moon_topos.py", line 18, in <module>
    moon_0_0  = moon.topos(0.0, 0.0)    # also fails with ("0.0 N", "0.0 E")
  File ".../anaconda2/anaconda/envs/scipyenv/lib/python2.7/site-packages/skyfield/jpllib.py", line 114, in topos
    assert self.code == 399
AssertionError

from skyfield.api import earth
tle = str.join("\n", [
    '1 99999U 14360DT  15001.50000000  .00000000  00000-0  00000-0 0  0010',
    '2 99999  90.0000 000.0000 0000000 000.0000  00.0000 15.00000000 00013'
])
print(earth.satellite(tle).epoch.utc_datetime())

Instead of "2015-01-01 12:00:00+00:00", I get "2014-12-31 12:00:00+00:00"

I have confirmed that this bug exists in 0.4. Installation of some dependencies is failing for me, so I do not know if this bug has been fixed in the master branch.

pyramid has a nice decorator called @reify to replace a method with its value once called:

http://docs.pylonsproject.org/projects/pyramid/en/latest/api/decorator.html

check the source, there is a link in the top right corner.

Seems like supporting fixed bodies (stars, extragalactic bodies, etc) would be pretty easy to implement and would be great functionality to have.

I think it can be confusing both for basic Python people like me and numpy users to see “max” in code examples actually be “numpy.max”. The scientific Python community defines a convention of importing numpy as np, so that using “np.min” has the advantages of being namespaced but still unambiguous.

The class that represents a geocentric position should gain a method by which the user can ask which geographic location is directly beneath the position, in latitude and longitude. It should probably, in fact, return a topos object!

Hi Brandon,

I would like to be obtain moon phase information similar to what PyEphem provides:

ephem.next_full_moon
ephem.Moon().phase

and similar please.

Thanks,

Bernard.

Downloading entire ephemerides takes up a lot of time and space. As it turns out, NASA's servers support HTTP/1.1 byte range requests. For example:

curl -o test.txt -r 120000000-130000000 ftp://ssd.jpl.nasa.gov/pub/eph/planets/ascii/de431/ascm06000.431

yields 10,000,001 bytes from the specified position in ascm06000.431.

Since these files (and the xsp and bsp files) are very well structured, it should be possible to pull (and cache) data on an as-needed basis, instead of having to download entire ephemerides.

This would be particularly useful when computing non-Terran satellite positions, whose ephemerides can be very large.

This has been separated from #79 since it only requires data which is already present in (most of) the ephemerides.

To generate a view of the disk of the moon from earth (not only the position), the libration data in the ephemerides is useful.

Is it possible to access that now through Skyfield? If not, a new method would be helpful - possibly something like body.at().position.libration() or body.at().libration(), giving an error when body isn't the moon [301].

In the mean time: - is the major contribution to the libration just the difference between the steady rotation of the moon and the non-steady increase in true anomaly due to the elliptical orbit? If so, I can approximate the libration for short periods of time just using geometry and a rotation period.

DE430 is fairly large and doesn't cover even large non-Earth satellites or major asteroids such as Ceres. Consider generating SPK files for "important" bodies on a year-to-year basis using http://ssd.jpl.nasa.gov/x/spk.html

Instead of downloading 100 years worth of data for a FEW bodies, download 1 years worth of data for MANY bodies.

Possibly even make skyfield intelligent enough to retrieve these files "on the fly" as needed (or some combination: eg, provide major planet files for 10 years and download on-the-fly for anything else the user requests).

I recently reinstalled OSX 10.11.3 and some pythons (IDLE, Anaconda). When I download DE405 from the Anaconda installation from a terminal, it takes a few seconds as always. However, when I try to do the same from my IDLE interface it starts throttling down to practically zero. After waiting 30 minutes, IDLE reached nearly 100% then gave the output shown below. The shell is frozen and I need to Force Quit from the OS. The download is not complete - a large, incomplete download file is shown in Finder.

I've never tried to use IDLE for Skyfield before so I can't say if this is new or not. Besides "don't use IDLE", any other suggestions? I'm OK with other pythons, but for some reason I like working in IDLE - the default UI downloaded with python at phython.org.

Ran this script:

from skyfield.api import load
from skyfield.positionlib import ICRF

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

data  = load('de405.bsp')

Starts fast, then starts to throttle down after roughly 30% (wasn't watching closely), so that after 30 minutes, the following happens. The lines started printing at about 2 per second. At 50% they slow to 1 per second. However, here near the end, the lines were printing at about 1 per minute.

[################################ ]  97% de405.bsp
[################################ ]  98% de405.bsp
[################################ ]  98% de405.bsp
[################################ ]  98% de405.bsp
[################################ ]  98% de405.bsp
[################################ ]  98% de405.bsp
[#################################]  99% de405.bsp
[#################################]  99% de405.bsp
[#################################]  99% de405.bsp

Traceback (most recent call last):
  File "/Users/me/Documents/fishing/Cosmo/make solar system v00.py", line 8, in <module>
    data  = load('de405.bsp')
  File "/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/skyfield/iokit.py", line 71, in __call__
    return load(filename, self.directory, verbose=self.verbose)
  File "/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/skyfield/iokit.py", line 176, in load
    download(url, path, verbose=verbose)
  File "/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/skyfield/iokit.py", line 248, in download
    raise IOError('error getting {0} - {1}'.format(url, e))
IOError: error getting ftp://ssd.jpl.nasa.gov/pub/eph/planets/bsp/de405.bsp - [Errno 54] Connection reset by peer
>>>

The ICRS class contains an undocumented function ecliptic_latlon which is supposed to return ecliptic latitude, longitude, and distance.

However, the results of this function don't appear plausible to me.

Here is an example, plotting the 1900-2049 orbit of Pluto projected in a plane perpendicular to the ecliptic:

import numpy as np
import matplotlib.pyplot as plt

from skyfield.api import pluto, JulianDate
from skyfield.functions import from_polar

jd = JulianDate(utc=(range(1900,2050),))
pluto_orbit = pluto(jd)

# top: using ecliptic_latlon - wrong?
plt.subplot(211)
lat, lon, r = pluto_orbit.ecliptic_latlon()
pos = from_polar(r.AU, lat.radians, lon.radians)
plt.scatter(pos[0], pos[1], linewidths=(0,), marker='.')

# bottom: rotating it myself - right?
plt.subplot(212)
from skyfield.constants import ASEC2RAD, T0
from skyfield.functions import rot_x
from skyfield.nutationlib import mean_obliquity

epsilon = mean_obliquity(T0) * ASEC2RAD
rot_to_ecliptic = rot_x(epsilon)

pos = rot_to_ecliptic.dot(pluto_orbit.position.AU)
plt.scatter(pos[1], pos[2], linewidths=(0,), marker='.')

plt.show()

Here is the output:

In the top plot (created using ecliptic_latlon) Pluto's orbit "bounces" on the ecliptic in a strange way.

In the bottm plot, (rotating the coordinates using code "stolen" from ecliptic_latlon), the orbit looks like an ellipse, as I would expect it to be.

Disclaimer: I have zero experience in astronomy, and these are my first attempts with Skyfield, so I it's very likely that i did something wrong. Maybe I messed up the polar to cartesian transformation?

Or maybe ecliptic_latlon is undocumented on purpose?

If the caller types

Topos('0 W', '90 N')

then the result at the moment is the only slightly informative

ValueError: your latitude string '0 W' does not end with either 'N' or 'S'

Instead, we should detect that the W is in fact one of the two possible endings that can be put after a longitude, and recommend that the user offer the attributes in reverse order. Better yet: instead of just saying “put them in reverse order,” we should print a message with the repr's of the two strings that the user has offered, so that (if they want) they can cut and paste from the error message directly back into their code, telling them that they can type

Topos('90 N', '0 W')

I see a few test failures on my Mac:
https://gist.github.com/cdeil/c5f748b62e6653fde78c

Let me know if there's something I can do to help track those down.

Hello,

Thanks you very much for your work!

I am translating pyephem code to skyfield and I have noticed pyephem calculation run at least at 30x faster than skyfield one.

My program wants to calculate all the artificial satellite positions for a given lat/lon/date using the NORAD TLEs

The time consuming math is in:
sat = earth.satellite(tle)
apparent = self.here(self.date).observe(sat)

real 0m31.959s
user 0m31.550s
sys 0m0.208s

compare with the pyephem:
tle.compute(self.here)

real 0m0.766s
user 0m0.532s
sys 0m0.064s

Are you aware of that? ..Or perhaps I am doing in the wrong way ...
For my app speed is critical.. Is there any way to speed up?

Nacho

a vectorized version of ephem is a great idea!

I have this simple az,el to ra,dec conversion, could you please help me to translate it using skyfield?

observer = ephem.Observer()
observer.lon = lon
observer.lat = lat
observer.elevation = alt
observer.date = utc
print observer.radec_of(azimuth, elevation)

I found a typo in http://rhodesmill.org/skyfield/planets.html: "Once an ephemeris file has been downloaded to your current directly,"

from skyfield.api import *
jd = JulianDate(utc=(2014,1,1))
print saturn(jd).position.km

The above yields: [ -1.02986183e+09 -9.95072753e+08 -3.66678097e+08]

Horizons' xyz for the same: -1.029862744489472E+09, -1.058816636837933E+09, 5.939641314965360E+07

Pretty good on the x and y, but way off on the z coordinate. I tried w/ several other times, and got similar results. For reference, the full HORIZONS computational data:

Ephemeris / MAIL_REQUEST Sat Sep 10 10:41:56 2011 Pasadena, USA / Horizons

Target body name: Saturn (699) {source: SAT317}
Center body name: Solar System Barycenter (0) {source: DE405}
Center-site name: BODY CENTER

Start time : A.D. 2014-Jan-01 00:00:00.0000 CT
Stop time : A.D. 2015-Jan-01 00:00:00.0000 CT
Step-size : 6 minutes

Center geodetic : 0.00000000,0.00000000,0.0000000 {E-lon(deg),Lat(deg),Alt(km)}
Center cylindric: 0.00000000,0.00000000,0.0000000 {E-lon(deg),Dxy(km),Dz(km)}
Center radii : (undefined)
Output units : KM-S
Output format : 03
Reference frame : ICRF/J2000.0
Output type : GEOMETRIC cartesian states
Coordinate systm: Ecliptic and Mean Equinox of Reference Epoch

JDCT , , X, Y, Z, VX, VY, VZ, LT, RG, RR,

$$SOE
2456658.500000000, A.D. 2014-Jan-01 00:00:00.0000, -1.029862744489472E+09, -1.058816636837933E+09, 5.939641314965360E+07, 6.398829343803826E+00, -6.760884214350837E+00, -1.365989925103413E-01, 4.930927951351853E+03, 1.478255010756676E+09, 3.791681692850611E-01,
[etc]

This problem may be part of jplephem, though. I did this:

import de421
from jplephem import Ephemeris
eph = Ephemeris(de421)
x, y, z = eph.position('saturn', 2456658.500000000);
print x,y,z

and got the same results: [ -1.02986226e+09] [ -9.95072340e+08] [ -3.66677908e+08]

NOTE: I ran HORIZONS again just now and the Z coordinate is now 5.939672039856304E+07 which is slightly different, but still way off from what skyfield/jplephem gives. HORIZONS now uses SAT360 for Saturn's position, not SAT317 as in my earlier run.

I have no idea what either is, but:

http://naif.jpl.nasa.gov/pub/naif/generic_kernels/spk/satellites/sat360.inp

suggests SAT360 relates to DE430, and:

http://naif.jpl.nasa.gov/pub/naif/generic_kernels/spk/satellites/a_old_versions/sat317.cmt

suggests SAT317 relates to DE421. What's odd is that SAT epidermedes are normally used for satellites, which makes me wonder if I'm confusing two different sats (Saturn and satellites).

repo:
https://github.com/hoylemd/python-skyfield/tree/keplerian_orbits

paper on the task:
http://ssd.jpl.nasa.gov/txt/aprx_pos_planets.pdf

When trying to observe mars multiple times, using Date arrays, I get this TypeError:

>>> from skyfield.api import earth, JulianDate, sun
>>> boston = earth.topos('42.3583 N', '71.0636 W')
>>> jd = JulianDate(utc=(2015, 10, range(1,5)))
>>> boston(jd).observe(sun).apparent().altaz()
(Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/home/dneise/anaconda/lib/python2.7/site-packages/skyfield/units.py", line 165, in __repr__
    return '<{0} {1}>'.format(type(self).__name__, self)
TypeError: __str__ returned non-string (type list)

I was actually trying to calculate the sunset for a certain location on the surface and wanted to interpolate it, by observing the sun using date array. Is there a better way?

Per the discussion brandon-rhodes/pyephem#48 over in the old PyEphem project, Skyfield should add support for planetary reference frames like IAU_JUPITER that express positions relative to the pole and equator of other planets.

ditto: http://stackoverflow.com/questions/36740731/convert-a-direction-vector-to-r-a-and-dec-in-skyfield

Within Skyfield, if I have a vector (x, y, z) from the center of the earth, how can I back-convert it to a point on the "Celestial Sphere" (i.e. Right Ascension and Declination)?

For example purposes, suppose I have a few points in the orbit of a telescope around the earth, and I want to calculate the exact RA and Dec of the direction normal to the orbit - for whatever reason.

Effects like aberration and gravitation can be neglected here - I'm just looking to transform a direction in ICRF to RA and Dec.

In pseudocode:

from skyfield.api import load
import numpy as np

eph = load('de421.bsp')
ts  = load.timescale()
now = ts.now()

vec     = np.array([3141, 2718, 5820], dtype=float)
nvec = vec / np.sqrt((vec**2).sum())  # normalize for the heck of it

earth = eph['earth']
evec  = earth.at(now).vector(vec)   # pseudocode

print "It's pointing toward: ", evec.radec()   # pseudocode

five_seconds = range(58, 58 + 5)
tup = (2015, 6, 30, 23, 59, five_seconds)
jd = JulianDate(utc=tup)

for string in jd.utc_jpl():
print(string)

Working with skyfield I observed a possible memory leak when one repeatedly calculates a satellite's altitude like in the satellite observation notebook.
Feel free to close this issue if I'm doing sth. wrong in the script below.

Steps to reproduce:

1. mkvirtualenv skyfield-test; pip install skyfield
2. Run ./test.py attached below
3. Watch the memory usage grow
4. Hit Ctrl+C and see the output of tracemalloc

The first line of the output (highest memory usage) will be for instance:
[...]/skyfield/timelib.py:602: size=6962 KiB, count=52419, average=136 B
This indicates a leak of def delta_t() in the following code:

def interpolate_delta_t(cache, tt):
    """Given TT, return interpolated Delta T, falling back to a formula."""
    def delta_t():  # TODO
        "Fake placeholder function, until I rewrite how the cache works."
    x, y = cache.run(delta_t)
    delta_t = interp(tt, x, y, nan, nan)
    missing = isnan(delta_t)
    if missing.any():
        tt = tt[missing]
        delta_t[missing] = delta_t_formula_morrison_and_stephenson_2004(tt)
    return delta_t

When replacing the complete function body as shown below, the memory usage is stable.

def interpolate_delta_t(cache, tt):
    return delta_t_formula_morrison_and_stephenson_2004(tt)

test.py script:

#!/usr/bin/env python3

import tracemalloc
from skyfield.api import load, JulianDate, Angle

tle = """ISS (ZARYA)
1 25544U 98067A   16036.52823010  .00009121  00000-0  14318-3 0  9990
2 25544  51.6441 359.5913 0006959  83.4035  25.9071 15.54440498984346
"""
lat = 50.000
lon = 10.000
starttime = (2016, 2, 11, 6, 0)

planets = load("de421.bsp")
earth = planets["earth"]
topos = earth.topos(Angle(degrees=lat), Angle(degrees=lon))
sat = earth.satellite(tle)
t = JulianDate(utc=starttime).tai
t_step = 1.0 / 24 / 60 / 60

tracemalloc.start()
try:
    while True:
        t += t_step
        alt = topos.at(tai=t).observe(sat).altaz()[0].degrees
except KeyboardInterrupt:
    snapshot = tracemalloc.take_snapshot()
    top_stats = snapshot.statistics('lineno')
    for stat in top_stats[:10]:
        print(stat)