### Instructions

When running the notebook the first time, make sure to run all cells before making changes in the notebook. Hit Shift + Enter to run the selected cell or, in the top menu, click on: `Kernel` > `Restart Kernel and Run All Cells...` to rerun the whole notebook. If you make any changes in a cell, rerun that cell.

If you make any changes in a coding cell, rerun the notebook by `Run` > `Run Selected Cell and All Below`


```python
# Importing libraries and dependencies
import sys
sys.path.append('python/')
import NGC5533_functions as funcNGC5533        # Functions for the components of NGC 5533
import NGC5533_traced_curves as dataNGC5533    # Traced data of Galaxy NGC 5533
import NGC5533_fitting as fitNGC5533           # Fitting parameters for best fit values
import dataPython as dp
import numpy as np
import matplotlib.pyplot as plt                # Plotting
%matplotlib inline
from ipywidgets import interactive, fixed, FloatSlider, HBox, Layout, Button, Label, Output, VBox   # Widget
from IPython.display import display, clear_output
from IPython.display import Javascript
import scipy.stats as stats
import warnings
warnings.filterwarnings("ignore")              # Ignore warnings
```

# Interactive Widget with Variable Components - Galaxy: NGC 5533

In the following activity, you can adjust the ingredients of a galaxy and investigate how the total velocity of stars and gases change. 

### Measured data points

Import the measured data points for the NGC 5533 galaxy. These are measured velocities of stars and gas in the galaxy as a function of radius. The data points were traced from Figure 4. in Noordermeer (2008). 


```python
data = dp.getXYdata_wXYerr('data/NGC5533/noord-120kpc-datapoints.txt')
r_dat = np.asarray(data['xx'])
v_dat = np.asarray(data['yy'])
v_err0 = np.asarray(data['ex'])
v_err1 = np.asarray(data['ey'])
```

### Components

The stars and gas in a galaxy can be categorized into components depending on the shape of the distribution. These components are then added in quadrature to calculate the total velocity or the rotation curve of a galaxy. In the case of NGC 5533, we can define a central supermassive black hole, a central bulge, a flat disk, a mostly atomic gas and a dark matter halo component. 

>__Total velocity__: <br>
    \begin{equation}
    v_{total}(r) = \sqrt{v_{gas}^2 + \Upsilon _{bulge} v_{bulge}^2 + \Upsilon _{disk} v_{disk}^2 + v_{halo}^2}
    \end{equation}<br>

Lets import the relevant functions or traced curves from the `NGC5533_functions.py` library below:


```python
def blackhole(r,Mbh):
    return funcNGC5533.bh_v(r,Mbh,load=False)

def bulge(r,bpref):
    return bpref*funcNGC5533.b_v(r,load=True)

def disk(r,dpref):
    return dpref*funcNGC5533.d_thief(r)

def gas(r,gpref):
    return gpref*funcNGC5533.g_thief(r)

def halo(r,rc,rho0):
    return funcNGC5533.h_v(r,rc,rho0,load=False)

# Total velocity containing all components
def total_all(r,Mbh,bpref,dpref,gpref,rc,rho0):
    total = np.sqrt(blackhole(r,Mbh)**2 
                    + bulge(r,bpref)**2 
                    + disk(r,dpref)**2
                    + gas(r,gpref)**2
                    + halo(r,rc,rho0)**2)
    return total

# Total velocity of baryonic or luminous matter (no dark matter component)
def total_bary(r,Mbh,bpref,dpref,gpref):
    total = np.sqrt(blackhole(r,Mbh)**2 
                    + bulge(r,bpref)**2 
                    + disk(r,dpref)**2
                    + gas(r,gpref)**2)
    return total
```

### Parameters

The scaling parameters for each components can be found by fitting the total velocity to the measured data points. Import these fitting parameters for our widget from the _NGC5533\_fitting.py_ library.  


```python
best_Mbh = fitNGC5533.f_Mbh
best_bpref = fitNGC5533.f_bpref
best_dpref = fitNGC5533.f_dpref
best_gpref = fitNGC5533.f_gpref
best_rc = fitNGC5533.f_rc
best_rho0 = fitNGC5533.f_rho00
```

### Define plotting function and sliders for widget

The interactive widget is defined with six adjustable sliders that controls the parameters of each component of a galaxy. 


