SymPy Coordinate systems.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Let's see what we can do ...\n",
    "import time\n",
    "import numpy as np\n",
    "from sympy import symbols, sin, cos, pi\n",
    "from sympy.diffgeom import Manifold, Patch, CoordSystem, Point\n",
    "from sympy.utilities.lambdify import lambdastr\n",
    "import sympy\n",
    "r = symbols('r', real=True, positive=True)\n",
    "theta, phi, x, y, z = symbols('theta phi x y z', real=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "m = Manifold('M', 3)\n",
    "patch = Patch('P', m)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "cart = CoordSystem('cart', patch)\n",
    "cyl_polar = CoordSystem('cyl_polar', patch)\n",
    "spherical = CoordSystem('spherical', patch)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "spherical.connect_to(cart, [r, theta, phi], [r * sin(theta) * cos(phi), r * sin(theta) * sin(phi), r * cos(theta)], inverse=False)\n",
    "cart.connect_to(spherical, [x, y, z], [sympy.sqrt(x**2+y**2+z**2), sympy.acos(z/(sympy.sqrt(x**2+y**2+z**2))), sympy.atan2(y, x)], inverse=False)\n",
    "cyl_polar.connect_to(cart, [r, theta, z], [r * cos(theta), r * sin(theta), z])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Matrix([\n",
       "[sin(theta)*cos(phi), r*cos(phi)*cos(theta), -r*sin(phi)*sin(theta)],\n",
       "[sin(phi)*sin(theta), r*sin(phi)*cos(theta),  r*sin(theta)*cos(phi)],\n",
       "[         cos(theta),         -r*sin(theta),                      0]])"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "spherical.jacobian(cart, [r, theta, phi])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Matrix([[r*sin(theta)*cos(phi)], [r*sin(phi)*sin(theta)], [r*cos(theta)]])\n"
     ]
    }
   ],
   "source": [
    "m = spherical.coord_tuple_transform_to(cart, [r, theta, phi])\n",
    "print(m)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Matrix([\n",
       "[        sqrt(x**2 + y**2 + z**2)],\n",
       "[acos(z/sqrt(x**2 + y**2 + z**2))],\n",
       "[                     atan2(y, x)]])"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "cart.coord_tuple_transform_to(spherical, [x, y, z])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Matrix([\n",
       "[(-x**2 + x*sqrt(x**2 + y**2) - y**2)/(x - sqrt(x**2 + y**2))],\n",
       "[                          -2*atan((x - sqrt(x**2 + y**2))/y)],\n",
       "[                                                           z]])"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "cart.coord_tuple_transform_to(cyl_polar, [x, y, z])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Matrix([\n",
       "[(-r**2*sin(phi)**2*sin(theta)**2 - r**2*sin(theta)**2*cos(phi)**2 + r*sqrt(r**2*sin(phi)**2*sin(theta)**2 + r**2*sin(theta)**2*cos(phi)**2)*sin(theta)*cos(phi))/(r*sin(theta)*cos(phi) - sqrt(r**2*sin(phi)**2*sin(theta)**2 + r**2*sin(theta)**2*cos(phi)**2))],\n",
       "[                                                                                                                                -2*atan((r*sin(theta)*cos(phi) - sqrt(r**2*sin(phi)**2*sin(theta)**2 + r**2*sin(theta)**2*cos(phi)**2))/(r*sin(phi)*sin(theta)))],\n",
       "[                                                                                                                                                                                                                                                    r*cos(theta)]])"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "cart.coord_tuple_transform_to(cyl_polar, spherical.coord_tuple_transform_to(cart, [r, theta, phi]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Matrix([\n",
       "[        sqrt(r**2*sin(phi)**2 + r**2*cos(phi)**2 + z**2)],\n",
       "[acos(z/sqrt(r**2*sin(phi)**2 + r**2*cos(phi)**2 + z**2))],\n",
       "[                           atan2(r*sin(phi), r*cos(phi))]])"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ct = cart.coord_tuple_transform_to(spherical, cyl_polar.coord_tuple_transform_to(cart, [r, phi, z]))\n",
    "ct"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Matrix([\n",
       "[            sqrt(r**2 + z**2)],\n",
       "[    acos(z/sqrt(r**2 + z**2))],\n",
       "[atan2(r*sin(phi), r*cos(phi))]])"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sympy.simplify(ct)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "atan2(r*sin(phi), r*cos(phi))"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sympy.atan2(r*sin(phi), r*cos(phi))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "atan2(sin(phi), cos(phi))"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sympy.refine(sympy.atan2(sin(phi), cos(phi)), sympy.Q.is_true(phi < 2*pi) & sympy.Q.is_true(phi > 0))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Matrix([\n",
       "[            sqrt(r**2 + z**2)],\n",
       "[    acos(z/sqrt(r**2 + z**2))],\n",
       "[atan2(r*sin(phi), r*cos(phi))]])"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sympy.trigsimp(ct)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sympy.diffgeom.rn as rn"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "CoordSystem(rectangular, Patch(origin, Manifold(R**3, 3)), (x, y, z))"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "rn.R3_r"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "CoordSystem(cylindrical, Patch(origin, Manifold(R**3, 3)), (rho, psi, z))"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "rn.R3_c"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "CoordSystem(cart, Patch(P, Manifold(M, 3)), (cart_0, cart_1, cart_2))"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "cart"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sympy.diffgeom.rn import R3_r, R3_c, R3_s, R3_origin, R3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "rho/r"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "R3_c.rho/R3_s.r"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "ll = []\n",
    "args = tuple(R3_r.transforms[R3_s][0])\n",
    "for i, c in enumerate(R3_r.coord_functions()):\n",
    "    ll.append((c, sympy.lambdify(args, R3_r.transforms[R3_s][1][i], modules=\"numpy\"),\n",
    "             lambdastr(c, R3_r.transforms[R3_s][1][i])))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Help on function <lambda> in module numpy:\n",
      "\n",
      "<lambda> lambda _Dummy_223, _Dummy_224, _Dummy_225\n",
      "    Created with lambdify. Signature:\n",
      "    \n",
      "    func(x, y, z)\n",
      "    \n",
      "    Expression:\n",
      "    \n",
      "    sqrt(_x**2 + _y**2 + _z**2)\n",
      "\n"
     ]
    }
   ],
   "source": [
    "help(ll[0][1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2.626e-01 2.615e-01\n"
     ]
    }
   ],
   "source": [
    "N = 256**3\n",
    "a, b, c = np.random.random(N), np.random.random(N), np.random.random(N)\n",
    "t1 = time.time()\n",
    "ll[0][1](a, b, c)\n",
    "t2 = time.time()\n",
    "np.sqrt(a**2 + b**2 + c**2)\n",
    "t3 = time.time()\n",
    "print (\"%0.3e %0.3e\" % (t2-t1, t3-t2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(x, <function <lambda> at 0x7fd67f9a51e0>, 'lambda x: (sqrt(_x**2 + _y**2 + _z**2))')\n"
     ]
    }
   ],
   "source": [
    "print (ll[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "mm = R3_r.transforms[R3_s][0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['C', 'D', 'H', 'LDLdecomposition', 'LDLsolve', 'LUdecomposition', 'LUdecompositionFF', 'LUdecomposition_Simple', 'LUsolve', 'QRdecomposition', 'QRsolve', 'T', '_LDLdecomposition', '__abs__', '__add__', '__array__', '__array_priority__', '__class__', '__delattr__', '__dict__', '__dir__', '__div__', '__doc__', '__eq__', '__format__', '__ge__', '__getattr__', '__getattribute__', '__getitem__', '__gt__', '__hash__', '__init__', '__init_subclass__', '__le__', '__len__', '__lt__', '__mathml__', '__matmul__', '__module__', '__mul__', '__ne__', '__neg__', '__new__', '__pow__', '__radd__', '__reduce__', '__reduce_ex__', '__repr__', '__rmatmul__', '__rmul__', '__rsub__', '__setattr__', '__setitem__', '__sizeof__', '__str__', '__sub__', '__subclasshook__', '__truediv__', '__weakref__', '_cache_eigenvects', '_cache_is_diagonalizable', '_cholesky', '_class_priority', '_diagonal_solve', '_diagonalize_clear_subproducts', '_eval_Abs', '_eval_add', '_eval_adjoint', '_eval_applyfunc', '_eval_as_real_imag', '_eval_atoms', '_eval_berkowitz_toeplitz_matrix', '_eval_berkowitz_vector', '_eval_col_del', '_eval_col_insert', '_eval_col_join', '_eval_col_op_add_multiple_to_other_col', '_eval_col_op_multiply_col_by_const', '_eval_col_op_swap', '_eval_conjugate', '_eval_det_bareiss', '_eval_det_berkowitz', '_eval_det_lu', '_eval_determinant', '_eval_diag', '_eval_diff', '_eval_echelon_form', '_eval_extract', '_eval_eye', '_eval_free_symbols', '_eval_get_diag_blocks', '_eval_has', '_eval_inverse', '_eval_is_Identity', '_eval_is_anti_symmetric', '_eval_is_diagonal', '_eval_is_echelon', '_eval_is_lower', '_eval_is_lower_hessenberg', '_eval_is_matrix_hermitian', '_eval_is_symbolic', '_eval_is_symmetric', '_eval_is_upper_hessenberg', '_eval_is_zero', '_eval_jordan_block', '_eval_matrix_mul', '_eval_matrix_mul_elementwise', '_eval_matrix_rmul', '_eval_ones', '_eval_permute_cols', '_eval_permute_rows', '_eval_pow_by_recursion', '_eval_row_del', '_eval_row_insert', '_eval_row_join', '_eval_row_op_add_multiple_to_other_row', '_eval_row_op_multiply_row_by_const', '_eval_row_op_swap', '_eval_rref', '_eval_scalar_mul', '_eval_scalar_rmul', '_eval_simplify', '_eval_tolist', '_eval_trace', '_eval_transpose', '_eval_values', '_eval_vec', '_eval_zeros', '_format_str', '_handle_creation_inputs', '_lower_triangular_solve', '_mat', '_matrix_pow_by_jordan_blocks', '_new', '_normalize_op_args', '_op_priority', '_permute_complexity_right', '_row_reduce', '_setitem', '_simplify', '_sympify', '_upper_triangular_solve', 'add', 'adjoint', 'adjugate', 'applyfunc', 'as_immutable', 'as_mutable', 'as_real_imag', 'atoms', 'berkowitz', 'berkowitz_charpoly', 'berkowitz_det', 'berkowitz_eigenvals', 'berkowitz_minors', 'charpoly', 'cholesky', 'cholesky_solve', 'cofactor', 'cofactorMatrix', 'cofactor_matrix', 'col', 'col_del', 'col_insert', 'col_join', 'col_op', 'col_swap', 'cols', 'columnspace', 'condition_number', 'conjugate', 'copy', 'copyin_list', 'copyin_matrix', 'cross', 'det', 'det_LU_decomposition', 'det_bareis', 'det_bareiss', 'diag', 'diagonal_solve', 'diagonalize', 'diff', 'doit', 'dot', 'dual', 'echelon_form', 'eigenvals', 'eigenvects', 'elementary_col_op', 'elementary_row_op', 'equals', 'evalf', 'exp', 'expand', 'extract', 'eye', 'fill', 'free_symbols', 'gauss_jordan_solve', 'get_diag_blocks', 'has', 'hstack', 'integrate', 'inv', 'inv_mod', 'inverse_ADJ', 'inverse_GE', 'inverse_LU', 'is_Identity', 'is_Matrix', 'is_MatrixExpr', 'is_anti_symmetric', 'is_diagonal', 'is_diagonalizable', 'is_echelon', 'is_hermitian', 'is_lower', 'is_lower_hessenberg', 'is_nilpotent', 'is_square', 'is_symbolic', 'is_symmetric', 'is_upper', 'is_upper_hessenberg', 'is_zero', 'jacobian', 'jordan_block', 'jordan_cell', 'jordan_cells', 'jordan_form', 'key2bounds', 'key2ij', 'left_eigenvects', 'limit', 'lower_triangular_solve', 'minor', 'minorEntry', 'minorMatrix', 'minor_submatrix', 'multiply', 'multiply_elementwise', 'n', 'norm', 'normalized', 'nullspace', 'ones', 'orthogonalize', 'permute', 'permuteBkwd', 'permuteFwd', 'permute_cols', 'permute_rows', 'pinv', 'pinv_solve', 'print_nonzero', 'project', 'rank', 'refine', 'replace', 'reshape', 'row', 'row_del', 'row_insert', 'row_join', 'row_op', 'row_swap', 'rows', 'rowspace', 'rref', 'shape', 'simplify', 'singular_values', 'solve', 'solve_least_squares', 'subs', 'table', 'tolist', 'trace', 'transpose', 'upper_triangular_solve', 'values', 'vec', 'vech', 'vstack', 'xreplace', 'zeros', 'zip_row_op']\n"
     ]
    }
   ],
   "source": [
    "print (dir(mm))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[_x, _y, _z]"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mm.values()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(_x, _y, _z)"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tuple(mm)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "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.6.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}