{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Pyephem"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Pyephem is a very useful tool for calculating positions of objects on the sky for an earth-based observer. We can simply calculate observabilities of objects at given RA, DEC positions, convert between epochs and coordinate systems and even track asteroids, solar-system bodies and extrasolar objects with given proper motions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import ephem\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import time"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We start by defining an observer, which is an object containing a set of lat/lon coordinates on the Earth, a date (defaults to now) and an observing epoch (defaults to 2000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'11:22'"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "obs = ephem.city('Santiago')\n",
    "ephem.localtime(obs.date).strftime(\"%H:%M\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's try with a predefined object:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "obj = ephem.Jupiter()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Nothing happens until we explicitly tell ephem to compute the position:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "33:57:58.8 61:52:25.1\n"
     ]
    }
   ],
   "source": [
    "obj.compute(obs)\n",
    "print obj.alt, obj.az"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can compute the next rise/set times using next_rising(), next_setting()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2016/9/9 13:07:29\n",
      "2016/9/8 21:50:25\n"
     ]
    }
   ],
   "source": [
    "obs.horizon = '20'#To take into account the cordillera, but this could also be your observing limits\n",
    "print obs.next_rising(obj)\n",
    "print obs.next_setting(obj)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "By defining the sun as an object, we can calculate the start and end of astronomical twilight:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2016/9/8 23:52:09\n",
      "2016/9/9 09:26:30\n"
     ]
    }
   ],
   "source": [
    "sol = ephem.Sun()\n",
    "obs.horizon = '-18'\n",
    "sol.compute(obs)\n",
    "print obs.next_setting(sol, use_center=True)\n",
    "print obs.next_rising(sol, use_center=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What about custom objects? We can define those with XEphem formatted strings http://www.clearskyinstitute.com/xephem/help/xephem.html#mozTocId468501"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-20:53:37.1 242:50:58.6\n",
      "2016/9/9 01:22:13 2016/9/9 10:57:14\n"
     ]
    }
   ],
   "source": [
    "ra, dec = \"00:41:30.347\",\"-09:15:45,14.46\"\n",
    "JO201 = ephem.readdb(\"JO201,f,\"+ra+\",\"+dec)\n",
    "JO201.compute(obs)\n",
    "print JO201.alt, JO201.az\n",
    "obs.horizon = '20'\n",
    "print obs.next_rising(JO201, use_center=True),obs.next_setting(JO201, use_center=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Coordinate transformations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Coordinate transformations with ephem are a breeze. We can work with Galactic lat/lon:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-71:58:51.1 114:59:11.8\n"
     ]
    }
   ],
   "source": [
    "JO201_gal = ephem.Galactic(JO201)\n",
    "print JO201_gal.lat, JO201_gal.lon"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Or we can change the epoch:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0:38:58.55 -9:32:11.8\n"
     ]
    }
   ],
   "source": [
    "JO201.compute(ephem.now(),epoch='1950')\n",
    "print JO201.a_ra, JO201.a_dec"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Try a quick example: Plotting the altitude over 24h"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "resolution = 0.1\n",
    "datelist = [ephem.date(obs.date + x*ephem.hour) for x in np.arange(0,24,resolution)]\n",
    "positions_JO201 = [0]*len(datelist)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "for i in range(len(datelist)):\n",
    "    obs.date = datelist[i]\n",
    "    JO201.compute(obs)\n",
    "    positions_JO201[i] = JO201.alt*(180./3.1415)\n",
    "\n",
    "plt.plot(np.arange(0,24,resolution),positions_JO201,'-k')\n",
    "plt.gca().set_ylim(0,90)\n",
    "plt.gca().set_xlim(0,24)\n",
    "plt.xticks(np.arange(0,25),rotation='vertical')\n",
    "plt.gca().set_xticklabels([ephem.localtime(d).strftime(\"%H:%M\") for d in datelist][::int(1/resolution)])\n",
    "plt.grid()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Calculations with angles"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also check the lunar distance right now, useful to know if we're in grey time"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "119:27:43.3\n"
     ]
    }
   ],
   "source": [
    "obs.date = ephem.now()\n",
    "moon = ephem.Moon()\n",
    "\n",
    "JO201.compute(obs)\n",
    "moon.compute(obs)\n",
    "\n",
    "print ephem.separation(JO201, moon)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "or for a given date/time:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "165:12:27.7\n"
     ]
    }
   ],
   "source": [
    "d = ephem.Date('2016/09/30 01:00:00')\n",
    "JO201.compute(d)\n",
    "moon.compute(d)\n",
    "\n",
    "print ephem.separation(JO201, moon)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## A bit of fun, projects for your own enjoyment"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Current planets visible:\n",
      "Mercury: Leo (6 deg E of sun)\n",
      "  Venus: Virgo \n",
      "Jupiter: Virgo (13 deg E of sun)\n"
     ]
    }
   ],
   "source": [
    "sol = ephem.Sun()\n",
    "planets = [ephem.Mercury(),ephem.Venus(),ephem.Mars(),ephem.Jupiter(),ephem.Saturn(),ephem.Moon()]\n",
    "\n",
    "home = ephem.city(\"Santiago\")\n",
    "\n",
    "sol.compute(home)\n",
    "\n",
    "print \"Current planets visible:\"\n",
    "\n",
    "for planet in planets:\n",
    "    planet.compute(home)\n",
    "    if planet.alt > ephem.degrees('10'):\n",
    "        print \"{:>7}: {}\".format(str(planet.name), str(ephem.constellation(planet)[1])),\n",
    "        dec = ephem.degrees((planet.dec+sol.dec)/2.)\n",
    "        sep = 360/(2*np.pi) * ((planet.ra - sol.ra)*np.cos(dec))#+ve: east\n",
    "        dir = 'E'\n",
    "        if sep < 0:\n",
    "            dir = 'W'\n",
    "        if np.abs(sep) < 20:\n",
    "            print \"({:.0f} deg {} of sun)\".format(np.abs(sep),dir)\n",
    "        else: print ''"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Upcoming ISS passes:\n",
      "14:40 on 08/09/16 (2 deg elev.)\n",
      "16:15 on 08/09/16 (47 deg elev.)\n",
      "17:52 on 08/09/16 (16 deg elev.)\n",
      "19:30 on 08/09/16 (4 deg elev.)\n",
      "21:08 on 08/09/16 (3 deg elev.)\n"
     ]
    }
   ],
   "source": [
    "rad2deg = 360/6.28318\n",
    "#http://spaceflight.nasa.gov/realdata/sightings/SSapplications/Post/JavaSSOP/orbit/ISS/SVPOST.html\n",
    "ISS = ephem.readtle(\"ISS\",\n",
    "    \"1 25544U 98067A   16246.52721698  .00016717  00000-0  10270-3 0  9006\",\n",
    "    \"2 25544  51.6415  32.1269 0002931 280.0688  80.0134 15.54400533 17019\")\n",
    "\n",
    "home = ephem.city(\"Santiago\")\n",
    "\n",
    "print \"Upcoming ISS passes:\"\n",
    "\n",
    "for i in range(5):\n",
    "    ISS.compute(home)\n",
    "    nextpass = home.next_pass(ISS)\n",
    "    print (ephem.localtime(nextpass[2])).strftime(\"%H:%M\"),\"on\",(ephem.localtime(nextpass[2])).strftime(\"%d/%m/%y\"),\"({:.0f} deg elev.)\".format(nextpass[3]*rad2deg)\n",
    "    home.date = nextpass[4]+0.0001"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import ephem\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import datetime\n",
    "\n",
    "planets = [ephem.Mercury(),ephem.Venus(),ephem.Mars(),ephem.Jupiter(),ephem.Saturn(),ephem.Moon()]\n",
    "\n",
    "ISS = ephem.readtle(\"ISS\",\n",
    "    \"1 25544U 98067A   16246.52721698  .00016717  00000-0  10270-3 0  9006\",\n",
    "    \"2 25544  51.6415  32.1269 0002931 280.0688  80.0134 15.54400533 17019\")\n",
    "\n",
    "home = ephem.city(\"Santiago\")\n",
    "\n",
    "fig = plt.figure()\n",
    "ax = fig.add_subplot(111,projection='polar')\n",
    "ax.set_rmax(90)\n",
    "\n",
    "for obj in planets:\n",
    "    obj.compute(home)\n",
    "    if obj.alt > ephem.degrees('0'):\n",
    "        theta,r = obj.az, rad2deg*(ephem.degrees('90')-ephem.degrees(obj.alt))\n",
    "        plt.polar(theta,r,'o',mew=0)\n",
    "        plt.text(theta,r,obj.name,color='k', horizontalalignment='left',verticalalignment='bottom',fontsize=10)\n",
    "        ax.set_theta_offset(1.5*np.pi)\n",
    "\n",
    "for i in range(3):\n",
    "    length = (ephem.date(home.next_pass(ISS)[4]) - ephem.date(home.next_pass(ISS)[0]))*24*3600\n",
    "    #print length\n",
    "    dates = [ephem.date(home.next_pass(ISS)[0]+ephem.second*x) for x in range(0,int(length),5)]\n",
    "    home.date = dates[0]\n",
    "    ISS.compute(home)\n",
    "    theta, r = ISS.az, rad2deg*(ephem.degrees('90')-ephem.degrees(ISS.alt))\n",
    "    plt.text(theta,r,ephem.localtime(home.date).strftime(\"%H:%M\"),color='k', horizontalalignment='right',verticalalignment='top',fontsize=10)\n",
    "    \n",
    "    for datep in dates:\n",
    "        home.date = datep\n",
    "        ISS.compute(home)\n",
    "        theta, r = ISS.az, rad2deg*(ephem.degrees('90')-ephem.degrees(ISS.alt))\n",
    "        \n",
    "        #print theta, r\n",
    "        plt.polar(theta, r, '.',c='grey',markersize=2)\n",
    "    home.date = ephem.date(home.next_pass(ISS)[4])+0.001\n",
    "\n",
    "\n",
    "ax.set_thetagrids(range(0,360,30),fontsize=10)\n",
    "ax.set_rgrids(range(15,105,15),[75,60,45,30,15,0],fontsize=10)\n",
    "ax.set_rlabel_position(0)\n",
    "ax.set_rlim(0,90)\n",
    "plt.grid(color=\"k\")\n",
    "plt.show()"
   ]
  },
  {
   "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.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}