Arbitrary Position Question

[dreamalligator]

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'>

[brandon-rhodes]

My guess is that what you want to do is to perform your average or combination operation on the raw x-y-z vectors of the stars you are interested in — maybe it would be enough just to add together the vectors to build the Star object that represents your average position? That would be a much more meaningful position, I think, than trying to deal with the polar coordinates RA and declination — because polar coordinates mean something different depending on where on the celestial sphere they are: one degree, say, of RA means “a full degree of the sky” down at the equator but “very nearly nothing” if the declination is pointing up near the pole.

I will look into adding a more basic kind of position than Star that lets you specify the position vector manually instead of having to create the object using RA and declination.

[dreamalligator]

Hi Brandon, thanks for looking into a basic position type!

Yes I totally agree for stars, solely representing a star by RA and declination isn't overly meaningful. I've found it needed lately to express positions in space such as if I were directing a telescope manually, or if I were working with some kind of arbitrary space bound such as the IAU constellation bounds (link, andromeda text file).

[dreamalligator]

Oh! I just thought of something. What do you think of adding an optional epoch field to Star class and the hypothetical general position class? If not specified, then set as None. For example, in the Hipparcos function, since we know that the Hipparcos catalog (I/239) uses epoch J1991.25 (JD2448349.0625 (TT)) all of the positions pulled from the catalog could automatically include the epoch.

If you think this is a good idea, would you choose ICRS, Julian years, or Julian date?

[brandon-rhodes]

The idea of representing boundaries is an interesting one, but one that tends to not involve vectors, because by definition a boundary is a geometric entity on a sphere. If you want to see some recent work I’ve done involving boundaries, see this IPython Notebook:

http://nbviewer.ipython.org/github/brandon-rhodes/astronomy-notebooks/blob/master/Interactive-Sky.ipynb

I could perhaps pull that over into the beginnings of a constellations module for Skyfield. Feel free to open an issue describing what constellation operations you need implemented first so we can discuss them.

The Star class is a vector and all Skyfield vectors use the ICRS, with proper motion having the epoch J2000. Hipparcos objects are adjusted to epoch J2000 on their way in to Skyfield to match this, and I think that this approach with other catalogs will be simpler than switching epochs later.

Meanwhile, I will try to make some progress next week on the idea of a basic position, and update this issue with any progress that I make.

[dreamalligator]

I think the epoch consistency is a good choice that resolves that need for Skyfield. I really like that notebook! If you do happen to implement constellation bounds in Skyfield, see the last couple paragraphs of my post about the IAU bounds. Pierre Barbier found that the IAU listed file bound_18.dat is close to the original but converted, truncated, and vertices were added to aid plotting.

[brandon-rhodes]

All right! I have now had a few minutes to look through your code sample and think about it.

Think about the idea of a visual barycenter: the “middle” of those three stars. The three stars are in a rough circle about halfway up the celestial sphere from the equator toward the north pole: at 3h, 13h, and 18h, they are all directions around the pole, but they are all within about ten degrees of the halfway point of 45°. Their visual barycenter ought therefore to lie somewhere up in the circle they create around the north pole — it should have a fairly high inclination, but probably be tipped down toward Vega because that is the star (at around 18h RA) that is lowest of the three and that will therefore drag the “middle” the farthest down in its own direction.

There are already a number of types in positionlib for representing pure positions. In this case you will want an Astrometric position: the same kind you get when you ask “where is this star tonight?” by asking the Earth to observe it. All you need is to sum the three position vectors of your stars, install that summed vector in a new Astrometric position, and ask it its RA and dec:

from skyfield.api import Star, earth
from skyfield.positionlib import Astrometric

algol = Star(ra_hours=( 3,  8, 10.1315),  dec_degrees=(40, 57, 20.332))
mizar = Star(ra_hours=(13, 23, 55.5),     dec_degrees=(54, 55, 31))
vega  = Star(ra_hours=(18, 36, 56.33635), dec_degrees=(38, 47, 01.2802))

stars = [algol, mizar, vega]
here = earth(utc=(2015, 1, 13))

vector_AU = sum(here.observe(star).position.AU for star in stars)
ra, dec, distance = Astrometric(vector_AU).radec()
print ra
print dec

This gives a quite reasonable answer for a visual barycenter — over towards Vega, and tipped down a bit from a pure 90° solution up at the pole:

18h 54m 54.68s
+78deg 18' 05.5"

I am going to tentatively close this issue, now that I have had time to work on your problem and realize that all you needed was the position classes already in positionlib — but please re-open it if you think there is still more work to do here!