Workspace
Ricardo Simão/

Star type, Temperature, Luminosity, a galactic 3D-map and more

0
Beta
Spinner
# Start coding here... 
Use the hipparcos mission data to plot the Hertzsprung Russell

diagram and from it estimate the mass, radius, luminosity, effective temperature. Also make a 3D-map of the bodies covered in this catalogue, map the density distribution, map the velocity profile, map the energy density, estimate the distribution of star-type per quadrant.

I am a novice and this is my first publicly available

data science project. As the project is big, I decided to split it into at least two parts, this being the first one.

Here I shall download the .csv.gz file from the ESA repository

untar it in a .csv file, get the relevant parameters, calculate the interest variables, estimate new variables from this calculated variables, plot the 3D-map given two distance limits, make some plots to understand the data

# libraries
import gzip
import urllib.request
import re
import pandas as pd
import numpy as np
import seaborn as sns
import matplotlib as mpl
import matplotlib.pyplot as plt
from matplotlib.ticker import AutoMinorLocator
import math
from scipy.optimize import curve_fit
import scipy.stats as stats
#----------------------------------------------------------------
### useful functions 
def poly(x, vec):
    '''Return polinomial function
    ______________________
    Takes: the number which the function is calculated(x: a float) and
    the vector of coefficients (vec: an array)
    ______________________
    Returns: the value of the polynimial at the point x (a float) '''
    temp = 0
    for coef in range(len(vec)):
        temp = temp + vec[coef]*x**coef
    return temp     
def point_distance(p_1, p_2):
    '''Calculate the distance between two points in a 2d environment
    ______________________
    Takes: two tuples containig (x, y) locations of each
    ______________________
    Returns: the distance between these points(a float) '''
    x_1, y_1 = p_1
    x_2, y_2 = p_2
    return np.sqrt((x_1 - x_2)**2 + (y_1 + y_2)**2)
def weigths(x , y, limits): 
    '''Calculate the weigths (weigthed probability) of the star types as the inverse 
    of the distanceto the star type lines in the pre-defined limits
    ______________________
    Takes: the coordinates of the point and the star type lines limits 
    ______________________
    Returns: a normalized vector of the distance of the point to each star type
    lines (a float), which may be interpreted as probabilities '''
    distances = []
    for limit in limits:
        x_il, x_sl, coef = limit
        if x >= x_il and x <= x_sl : 
            distances.append(1/max(abs(y - poly(x, vec=coef)), 0.01))
        elif x < x_il:
            distances.append(1/max(point_distance((x, y), (x_il, poly(x_il, vec=coef))), 0.01))
        elif x > x_sl:
            distances.append(1/max(point_distance((x, y), (x_sl, poly(x_sl, vec=coef))), 0.01))
    return distances/np.sum(distances)
#----------------------------------------------------------------------

Libraries and some useful functions loaded, let us write some functions to download and process the data

####### functions to import the data from the ESA website and preprocessing ########
def web_download(url_file, filename):
    """Download the file named filename
       ____________
       Args (str): the url addres of the file       
       ____________
       Returns (str): name pf the dowloaded data 
    """
    try:
        urllib.request.urlretrieve(url_file, filename)
        return(filename)
    except:
        print('NonExistentFile')
        return False 
def gz_data_extraction(url, output):
    """Extract a .gz data from the web apri uncompress it
       ____________
       Args (str): its url adress 
       ____________
       Returns (str): None, bu it saves the file in another file
    """
    # Download archive
    try:
        # Read the file inside the .gz archive located at url
        with urllib.request.urlopen(url) as response:
            with gzip.GzipFile(fileobj=response) as uncompressed:
                file_content = uncompressed.read()
        # write to file in binary mode 'wb'
        with open(output, 'wb') as f:
            f.write(file_content)
            return 0
    except Exception as e:
        print(e)
        return None
def dat_to_dataframe(filename):
    """Eliminate the empty spaces of each line and import the data as a dataframe
       ____________
       Args (str): the name of the file  
       ____________
       Returns (str): returns a dataframe containing the data in the file
    """
    file = open(filename, "r")
    file_lines = file.readlines()
    new = []
    for i in range(len(file_lines)):
        temp_line = re.sub("\s+", ",", file_lines[i].strip())
        values = temp_line.split(',')
        #according to the documentation, the lines beyond the 27-th are entries of a matrix, which we will not use
        new.append(values[0:27])
    file.close()
    # data taken from the Read me file 
    names = ['HIP', 'Sn', 'So', 'Nc', 'RArad', 'DErad', 'Plx', 'pmRA', 'pmDE', 'e_RArad', \
    'e_DErad', 'e_Plx', 'e_pmRA', 'e_pmDE', 'Ntr', 'F2', 'F1', 'var', 'ic', 'Vmag', \
    'e_Hpmag', 'sHp', 'VA', 'B-V', 'e_B-V', 'V-I', 'UW']
    return  pd.DataFrame(new, columns = names)
def preprocessing():
    """bundle togheter the preprocessing steps
       ____________
       Args (str): None
       ____________
       Returns (str): a dataframe of processed data 
    """
    # main file download
    gz_data_extraction('https://cdsarc.cds.unistra.fr/ftp/I/311/hip2.dat.gz', 'hip2.dat')
    print('gz downloaded from the web')
    # read me download 
    web_download('https://cdsarc.cds.unistra.fr/ftp/I/311/ReadMe', "ReadMe_hipparcos_2")
    print('Read me downloaded from the web')
    df = dat_to_dataframe("./hip2.dat")
    print('data transformed to dataframe')
    df = df[['B-V','Vmag', 'Plx', 'RArad', 'DErad']] 
    df = df.astype(float)
    print('B-V','Vmag', 'Plx', 'RArad', 'DErad features selected')
    temp_size = df.shape[0]
    df = df[abs(df['Plx']) >= 0.0001] 
    per_of_del_data = (1 - df.shape[0]/temp_size)*100
    print("Far away objects were deleted (Plx approx 0). A total of {var}% of the data was lost due this procedure.".format(var=per_of_del_data))
    df = df.reset_index()
    return df
