Interactives/Synthetic_Galaxy_HI_Spectra.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Synthetic HI Spectra for Spiral Galaxies\n",
    "\n",
    "### By Juan Cabanela (Minnesota State University Moorhead)\n",
    "\n",
    "This program is designed to allow you to model the neutral hydrogen line spectrum (the so-called 21-cm line) of a model spiral galaxy as viewed from outside the galaxy.  In other words, the model radio spectrum of spiral galaxies.  We do this for\n",
    "1. Unresolved \"single-dish\" Arecibo Observatory spectra where we assume the entire galaxy falls within the radio beam (assuming a 145 arcsecond beam)with a 5 km/s velocity resolution\n",
    "2. Resolved aperture synthesis spectra (as might be produced by the VLA) in which we can resolve the galaxy, essentially getting the spectrum for every point within the galaxy.\n",
    "\n",
    "This code generates spectra for model spiral galaxies with the HI gas density profile (the density versus radius) and rotation curve initially assumed to be similar to the Milky Way. However, you can sketch any rotation curve or HI gas density profile on those plots and the modelled unresolved \"single-dish\" spectrum and the resolved \"aperture synthesis\" velocity map are updated.  You can also modify the redshift $cz$ or gas velocity dispersion by sliding thei respective sliders."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Author: Juan Cabanela\n",
    "#\n",
    "# Originally written in perl (using pgplot libraries) in 2002, this could was ported to a python script \n",
    "# using matplotlib in 2013 (roughly) and in 2017 it was moved to a Jupyter notebook.  The current version\n",
    "# was ported in summer 2019 to allow use of ipywidget widgets as controls for the l-v diagram generator."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import necessary libraries\n",
    "import matplotlib.cm as cm\n",
    "import matplotlib as mpl\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "import random\n",
    "import string\n",
    "import numpy as np\n",
    "import io\n",
    "import sys\n",
    "import bqplot as bq\n",
    "from ipywidgets import Layout\n",
    "from ipywidgets import widgets\n",
    "from IPython.display import display\n",
    "\n",
    "# Import specific functions designed for galaxy model building\n",
    "from galaxyparam import *"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Configure some basic settings for this notebook here\n",
    "#\n",
    "# You will need to set up a few things before running this Jupyter notebook (the defaults are fine for the first run):\n",
    "#\n",
    "# Changing the galaxy model: You may want to edit galaxy_datafile with the filename of the file describing the \n",
    "#    galaxy.  That file must be in the same directory with the IPython notebook.  That file is a 3 column (where \n",
    "#    the columns are comma-seperated) text file.  The three columns are:\n",
    "#\n",
    "#    1. radius (in kpc): Radius of this data point in kpc\n",
    "#    2. rot_vel (in km/s): Rotational velocity of this data point in km/s\n",
    "#    3. density (in atoms/cm^3): Assumed gas density at this point (in atoms/cm^3)\n",
    "#    The 'header' row is required to be \"radius,rot_vel,density\" or this script will not properly read the file.\n",
    "#\n",
    "#   After running the default realistic Milky Way model, you will create a new galaxy model (in a different file\n",
    "#   than the original galaxy_description.csv) to represent a Keplerian (~ r^{-1/2}) rotation curve.  \n",
    "#   You will need to set galaxy_datafile to point to the file containing the new model.  Don't worry about \n",
    "#   evenly spacing data points in radius.  The part of the script below that reads this datafile will \n",
    "#   interpolate between your radial points with a spline fit, so it is pretty robust as smoothly fitting \n",
    "#   between your data points.\n",
    "#\n",
    "# Changing the image filename: Set up the image_filename for the file you will want to save the final image.\n",
    "#    It must end with .jpg, .png, or .pdf.\n",
    "\n",
    "# Initialize filenames and parameters related to data import\n",
    "galaxy_datafile = \"data/galaxy_description.csv\"  # File containing rotation curve and density data\n",
    "HIradius = 25    # Assumed radius of HI disk\n",
    "dens_drop = 0.01  # Drop in density per kpc of radius (in atoms/cm^3)\n",
    "# Load the Galaxy Data (go out sqrt(2) times HI radius to allow for square computational grid later)\n",
    "(rad_raw, rotvel_raw, density_raw) = LoadData(galaxy_datafile, np.sqrt(2)*HIradius, dens_drop) \n",
    "\n",
    "image_directory = \"tmp\"  # Name of temporary directory to store images in\n",
    "\n",
    "randkey = ''.join(random.choice(string.ascii_letters) for i in range(5) )\n",
    "basename = \"Synthetic_HISpectra\"\n",
    "image_filename = \"{0}_{1}.png\".format(basename, randkey)     # Name of the final PNG image file to generate\n",
    "\n",
    "#\n",
    "# Constants related to the spline fit routine\n",
    "#\n",
    "dr_spline = 0.1    # Step size in kpc for spline fit (global, so it can be used within functions)\n",
    "\n",
    "# Make the initial spline fit to the rotation curve and gas density curve data\n",
    "rad_spline, rotvel_spline, density_spline = spline_curves(rad_raw, rotvel_raw, density_raw, dr_spline)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "############################################\n",
    "# Initialize Constants and other variables #\n",
    "############################################\n",
    "deg2rad = np.pi/180  # Radians per degree\n",
    "\n",
    "# Constants related to computing line of sight velocity\n",
    "# of different points in the galaxy's disk.\n",
    "radial_step = 0.05   # Step size [in kpc]\n",
    "theta_step = 1       # Step size [in degrees ]\n",
    "vel_step = 2         # Velocity bin size [in km/s]\n",
    "\n",
    "# Additional constants defined for modeling the HI \n",
    "# spectrum of a galaxy as seen from the outside.\n",
    "incl = 75.0     # Inclination of galaxy between 0 and 85 [degrees] \n",
    "cz   = 5000.0   # Redshift of galaxy in km/s\n",
    "\n",
    "H_0 = 72      # Assumed value of Hubble's constant (in km/s/Mpc)\n",
    "\n",
    "dx = 0.5           # Stepsize for integrating over galaxy (x,y) grid (in Kpc)\n",
    "dz = 1             # Assumed thickness of galaxy here (in Kpc)\n",
    "vel_sigma = 17     # Assumed FWHM of gas velocity distribution width (in km/s)  \n",
    "\n",
    "dv = 5              # Channel width of unresolved dish spectra [Arecibo@21cm] (in km/s)\n",
    "rstep = 0.05        # Stepsize in radius (kpc) for fit of rotation & density curve\n",
    "velocity_step = 20  # Velocity step size in contour plots\n",
    "dens_step = 0.1     # Density step size in contour plots\n",
    "half_width = 1000   # Allowed velocity difference from cz (in km/s)\n",
    "maxdens = 1.0       # Maximum allowed gas density (in atoms/cm^2)\n",
    "maxvel = 400        # Maximum rotational velocity allowed (in km/s)\n",
    "\n",
    "# Constants related to converting gas density from atoms/cm^3 to\n",
    "# solar masses/Kpc^3.\n",
    "vol_factor = (3.09e21)*(3.09e21)*(3.09e21)  # cm^3/Kpc^3\n",
    "mass_factor = 1.67e-27                      # kg/atom\n",
    "solar_factor = 1.99e30                      # kg/M_solar\n",
    "factor = vol_factor*mass_factor/solar_factor  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "##########################\n",
    "# FUNCTION LIBRARY BELOW #\n",
    "##########################"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def check_beam(cz, H_0, radius, beam_fwhm = 145):\n",
    "    # Compute size of beam in Kpc at galaxy and confirm galaxy unresolved.\n",
    "    # Assumed default beam size (beam_fwhm) is Arecibo (in arcseconds).\n",
    "    global cz_notice\n",
    "    \n",
    "    # Assuming a sharp edged radio beam, which is not realistic...\n",
    "    beam_galaxy = (beam_fwhm/206265.)*(cz/H_0)*1000.\n",
    "    \n",
    "    # Check this beam radius versus galaxy radius and report back\n",
    "    if (beam_galaxy < radius): \n",
    "        cz_notice.value = \"<p style=\\\"color: red; font-size: smaller;\\\"><B>WARNING:</B> Galaxy is resolved by Arecibo beam!!!</P>\"\n",
    "    else:\n",
    "        cz_notice.value = \"\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def build_galaxy_velcube(rotvel_prof, density_prof, HIradius, dx, dr_spline):\n",
    "    #\n",
    "    # This function builds a popultion of noisy velocity cubes (velocity \n",
    "    # versus position on sky) for a galaxy, along with a (simple) \n",
    "    # density cube.  It returns that data.\n",
    "    #\n",
    "    # Refactored the previous approach of looping over all the \n",
    "    # x,y positions in the galaxy's plane and converting them individually\n",
    "    # into polar (r, theta) positions which allow rotational velocity\n",
    "    # and gas density lookup.  The major change was to use array \n",
    "    # arithmetic and thus speed up the computations immensely versus \n",
    "    # inefficient loops.\n",
    "    #\n",
    "    # This was done during the migration to an ipywidgets interface in\n",
    "    # summer 2019.\n",
    "\n",
    "    # Build arrays of positions in the galaxy in cartesian coordinates\n",
    "    # and then in polar coordinates.\n",
    "    radspace = np.linspace(-HIradius, HIradius, int(((2*HIradius)/dx)+1))\n",
    "    x, y_gal = np.meshgrid(radspace, radspace)\n",
    "    (r, theta) = RotCoord(x, y_gal)\n",
    "    # Compute position on plane in the sky (assuming thin disk, so no z component)\n",
    "    y = y_gal*np.cos(incl*deg2rad)\n",
    "\n",
    "    # Retrieve HI density and rotational velocities for all galaxy positions\n",
    "    R_idx = np.floor(r/dr_spline).astype(np.int64)\n",
    "    dens_0 = density_prof[R_idx] # Faceon density at each point\n",
    "    vel  = rotvel_prof[R_idx]  # Rotational velocity at each point\n",
    "\n",
    "    # Project these rotational velocities to line of sight velocities.\n",
    "    # Assumes positive y-axis pointed toward the earth if the inclination \n",
    "    # angle i was 90 degrees.\n",
    "    v_rad = vel*np.cos(theta)*np.sin(incl*deg2rad)\n",
    "\n",
    "    # Generate population of n_copies copies of the velocities with \n",
    "    # a Gaussian scatter of velocities added to them.  Also make copies\n",
    "    # of the gas densities and use that to compute the HI masses involved.\n",
    "    n_copies = 20\n",
    "    v_rad_pop = np.array([v_rad]*n_copies)\n",
    "    v_rad_pop += np.random.normal(loc=0, scale=vel_sigma, size=v_rad_pop.shape)\n",
    "    dens_0_pop = np.array([dens_0]*n_copies)\n",
    "    \n",
    "    return (x, y, v_rad_pop, dens_0_pop)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def unresolved_spectra(vrad_pop, dens_pop, cz):\n",
    "    global H_0, half_width, dv\n",
    "    \n",
    "    #####################################\n",
    "    # GENERATING THE UNRESOLVED SPECTRA #\n",
    "    #####################################\n",
    "\n",
    "    # Construct numpy array for single dish spectram storage\n",
    "    chvel = np.arange(cz-half_width, cz+half_width+dv, dv)\n",
    "    MassHisto = np.zeros_like(chvel)\n",
    "\n",
    "    # Determine the unresolved spectra velocity bin (in array form)\n",
    "    # and then sum of the mass contributions.\n",
    "    velbin1D_pop = np.floor((vrad_pop + half_width)/dv).astype(np.int64)\n",
    "    vels1D = velbin1D_pop.reshape(velbin1D_pop.size)\n",
    "    \n",
    "    # Determine mass distribution in the galaxy\n",
    "    dens = dens_pop.reshape(dens_pop.size)\n",
    "    mass = factor*dens*dx*dx*dz   # mass of these particular copies\n",
    "\n",
    "    # Record these density contributions to approprate velbin and lonbin of dens_map\n",
    "    for i in range(vels1D.size):\n",
    "            MassHisto[vels1D[i]] += mass[i]\n",
    "    total_mass= np.sum(mass)   # Track the total gas mass in solar masses\n",
    "\n",
    "    # At this point, we have a Mass histogram with the amount of HI mass in each\n",
    "    # velocity bin (chvel).  We convert this into a flux array by inverting the \n",
    "    # flux to hydrogen mass equation:\n",
    "    #    HI_Mass (M_sun) = (235600)*(Dist[Mpc]^2)*Integral(F*delV Jy-km/s)\n",
    "    # such that (assuming a boxcar profile for the gas in each velocity channel)\n",
    "    #    F = HI_Mass/(delV*235600*Dist^2) \n",
    "    mass2Jy = (dv*235600.*(cz/H_0)*(cz/H_0)) # conversion in Msolar/Jy-km/s\n",
    "    unresolved_flux = MassHisto/mass2Jy\n",
    "    \n",
    "    return(chvel, unresolved_flux)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def resolved_spectra(cz, vrad_pop, dens_pop, velocity_step):\n",
    "    global incl, deg2rad\n",
    "    \n",
    "    ##########################################\n",
    "    # GENERATING THE RESOLVED VELOCITY FIELD #\n",
    "    ##########################################\n",
    "\n",
    "    # To determine the resolved spectrum, we take the previous\n",
    "    # randomized copies of the velocity field and take the average\n",
    "    # along the number of copies to get an average velocity\n",
    "    # field.  Then add the systemic redshift\n",
    "    vrad_avg = cz + np.average(vrad_pop, axis=0)\n",
    "\n",
    "    # Store the density information (maybe for later use)\n",
    "    resolved_dens= np.average(dens_pop, axis=0)  # Sum up the densities\n",
    "\n",
    "    # Determine the values of the velocity contours to use.\n",
    "    # Set the velocity limits, rounded to velocity_step\n",
    "    min_resolved_vel = velocity_step*np.floor(np.min(vrad_avg)/velocity_step)\n",
    "    max_resolved_vel = velocity_step*np.ceil(np.max(vrad_avg)/velocity_step)\n",
    "    if (max_resolved_vel - min_resolved_vel)>200:\n",
    "        vel_levels = np.linspace(min_resolved_vel, max_resolved_vel, 11)\n",
    "    else: # Space every 20 km/s\n",
    "        vel_levels = np.arange(min_resolved_vel, max_resolved_vel, 20)\n",
    "\n",
    "    # Set the density limit rounded to 0.1 atom/cm^3\n",
    "    max_resolved_dens = dens_step*np.ceil(np.max(resolved_dens)/dens_step)\n",
    "    dens_levels = np.linspace((dens_step/2), max_resolved_dens, int(np.ceil(max_resolved_dens/(dens_step/2))) )\n",
    "\n",
    "    return(vrad_avg, vel_levels, resolved_dens, dens_levels)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "##\n",
    "## Functions to build plots in bqplot and matplotlib (including printable ones)\n",
    "##\n",
    "\n",
    "# Set up graph width and height\n",
    "graph_height = '250px'\n",
    "graph_width = '400px'\n",
    "\n",
    "\n",
    "def rotvel_plot(rad_spline, rotvel_spline):\n",
    "    global rotcurve_spline, rotcurve_handdraw\n",
    "    rotcurve_sc_x = bq.LinearScale()\n",
    "    rotcurve_sc_y = bq.LinearScale()\n",
    "    rotcurve_sc_y.max = maxvel\n",
    "    rotcurve_sc_y.min = 0\n",
    "    rotcurve_ax_x = bq.Axis(scale=rotcurve_sc_x, label='Galactocentric Radius (kpc)')\n",
    "    rotcurve_ax_y = bq.Axis(scale=rotcurve_sc_y, orientation='vertical', label='Rotational Velocity (km/s)')\n",
    "    rotcurve_ax_y.label_offset = '3.5em'\n",
    "    rotcurve_spline = bq.Lines(x=rad_spline, y=rotvel_spline, scales={'x': rotcurve_sc_x, 'y': rotcurve_sc_y},\n",
    "                             colors=['Blue'], labels=['cubic spline fit'], display_legend=False)\n",
    "    rotcurve_handdraw = bq.interacts.HandDraw(lines=rotcurve_spline)\n",
    "    \n",
    "    return (bq.Figure(marks=[rotcurve_spline], axes=[rotcurve_ax_x, rotcurve_ax_y], \n",
    "                      title='Galaxy Rotation Curve', legend_location='bottom', \n",
    "                      interaction=rotcurve_handdraw,\n",
    "                      layout=widgets.Layout(width=graph_width, height=graph_height, \n",
    "                                            fig_margin='0px 0px 0px 0px')))\n",
    "\n",
    "\n",
    "def densprof_plot(rad_spline, density_spline):\n",
    "    global densprof_spline, densprof_handdraw\n",
    "    densprof_sc_x = bq.LinearScale()\n",
    "    densprof_sc_y = bq.LinearScale()\n",
    "    densprof_sc_y.max = maxdens\n",
    "    densprof_sc_y.min = 0.0\n",
    "    densprof_ax_x = bq.Axis(scale=densprof_sc_x, label='Galactocentric Radius (kpc)')\n",
    "    densprof_ax_y = bq.Axis(scale=densprof_sc_y, orientation='vertical', label='HI Gas Density (atom/cm^3)')\n",
    "    densprof_ax_y.label_offset = '3.5em'\n",
    "    densprof_spline = bq.Lines(x=rad_spline, y=density_spline, scales={'x': densprof_sc_x, 'y': densprof_sc_y},\n",
    "                             colors=['Blue'], labels=['cubic spline fit'], display_legend=False)\n",
    "    densprof_handdraw = bq.interacts.HandDraw(lines=densprof_spline)\n",
    "    return(bq.Figure(marks=[densprof_spline], axes=[densprof_ax_x, densprof_ax_y], \n",
    "                     title='Galaxy Gas Density Profile', legend_location='top', \n",
    "                     interaction=densprof_handdraw,\n",
    "                     layout=widgets.Layout(width=graph_width, height=graph_height,\n",
    "                                           fig_margin='0px 0px 0px 0px')))\n",
    "\n",
    "\n",
    "def spectra_plots(x, y, vel, flux, vrad_map, vel_levels, dens_map, dens_levels):\n",
    "    # \n",
    "    # Generate plot of HI spectra (both 1D and resolved) as an image\n",
    "\n",
    "    # Clear any existing figure and create offline figure\n",
    "    _ = plt.clf()\n",
    "    _ = plt.ioff()\n",
    "\n",
    "    # Adjust plot defaults\n",
    "    mpl.rc('axes',titlesize='large')\n",
    "    mpl.rc('axes',labelsize='medium')\n",
    "    mpl.rc('xtick',labelsize='medium')\n",
    "    mpl.rc('ytick',labelsize='medium')\n",
    "\n",
    "    fig = plt.figure(1, facecolor='w', edgecolor='k')\n",
    "    _ = fig.set_size_inches(8, 4)\n",
    "\n",
    "    # Set up the plot grid\n",
    "    ax1 = plt.subplot2grid((1,2), (0,0))\n",
    "    ax2 = plt.subplot2grid((1,2), (0,1))\n",
    "\n",
    "    # Set up spacing being plots\n",
    "    plt.subplots_adjust(wspace=0.3)\n",
    "    plt.subplots_adjust(hspace=0.7)\n",
    "\n",
    "    #\n",
    "    # Plot the unresolved spectral profile\n",
    "    #\n",
    "    ax1.set_xlabel('Heliocentric Velocity (km/s)', fontsize='medium')\n",
    "    ax1.set_ylabel('Flux (Jy-km/s)', fontsize='medium')\n",
    "    unresolved_title = \"Unresolved HI Spectra (i: {:02d}$^\\circ$)\".format(int(incl))\n",
    "    ax1.set_title(unresolved_title, fontsize='large')\n",
    "    ax1.plot(vel, flux, 'b-')\n",
    "    \n",
    "    #\n",
    "    # Plot the resolved velocity field\n",
    "    #\n",
    "    ax2.set_xlabel('X (kpc)', fontsize='medium')\n",
    "    ax2.set_ylabel('Y (kpc)', fontsize='medium')\n",
    "    ax2.set_xlim(-HIradius, HIradius)\n",
    "    ax2.set_ylim(-HIradius, HIradius)\n",
    "    resolved_title = \"Resolved HI Velocity Field (i: {:02d}$^\\circ$)\".format(int(incl))\n",
    "    ax2.set_title(resolved_title, fontsize='large')\n",
    "\n",
    "    #imshow(resolved_vel, interpolation='bilinear',cmap=gray)\n",
    "    cm_name = 'Greys'\n",
    "    CSF = ax2.contourf(x, y, dens_map, levels=dens_levels, alpha=1, cmap = cm.get_cmap(cm_name), extend='max')\n",
    "    CS = ax2.contour(x, y, vrad_map, levels=vel_levels, extend='neither', alpha=0.75, colors='k')\n",
    "\n",
    "    # Label levels with specially formatted floats\n",
    "    CS.levels = [nf(val) for val in CS.levels ]\n",
    "    ax2.clabel(CS, CS.levels, inline=True, fmt='%r', fontsize=12)\n",
    "\n",
    "    # Display the plot to a PNG memory buffer\n",
    "    buf = io.BytesIO()\n",
    "    _ = fig.savefig(buf, format='png', dpi=150)\n",
    "    buf.seek(0)\n",
    "    pngimg = buf.read()\n",
    "    buf.close()\n",
    "    plt.clf()\n",
    "\n",
    "    return(pngimg)\n",
    "\n",
    "\n",
    "def make_printable_plot(b=None):\n",
    "    global imgfile_name\n",
    "    global rad_spline, rotvel_spline, density_spline\n",
    "    global x, y, vel, flux, vrad_map, vel_levels, dens_map, dens_levels\n",
    "\n",
    "    ######################\n",
    "    #  CREATE THE PLOTS  #\n",
    "    ###################### \n",
    "\n",
    "    # Adjust plot defaults\n",
    "    mpl.rc('axes',titlesize='medium')\n",
    "    mpl.rc('axes',labelsize='small')\n",
    "    mpl.rc('xtick',labelsize='small')\n",
    "    mpl.rc('ytick',labelsize='small')\n",
    "\n",
    "    fig = plt.figure(1, facecolor='w', edgecolor='k')\n",
    "    fig.set_size_inches(14,10)\n",
    "    plt.subplots_adjust(wspace=0.3)\n",
    "    plt.subplots_adjust(hspace=0.4)\n",
    "\n",
    "    #\n",
    "    # (1) Plot the rotation curve\n",
    "    #\n",
    "    plt.subplot(221) \n",
    "    plt.xlabel('Galactocentric Radius (kpc)')\n",
    "    plt.ylabel('Rotational Velocity (km/s)')\n",
    "    plt.title('Galaxy Rotation Curve')\n",
    "    plt.ylim(0, maxvel)\n",
    "    plt.plot(rad_spline, rotvel_spline, 'b-', label='spline fit')\n",
    "    \n",
    "    #\n",
    "    # (2) Plot the density profile\n",
    "    #\n",
    "    plt.subplot(222) \n",
    "    plt.xlabel('Galactocentric Radius (kpc)')\n",
    "    plt.ylabel('HI Gas Density (atom/cm^3)')\n",
    "    plt.title('Galaxy Gas Density Profile')\n",
    "    plt.ylim(0, maxdens)\n",
    "    plt.plot(rad_spline, density_spline, 'b-', label='spline fit')\n",
    "    \n",
    "    #\n",
    "    # (3) Plot the unresolved spectral profile\n",
    "    #\n",
    "    plt.subplot(223) \n",
    "    plt.xlabel('Heliocentric Velocity (km/s)')\n",
    "    plt.ylabel('Flux (Jy-km/s)')\n",
    "    unresolved_title = \"Unresolved HI Spectra (i: {:02d}$^\\circ$)\".format(int(incl))\n",
    "    plt.title(unresolved_title)\n",
    "    plt.plot(vel, flux, 'b-')\n",
    "\n",
    "    #\n",
    "    # (4) Plot the resolved velocity field\n",
    "    #\n",
    "    plt.subplot(224) \n",
    "    plt.xlabel('X (kpc)')\n",
    "    plt.ylabel('Y (kpc)')\n",
    "    plt.xlim(-HIradius, HIradius)\n",
    "    plt.ylim(-HIradius, HIradius)\n",
    "    resolved_title = \"Resolved HI Velocity Field (i: {:02d}$^\\circ$)\".format(int(incl))\n",
    "    plt.title(resolved_title)\n",
    "\n",
    "    #imshow(resolved_vel, interpolation='bilinear',cmap=gray)\n",
    "    cm_name = 'Greys'\n",
    "    CSF = plt.contourf(x, y, dens_map, levels=dens_levels, alpha=1, cmap = cm.get_cmap(cm_name), extend='max')\n",
    "    CS = plt.contour(x, y, vrad_map, levels=vel_levels, extend='neither', alpha=0.75, colors='k')\n",
    "\n",
    "    # Label levels with specially formatted floats\n",
    "    CS.levels = [nf(val) for val in CS.levels ]\n",
    "    plt.clabel(CS, CS.levels, inline=True, fmt='%r', fontsize=12)\n",
    "\n",
    "    # Save this image as a figure too\n",
    "    plt.savefig(imgfile_name.value, dpi=300)\n",
    "    \n",
    "    # Create the image file and create a link if the filename is legitimate\n",
    "    filetypes = [\"jpg\", \"png\", \"pdf\"]\n",
    "    img_type = imgfile_name.value[-3:]\n",
    "    if img_type in filetypes:\n",
    "        # Store images in a temporary directory\n",
    "        true_loc = \"{0}/{1}\".format(image_directory,imgfile_name.value)\n",
    "        plt.savefig(true_loc, dpi=300)\n",
    "        savelink.value=\"<a href=\\\"{}\\\">Click to download!</a>\".format(\"DisplayRecentImage.ipynb\")\n",
    "    else:\n",
    "        savelink.value=\"<b>Illegal image type ({})!</b>\".format(img_type)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#\n",
    "# Functions for handling interactive elements\n",
    "#\n",
    "\n",
    "def replot(b=None):\n",
    "    global HIradius, dx, dr_spline, cz, velocity_step, incl\n",
    "    global rad_raw, rotvel_raw, density_raw\n",
    "    global rad_spline, rotvel_spline, density_spline\n",
    "    global x, y, vel, flux, vrad_map, vel_levels, dens_map, dens_levels\n",
    "    \n",
    "    # Update the spectra plot information\n",
    "    x, y, vrad_pop, dens_pop = build_galaxy_velcube(rotvel_spline, density_spline, HIradius, dx, dr_spline)\n",
    "    vel, flux = unresolved_spectra(vrad_pop, dens_pop, cz)\n",
    "    vrad_map, vel_levels, dens_map, dens_levels = resolved_spectra(cz, vrad_pop, dens_pop, velocity_step)\n",
    "\n",
    "    # Update spectra plots\n",
    "    spectra.value = spectra_plots(x, y, vel, flux, vrad_map, vel_levels, dens_map, dens_levels)\n",
    "    \n",
    "\n",
    "def cz_changed(change):\n",
    "    global cz, H_0, HIradius\n",
    "    \n",
    "    # Update the cz value\n",
    "    cz = change.new\n",
    "    \n",
    "    # Check the beam size\n",
    "    check_beam(cz, H_0, HIradius)\n",
    "    \n",
    "    # Replot\n",
    "    replot()    \n",
    "\n",
    "def incl_changed(change):\n",
    "    global incl\n",
    "    \n",
    "    # Update the incl value\n",
    "    incl = change.new\n",
    "    \n",
    "    # Replot\n",
    "    replot()  \n",
    "\n",
    "def sigma_changed(change):\n",
    "    global vel_sigma\n",
    "    \n",
    "    # Update the sigma value\n",
    "    vel_sigma = change.new\n",
    "    \n",
    "    # Replot\n",
    "    replot()    \n",
    "    \n",
    "def update_rotcurve(change):\n",
    "    global rotcurve_spline, rad_spline, rotvel_spline, density_spline\n",
    "    \n",
    "    # Force any negative points to zero\n",
    "    # (there doesn't seem to be a way to update handdrawn line)\n",
    "    rotcurve_spline.y[rotcurve_spline.y<0] = 0\n",
    "    \n",
    "    # Update spline based on handdraw action on rotcurve_spline\n",
    "    rotvel_spline = rotcurve_spline.y\n",
    "    \n",
    "    # Replot\n",
    "    replot()\n",
    "\n",
    "def update_densprof(change):\n",
    "    global densprof_spline, rad_spline, rotvel_spline, density_spline\n",
    "\n",
    "    # Force any negative points to zero\n",
    "    # (there doesn't seem to be a way to update handdrawn line)\n",
    "    densprof_spline.y[densprof_spline.y<0] = 0\n",
    "\n",
    "    # Update spline based on handdraw action on densprof_spline\n",
    "    density_spline = densprof_spline.y\n",
    "    \n",
    "    # Replot\n",
    "    replot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "hide-output"
    ]
   },
   "outputs": [],
   "source": [
    "%%capture --no-stderr --no-stdout \n",
    "\n",
    "##########################################\n",
    "#  GET INITIAL SPECTRA AND VELOCITY MAP  #\n",
    "##########################################\n",
    "\n",
    "# Generate a population of radial velocity and density maps\n",
    "x, y, vrad_pop, dens_pop = build_galaxy_velcube(rotvel_spline, density_spline, HIradius, dx, dr_spline)\n",
    "\n",
    "# Compute the unresolved 1-D spectra\n",
    "vel, flux = unresolved_spectra(vrad_pop, dens_pop, cz)\n",
    "\n",
    "# Compute 2-D velocity map\n",
    "vrad_map, vel_levels, dens_map, dens_levels = resolved_spectra(cz, vrad_pop, dens_pop, velocity_step)\n",
    "\n",
    "# Generate plots\n",
    "# (1) Generate plot of plot of the rotation curve of the Galaxy\n",
    "# (2) Draw the density profile of the Galaxy\n",
    "# (3) Using matplotlib to generate the unresolved and resolved spectra plots\n",
    "#\n",
    "rot_curve = rotvel_plot(rad_spline, rotvel_spline)\n",
    "density_profile = densprof_plot(rad_spline, density_spline)\n",
    "\n",
    "#\n",
    "# Generate PNG to show\n",
    "#\n",
    "spectra2show = spectra_plots(x, y, vel, flux, vrad_map, vel_levels, dens_map, dens_levels)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "lines_to_next_cell": 0
   },
   "outputs": [],
   "source": [
    "#\n",
    "# Create Control Sliders/Fields/Buttons\n",
    "#\n",
    "sigma_label = widgets.Label(value='$\\sigma_V$ (km/s): ')\n",
    "sigma_slider = widgets.IntSlider(min=0, max=50, step=1, value=vel_sigma,\n",
    "                                 continuous_update=False, orientation='horizontal',\n",
    "                                 readout=True, readout_format='02d',\n",
    "                                 layout=widgets.Layout(height='40px', min_height='40px', max_height='100px',\n",
    "                                                       min_width='50px', width='200px', max_width='300px',\n",
    "                                                       overflow='hidden'))\n",
    "sigma_control = widgets.HBox([sigma_label, sigma_slider], \n",
    "                             layout=widgets.Layout(align_content='center', align_items='center'))\n",
    "\n",
    "incl_label = widgets.Label(value='Inclination (deg): ')\n",
    "incl_slider = widgets.IntSlider(min=0, max=80, step=1, value=incl,\n",
    "                                 continuous_update=False, orientation='horizontal',\n",
    "                                 readout=True, readout_format='04d',\n",
    "                                 layout=widgets.Layout(height='40px', min_height='40px', max_height='100px',\n",
    "                                                       min_width='50px', width='200px', max_width='300px',\n",
    "                                                       overflow='hidden'))\n",
    "incl_control = widgets.HBox([incl_label, incl_slider], \n",
    "                             layout=widgets.Layout(align_content='center', align_items='center'))\n",
    "\n",
    "cz_label = widgets.Label(value='$cz$ (km/s): ')\n",
    "cz_slider = widgets.IntSlider(min=0, max=20000, step=10, value=cz,\n",
    "                                 continuous_update=False, orientation='horizontal',\n",
    "                                 readout=True, readout_format='04d',\n",
    "                                 layout=widgets.Layout(height='40px', min_height='40px', max_height='100px',\n",
    "                                                       min_width='50px', width='200px', max_width='300px',\n",
    "                                                       overflow='hidden'))\n",
    "cz_control = widgets.HBox([cz_label, cz_slider], \n",
    "                             layout=widgets.Layout(align_content='center', align_items='center'))\n",
    "cz_notice = widgets.HTML(value='', layout=widgets.Layout(height='40px', min_height='40px', max_height='100px',\n",
    "                                                       min_width='50px', width='200px', max_width='300px'))\n",
    "\n",
    "imgfile_label = widgets.HTML(value='Image Filename: ', \n",
    "                             layout=widgets.Layout(min_width='70px', width='100px', max_width='200px'))\n",
    "imgfile_name = widgets.Text(value=image_filename, \n",
    "                             layout=widgets.Layout(min_width='100px', width='200px', max_width='300px'))\n",
    "imgfile_control = widgets.HBox([imgfile_label, imgfile_name], \n",
    "                               layout=widgets.Layout(align_content='center', align_items='center'))\n",
    "\n",
    "savefile_button = widgets.Button(description='Save Plots', disabled=False,\n",
    "                                 tooltip='Click to save image file',\n",
    "                                 layout=widgets.Layout(min_width='100px', width='150px', max_width='200px'))\n",
    "\n",
    "savelink = widgets.HTML(value='')\n",
    "\n",
    "#\n",
    "# Generate image widget to display spectra\n",
    "#\n",
    "spectra = widgets.Image(value=spectra2show, format='png', width=800, height=400)\n",
    "\n",
    "#\n",
    "# Display the plots\n",
    "#\n",
    "\n",
    "settings_top = widgets.HBox([sigma_control, imgfile_control, savefile_button, savelink], \n",
    "                               layout=widgets.Layout(align_content='center', align_items='center',\n",
    "                                                     margin='0px', width='850px', \n",
    "                                                     overflow='hidden'))\n",
    "settings_middle = widgets.HBox([cz_control, cz_notice, incl_control], \n",
    "                               layout=widgets.Layout(align_content='center', align_items='center',\n",
    "                                                     margin='0px', width='850px', \n",
    "                                                     overflow='hidden'))\n",
    "\n",
    "settings_bottom = widgets.HBox([rot_curve, density_profile], \n",
    "                            layout=widgets.Layout(margin='0px', width='850px', \n",
    "                                                  overflow='hidden'))\n",
    "settings = widgets.VBox([settings_top, settings_middle, settings_bottom], \n",
    "                        layout=widgets.Layout(align_content='center', align_items='center',\n",
    "                                              margin='0px', width='850px', \n",
    "                                              overflow='hidden'))\n",
    "MainDisplay = widgets.VBox([settings, spectra], \n",
    "                           layout=widgets.Layout(align_content='center', align_items='center',\n",
    "                                                 margin='0px', width='850px', \n",
    "                                                 overflow='hidden'))\n",
    "display(MainDisplay)\n",
    "\n",
    "#\n",
    "# Direct functions to be called if various widgets change\n",
    "#\n",
    "sigma_slider.observe(sigma_changed, 'value')\n",
    "cz_slider.observe(cz_changed, 'value')\n",
    "incl_slider.observe(incl_changed, 'value')\n",
    "rotcurve_spline.observe(update_rotcurve, names=['y'])\n",
    "densprof_spline.observe(update_densprof, names=['y'])\n",
    "savefile_button.on_click(make_printable_plot)\n",
    "\n",
    "# Still need to add Warning about beam size"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}