# Plotting Astronomical Data Using Python

I have been trying to create a graph which shows the distances between the stars and the earth, yet I could not obtain the desired graph as seen below (click to zoom in): My graph:

As @uhoh suggested, I overlaid one graph on top of the other. I apologize for posting this confusing graph

The combined graph denotes that the positions of the curves in my graph totally differ from those in the expected graph.

I followed this guide to calculate the position functions of the stars.

My Question: Did I made any mistakes that resulted in the anomalies?

Annotations for the Code:

• Name - Name of the star

• RA - Right Ascension in degrees, ICRS coord. (J2000)

• Dec - Declination in degrees, ICRS coord. (J2000)

• pm_mRA - Proper motion in right ascension, in miliarcseconds per year

• pm_mdec - Proper motion in declination, in miliarcseconds per year

• vr - radial velocity in kilometers per second, a positive value means that the star is moving away from us

• mparallax - parallax of the star in miliarcseconds

• d - distance between the star and the earth

My Code:

def parseTextFile(file_name, delimiter=",", header=0):
""" Parse a text file to a list. The file contents are delimited and have a header. """

with open(file_name) as f:

next(f)

data = []

# Parse file contents
for line in f:

# Remove the newline char
line = line.replace('\n', '').replace('\r', '')

# Split the line by the delimiter
line = line.split(delimiter)

# Strip whitespaces from individual entries in the line
for i, entry in enumerate(line):
line[i] = entry.strip()

# Add the contents of the line to the data list
data.append(line)

return data

fig = plt.figure()

#time span
time = np.arange(-60000,100000,10)
count = 1

xdic = {}
ydic = {}
zdic = {}

#multiple lines of data

name = str(star[0])
RA = float(star[1])
Dec = float(star[2])
pm_mRA = float(star[3])
pm_mDec = float(star[4])
vr = float(star[5])
mparallax = float(star[6])

pm_RA = pm_mRA * 0.001
pm_Dec = pm_mDec * 0.001
d = 1 / (mparallax * 0.001)

#Transverse velocities
vta = pm_RA * d * 4.740
vtd = pm_Dec * d * 4.740

#Linear velocities
vx = vr * np.cos(Dec) * np.cos(RA) - vta * np.sin(RA) - vtd * np.sin(Dec) * np.cos(RA)
vy = vr * np.cos(Dec) * np.sin(RA) + vta * np.cos(RA) - vtd * np.sin(Dec) * np.sin(RA)
vz = vr * np.sin(Dec) + vtd * np.cos(Dec)

#unit conversion from km/s to pc/year
vx_pcyr = vx / 977780
vy_pcyr = vy / 977780
vz_pcyr = vz / 977780

#initial positions
xi = d * np.cos(Dec) * np.cos(RA)
yi = d * np.cos(Dec) * np.sin(RA)
zi = d * np.sin(Dec)

#position functions
x = xi + vx_pcyr * time
y = yi + vy_pcyr * time
z = zi + vz_pcyr * time

distance = np.sqrt(x ** 2 + y ** 2 + z ** 2)

ax.plot(time,distance,label=name)

ax.set_xlabel('Time (Year)')
ax.set_ylabel('Distance (pc)')
ax.legend()
plt.show()

• Your question is well written, but I wonder what kind of answer you are expecting? What makes you think your plot might be incorrect? Try this: plot using the same units and scale as the original (0.0 to 10.0 light years) using the correct conversion factor and ax.set_ylim(0.0, 10.0) then save a copy with a transparent background with plt.savefig('mytransparentplot', transparent=True) and overlay it on the original (in powerpoint or any image manipulating tool) and stretch it and see if it overlays perfectly. Or print out some distances numerically and ask how to check those
– uhoh
Jan 5 '20 at 1:25
• @uhoh Thank you for your reply. I have updated my question which remains unsolved. Jan 5 '20 at 2:00
• did you convert the units of your plot from parsec to light years before overlaying?
– uhoh
Jan 5 '20 at 2:20
• If you solve the problem it's perfectly fine to post a short answer and to click "accept". No need to close it, the problem has been solved. Welcome to Stack Exchange!
– uhoh
Jan 5 '20 at 2:26
• I'm voting to close this question as off-topic because this is a programming question and does not invoke or require astronomical information. Jan 6 '20 at 18:47

I forgot changing the unit of the distance from parsec to light year. This is a simple unit conversion error which should have been avoided.

Ultimate Graph:

Data Used:

