Prerequisites#

These exercises assume:

1. Basic knowledge of python coding

2. Basic knowledge of RA/DEC cooridinates for finding star clusters

3. Basic knowledge of how color magnitude diagrams are generated. More information on Star Clusters from this website HR Diagrams and Cluster Evolution.

Introduction - Data Reduction#

In this activity, we will demonstrate how code can be used to query data from Gaia Data Release 3 to generate a color magnitude diagram. However, when it comes to extracting data from any database, it is important to be sure you are evaluating the data from your source. Often, data may have some values that do not make sense, such as a magnitude of 40 for a star. That value would be too faint for our telescopes, therefore that would skew the data. Second, sometimes measurements can’t be made and values such as -999 is point in its place. This would again skew the data, giving values that does not make sense.

We are going to look for stars at a specific cooridinate and do a radial search around that cooridinate. This poses another problem. We will be targeting stars both in our cluster and outside the cluster. This makes it difficult to generate the HR Diagram because this only works for stars in the same cluster which are at the same distance away. To identify stars in the cluster, will we target stars with the same proper motions. Using this fact, we can filter out the field stars and target the cluster starts. Once we are able to filter out those field stars, one can generate a HR Diagram of the cluster and begin analysis on the cluster. This activity will take you through these steps and validate with an animation that the cluster stars do move together with the same proper motion, as opposed to the field stars which vary in velocity.

Set up your Machine#

It is important to run the following installs for the necessary libraries. Whether you are using a virtual machine like Colab, or a non-virtual machine like jupyter notebooks, you must install and set up the following. Copy and paste the code into your terminal.

Specifically for Google Colab: When running the cell above, Notice that the infinity symbol at the top of the Colab turns gray, and turns back to yellow when it finishes. While it is running, click "Show code" to see what it is doing. The "install" statements are installing software on your virtual machine, and the "import" statements are importing functionality that you will need to run the code below. If you leave the Colab, Google might reallocate the virtual machine you are using. In this case, when you return, you would get a new virtual machine, which does not have the software installed. If something breaks, the infinity symbol at the top of the Colab will turn red and an error message will be displayed below the corresponding cell.

#@title Installs and auxiliary functions
%pip install numpy
%pip install astroquery
%pip install astropy
%pip install pandas
%pip install ipywidgets
%pip install matplotlib
from astroquery.gaia import Gaia
from astropy.table import Table
import astropy.units as u
from matplotlib import pyplot as plt
import numpy as np
from matplotlib.animation import FuncAnimation
from IPython.display import HTML, display
import pandas as pd
import ipywidgets as widgets
from matplotlib import rc
from astroquery.vizier import Vizier

