4
$\begingroup$

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): enter image description here My graph: enter image description here

As @uhoh suggested, I overlaid one graph on top of the other. enter image description here 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?

Thanks in advance.

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:

        # Skip the header
        for i in range(header):
            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()
ax = fig.add_subplot() 

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

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

#multiple lines of data
for star in parseTextFile(file_name, header=1):

    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()
$\endgroup$
5
  • $\begingroup$ 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 $\endgroup$
    – uhoh
    Commented Jan 5, 2020 at 1:25
  • 1
    $\begingroup$ @uhoh Thank you for your reply. I have updated my question which remains unsolved. $\endgroup$
    – Leo Liu
    Commented Jan 5, 2020 at 2:00
  • $\begingroup$ did you convert the units of your plot from parsec to light years before overlaying? $\endgroup$
    – uhoh
    Commented Jan 5, 2020 at 2:20
  • 1
    $\begingroup$ 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! $\endgroup$
    – uhoh
    Commented Jan 5, 2020 at 2:26
  • 2
    $\begingroup$ I'm voting to close this question as off-topic because this is a programming question and does not invoke or require astronomical information. $\endgroup$ Commented Jan 6, 2020 at 18:47

1 Answer 1

2
$\begingroup$

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: enter image description here

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


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:

        # Skip the header
        for i in range(header):
            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()
ax1 = fig.add_subplot(211, projection='3d')
ax2 = fig.add_subplot(212) 

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

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


for star in parseTextFile(file_name, header=1):

    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.add_patch(patches.Rectangle((-60000, 0.03), w_oort,h_oort,color='slateblue',alpha=0.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.

$\endgroup$
3
  • 2
    $\begingroup$ 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. $\endgroup$
    – uhoh
    Commented Jan 5, 2020 at 2:49
  • 1
    $\begingroup$ It's also okay to click "accept" on your own answer $\endgroup$
    – uhoh
    Commented Jan 5, 2020 at 12:58
  • 1
    $\begingroup$ Okay, but the system tells me that I can accept my answer tomorrow. $\endgroup$
    – Leo Liu
    Commented Jan 5, 2020 at 13:14

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .