test_ephem.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Astropy implementation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 188,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from astropy.coordinates import solar_system_ephemeris, EarthLocation, get_body\n",
    "from astropy.time import Time\n",
    "from numpy import pi as PI\n",
    "from numpy import arccos, cos, sin, sqrt\n",
    "from astropy.coordinates import AltAz\n",
    "\n",
    "MIN_VENUS_SEP = 2.0 * PI / 180.0\n",
    "MIN_MARS_SEP = 2.0 * PI / 180.0\n",
    "MIN_JUPITER_SEP = 2.0 * PI / 180.0\n",
    "MIN_SATURN_SEP = 2.0 * PI / 180.0\n",
    "MIN_NEPTUNE_SEP = 2.0 * PI / 180.0\n",
    "MIN_URANUS_SEP = 2.0 * PI / 180.0\n",
    "MIN_CERES_SEP = 2.0 * PI / 180.0\n",
    "def angdist(ra1,dec1,ra2,dec2):\n",
    "    \"\"\"\n",
    "    Computes the angular distance between two objects\n",
    "\n",
    "    Args:\n",
    "        ra1: float (R.A for one of the objects, degrees)\n",
    "        ra2: float (R.A. for the other object, degrees)\n",
    "        dec1: float (DEC for one of the objects, degrees)\n",
    "        dec2: float (DEC for the other object, degrees)\n",
    "\n",
    "    Returns:\n",
    "        float, Angular distance in radians between the objects\n",
    "    \"\"\"\n",
    "    angdist = arccos(sin(dec1*180/PI)*sin(dec2*180/PI)+ \\\n",
    "                     cos(dec1*180/PI)*cos(dec2*180/PI)*cos((ra1-ra2)*180/PI))\n",
    "    return angdist\n",
    "\n",
    "\n",
    "def avoidObject(datetime, ra0, dec0):\n",
    "    \"\"\"\n",
    "    Checks whether all the objects on the list are far enough away from\n",
    "    the input coordinates.\n",
    "    The current list has: Venus, Mars, Jupiter, Saturn, Neptune, Uranus;\n",
    "    the Moon is treated separately.\n",
    "\n",
    "    Args:\n",
    "        datetime: astropy Time; should have timezone info\n",
    "        ra0: float (apparent or observed, degrees)\n",
    "        dec0: float (apparent or observed, degrees)\n",
    "\n",
    "    Returns:\n",
    "        bool, True if all objects on the list are far enough away\n",
    "    \"\"\"\n",
    "\n",
    "    #KPNO's location (slightly different than in kpno.py)\n",
    "    location = EarthLocation.of_site('Kitt Peak National Observatory')\n",
    "    with solar_system_ephemeris.set('de432s'):\n",
    "        venus = get_body('venus', datetime, location)\n",
    "        mars = get_body('mars', datetime, location)\n",
    "        jupiter = get_body('jupiter', datetime, location)\n",
    "        saturn = get_body('saturn', datetime, location)\n",
    "        uranus = get_body('uranus', datetime, location)\n",
    "        neptune = get_body('neptune', datetime, location)\n",
    "    \n",
    "    dist_venus = angdist(venus.ra.deg,venus.dec.deg,ra0,dec0)\n",
    "    if dist_venus < MIN_VENUS_SEP:\n",
    "        return False\n",
    "    dist_mars = angdist(mars.ra.deg,mars.dec.deg,ra0,dec0)\n",
    "    if dist_mars < MIN_MARS_SEP:\n",
    "        return False\n",
    "    dist_jupiter = angdist(jupiter.ra.deg,jupiter.dec.deg,ra0,dec0)\n",
    "    if dist_jupiter < MIN_JUPITER_SEP:\n",
    "        return False\n",
    "    dist_saturn = angdist(saturn.ra.deg,saturn.dec.deg,ra0,dec0)\n",
    "    if dist_saturn < MIN_SATURN_SEP:\n",
    "        return False\n",
    "    dist_neptune = angdist(neptune.ra.deg,neptune.dec.deg,ra0,dec0)\n",
    "    if dist_neptune < MIN_NEPTUNE_SEP:\n",
    "        return False\n",
    "    dist_uranus = angdist(uranus.ra.deg,uranus.dec.deg,ra0,dec0)\n",
    "    if dist_uranus < MIN_URANUS_SEP:\n",
    "        return False\n",
    "\n",
    "    # If still here, return True\n",
    "    return True\n",
    "\n",
    "def moonLoc (datetime, ra0, dec0):\n",
    "    \"\"\"\n",
    "    Returns the distance to the Moon if RA and DEC as well as alt, az.\n",
    "\n",
    "    Args:\n",
    "        datetime: astropy Time; should have timezone info\n",
    "        ra0: float (apparent or observed, degrees)\n",
    "        dec0: float (apparent or observed, degrees)\n",
    "\n",
    "    Returns:\n",
    "        float, distance from the Moon (degrees)\n",
    "        float, Moon altitude (degrees)\n",
    "        float, Moon azimuth (degrees)\n",
    "    \"\"\"\n",
    "    \n",
    "    location = EarthLocation.of_site('Kitt Peak National Observatory')\n",
    "    aa = AltAz(location=location, obstime=datetime)\n",
    "    with solar_system_ephemeris.set('de432s'):\n",
    "        moon = get_body('moon', datetime, location)\n",
    "    moondist = angdist(moon.ra.deg,moon.dec.deg,ra0,dec0)\n",
    "    moon = moon.transform_to(aa)\n",
    "    return moondist*180.0/PI, moon.alt.deg, moon.az.deg\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 189,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "times = Time(\"2020-09-22 23:22\")\n",
    "time = Time(times, scale='utc', format='isot')\n",
    "ra = 30.0\n",
    "dec = 20.0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 190,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 4 µs, sys: 1 µs, total: 5 µs\n",
      "Wall time: 9.06 µs\n",
      "Separation angle 134.003263289\n",
      "Moon Altitude  34.3016855459\n",
      "Moon Azimuth 164.148457207\n"
     ]
    }
   ],
   "source": [
    "%time\n",
    "mloc = moonLoc(time,ra,dec)\n",
    "print 'Separation angle', mloc[0]\n",
    "print 'Moon Altitude ', mloc[1]\n",
    "print 'Moon Azimuth', mloc[2]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### PyEphem implementation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 191,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import ephem\n",
    "from datetime import datetime\n",
    "from numpy import pi as PI\n",
    "from surveysim.kpno import mayall\n",
    "\n",
    "MIN_VENUS_SEP = 2.0 * PI / 180.0\n",
    "MIN_MARS_SEP = 2.0 * PI / 180.0\n",
    "MIN_JUPITER_SEP = 2.0 * PI / 180.0\n",
    "MIN_SATURN_SEP = 2.0 * PI / 180.0\n",
    "MIN_NEPTUNE_SEP = 2.0 * PI / 180.0\n",
    "MIN_URANUS_SEP = 2.0 * PI / 180.0\n",
    "MIN_CERES_SEP = 2.0 * PI / 180.0\n",
    "\n",
    "def avoidObject(datetime, ra0, dec0):\n",
    "    \"\"\"\n",
    "    Checks whether all the objects on the list are far enough away from\n",
    "    the input coordinates.\n",
    "    The current list has: Venus, Mars, Jupiter, Saturn, Neptune, Uranus;\n",
    "    the Moon is treated separately.\n",
    "    Args:\n",
    "        datetime: datetime object; should have timezone info\n",
    "        ra0: float (apparent or observed, degrees)\n",
    "        dec0: float (apparent or observed, degrees)\n",
    "    Returns:\n",
    "        bool, True if all objects on the list are far enough away\n",
    "    \"\"\"\n",
    "\n",
    "    ra = PI * ra0/180.0\n",
    "    dec = PI * dec0/180.0\n",
    "    \n",
    "    dt = ephem.Date(datetime)\n",
    "    gatech = ephem.Observer()\n",
    "    gatech.lon, gatech.lat = mayall.west_lon_deg, mayall.lat_deg\n",
    "    gatech.date = dt\n",
    "    gatech.epoch = dt\n",
    "\n",
    "    venus = ephem.Venus()\n",
    "    venus.compute(gatech)\n",
    "    if ephem.separation(venus, (ra, dec)) < MIN_VENUS_SEP:\n",
    "        return False\n",
    "    mars = ephem.Mars()\n",
    "    mars.compute(gatech)\n",
    "    if ephem.separation(mars, (ra, dec)) < MIN_MARS_SEP:\n",
    "        return False\n",
    "    #ceres = ephem.Ceres()\n",
    "    #ceres.compute(gatech)\n",
    "    #if ephem.separation(ceres, (ra, dec)) < MIN_CERES_SEP:\n",
    "    #    return False\n",
    "    jupiter = ephem.Jupiter()\n",
    "    jupiter.compute(gatech)\n",
    "    if ephem.separation(jupiter, (ra, dec)) < MIN_JUPITER_SEP:\n",
    "        return False\n",
    "    saturn = ephem.Saturn()\n",
    "    saturn.compute(gatech)\n",
    "    if ephem.separation(saturn, (ra, dec)) < MIN_SATURN_SEP:\n",
    "        return False\n",
    "    neptune = ephem.Neptune()\n",
    "    neptune.compute(gatech)\n",
    "    if ephem.separation(neptune, (ra, dec)) < MIN_NEPTUNE_SEP:\n",
    "        return False\n",
    "    uranus = ephem.Uranus()\n",
    "    uranus.compute(gatech)\n",
    "    if ephem.separation(uranus, (ra, dec)) < MIN_URANUS_SEP:\n",
    "        return False\n",
    "\n",
    "    # If still here, return True\n",
    "    return True\n",
    "\n",
    "def moonLoc (datetime, ra0, dec0):\n",
    "    \"\"\"\n",
    "    Returns the distance to the Moon if RA and DEC as well as alt, az.\n",
    "    Args:\n",
    "        datetime: datetime object; should have timezone info\n",
    "        ra0: float (apparent or observed, degrees)\n",
    "        dec0: float (apparent or observed, degrees)\n",
    "    Returns:\n",
    "        float, distance from the Moon (degrees)\n",
    "        float, Moon altitude (degrees)\n",
    "        float, Moon azimuth (degrees) \n",
    "    \"\"\"\n",
    "\n",
    "    dt = ephem.Date(datetime)\n",
    "    gatech = ephem.Observer()\n",
    "    gatech.lon, gatech.lat = mayall.west_lon_deg, mayall.lat_deg\n",
    "    gatech.date = dt\n",
    "    gatech.epoch = dt\n",
    "\n",
    "    moon = ephem.Moon()\n",
    "    moon.compute(gatech)\n",
    "    ra = PI * ra0/180.0\n",
    "    dec = PI * dec0/180.0\n",
    "    moondist = ephem.separation(moon, (ra, dec))\n",
    "    moondist*180.0/PI\n",
    "    \n",
    "    return moondist*180.0/PI, (moon.alt)*180.0/PI, (moon.az)*180.0/PI"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 192,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2020-09-22 23:22:00\n"
     ]
    }
   ],
   "source": [
    "time = time.datetime\n",
    "print time"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 193,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Separation angle 137.81981246\n",
      "Moon Altitude  -80.3756444167\n",
      "Moon Azimuth 18.4890676587\n",
      "CPU times: user 188 µs, sys: 25 µs, total: 213 µs\n",
      "Wall time: 197 µs\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "mloc = moonLoc(time, ra, dec)\n",
    "print 'Separation angle', mloc[0]\n",
    "print 'Moon Altitude ', mloc[1]\n",
    "print 'Moon Azimuth', mloc[2]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The results don't agree"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Implementation using tables"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 194,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from astropy.table import Table\n",
    "from astropy.time import Time\n",
    "from numpy import pi as PI\n",
    "from numpy import arccos, cos, sin, sqrt\n",
    "from scipy.interpolate import interp1d\n",
    "from numpy import argmin, fabs\n",
    "MIN_SEP = 2. * PI/180.\n",
    "def angdist(ra1,dec1,ra2,dec2):\n",
    "    \"\"\"\n",
    "    Computes the angular distance between two objects\n",
    "\n",
    "    Args:\n",
    "        ra1: float (R.A for one of the objects, degrees)\n",
    "        ra2: float (R.A. for the other object, degrees)\n",
    "        dec1: float (DEC for one of the objects, degrees)\n",
    "        dec2: float (DEC for the other object, degrees)\n",
    "\n",
    "    Returns:\n",
    "        float, Angular distance in radians between the objects\n",
    "    \"\"\"\n",
    "    angdist = arccos(sin(dec1*180/PI)*sin(dec2*180/PI)+ \\\n",
    "                     cos(dec1*180/PI)*cos(dec2*180/PI)*cos((ra1-ra2)*180/PI))\n",
    "    return angdist\n",
    "\n",
    "\n",
    "def avoidObject(datetime, ra0, dec0):\n",
    "    \"\"\"\n",
    "    Checks whether all the objects on the list are far enough away from\n",
    "    the input coordinates.\n",
    "    The current list has: Venus, Mars, Jupiter, Saturn, Neptune, Uranus;\n",
    "    the Moon is treated separately.\n",
    "\n",
    "    Args:\n",
    "        datetime: astropy Time; should have timezone info\n",
    "        ra0: float (apparent or observed, degrees)\n",
    "        dec0: float (apparent or observed, degrees)\n",
    "\n",
    "    Returns:\n",
    "        bool, True if all objects on the list are far enough away\n",
    "    \"\"\"\n",
    "\n",
    "    table = Table.read('/Users/javiers/gitrepos/downloads/surveysim/py/surveysim/data/ephemeris_table.fits.gz')\n",
    "    ind_date = argmin(fabs(table['MJD']-datetime.mjd))\n",
    "    k=5\n",
    "    if(ind_date+k>=len(table)):\n",
    "        k=len(table)-ind_date\n",
    "    bodies = ['venus','mars','jupiter','saturn','uranus','neptune']\n",
    "    for body in bodies:\n",
    "        spline_ra = interp1d(table['MJD'][ind_date-k:ind_date+k],table['RA_'+body][ind_date-k:ind_date+k], kind='cubic')\n",
    "        spline_dec = interp1d(table['MJD'][ind_date-k:ind_date+k], table['DEC_'+body][ind_date-k:ind_date+k], kind='cubic')\n",
    "        ra = spline_ra(datetime.mjd)\n",
    "        dec = spline_dec(datetime.mjd)\n",
    "        dist = angdist(ra,dec,ra0,dec0)\n",
    "        if(dist < MIN_SEP):\n",
    "            return False\n",
    "    return True\n",
    "\n",
    "def moonLoc (datetime, ra0, dec0):\n",
    "    \"\"\"\n",
    "    Returns the distance to the Moon if RA and DEC as well as alt, az.\n",
    "\n",
    "    Args:\n",
    "        datetime: astropy Time; should have timezone info\n",
    "        ra0: float (apparent or observed, degrees)\n",
    "        dec0: float (apparent or observed, degrees)\n",
    "\n",
    "    Returns:\n",
    "        float, distance from the Moon (degrees)\n",
    "        float, Moon altitude (degrees)\n",
    "        float, Moon azimuth (degrees)\n",
    "    \"\"\"\n",
    "    table = Table.read('/Users/javiers/gitrepos/downloads/surveysim/py/surveysim/data/moon_ephemeris.fits.gz')\n",
    "    ind_date = argmin(fabs(table['MJD']-datetime.mjd))\n",
    "    k=5\n",
    "    if(ind_date+k>=len(table)):\n",
    "        k=len(table)-ind_date\n",
    "    spline_ra = interp1d(table['MJD'][ind_date-k:ind_date+k],table['RA'][ind_date-k:ind_date+k])\n",
    "    spline_dec = interp1d(table['MJD'][ind_date-k:ind_date+k],table['DEC'][ind_date-k:ind_date+k])\n",
    "    spline_alt = interp1d(table['MJD'][ind_date-k:ind_date+k],table['ALT'][ind_date-k:ind_date+k])\n",
    "    spline_az = interp1d(table['MJD'][ind_date-k:ind_date+k],table['AZ'][ind_date-k:ind_date+k])\n",
    "    ra = spline_ra(datetime.mjd)\n",
    "    dec = spline_dec(datetime.mjd)\n",
    "    alt = spline_alt(datetime.mjd)\n",
    "    az = spline_az(datetime.mjd)\n",
    "    moondist = angdist(ra,dec,ra0,dec0)\n",
    "    return moondist*180.0/PI, alt, az\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 195,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "times = Time(\"2020-09-22 23:22\")\n",
    "time = Time(times, scale='utc', format='isot')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 196,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Separation angle 133.041651396\n",
      "Moon Altitude  30.6975848332\n",
      "Moon Azimuth 166.121826916\n",
      "CPU times: user 87.4 ms, sys: 6.76 ms, total: 94.1 ms\n",
      "Wall time: 94.4 ms\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "mloc = moonLoc(time, ra, dec)\n",
    "print 'Separation angle', mloc[0]\n",
    "print 'Moon Altitude ', mloc[1]\n",
    "print 'Moon Azimuth', mloc[2]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}