#-----------------------------------------------------------------------
To check the data processing, let me define some easy to use data visualization functions to check the sanity of the data
#-----------------------------------------------------------------------
### plot histograms 
def plot_histograms(df, variables, n_rows, n_cols) : 
    """Plot histograms of the data
       ____________
       Args (str): the dataframe (df), the variables to plot (variables) the number of rows
       and columns to plot (n_rows, n_columns)
       ____________
       Returns (str): None
    """
    fig=plt.figure()
    for i, var_name in enumerate(variables):
        ax=fig.add_subplot(n_rows,n_cols,i+1)
        df[var_name].hist(bins=80,ax=ax, color='black')
        ax.set_title(var_name+" Distribution")
    fig.tight_layout()  # Improves appearance a bit.
    plt.show()    
    return None
def print_NaNs(df) : 
    """Plot histograms of the relative number of NaNs. The plot presents a horizintal line 
    at the heigth 5% so that it is ok the delete the data if the bar is below this line
       ____________
       Args: the dataframe (df)
       ____________
       Returns : None
    """
    NaSum = df.isna().mean().values  
    p = sns.barplot(y=NaSum, x=df.columns, color='black')
    p.axhline(y=0.05, linestyle = '--', color = 'red' )
    plt.xticks(rotation=90)
    plt.show()
    return None
#-----------------------------------------------------------------------
Now, functions to calculate the features that can be taken from the data alone. Also is present a function to plot the HR diagram, a function to interpolate the typical region of each star type and a function to get the coefficients of a polynomial interpolation to those regions.
#-----------------------------------------------------------------------
### Add calculable properties
def add_distance(df) : 
    """Calculate the distances from the object using the formula 
    D = 1/plx * 1000 * 3,26156 (in ly) and include 
    in the dataframe
       ____________
       Args: the dataframe (df)
       ____________
       Returns: the dataframe with the additinal column
    """
    df['Distance'] = abs(3.26156*1000/np.array(df["Plx"], dtype=float))
    return df 
def add_cartesian_coordinates(df) :
    """Calculate the cartesian coordinates of the objects for posterior plotting
       ____________
       Args: the dataframe (df)
       ____________
       Returns: the dataframe with three additinal columns (the x,y,z coordinates)
    """
    r = np.array(df['Distance'])
    cos_theta =  np.array(np.cos(df['DErad']))
    sin_theta =  np.array(np.sin(df['DErad']))
    cos_phi = np.array(np.cos((df['RArad'])))
    sin_phi = np.array(np.sin(df['RArad']))
    df['x_pos'] = r * cos_theta * sin_phi
    df['y_pos'] = r * sin_theta * sin_phi
    df['z_pos'] = r * cos_phi
    return df
def add_abs_magnitude(df) :
    """Calculate the absolute magnitude using the formula  M = m + 5 - 5LOG_10(d(pc))
    notice the distance in parsecs
       ____________
       Args: the dataframe (df)
       ____________
       Returns: the dataframe with the additinal column, the absolute magnitude
    """
    df['AbsMag'] = np.array(df["Vmag"]) + 2.43287314 -5*np.log10(np.array(df["Distance"])/3.26156)
    return df
#-------------------------------------------------------------------------
### Plot the HR-Diagram
def HR_diagram(Magnitude_B_V, Magnitude_absoluta) :
    """Plot the HR-diagram 
       ____________
       Args: the B_V magnitude feature and the absolute magnitude feature
       ____________
       Returns: None
    """
    data = pd.DataFrame({'X': Magnitude_B_V, 'Y': Magnitude_absoluta})
    F1 = {"family": "serif", "weight": "bold", "color": "black", "size": 18}
    F2 = {"family": "serif", "weight": "bold", "color": "black", "size": 20}
    norm = mpl.colors.Normalize(vmin=Magnitude_B_V.min(), vmax=Magnitude_B_V.max())
    fig, ax = plt.subplots(figsize=(10, 8))
    ax.scatter(data['X'], data['Y'], c=data['X'],  cmap = "YlOrRd",\
        marker = '.', s = 2, linewidth = 0.0001, norm=norm)
    plt.xlabel("Color Index (B-V)", fontdict = F1)
    plt.ylabel("Absolute Magnitude (Mv)", fontdict = F1)
    plt.title("Hertzsprung-Russell diagram", fontdict = F2)
    plt.xlim(Magnitude_B_V.min(), Magnitude_B_V.max())
    plt.ylim(Magnitude_absoluta.max(), Magnitude_absoluta.min()) 
    ax.xaxis.set_minor_locator(AutoMinorLocator())
    ax.yaxis.set_minor_locator(AutoMinorLocator())
    Color_BC = plt.gca()
    Color_BC.set_facecolor("darkblue")
    plt.tight_layout()
    plt.show()
