Last active 1 year ago

mga's Avatar mga revised this gist 1 year ago. Go to revision

5 files changed, 184 insertions

README.md(file created)

@@ -0,0 +1,3 @@
1 + ## Dataset
2 +
3 + https://www.kaggle.com/datasets/konivat/hipparcos-star-catalog

distribution_plot.py(file created)

@@ -0,0 +1,50 @@
1 + import os
2 + import pandas as pd
3 + import numpy as np
4 + import matplotlib.pyplot as plt
5 +
6 + # Open the file "./hipparcos-voidmain.csv" for reading
7 + f = open("./hipparcos-voidmain.csv", "r")
8 +
9 + dataset = pd.read_csv(f, sep=",", header=0)
10 +
11 + # Get the data we want to plot
12 + # Apparent magnitude in the V band (visible light)
13 + apparentMag = np.array(dataset["Vmag"], dtype=float)
14 + # B-V color index
15 + bmv = np.array(dataset["B-V"], dtype=float)
16 + # Paralax
17 + paralax = np.array(dataset["Plx"].notna().notnull(), dtype=float)
18 + paralax = paralax[paralax > 0.1]
19 + # Compute distance only if parallax is not zero
20 + distance = 1000 / paralax # unit: parsec
21 + # Absolute magnitude
22 + absoluteMag = apparentMag - 5 * np.log10(distance / 10)
23 +
24 + # Plot the KDE of apparent magnitude, B-V color index, distance and absolute magnitude
25 + figure, axis = plt.subplots(2, 2, figsize=(10, 10))
26 + axis[0, 0].set_title("KDE of apparent magnitude")
27 + axis[0, 0].set_xlabel("Apparent magnitude")
28 + axis[0, 0].set_ylabel("Density")
29 +
30 + axis[0, 1].set_title("KDE of B-V color index")
31 + axis[0, 1].set_xlabel("B-V color index")
32 + axis[0, 1].set_ylabel("Density")
33 +
34 + axis[1, 0].set_title("KDE of distance")
35 + axis[1, 0].set_xlabel("Distance")
36 + axis[1, 0].set_ylabel("Density")
37 +
38 + axis[1, 1].set_title("KDE of absolute magnitude")
39 + axis[1, 1].set_xlabel("Absolute magnitude")
40 + axis[1, 1].set_ylabel("Density")
41 +
42 + axis[0, 0].hist(apparentMag, bins=100, density=True)
43 + axis[0, 1].hist(bmv, bins=100, density=True)
44 + axis[1, 0].hist(paralax, bins=100, density=True)
45 + axis[1, 1].hist(absoluteMag, bins=100, density=True)
46 +
47 + plt.show()
48 +
49 + # Close the file
50 + f.close()

hr_plot.py(file created)

@@ -0,0 +1,47 @@
1 + import os
2 + import pandas as pd
3 + import numpy as np
4 + import matplotlib.pyplot as plt
5 +
6 + import utilities
7 +
8 + # Open the file "/hipparcos-voidmain.csv" for reading
9 + f = open("./hipparcos-voidmain.csv", "r")
10 +
11 + dataset = pd.read_csv(f, sep=",", header=0)
12 +
13 + # Isolate the data we want to plot
14 + inputData = np.array(dataset[["B-V", "e_B-V", "Vmag", "Plx", "e_Plx"]], dtype=float)
15 +
16 + # Filter the data
17 + # Get rid of the NaN values
18 + inputData = inputData[~np.isnan(inputData).any(axis=1)]
19 + # Apply the following conditions:
20 + # - e_B-V should be smaller than 0.1
21 + # - Plx should be larger than 0.1
22 + # - e_Plx should be smaller than 5
23 + inputData = inputData[(inputData[:, 1] < 0.1) & (inputData[:, 3] > 0.1) & (inputData[:, 4] < 5)]
24 +
25 + # Compute the distance
26 + distance = 1000 / inputData[:, 3] # unit: parsec
27 + # Compute the absolute magnitude
28 + absoluteMag = inputData[:, 2] - 5 * np.log10(distance / 10)
29 + bmv = inputData[:, 0]
30 +
31 + # Map B-V to RGB
32 + rgb = np.array([utilities.bv2rgb(item) for item in bmv])
33 +
34 + # Plot the Hertzsprung-Russell diagram
35 + plt.figure(figsize=(10, 10))
36 + plt.scatter(bmv, absoluteMag, s=0.05, c=rgb, alpha=0.75)
37 + plt.xlim(-0.5, 2.5)
38 + plt.ylim(-15, 15)
39 + plt.xlabel("B-V")
40 + plt.ylabel("Vmag")
41 + plt.title("Hertzsprung-Russell diagram")
42 + plt.gca().set_facecolor((0, 0, 0))
43 + plt.gca().invert_yaxis()
44 + plt.show()
45 +
46 + # Close the file
47 + f.close()

