daily pastebin goal
62%
SHARE
TWEET

Untitled

a guest Apr 16th, 2018 51 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. {
  2.  "cells": [
  3.   {
  4.    "cell_type": "code",
  5.    "execution_count": 4,
  6.    "metadata": {
  7.     "ExecuteTime": {
  8.      "end_time": "2018-04-08T17:04:37.142633Z",
  9.      "start_time": "2018-04-08T17:04:37.119296Z"
  10.     }
  11.    },
  12.    "outputs": [
  13.     {
  14.      "name": "stdout",
  15.      "output_type": "stream",
  16.      "text": [
  17.       "Using matplotlib backend: MacOSX\n"
  18.      ]
  19.     }
  20.    ],
  21.    "source": [
  22.     "import matplotlib as mpl\n",
  23.     "from mpl_toolkits.mplot3d import Axes3D, proj3d\n",
  24.     "import numpy as np\n",
  25.     "import matplotlib.pyplot as plt\n",
  26.     "%matplotlib auto"
  27.    ]
  28.   },
  29.   {
  30.    "cell_type": "code",
  31.    "execution_count": 5,
  32.    "metadata": {
  33.     "ExecuteTime": {
  34.      "end_time": "2018-04-08T17:04:38.420881Z",
  35.      "start_time": "2018-04-08T17:04:38.126231Z"
  36.     }
  37.    },
  38.    "outputs": [],
  39.    "source": [
  40.     "def circle(p0, p1, ax, **kwargs):\n",
  41.     "    \"\"\"Draw circle passing through p0 and p1\"\"\"\n",
  42.     "    v = np.cross(p0, p1)\n",
  43.     "    u = np.cross(p0, v)\n",
  44.     "    p0 = np.array(p0)/np.sqrt(np.dot(p0, p0))\n",
  45.     "    u = u/np.sqrt(np.dot(u, u))\n",
  46.     "    theta = np.linspace(0, 2*np.pi, 100)\n",
  47.     "    r = np.outer(p0, np.cos(theta)) + np.outer(u, np.sin(theta))\n",
  48.     "    ax.plot(r[0], r[1], r[2], **kwargs)\n",
  49.     "\n",
  50.     "def circle_pole(p0, ax, **kwargs):\n",
  51.     "    \"\"\"Draw great circle with pole at p0\"\"\"\n",
  52.     "    rand = np.random.uniform(size=3)\n",
  53.     "    rand /= np.sqrt(np.dot(rand, rand))\n",
  54.     "    u = np.cross(p0, rand)\n",
  55.     "    v = np.cross(p0, u)\n",
  56.     "    u /= np.sqrt(np.dot(u, u))\n",
  57.     "    v /= np.sqrt(np.dot(v, v))\n",
  58.     "    theta = np.linspace(0, 2*np.pi, 100)\n",
  59.     "    r = np.outer(u, np.cos(theta)) + np.outer(v, np.sin(theta))\n",
  60.     "    ax.plot(r[0], r[1], r[2], **kwargs)\n",
  61.     "\n",
  62.     "def surface_arc(p0, radius, ax, start=0, stop=2*np.pi, label=None, **kwargs):\n",
  63.     "    \"\"\"Draw an arc around `p0` with `radius` from angle `start` to `stop`, where\n",
  64.     "    start and stop are measured CCW from +Alt when viewing from outside sphere.\"\"\"\n",
  65.     "    # Almost the same as circle_pole.\n",
  66.     "    # First find a vector perp to vertical and to p0\n",
  67.     "    u = np.cross(p0, [0,0,1])\n",
  68.     "    u /= np.sqrt(np.dot(u, u))\n",
  69.     "    # Next find a vector perp to u and p0\n",
  70.     "    v = np.cross(u, p0)\n",
  71.     "    v /= np.sqrt(np.dot(v, v))\n",
  72.     "    \n",
  73.     "    theta = np.linspace(start, stop, 100)\n",
  74.     "    # Tricky bit is center of circle is not (0,0,0), but point, and radius != 1\n",
  75.     "    r = np.outer(radius*v, np.cos(theta)) + np.outer(radius*u, np.sin(theta))\n",
  76.     "    r += np.array(p0)[:,None]\n",
  77.     "    ax.plot(r[0], r[1], r[2], **kwargs)\n",
  78.     "\n",
  79.     "    # return label at midway point along arc\n",
  80.     "    if label:\n",
  81.     "        theta = 0.5*(start+stop)\n",
  82.     "        r = radius*v*np.cos(theta) + radius*u*np.sin(theta)\n",
  83.     "        r += np.array(p0)\n",
  84.     "        return r, ax.annotate(label, xy=(0, 0))\n",
  85.     "\n",
  86.     "def solve_triangle(alt, az, lat):\n",
  87.     "    # correct for lat>0 and HA<0\n",
  88.     "    c90malt = np.cos(np.pi/2-alt)\n",
  89.     "    c90mlat = np.cos(np.pi/2-lat)\n",
  90.     "    s90malt = np.sin(np.pi/2-alt)\n",
  91.     "    s90mlat = np.sin(np.pi/2-lat)\n",
  92.     "    caz = np.cos(az)\n",
  93.     "    saz = np.sin(az)\n",
  94.     "    c90mdec = c90malt*c90mlat + s90malt*s90mlat*caz\n",
  95.     "    dec = np.pi/2 - np.arccos(c90mdec)  # This arccos is unambiguous\n",
  96.     "\n",
  97.     "    # arcsin, arccos would be ambiguous here, so use arctan2\n",
  98.     "    s90mdec = np.sin(np.pi/2-dec)\n",
  99.     "    smq = saz/s90mdec*s90mlat\n",
  100.     "    cmq = (c90mlat-c90malt*c90mdec)/s90malt*s90mdec\n",
  101.     "    q = -np.arctan2(smq, cmq)\n",
  102.     "\n",
  103.     "    # arcsin here is ambiguous, use arctan2\n",
  104.     "    smha = saz/s90mdec*s90malt\n",
  105.     "    cmha = (c90malt-c90mlat*c90mdec)/(s90mlat*s90mdec)\n",
  106.     "    ha = -np.arctan2(smha, cmha)\n",
  107.     "\n",
  108.     "    print(\"alt = \", alt*180/np.pi)\n",
  109.     "    print(\"az = \", az*180/np.pi)\n",
  110.     "    print(\"lat = \", lat*180/np.pi)\n",
  111.     "    print(\"ha = \", ha*180/np.pi)\n",
  112.     "    print(\"dec = \", dec*180/np.pi)\n",
  113.     "    print(\"q = \", q*180/np.pi)\n",
  114.     "\n",
  115.     "    return dict(ha=ha, dec=dec, q=q)"
  116.    ]
  117.   },
  118.   {
  119.    "cell_type": "code",
  120.    "execution_count": 6,
  121.    "metadata": {
  122.     "ExecuteTime": {
  123.      "end_time": "2018-04-08T17:04:39.318415Z",
  124.      "start_time": "2018-04-08T17:04:38.797974Z"
  125.     }
  126.    },
  127.    "outputs": [
  128.     {
  129.      "name": "stdout",
  130.      "output_type": "stream",
  131.      "text": [
  132.       "alt =  65.6\n",
  133.       "az =  231.36\n",
  134.       "lat =  19.800000000000004\n",
  135.       "ha =  18.8668475078\n",
  136.       "dec =  3.77165715688\n",
  137.       "q =  47.5581097468\n"
  138.      ]
  139.     }
  140.    ],
  141.    "source": [
  142.     "# constants\n",
  143.     "zenith = [0,0,1]\n",
  144.     "nadir = [0,0,-1]\n",
  145.     "horizon_north = [0,1,0]\n",
  146.     "horizon_east = [1,0,0]\n",
  147.     "\n",
  148.     "# inputs and derived\n",
  149.     "lat = 19.8*np.pi/180\n",
  150.     "NCP = [0,np.cos(lat),np.sin(lat)]\n",
  151.     "SCP = [0,-np.cos(lat),-np.sin(lat)]\n",
  152.     "alt = 65.6*np.pi/180\n",
  153.     "az = 231.36*np.pi/180\n",
  154.     "point = [np.cos(alt)*np.sin(az),\n",
  155.     "         np.cos(alt)*np.cos(az),\n",
  156.     "         np.sin(alt)]\n",
  157.     "triangle = solve_triangle(alt, az, lat)\n",
  158.     "q = triangle['q']\n",
  159.     "ha = triangle['ha']\n",
  160.     "dec = triangle['dec']\n",
  161.     "\n",
  162.     "fig = plt.figure(figsize=(8, 7))\n",
  163.     "ax = fig.gca(projection='3d')\n",
  164.     "ax.view_init(10,0)\n",
  165.     "\n",
  166.     "# plot unit sphere\n",
  167.     "phi, theta = np.mgrid[0.0:np.pi:20j, 0.0:2.0*np.pi:20j]\n",
  168.     "r = 1\n",
  169.     "x = r*np.sin(phi)*np.cos(theta)\n",
  170.     "y = r*np.sin(phi)*np.sin(theta)\n",
  171.     "z = r*np.cos(phi)\n",
  172.     "ax.plot_surface(x, y, z,\n",
  173.     "                rstride=1, cstride=1, \n",
  174.     "                color='c', alpha=0.1, linewidth=0)\n",
  175.     "\n",
  176.     "# Plot hour angle=0 circle\n",
  177.     "circle(zenith, horizon_north, ax, c='k')\n",
  178.     "ax.set_xlabel(\"x\")\n",
  179.     "ax.set_ylabel(\"y\")\n",
  180.     "ax.set_zlabel(\"z\")\n",
  181.     "\n",
  182.     "# Plot Zenith/nadir point\n",
  183.     "ax.scatter(*zenith, c='r')\n",
  184.     "ax.scatter(*nadir, c='r')\n",
  185.     "\n",
  186.     "# Plot NCP/SCP point\n",
  187.     "ax.scatter(*NCP, c='b')\n",
  188.     "ax.scatter(*SCP, c='b')\n",
  189.     "\n",
  190.     "# Plot horizon\n",
  191.     "circle(horizon_north, horizon_east, ax, c='r')\n",
  192.     "\n",
  193.     "# Plot *\n",
  194.     "ax.scatter(*point, c='k')\n",
  195.     "\n",
  196.     "# Plot meridian.  Goes through zenith and point\n",
  197.     "circle(zenith, point, ax, c='r')\n",
  198.     "\n",
  199.     "# Plot hour circle through NCP and point\n",
  200.     "circle(NCP, point, ax, c='b')\n",
  201.     "\n",
  202.     "# Plot equator.\n",
  203.     "circle_pole(NCP, ax, c='b')\n",
  204.     "\n",
  205.     "# Label the hour angle\n",
  206.     "if ha > 0:\n",
  207.     "    halabel = surface_arc(NCP, 0.15, ax, start=0, stop=-ha, label=\"HA\")\n",
  208.     "else:\n",
  209.     "    halabel = surface_arc(NCP, 0.15, ax, start=0, stop=-ha, label=\"-HA\")\n",
  210.     "\n",
  211.     "# Label the parallactic angle\n",
  212.     "if q > 0:\n",
  213.     "    qlabel = surface_arc(point, 0.15, ax, start=0, stop=q, label=\"q\")\n",
  214.     "else:\n",
  215.     "    qlabel = surface_arc(point, 0.15, ax, start=0, stop=q, label=\"-q\")\n",
  216.     "\n",
  217.     "# Add labels\n",
  218.     "labels = []\n",
  219.     "labels.append((NCP, ax.annotate(\"NCP\", xy=(0,0))))\n",
  220.     "labels.append((SCP, ax.annotate(\"SCP\", xy=(0,0))))\n",
  221.     "labels.append((zenith, ax.annotate(\"zenith\", xy=(0,0))))\n",
  222.     "labels.append((nadir, ax.annotate(\"nadir\", xy=(0,0))))\n",
  223.     "labels.append((point, ax.annotate(\"*\", xy=(0,0))))\n",
  224.     "labels.append(halabel)\n",
  225.     "labels.append(qlabel)\n",
  226.     "\n",
  227.     "def update_position(labels, fig, ax):\n",
  228.     "    for p, label in labels:\n",
  229.     "        x, y, _ = proj3d.proj_transform(p[0], p[1], p[2], ax.get_proj())\n",
  230.     "        label.set_x(x)\n",
  231.     "        label.set_y(y)\n",
  232.     "    fig.canvas.draw()\n",
  233.     "update_position(labels, fig, ax)\n",
  234.     "fig.canvas.mpl_connect('motion_notify_event', lambda event: update_position(labels, fig, ax))\n",
  235.     "pass"
  236.    ]
  237.   },
  238.   {
  239.    "cell_type": "code",
  240.    "execution_count": null,
  241.    "metadata": {},
  242.    "outputs": [],
  243.    "source": []
  244.   }
  245.  ],
  246.  "metadata": {
  247.   "anaconda-cloud": {},
  248.   "kernelspec": {
  249.    "display_name": "Python 3",
  250.    "language": "python",
  251.    "name": "python3"
  252.   },
  253.   "language_info": {
  254.    "codemirror_mode": {
  255.     "name": "ipython",
  256.     "version": 3
  257.    },
  258.    "file_extension": ".py",
  259.    "mimetype": "text/x-python",
  260.    "name": "python",
  261.    "nbconvert_exporter": "python",
  262.    "pygments_lexer": "ipython3",
  263.    "version": "3.6.2"
  264.   }
  265.  },
  266.  "nbformat": 4,
  267.  "nbformat_minor": 1
  268. }
RAW Paste Data
Top