{
  "cells": [
    {
      "cell_type": "markdown",
      "source": [
        "# Radio Telescope - Example of Loading Spectrum Data"
      ],
      "metadata": {
        "id": "TA3dSyKj1019"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "## Load the Data from a .rad Text File"
      ],
      "metadata": {
        "id": "d0U8GYPu18LH"
      }
    },
    {
      "cell_type": "code",
      "execution_count": 1,
      "metadata": {
        "id": "NEmH7StEx4bT"
      },
      "outputs": [],
      "source": [
        "import numpy as np\n",
        "from scipy import stats\n",
        "import matplotlib.pyplot as plt"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "SS5gGl5Rx4bU"
      },
      "outputs": [],
      "source": [
        "#%matplotlib notebook"
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "Before you run this code, open the file 'Test.rad' using a text editor.  Make sure you are familiar with what is in the file: (more explanation in that file)\n",
        "* rows: different times\n",
        "* columns (after the first few): different frequencies"
      ],
      "metadata": {
        "id": "1leUFy7oyyz2"
      }
    },
    {
      "cell_type": "code",
      "execution_count": 27,
      "metadata": {
        "scrolled": true,
        "id": "tG5EVa4Ax4bV"
      },
      "outputs": [],
      "source": [
        "file=\"Test.rad\"                                              # edit in file name\n",
        "time= np.loadtxt(file, dtype='str', usecols=0, comments='*')     # reads time stamp as a string (for now)\n",
        "freq0= np.loadtxt(file, dtype='float', usecols=5, comments='*')       # reads the frequency of the first element in the spectrum\n",
        "freqstep= np.loadtxt(file, dtype='float', usecols=6, comments='*')     # reads the frequency step\n",
        "nspec = np.loadtxt(file, dtype='int',usecols=8, comments='*')          # reads the number of elements in each spectrum\n",
        "endspec=9+nspec[0]\n",
        "mylist2=np.arange(9,endspec,1,dtype=int)            # spectrum is written starting with column 9 and extending for nspec elements\n",
        "data= np.loadtxt(file, dtype='float', usecols=mylist2, comments='*')   # data is a two dimensional array\n",
        "vlsr= np.loadtxt(file, dtype='float', usecols=endspec+1, comments='*') # VLSR (velocity of local standard of rest)\n"
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "Freq0 and freqstep should ideally be the same value at every time step - if you change frequency settings, you should save and start a new data file.  Next step is to generate an array containing the list of the frequencies of all the channels (the x-axis for your spectrum plots)."
      ],
      "metadata": {
        "id": "vEDsT0LjKaY2"
      }
    },
    {
      "cell_type": "code",
      "execution_count": 3,
      "metadata": {
        "id": "TRG7_Eawx4bV"
      },
      "outputs": [],
      "source": [
        "freqval = np.linspace(freq[0],freq[0]+(nspec[0]-1)*freqstep[0],nspec[0])"
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "Add a code cell below.  Check the shape of the 2D data array using data.shape.  Then check the length of the time array and the length of the freqval array.  Out of the rows and columns in data, which direction matches up to frequency, and which to time?  Document your conclusions here."
      ],
      "metadata": {
        "id": "wiULqJwt1RYH"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "## Plot the Data"
      ],
      "metadata": {
        "id": "WaJKfbsp2B_I"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "First plot the data from a single timestamp."
      ],
      "metadata": {
        "id": "kDiyiG9Z_xrR"
      }
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "scrolled": false,
        "id": "OFpO8eVGx4bV"
      },
      "outputs": [],
      "source": [
        "plt.figure()\n",
        "ind = 7       # choosing an arbitrary timestamp for which to plot the spectrum\n",
        "#data = data - 1550     # subtracting off the \"plateau\" emission just to make the plot look nice; not necessary\n",
        "plt.ylim((1500,2300))          # setting up y-limits to make the plot look nice; not necessary\n",
        "plt.plot(freqval,data[ind])\n",
        "plt.title(\"Radio Telescope spectrum\")\n",
        "plt.xlabel('Frequency (MHz)')\n",
        "plt.ylabel(\"Intensity (arbitrary units)\")"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "CmFrnWc-x4bW"
      },
      "outputs": [],
      "source": [
        "avgspec = np.mean(data,axis=0)  # average all the rows together in the data array (why?)\n",
        "plt.figure()\n",
        "plt.ylim((1500,2300))          # setting up y-limits to make the plot look nice; not necessary\n",
        "plt.plot(freqval,avgspec)\n",
        "plt.title(\"Radio Telescope spectrum\")\n",
        "plt.xlabel('Frequency (MHz)')\n",
        "plt.ylabel(\"Intensity (arbitrary units)\")"
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "What is the difference between the first and second plot above?  Why are they different?  What does the first line of code above the second plot do? (the avgspec line)"
      ],
      "metadata": {
        "id": "RWVePi8-LIFN"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "## Measuring $v_\\textrm{rot}$\n",
        "\n",
        "To measure the galactic rotation speed $v_\\textrm{rot}$, you need to estimate the *maximum* frequency at which there is emission above the baseline level.  Carry out the following steps to get an idea for how you will measure the rotation speed of the galaxy.  You'll eventually collect data (like the above spectrum graph) pointing in many different directions, to make a graph of rotation speed vs. distance from the Galactic Center.\n",
        "\n",
        "1) Estimate this maximum frequency by eye (and start thinking about other more quantitative ways you could measure it later).  What did you get?  This is $f_{obs}$.\n",
        "\n",
        "2) Use $f_{obs}$ to calculate $v_{obs}$.\n",
        "\n",
        "3) Following the explanation in the intro packet, convert this $v_{obs}$ to $v_{los}$.  Since this is the maximum value for $v_{los}$ at which there is emission, we call this the tangential velocity $v_\\textrm{rot}$."
      ],
      "metadata": {
        "id": "SKxfjjCkLy3b"
      }
    },
    {
      "cell_type": "code",
      "source": [],
      "metadata": {
        "id": "1cxC29wtNgOk"
      },
      "execution_count": null,
      "outputs": []
    }
  ],
  "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.9.13"
    },
    "colab": {
      "provenance": []
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}