def HR_diagram_regions(sc_x, sc_y) :
    """Plot the HR-diagram, plot the limits of the regions considered to each star type 
    in the HR diagram, plot a polynomial regression  of each subset considering the limits
    and return those regressions cefficients
       ____________
       Args: The B_V magnitude and the Absolute magnitude (the same used to plot the HR diagram)
       ____________
       Returns: 7 arrays, each containing the parameters of the polynomial coefficients to fit
    the HR diagram and define the Star type regions
    """
    ### plot the points
    data = pd.DataFrame({'X': sc_x, 'Y': sc_y})
    F1 = {"family": "serif", "weight": "bold", "color": "black", "size": 18}
    F2 = {"family": "serif", "weight": "bold", "color": "black", "size": 20}
    norm = mpl.colors.Normalize(vmin=sc_x.min(), vmax=sc_x.max())
    fig, ax = plt.subplots(figsize=(10, 8))
    ax.scatter(data['X'], data['Y'], c=data['X'],  cmap = "YlOrRd",\
        marker = '.', s = 2, linewidth = 0.0001, norm=norm)
    plt.xlabel("Color Index (B-V)", fontdict = F1)
    plt.ylabel("Absolute Magnitude (Mv)", fontdict = F1)
    plt.title("Hertzsprung-Russell diagram", fontdict = F2)
    plt.xlim(sc_x.min(), sc_x.max())
    plt.ylim(sc_y.max(), sc_y.min()) 
    ax.xaxis.set_minor_locator(AutoMinorLocator())
    ax.yaxis.set_minor_locator(AutoMinorLocator())
    Color_BC = plt.gca()
    Color_BC.set_facecolor("darkblue")
    plt.tight_layout()
    # useful polynomials 
    def func_2(x, a_0, a_1, a_2):
        return a_0 + a_1*x + a_2*x**2
    def func_3(x, a_0, a_1, a_2, a_3):
        return a_0 + a_1*x + a_2*x**2 + a_3*x**3
    def func_4(x, a_0, a_1, a_2, a_3, a_4):
        return a_0 + a_1*x + a_2*x**2 + a_3*x**3 + a_4*x**4
    def func_10(x, a_0, a_1, a_2, a_3, a_4, a_5, a_6, a_7, a_8, a_9, a_10):
        return a_0 + a_1*x + a_2*x**2 + a_3*x**3 + a_4*x**4 + a_5*x**5 + a_6*x**6 +\
            a_7*x**7 + a_8*x**8+ a_9*x**9 + a_10*x**10
    ### fitting degree 2 polynomial on the white dwarfs region
    f_test = lambda x: 5*x + 6
    test = pd.DataFrame({'X' : np.linspace(-0.4, 0.8, num = 50)})
    test['Y'] = test['X'].map(f_test)
    #sns.lineplot(x=test['X'], y=test['Y'],color='red', ax=ax)  
    new_df=pd.DataFrame({'x_0': sc_x, 'y_0':sc_y})
    new_df['limit'] = new_df['x_0'].map(f_test)
    new_df = new_df[new_df['y_0']>new_df['limit']]
    #sns.scatterplot(x=new_df['x_0'], y=new_df['y_0'], marker='*', color='gold')
    popt, pcov = curve_fit(func_2, new_df['x_0'], new_df['y_0'])
    test['Z'] = test['X'].map(lambda x: popt[0] + popt[1]*x + popt[2]*x**2)
    sns.lineplot(x=test['X'], y=test['Z'],color='white', ax=ax)  
    popt_wd  = popt.copy()
    ### fitting degree 4 polynomial on the brown dwarfs region
    test = pd.DataFrame({'X' : np.linspace(1.35, 2.5, num = 50)})
    test['Y_0'] = test['X'].map(lambda x: 15*x - 16.5)
    test['Y_1'] = test['X'].map(lambda x: 4.5*x - 4)
    #sns.lineplot(x=test['X'], y=test['Y_0'],color='red', ax=ax)  
    #sns.lineplot(x=test['X'], y=test['Y_1'],color='red', ax=ax)  
    new_df=pd.DataFrame({'x_0': sc_x, 'y_0':sc_y})
    new_df['limit_0'] = new_df['x_0'].map(lambda x: 15*x - 16.5)
    new_df['limit_1'] = new_df['x_0'].map(lambda x: 4.5*x - 4)
    new_df = new_df[new_df['y_0']<new_df['limit_0']]
    new_df = new_df[new_df['y_0']>new_df['limit_1']]
    new_df = new_df[new_df['x_0']>1.35]
    #sns.scatterplot(x=new_df['x_0'], y=new_df['y_0'], marker='*', color='gold')
    #sns.regplot(x=new_df['x_0'], y=new_df['y_0'], order = 4, ax= ax)
    popt, pcov = curve_fit(func_4, new_df['x_0'], new_df['y_0'])
    test['Z'] = test['X'].map(lambda x:  popt[4]*x**4 + popt[3]*x**3 + popt[2]*x**2 + popt[1]*x + popt[0])
    test = test[(test['X'] > 1.35)]
    test = test[(test['X'] < 2.5)]
    sns.lineplot(x=test['X'], y=test['Z'],color='white', ax=ax)  
    popt_bd  = popt.copy()
    # fitting the subdwarf branch with 3th degree polynomial
    test = pd.DataFrame({'X' : np.linspace(-0.4, 3.5, num = 50)})
    test['Y_0'] = test['X'].map(lambda x: 5*x + 6)
    test['Y_1'] = test['X'].map(lambda x: 6.5*x )
    #sns.lineplot(x=test['X'], y=test['Y_0'],color='red', ax=ax) 
    #sns.lineplot(x=test['X'], y=test['Y_1'],color='red', ax=ax) 
    new_df=pd.DataFrame({'x_0': sc_x, 'y_0':sc_y})
    new_df['limit_0'] = new_df['x_0'].map(lambda x: 5*x + 6)
    new_df['limit_1'] = new_df['x_0'].map(lambda x: 6.5*x)
    new_df = new_df[new_df['y_0']<new_df['limit_0']]
    new_df = new_df[new_df['y_0']>new_df['limit_1']]
    #sns.scatterplot(x=new_df['x_0'], y=new_df['y_0'], marker='*', color='gold')
    #sns.regplot(x=new_df['x_0'], y=new_df['y_0'], order = 3, ax= ax)
    popt, pcov = curve_fit(func_3, new_df['x_0'], new_df['y_0'])
    test['Z'] = test['X'].map(lambda x:  popt[3]*x**3 + popt[2]*x**2 + popt[1]*x + popt[0])
    test = test[(test['X'] < 1.4)]
    sns.lineplot(x=test['X'], y=test['Z'],color='white', ax=ax)  
    popt_sd  = popt.copy()
    # fitting the Main sequence stars to a 10-th order polynomial
    test = pd.DataFrame({'X' : np.linspace(-0.4, 3.5, num = 50)})
    test['Y_0'] = test['X'].map(lambda x: 6.5*x )
    test['Y_1'] = test['X'].map(lambda x: 6.5*x - 4 )
    #sns.lineplot(x=test['X'], y=test['Y_0'],color='red', ax=ax)  
    #sns.lineplot(x=test['X'], y=test['Y_1'],color='red', ax=ax) 
    new_df=pd.DataFrame({'x_0': sc_x, 'y_0':sc_y})
    new_df['limit_0'] = new_df['x_0'].map(lambda x: 6.5*x)
    new_df['limit_1'] = new_df['x_0'].map(lambda x: 6.5*x - 4)
    new_df = new_df[new_df['y_0']<new_df['limit_0']]
    new_df = new_df[new_df['y_0']>new_df['limit_1']]
    new_df = new_df[new_df['x_0']>-0.3]
    #sns.scatterplot(x=new_df['x_0'], y=new_df['y_0'], marker='*', color='gold')
    #sns.regplot(x=new_df['x_0'], y=new_df['y_0'], line_kws={'color':'black'}, scatter_kws={'s':1}, order=10)
    popt, pcov = curve_fit(func_10, new_df['x_0'], new_df['y_0'])
    test['Z'] = test['X'].map(lambda x: popt[10]*x**10 + popt[9]*x**9 + popt[8]*x**8 + popt[7]*x**7 + \
        popt[6]*x**6 + popt[5]*x**5 + popt[4]*x**4 + popt[3]*x**3 + \
        popt[2]*x**2 + popt[1]*x + popt[0])
    test = test[test['X']>-0.3]
    test = test[test['X']<1.8]
    sns.lineplot(x=test['X'], y=test['Z'],color='white', ax=ax)     
    popt_ps  = popt.copy()
    # fitting the subgiant branch with 2th degree polynomial
    test = pd.DataFrame({'X' : np.linspace(0.36, 2, num = 50)})
    test['Y_0'] = test['X'].map(lambda x: -5*(x-1.3)**2 + 3)
    test['Y_1'] = test['X'].map(lambda x: x - 2 )
    #sns.lineplot(x=test['X'], y=test['Y_0'],color='red', ax=ax) 
    #sns.lineplot(x=test['X'], y=test['Y_1'],color='red', ax=ax)
    new_df=pd.DataFrame({'x_0': sc_x, 'y_0':sc_y})
    new_df['limit_0'] = new_df['x_0'].map(lambda x: -5*(x-1.3)**2 + 3)
    new_df['limit_1'] = new_df['x_0'].map(lambda x: x - 2 )
    new_df = new_df[new_df['y_0']<new_df['limit_0']]
    new_df = new_df[new_df['y_0']>new_df['limit_1']]
    #sns.scatterplot(x=new_df['x_0'], y=new_df['y_0'], marker='*', color='gold')
    #sns.regplot(x=new_df['x_0'], y=new_df['y_0'], line_kws={'color':'black'}, scatter_kws={'s':1}, order=2)
    popt, pcov = curve_fit(func_2, new_df['x_0'], new_df['y_0'])
    test['Z'] = test['X'].map(lambda x: popt[2]*x**2 + popt[1]*x + popt[0])
    sns.lineplot(x=test['X'], y=test['Z'],color='white', ax=ax)
    popt_sG  = popt.copy()
    # fitting the giant branch with 4th degree polynomial
    test = pd.DataFrame({'X' : np.linspace(-0.4, 3.2, num = 50)})
    test['Y_0'] = test['X'].map(lambda x: 0.5*x - 6)
    test['Y_1'] = test['X'].map(lambda x: 0.5*x - 3)
    new_df=pd.DataFrame({'x_0': sc_x, 'y_0':sc_y})
    new_df['limit_0'] = new_df['x_0'].map(lambda x: 0.5*x - 6)
    new_df['limit_1'] = new_df['x_0'].map(lambda x: 0.5*x - 3)
    new_df = new_df[new_df['y_0']>new_df['limit_0']]
    new_df = new_df[new_df['y_0']<new_df['limit_1']]
    #sns.scatterplot(x=new_df['x_0'], y=new_df['y_0'], c='red', marker = '*', ax=ax) 
    #sns.lineplot(x=test['X'], y=test['Y_0'],color='red', ax=ax) 
    #sns.lineplot(x=test['X'], y=test['Y_1'],color='red', ax=ax) 
    #sns.regplot(x=new_df['x_0'], y=new_df['y_0'], order = 4)
    #sns.lineplot(x=test['X'], y=test['Y_2'],color='red', ax=ax)
    popt, pcov = curve_fit(func_4, new_df['x_0'], new_df['y_0'])
    test['Z'] = test['X'].map(lambda x: popt[4]*x**4 + popt[3]*x**3 + popt[2]*x**2 + popt[1]*x + popt[0]  )
    test = test[(test['X'] > 0.)]
    #test = test[(test['X'] < 1.9)]
    sns.lineplot(x=test['X'], y=test['Z'],color='white', ax=ax)
    popt_G  = popt.copy()
    # fitting a 2-th polynomial to the Bright Giant Ib branch 
    test = pd.DataFrame({'X' : np.linspace(-0.4, 3.5, num = 50)})
    test['Y_0'] = test['X'].map(lambda x: 0.5*x - 7.)
    test['Y_1'] = test['X'].map(lambda x: 0.5*x - 10.)
    #sns.lineplot(x=test['X'], y=test['Y_0'],color='red', ax=ax) 
    #sns.lineplot(x=test['X'], y=test['Y_1'],color='red', ax=ax) 
    #sns.scatterplot(x=new_df['x_0'], y=new_df['y_0'], c='red', marker = '*', ax=ax) 
    new_df=pd.DataFrame({'x_0': sc_x, 'y_0':sc_y})
    new_df['limit_0'] = new_df['x_0'].map(lambda x: 0.5*x - 7)
    new_df['limit_1'] = new_df['x_0'].map(lambda x: 0.5*x - 10)
    new_df = new_df[new_df['y_0']<new_df['limit_0']]
    new_df = new_df[new_df['y_0']>new_df['limit_1']]
    #sns.regplot(x=new_df['x_0'], y=new_df['y_0'], order = 4, line_kws={'color':'black'})
    popt, pcov = curve_fit(func_4, new_df['x_0'], new_df['y_0'])
    f_test = lambda x: popt[2]*x**2 + popt[1]*x + popt[0]  
    test['Z'] = test['X'].map(f_test)
    test = test[(test['X'] > -0.4)]
    test = test[(test['X'] < 3.0)]
    sns.lineplot(x=test['X'], y=test['Z'],color='white', ax=ax)
    popt_SGb  = popt.copy() 
    # fitting a 2-th polynomial to the Bright Giant Ia branch 
    test = pd.DataFrame({'X' : np.linspace(-0.4, 3.5, num = 50)})
    test['Y_0'] = test['X'].map(lambda x: 0.5*x - 10)
    #sns.lineplot(x=test['X'], y=test['Y_0'],color='red', ax=ax) 
    new_df=pd.DataFrame({'x_0': sc_x, 'y_0':sc_y})
    new_df['limit_0'] = new_df['x_0'].map(lambda x: 0.5*x - 10)
    new_df = new_df[new_df['y_0']<new_df['limit_0']]
    #sns.lineplot(x=test['X'], y=test['Y_1'],color='red', ax=ax) 
    #sns.scatterplot(x=new_df['x_0'], y=new_df['y_0'], c='red', marker = '*', ax=ax) 
    #sns.regplot(x=new_df['x_0'], y=new_df['y_0'], order = 4) 
    popt, pcov = curve_fit(func_4, new_df['x_0'], new_df['y_0']) 
    test['Z'] = test['X'].map(lambda x: popt[4]*x**4 + popt[3]*x**3 + popt[2]*x**2 + popt[1]*x + popt[0] )
    test = test[(test['X'] > -0.4)]
    test = test[(test['X'] < 3.5)]
    sns.lineplot(x=test['X'], y=test['Z'],color='white', ax=ax)
    popt_SGa  = popt.copy()
    plt.show()
    return popt_wd, popt_ps, popt_bd, popt_sd, popt_sG, popt_G, popt_SGb, popt_SGa