def interactive_display_astropy_table(astropy_table):
    """
    Interactive illustration of an Astropy table in a Jupyter Notebook using Pandas.

    Parameters:
    astropy_table (astropy.table.Table): The Astropy table to be displayed.
    """
    # Convert Astropy Table to Pandas DataFrame
    df = astropy_table.to_pandas()

    # Create widgets for interactive display
    search_widget = widgets.Text(
        description='Search:',
        placeholder='Type here...',
        style={'description_width': 'initial'}
    )

    dropdown_widget = widgets.Dropdown(
        options=['All'] + list(df.columns),
        description='Filter by:',
        style={'description_width': 'initial'}
    )

    output = widgets.Output()

    def update_table(change):
        with output:
            output.clear_output()
            if search_widget.value:
                if dropdown_widget.value == 'All':
                    mask = df.apply(lambda row: row.astype(str).str.contains(search_widget.value, case=False).any(), axis=1)
                else:
                    mask = df[dropdown_widget.value].astype(str).str.contains(search_widget.value, case=False)
                display(df[mask])
            else:
                display(df)

    search_widget.observe(update_table, names='value')
    dropdown_widget.observe(update_table, names='value')

    # Display the widgets and the table
    display(widgets.VBox([widgets.HBox([search ```_widget, dropdown_widget]), output]))
    update_table(None)

Query Gaia#

In the following exercise, you’ll query gaia databases to retrieved the information we need thoughout this activity.

Photometry & CMD#

We will now plot the Color-Magnitude Diagram for the photometry data you extracted from the star cluster you searched for in Gaia. You will notice the standard clustering of points indicating the HR Diagram for star clusters, however, you will also notice a fair amount of stars in the diagram that are far away from this relation. These are known as field stars. These are the stars in the same region of the sky as the cluster, but which are not actually part of the cluster.

Field Star Removal with Astrometry Plots#

Next we need to remove the field stars from our HR plot. We identify cluster stars through their proper motions. These proper motions in the RA and DEC direction will all be in the same if they are part of the cluster, as those stars are all moving (to first order) in the direction and at the velocity that the cluster as a whole is moving.

Visulization of Cluster Motion#

The code below will show how cluster members move within a star field. We used data from above to identify star cluster members through their proper motion. Here, this shows all of the star proper motions.

def plot_cluster_dynamic(_table, _cluster_table, _field_table, _mf, _speed):
    """
    Plot the cluster dynamically based on the proper motion.

    Parameters:
    _table (astropy.table.Table): The table containing the data.
    _speed (float): The speed factor for the animation. sec/yr
    """

    # Append the data
    ra = _table["ra"]
    dec = _table["dec"]
    pmra = _table["pmra"]
    pmdec = _table["pmdec"]

    # range of ra/dec
    _ra_range = max(ra) - min(ra)
    _dec_range = max(dec) - min(dec)
    _min_plot_ra = min(ra) - _mf * _ra_range
    _max_plot_ra = max(ra) + _mf * _ra_range
    _min_plot_dec = min(dec) - _mf * _dec_range
    _max_plot_dec = max(dec) + _mf * _dec_range

    # Set up the plot
    fig, ax = plt.subplots(figsize=(10, 10))
    ax.set_xlim(_min_plot_ra, _max_plot_ra)
    ax.set_ylim(_min_plot_dec, _max_plot_dec)
    ax.set_xlabel("RA (degrees)")
    ax.set_ylabel("Dec (degrees)")
    ax.set_title("Cluster with Proper Motion Animation")

    # Create the scatter plot
    cluster_ra = _cluster_table["ra"]
    cluster_dec = _cluster_table["dec"]
    ax.scatter(cluster_ra, cluster_dec, s=3, color="red")
    field_ra = _field_table["ra"]
    field_dec = _field_table["dec"]
    ax.scatter(field_ra, field_dec, s=1, color="grey")

    _pm_factor = 0.001 / 3600  # milli-arsec per year to degree per year

    _delta_ra_cluster = _cluster_table["pmra"] * _pm_factor * _speed
    _delta_dec_cluster = _cluster_table["pmdec"] * _pm_factor * _speed
    _delta_ra_field = _field_table["pmra"] * _pm_factor * _speed
    _delta_dec_field = _field_table["pmdec"] * _pm_factor * _speed

    # Update function for the animation
    def update(_i):
        nonlocal cluster_ra, cluster_dec, field_ra, field_dec, ax
        ax.clear()
        ax.set_xlim(_min_plot_ra, _max_plot_ra)
        ax.set_ylim(_min_plot_dec, _max_plot_dec)
        ax.set_xlabel("RA (degrees)")
        ax.set_ylabel("Dec (degrees)")
        ax.set_title("Cluster with Proper Motion Animation")
        ax.scatter(cluster_ra, cluster_dec, s=3, color="red")
        ax.scatter(field_ra, field_dec, s=1, color="grey")
        cluster_ra += _delta_ra_cluster
        cluster_dec += _delta_dec_cluster
        field_ra += _delta_ra_field
        field_dec += _delta_dec_field

    # Animate the plot
    anim = FuncAnimation(fig, update, frames=100, interval=100, repeat=True)
    # save as gif
    # anim.save("cluster.gif", fps=10)
    plt.close()
    return anim


# Example usage
anim = plot_cluster_dynamic(table, cluster_table, field_table, 0.5, 2000)

# equivalent to rcParams['animation.html'] = 'html5'
rc('animation', html='html5')
anim

Calculate T_eff and Luminosity for cluster stars#

The code below will calculate the temperature and lumonisity of the cluster members to identify their characteristics. They will be plotted as well. The process to calculate these features are identified and commented out in the code.