position_plot.py(file created)

@@ -0,0 +1,37 @@
1 + import os
2 + import pandas as pd
3 + import numpy as np
4 + import matplotlib.pyplot as plt
5 +
6 + import utilities
7 +
8 + # Open the file "./hipparcos-voidmain.csv" for reading
9 + f = open("./hipparcos-voidmain.csv", "r")
10 +
11 + dataset = pd.read_csv(f, sep=",", header=0)
12 +
13 + # Get the position data
14 + x = np.array(dataset["RAdeg"], dtype=float)
15 + y = np.array(dataset["DEdeg"], dtype=float)
16 +
17 + # Get the magnitude data (will be use to determine the size of the point)
18 + apparentMag = np.array(dataset["Vmag"], dtype=float)
19 + # Convert apparentMag to size
20 + size = 10 ** (apparentMag / -2.5)
21 +
22 + # Get the color data (will be use to determine the color of the point)
23 + bmv = np.array(dataset["B-V"], dtype=float)
24 + # Map B-V to RGB
25 + rgb = np.array([utilities.bv2rgb(bv) for bv in bmv])
26 +
27 + # Plot the stars position distribution
28 + plt.scatter(x, y, s=size, c=rgb, alpha=0.75)
29 + plt.title("Stars position")
30 + plt.xlabel("RAdeg")
31 + plt.ylabel("DEdeg")
32 + plt.gca().set_facecolor('black')
33 + plt.gca().invert_xaxis()
34 + plt.show()
35 +
36 + # Close the file
37 + f.close()

utilities.py(file created)

@@ -0,0 +1,47 @@
1 + # Get the rgb color from B-V color
2 + # More info: https://wikipedia.org/wiki/Color_index
3 + # Function from https://stackoverflow.com/questions/21977786/star-b-v-color-index-to-apparent-rgb-color
4 + # By AymericG
5 + def bv2rgb(bv):
6 + if bv < -0.40: bv = -0.40
7 + if bv > 2.00: bv = 2.00
8 + r = 0.0
9 + g = 0.0
10 + b = 0.0
11 + if -0.40 <= bv<0.00:
12 + t=(bv+0.40)/(0.00+0.40)
13 + r=0.61+(0.11*t)+(0.1*t*t)
14 + elif 0.00 <= bv<0.40:
15 + t=(bv-0.00)/(0.40-0.00)
16 + r=0.83+(0.17*t)
17 + elif 0.40 <= bv<2.10:
18 + t=(bv-0.40)/(2.10-0.40)
19 + r=1.00
20 + if -0.40 <= bv<0.00:
21 + t=(bv+0.40)/(0.00+0.40)
22 + g=0.70+(0.07*t)+(0.1*t*t)
23 + elif 0.00 <= bv<0.40:
24 + t=(bv-0.00)/(0.40-0.00)
25 + g=0.87+(0.11*t)
26 + elif 0.40 <= bv<1.60:
27 + t=(bv-0.40)/(1.60-0.40)
28 + g=0.98-(0.16*t)
29 + elif 1.60 <= bv<2.00:
30 + t=(bv-1.60)/(2.00-1.60)
31 + g=0.82-(0.5*t*t)
32 + if -0.40 <= bv<0.40:
33 + t=(bv+0.40)/(0.40+0.40)
34 + b=1.00
35 + elif 0.40 <= bv<1.50:
36 + t=(bv-0.40)/(1.50-0.40)
37 + b=1.00-(0.47*t)+(0.1*t*t)
38 + elif 1.50 <= bv<1.94:
39 + t=(bv-1.50)/(1.94-1.50)
40 + b=0.63-(0.6*t*t)
41 + return (r, g, b)
42 +
43 + # Get the temperature from B-V color
44 + # More info: https://wikipedia.org/wiki/Color_index
45 + def bv2temp(bv):
46 + temp=4600.0*(1.0/(0.92*bv+1.7)+1.0/(0.92*bv+0.62))
47 + return temp
Newer Older