{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "# Determining and plotting the altitude/azimuth of a celestial object\n",
        "\n",
        "This notebook is from astropy examples repo\n",
        "\n",
        "This example demonstrates coordinate transformations and the creation of\n",
        "visibility curves to assist with observing run planning.\n",
        "\n",
        "In this example, we make a `~astropy.coordinates.SkyCoord` instance for M33.\n",
        "The altitude-azimuth coordinates are then found using\n",
        "`astropy.coordinates.EarthLocation` and `astropy.time.Time` objects.\n",
        "\n",
        "This example is meant to demonstrate the capabilities of the\n",
        "`astropy.coordinates` package. For more convenient and/or complex observation\n",
        "planning, consider the [astroplan](https://astroplan.readthedocs.org/)\n",
        "package.\n",
        "\n",
        "\n",
        "*By: Erik Tollerud, Kelle Cruz*\n",
        "\n",
        "*License: BSD*\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Let's suppose you are planning to visit picturesque Bear Mountain State Park\n",
        "in New York, USA. You're bringing your telescope with you (of course), and\n",
        "someone told you M33 is a great target to observe there. You happen to know\n",
        "you're free at 11:00 pm local time, and you want to know if it will be up.\n",
        "Astropy can answer that.\n",
        "\n",
        "Import numpy and matplotlib. For the latter, use a nicer set of plot\n",
        "parameters and set up support for plotting/converting quantities.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 2,
      "metadata": {
        "collapsed": false
      },
      "outputs": [
        {
          "data": {
            "text/plain": [
              "<astropy.visualization.units.quantity_support.<locals>.MplQuantityConverter at 0x7fc004e85d20>"
            ]
          },
          "execution_count": 2,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "import matplotlib.pyplot as plt\n",
        "import numpy as np\n",
        "\n",
        "from astropy.visualization import astropy_mpl_style, quantity_support\n",
        "\n",
        "plt.style.use(astropy_mpl_style)\n",
        "quantity_support()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Import the packages necessary for finding coordinates and making\n",
        "coordinate transformations\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 3,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import astropy.units as u\n",
        "from astropy.coordinates import AltAz, EarthLocation, SkyCoord\n",
        "from astropy.time import Time"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "`astropy.coordinates.SkyCoord.from_name` uses Simbad to resolve object\n",
        "names and retrieve coordinates.\n",
        "\n",
        "Get the coordinates of M33:\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 4,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "m33 = SkyCoord.from_name('M31')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Use `astropy.coordinates.EarthLocation` to provide the location of Bear\n",
        "Mountain and set the time to 11pm EDT on 2012 July 12:\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "bear_mountain = EarthLocation(lat=41.3*u.deg, lon=-74*u.deg, height=390*u.m)\n",
        "utcoffset = -4*u.hour  # Eastern Daylight Time\n",
        "time = Time('2012-7-12 23:00:00') - utcoffset"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "`astropy.coordinates.EarthLocation.get_site_names` and\n",
        "`~astropy.coordinates.EarthLocation.get_site_names` can be used to get\n",
        "locations of major observatories.\n",
        "\n",
        "Use `astropy.coordinates` to find the Alt, Az coordinates of M33 at as\n",
        "observed from Bear Mountain at 11pm on 2012 July 12.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "m33altaz = m33.transform_to(AltAz(obstime=time,location=bear_mountain))\n",
        "print(f\"M33's Altitude = {m33altaz.alt:.2}\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "This is helpful since it turns out M33 is barely above the horizon at this\n",
        "time. It's more informative to find M33's airmass over the course of\n",
        "the night.\n",
        "\n",
        "Find the alt,az coordinates of M33 at 100 times evenly spaced between 10pm\n",
        "and 7am EDT:\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "midnight = Time('2012-7-13 00:00:00') - utcoffset\n",
        "delta_midnight = np.linspace(-2, 10, 100)*u.hour\n",
        "frame_July13night = AltAz(obstime=midnight+delta_midnight,\n",
        "                          location=bear_mountain)\n",
        "m33altazs_July13night = m33.transform_to(frame_July13night)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "convert alt, az to airmass with `~astropy.coordinates.AltAz.secz` attribute:\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "m33airmasss_July13night = m33altazs_July13night.secz"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Plot the airmass as a function of time:\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "plt.plot(delta_midnight, m33airmasss_July13night)\n",
        "plt.xlim(-2, 10)\n",
        "plt.ylim(1, 4)\n",
        "plt.xlabel('Hours from EDT Midnight')\n",
        "plt.ylabel('Airmass [Sec(z)]')\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Use  `~astropy.coordinates.get_sun` to find the location of the Sun at 1000\n",
        "evenly spaced times between noon on July 12 and noon on July 13:\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "from astropy.coordinates import get_sun\n",
        "\n",
        "delta_midnight = np.linspace(-12, 12, 1000)*u.hour\n",
        "times_July12_to_13 = midnight + delta_midnight\n",
        "frame_July12_to_13 = AltAz(obstime=times_July12_to_13, location=bear_mountain)\n",
        "sunaltazs_July12_to_13 = get_sun(times_July12_to_13).transform_to(frame_July12_to_13)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Do the same with `~astropy.coordinates.get_body` to find when the moon is\n",
        "up. Be aware that this will need to download a 10MB file from the internet\n",
        "to get a precise location of the moon.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "from astropy.coordinates import get_body\n",
        "\n",
        "moon_July12_to_13 = get_body(\"moon\", times_July12_to_13)\n",
        "moonaltazs_July12_to_13 = moon_July12_to_13.transform_to(frame_July12_to_13)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Find the alt,az coordinates of M33 at those same times:\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "m33altazs_July12_to_13 = m33.transform_to(frame_July12_to_13)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Make a beautiful figure illustrating nighttime and the altitudes of M33 and\n",
        "the Sun over that time:\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "plt.plot(delta_midnight, sunaltazs_July12_to_13.alt, color='r', label='Sun')\n",
        "plt.plot(delta_midnight, moonaltazs_July12_to_13.alt, color=[0.75]*3, ls='--', label='Moon')\n",
        "plt.scatter(delta_midnight, m33altazs_July12_to_13.alt,\n",
        "            c=m33altazs_July12_to_13.az.value, label='M33', lw=0, s=8,\n",
        "            cmap='viridis')\n",
        "plt.fill_between(delta_midnight, 0*u.deg, 90*u.deg,\n",
        "                 sunaltazs_July12_to_13.alt < -0*u.deg, color='0.5', zorder=0)\n",
        "plt.fill_between(delta_midnight, 0*u.deg, 90*u.deg,\n",
        "                 sunaltazs_July12_to_13.alt < -18*u.deg, color='k', zorder=0)\n",
        "plt.colorbar().set_label('Azimuth [deg]')\n",
        "plt.legend(loc='upper left')\n",
        "plt.xlim(-12*u.hour, 12*u.hour)\n",
        "plt.xticks((np.arange(13)*2-12)*u.hour)\n",
        "plt.ylim(0*u.deg, 90*u.deg)\n",
        "plt.xlabel('Hours from EDT Midnight')\n",
        "plt.ylabel('Altitude [deg]')\n",
        "plt.show()"
      ]
    }
  ],
  "metadata": {
    "kernelspec": {
      "display_name": "Python 3",
      "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.10.12"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}