def find_HR_diagram(sc_x, sc_y) :
    """Based on the plot in the HR_diagram_regions(sc_x, sc_y) function, find the coefficients 
    of the fitted polynomials without ploting anything
       ____________
       Args: The B_V magnitude and the Absolute magnitude (the same used to plot the HR diagram)
       ____________
       Returns: 7 arrays, each containing the parameters of the polynomial coefficients to fit
    the HR diagram and define the Star type regions
    """
    # useful polynomials 
    def func_2(x, a_0, a_1, a_2):
        return a_0 + a_1*x + a_2*x**2
    def func_3(x, a_0, a_1, a_2, a_3):
        return a_0 + a_1*x + a_2*x**2 + a_3*x**3
    def func_4(x, a_0, a_1, a_2, a_3, a_4):
        return a_0 + a_1*x + a_2*x**2 + a_3*x**3 + a_4*x**4
    def func_10(x, a_0, a_1, a_2, a_3, a_4, a_5, a_6, a_7, a_8, a_9, a_10):
        return a_0 + a_1*x + a_2*x**2 + a_3*x**3 + a_4*x**4 + a_5*x**5 + a_6*x**6 +\
            a_7*x**7 + a_8*x**8+ a_9*x**9 + a_10*x**10
    ### fitting degree 2 polynomial on the white dwarfs region 8, 5, 0
    f_test = lambda x: 5*x + 6
    new_df=pd.DataFrame({'x_0': sc_x, 'y_0':sc_y})
    new_df['limit'] = new_df['x_0'].map(f_test)
    new_df = new_df[new_df['y_0']>new_df['limit']]
    popt, pcov = curve_fit(func_2, new_df['x_0'], new_df['y_0'])
    popt_wd  = popt.copy()
    ### fitting degree 4 polynomial on the brown dwarfs region
    new_df=pd.DataFrame({'x_0': sc_x, 'y_0':sc_y})
    new_df['limit_0'] = new_df['x_0'].map(lambda x: 15*x - 16.5)
    new_df['limit_1'] = new_df['x_0'].map(lambda x: 4.5*x - 4)
    new_df = new_df[new_df['y_0']<new_df['limit_0']]
    new_df = new_df[new_df['y_0']>new_df['limit_1']]
    new_df = new_df[new_df['x_0']>1.35]
    popt, pcov = curve_fit(func_4, new_df['x_0'], new_df['y_0'])
    popt_bd  = popt.copy()
    # fitting the subdwarf branch with 3th degree polynomial
    new_df=pd.DataFrame({'x_0': sc_x, 'y_0':sc_y})
    new_df['limit_0'] = new_df['x_0'].map(lambda x: 5*x + 6)
    new_df['limit_1'] = new_df['x_0'].map(lambda x: 6.5*x)
    new_df = new_df[new_df['y_0']<new_df['limit_0']]
    new_df = new_df[new_df['y_0']>new_df['limit_1']]
    popt, pcov = curve_fit(func_3, new_df['x_0'], new_df['y_0'])
    popt_sd  = popt.copy()
    # fitting the Main sequence stars to a 10-th order polynomial
    new_df=pd.DataFrame({'x_0': sc_x, 'y_0':sc_y})
    new_df['limit_0'] = new_df['x_0'].map(lambda x: 6.5*x)
    new_df['limit_1'] = new_df['x_0'].map(lambda x: 6.5*x - 4)
    new_df = new_df[new_df['y_0']<new_df['limit_0']]
    new_df = new_df[new_df['y_0']>new_df['limit_1']]
    new_df = new_df[new_df['x_0']>-0.3]
    new_df = new_df[new_df['x_0']<1.8]
    popt, pcov = curve_fit(func_10, new_df['x_0'], new_df['y_0'])
    popt_ps  = popt.copy()
    # fitting the subgiant branch with 2th degree polynomial
    new_df=pd.DataFrame({'x_0': sc_x, 'y_0':sc_y})
    new_df['limit_0'] = new_df['x_0'].map(lambda x: -5*(x-1.3)**2 + 3)
    new_df['limit_1'] = new_df['x_0'].map(lambda x: x - 2 )
    new_df = new_df[new_df['y_0']<new_df['limit_0']]
    new_df = new_df[new_df['y_0']>new_df['limit_1']]
    popt, pcov = curve_fit(func_2, new_df['x_0'], new_df['y_0'])
    popt_sG  = popt.copy()
    # fitting the giant branch with 4th degree polynomial
    new_df=pd.DataFrame({'x_0': sc_x, 'y_0':sc_y})
    new_df['limit_0'] = new_df['x_0'].map(lambda x: 0.5*x - 6)
    new_df['limit_1'] = new_df['x_0'].map(lambda x: 0.5*x - 3)
    new_df = new_df[new_df['y_0']>new_df['limit_0']]
    new_df = new_df[new_df['y_0']<new_df['limit_1']]
    new_df = new_df[new_df['x_0']>-0.05]
    new_df = new_df[new_df['x_0']<3.5]
    popt, pcov = curve_fit(func_4, new_df['x_0'], new_df['y_0'])
    popt_G  = popt.copy()
    # fitting a 2-th polynomial to the Bright Giant Ib branch 
    new_df=pd.DataFrame({'x_0': sc_x, 'y_0':sc_y})
    new_df['limit_0'] = new_df['x_0'].map(lambda x: 0.5*x - 7)
    new_df['limit_1'] = new_df['x_0'].map(lambda x: 0.5*x - 10)
    new_df = new_df[new_df['y_0']<new_df['limit_0']]
    new_df = new_df[new_df['y_0']>new_df['limit_1']]
    popt, pcov = curve_fit(func_4, new_df['x_0'], new_df['y_0'])
    popt_SGb  = popt.copy()
    # fitting a 2-th polynomial to the Bright Giant Ia branch 
    new_df=pd.DataFrame({'x_0': sc_x, 'y_0':sc_y})
    new_df['limit_0'] = new_df['x_0'].map(lambda x: 0.5*x - 10)
    new_df = new_df[new_df['y_0']<new_df['limit_0']]
    popt, pcov = curve_fit(func_4, new_df['x_0'], new_df['y_0']) 
    popt_SGa  = popt.copy()
    return popt_wd, popt_bd, popt_sd, popt_ps, popt_sG, popt_G, popt_SGb, popt_SGa
