{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Muistiinpanoja erityisestä suhteellisuusteoriasta"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import sympy as sp\n",
    "from numpy.linalg import eig\n",
    "from sympy import *"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "phi=symbols('phi')\n",
    "rotation=Matrix([[cos(phi), -sin(phi)],\n",
    "                 [sin(phi), cos(phi)]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x = sp.symbols(\"x\")\n",
    "sp.simplify(cos(x)**2 + sin(x)**2)  # returns 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "x"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "cos(acos(x))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Minkowskin metriikan matriisi\n",
    "eta_matrix=Matrix([[1,0,0,0],[0,-1,0,0],[0,0,-1,0],[0,0,0,-1]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Matrix([\n",
       "[1, 0, 0, 0],\n",
       "[0, 1, 0, 0],\n",
       "[0, 0, 1, 0],\n",
       "[0, 0, 0, 1]])"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Tarkistetaan matriisitulo...pitäisi saada yksikkömatriisi\n",
    "eta_matrix*eta_matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Nelipaikka\n",
    "x,y,z,ct = sp.symbols(\"x y z ct\")  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "Xv=Matrix([ct,x,y,z])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Matrix([[ct**2 - x**2 - y**2 - z**2]])"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#Pistetulo Minkowskin avaruudessa...pitäisi olla sama kaikissa inertiaalikoordinaatistoissa\n",
    "Transpose(Xv)*eta_matrix*Xv"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "gamma,beta,v,c = sp.symbols(\"gamma beta v c\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "beta=v/c\n",
    "gamma=1/sqrt(1-Pow(beta,2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Lorentzin muunnos, kun nopeus x-suuntaan\n",
    "BoostX=Matrix([[gamma,-beta*gamma,0,0],[-beta*gamma,gamma,0,0],[0,0,1,0],[0,0,0,1]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Matrix([[1/sqrt(1 - v**2/c**2), -v/(c*sqrt(1 - v**2/c**2)), 0, 0], [-v/(c*sqrt(1 - v**2/c**2)), 1/sqrt(1 - v**2/c**2), 0, 0], [0, 0, 1, 0], [0, 0, 0, 1]])\n"
     ]
    }
   ],
   "source": [
    "print(BoostX)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Muunnettu nelipaikka\n",
    "Xvp=BoostX*Xv"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Matrix([[ct**2 - x**2 - y**2 - z**2]])"
      ]
     },
     "execution_count": 46,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "simplify(Transpose(Xvp)*eta_matrix*Xvp)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "vx,vy,vz,c=sp.symbols(\"vx vy vz c\")\n",
    "betav=Matrix([vx,vy,vz])/c"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Matrix([[vx/c], [vy/c], [vz/c]])\n"
     ]
    }
   ],
   "source": [
    "print(betav)\n",
    "betasqr=betav.dot(betav)\n",
    "gamma=1/sqrt(1-betav.dot(betav))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "vx**2/c**2 + vy**2/c**2 + vz**2/c**2\n"
     ]
    }
   ],
   "source": [
    "print(betasqr)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "vx**2/c**2\n"
     ]
    }
   ],
   "source": [
    "print(pow(betav[0],2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Kaikenlaisia Lorentzin muunnoksia eri suuntiin\n",
    "BoostX=Matrix([[gamma,-betav[0]*gamma,0,0],[-betav[0]*gamma,gamma,0,0],[0,0,1,0],[0,0,0,1]])\n",
    "BoostY=Matrix([[gamma,0,-betav[1]*gamma,0],[0,1,0,0],[-betav[1]*gamma,0,gamma,0],[0,0,0,1]])\n",
    "BoostZ=Matrix([[gamma,0,0,-betav[2]*gamma],[0,1,0,0],[0,0,1,0],[-betav[2]*gamma,0,0,gamma]])\n",
    "BoostGeneral=Matrix([[gamma,-betav[0]*gamma,-betav[1]*gamma,-betav[2]*gamma],[-betav[0]*gamma,1+(gamma-1)*pow(betav[0],2)/betasqr,(gamma-1)*betav[0]*betav[1],(gamma-1)*betav[0]*betav[2]],[-betav[1]*gamma,(gamma-1)*betav[0]*betav[1],1+(gamma-1)*pow(betav[1],2)/betasqr,(gamma-1)*betav[1]*betav[2]],[-betav[2]*gamma,(gamma-1)*betav[0]*betav[2],(gamma-1)*betav[1]*betav[2],1+(gamma-1)*pow(betav[2],2)/betasqr]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 79,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Matrix([\n",
       "[      1/sqrt(1 - vx**2/c**2 - vy**2/c**2 - vz**2/c**2),                                                          -vx/(c*sqrt(1 - vx**2/c**2 - vy**2/c**2 - vz**2/c**2)),                                                          -vy/(c*sqrt(1 - vx**2/c**2 - vy**2/c**2 - vz**2/c**2)),                                                          -vz/(c*sqrt(1 - vx**2/c**2 - vy**2/c**2 - vz**2/c**2))],\n",
       "[-vx/(c*sqrt(1 - vx**2/c**2 - vy**2/c**2 - vz**2/c**2)), 1 + vx**2*(-1 + 1/sqrt(1 - vx**2/c**2 - vy**2/c**2 - vz**2/c**2))/(c**2*(vx**2/c**2 + vy**2/c**2 + vz**2/c**2)),                                              vx*vy*(-1 + 1/sqrt(1 - vx**2/c**2 - vy**2/c**2 - vz**2/c**2))/c**2,                                              vx*vz*(-1 + 1/sqrt(1 - vx**2/c**2 - vy**2/c**2 - vz**2/c**2))/c**2],\n",
       "[-vy/(c*sqrt(1 - vx**2/c**2 - vy**2/c**2 - vz**2/c**2)),                                              vx*vy*(-1 + 1/sqrt(1 - vx**2/c**2 - vy**2/c**2 - vz**2/c**2))/c**2, 1 + vy**2*(-1 + 1/sqrt(1 - vx**2/c**2 - vy**2/c**2 - vz**2/c**2))/(c**2*(vx**2/c**2 + vy**2/c**2 + vz**2/c**2)),                                              vy*vz*(-1 + 1/sqrt(1 - vx**2/c**2 - vy**2/c**2 - vz**2/c**2))/c**2],\n",
       "[-vz/(c*sqrt(1 - vx**2/c**2 - vy**2/c**2 - vz**2/c**2)),                                              vx*vz*(-1 + 1/sqrt(1 - vx**2/c**2 - vy**2/c**2 - vz**2/c**2))/c**2,                                              vy*vz*(-1 + 1/sqrt(1 - vx**2/c**2 - vy**2/c**2 - vz**2/c**2))/c**2, 1 + vz**2*(-1 + 1/sqrt(1 - vx**2/c**2 - vy**2/c**2 - vz**2/c**2))/(c**2*(vx**2/c**2 + vy**2/c**2 + vz**2/c**2))]])"
      ]
     },
     "execution_count": 79,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "BoostGeneral.subs(vx,0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Matrix([[1/sqrt(1 - vx**2/c**2 - vy**2/c**2 - vz**2/c**2), -vx/(c*sqrt(1 - vx**2/c**2 - vy**2/c**2 - vz**2/c**2)), -vy/(c*sqrt(1 - vx**2/c**2 - vy**2/c**2 - vz**2/c**2)), -vz/(c*sqrt(1 - vx**2/c**2 - vy**2/c**2 - vz**2/c**2))], [-vx/(c*sqrt(1 - vx**2/c**2 - vy**2/c**2 - vz**2/c**2)), 1 + vx**2*(-1 + 1/sqrt(1 - vx**2/c**2 - vy**2/c**2 - vz**2/c**2))/(c**2*(vx**2/c**2 + vy**2/c**2 + vz**2/c**2)), vx*vy*(-1 + 1/sqrt(1 - vx**2/c**2 - vy**2/c**2 - vz**2/c**2))/c**2, vx*vz*(-1 + 1/sqrt(1 - vx**2/c**2 - vy**2/c**2 - vz**2/c**2))/c**2], [-vy/(c*sqrt(1 - vx**2/c**2 - vy**2/c**2 - vz**2/c**2)), vx*vy*(-1 + 1/sqrt(1 - vx**2/c**2 - vy**2/c**2 - vz**2/c**2))/c**2, 1 + vy**2*(-1 + 1/sqrt(1 - vx**2/c**2 - vy**2/c**2 - vz**2/c**2))/(c**2*(vx**2/c**2 + vy**2/c**2 + vz**2/c**2)), vy*vz*(-1 + 1/sqrt(1 - vx**2/c**2 - vy**2/c**2 - vz**2/c**2))/c**2], [-vz/(c*sqrt(1 - vx**2/c**2 - vy**2/c**2 - vz**2/c**2)), vx*vz*(-1 + 1/sqrt(1 - vx**2/c**2 - vy**2/c**2 - vz**2/c**2))/c**2, vy*vz*(-1 + 1/sqrt(1 - vx**2/c**2 - vy**2/c**2 - vz**2/c**2))/c**2, 1 + vz**2*(-1 + 1/sqrt(1 - vx**2/c**2 - vy**2/c**2 - vz**2/c**2))/(c**2*(vx**2/c**2 + vy**2/c**2 + vz**2/c**2))]])\n"
     ]
    }
   ],
   "source": [
    "print(BoostGeneral)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 68,
   "metadata": {},
   "outputs": [],
   "source": [
    "Xvp=BoostGeneral*Xv"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 72,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Matrix([[(c**2*(vx**2 + vy**2 + vz**2)**2*(-c*ct + vx*x + vy*y + vz*z)**2 - (-c**2*x*(-vx**2*(sqrt((c**2 - vx**2 - vy**2 - vz**2)/c**2) - 1) + sqrt((c**2 - vx**2 - vy**2 - vz**2)/c**2)*(vx**2 + vy**2 + vz**2)) + c*ct*vx*(vx**2 + vy**2 + vz**2) + vx*(sqrt((c**2 - vx**2 - vy**2 - vz**2)/c**2) - 1)*(vy*y + vz*z)*(vx**2 + vy**2 + vz**2))**2 - (-c**2*y*(-vy**2*(sqrt((c**2 - vx**2 - vy**2 - vz**2)/c**2) - 1) + sqrt((c**2 - vx**2 - vy**2 - vz**2)/c**2)*(vx**2 + vy**2 + vz**2)) + c*ct*vy*(vx**2 + vy**2 + vz**2) + vy*(sqrt((c**2 - vx**2 - vy**2 - vz**2)/c**2) - 1)*(vx*x + vz*z)*(vx**2 + vy**2 + vz**2))**2 - (-c**2*z*(-vz**2*(sqrt((c**2 - vx**2 - vy**2 - vz**2)/c**2) - 1) + sqrt((c**2 - vx**2 - vy**2 - vz**2)/c**2)*(vx**2 + vy**2 + vz**2)) + c*ct*vz*(vx**2 + vy**2 + vz**2) + vz*(sqrt((c**2 - vx**2 - vy**2 - vz**2)/c**2) - 1)*(vx*x + vy*y)*(vx**2 + vy**2 + vz**2))**2)/(c**2*(vx**2 + vy**2 + vz**2)**2*(c**2 - vx**2 - vy**2 - vz**2))]])"
      ]
     },
     "execution_count": 72,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "vx=0\n",
    "vz=0\n",
    "simplify(Transpose(Xvp)*eta_matrix*Xvp)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Field transformation under Lorentz transformation\n",
    "* Lorentz force under these transformations\n",
    "* Maxwell's equations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "eta_matrix=Matrix([[1,0,0,0],[0,-1,0,0],[0,0,-1,0],[0,0,0,-1]])\n",
    "x,y,z,ct = sp.symbols(\"x y z ct\")\n",
    "Xv=Matrix([ct,x,y,z])\n",
    "Transpose(Xv)*eta_matrix*Xv\n",
    "gamma,beta,v,c = sp.symbols(\"gamma beta v c\")\n",
    "beta=v/c\n",
    "gamma=1/sqrt(1-Pow(beta,2))\n",
    "gamma=1/sqrt(1-Pow(beta,2))\n",
    "BoostX=Matrix([[gamma,-beta*gamma,0,0],[-beta*gamma,gamma,0,0],[0,0,1,0],[0,0,0,1]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "Xvp=BoostX*Xv"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Matrix([\n",
       "[ct/sqrt(1 - v**2/c**2) - v*x/(c*sqrt(1 - v**2/c**2))],\n",
       "[x/sqrt(1 - v**2/c**2) - ct*v/(c*sqrt(1 - v**2/c**2))],\n",
       "[                                                   y],\n",
       "[                                                   z]])"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Xvp"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Field tensor"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "Bx,By,Bz,Ex,Ey,Ez = sp.symbols(\"Bx By Bz Ex Ey Ez\")\n",
    "Fmat=Matrix([[0,-Ex/c,-Ey/c,-Ez/c],[Ex/c,0,-Bz,By],[Ey/c,Bz,0,-Bx],[Ez/c,-By,Bx,0]])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Transform the field tensor into new frame via boost along x."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "Fmatp=simplify(BoostX*Fmat*BoostX)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Matrix([\n",
       "[                                        0,                                                    -Ex/c,                      (Bz*v - Ey)/(c*sqrt(1 - v**2/c**2)),                    -(By*v + Ez)/(c*sqrt(1 - v**2/c**2))],\n",
       "[                                     Ex/c,                                                        0, sqrt((c**2 - v**2)/c**2)*(-Bz*c**2 + Ey*v)/(c**2 - v**2), sqrt((c**2 - v**2)/c**2)*(By*c**2 + Ez*v)/(c**2 - v**2)],\n",
       "[(-Bz*v + Ey)/(c*sqrt((c**2 - v**2)/c**2)),  sqrt((c**2 - v**2)/c**2)*(Bz*c**2 - Ey*v)/(c**2 - v**2),                                                        0,                                                     -Bx],\n",
       "[      (By*v + Ez)/(c*sqrt(1 - v**2/c**2)), sqrt((c**2 - v**2)/c**2)*(By*c**2 + Ez*v)/(-c**2 + v**2),                                                       Bx,                                                       0]])"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Fmatp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Ex/c and Ex/c\n"
     ]
    }
   ],
   "source": [
    "i=1\n",
    "j=0\n",
    "print(Fmat[i,j],\"and\",Fmatp[i,j])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "Bxp,Byp,Bzp,Exp,Eyp,Ezp = sp.symbols(\"Bxp Byp Bzp Exp Eyp Ezp\")\n",
    "Byp=Fmatp[1,3];Bzp=Fmatp[2,1];Bxp=Fmatp[3,2];Exp=c*Fmatp[1,0];Eyp=c*Fmatp[2,0];Ezp=c*Fmatp[3,0];"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Ex"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Exp"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Checking that tensor is anti-symmetric...in case of typos"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "i=3\n",
    "j=1\n",
    "simplify(Fmatp[i,j]+Fmatp[j,i])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Checking that tensor is anti-symmetric...in case of typos"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "simplify(Fmat[i,j]+Fmat[j,i])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It seems that these results are consistent with the results given in Wikipedia https://en.wikipedia.org/wiki/Classical_electromagnetism_and_special_relativity\n",
    "Some sources seem to have differences in signs and sometimes factors of c^2 which might be due to c=1 assumptions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [],
   "source": [
    "vx,vy,vz = sp.symbols(\"vx vy vz\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [],
   "source": [
    "Bpv=Matrix([Bxp,Byp,Bzp]);\n",
    "Epv=Matrix([Exp,Eyp,Ezp]);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1/sqrt(1 - v**2/c**2)"
      ]
     },
     "execution_count": 35,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "gamma"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Relativistinen nopeuksien summaus\n",
    "vp=Matrix([(vx+v)/(1+v*vx/Pow(c,2)),vy/gamma/(1+v*vx/Pow(c,2)),vz/gamma/(1+v*vx/Pow(c,2))])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 80,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Lorentz force/q in the new frame"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [],
   "source": [
    "Fp=Epv+vp.cross(Bpv)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Ex + vy*sqrt((c**2 - v**2)/c**2)*sqrt(1 - v**2/c**2)*(Bz*c**2 - Ey*v)/((1 + v*vx/c**2)*(c**2 - v**2)) - vz*sqrt((c**2 - v**2)/c**2)*sqrt(1 - v**2/c**2)*(By*c**2 + Ez*v)/((1 + v*vx/c**2)*(c**2 - v**2))"
      ]
     },
     "execution_count": 38,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Fp[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [],
   "source": [
    "#original frame\n",
    "vv=Matrix([vx,vy,vz]);Ev=Matrix([Ex,Ey,Ez]);Bv=Matrix([Bx,By,Bz]);\n",
    "Forig=Ev+vv.cross(Bv)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Matrix([\n",
       "[-By*vz + Bz*vy + Ex],\n",
       "[ Bx*vz - Bz*vx + Ey],\n",
       "[-Bx*vy + By*vx + Ez]])"
      ]
     },
     "execution_count": 40,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Forig"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Matrix([\n",
       "[vx],\n",
       "[vy],\n",
       "[vz]])"
      ]
     },
     "execution_count": 42,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "vv"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "See https://en.wikipedia.org/wiki/Lorentz_force\n",
    "http://www.sciencebits.com/node/176"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}