Name,RA(deg),Dec(deg),pm_RA(mas/yr),pm_Dec(mas/yr),Vr(km/s),parallax(mas)
Wolf 359,164.120271,+07.014658,-3842.0,-2725.0,19.321,418.3
Proxima Centauri,217.42895219,-62.67948975,-3775.75,765.54,-22.40,768.13
Alpha Centauri,219.900850,-60.835619,-3608,686,-22.3,742
Barnard's star,269.45207511,+04.69339088,-798.58,10328.12,-110.51,548.31
Gliese 445,176.92240640,+78.69116300,743.61,481.40,-111.65,186.86
Luhman 16A,150.8218675,-53.319405556,-2754.77,358.72,23.1,500.51
Sirius,101.28715533,-16.71611586,-546.01,-1223.07,-5.50,379.21
Lalande 21185,165.83414166,+35.96988004,-580.27,-4765.85,-84.69,392.64
Ross 248,355.479122,+44.177994,115.10,-1592.77,-77.715,316.7
Gliese 65,24.756054,-17.950569,3321,562,29,373.70

Code:

%matplotlib qt
import numpy as np
from mpl_toolkits.mplot3d import axes3d
from matplotlib import pyplot as plt
from matplotlib import patches

""" Parse a text file to a list. The file contents are delimited and have a header. """

with open(file_name) as f:

next(f)

data = []

# Parse file contents
for line in f:

# Remove the newline char
line = line.replace('\n', '').replace('\r', '')

# Split the line by the delimiter
line = line.split(delimiter)

# Strip whitespaces from individual entries in the line
for i, entry in enumerate(line):
line[i] = entry.strip()

# Add the contents of the line to the data list
data.append(line)

return data

if __name__ == "__main__":

file_name = 'C:\\Users\\The Wings of Dream\\Desktop\\UWO-PA-Python-Course\\Lecture 5\\Homework 2\\star_data.txt'

#Program Begin:

fig = plt.figure()

time = np.arange(-60000,100000,10)
count = 1

xdic = {}
ydic = {}
zdic = {}

name = str(star[0])
RA = float(star[1])
Dec = float(star[2])
pm_mRA = float(star[3])
pm_mDec = float(star[4])
vr = float(star[5])
mparallax = float(star[6])

pm_RA = pm_mRA * 0.001
pm_Dec = pm_mDec * 0.001
d = 1 / (mparallax * 0.001)

vta = pm_RA * d * 4.740
vtd = pm_Dec * d * 4.740

vx = vr * np.cos(Dec) * np.cos(RA) - vta * np.sin(RA) - vtd * np.sin(Dec) * np.cos(RA)
vy = vr * np.cos(Dec) * np.sin(RA) + vta * np.cos(RA) - vtd * np.sin(Dec) * np.sin(RA)
vz = vr * np.sin(Dec) + vtd * np.cos(Dec)

vx_pcyr = vx / 977780
vy_pcyr = vy / 977780
vz_pcyr = vz / 977780

xi = d * np.cos(Dec) * np.cos(RA)
yi = d * np.cos(Dec) * np.sin(RA)
zi = d * np.sin(Dec)

x = xi + vx_pcyr * time
y = yi + vy_pcyr * time
z = zi + vz_pcyr * time

xdic['x'+str(count)] = x
ydic['y'+str(count)] = y
zdic['z'+str(count)] = z

distance = np.sqrt(x ** 2 + y ** 2 + z ** 2) * 3.26156

ax1.plot(xdic['x'+str(count)],ydic['y'+str(count)],zdic['z'+str(count)])
ax2.plot(time,distance,label=name)

count = count + 1

w_oort, h_oort = 160000, 3.2
ax2.annotate('Oort Cloud', xy=(15000,1.6), size=12)

plt.axvline(x=0,color='gray',linestyle='--',linewidth='0.5')

#plotting constraints
ax2.set_ylim(0.0, 10.0)
ax2.set_xlim(-60000, 100000)

ax1.set_xlabel('x axis')
ax1.set_ylabel('y axis')
ax1.set_zlabel('z axis')
ax1.title.set_text('Motion of Stars in Space')
ax2.title.set_text('Distance-Time')
ax2.set_xlabel('Time (Year)')
ax2.set_ylabel('Distance (Light Year)')
ax2.legend()
plt.show()


ax1 subplot will give you a 3d parametricplot showing the motion of the stars.

• Happy ending to an interesting story! I think the plot and the original data is really interesting, and adding and accepting an answer rather than closing it will probably allow more people to see it.
– uhoh
Jan 5 '20 at 2:49