#--------------------------------------------------------------------------
Now, add the features that we estimate through physical arguments and the star type approximated by the distance the obct is from the characteristic regions in the HR diagram
#--------------------------------------------------------------------------
### Spectral type measurements, beggining of the estimated features
def add_luminosity(df) : 
    """Estimate the luminosity using L/Ls = 2.511886432**(Ms-M), L1_sol = 1, M1_sol = 4.83
       ____________
       Args: the dataframe (df)
       ____________
       Returns: the dataframe with the additinal column, the absolute magnitude
    """
    df['EstimatedLuminosity'] = 2.511886432**(4.83-np.array(df['AbsMag']))
    return df
def add_mass_estimative(df) : 
    """Estimate the mass using mass-luminosity estimatives like 
    class 1 : M < 0.43 M_sun -> L/Ls = 0.23(M/Ms)**2.3
    class 2 : 2 M_sun > M > 0.43 M_sun
    class 3 : 55 M_sun > M > 2 M_sun
    class 4 : M > 55 M_sun
       ____________
        Args: the dataframe (df)
       ____________
        Returns: the dataframe with the additinal column, the absolute magnitude
    """
    def mass_class_1(L) : 
        # L < 0.033
        # L/L_sun = 0.23(M/M_sun)**2.3
        mass_estimative = (4.347826087 * L)**(0.434782609)
        return mass_estimative
    def mass_class_2(L) : 
        #  16 > L > 0.033
        # L/L_sun = (M/M_sun)**4
        mass_estimative = L**(0.25)
        return mass_estimative
    def mass_class_3(L) : 
        # 1 727 418 > L > 16
        # L/L_sun = 1.4 (M/M_sun)**3.5
        mass_estimative = (0.714285714 * L)**(0.285714286)
        return mass_estimative
    def mass_class_4(L) : 
        # L > 1 727 418
        # L/L_sun = 32000 (M/M_sun)
        mass_estimative = 0.00003125 * L
        return mass_estimative
    mass_array = [] 
    for L in df['EstimatedLuminosity'] : 
        if L < 0.033 :
            mass_array.append(mass_class_1(L))
        elif L < 16 :
            mass_array.append(mass_class_2(L))   
        elif L < 1727418 :
            mass_array.append(mass_class_3(L))
        else :
            mass_array.append(mass_class_4(L))  
    df['EstimatedMass'] = mass_array
    return df
