Great circles

Great Circles.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-04-08T17:04:37.142633Z",
     "start_time": "2018-04-08T17:04:37.119296Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Using matplotlib backend: MacOSX\n"
     ]
    }
   ],
   "source": [
    "import matplotlib as mpl\n",
    "from mpl_toolkits.mplot3d import Axes3D, proj3d\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib auto"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-04-08T17:04:38.420881Z",
     "start_time": "2018-04-08T17:04:38.126231Z"
    }
   },
   "outputs": [],
   "source": [
    "def circle(p0, p1, ax, **kwargs):\n",
    "    \"\"\"Draw circle passing through p0 and p1\"\"\"\n",
    "    v = np.cross(p0, p1)\n",
    "    u = np.cross(p0, v)\n",
    "    p0 = np.array(p0)/np.sqrt(np.dot(p0, p0))\n",
    "    u = u/np.sqrt(np.dot(u, u))\n",
    "    theta = np.linspace(0, 2*np.pi, 100)\n",
    "    r = np.outer(p0, np.cos(theta)) + np.outer(u, np.sin(theta))\n",
    "    ax.plot(r[0], r[1], r[2], **kwargs)\n",
    "\n",
    "def circle_pole(p0, ax, **kwargs):\n",
    "    \"\"\"Draw great circle with pole at p0\"\"\"\n",
    "    rand = np.random.uniform(size=3)\n",
    "    rand /= np.sqrt(np.dot(rand, rand))\n",
    "    u = np.cross(p0, rand)\n",
    "    v = np.cross(p0, u)\n",
    "    u /= np.sqrt(np.dot(u, u))\n",
    "    v /= np.sqrt(np.dot(v, v))\n",
    "    theta = np.linspace(0, 2*np.pi, 100)\n",
    "    r = np.outer(u, np.cos(theta)) + np.outer(v, np.sin(theta))\n",
    "    ax.plot(r[0], r[1], r[2], **kwargs)\n",
    "\n",
    "def surface_arc(p0, radius, ax, start=0, stop=2*np.pi, label=None, **kwargs):\n",
    "    \"\"\"Draw an arc around `p0` with `radius` from angle `start` to `stop`, where\n",
    "    start and stop are measured CCW from +Alt when viewing from outside sphere.\"\"\"\n",
    "    # Almost the same as circle_pole.\n",
    "    # First find a vector perp to vertical and to p0\n",
    "    u = np.cross(p0, [0,0,1])\n",
    "    u /= np.sqrt(np.dot(u, u))\n",
    "    # Next find a vector perp to u and p0\n",
    "    v = np.cross(u, p0)\n",
    "    v /= np.sqrt(np.dot(v, v))\n",
    "    \n",
    "    theta = np.linspace(start, stop, 100)\n",
    "    # Tricky bit is center of circle is not (0,0,0), but point, and radius != 1\n",
    "    r = np.outer(radius*v, np.cos(theta)) + np.outer(radius*u, np.sin(theta))\n",
    "    r += np.array(p0)[:,None]\n",
    "    ax.plot(r[0], r[1], r[2], **kwargs)\n",
    "\n",
    "    # return label at midway point along arc\n",
    "    if label:\n",
    "        theta = 0.5*(start+stop)\n",
    "        r = radius*v*np.cos(theta) + radius*u*np.sin(theta)\n",
    "        r += np.array(p0)\n",
    "        return r, ax.annotate(label, xy=(0, 0))\n",
    "\n",
    "def solve_triangle(alt, az, lat):\n",
    "    # correct for lat>0 and HA<0\n",
    "    c90malt = np.cos(np.pi/2-alt)\n",
    "    c90mlat = np.cos(np.pi/2-lat)\n",
    "    s90malt = np.sin(np.pi/2-alt)\n",
    "    s90mlat = np.sin(np.pi/2-lat)\n",
    "    caz = np.cos(az)\n",
    "    saz = np.sin(az)\n",
    "    c90mdec = c90malt*c90mlat + s90malt*s90mlat*caz\n",
    "    dec = np.pi/2 - np.arccos(c90mdec)  # This arccos is unambiguous\n",
    "\n",
    "    # arcsin, arccos would be ambiguous here, so use arctan2\n",
    "    s90mdec = np.sin(np.pi/2-dec)\n",
    "    smq = saz/s90mdec*s90mlat\n",
    "    cmq = (c90mlat-c90malt*c90mdec)/s90malt*s90mdec\n",
    "    q = -np.arctan2(smq, cmq)\n",
    "\n",
    "    # arcsin here is ambiguous, use arctan2\n",
    "    smha = saz/s90mdec*s90malt\n",
    "    cmha = (c90malt-c90mlat*c90mdec)/(s90mlat*s90mdec)\n",
    "    ha = -np.arctan2(smha, cmha)\n",
    "\n",
    "    print(\"alt = \", alt*180/np.pi)\n",
    "    print(\"az = \", az*180/np.pi)\n",
    "    print(\"lat = \", lat*180/np.pi)\n",
    "    print(\"ha = \", ha*180/np.pi)\n",
    "    print(\"dec = \", dec*180/np.pi)\n",
    "    print(\"q = \", q*180/np.pi)\n",
    "\n",
    "    return dict(ha=ha, dec=dec, q=q)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-04-08T17:04:39.318415Z",
     "start_time": "2018-04-08T17:04:38.797974Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "alt =  65.6\n",
      "az =  231.36\n",
      "lat =  19.800000000000004\n",
      "ha =  18.8668475078\n",
      "dec =  3.77165715688\n",
      "q =  47.5581097468\n"
     ]
    }
   ],
   "source": [
    "# constants\n",
    "zenith = [0,0,1]\n",
    "nadir = [0,0,-1]\n",
    "horizon_north = [0,1,0]\n",
    "horizon_east = [1,0,0]\n",
    "\n",
    "# inputs and derived\n",
    "lat = 19.8*np.pi/180\n",
    "NCP = [0,np.cos(lat),np.sin(lat)]\n",
    "SCP = [0,-np.cos(lat),-np.sin(lat)]\n",
    "alt = 65.6*np.pi/180\n",
    "az = 231.36*np.pi/180\n",
    "point = [np.cos(alt)*np.sin(az),\n",
    "         np.cos(alt)*np.cos(az),\n",
    "         np.sin(alt)]\n",
    "triangle = solve_triangle(alt, az, lat)\n",
    "q = triangle['q']\n",
    "ha = triangle['ha']\n",
    "dec = triangle['dec']\n",
    "\n",
    "fig = plt.figure(figsize=(8, 7))\n",
    "ax = fig.gca(projection='3d')\n",
    "ax.view_init(10,0)\n",
    "\n",
    "# plot unit sphere\n",
    "phi, theta = np.mgrid[0.0:np.pi:20j, 0.0:2.0*np.pi:20j]\n",
    "r = 1\n",
    "x = r*np.sin(phi)*np.cos(theta)\n",
    "y = r*np.sin(phi)*np.sin(theta)\n",
    "z = r*np.cos(phi)\n",
    "ax.plot_surface(x, y, z,\n",
    "                rstride=1, cstride=1, \n",
    "                color='c', alpha=0.1, linewidth=0)\n",
    "\n",
    "# Plot hour angle=0 circle\n",
    "circle(zenith, horizon_north, ax, c='k')\n",
    "ax.set_xlabel(\"x\")\n",
    "ax.set_ylabel(\"y\")\n",
    "ax.set_zlabel(\"z\")\n",
    "\n",
    "# Plot Zenith/nadir point\n",
    "ax.scatter(*zenith, c='r')\n",
    "ax.scatter(*nadir, c='r')\n",
    "\n",
    "# Plot NCP/SCP point\n",
    "ax.scatter(*NCP, c='b')\n",
    "ax.scatter(*SCP, c='b')\n",
    "\n",
    "# Plot horizon\n",
    "circle(horizon_north, horizon_east, ax, c='r')\n",
    "\n",
    "# Plot *\n",
    "ax.scatter(*point, c='k')\n",
    "\n",
    "# Plot meridian.  Goes through zenith and point\n",
    "circle(zenith, point, ax, c='r')\n",
    "\n",
    "# Plot hour circle through NCP and point\n",
    "circle(NCP, point, ax, c='b')\n",
    "\n",
    "# Plot equator.\n",
    "circle_pole(NCP, ax, c='b')\n",
    "\n",
    "# Label the hour angle\n",
    "if ha > 0:\n",
    "    halabel = surface_arc(NCP, 0.15, ax, start=0, stop=-ha, label=\"HA\")\n",
    "else:\n",
    "    halabel = surface_arc(NCP, 0.15, ax, start=0, stop=-ha, label=\"-HA\")\n",
    "\n",
    "# Label the parallactic angle\n",
    "if q > 0:\n",
    "    qlabel = surface_arc(point, 0.15, ax, start=0, stop=q, label=\"q\")\n",
    "else:\n",
    "    qlabel = surface_arc(point, 0.15, ax, start=0, stop=q, label=\"-q\")\n",
    "\n",
    "# Add labels\n",
    "labels = []\n",
    "labels.append((NCP, ax.annotate(\"NCP\", xy=(0,0))))\n",
    "labels.append((SCP, ax.annotate(\"SCP\", xy=(0,0))))\n",
    "labels.append((zenith, ax.annotate(\"zenith\", xy=(0,0))))\n",
    "labels.append((nadir, ax.annotate(\"nadir\", xy=(0,0))))\n",
    "labels.append((point, ax.annotate(\"*\", xy=(0,0))))\n",
    "labels.append(halabel)\n",
    "labels.append(qlabel)\n",
    "\n",
    "def update_position(labels, fig, ax):\n",
    "    for p, label in labels:\n",
    "        x, y, _ = proj3d.proj_transform(p[0], p[1], p[2], ax.get_proj())\n",
    "        label.set_x(x)\n",
    "        label.set_y(y)\n",
    "    fig.canvas.draw()\n",
    "update_position(labels, fig, ax)\n",
    "fig.canvas.mpl_connect('motion_notify_event', lambda event: update_position(labels, fig, ax))\n",
    "pass"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "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.6.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}