API tools faq
paste
Guest User

Untitled

a guest
Apr 16th, 2018
82
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
text 8.43 KB | None | 0 0
raw download clone embed print report
  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. }
Add Comment
Please, Sign In to add comment
Public Pastes
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.  OK, I Understand
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!