def add_radius_estimative(df) :
    """Estimate the radius of the object. For now, assuming that is a main sequence star
       ____________
       Args: the dataframe (df)
       ____________
       Returns: the dataframe with the additinal column, the estimated radius 
    """
    def main_sequence(M):
        if M < 1 : 
            R = M**0.57
        else :
            R = M**0.8
        return R    
    radius_estimates = []
    for M in df['EstimatedMass']: 
        radius_estimates.append(main_sequence(M))
    df['EstimatedRadius'] = radius_estimates
    return df
def add_effective_temperature_estimative(df) :
    """Estimate the effective temperatue using the stefan boltzman law: L = 4 pi R**2 cB T**4
       ____________
       Args: the dataframe (df)
       ____________
       Returns: the dataframe with the additinal column, the estimated effective temperatue
    """
    R_sun = 696340
    cB = 0.0567
    df['EstimatedTemperature'] = (np.array(df['EstimatedLuminosity']*3.846 * 10**26) / (4 * math.pi * (R_sun * np.array(df['EstimatedRadius']))**2 * cB))**(0.25)
    return df
def add_star_type_estimation(df) : 
    """Estimate the star type by its position in the HR diagram
       ____________
       Args: the dataframe (df)
       ____________
       Returns: the dataframe with the additinal column, the estimative of the star type
    """
    #popt_wd, popt_bd, popt_sd, popt_ps, popt_sG, popt_G, popt_SGb, popt_SGa = figure_test(temp['B-V'], temp['AbsMag'])
    popt_wd, popt_bd, popt_sd, popt_ps, popt_sG, popt_G, popt_SGb, popt_SGa =\
        find_HR_diagram(df['B-V'], df['AbsMag'])
    test = pd.DataFrame({'X' : np.linspace(df['B-V'].min(), df['B-V'].max(), num = 100)})
    test['wd_curve']= test['X'].map(lambda x: poly(x, vec=popt_wd))
    test['bd_curve']= test['X'].map(lambda x: poly(x, vec=popt_bd))
    test['swd_curve']= test['X'].map(lambda x: poly(x, vec=popt_sd))
    test['PS_curve']= test['X'].map(lambda x: poly(x, vec=popt_ps))
    test['sG_curve']= test['X'].map(lambda x: poly(x, vec=popt_sG))
    test['G_curve']= test['X'].map(lambda x: poly(x, vec=popt_G))
    test['SGb_curve']= test['X'].map(lambda x: poly(x, vec=popt_SGb))
    test['SGa_curve']= test['X'].map(lambda x: poly(x, vec=popt_SGa))
    wd_limits = (-0.4, 0.8, popt_wd)
    bd_limits = (1.4, 2.5, popt_bd) 
    swd_limits = (-0.4, 1.4, popt_sd)
    PS_limits = (-0.3, 1.8, popt_ps)
    sG_limits = (0.4, 2, popt_sG)
    G_limits = (0., 3.5, popt_G)
    SGb_limits = (-0.4, 3.5, popt_SGb) 
    SGa_limits = (-0.4, 3.5, popt_SGa)
    limits = [wd_limits, bd_limits, swd_limits, PS_limits, sG_limits, G_limits, SGa_limits, SGb_limits]
    p_wd, p_bd, p_sd, p_PS, p_sG, p_G, p_SGb, p_SGa, t = [], [], [], [], [], [], [], [], []
    StarType = ['white_dwarf', 'brown_dwarf','sub_dwarf', 'MS', 'sub_giant',\
        'giant', 'super_giant_b', 'super_giant_a']
    for index in range(df.shape[0]):
        temp = weigths(df.loc[index, 'B-V'], df.loc[index,'AbsMag'], limits)
        p_wd.append(temp[0])
        p_bd.append(temp[1])
        p_sd.append(temp[2])
        p_PS.append(temp[3])
        p_sG.append(temp[4])
        p_G.append(temp[5])
        p_SGb.append(temp[6])
        p_SGa.append(temp[7])
        i = np.where(temp==temp.max())[0][0]
        t.append(StarType[i])
    temp = pd.DataFrame({'p_wd':p_wd, 'p_bd':p_bd, 'p_sd':p_sd, 'p_PS':p_PS, 'p_sG':p_sG,\
        'p_G':p_G, 'p_SGb':p_SGb, 'p_SGa':p_SGa, 'StarType':t})   
    df = pd.concat([df, temp], axis=1)
    return df   