```python
# Plotting function
def f(Mbh,bpref,dpref,gpref,rc,rho0):
        
    # Define radius
    r = np.linspace(np.min(r_dat),np.max(r_dat),1000)
    
    # Plot
    plt.figure(figsize=(11,7))
    plt.xlim(0,np.max(r_dat)+0.2)
    plt.ylim(0,np.max(v_dat)+100)
    
    plt.errorbar(r_dat,v_dat,yerr=v_err1,fmt='bo',label='Data')     # Measured data points
    plt.plot(r,blackhole(r,Mbh),label=("Central Supermassive Black Hole"),color='black') # Black hole component
    plt.plot(r,bulge(r,bpref),label=("Bulge"),color='orange')       # Bulge component
    plt.plot(r,disk(r,dpref),label=("Disk"),color='purple')         # Disk component 
    plt.plot(r,gas(r,gpref),label=("Gas"),color='blue')             # Gas component
    plt.plot(r,halo(r,rc,rho0),label=("Halo"),color='green')        # Dark matter halo component
    plt.plot(r,total_all(r,Mbh,bpref,dpref,gpref,rc,rho0),label=("Total Curve"),color='red')    # Total velocity with dark matter
    plt.plot(r,total_bary(r,Mbh,bpref,dpref,gpref),label=("Luminous Matter"),linestyle='--')    # Total velocity without dark matter
    plt.fill_between(r,dataNGC5533.greyb_bottom(r),dataNGC5533.greyb_top(r),color='#dddddd',label="Confidence Band")    # Confidence band
    plt.title("Interactive Rotation Curve - Galaxy: NGC 5533",fontsize=16)
    plt.xlabel("Radius (kpc)")
    plt.ylabel("Velocity (km/s)")
    plt.legend(bbox_to_anchor=(1,1), loc="upper left")              # Put legend outside of the plot
    
    # Chi squared and reduced chi squared
    # Residuals
    residuals = v_dat - total_all(r_dat,Mbh,bpref,dpref,gpref,rc,rho0)
    # Error
    error = np.sqrt(v_err1**2 + dataNGC5533.band**2)
    # Chi squared
    chisquared = np.sum(residuals**2/error**2)
    # Degrees of freedom
    dof = len(r_dat) - 6                 # number of degrees of freedom = number of observed data - number of fitting parameters
    # Reduced chi squared
    reducedchisquared = chisquared / dof 
    
    # Annotation
    props = dict(boxstyle='round', facecolor='white', alpha=0.5)
    plt.text(80,373,r"$\chi^2$: {:.5f}".format(chisquared)+'\n'+r"Reduced $\chi^2$: {:.5f}".format(reducedchisquared),bbox=props,size=10)
    props = dict(boxstyle='round', facecolor='white', alpha=0.5)
    plt.annotate('Data source: E. Noordermeer. The rotation curves of flattened Sérsic bulges. MNRAS,385(3):1359–1364, Apr 2008',
            xy=(0, 0), xytext=(0,5),
            xycoords=('axes fraction', 'figure fraction'),
            textcoords='offset points',
            size=10, ha='left', va='bottom')
    
    plt.show()
```


```python
# Widget appearance 
style = {'description_width': 'initial'}
layout = {'width':'600px'}

# Define widget sliders
# Mass of central supermassive black hole
Mbh = FloatSlider(min=0, max=5e9, step=1e8, 
                value=best_Mbh, 
                description='Black Hole Mass [$M_{\odot}$]', 
                readout_format='.2e', 
                orientation='horizontal', 
                style=style, layout=layout)

# Bulge prefactor (multiplier)
bpref = FloatSlider(min=0, max=5, step=0.1, 
                    value=round(best_bpref,1), 
                    description='Bulge Prefactor', 
                    readout_format='.2f', 
                    orientation='horizontal', 
                    style=style, layout=layout)

# Disk prefactor (multiplier)
dpref = FloatSlider(min=0, max=5, step=0.1, 
                    value=round(best_dpref,1), 
                    description='Disk Prefactor', 
                    readout_format='.2f', 
                    orientation='horizontal', 
                    style=style, layout=layout)

# Gas prefactor (multiplier)
gpref = FloatSlider(min=0, max=5, step=0.1, 
                    value=round(best_gpref,1), 
                    description='Gas Prefactor', 
                    readout_format='.2f', 
                    orientation='horizontal', 
                    style=style, layout=layout)

# Core radius of dark matter halo
rc = FloatSlider(min=0, max=5, step=0.1, 
                 value=round(best_rc,1), 
                 description='Halo Core Radius [kpc]', 
                 readout_format='.2f', 
                 orientation='horizontal', 
                 style=style, layout=layout)
#rc = fixed(best_rc)      # If we need the core radius to be fixed

# Central mass density of dark matter halo
rho0 = FloatSlider(min=0, max=1e9, step=1e7, 
                    value=best_rho0, 
                    description=r'Halo Central Mass Density [$M_{\odot} / kpc^3$]', 
                    readout_format='.2e', 
                    orientation='horizontal', 
                    style=style, layout=layout)

# Interactive widget
def interactive_plot(f):
    interact = interactive(f, Mbh = Mbh, 
                           bpref = bpref, 
                           dpref = dpref, 
                           gpref = gpref,
                           rc = rc,
                           rho0 = rho0,
                           continuous_update=False)
    return interact

# Button to revert back to Best Fit
button = Button(
    description="Best Fit",
    button_style='warning',
    icon='check')
out = Output()

# Define what happens when we click on the "Best Fit" button
def on_button_clicked(_):
    Mbh.value = best_Mbh
    bpref.value = round(best_bpref,1)
    dpref.value = round(best_dpref,1)
    gpref.value = round(best_gpref,1)
    rc.value = round(best_rc,1)
    rho0.value = best_rho0

button.on_click(on_button_clicked)
```

