Rotation example

rotation_example

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Populating the interactive namespace from numpy and matplotlib\n"
     ]
    }
   ],
   "source": [
    "%pylab inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 74,
   "metadata": {},
   "outputs": [],
   "source": [
    "def Mrot(alpha_pole, delta_pole, phi1_0):\n",
    "    \"\"\"\n",
    "    Computes the rotation matrix to coordinates aligned with a pole\n",
    "    where alpha_pole, delta_pole are the poles in the original coorindates\n",
    "    and phi1_0 is the zero point of the azimuthal angle, phi_1, in the new coordinates\n",
    "\n",
    "    Critical: All angles must be in degrees!\n",
    "    \n",
    "    \"\"\"\n",
    "\n",
    "    alpha_pole *= np.pi / 180.\n",
    "    delta_pole = (90. - delta_pole) * np.pi / 180.\n",
    "    phi1_0 *= np.pi / 180.\n",
    "\n",
    "    M1 = np.array([[np.cos(alpha_pole),\n",
    "                    np.sin(alpha_pole), 0.],\n",
    "                   [-np.sin(alpha_pole),\n",
    "                    np.cos(alpha_pole), 0.], [0., 0., 1.]])\n",
    "\n",
    "    M2 = np.array([[np.cos(delta_pole), 0., -np.sin(delta_pole)], [0., 1., 0.],\n",
    "                   [np.sin(delta_pole), 0.,\n",
    "                    np.cos(delta_pole)]])\n",
    "\n",
    "    M3 = np.array([[np.cos(phi1_0), np.sin(phi1_0), 0.],\n",
    "                   [-np.sin(phi1_0), np.cos(phi1_0), 0.], [0., 0., 1.]])\n",
    "\n",
    "    return np.dot(M3, np.dot(M2, M1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Can use this manually like this"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 98,
   "metadata": {},
   "outputs": [],
   "source": [
    "def phi12_rotmat(alpha, delta, R_phi12_radec):\n",
    "    '''\n",
    "    \n",
    "    Author: Denis Erkal\n",
    "    Converts coordinates (alpha,delta) to ones defined by a rotation matrix\n",
    "    R_phi12_radec, applied on the original coordinates\n",
    "\n",
    "    Critical: All angles must be in degrees\n",
    "    '''\n",
    "\n",
    "    vec_radec = np.array([\n",
    "        np.cos(alpha * np.pi / 180.) * np.cos(delta * np.pi / 180.),\n",
    "        np.sin(alpha * np.pi / 180.) * np.cos(delta * np.pi / 180.),\n",
    "        np.sin(delta * np.pi / 180.)\n",
    "    ])\n",
    "\n",
    "    vec_phi12 = np.zeros(np.shape(vec_radec))\n",
    "\n",
    "    vec_phi12[0] = np.sum(R_phi12_radec[0][i] * vec_radec[i] for i in range(3))\n",
    "    vec_phi12[1] = np.sum(R_phi12_radec[1][i] * vec_radec[i] for i in range(3))\n",
    "    vec_phi12[2] = np.sum(R_phi12_radec[2][i] * vec_radec[i] for i in range(3))\n",
    "\n",
    "    vec_phi12 = vec_phi12.T\n",
    "\n",
    "    vec_phi12 = np.dot(R_phi12_radec, vec_radec).T\n",
    "\n",
    "    phi1 = np.arctan2(vec_phi12[:, 1], vec_phi12[:, 0]) * 180. / np.pi\n",
    "    phi2 = np.arcsin(vec_phi12[:, 2]) * 180. / np.pi\n",
    "\n",
    "    return [phi1, phi2]\n",
    "\n",
    "def pmphi12(alpha, delta, mu_alpha_cos_delta, mu_delta, R_phi12_radec):\n",
    "    '''\n",
    "    Author: Denis Erkal\n",
    "\n",
    "    Converts proper motions (mu_alpha_cos_delta,mu_delta) to those in coordinates defined by the rotation matrix, R_phi12_radec, applied to the original coordinates\n",
    "\n",
    "    Critical: All angles must be in degrees\n",
    "    '''\n",
    "\n",
    "    k_mu = 4.74047\n",
    "\n",
    "    phi1, phi2 = phi12_rotmat(alpha, delta, R_phi12_radec)\n",
    "\n",
    "    r = np.ones(len(alpha))\n",
    "\n",
    "    vec_v_radec = np.array([\n",
    "        np.zeros(len(alpha)), k_mu * mu_alpha_cos_delta * r,\n",
    "        k_mu * mu_delta * r\n",
    "    ]).T\n",
    "\n",
    "    worker = np.zeros((len(alpha), 3))\n",
    "\n",
    "    worker[:, 0] = (\n",
    "        np.cos(alpha * np.pi / 180.) * np.cos(delta * np.pi / 180.) *\n",
    "        vec_v_radec[:, 0] - np.sin(alpha * np.pi / 180.) * vec_v_radec[:, 1] -\n",
    "        np.cos(alpha * np.pi / 180.) * np.sin(\n",
    "            delta * np.pi / 180.) * vec_v_radec[:, 2])\n",
    "\n",
    "    worker[:, 1] = (\n",
    "        np.sin(alpha * np.pi / 180.) * np.cos(delta * np.pi / 180.) *\n",
    "        vec_v_radec[:, 0] + np.cos(alpha * np.pi / 180.) * vec_v_radec[:, 1] -\n",
    "        np.sin(alpha * np.pi / 180.) * np.sin(\n",
    "            delta * np.pi / 180.) * vec_v_radec[:, 2])\n",
    "\n",
    "    worker[:, 2] = (np.sin(delta * np.pi / 180.) * vec_v_radec[:, 0] +\n",
    "                    np.cos(delta * np.pi / 180.) * vec_v_radec[:, 2])\n",
    "\n",
    "    worker2 = np.zeros((len(alpha), 3))\n",
    "\n",
    "    worker2[:, 0] = np.sum(\n",
    "        R_phi12_radec[0][axis] * worker[:, axis] for axis in range(3))\n",
    "    worker2[:, 1] = np.sum(\n",
    "        R_phi12_radec[1][axis] * worker[:, axis] for axis in range(3))\n",
    "    worker2[:, 2] = np.sum(\n",
    "        R_phi12_radec[2][axis] * worker[:, axis] for axis in range(3))\n",
    "\n",
    "    worker[:, 0] = (np.cos(phi1 * np.pi / 180.) * np.cos(\n",
    "        phi2 * np.pi / 180.) * worker2[:, 0] + np.sin(\n",
    "            phi1 * np.pi / 180.) * np.cos(phi2 * np.pi / 180.) * worker2[:, 1]\n",
    "                    + np.sin(phi2 * np.pi / 180.) * worker2[:, 2])\n",
    "\n",
    "    worker[:, 1] = (-np.sin(phi1 * np.pi / 180.) * worker2[:, 0] +\n",
    "                    np.cos(phi1 * np.pi / 180.) * worker2[:, 1])\n",
    "\n",
    "    worker[:, 2] = (-np.cos(phi1 * np.pi / 180.) * np.sin(\n",
    "        phi2 * np.pi / 180.) * worker2[:, 0] - np.sin(\n",
    "            phi1 * np.pi / 180.) * np.sin(phi2 * np.pi / 180.) * worker2[:, 1]\n",
    "                    + np.cos(phi2 * np.pi / 180.) * worker2[:, 2])\n",
    "\n",
    "    mu_phi1_cos_delta = worker[:, 1] / (k_mu * r)\n",
    "    mu_phi2 = worker[:, 2] / (k_mu * r)\n",
    "\n",
    "    return mu_phi1_cos_delta, mu_phi2\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 99,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 0.5019448  -0.05098569 -0.86339555]\n",
      " [ 0.1010563   0.99488071  0.        ]\n",
      " [ 0.85897558 -0.08725156  0.50452762]]\n"
     ]
    }
   ],
   "source": [
    "Rotmat = Mrot(354.2, 30.3, 0.) ## Poles taken from Shipp+18\n",
    "print(Rotmat)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 100,
   "metadata": {},
   "outputs": [],
   "source": [
    "alpha = np.array([-6.3, 3.2])\n",
    "delta = np.array([-59.7, -59.4])\n",
    "pmracosdec = np.array([1, 2])\n",
    "pmdec = np.array([1, 2])\n",
    "\n",
    "## Convert using Rotmat\n",
    "newcoo = phi12_rotmat(alpha, delta, Rotmat)\n",
    "newPM = pmphi12(alpha, delta, pmracosdec, pmdec, Rotmat)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Or you can use astropy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 101,
   "metadata": {},
   "outputs": [],
   "source": [
    "from astropy.coordinates import frame_transform_graph\n",
    "from astropy.coordinates.matrix_utilities import rotation_matrix, matrix_product, matrix_transpose\n",
    "import astropy.coordinates as coord\n",
    "import astropy.units as u\n",
    "\n",
    "class StreamFrame(coord.BaseCoordinateFrame):\n",
    "    \"\"\"\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    representation : `BaseRepresentation` or None\n",
    "        A representation object or None to have no data (or use the other keywords)\n",
    "    Lambda : `Angle`, optional, must be keyword\n",
    "        The longitude-like angle corresponding to Sagittarius' orbit.\n",
    "    Beta : `Angle`, optional, must be keyword\n",
    "        The latitude-like angle corresponding to Sagittarius' orbit.\n",
    "    distance : `Quantity`, optional, must be keyword\n",
    "        The Distance for this object along the line-of-sight.\n",
    "    pm_Lambda_cosBeta : :class:`~astropy.units.Quantity`, optional, must be keyword\n",
    "        The proper motion along the stream in ``Lambda`` (including the\n",
    "        ``cos(Beta)`` factor) for this object (``pm_Beta`` must also be given).\n",
    "    pm_Beta : :class:`~astropy.units.Quantity`, optional, must be keyword\n",
    "        The proper motion in Declination for this object (``pm_ra_cosdec`` must\n",
    "        also be given).\n",
    "    radial_velocity : :class:`~astropy.units.Quantity`, optional, must be keyword\n",
    "        The radial velocity of this object.\n",
    "\n",
    "    \"\"\"\n",
    "    default_representation = coord.SphericalRepresentation\n",
    "    default_differential = coord.SphericalCosLatDifferential\n",
    "\n",
    "    frame_specific_representation_info = {\n",
    "        coord.SphericalRepresentation: [\n",
    "            coord.RepresentationMapping('lon', 'Lambda'),\n",
    "            coord.RepresentationMapping('lat', 'Beta'),\n",
    "            coord.RepresentationMapping('distance', 'distance')],\n",
    "        coord.SphericalCosLatDifferential: [\n",
    "            coord.RepresentationMapping('d_lon_coslat', 'pm_Lambda_cosBeta'),\n",
    "            coord.RepresentationMapping('d_lat', 'pm_Beta'),\n",
    "            coord.RepresentationMapping('d_distance', 'radial_velocity')],\n",
    "        coord.SphericalDifferential: [\n",
    "            coord.RepresentationMapping('d_lon', 'pm_Lambda'),\n",
    "            coord.RepresentationMapping('d_lat', 'pm_Beta'),\n",
    "            coord.RepresentationMapping('d_distance', 'radial_velocity')]\n",
    "    }\n",
    "\n",
    "    frame_specific_representation_info[coord.UnitSphericalRepresentation] = \\\n",
    "        frame_specific_representation_info[coord.SphericalRepresentation]\n",
    "    frame_specific_representation_info[coord.UnitSphericalCosLatDifferential] = \\\n",
    "        frame_specific_representation_info[coord.SphericalCosLatDifferential]\n",
    "    frame_specific_representation_info[coord.UnitSphericalDifferential] = \\\n",
    "        frame_specific_representation_info[coord.SphericalDifferential]\n",
    "\n",
    "SGR_MATRIX = Rotmat\n",
    "\n",
    "@frame_transform_graph.transform(coord.StaticMatrixTransform, coord.ICRS, StreamFrame)\n",
    "def galactic_to_sgr():\n",
    "    \"\"\" Compute the transformation matrix from Galactic spherical to\n",
    "        heliocentric Stream coordinates.\n",
    "    \"\"\"\n",
    "    return SGR_MATRIX\n",
    "\n",
    "\n",
    "@frame_transform_graph.transform(coord.StaticMatrixTransform, StreamFrame, coord.ICRS)\n",
    "def sgr_to_galactic():\n",
    "    \"\"\" Compute the transformation matrix from heliocentric Stream coordinates to\n",
    "        spherical Galactic.\n",
    "    \"\"\"\n",
    "    return matrix_transpose(SGR_MATRIX)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 102,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<Galactic Coordinate: (l, b) in deg\n",
       "    [(319.65555678, -54.869743  ), (311.96572016, -57.04543154)]\n",
       " (pm_l_cosb, pm_b) in mas / yr\n",
       "    [(-0.35399007, -1.36919357), (-1.37411612, -2.47220648)]>"
      ]
     },
     "execution_count": 102,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "coo = coord.FK5(ra=alpha*u.deg, dec=delta*u.deg, pm_ra_cosdec=pmracosdec*u.mas/u.yr, pm_dec=pmdec*u.mas/u.yr)\n",
    "coo = coo.transform_to(coord.Galactic)\n",
    "coo"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 103,
   "metadata": {},
   "outputs": [],
   "source": [
    "newcoo_apy = coo.transform_to(StreamFrame)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 104,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "('phi1 (along the stream)', array([-0.25226143,  4.56739172]))\n",
      "('phi2 (perpedincular to the stream)', array([-0.00095034, -0.01003025]))\n",
      "('pmphi1cosphi2 (along the stream)', array([0.99243716, 2.25180317]))\n",
      "('pmphi2 (perpedincular to the stream)', array([1.00750607, 1.71154388]))\n"
     ]
    }
   ],
   "source": [
    "print(\"phi1 (along the stream)\", newcoo[0])\n",
    "print(\"phi2 (perpedincular to the stream)\", newcoo[1])\n",
    "print(\"pmphi1cosphi2 (along the stream)\", newPM[0])\n",
    "print(\"pmphi2 (perpedincular to the stream)\", newPM[1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 105,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<StreamFrame Coordinate: (Lambda, Beta) in deg\n",
       "    [(359.74773989, -0.00095345), (  4.56739304, -0.01003266)]\n",
       " (pm_Lambda_cosBeta, pm_Beta) in mas / yr\n",
       "    [(0.99243702, 1.00750621), (2.25180293, 1.71154421)]>"
      ]
     },
     "execution_count": 105,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "newcoo_apy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python [conda env:anaconda]",
   "language": "python",
   "name": "conda-env-anaconda-py"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.15"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}