Support for CIRS coordinates

[jrs65]

Here's a patch for CIRS coordinates adding a cirs option to .radec, and adding the C and CT matrices to Time objects. There's also a set of tests. Some are internal, one is testing against an external calculation using SOFA (this agrees at the 10 micro-arc-second level.

Please let me know if there's any more changes or suggestions that you have.

Fixes #104

[brandon-rhodes]

Thanks for putting this together! I have a few quick comments; once you've responded, I suspect that this can be merged!

[brandon-rhodes]

Since this is a different sort of coordinate, I think it would be clearer to have a separate method cirs_radec() — since so far the pattern with this class is to add a new method for each new coordinate, rather than a single method with many flags.

[jrs65]

Sounds good. Although I have another suggestion I'll put in another comment.

[brandon-rhodes]

Thanks the the comment and for the suggestion that there's a different way forward here (making objects for coordinates, not just methods). But since I suspect that many Skyfield users are well served by the distinction between positions (objects) and coordinates (methods), let's move forward with the plan to make this a separate method. Since some folks might want the xyz coordinates, feel free to break this into two methods cirs_xyz() and cirs_radec() like ecliptic and galactic coordinates do.

I look forward to merging the result — thanks again!

[brandon-rhodes]

I'm a bit surprised than an ephemeris is needed here and that this can't simply be computed using Earth-local vectors — maybe after this lands I'll need to add some print() statements to watch the match work out.

[brandon-rhodes]

Add a newline at the end, please — thanks!

[jrs65]

Your comment about using a .cirs_radec() method got me thinking a little bit more about the semantics of how to do this.

I think an alternative approach which is maybe more consistent with the idea of the equinox-based and CIRS representing frames in their own right would be to add them as specific frame classes. Something like this approximate sketch:

class EquinoxSystem(ICRF):
    ...

class CIRS(ICRF):
    ...

And then add two methods to ICRF that transform the instance into the new frame:

def equinox(self, time):
    # Multiply by time.M to get pos etc...
    return EquinoxSystem(pos)  

def cirs(self, time):
    # Multiply by time.C to get pos etc...
    return CIRS(pos)

and modify the radec method to be like:

def radec(self, epoch=None):

    # For backwards compatibility essentially does what it always did
    if epoch is not None:
        return self.equinox(epoch).radec()

    position_au = self.position.au

    r_au, dec, ra = to_polar(position_au)
    return (Angle(radians=ra, preference='hours'),
                 Angle(radians=dec, signed=True),
                 Distance(r_au))

Overall this would allow making calls like:

# CIRS ...
ra, dec, d = earth.at(st).observe(star).cirs(st).radec()

# Equinox ...
ra, dec, d = earth.at(st).observe(star).equinox(st).radec()
# ... equivalent to the existing ...
ra, dec, d = earth.at(st).observe(star).radec(st)

It's certainly a little bit more complex, but I think it's more conceptually clear, and easy to make backwards compatible, i.e. if a user doesn't care/know about the distinction it just does what they'd always thought it should do.

[brandon-rhodes]

That's a great question, and I should add developer docs addressing this. The thinking behind the design is this:

With previous libraries, I found that new users get confused very quickly about when two different sets of numbers represent truly different positions in the sky (like the astrometric vs the apparent position for a planet; truly two different locations amongst the stars) and when different numbers are simply different ways of naming a single position.

So Skyfield's answer is to never create two different Python objects for the same position. Instead, any given place-in-the-sky gets a single position object and only ever a single position object, with as many methods as we want for turning those into numbers.

So, by this scheme, CORS coordinates must be merely one more of the many methods on the single position object that, whatever the coordinates used to represent it, is exactly and precisely "one place in the sky".

Thanks for asking, it's a very important design principle that I had been needing to write up! Maybe this weekend I can paste this explanation into a file. :)

[brandon-rhodes]

All right! A talk I gave a couple of weeks ago at PyBay gave me practice describing the Skyfield approach that avoids creating separate objects for separate coordinates, which I've now made a part of the official Design document:

b4e3b51

Now that I've better surfaced that rationale, I'm going to officially put this pull request in "Request changes" mode so that it drops temporarily off of my GitHub dashboard. As soon as you're ready to make the quick pivot to a separate method, simply request review and I think we'll be ready to merge.

Thanks for making this new functionality available in Skyfield!

[brandon-rhodes]

Once you've swapped the code over to its own method, simply re-submit and I should be able to merge!

[jrs65]

@brandon-rhodes I've updated this PR with your suggestions.

Unfortunately the build failed for Python 3.6 though I'm not sure that this failure has anything to do with my changes.

[brandon-rhodes]

You're right, the error has nothing to do with your diff! Hopefully I'll have a few minutes to review it in the morning :)

[brandon-rhodes]

Only one real quick comment and I think we'll be done with the review ­— thanks for this update to the PR!

[brandon-rhodes]

What if the opposite is true — what if epoch is indeed None?

[brandon-rhodes]

(Oh — and, I have a vague feeling I should really have named this parameter equinox this whole time. What do you think? Does that sound more correct? Should we go ahead and name it that here?

[jrs65]

Good point, should have thought about this more carefully. I've made it throw the same ValueError if it gets None.

As for the parameter name ,epoch is probably the correct terminology.

[brandon-rhodes]

My worry is this wikipedia article (which could of course be half incorrect):

"Equinox is often confused with epoch, with the difference between the two being that the equinox addresses changes in the coordinate system, while the epoch addresses changes in the position of the celestial body itself"

https://en.wikipedia.org/wiki/Equinox_(celestial_coordinates)

But it's okay for now to simply make this match the rest of the API.

[jrs65]

Hmmmm. I've never come across this distinction before. That nomenclature might make sense for traditional equinox based coordinates, though probably maintaining consistency for the moment is more important.

For CIO based coordinates (e.g. the CIRS) using equinox might be very confusing as it's a coordinate system designed to be separated from the actual equinoxes!

[jrs65]

One thing I've noticed just now: the definitions in to_polar are unconventional and don't match the documentation. I'll create a separate issue for it.

[brandon-rhodes]

Merged! Thanks for this contribution :)