P4.6.3The Yale Bright Star catalog
Write a program to create an image of a constellation using the data from the Yale Bright Star Catalog.
Create a class, Star, to represent a star with attributes for its name, magnitude and position in the sky, parsed from the file bsc5.dat which forms part of the catalog. Implement a method for this class which converts the star's position on the celestial sphere as (Right Ascension: \(\alpha\), Declination: \(\delta\)) to a point in a plane, \((x,y)\), for example using the orthographic projection about a central point \((\alpha_0, \delta_0)\):
Suitably-scaled projected, star positions can be output to an SVG image as circles (with a larger radius for brighter stars). For example, the line
<circle cx="200" cy="150" r="5" stroke="none" fill="#ffffff"/>
represents a white circle of radius 5 pixels, centred on the canvas at (200, 150).
Hints: You will need to convert the right ascension from (hr, min, sec) and the declination from (deg, min, sec) to radians - use the data corresponding to "equinox J2000, epoch 2000.0" in each line of bsc5.dat. Let the user select the constellation from the command line using its three-letter abbreviation (e.g. 'Ori' for Orion): this is given as part of the star name in the catalog. Don't forget that star magnitudes are smaller for brighter stars. If you are using the orthographic projection suggested above, choose \((\alpha_0, \delta_0)\) to be the mean of \((\alpha, \delta)\) for stars in the constellation.
Show solution
Here is one solution.
import sys
import math
try:
constellation = sys.argv[1]
except IndexError:
print(
f"usage:\n"
f"python {sys.argv[0]} <constellation>\n"
f"where <constellation> is the three-letter abbreviation for a"
f" constellation name (e.g. Ori, Lyr, UMa, ...)"
)
sys.exit(1)
class Star:
"""
A class representing a star, with name, magnitude, right ascension (ra),
declination (dec) and projected co-ordinates x,y calculated by a class
method.
"""
def __init__(self, name, mag, ra, dec):
"""
Initializes the star object with its name and magnitude, and position
as right ascension (ra) and declination (dec), both in radians.
"""
self.name = name
self.mag = mag
self.ra = ra
self.dec = dec
def project_orthographic(self, ra0, dec0):
"""
Calculates, stores and returns the projected co-ordinates (x, y) of
this star's position using an orthographic projection about the
angular position (ra0, dec0).
"""
delta_ra = self.ra - ra0
self.x = math.cos(self.dec) * math.sin(delta_ra)
self.y = (math.sin(self.dec) * math.cos(dec0) - math.cos(self.dec)
* math.cos(delta_ra) * math.sin(dec0))
return self.x, self.y
stars = []
with open("bsc5.dat", "r") as fi:
for line in fi.readlines():
if line[11:14] != constellation:
# We have no interest of stars which do not belong to constellation
continue
name = line[4:14]
try:
mag = float(line[102:107])
except ValueError:
# some stars do not have magnitudes: ignore these entries
continue
# Right ascension (hrs, mins, secs), equinox J2000, epoch 2000.0
ra_hrs, ra_min, ra_sec = [
float(x) for x in (line[75:77], line[77:79], line[79:83])
]
# Declination (hrs, mins, secs), equinox J2000, epoch 2000.0
dec_deg, dec_min, dec_sec = [
float(x) for x in (line[83:86], line[86:88], line[88:90])
]
# Convert both RA and declination to radians
ra = math.radians((ra_hrs + ra_min / 60 + ra_sec / 3600) * 15.0)
# NB in the Southern Hemisphere be careful to subtract the minutes and
# seconds from the (negative) degrees.
sgn = math.copysign(1, dec_deg)
dec = math.radians(dec_deg + sgn * dec_min / 60 + sgn * dec_sec / 3600)
# Create a new Star object and add it to the list of stars
stars.append(Star(name, mag, ra, dec))
n = len(stars)
if n == 0:
print(f"Constellation {constellation} not found.")
sys.exit(1)
else:
print(f"Found {n} stars in the constellation {constellation}")
# Now calculate the projected co-ordinates of each star, (x,y), finding the
# maximum and minimum values and hence the aspect ratio for our image.
# The "centre" of the constellation in RA, dec:
ra0 = sum([star.ra for star in stars]) / n
dec0 = sum([star.dec for star in stars]) / n
x, y = [None] * n, [None] * n
for i, star in enumerate(stars):
# Orthographic projection (ra, dec) -> (x1, y1)
x[i], y[i] = star.project_orthographic(ra0, dec0)
xmin, xmax, ymin, ymax = min(x), max(x), min(y), max(y)
aspect_ratio = (xmax - xmin) / (ymax - ymin)
# The stars will be output on a canvas of dimensions width x height, with
# some extra padding around the outside of the "paintable" area.
padding = 50
height = 500
width = int(height * aspect_ratio)
# Write the SVG image file for the constellation
with open(f"{constellation}.svg", "w") as f:
print('<?xml version="1.0" encoding="utf-8"?>', file=f)
print('<svg xmlns="http://www.w3.org/2000/svg"', file=f)
print(' xmlns:xlink="http://www.w3.org/1999/xlink"', file=f)
print(
f' width="{width + 2 * padding:d}" height="{height + 2 * padding:d}"'
f' style="background: #000000">',
file=f,
)
for star in stars:
rx = (star.x - xmin) / (xmax - xmin)
ry = (star.y - ymin) / (ymax - ymin)
cx = padding + (1 - rx) * (width - 2 * padding)
cy = padding + (1 - ry) * (height - 2 * padding)
print(
f'<circle cx="{cx:.1f}" cy="{cy:.1f}" r="{max(1, 5 - star.mag):.1f}"'
' stroke="none" fill="#ffffff" name="{star.name}"/>',
file=f,
)
print("</svg>", file=f)
In the above code, we first retrieve the desired constellation abbreviation as a the command line argument; if none is supplied, exit with a usage message.
Next comes the Star class definition, which implements an orthographic projection method, generating and storing coordinates x,y from the star's right ascension and declination.
The star objects themselves are generated by parsing the bsc5.dat file: we convert the angular positions to radians first.
After a quick check that we actually retrieved any stars for the requested constellation, we call the projection method: this returns the coordinates x,y as well as storing them as attributes to each Star object: we need to find the maximum and minimum values of x and y in order to map them to the coordinates of the final image.
To create an SVG image of the constellation we need the coordinates on the SVG "canvas" corresponding to the coordinates of the star on the projected plane. We can easily transform (x, y) to (rx, ry) where rx and ry are between 0 and 1:
rx = (x - xmin)/(xmax - xmin)
ry = (y - ymin)/(ymax - ymin)
Then, if the canvas has width width and height height, we could map (rx, ry) to the canvas coordinates as:
cx = rx * width
cy = ry * height
However, this does not allow for any "padding" around the edges of the constellation image (if the centres of the circles at the extreme edges of the constellation are at the edges of the image, they won't appear as full circles). We can fix this by allowing for an amount of space padding all around the image (so that its total width is width + 2*padding by height + 2*padding), and calculating the canvas position as:
cx = padding + rx * (width - 2*padding)
cy = padding + ry * (height - 2*padding)
There is a further correction to make: because of the nature of the projection, the constellation will look back-to-front with x increasing to the right; this is fixed by replacing rx with 1-rx. Furthermore, because the vertical coordinate of an SVG image increases_downwards_, the constellation will appear upside-down unless we replace ry with 1-ry.
The circles' radii are calculated as max(1, 5 - star.mag) — the brightest star in the night sky, Sirius, has magnitude -1.5 and so is depicted as a circle of radius 6.5; other stars are rendered as smaller circles for increasing star.mag, but we ensure a minimum radius of 1.
An example output image for the constellation Orion is given below.