def add_main_sequence_type(df) :
    """Estimate the main sequence type of a main sequence star based on the estimative of the 
    effective temperature
       ____________
       Args: the dataframe (df)
       ____________
       Returns: the dataframe with the StarType column changed to have main sequence star types
    """
    for index in range(df.shape[0]):
        if  df.loc[index, 'StarType'] == 'MS' :
            if df.loc[index, 'EstimatedTemperature'] >= 25000:
                df.loc[index, 'StarType'] = 'MS_O_type'
            elif df.loc[index, 'EstimatedTemperature'] >= 11000:
                df.loc[index, 'StarType'] = 'MS_B_type'
            elif df.loc[index, 'EstimatedTemperature'] >= 7500:
                df.loc[index, 'StarType'] = 'MS_A_type'
            elif df.loc[index, 'EstimatedTemperature'] >= 6000:
                df.loc[index, 'StarType'] = 'MS_F_type'
            elif df.loc[index, 'EstimatedTemperature'] >= 5000:
                df.loc[index, 'StarType'] = 'MS_G_type'
            elif df.loc[index, 'EstimatedTemperature'] >= 3500:
                df.loc[index, 'StarType'] = 'MS_K_type'
            else:     
                df.loc[index, 'StarType'] = 'MS_M_type'
    return df         
#-------------------------------------------------------------------------
### Add properties 
def add_properties(df) : 
    """bundle up the properties wich can be directely derived from the data, 
    such as the absolute magnitude and the distances, or estimatives made with 
    physical arguments, sucha as scaling relations 
       ____________
       Args: the dataframe (df)
       ____________
       Returns: a new dataframe with the new features appended
    """
    df = add_distance(df)
    print('Distances added')
    df = add_cartesian_coordinates(df)
    print('3-dimensional cartesian coordinates added')
    df = add_abs_magnitude(df)
    print('Absolute magnitude added')
    aux = df.shape[0]
    df = df[abs(df['B-V']) > 0.000001]
    print("More data was deleted due poin float errors. An amount of {var}% was deleted".format(var=round((1-df.shape[0]/aux)*100,4)))
    df = df.reset_index()
    df = add_luminosity(df)
    print('Luminosity estimative added')
    df = add_mass_estimative(df)
    print('Mass estimative added')
    df = add_star_type_estimation(df)
    print('Gross star classification estimative added')
    df = add_radius_estimative(df)
    print('Radius estimative added')
    df = add_effective_temperature_estimative(df)
    print('Effective temperatue estimative added')
    df = add_main_sequence_type(df)
    print('Main sequence stars classification estimative added')
    return df