Using the sliders below, answer the following questions: 

<div class="alert-info">Activity 1)</div>

>Can you get the Reduced $\chi^2$ close to zero?

<div class="alert-info">Activity 2)</div>

>Can you get the Reduced $\chi^2$ close to one?

<div class="alert-info">Activity 3)</div>

>Remove the Dark Matter by changing both Dark Matter parameters (halo core radius and halo central mass density) to zero using the sliders. Are you able to scale the other components (central black hole, bulge, disk, gas) to have the total curve aligned with the measured data points (Reduced $\chi^2$ close to one)?

In a strictly statistical sense, a Reduced $\chi^2$ of zero represents a perfect fit. However, a Reduced $\chi^2$ much less than 1 often indicates that more parameters were used than represented by the physical system. We would be simply connecting the measured data points. <br>
In general, the gas component of a galaxy can be measured with great accuracy. It is made of mostly atomic gas. In most cases, the measured gas component is multiplied by 1.33 to account for helium. In spite of that, the gas component can also be adjusted in this widget. Here, we allow you to scale each component as you please. See how 'good' of a fit you can get!


```python
#NBVAL_IGNORE_OUTPUT
#Because the figure doesn't save to the repository correctly.
# Widget
VBox([button,out,interactive_plot(f)])
```


    VBox(children=(Button(button_style='warning', description='Best Fit', icon='check', style=ButtonStyle()), Outp…


### Slider key

The prefactor for the **black hole** is the square root of its mass. 
Because it is more physically meaningful and not computationally intensive, we have you control the mass of the black hole rather than the prefactor itself.
The units of the black hole's mass are solar masses (multiply the value by the mass of our Sun. Note that 'e+9' equates to a billion--that's a billion solar masses!).

The prefactors for the **bulge** and the **disk** relate to the mass-to-light ratio.
A prefactor of one means that you are using the M/L we generally expect from stars. Increasing the prefactor means increasing the expected mass for a given amount of light. 

The halo central mass density behaves as a prefactor for the **Dark Matter "halo"** in and around the galaxy. 
Its units are solar masses per cubic kiloparsec and, in this fit, is on the scale of hundreds of millions.
This represents how much Dark Matter we think there is overall in the center of a galaxy. 

The **gas** prefactor is similar to the bulge and disk prefactors. 
Rather than being for a stellar M/L, it represents the amount of radiation emitted by interstellar gas. 
A value of 1 means using the expected amount of radiation. Increasing the prefactor means increasing the expected mass for a given amount of radiation (note that this is why the gas component is kept static during fitting, as increasing or decreasing the gas prefactor would be digressing from the known physical laws of H1 emissions).

### References: 

>Jimenez, Raul, Licia Verde, and S. Peng Oh. **Dark halo properties from rotation curves.** _Monthly Notices of the Royal Astronomical Society_ 339, no. 1 (2003): 243-259. https://doi.org/10.1046/j.1365-8711.2003.06165.x. <br><br>
>Noordermeer, E., &amp; Van Der Hulst, J. M. (2007). **The stellar mass distribution in early-type disc galaxies: Surface Photometry and bulge–disc decompositions.** Monthly Notices of the Royal Astronomical Society, 376(4), 1480–1512. https://doi.org/10.1111/j.1365-2966.2007.11532.x <br><br>
>Noordermeer, E. (2008), **The rotation curves of flattened Sérsic bulges**. Monthly Notices of the Royal Astronomical Society, 385: 1359-1364. https://doi.org/10.1111/j.1365-2966.2008.12837.x <br><br>
>Taylor, John Robert. __An Introduction to Error Analysis: The Study of Uncertainties in Physical Measurements.__ 2nd ed. United States of America: University Science Books, 1997.

ce9d5a8cae98d2be6a3f9628193a8c73c97380c7
ipynb
Jupyter Notebook
binder/05_Widget_NGC5533_All_Components.ipynb
villano-lab/galactic-spin-W1
d95c706ccbc347f9bc61bb7c96b1314460bc2d0f
CC-BY-4.0
2022-03-22T04:00:17.000Z
2022-03-22T04:00:17.000Z
binder/05_Widget_NGC5533_All_Components.ipynb
villano-lab/galactic-spin-W1
d95c706ccbc347f9bc61bb7c96b1314460bc2d0f
CC-BY-4.0
2021-11-05T18:17:19.000Z
2022-02-19T20:35:05.000Z
binder/05_Widget_NGC5533_All_Components.ipynb
villano-lab/galactic-spin-W1
d95c706ccbc347f9bc61bb7c96b1314460bc2d0f
CC-BY-4.0
Qwen/Qwen-72B
1. YES
2. YES