• API
• FAQ
• Tools
• Archive
daily pastebin goal
14%
SHARE
TWEET

# Untitled

a guest Apr 16th, 2018 52 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,
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,
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",
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",
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,
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",
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",
183.     "ax.scatter(*zenith, 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",
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",
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,
242.    "outputs": [],
243.    "source": []
244.   }
245.  ],
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
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.

Top