#-------------------------------------------------------------------------
Plot the 3D-map and some plotting functions to make a preliminary analysis of the data
### 3d neighborhood map according to the data
def plot_map(df, lower_distance, upper_distance) : 
    """Plot a 3D-map (which is navigable using the right-left mouse buttons) with 
    color defined by star type and temperature, size define by the star radius and the 
    map centered in the sun
       ____________
       Args: the dataframe (df) and two distance limits, will be plotted all stars
    distant beteween the estipulated limits
       ____________
       Returns: None
    """
    df.sort_values(by='Distance', ascending=True, inplace=True, ignore_index=True)
    temp = df[df['Distance'] < upper_distance]
    temp = temp[temp['Distance'] > lower_distance]
    temp.reset_index(drop=True, inplace=True)
    CMap = []
    Marker = []
    all_bool = []
    for index in range(temp.shape[0]): 
        if temp.loc[index, 'StarType'] == 'white_dwarf':
            CMap.append('lightgray')
            Marker.append('v')
            all_bool.append(True)
        elif temp.loc[index, 'StarType'] == 'brown_dwarf':
            CMap.append('darkorange')
            Marker.append('s')
            all_bool.append(True)  
        elif temp.loc[index, 'StarType'] == 'sub_dwarf':
            CMap.append('cadetblue')
            Marker.append('p')
            all_bool.append(True)  
        elif temp.loc[index, 'StarType'] == 'MS_O_type':
            CMap.append('blue')
            Marker.append('o')
            all_bool.append(True)   
        elif temp.loc[index, 'StarType'] == 'MS_B_type':
            CMap.append('aqua')
            Marker.append('o')
            all_bool.append(True)       
        elif temp.loc[index, 'StarType'] == 'MS_A_type':
            CMap.append('lightseagreen')
            Marker.append('o')
            all_bool.append(True)    
        elif temp.loc[index, 'StarType'] == 'MS_F_type':
            CMap.append('yellow')
            Marker.append('o')
            all_bool.append(True)     
        elif temp.loc[index, 'StarType'] == 'MS_G_type':
            CMap.append('orange')
            Marker.append('o')
            all_bool.append(True)     
        elif temp.loc[index, 'StarType'] == 'MS_K_type':
            CMap.append('orangered')
            Marker.append('o')
            all_bool.append(True)    
        elif temp.loc[index, 'StarType'] == 'MS_M_type':
            CMap.append('red')
            Marker.append('o')
            all_bool.append(True)
        elif temp.loc[index, 'StarType'] == 'sub_giant':
            CMap.append('coral')
            Marker.append('D')
            all_bool.append(True)
        elif temp.loc[index, 'StarType'] == 'giant':
            CMap.append('crimson')
            Marker.append('X')
            all_bool.append(True)   
        elif temp.loc[index, 'StarType'] == 'super_giant_b':
            CMap.append('thistle')
            Marker.append('P')
            all_bool.append(True)      
        elif temp.loc[index, 'StarType'] == 'super_giant_a':
            CMap.append('lavenderblush')
            Marker.append('*')
            all_bool.append(True)
        else:
            all_bool.append(False)     
    temp = temp[all_bool]
    fig = plt.figure()
    ax = fig.add_subplot(projection='3d')
    ax.scatter(0.0, 0.0, 0.0, s=1, marker='*', c='white')
    ax.scatter(temp['x_pos'], temp['y_pos'], temp['z_pos'], s=temp['EstimatedRadius'], \
        alpha = 1, c=CMap, marker = 'o')
    ax.grid(False)
    ax.set_xticks([])
    ax.set_yticks([])
    ax.set_zticks([])
    plt.axis('off')
    plt.grid(b=None)
    Color_BC = plt.gca()
    Color_BC.set_facecolor("black")
    plt.show()
def plot_bexenplots_degree_distance(df, feat_to_plot, inf_distance, sup_distance):
    """Plot boxenplots of a object feature in a distance interval
       ____________
       Args: the dataframe (df), the feature (to plot) name and two distance limits, will be plotted all stars
    distant between the estipulated limits
       ____________
       Returns: None
    """
    df.sort_values(by='Distance', ascending=True, inplace=True, ignore_index=True)
    temp = df[df['Distance'] <= sup_distance]
    temp = temp[temp['Distance'] >= inf_distance]
    fig, ax = plt.subplots(3,1, figsize=(16, 14))
    plt.subplots_adjust(wspace=.3, hspace=.3)
    g_0 = sns.boxenplot(data=temp, y=feat_to_plot, x = 'RArad', ax=ax[0], color = 'yellow')
    g_1 = sns.boxenplot(data=temp, y=feat_to_plot, x = 'DErad', ax=ax[1], color = 'orange')
    g_2 = sns.boxenplot(data=temp, y=feat_to_plot, x = 'Distance', ax=ax[2], color = 'violet')
    ax[0].set_xticks([-math.pi, -math.pi/2 ,0, math.pi/2, math.pi, 3*math.pi/2], ['180', '270','0', '90', '180', '270'])
    ax[1].set_xticks([-math.pi, -math.pi/2 ,0, math.pi/2, math.pi], ['180', '270','0', '90', '180'])
    ax[0].set_xlabel('Degrees')
    ax[1].set_xlabel('Degrees')
    ax[2].set_xlabel('Distance (ly)')
    plt.show()
    return None
def plot_lineplots_degree_distance(df, feat_to_plot, inf_distance, sup_distance):
    """Plot lineplots of a numerical feature in a distance interval
       ____________
       Args: the dataframe (df), the feature (to plot) name and two distance limits, will be plotted all stars
    distant between the estipulated limits
       ____________
       Returns: None
    """
    df.sort_values(by='Distance', ascending=True, inplace=True, ignore_index=True)
    temp = df[df['Distance'] <= sup_distance]
    temp = temp[temp['Distance'] >= inf_distance]
    fig, ax = plt.subplots(3,1, figsize=(16, 14))
    plt.subplots_adjust(wspace=.3, hspace=.3)
    g_0 = sns.lineplot(data=temp, y=feat_to_plot, x = 'RArad', ax=ax[0], color = 'yellow')
    g_1 = sns.lineplot(data=temp, y=feat_to_plot, x = 'DErad', ax=ax[1], color = 'orange')
    g_2 = sns.lineplot(data=temp, y=feat_to_plot, x = 'Distance', ax=ax[2], color = 'violet')
    ax[0].set_xticks([-math.pi, -math.pi/2 ,0, math.pi/2, math.pi, 3*math.pi/2], ['180', '270','0', '90', '180', '270'])
    ax[1].set_xticks([-math.pi, -math.pi/2 ,0, math.pi/2, math.pi], ['180', '270','0', '90', '180'])
    ax[0].set_xlabel('Degrees')
    ax[1].set_xlabel('Degrees')
    ax[2].set_xlabel('Distance (ly)')
    plt.show()
    return None
#-------------------------------------------------------------------------



  • AI Chat
  • Code