Today, I am illustrating all of the combinations of computing and plotting a constellation with an ancient epoch and/or year using PyEphem. In this notebook, I wanted to double-check my observations in proper-motion.ipynb (nbviewer) with some greater detail. There is a lot more code you'll notice compared to the first example. This is from graphing all of the segments, and that the code accounts for wrapping about 0h-to-24h right ascension. It linearly interpolates for the bound's intersection and splits the segments. There aren't any new conclusions here, but it is definitely easier to see the constellation.
%pylab inline
import ephem
from shapely.geometry import LineString, Point
UMa = {
# star proper name: (abbreviation, hipparcos #),
'dubhe': ('αUMa', 54061),
'merak': ('βUMa', 53910),
'phecda':('γUMa', 58001),
'megrez':('δUMa', 59774),
'alioth':('εUMa', 62956),
'mizar': ('ζUMa', 65378),
'alcor': ('80UMa',65477),
'alcaid':('ηUMa', 67301),
}
def get_hip(name):
for const_name in UMa.keys():
if(name.lower() == const_name.lower()):
return UMa[const_name][1]
constellationship = 'UMa 7 67301 65378 65378 62956 62956 59774 59774 54061 54061 53910 53910 58001 58001 59774'
BLUE = '#268bd2'
BASE03 = '#002b36'
def plot_line(ax, ob):
x, y = ob.xy
ax.plot(x, y, color=BLUE, linewidth=3, solid_capstyle='round', zorder=1)
def plot_pt(ax, ob):
x, y = ob.xy
ax.plot(x, y, 'o',color=BASE03, zorder=1)
def linear_interpolate(x0,y0,x1,y1,x=None,y=None):
# if not passed any x or y, assume want bisection pt
if(x==None and y==None):
x=(x0+x1)/2.
y=(y0+y1)/2.
elif(y==None):
y=y0+(y1-y0)*(x-x0)/(x1-x0)
elif(x==None):
x=x0+(x1-x0)*(y-y0)/(y1-y0)
return x,y
def const(const,years=None,epochs=None,title=None,lines=None,plot_pts=False):
tau = 2.0 * pi
degree = tau / 360.0
hour = tau / 24.0
this_year = ephem.now().tuple()[0]
if(years is None and epochs is None):
years = [this_year]
epochs = [this_year]
elif(years is None):
years = [this_year]
elif(epochs is None):
epochs = [this_year]
n = len(years)*len(epochs)
y = 0
e = 0
fig = pyplot.figure(1,figsize(8*n,5),dpi=180)
for year in years:
y += 1
for epoch in epochs:
e += 1
s = []
for star in const.keys():
s.append(ephem.star(star.capitalize()))
for i in range(len(s)):
if(year is None and epoch is None):
s[i].compute(this_year,epoch=this_year)
elif(epoch is None):
s[i].compute(year,epoch=this_year)
elif(year is None):
s[i].compute(this_year,epoch=epoch)
else:
s[i].compute(year,epoch=epoch)
ra_list = [star.a_ra / hour for star in s]
dec_list = [star.a_dec / degree for star in s]
mag_array = np.array([star.mag for star in s])
size_array = (5 - mag_array) ** 1.5 * 4
ax = fig.add_subplot(1,n,e)
if(title is None):
this_title = 'year={}, epoch={}'.format(year,epoch)
else:
this_title = '{}: year={}, epoch={}'.format(title,year,epoch)
pyplot.title(this_title)
if(lines!=None):
segs = lines.split()
segs.remove(segs[0])
#ns = segs[0]
segs.remove(segs[0])
for i in range(len(segs)):
segs[i] = int(segs[i])
for i in range(0,len(segs),2):
star1 = segs[i]
star2 = segs[i+1]
pt1=pt2=(0.,0.)
for j in range(len(s)):
if(get_hip(s[j].name) == star1):
pt1 = (ra_list[j],dec_list[j])
if(get_hip(s[j].name) == star2):
pt2 = (ra_list[j],dec_list[j])
# gotta check if the lines cross the bounds
p = .25 # NOTE: using 25%
ub_x = 24. # upper bound (RA)
lb_x = 0. # lower bound (RA)
if((pt1[0]>ub_x*(1.-p)) and (pt2[0]<(lb_x+(ub_x*p)))):
if(lines!=None):
#print('crosses higher to lower: ({},{})'.format(pt1,pt2))
# amend left
est_pt1 = linear_interpolate(pt1[0],pt1[1],pt2[0]+ub_x,pt2[1],ub_x)
# amend right
est_pt2 = (lb_x,est_pt1[1])
l1 = LineString([est_pt1,pt1])
plot_line(ax,l1)
l2 = LineString([pt2,est_pt2])
plot_line(ax,l2)
if(plot_pts):
plot_pt(ax,l1)
plot_pt(ax,l2)
axis([0.0,24.0,-90.0,90.0])
elif((pt2[0]>ub_x*(1.-p)) and (pt1[0]<(lb_x+(ub_x*p)))):
if(lines!=None):
#print('crosses lower to higher: ({},{})'.format(pt1,pt2))
# amend left
est_pt1 = linear_interpolate(pt1[0]+ub_x,pt1[1],pt2[0],pt2[1],ub_x)
# amend right
est_pt2 = (lb_x,est_pt1[1])
l1 = LineString([est_pt2,pt1])
plot_line(ax,l1)
l2 = LineString([pt2,est_pt1])
plot_line(ax,l2)
if(plot_pts):
plot_pt(ax,l1)
plot_pt(ax,l2)
axis([0.0,24.0,-90.0,90.0])
else:
l = LineString([pt1,pt2])
if(lines!=None):
plot_line(ax,l)
if(plot_pts):
plot_pt(ax,l)
gca().xaxis.grid(True)
gca().yaxis.grid(True)
gca().invert_xaxis()
return s
Populating the interactive namespace from numpy and matplotlib
The following depictions are made from calling the above function const
for each combination.
This computes the current position of the Big Dipper stars and uses coordinates relative to the direction the earth's axis is pointing currently.
ancient_e = const(UMa,lines=constellationship,plot_pts=True,title='Big Dipper')
#savefig('uma.png')
The following computes the current position of the Big Dipper stars but uses coordinates relative to the direction the earth's axis was pointing at year -100000. This is rather useless unless maybe you are an ancient alien on earth pondering the future and you are planning on getting in your spaceship and staying in the same celestial spot until 2014.
ancient_e = const(UMa,epochs=['-100000'],lines=constellationship,plot_pts=True,title='Big Dipper')
#savefig('uma_epoch.png')
This call to const
computes the positon of the Big Dipper as it was in the year -100000. It returns coordinates for the orientation of the earth that year. At the time of writing, I believe this is the representation I want regarding an ancient ancient observer looking at the night sky.
ancient_ey = const(UMa,years=['-100000'],epochs=['-100000'],lines=constellationship,plot_pts=True,title='Big Dipper')
#savefig('uma_epoch_and_year.png')
This gets the positon of the Big Dipper in the year -100000 but expressed in 2014 coordinates. To me this is useless, but am I wrong? This is what is commonly displayed by many popular stellar programs.
ancient_y = const(UMa,years=['-100000'],lines=constellationship,plot_pts=True,title='Big Dipper')
#savefig('uma_year')
ancient_y = const(UMa,years=['-100000','2014','100000'],lines=constellationship,plot_pts=True,title='Big Dipper')
#savefig('uma_changing_year.png')
ancient_y = const(UMa,epochs=['-100000','2014','100000'],lines=constellationship,plot_pts=True,title='Big Dipper')
#savefig('uma_changing_epoch.png')
ancient_y = const(UMa,years=['-100000','2014','100000'],epochs=['-100000','2014','100000'],lines=constellationship,plot_pts=True,title='Big Dipper')
#savefig('uma_epoch_year_iterations.png')
I believe that the above depiction is the one I want regarding an ancient observer seeing the night sky in his perspective, and a futuristic observer seeing the night sky in his own perspective. This is rather hard to comprehend sometimes, so I would absolutely love some community input on this subject.
We find out later that this is indeed the correct methodology. That is, to vary both epoch and year together. However, the IAU precession model PyEphem inherits from LibAstro is only valid for a few hundred years. That means that although the methodology above is fine, the depictions are broken when using epochs far in the future or far in the past. We're in need of a long-term precession model to handle many thousands of years. In far-flung-proper-motion.ipynb (nbviewer) I find out why this happens, and in vondrak.ipynb (nbviewer) I implement a long-term precession model that works for hundreds of thousands of years.