jobovy · December 22, 2020 21:14 · salmanhiro · Jul 3, 2025 · jobovy · Jul 14, 2025
diff --git a/align-stream-phi12.ipynb b/align-stream-phi12.ipynb
 {
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "/Users/bovy/tmp/galpy-spheredf/galpy/util/bovy_conversion.py:6: FutureWarning: galpy.util.bovy_conversion is being deprecated in favor of galpy.util.conversion; all functions in there are the same; please switch to the new import, because the old import will be removed in v1.9\n",
      "  warnings.warn('galpy.util.bovy_conversion is being deprecated in favor of galpy.util.conversion; all functions in there are the same; please switch to the new import, because the old import will be removed in v1.9',FutureWarning)\n",
      "\n",
      "/Users/bovy/tmp/galpy-spheredf/galpy/util/bovy_coords.py:6: FutureWarning: galpy.util.bovy_coords is being deprecated in favor of galpy.util.coords; all functions in there are the same; please switch to the new import, because the old import will be removed in v1.9\n",
      "  warnings.warn('galpy.util.bovy_coords is being deprecated in favor of galpy.util.coords; all functions in there are the same; please switch to the new import, because the old import will be removed in v1.9',FutureWarning)\n",
      "\n",
      "Populating the interactive namespace from numpy and matplotlib\n"
     ]
    }
   ],
   "source": [
    "import numpy\n",
    "from astropy import units\n",
    "from galpy.potential import LogarithmicHaloPotential\n",
    "from galpy.orbit import Orbit\n",
    "from galpy.util import conversion, coords\n",
    "from galpy.util import plot as galpy_plot\n",
    "galpy_plot.start_print(axes_labelsize=17.,text_fontsize=12.,\n",
    "                     xtick_labelsize=15.,ytick_labelsize=15.)\n",
    "from streamtools.df import streamspraydf\n",
    "from galpy.util import _rotate_to_arbitrary_vector\n",
    "%pylab inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Approximately aligning a stream horizontally\n",
    "\n",
    "In this brief notebook, I illustrate a simple way to transform the sky coordinates of a stream of stars from RA and Dec to $(\\phi_1,\\phi_2)$ such that the stream lies approximately along a constant value of $\\phi_2$. I do this by rotating the sky coordinate frame such that the stream progenitor's angular momentum measured with respect to a non-moving observer at the Sun aligns with the $z$ axis (non-moving in Galactocentric coordinates). I use the example from the ``streamtools`` package.\n",
    "\n",
    "First we set up the stream as in the ``streamtools`` example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "o= Orbit([1.56148083,0.35081535,-1.15481504,0.88719443,-0.47713334,0.12019596])\n",
    "lp= LogarithmicHaloPotential(normalize=1.,q=0.9)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "spdf= streamspraydf(2*10.**4.*units.Msun,progenitor=o,pot=lp,tdisrupt=4.5*units.Gyr)\n",
    "spdft= streamspraydf(2*10.**4.*units.Msun,progenitor=o,pot=lp,leading=False,tdisrupt=4.5*units.Gyr)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then we generate a leading and trailing stream:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "numpy.random.seed(4)\n",
    "RvR,dt= spdf.sample(n=100,returndt=True,integrate=True)\n",
    "RvRt,dt= spdft.sample(n=100,returndt=True,integrate=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We load the stream particles into an ``Orbit`` instance such that we can easily convert coordinates:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "streaml= Orbit(RvR.T)\n",
    "streamt= Orbit(RvRt.T)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following is then the stream in RA and Dec.:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "galpyWarning: Method ra(.) requires ro to be given at Orbit initialization or at method evaluation; using default ro which is 8.000000 kpc\n",
      "galpyWarning: Method dec(.) requires ro to be given at Orbit initialization or at method evaluation; using default ro which is 8.000000 kpc\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAp0AAAF6CAYAAABbfOOZAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3dT2wb6X3/8Q8pW/KpkKwtEGwAF+Aenl72EK3dgofkp8AyilwaQOstb1psUslNsgFctLDgQ4EFUkCxsW2NZJPAarxGdGPsVQ8t0CaWscz2MEC6dgsEPTwHE+geigBdWbrFki3xd3hmzNFoSJHSzJCceb+ABT1/JI5mJfLD7/Ov1Gq1BAAAAKSpPOgLAAAAQP4ROgEAAJA6QicAAABSR+gEAABA6gidAAAASB2hEwAAAKk7NegLiPPHf/zHrS9+8YuDvoyXtra2NDU1NejLGAncq/5wv/rD/eod96o/3K/+cL/6U7T79d///d+fW2t/P7p/KEPnF7/4Ra2vrw/6Ml6an58fqusZZtyr/nC/+sP96h33qj/cr/5wv/pTtPtljPmfuP00rwMAACB1hE4AAACkjtDZg1qtNuhLGBncq/5wv/rD/eod96o/3K/+cL/6w/1ySsO49vr8/HyrSH0fAAAA8sIY88haez66n0onAAAAUkfoBAAAQOoInQAAAEgdoRMAAACpI3QCAAAgdYROAAAApI7QCQAAgNQN5drrAID0eJ60tib99rfSF74gfelL0uamNDsrVasHz2s0pOlpdzx4jJ4X/p6StLAQf7zRiP/ak/4s4WtM+vsDSE4qodMYU5F0WdK2JFlrV0PHrklqSqpI2rDWPk7jGgCgSDqFuuh+z5O+8hXpxYuDX18uSxMT0sOH7fMuXpR2dqT9falUklqt+PPW1qSf/rT9PT/80D1ncB3B99rdlcbH21/by8/R7edaW3PP9eKFu8ZyWTp9WnrnHRd8pcNfm1b4BXC0xEOnHzhvWGvf8rcfGWM+tdY+Nsbck7QSBE1jzANJl5K+BgAYRccNRJ1CXdz+tbXDgVNyoW13tx0WGw23vb/vjgeL14XPk9z3f/asfVySnj8/GDqD77W3d/A5oj9DECL39tz13rolXb3a+eeKPu/+vgvJt2+771MquZ81+NrgeoPvd+sW1VEgS2lUOm9LuhHavmit3fb/PReEUV/TGDNnrd1I4ToAYOA6NTvHNV3HBaxedAp1cfs7KZfd887Ouu3ZWbcdVDrL5fZjcF7w/aOrKZfL7e8T/l7BzxY+FtyLaIjc3ZU++qj7zxV93qAa22q54Cu5f4d/9uD77exI3/mOOx4OpTdvSv/7v9I3vyktLbWvj+oocHKJhk5jzKRcsHxZvQwCpzFmTq5ZPWxbrtJJ6AQw8qIB8ze/kb71rXa18O5d6eOP3b/DTdflsjQ25sJQtOLYi06hrtP+u3fdvlJJ+tM/lb72tcMVv2rVBbGj+nSGg6nkfo4f//jgtYe/V1xwi4bIUsl93zfflP7937v/XKdOueb0L31J+s//bFdKx8YOVjqDrw2+rlw+eL/X1qQ7d9ph9de/bl/fu++6c6PdCgiiQH+SrnRWJG37AXPS337sVzInY87flHQh4WsAgMx5ngsgu7tu+86ddqgJRCtuwbHgsVxuB65oNbCbTqGu0/6PP+4tMFWrRx+PBtNO37Pb94oLkUFV+PXXe/+5JPd1wX7p8Dnh6w1XlqV24AzcuSM9ftzujrCzc7BbwXGq0kCRpRE6Jelp0GTu9+l8S9LZXr/J1taW5ufnX27XajXVarVELxQAjtLr4BzJbYdDSzTASAebnaNN1xMTJ+tj2CnUxe0/Kkwm8bz9fo9OIbKfnytuf7efPRxopYOVTkl69VXp00/b22NjB7sVdOujGtbr4ChgVNXrddXr9fCuV+LOSzp0bkuajIxIb0q6Iuk/Ys6fjvsmU1NTWl9fT/jSAOCgbm/+/QzOqVbd9zh9ul3pPH3aPQYhplw+2Ozca4WwKJIMwsd9zl/96mCfztdfl37xC/fhoFyWPvjgYLeCTn1Uw6K/L50GRwGjLFocNMZ8Hnde0qGzKX+apMi+iqQHim9ij/bzBIBUBX0v7949OLo5/Obfz+CcILw0Ggf7dEqd564cRMhCd9Wq9E//dHBfP90W4kR/XzoNjgKKINHQaa1t+oOJwiYlNa21G8aYaBN7RW60OwBkotNI6eibf7+Dc6TOTdkYXf0270dFf186DY4CiiCNKZNuRqZBOi8pmCZpwxgzE2p+rzBdEoCkHLUqjtR5pHT0zb/fwTlAnLjfl7jBUf2gTyhGVakVnegsAcaYG5KeSHpNUj00GfykpOty/TsvhI+Fzc/Pt+jTCaBXcaviTEy4Udrd+mqOjUnf+EbngAoMm35WdwIGxRjzyFp7Pro/lWUwrbXLHfZvSwqO3U/juQEUS6fVaTr1l6NSiVHWbeQ8FVAMu1RCJwBkpdPqNN36yzGIB6OqU59iKqAYBYROAAORVFUmOrH4174mfeELNJkjnzpV6vudO/QoVE2RBkIngEwF/S/v3GlPVxTX97JXNJejaOIq9Uetb98PqqZIC6ETQGaiS0VKbuLtmzcPz4/YD5rLUXRJfvhKumoKBAidADKztnYwcAb++Z9dIOWNDTi+pD58JVk1BcLKg74AAPnjedLKinvsRavlqikABi+omn7vezStI1lUOgGcWHjQgdS5P9jCglt6cnfXrWUtucA5MUE1BRgm3aqmq6uuT/arr0rXrhFK0TtCJ4ATiQ46ePvtdn+wZ89ck3p4NZ+PPz4YUBkABIyO1VXpypX29r/8i/TJJ+7f/C3jKIROAF0dNXVKdNCB5KYu2ttzVcwPPzw4fVG0gsIbFDA6Pvro4PaLF+6D5c9+1p627J13mLIM8ejTCaCjoIr5N3/jHuP6aAaDDsbG3OPCgnvTKZXc8b09+msCefHmmwe3T/mlq+CD586OdPu2e1341rd679eNYqDSCeCQoLr52WdHT53SaaqWoPLB6FcgP5aW3GO4T6fk/t6DpWhbLfe3f/u2a+n4xjeofMIptaJrxw2B+fn51vr6+qAvAyikcB/NU6fcG8jeXv+TRLOiCVAcwaIPH34oPX9+cFnaUkk6c8a9fkj06S4CY8wja+356H4qnQAOCPfRlKTFRencuf7fGJiwHSiO4O99YcGFz2CWinDlM9r387gfaDG6CJ0AXvI816Q+Nua2gz6avBkA6EVc+AyWu5XaH2j39912OJBS9cw/QidQQEFTmNQOldFm9cVFAieA4wmHz3Bzelylc2zsYDil6plfhE6gYKLrn9+92547M9ysfu4cL/wATibazSY86FBqD1j8x39krfciIHQCORcd0NNouI7+geBFnvWWAaQtbp5ez2O2i6IgdAI5FjSZ7+y4JqwPPnAv6KdPtyudwYt8p6mPACBN/bz2MCvGaCN0AjkSXQP9vffac+ft70vf+Y5bsq7RONynU2LEOYDB6OW1J7rkLn0/Rw+hE8iJ1VUXKvf3XSWz1XId88Pz5e3vu8B5/Tov1gBGS3TJ3WClMyqfo4PQCeSA50nvvutCpuSa00slFzhLpfaSlBMT9JcCMJqi/c6np6l8jhpCJ5ADjUZ71Lnk+m+OjbUnXr51S9rcpBoAYHRF+37GVT55fRtuhE4gB2ZnXRUzPGDo9ddpdgKQL9G+n8y4MVoIncAIiJvMPazT6E/CJoC86jbqnVHuw4nQCQyxIGzeudOeW/OnP3Uj0OOCJy+uAIok7nWPUe7Di9AJDJngE/r0tHT1anvKo8CLFy6I8iIKAIcxyn14JR46jTFL/j9/LumspCvW2uXQ8WuSmpIqkjastY+TvgZgVIU/oUtuiqNw4Az89rfZXhcAjIroKPftbenLX3avp+PjbtlfgudgpFHpnJR0Q9JtuXB5KThgjLknaSUImsaYB+HjQNGtrR2ubErt6Y8CX/hCttcFAKMi3Ndzelr69rfbs3vs7NBSNEjlFL7ntqQpSVPW2testc3QsblIZbNpjJlL4RqAkeN50ocfxgfOr3/dfUIvldzjwsJgrhEARkG16hbB2Nw8OJ0cBiuVPp3W2u3oPj9cNiO7t+UqnRtpXAcwSuLm2pRcyLx2zf1HnyQA6F24qV1yq7VFP7Qz0j07qYROv1/nU0kXJNX96uZkzKmb/jlA4UX7IcVN6M4LIgD0rlp1gbLTlHOMdM9WGqFzI9Skft8Y88QY84bcoKKebG1taX5+/uV2rVZTrVZL+DKBIeF/zK5OT+vh22fU0P/T7MIf8MIHAAnoNp0cqxolo16vq16vh3e9Ende4qEz0odTck3ofyZX+YyajvseU1NTWl9fT/rSgOETfMze2ZH291Utl1WdmJAWHkrilQ8A0hRtYWJVo+OJFgeNMZ/HnZfoQCJjTMUYsxXZ3ZT0mlz4jGtij4ZUIL88T1pZcY9S+2P2/r7b3t8/OLEcACA1wUj3732v3bQefZlGctJoXl+ObE9KemKt3TDGRJvYK3JTKwH5F9d5KPiY7Vc6VS7zcRsAMhRufqePZ7oSrXT6Tesvq5nGmElJFWvtqr9rwxgzE/qSirWWkesohk6dhx4+lP72b6Xbt90jr3IAMBCdVjNCMtKodK76qw5Jrlk9PPn7oqTrxpiK3Kj1xRSeHxhOnToPsWg6AAwF+nimK42BRNuSbnY5FjS/30/6uYGh4nmH5+kIlslgQjgAGDq8TKcrlXk6gcLzPPeKFcxIfPdue8FfXsUAYGj18jLNhPLHQ+gE0tBoSM+ft7eZAA4AcoHBRseXxtrrAGZn3XprAToHAUAuMNjo+Kh0AkmJtrd0W3sNADCSGGx0fIRO4KSCAUN370ovXhxsbyFoAkCuMNjo+AidwEkEnXuePZNaLbeP/psAkGvUFI6HPp3ASaytHQycpRLtLQAAxCB0Asflea5JPQicp05JV64wlBEAgBg0rwPH1Wi4PpySq3D++Z9LP/nJQC8JAIBhRaUTOK5gCOPYmHTmjBuhDgAAYlHpBI6LIYwAAPSM0AmcBEMYAQDoCc3rAAAASB2hEwAAIAOeJ62suMcionkdAAAgZcFaIsHymUWcXY9KJxBV9I+iAIDENRoucO7ttReuKxoqnYDntUegS3wUBQAkLphlL3h7KeLCdYROFFu0vePttw9/FCV0AgBOiFn2CJ0oknBFs1p12++9J+3sSPv7bg313/6Wj6IAgFQUfZY9QieKIVrRvHVLunq1HTglt4b6v/6r9IMfSJubxf0oCgBACgidKIZoD+6PPnKP+/tu3fRWy5334oULnNevD/RyAQDIG0avoxjC66SPj0tvvtnePn1amphoH6NJHQCAxFHpRDEEPbjX1tz2668f7NEtFbt3NwAAKSN0olh+9jPXrP6zn7nQGW5GJ2wCAJAamtdRHMzMCwDAwKRa6TTGzEmatNbeD+27JqkpqSJpw1r7OM1rAF6anpbKZTdoiL6bAABkKu1K5w1JZ4MNY8w9uaB531p70z8OpM/z3BRJe3sueN66RXM6AAAZSq3S6Vc5m5Hdc9bat0LbTWPMnLV2I63rQAFFl7VsNKRf/9pN/t5quSmSNjcHeIEAABRPms3rk5KeBhsdQui2pEuSCJ1IRngS+LExFzCfP29PAC9Jp07RtA4AGBnRBfVGVSqh0xhz2Vp73xhzKbR7MubUTUkX0rgGFFR4sFB4paFAqSS9885o/9UCAAojuqDew4ej+xaWeJ9OY8ykXAUz6mzMPiBZwWChctlN+j4+7v4tucczZ6SFhcFeIwAAPcrTxCtpVDr/zFq7GrP/acy+6bhvsLW1pfn5+ZfbtVpNtVotoctDbkUHC/3wh24S+EbDhVHWUwcAjJhgQb2g0jmMvcPq9brq9Xp41ytx5yUaOo0xFUmfdji8rfgm9mg/T01NTWl9fT3JS0MRBB8Hg/XUNzddwCRkAgBGVLCg3jD36YwWB40xn8edl3Slc0ZSxR80JEnnJZ01xshau2qMiTaxVyTdTvgaUFSj8HEQAIA+5aV+kmjoDE8CL0nGmAuSHoSa2zeMMTOhCeErTJeExIzCx0EAAAoqzXk6r0mak6t8PvUD6aKk634z/AV/G+hfdP6I8HZ4PXUAADAUUgud/opDNyP7tiUt+5v3D30R0AvPc+Hy+XM3Qv2HP3QDiPIwnwQAADmV9jKYwMl5nrSy4h4laW3NBcxWyz3euZOf+SQAAMipNFckAk4ublbcqFdflX7zGwYQAQAwxAidGG5xs+J+6UsHz/na16Rr1xhABADAECN0YrjFTYPUaLjJ3/f33SPzcQIAMPQInRhunaZBmpigOR0AgBFC6MTwi1YxmY8TAICRQ+jE8InOwRmH5nQAAEYKoRPDJW60OuESAIC+9VLDyRKhE8MlbrT6MPylAAAwQoaxhsPk8BguwWj1sTEGCQEAcExxNZxBo9KJ4RG0A9y65aZBGpb2AAAARkzcjIODRujEcBjGdgAAAEbUME70QujEcGg0pJ0dN+H7zg59OQEAOKFhm+iFPp0YDtPTLnBK7nF6erDXAwAAEkWlE4MV9OP87LPDS1sCAIDcIHRicJaXpfffd0FzbMz9VyoNT49nAACQGEInBuNP/kT65S/b23t77nFxUVpYGK5OKAAA4MTo04nsLS8fDJyBVks6d47ACQBADhE6kb319fj9ExM0qwMAkFOETmTL86RXXz24r1SS/uIvmJsTAIAco08nshOeAH5szE2L9Id/KH3/+4RNAAByjkonshNeCFaSrl6VfvUrAicAAAVA6ER2goVgx8aYFgkAgIKheR3ZGcaFYAEAQCYIncjWsC0ECwAAMkHzOgAAAFKXeKXTGDMpaUnStqRLkm5bazdCx69JakqqSNqw1j5O+hoAAAAwXNJoXr9urV2WJGPMhqQnxpgpa+22MeaepJUgaBpjHsgFUwAAAORYGs3rS8aYOUmy1jb9fRX/cS5S2WwG5wIAACC/0gidbwTN6caYIGwG4bIZOTdoggcAAECOJR46Q9VNSboiadlauy1pMub0TbWroAAAAMipVKZM8iucl+UC5Yq/+2waz4Uh5XnMxwkAAF5KJXT61c6bfvh8ZIx5Q9LTmFOn475+a2tL8/PzL7drtZpqtVoal4o0hNdYHx93E8ITPAEAyKV6va56vR7e9UrcealMmeQ3p8ta2zTGbEu6LumB4pvYo/08NTU1pfX19aQvDVkJr7G+u+u2CZ0AAORStDhojPk87rxE+3T6g4W2Yg5N+oOLok3sFbkwijxhjXUAABCR9ECipqTlyL6KpHv+vzeMMTPhY+GJ4zGiPE9aWXGPUnuN9e99j6Z1AAAgKeHmdb85/bG/6tC2pDckLYaC5aKk635fzwv+NkbZ6qr07ruuKf30aemdd6SFBdZYBwAAB5Rardagr+GQ+fn5Fn06R4DnSV/5ivTiRXtfqSSdOUOFEwCAgjLGPLLWno/uT2NyeBRFo+EqnGGtVnvwEAAAgC+VKZNQENvbLmQGxsbcI4OHAABABKETx+N50j/8Q3u7VJIWF6Vz55gQHgAAHELoxPFEm9bHxtoDiAAAACLo04njmZ2VJiakctmNWv/RjwicAACgIyqdOJ5gLk7WVwcAAD0gdOL4mIsTAAD0iOZ19C668hAAAECPqHSiN54nffWrbg7O8XHp44+pcgIAgJ5R6URv1taknR03L+fOjtsGAADoEaETAAAAqSN0ojcLC65ZvVRyjwsLg74iAAAwQujTid5Uq256JKZIAgAAx0DoRGeedzBkMkUSAAA4JkIn4nmedPGiGzRULrsVh5aWBn1VAABgRNGnE/EaDRc49/elFy+kd99lfk4AAHBshE4c5nnSZ5+5QUOBvT0XRAEAAI6B5nUcFDSr7+660Dk25ubmnJhwfTsBAACOgdCJgxoN6dkzFzRLJenKFencOUasAwCAEyF0os3zpH/7Nxc4Jff4e78nXb8+2OsCAAAjj9AJJ2hW/93vDu7/r/8azPUAAIBcYSARnEbD9eOMevPNzC8FAADkD5VOONPTbj7OVss9zsxI3/wmc3MCAIBEEDrhmtavXnXTIjERPAAASAHN62g3re/vu0rn5uagrwgAAOQMoRNuOqTxcTcn5/g483ECAIDEJd68boyZlBS0zV6QtGKtfRw6fk1SU1JF0kb4GAakWpUePnQVT+bjBAAAKUijT+cNa+0VSTLGVCQ9Msa8Ya1tGmPuKRRCjTEPJF1K4RrQr2qVsAkAAFKTaPO6HzKfBNvW2qZcVfOyv2suUtlsGmPmkrwGAAAADJ+k+3ROSroRs3/aD5fNyP5tUekEAADIvURDp1/FfCOye0bSA7lAGrUp17cTAAAAOZb46PXIoKElucFCG5LOJv1cAAAAGA2pTQ7vj2J/y1obNJ8/jTltOu5rt7a2ND8//3K7VqupVqslf5EAAAA4kXq9rnq9Ht71Stx5aa5IdEPSW6HtbcU3sUf7eWpqakrr6+tpXRcAAAASEi0OGmM+jzsvlcnh/bk4b1hrt/3tmQ5N7BW5/p7IiudJKyvuEQAAICNpTA5/WdJjSU/9JvaKpPP+vg0/gAb9Pit+GEUWPE+6eNEteTk+7iaEZ25OAACQgURDpz9P572YQ0G/zkVJ1/3zLvjbyMramvTsmVtffXfXrUBE6AQAABlINHT6k8GXuhzflrTsb95P8rlxBM+T7t51gVNy66yzxjoAAMhIKn06MYQaDenFC/fvUkn6xjeocgIAgMykOXodw8LzpM8+c9VNyfXnXFgY7DUBAIBCIXTm3eqq9O1vS/v70qlT0uKiC5xUOQEAQIYInXnleW7g0OqqC5yS9Py5eyRwAgCAjBE68yiYGul3vxv0lQAAAEhiIFE+NRpuaqSo06fpywkAAAaC0JlH09PtqZECf/RH0q9+RdM6AAAYCEJnHm1uummRAqdPS7duETgBAMDAEDrzaHraTY9UKrkR6x98QOAEAAADRejMG8+Trl5tT5H0ox9JS0uDvioAAFBwjF7Pm0ZD2tlpT5O0uTnQywEAAJCodObP9HQ7cO7vu20AAIABI3TmzeamVPb/t5bLVDoBAMBQoHk9LzzPNa1PT0sTE9LurltjfXZ20FcGAABA6MyFYAWiIGjeuuUqnLOzjFoHAABDgdCZB42GC5x7e+5xc1O6fn3QVwUAAPASfTrzYHbWVTjHxmhSBwAAQ4lKZx5Uq9LDh67iSZM6AAAYQoTOvKhWCZsAAGBo0bwOAACA1BE6AQAAkDpCJwAAAFJH6AQAAEDqCJ2jyvOklRX3CAAAMOQYvT6KoisQPXzIyHUAADDUUgmdxpgbkh5Yazci+69JakqqSNqw1j5O4/lzzfOk996Tdnak/X0XPBsNQicAABhqiYZOY8ycpBlJlyU9iBy7J2klCJrGmAeSLiX5/LkXVDifPZNaLalUYgUiAAAwEhLt02mt3bDW3pSrZkbNRSqbTT+koleNRjtwSi503rpFlRMAAAy9TAYS+eEyGkS3RaWzP9PT7cApuX9vbg7uegAAAHqU1ej1yZh9m3J9O9GrzU1X3QycOkXTOgAAGAlZhc6zGT1Pvs3OSmfOSOWyC5wffEDTOgAAGAlZTZn0NGbfdEbPnR/VqpseqdFwAZTACQAARkRWoXNb8U3scQOOtLW1pfn5+ZfbtVpNtVotpUsbMdUqYRMAAAyNer2uer0e3vVK3HmZhE5r7YYxJtrEXpF0O+78qakpra+vp39hAAAAOJFocdAY83nceVkug7lhjJkJbVeik8cDAAAgn5KeHH5GUk3SnKSzxpi6P2+nJC1Kum6MqUi64G8DAACgABINnf7k748lLccc2w7tv5/k8wIAAGC4Zdm8DgAAgIIidAIAACB1hM5h5XnSyop7BAAAGHFZzdOJfniedPGitLsrjY+7CeGZmxMAAIwwKp3DqNFwgXNvzz02GoO+IgAAgBOh0jmMpqfd+uqtlqt0zs4O+ooAAABOhErnsPE86epVV+Usl6Vbt2haBwAAI4/QOWyCpvX9fVfp3Nwc9BUBAACcGKFz2MzOuib1sTGa1gEAQG7Qp3OYeJ6rdN665Sqcs7M0rQMAgFwgdA4LpkkCAAA5RvP6sGCaJAAAkGNUOoeB50mffeb6cUr05QQAALlD6By01VXp3XddhfP0aWlxUVpYoGkdAADkCqFzkDxP+s53pBcv3Pbz59K5cwROAACQO/TpHKRGw83HGRgbo1kdAADkEqFzkGZnpYkJqVRyqw/95V9S5QQAALlE6BykalX67ndd6Gy1pB/+0DW5AwAA5Ayhc5BWV6X3328vebmzw1RJAAAglwidgxIMIqJPJwAAKABC56BEBxGVy9IHH9CnEwAA5BKhc1CCQUTlsnTqlPSTn0hLS4O+KgAAgFQwT+egVKtuffVGwwVQKpwAACDHCJ2DVK0SNgEAQCHQvA4AAIDUEToBAACQOkInAAAAUpd5n05jzDVJTUkVSRvW2sdZXwMAAACylWnoNMbck7QSBE1jzANJl7K8BgAAAGQv6+b1uUhls2mMmcv4GgAAAJCxzEKnHy6bkd3bKkql0/Okb33L/ed5g74aAACATGXZvD4Zs29T0oUMr2EwPE/66lelnR23/dOfSp98whydAACgMLJsXj+b4XMNl0ajHTgl6cULaW1tYJcDAACQtSwrnU9j9k3Hnbi1taX5+fmX27VaTbVaLa3rSt/srDQ2Ju3tDfpKAAAAElWv11Wv18O7Xok7L8vQua34JvZoP09NTU1pfX09/SvKgue5Sudf/ZX0d38n7e9L4+PSwsKgrwwAAODEosVBY8zncedlFjqttRvGmGgTe0XS7ayuIXOeJ128KO3uuqD54x9Lm5uu8kl/TgAAUCBZT5m0YYyZCW1XrLUbGV9DdoK+nHt77nFzU7p+ncAJAAAKJ+sViRYlXTfGVORGrS9m/PzZ8Tzp1792zemSe5yO7cIKAACQe5mGTmvttqRlf/N+ls+dqaBZ/dmz9r5y2VU6AQAACijr5vViaDRcP85Wy22XStLEhOvLCQAAUECEzjTMzrqBQ+Wymyrp61+XHj6kLycAACgsQmcaqlXpu991/97fl37xi8FeDwAAwIAROtPgedLf/70LnK2WG7neaAz6qmf1d20AAA38SURBVAAAAAaG0Jk0z5Pee+/g6kNjY/TnBAAAhZb1lEn5Foxa39lxFc5SyQXODz6gPycAACg0Kp1JCkat7++7QUSXLkmffCItLQ36ygAAAAaK0JkUz5M++8xVN0sl6dQp18xOhRMAAIDm9USEm9X399vBEwAAAJKodCYj3Kwuuf6cL14wYh0AAMBH6EzC7KxrTg+Uy25yeEasAwAASCJ0JuM3v3GVTcmFz6UlViACAAAIIXSelOdJ777bnpdzf186d47ACQAAEELoPKlG4+BE8OUyzeoAAAARhM6Tmp2VJiZc2Dx9WvrRj6hyAgAARDBlUhLefts9LiwQOAEAAGIQOk9ieVl6/303RdKZMy50AgAA4BCa149rdVW6edMNHGq1pGfPmJcTAACgA0Lncd25c3gfA4gAAABiETqP69VXD25/+cv05wQAAOiA0Hlc1661VyE6dUr6/vcHez0AAABDjIFEx1WtSp984vpxzs5S5QQAAOiC0HkS1SphEwAAoAeEzuPwPGltzf2buTkBAACOROjsl+dJX/2qtLPjtj/80DWxEzwBAAA6YiBRvxoNaXe3vf38OfNzAgAAHCGVSqcx5oakB9bajcj+a5KakiqSNqy1j9N4/lTNzkrj4+1K5+nTzM8JAABwhERDpzFmTtKMpMuSHkSO3ZO0EgRNY8wDSZeSfP5MVKvSxx/TpxMAAKAPiYZOv7K5YYyJC5Nz1tq3QttNY8xctBo6Ehi1DgAA0JdM+nT6FdBmZPe2RrHSCQAAgL5lNZBoMmbfplzfTgAAAORcVlMmne3n5K2tLc3Pz7/crtVqqtVqiV8UAAAATqZer6ter4d3vRJ3Xlah82nMvulOJ09NTWl9fT3FyzkGz2PJSwAAgIhocdAY83nceV1DpzFmSdIbRzzXDWtttL9m1Lbim9iP+rrhsLwsvf++1GpJZ85IDx8SPAEAAPrQNXRaa1eTeBJr7YYxJtrEXpF0O4nvn6rVVenmzfb2zg4rEAEAAPQpyxWJNowxM6HtykhMl/TRR4f3MRk8AABAX5KeHH5GUk3SnKSzxpi6tTYoEy5Kum6MqUi64G8PvzfflH75y/b2X/81VU4AAIA+JT05/GNJjyUtxxzbDu2/n+TzpmppyT1+9JELoME2AAAAepbV6PXRtrRE2AQAADiBLPt0AgAAoKAInUfxPGllxT0CAADgWGhe78bzpIsXpd1daXyc+TkBAACOiUpnN42GC5x7e+6x0Rj0FQEAAIwkQmc3s7Ouwjk25h6ZnxMAAOBYaF7vplp1TeqsuQ4AAHAihM6jVKuETQAAgBOieR0AAACpI3QCAAAgdYROAAAApI7QCQAAgNQROgEAAJA6QmeA5S4BAABSw5RJkguas7PS8+fS6dNuXk6mSQIAAEgMlU5JWltzy1y2Wu5xbW3QVwQAAJArhE4AAACkjtApSQsL0sSEVCq5x4WFQV8RAABArtCn0/NcH84f/EDa3GSNdQAAgBQUO3R6nnTxouvHOT4uPXxI4AQAAEhBsZvXGw0XOPf23GOjMegrAgAAyKVih87ZWVfhHBtzj7Ozg74iAACAXCp283q16prUGw36cgIAAKSo2KFTckGTsAkAAJCqYjevAwAAIBOETgAAAKSO0AkAAIDUJdqn0xgzKWnJ37wgacVa+zh0/JqkpqSKpI3wMQAAAORX0gOJblhrr0iSMaYi6ZEx5g1rbdMYc0+hEGqMeSDpUsLPDwAAgCGUWPO6HzKfBNvW2qZcVfOyv2suUtlsGmPmknp+AAAADK8k+3ROSroRs3/aD5fNyP5tUekEAAAohMRCp1/FfCOye0bSA7lAGrUp17fzkK2tLc3Pz7/8r16vJ3WZxzLo5x8l3Kv+cL/6w/3qHfeqP9yv/nC/+pP3+1Wv1w/kNkmvxJ2X6Oj1yKChJbnBQhuSzvbzfaamprS+vv7yv1qtluRl9i3vvyxJ4l71h/vVH+5X77hX/eF+9Yf71Z+8369arXYgt0n6PO68VKZM8kexv2WtDZrPn8acNp3GcwMAAGD4dB297lcro03mUTf8QUMH9kl6K7S9rfgm9ujXAQAAIIdKrVYr0W/oz8V5PwiixpgZa+1jY8yWtXYqdN49Sbf95vfo9/g/Sf+T6IWdzCvqUCrGIdyr/nC/+sP96h33qj/cr/5wv/pTtPv1B9ba34/uTHpy+MuSHkt66jexVySd9/dtBAHUP70SFzglKe5CAQAAMLoSC53+PJ33Yg4F/ToXJV33z7vgbw+cMeaGpAfhAMzKSvHi7lXk+JykSWvt/dC+Qt4rqfP98v8GLst1O5G1djV0jPsV/7cYdNF5HDle2PsFpOGo97/Qebzeo2+JhU6/Ob3U5fi2pGV/836n87Li/8HMyL35P4gcZmWlkCPuVdgNSbdDX1e4eyV1v1/+79MNa+1b/vYjY8ynfhcU7tfh368la+3N0Lk3/Pu1XdT7JZ3sg3HRgkG3e0WBIVbH97/oeeL1vpffIYoMIamMXh8F1toN/83swB8SKysd1ulehXVYAKBw90o68n7dVuiFWtLF0D3ifh0WfdN6ovb8voW8X74b1tqb/n1blvTQf+0K3vw3rLX3/eMvF+3odizHOt6rbseKeK96eP8LzuP1vq3b71AldHxV0hVjzIx/rHC/X1KBQ2cXrKx0PJMKTY3FvTrM/0Q8F24e9lsAuF+dnfWb3QOX/KpwYe/XCT8YFyoYdLtXFBhidXz/izmv8K/3PfwOUWSIIHRGJLmyUlEYYy6H+/X4uFeHVSRtG2PmjDGXjTHXQi8y3K94i5KW/G4I19TuolPk+3WsD8YFDQbdQhQFhogj3v8k8Xof0e13iCJDDEJnjKRWVioC/w9rO+YQ9+qw4AX4abhJxf+0zP2K4f8t/lztF/fgHhb2fp3gg3HhgkG3e0WBIV6X9z9e7yOO+B2iyBCD0NkFKyv15M86jGbnXh22LTfa80CTiqQr4n7FMsbclusT9ZqkVUkP/D5Rhb5fx/xgXORgIOlwiKLA0FnM+5/E6/0hXX6HKDLEIHR2x8pKXfh/PJ92OMy9Oqypw1WCYOQi9yvCD5dPglGz/ojaZbmQzv1S3x+MCxsMpI4hqtOxQt8r34H3P17vu4v5HaLIEIPQ2YHff+xGqA/GTIdPwBV1n0Yoz2YkzfnNBtfkFgK4ZIxZ4l4d5oen6AvzpKQm9ytWRYfftFYlN+Jd3C+pvw/GRQ8G0XvV7Vih71Xc+594vT9K9HeIIkMMQmcME1lZyf+DO+8f3gimPPB1XFkp74Img9B0EU25vlLBPGTcq8NuRkYonld7dCP366ANSbXIvjlxvyT1/8G4yMGgQ4jqeKzg9yr2/Y/X+846/A5RZIiR6DKYo8T/46jJvYmdNcbUrbU3zYiurJSmTvcqcs41/3jFGPPUH91YuHsldb9f1tpl4yY4r0h6TdJiaNJl7lfoflk3AfyKP2VSMC1JM9RcVcj7JZ1oyeGelyPOi2736gT3MZd6eP8LzuP13nfE79BNY0x4BPt5tauhhfv9kqRSq9Ua9DUAAHoUnRsw5JK1dsN/47su6T/k3vzrwRtbt2N51O1eyVXqjnUfAenov0X/nOBD82sq8N9igNAJAACA1NGnEwAAAKkjdAIAACB1hE4AAACkjtAJAACA1BE6AQAAkDpCJwBkzJ/bL27/pD+P6420nwsAssaUSQAKJzRH3jW5lY+ClUCm5VYNuZ3WnHn+ilTN0KIA0eMVSY+stVMJPd+MJBVhDkAAw43QCaCwjDGPJC2HVwLxA+k/ygXPxFcI8deqXj3inCfW2tcSfM5r0VXEACBrNK8DQIi/fvKipHt+AE2Mv3zgIJa6e0wzO4BBI3QCQIQfPH8u1wSfpNc6Naunya/YXjryRABI0alBXwAADKkHCoVOv+q5JOmxpBlJG+F+kn4VsynprKRPo30o/a/fjnsi/2uD85/GHO/43KFjTUmV0DXoqGZ8AMgSoRMA4m3LBbzAPUlv+VXQDb8/6BuSZIy5J2klFAQvqx0iA+clbUafJOZrZ+SHxl6e2z92xVrb9J+30qn/pjFmhgFFAAaF5nUAiPey4uiPOA+a3QOfGmMu+6PNo2Huhr8/LKhCvhT3tTEV0o7P7f97LtRk35RU6/DzPNHhMAsAmaHSCQDxKmpXK2ckbQcB0PdE7WpotHp4Kabv5qQiodP/2qP6eHZ7bulwk/2h5vnQeYkOjAKAfhA6ASDeBUmfhrabkSmUNiQ3BZIiYS6FwUKxz+37uTFmzj8+J+l2ws8NAImgeR0A4l2WtOz/OxjAc4A/iGdDrip6lLhK4+Mevrbbc0tuwFMl6Edqrb3f4ft0HMgEAFkgdAJAhD+4ZznoRxlUGcP9NEMr/TTl5sGcCx2rxPTpDEaXvxT62vD3nYuc0/G5fRestavW2vtHTGb/mo5uygeA1LAiEYDC6bIMpuTCWewymP6a6E/k+k02I1MmhY9txwVAY8wNa+1yzP7olEmPJN231r511HP7zfvLalcxP5V0I9rEb4y5ba290vGmAEDKCJ0AkJGkg59f/bwcTJHkh+nzcqHzjci5hE4AA0XzOgBk50lMs/tJLEt6OQG8tTaosNbDTfB+k/2DmK8HgMwQOgEgI35FMsk10B/JjViPivbfnOkywAgAMkHzOgBkKDTKPJFBPf73OysXMif9f28E3z804ImViAAMFKETADJmjLmcVeUxNIcnAAwUoRMAAACpo08nAAAAUkfoBAAAQOoInQAAAEgdoRMAAACpI3QCAAAgdYROAAAApO7/A+F0zA514xBPAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 864x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "figsize(12,6)\n",
    "plot(streaml.ra(),streaml.dec(),'r.')\n",
    "plot(streamt.ra(),streamt.dec(),'b.')\n",
    "xlabel(r'$\\mathrm{RA}\\,(\\mathrm{deg})$')\n",
    "xlabel(r'$\\mathrm{Dec}\\,(\\mathrm{deg})$')\n",
    "gca().set_aspect('equal');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we compute the conversion matrix, first defining a function that computes the relative angular momentum between two orbits:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def relative_angmom(orb,ref_orb):\n",
    "    return numpy.array([(orb.y()-ref_orb.y())*(orb.vz()-ref_orb.vz())\n",
    "                        -(orb.z()-ref_orb.z())*(orb.vy()-ref_orb.vy()),\n",
    "                        (orb.z()-ref_orb.z())*(orb.vx()-ref_orb.vx())\n",
    "                        -(orb.x()-ref_orb.x())*(orb.vz()-ref_orb.vz()),\n",
    "                        (orb.x()-ref_orb.x())*(orb.vy()-ref_orb.vy())\n",
    "                        -(orb.y()-ref_orb.y())*(orb.vx()-ref_orb.vx())])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "and we use that the compute the relative angular momentum of the stream progenitor wrt a non-moving Sun (``Orbit()`` is the Sun). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "sun= Orbit([0,0,0,0,0,0],radec=True)\n",
    "sun.vxvv[0][1]= 0.\n",
    "sun.vxvv[0][2]= 0.\n",
    "sun.vxvv[0][4]= 0.\n",
    "sun.turn_physical_off()\n",
    "o.turn_physical_off()\n",
    "angmom_vec= relative_angmom(o,sun)\n",
    "rot= _rotate_to_arbitrary_vector(numpy.atleast_2d(angmom_vec),[0.,0.,1.])[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can then use this ``rot`` matrix to convert sky coordinates as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "galpyWarning: Method ra(.) requires ro to be given at Orbit initialization or at method evaluation; using default ro which is 8.000000 kpc\n",
      "galpyWarning: Method dec(.) requires ro to be given at Orbit initialization or at method evaluation; using default ro which is 8.000000 kpc\n"
     ]
    }
   ],
   "source": [
    "phi12l= coords.radec_to_custom(streaml.ra().to_value(u.deg),streaml.dec().to_value(u.deg),T=rot.T,degree=True)\n",
    "phi12t= coords.radec_to_custom(streamt.ra().to_value(u.deg),streamt.dec().to_value(u.deg),T=rot.T,degree=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "and the resulting stream in $(\\phi_1,\\phi_2)$ is:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAtcAAAB1CAYAAACFxc2RAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAfaklEQVR4nO3dX2wcV/UH8O/uxg68gJMAqlqpv3aD9jyhX+P8QX6qQxwqKrWA09ZvQRTiVG1BSPx+iSJUKVIFxgFEBBXFTmlF3tw0oQj9qGic1oBUSyV/KqEgHQm7gBCqADduX2jiNvt7uPd6796dmZ1dz3q9u9+PFNm7M7szvtnZOXv2zLm5crkMIiIiIiJau3y7d4CIiIiIqFswuCYiIiIiygiDayIiIiKijDC4JiIiIiLKCINrIiIiIqKMMLgmIiIiIsrIpnbvQJY+/elPl2+77ba27sO1a9ewZcuWtu5DN+A4ZoPjmA2OYzY4jtngOGaD45iNXh3Hq1ev/ltVPx61rKuC69tuuw3nzp1r6z6Mjo62fR+6AccxGxzHbHAcs1E1jvPzwNwcMDwMDA21c7c6Dl+P2eA4ZqNXx1FE/hq3rKuCayIi6gDz88C+fcCNG0B/P3DhAgNsIuoarLkmIupg8/PAxIT52THm5kxg/cEH5ufcXOu2lTRAHTl4RLTRMXOdsbGxsXbvQlfgOGaD45iN9RrHNJUSbp1t24ArV4BnnzUx6qZNwJe/DBw8mJAEbmcpxvw8jm/eDExPA3/7m9lhwGSuh4frPrap/U7KkMctC7cVDvhbbwG33FJnoFuLx3U2OI7Z4DjWypXL5XbvQ2ZGR0fLvVj3Q0Sdy4/dvvENE+sVCsDDD9fGb/PzwN69wPXr0c+VywEf+lBMlUW9QLMVQXf4x12/Dty8CeTzJrh2fyQQvf35eeD06coniP5+4ORJYGkpel/97S0tmSD+1Cnz2EIBePJJ4Ngxs+7EBPDEE9XLhocrY1QoAPfeC/z618D775v99vX1AU89Fb8vRNTVROSSqu6KWpY6cy0idwHYCmDA3rUMYFFV/7LmPSQi6hBJcWgY28XFikAlpnSxHGDit3LZxHtTU8DPf14dA58+HR9YA+axrspiaCjY2ahSjKEhk0l+7DGz8c2bazO4LmPrdjptEOkH8/m82a4LUG/eNLdvv716EPzg2QXk771n/jDA/PGPPWZux2WiXQCfy5nt5m31Y5ghHx4297nthmP0wQfAiy/G/30rK8Cjj1aeO6u6cf9FsmNHc2PfyLZ4USlR5hKDaxE5AGA/gDKAHIAFmKAaALYDeEhEYJefV1WmjYmoo4UBsB+7nT4NPPecSWTWi+3yeROr+rHi179eCY6ffdYkbt36oZpAOYVczoshp6eBxx83QaLbkTCYnJ8367z/vnmC69cr9c/+H+M89xzw6qvpdsgPVMtls3O+QqE2oPWD53ze7JcLrHM58xgXpPuDMz8PHD9evb/uUwoAHDpUG5wODZn/wDC47O+vDuiTuOdv9D8qLqidnwfuvtsE7qG4sY97wTrT08DZs8CBA8D4eGX9t94yWXn3jQAvKiXKTGRwLSI7AIwAuKSqj6R5IhHZJyL/A2BWVd/IagdFZNz++jxM5vywqh7N6vmJiJwwtnnmGeB3vzO/79tXHXOF8ZSLEf3k7PXrJnZ1wbYfM62smPgmKrAuFMzPMNl68KCJlfzH5PMm7qyquca8CVL9oHlpqTaYnJioBIjuyfyAN9y5RoLIMDN8zz3AL39ZCbQffrg6oA0z3OVydYD7+c8Dn/tcpXbG/4AQZqz9x5XLJkMetc9DQ9EBt/8patMmYPt24E9/qn5soVAJ9tPUjTtJ5TmnT0cH1kD02Id1Qs8+W73O9DRw+LD5/eWXgYUF8yHLfVUSPjfQ2kw2M+XUI2qCa5utXlDV7zXyRKp6AcAFEdkhIp9R1Vcy2scBAJMApgAswmTSiYgSuYTdXXcBAwPpzudhbPP+++a+22838YefRI2rMvAz12Gs6GIxwJTs3nKLWcfFsK5mOq6seGgIePppU41w86Z5jh//OGLdibnqwNhlicNgcnjYZLWvXzfrPPVUdcAbZq4bCSLDzDAA/OY3laDS1cX46/mF5/l8JcDO54E9e0zm9VOfqv2A4D4I5PPAyIj5T//hDytZ+7T77PZnaMjsn7/ve/ea7eRywP33A0eOmPsbDRbjynPqiRp791zOykr18509W73+uXO1gbV77m3bmmuP6AfMbp/Ci0HdsuFhs499fY1l+ok6TE1wrapno1ZMS1WvrOXxEZYBbLHPvVxnXSLqYlGJr6hvxcOEXeKFfp633oq+f3jYJDBdYJzPmwA4rsrA1Vz7saIrJ/ZLaAFTV510EWMoKr6M3GEXNOfz1UGzL640Ivxjmq37DYP5qG2F67k/Lhw8F6BFfUDwM+THj5vlX/jC2rKk4XZefTV+3xsRVevtHDxoss8rK+b/7ZvfBN59t7Is3Jb/iQ4wQav/fAcOmAPAGR0FfvCD6m8r9uwxL8xmgn4/C18omAPN1UydPGnqoG7cMAfPJz9ZCexv3ABOnDAfUOpdwMAsN3Wicrlc91+pVPpqinXuKpVKH0nzfI38K5VK42nX/eIXv1gmos7z2mvl8ne+Y34mrfPhD5fLhYL5+dpr5t/mzS69WS7395v7PvvZyn3uX6FgtpHkkUeqH5PPV/bpkUfK5Vwu/XOl/dvS/O1NadkTr6O0f0On/a1J+9vo3/Laa+bF+cgj0Y+ZmjIHxNRU5XahYF7MmzdXHhN1gNXzne+Y9QHzfP4BsmdP7UEYHlybN1e2NzVV+bub2Rd/PDrptUAdq1QqXSzHxKOpWvGJyCFVPRWz7E6YeujnARQBTGbZQcSruX4bwG4AM6p6OWrdvXv3lv357cfGxth/kWgDifoGeds2U5e8smKSX/fdV91G2D0mqqsaAHzrW9XlGt/+tnlOl7l296fJXLsSVleV8JOfmEyxW8ZJBakrJF1Q2Ui2OClzfc89yd1W3AWuruzHlQH19wNf+lJ8C8Wkvyf8toMHKWVoZmYGMzMzq7evXr36V1W9I2rdtMH1AZjAeQzAEoAp1xlERJ4G8IKtuYaIfFVVn1nrH+Ftu6iqi97tBQA7o0pE2OeaaOMKO7O5OmQg+qK+/n5TT+z3fvbP3RcumPX867n6+6s7zDVac+32k99UE6UUV3MNVGqswwC6XDalIu6A9i9OKBRMdxdXL1UvSPbfWHI58xzuwtbdu2vrt4gykkWf6yLMxYSH7O0R76LF3QD87h3Xmt7TCH5gbS0DeAjAdJbbIaLmJF3P5AtbCNezsmKCY/cYwJxzb7+9+vlffTW6E9n4eCXr3Iiw1DbtMqKeFNVxxZmbq2387u4Pv77yM84HD1ZfUJp00PlvLK51DmAC+Ndfr3S/4YFL6yhtcL0Q9LC+IiKj9veyqr7rLUsMrm2Zx84625tU1UURKcK0A9ziLVuE6bFNROvMBdJXrwK//S3wsY+ZDmVu+m3XWjgq2eSuvQpbCPvdMnx9feZ6rN//vvqcG54jGfASbVBxB2dUMB51lW6zLR/9cpSwg0oa9TIG4VdY9XqNU89JG1xvF5E7XC21iHwGpuc0YDt5eIoAYtvwqWqjGeewp/UAzGQ2RJShuPPFW2+ZGugdO0xy6T//qTzm73+v/O7P3RHVbCBsIexqrJ96yrTf/f73K98e33efaSQwNJSiMwYRdb5mPyVHtXx86aX4Dir1+GUmURkDoHZG0XB2KGbKe17a4HoawBkRGYHJTD8PYFZEvgtg0k4eMw0z8UxmZSE2e+2mW4f9vdhEgE5EEfzrgFzXrP5+4Ec/Ar72teqWuIVC/KR1ru+zfx6KOp9FtRB256C4zmnMTBNRoqi2ic1mkv0yk6iMAVDdstDVrjkrK2bb/hT2UU3rqauluqCxHtsx5CiAi1lezGifewCAq5zcDlsyErUuL2ikXtboN5N+ggaoroP+xCeAf/2rNpju66uekRowmebx8UrfZmaZiahjrTVz7S7U9GejyudN33l2L+kqWVzQ6J5oFMCyqr4iIh+FrbdW1TcBpJomvVG2K8iJVjw3UafyM85ushL//f3UqUobubC8wwXefoLGXQPk/POftdt03TuWlio118Ui8N3vrm1ODSKiDSOqzCTMGISTIX3qU5XMBgBMTVU/582bjc3GCUTX6TFz0THStuLbAeAUgIsA/qyq3/fuL6vqGy3dy5SYuaZe4BIr/jTbuVxtB46+PlPPHJZ3bN5svjUFqhM0Kyu1Fxbu2QPcemt1UE5ERDFcs3yX6QAaz1xPTwOPPmrekPv6qnuS+hl0BtttlUXmesQ9gYjsc3eqqusasiGCa6Ju4men3ezTO3aYEj8XWAPRnTYAE2yfPVv97SRQSaAcO1adgPnjH837uQvSN29mi1giooYMDVXXfDdacz0/X/1GfOMG8LOfVb5mfO89M3X8b34T3wecWe+2SxtcX0pYtjVhGRE1ILzA0E9+OLmcKQMMf/rLczkTHB84UCn/cPyLDf3rgMJvN5mpJiJqwlquwp6bq82Y3Hor8MYbJrgul4Ff/cr8vHnTnCSOHzdv9q5G0M9ynzzJrHcbpA2ud4rIRdvPevU0LiJ3APhkK3aMqBvVm/0v7gJDn2tXNzJi3k+vXDHdn9w1NydPVidKXMAc1lxHYWcOIqI2Gh42pSDuRNDXZ/qS3nKLqeV2M1zm82b5zZvA+fPAyy9XZsF0s1S6rLebXODGDZP1/r//MycMv1SF2e1MNdKK7wURKQNYFpH9MBPBlFX1npbtHVEX8YPnQgF4+OHqQNe/wDCJy0ofP155bNJkZgyYiYg6hLvaPOorRH9K+JMnTd3f7GxtjWChYE4UhYLJvrivNvN5k/V2J5nr1yvtBV2deD4P3H9/ZaIBakpDrfhsvfWgvXlZVS+0ZK+axAsaaSObmACeeKLyvpbLAR/6UHXiwL/A0E0Tns9XkhWbNgFf+QpLNoiIek5ULXV4dbu7WGZpCfjb30zrKNcSavdu4OLFShDe12faPp0+Dfz0p9Xb6u/nZDh1ZNaKzwbTGyqgJsrKWr4VS/PYcPrvcCbDpA5Q/u98ryMi6kHh15D+ScP1ZQ2nafez3V/5irly3WWon3rKrOu3EXTCaeOnp02m/MAB0+OVEjUUXEcRkRlVHctiZ+ps5wiARZjp1WdV9XKrt0m9w88aR118HbW+H/imeWw4/ff779fOZBj13hn1OxERUWLdX5ixcRfhhJmagweBZ54xJyXHnzZ+eho4fNj8/vLL5icD7EQ1wbWd0nwcwFKwKK4ryEDM/ZkRkTMAJlxALSLnAexv9Xapd/j1zvV6/YeB+Je+lP6xSdN/ExERZSoqYxN1Yc7vfhd/5fvZs9Xrnz2bHFzPz5sLJ//xD5Mt78FAPC5zvdPOughgdbIYqOoVfyVbg/1263Zv1YiqPujdXhSREVWdXYdtUw9wJRsuYPazyaEwEAfSP9bhRYZERLRhJJ2UDhyoZKzd7Tjz88Ddd1cmWHj9dfOzxwLsqOB6QlXfCe7boqqvhCuq6gUR+Uxrds0QkRGYchDfMkzmmsE1NcxNBw5UPpxHfXsWJwzEDx5kJpqIiLqUC4zT1FzPzdXOXFYv092FaoLriMC63aLKTpYA7F7vHaHO4wJp903Xjh3A449Xjv3nnjOTabkAO01gHBeIM6gmIqKuND6eLkB2fbr9ANvPdPdIP+1GJpF5W1Wrpjm3k8jsBFCT1c5Q6hkgr127htHR0dXbY2NjGBtr+bWW1GZxM71GzXLo+us79Wqk47Csg4iIKDA0ZNr7RdVcN9o5INTmwHxmZgYzMzP+XR+LWzdVcK2q3xOR523ttSvRKML0um519BpV070tasUtW7aAfa57h8tKu9kJCwXg3nuBl14yFz3n87XfToWzyubz6WqkiYiIKIWhIeAXv6i9v5HOAaHpafO1cziz5DoKE7Yi8u+4dVO34lPVh0TkTlRPIvNm0mMysozo0pCwDpu6TJqpwl3PaMAccy++WFmnXDYBtz/jYV+f+bmyYpb95CfMQBMREbVcI50DfPPzwGOPVVoFupklN/C07Y1OIvMmgPUIqP1tzopIWBpSBDC1nvtB68Mv6fjGN+K/PXIfgOMmGHVThJ88aWZ/9bsLucdvsGORiIioezXSOcA3N1f9tXOhYB6/1jKTForqc30AwCVV/UszT2iz2ztUNcv6jFkRGfQmjimyDV/38Y+TfN5knG/ejP72yP8AvGkT8LnPVcpBCgXg4YeTpwjfIMcfERFR72jmgqXhYZMtC2eWnJhovsykxaK6hZwVkX22Bd5s2iDbXtz4IICFjANrADgE4JiIFGG6hBzK+PlpnUV9k+OXY5XL5hjK5aK/PYr6ALxBvx0iIiKiZsVlvJstM1kHuXLc9+oAROR/YUowlmHa3y3a3wFTB12EuVpyAMCfVfX7Ld3bOkZHR8u8oHHjqlfyEX7Dc/IksLTEYJmIiIgitDGrJiKXVHVX1LLEmmtV/Z59go8C2AUTTG+3i5cBXFHVCxnuK3Wp6Wng0UdNmUc+bzLTYclHs+VYRERE1IM2aF/ctK343gFwwf4jShTVd/rRRytdO1zbvEKh9pucDXqcEBEREaXSULcQonqiLt6dm6tuhwcA990H7NnDDDURERF1FwbXtCZhljqqR7x/oS9gunscOcKgmoiIiLoPg2tKlGYiFz9LHXXx7tAQ8OqrZjZFILlFHhEREVEnazi4FpGPwPSZfsO77y7/NnWHev3Zo7LUx45FX5TIWmoiIiLqBflGVrYTzCwDeEVElkTkv+2id0Xkg4SHUgeKCp59LksdXpg4NGSCbAbTRERE1GsaCq4BjKhqXlW3wkzm8i2btV7EOk+LTq0XFzw7rnXek09uqFlHiYiIiNqm0bIQN/04bED9kIgcEpEcgPjZaKgjpek7zXIPIiIioopGg+u3bWnIQwAOqeq7qnpKRPYB2Jb97hkiMm5/fR7AVgCHVfVoq7ZHFQyeiYiIiNJrqCxEVc/CTIH+vKq+691/AcC+jPfNNwBgCsA1AOft71TH/DwwMWF+EhEREVHrpc5c2y4hI/bmpXC5ql7JaqciLAPYYrez3MLtdI16nT6IiIiIKHupgmsRuRPAGZis9QCAXSJyHrY0pIX7t4pBdWOiOn0wuCYiIiJqrcjgWkTuUNW/eHcdUNVdwTrjMAH3Pa3bvZrtvQ3TpWRGVS/XeUhPi5rMhYiIiIhaK1cu1zb5sJnqQZgOIDkAW1T1mYj19gG4M2pZlkSkaLuTuNsLAHaG2ey9e/eWt2zZsnp7bGwMY2Njrdy1DS1pdkUiIiIiSmdmZgYzMzOrt69evfpXVb0jat3I4Dpkg+gBe7MMYNHNyCgiX211cB2xP5cATKnqtH//6Oho+dy5c+u5K0RERETUY0TkUljV4aS9oPGiqr7jPeGdtiXfVgCDIrKoqq9ksK81RKQI4JKqbvHuXgSwvRXb2wiYcSYiIiLqTKmCa1V9R0RGVfWcvf0m7IyMNsh+0/4EgLJbL0NhT+sBAAsZb2NDYJcPIiIios7VyCQyb4rI0wAm3cWOInIHgF22/3VLpj9X1UURcSUpsL8Xw5KQbsEuH0RERESdK3VwrapXRGQawAsisgOm9/QiWjt5jDMtIkfs79sB7F+HbbYFu3wQERERda6Gpj+3E8XsEpGP2tvv1HlIJmxXkBPrsa12GxoypSCsuSYiIiLqPA0F1856BdWdoBUXHw4NMagmIiIi6kT5du9AJ3MXHz7xhPk5P4+qHojUPI5jNjiO2eA4ZoPjmA2OYzY4jtngONZicL0GURcf8kWWDY5jNjiO2eA4ZoPjmA2OYzY4jtngONZicL0G7uLDQoEXHxIRERFRkzXXZPDiQyIiIiLypZr+vFOIyL8A/LXNu/ExAP9u8z50A45jNjiO2eA4ZoPjmA2OYzY4jtno1XH8L1X9eNSCrgquiYiIiIjaiTXXREREREQZYc11A0RkEsB5VZ317hsAMA4zY+UAgMvB8iMwM1kWAcyq6uX13WsiiuMdvwCwG8BE1DEqIiMABlT1Be8+HttERFSDwXUK9sQ6COABAOeDxeOqesJbd1JELqrqsoicgXeyFpHz6OKp29OoF8wkBSwMZiqSxnEtY9yDJlX1MACISBHAJRHZqaqL4XoAptwNHtvVUrzmijDvn8sAoKrT3jK+Hq2UxzUTOXUEY7UfwFTaseI4ViSNI88zdZTLZf5L+a9UKp0vlUoj4X3B7fFSqTRof78WLJsKH99r/0ql0pT3e7FUKl0rlUpFe/uMG7twbJOW9eK/OuPY1Bj32j87NkeC+y5F3Ddix23cu4/HdvD3B+Pqv+aKpVLpTDDG7j2Sr8f04xi+LidLpdIAxzFyHCeDcSynGSuOY0PjyPNMwj/WXK/dVlsu4uxX1cs22x1mv9ynv55ks1cL7rbNDi7CZLQAYCT4dLtox7Hesp6SNI5rHONeMwCTkQ5ti1jvbXeDx3a1FK+5KXhZfwD7vNcgX49WinEMX18LMFlBgOMYGnd/v/ctVJqx4jhWixxHnmfqY3C9dodgXoCX7NcgR+39AxHrLqFygPei2GAmKWBhMFMjKShsaoyz3b3OYN/8dwZ3D8Ir/RKRB/w6a4vHdrWk19wAzIl29St5VV0G+CElQr0Pe0zkpLfTK19wx+UizzMNixxH8DxTF4PrNbIn6OdRebG5F+DWtu3UBlUnmEkKWBjMeJLGcQ1j3JOCGsFxmNpAv6ZwOeJhPLY9dV5zRQDLIjIiIg+IyBEvg8XXoyfFhz0mclIKrpk4DOCo/VDH80wD4saR55n6GFyvkYhMwVwUtR3ANIDzIjII72tkT/h1c89JCGaSAhYGM4GkoLDJMe5pNpB+UFX97MpDfsbVw2M7kPCacyfUt1X1BXvx96TNgvH1GEhxXDORk5KIFO2HkCLMuRngeaZhMePI80wdDK7XwAbRC+7Tne06cBTmE17cp+Tw65KeFBHMJAUsDGZixASFccs4jvEmATzobtjg72LMujy2Y0S85pZhWhhW1V/CvEfy9Rgj6rhmIqcxqrpoP8wdhekCVHXthIfnmQQx47iK55loDK7XpojaE+o0AMR8giuitpVfr6oKZpAcsDCYiReOY9IyjmMEm5WZ9GqBB2G+4hyxZQxHAOyCqcsc57GdKHzNueMXwX1F8PWYJPywx0ROA/wA0I7ZMoBj4HmmIQnj6ON5JgKD67WZBTAW3DeCypXxs/ZN0SnGfM3cU6KCmaSAhcFMtJigMHYZx7GWiDwA4DKAt0VkwI7hLlfC4P7BnBjOe/2ZeWwHYl5z7uIn3wCARb4eo8Uc10zkpGRr+q9FLBrgeSa9pHH01uF5JgaD6xREZNBepT0CUy94BFi96n3CThwzbuuOXLE/YC5AGbMX8kza2z0tLpixi5MCFgYznqRxXMMY9xRb+nEG5k3/mv13CUEQY4/3EQCH7dgCPLar1HnNnQjacO0CExCREsaRiZz0FlG52NNxxzrA80xaiePI80yyXLlcbvc+UI8Ie2N69qvqrP0K6hiAP8DM+DQTzE4WuazXJI0jzBtiU2NM1Ix6x7VdZ9Kusx08riOleH8chAmw3TqLWt3ZhuNoSWVW5WWYrhbnXUtNnmfSixvHtZzLewWDayIiIiKijLAshIiIiIgoIwyuiYiIiIgywuCaiIiIiCgjDK6JiHqI1/EkvH/Adj6abPW2iIi6GYNrIqIOYaciLtZfM/bxIzDts2rY1qJTAMabff4Ii0FLLiKirsfgmoioc7iJWZpVTHq8XRY1fXFTbPutkborEhF1EQbXREQ9wE6G046JHC6zPISIegmDayKi3rB9jVnvptiJTvav93aJiNplU7t3gIiIktkZz3YB2Goz0C80Eijbxy/HLDuCSh12TUmIfey4XWcQwGwwo904zMygRftzKwCo6nTa/SMi6ibMXBMRbVAiMigiD9iLDQdUdVpVTwA4YwPbtHYBWIp4/jMwwfKsm0obNjj2nAEwbdc5AeBUsOwFO7X0IkxN93RUYM0LG4moVzC4JiLagGwwetgGrqFFeBcK2jZ6SV1EXFbZf/4izAWSq91D/N/tOiP2fj/rfdGroR7xMuiLAMZitr+A2qCdiKgrsSyEiGhjOgXgwZhlq1lrGwAfBvAHACcS1g/LSAYj7gsNAlh2Qba1gEqJSVhqEtdpZNnfZyKibsbgmohog7ElH6tt92yW2Q+Ed8EE1FDVWbu8VcHrolcyAlR3HHleREbs8hGYPtlERD2NZSFERBuTnxUe9C4iHIepk26k80dU5vgyTLlIEncRYxWv3vs8gKItE7kcU8ICu+3ICyqJiLoNg2siog3G1jgvenXUW4HVEpAHARxq8CldNw9/G4swPahX7w/KP1wbPQTr+MH2bnsB4wtBdju0HfVLUIiIugLLQoiINqZ9AI6JyB8A7LQZa6hqwz2jbelIzeNU9UEROSIifiu+ARE5o6oPuu2JyKSILNjli94Fjgv2fnf7IoDJqKx6O3psExG1A4NrIqINyAawR23WeLlOZjiNyJps217Pl4tY52h4n6vzVtXt9rbrxX0GwM417isRUcdiWQgR0cY2mBRYe6Ui+8OyjsBCnXZ9jToKYLWftaq6DwAzfumI3afzGW6XiGhDY+aaiKiD2YC2blZbVU/Y2Rjj2vU16hJMh5DwIsbt8IJumA8HWW2TiGjDY3BNRNQ7FkWkmEX9s6pOi8gDthZ8EabsZCtMzfUysHrx41rLWYiIOkquXC63ex+IiGid2OnU41rmZb2tkQxqxYmIOgqDayIiIiKijPCCRiIiIiKijDC4JiIiIiLKCINrIiIiIqKMMLgmIiIiIsoIg2siIiIioowwuCYiIiIiygiDayIiIiKijPw/QfyC0L3gupIAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 864x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "figsize(12,6)\n",
    "plot((phi12l[:,0]+180.) % 360.,phi12l[:,1],'r.')\n",
    "plot((phi12t[:,0]+180.) % 360.,phi12t[:,1],'b.')\n",
    "gca().set_aspect('equal')\n",
    "ylim(8.,-8.)\n",
    "xlabel(r'$\\phi_1\\,(\\mathrm{deg})$')\n",
    "ylabel(r'$\\phi_2\\,(\\mathrm{deg})$');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Much better aligned!"
   ]
  }
 ],
 "metadata": {
  "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.7.3"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
 }