{"cells":[{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true,"collapsed":true},"cell_type":"code","source":"# This Python 3 environment comes with many helpful analytics libraries installed\n# It is defined by the kaggle/python docker image: https://github.com/kaggle/docker-python\n# For example, here's several helpful packages to load in \n\nimport numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\nimport matplotlib.pyplot as plt\n\n# Input data files are available in the \"../input/\" directory.\n# For example, running this (by clicking run or pressing Shift+Enter) will list the files in the input directory\n\nimport os\nprint(os.listdir(\"../input\"))\n\n# Any results you write to the current directory are saved as output.","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"collapsed":true,"_uuid":"a99a4ec78a408ed07d7f56b38437911ba5cba1a8"},"cell_type":"code","source":"#v2: moved most inline math expressions into equation environments\n#v3: added discussion of cylindrical and spherical coordinates, corrected formulas involving charge sign\n#v4: typos corrected","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"3e50f025bbfdd2394792825c6f67ee4cdaacc575"},"cell_type":"markdown","source":"I collect some details concerning the geometry of the helix for charged particles in a constant magnetic field and discuss why the coordinate transformations used in the DBSCAN Benchmark kernel https://www.kaggle.com/mikhailhushchyn/dbscan-benchmark are useful."},{"metadata":{"_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","collapsed":true,"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","trusted":false},"cell_type":"markdown","source":"**Cartesian coordinates**\n\nWe fix the following notation:\n1. Particle rest mass $m_0$ and electric charge $q$.\n2. $m=\\gamma m_0$ with relativistic factor $\\gamma$ (considered constant).\n3. Coordinate system $(x,y,z)$.\n4. Magnetic field $\\vec{B}=Be_z$.\n5. Initial coordinates $x_0=y_0=z_0=0$ (for simplicity, see comment below).\n6. Initial velocities $v_{x0}, v_{y0}, v_{z0}$ and initial momenta $p_{x0}, p_{y0}, p_{z0}$ with $p_{i0}=m v_{i0}$ for $i=x,y,z$.\n7. Longitudinal momentum $p_\\parallel = p_{z0}$, transversal momentum \n\\begin{equation*}\np_\\perp = \\sqrt{p_{x0}^2+p_{y0}^2},\n\\end{equation*}\nand total momentum \n\\begin{equation*}\np=\\sqrt{p_\\parallel^2+p_\\perp^2}.\n\\end{equation*}\n\nWe define:\n\\begin{align*}\n\\omega &= \\frac{qB}{m}\\\\\nR&= \\frac{p_\\perp}{qB}\n\\end{align*}\nand an angle $\\phi_0$ by\n\\begin{align*}\n\\cos\\phi_0&=-\\frac{p_{y0}}{p_\\perp}\\\\\n\\sin\\phi_0&=\\frac{p_{x0}}{p_\\perp}.\n\\end{align*}\nThe system of differential equations\n\\begin{align*}\nm\\ddot{x} &= qB\\dot{y}\\\\\nm\\ddot{y} &= -qB\\dot{x}\\\\\nm\\ddot{z} &= 0\n\\end{align*}\nhas the helix solution\n\\begin{align*}\nx(t) &= R(\\cos(\\phi_0-\\omega t)-\\cos(\\phi_0))\\\\\ny(t) &= R(\\sin(\\phi_0-\\omega t)-\\sin(\\phi_0))\\\\\nz(t) &= R\\frac{p_\\parallel}{p_\\perp}\\omega t.\n\\end{align*}\nIf the initial coordinates are non-zero, we just add $x_0$, $y_0$ and $z_0$ to these expressions.\n\nProjected onto the $x$-$y$-plane, the helix forms a circle of radius $|R|$ with center at the point \n\\begin{equation*}\n-R(\\cos\\phi_0,\\sin\\phi_0).\n\\end{equation*}\nThe motion in $z$-direction is of constant velocity."},{"metadata":{"_uuid":"b66f4aa9bee7dab0a28e1dcf0b5d5162e55e6992"},"cell_type":"markdown","source":"Using\n\\begin{align*}\n\\cos^2(\\alpha)+\\sin^2(\\alpha)&=1\\\\\n\\cos(\\alpha)\\cos(\\beta)+\\sin(\\alpha)\\sin(\\beta)&=\\cos(\\alpha-\\beta),\n\\end{align*}\nwe have\n\\begin{align*}\nx^2+y^2 &= 2R^2(1-\\cos(\\omega t))\\\\\n&\\approx R^2\\left(\\omega^2t^2-\\frac{1}{12}\\omega^4t^4  \\right),\n\\end{align*}\nwhere the approximation in the second line follows from Taylor's formula and is valid for small $\\omega t$.\n\nWe define a new coordinate\n\\begin{equation*}\nz_2=\\frac{z}{\\sqrt{x^2+y^2}}=\\pm\\frac{p_\\parallel}{p_\\perp}\\frac{\\omega t}{\\sqrt{2(1-\\cos(\\omega t))}},\n\\end{equation*}\nwhich is independent of $R$ and $\\phi_0$ and the sign $\\pm$ is the sign of the charge $q$.\n\nFor small $\\omega t$ we get\n\\begin{align*}\nz_2&\\approx \\frac{p_\\parallel}{p_\\perp}\\frac{1}{\\sqrt{1-\\frac{1}{12}\\omega^2t^2}}\\\\\n&\\approx\\frac{p_\\parallel}{p_\\perp}\\left(1+\\frac{1}{24}\\omega^2t^2\\right).\n\\end{align*}\nHence up to first order in $\\omega t$ the coordinate $z_2$ is constant for each helix and given by\n\\begin{equation*}\nz_2\\approx \\frac{p_\\parallel}{p_\\perp}.\n\\end{equation*}"},{"metadata":{"_uuid":"4284fe60b5e7eab6a50817e380f333a1569ada82"},"cell_type":"markdown","source":"We set\n\\begin{align*}\nx_2 &= \\frac{x}{\\sqrt{x^2+y^2+z^2}}\\\\\ny_2 &= \\frac{y}{\\sqrt{x^2+y^2+z^2}}.\n\\end{align*}\nThen the expression\n\\begin{align*}\nx_2^2+y_2^2 &= \\frac{x^2+y^2}{x^2+y^2+z^2}\\\\\n&=\\frac{1}{1+z_2^2}\n\\end{align*}\nis independent of $R$ and $\\phi_0$.\n\nFor small $\\omega t$, where $z_2\\approx \\frac{p_\\parallel}{p_\\perp}$, the point $(x_2,y_2)$ lies approximately on a circle of radius\n\\begin{equation*}\n\\frac{1}{\\sqrt{1+\\frac{p_\\parallel^2}{p_\\perp^2}}} = \\frac{p_\\perp}{\\sqrt{p_\\parallel^2+p_\\perp^2}} = \\frac{p_\\perp}{p},\n\\end{equation*}\nwhere $p$ is the total momentum.\n\nUsing Taylor expansion for small $\\omega t$ we can calculate\n\\begin{align*}\nx(t)&\\approx R\\sin(\\phi_0)\\omega t\\\\\ny(t)&\\approx -R\\cos(\\phi_0)\\omega t.\n\\end{align*}\nThis implies\n\\begin{align*}\n\\lim_{t\\rightarrow 0}x_2(t)&=\\sin(\\phi_0)\\frac{p_\\perp}{p}\\\\\n\\lim_{t\\rightarrow 0}y_2(t)&=-\\cos(\\phi_0)\\frac{p_\\perp}{p}.\n\\end{align*}\nWe see that the (idealized) initial point $(x_2,y_2)$ on the circle of radius $\\frac{p_\\perp}{p}$ depends on the angle $\\phi_0$, even though $x_0=y_0=0$. Hence in the coordinates $x_2,y_2$ the different helices get separated, depending on the angle $\\phi_0$."},{"metadata":{"_uuid":"6f4c9629f8910234ecb3a8454f22cfba1f0f89c6"},"cell_type":"markdown","source":"**Cylindrical coordinates**\n\nCylindrical coordinates $(r,\\phi,z)$, with $r\\geq 0, \\phi\\in[0,2\\pi]$, are given by\n\\begin{align*}\nx &= r\\cos\\phi\\\\\ny &= r\\sin\\phi\\\\\nz &= z.\n\\end{align*}\nWe have  already calculated\n\\begin{align*}\nr&=\\sqrt{x^2+y^2}\\\\\n&=\\sqrt{2R^2(1-\\cos(\\omega t))}\\\\\n&=\\left|2R\\sin\\left(\\frac{\\omega t}{2}\\right)\\right|,\n\\end{align*}\nwhere we used in the final step the identity \n\\begin{equation*}\n\\cos(\\omega t) = \\cos^2\\left(\\frac{\\omega t}{2}\\right)-\\sin^2\\left(\\frac{\\omega t}{2}\\right).\n\\end{equation*}\n\nUsing the identities\n\\begin{align*}\n\\cos \\alpha-\\cos \\beta &=-2\\sin {\\frac  {\\alpha+\\beta}{2}}\\sin {\\frac  {\\alpha-\\beta}{2}}\\\\\n\\sin \\alpha-\\sin \\beta &=2\\cos \\frac{\\alpha+\\beta}{2}\\sin \\frac{\\alpha-\\beta}{2}\n\\end{align*}\nwe can write\n\\begin{align*}\nx(t) &= 2R\\sin\\left(\\phi_0-\\frac{1}{2}\\omega t\\right)\\sin(\\omega t)\\\\\ny(t) &= -2R\\cos\\left(\\phi_0-\\frac{1}{2}\\omega t\\right)\\sin(\\omega t).\n\\end{align*}\nIt follows that\n\\begin{align*}\n\\tan\\phi &= \\frac{y}{x}\\\\\n&=\\frac{-\\cos\\left(\\phi_0-\\frac{1}{2}\\omega t\\right)}{\\sin\\left(\\phi_0-\\frac{1}{2}\\omega t\\right)}\\\\ \n&=\\frac{\\sin\\left(\\phi_0-\\frac{1}{2}\\omega t-\\frac{\\pi}{2}\\right)}{\\cos\\left(\\phi_0-\\frac{1}{2}\\omega t-\\frac{\\pi}{2}\\right)}\\\\\n&=\\tan\\left(\\phi_0-\\frac{1}{2}\\omega t-\\frac{\\pi}{2}\\right).\n\\end{align*}\nHence the helix is given in cylindrical coordinates by\n\\begin{align*}\nr(t)&=\\left|2R\\sin\\left(\\frac{\\omega t}{2}\\right)\\right|\\\\\n\\phi(t)&=\\phi_0-\\frac{1}{2}\\omega t-\\frac{\\pi}{2}\\\\\nz(t)&=R\\frac{p_\\parallel}{p_\\perp}\\omega t.\n\\end{align*}\nNotice that strictly speaking the angle $\\phi$ is only defined if $r\\neq 0$. However, the expression for $\\phi(t)$ above can still be used.\n\nUp to terms of second order in $\\omega t$ we get\n\\begin{align*}\nr(t)&\\approx \\left|R\\omega t\\right|\\\\\n\\phi(t)&=\\phi_0-\\frac{1}{2}\\omega t-\\frac{\\pi}{2}\\\\\nz(t)&=R\\frac{p_\\parallel}{p_\\perp}\\omega t.\n\\end{align*}"},{"metadata":{"_uuid":"45eeeac94bf79e40e245dc76f3703a7ac716f51c"},"cell_type":"markdown","source":"**Spherical coordinates**\n\nSpherical coordinates $(r,\\phi,\\theta)$, with $r\\geq 0, \\phi\\in[0,2\\pi], \\theta\\in[0,\\pi]$, are given by\n\\begin{align*}\nx &= r\\sin\\theta\\cos\\phi\\\\\ny &= r\\sin\\theta\\sin\\phi\\\\\nz &= r\\cos\\theta.\n\\end{align*}\nUsing our calculations above we can describe the helix by\n\\begin{align*}\nr(t)&=|R|\\sqrt{\\left(\\frac{p_\\parallel}{p_\\perp}\\right)^2\\omega^2t^2+4\\sin^2\\left(\\frac{\\omega t}{2}\\right)}\\\\\n\\phi(t)&=\\phi_0-\\frac{1}{2}\\omega t-\\frac{\\pi}{2}\\\\\n\\theta(t)&=\\arctan\\left(\\frac{p_\\perp}{p_\\parallel}\\frac{2\\sin\\left(\\frac{\\omega t}{2}\\right)}{\\omega t}\\right).\n\\end{align*}\nUp to terms of second order in $\\omega t$ we get\n\\begin{align*}\nr(t)&\\approx \\frac{p}{p_\\perp}|R\\omega t|\\\\\n\\phi(t)&=\\phi_0-\\frac{1}{2}\\omega t-\\frac{\\pi}{2}\\\\\n\\theta(t)&\\approx\\arctan\\left(\\frac{p_\\perp}{p_\\parallel}\\right).\n\\end{align*}\nIn particular, the angle $\\theta$ is approximately constant.\n"},{"metadata":{"_uuid":"116f83502b3c19957b94d57a34e3bc1be9402114"},"cell_type":"markdown","source":"**Visualization**\n\nLet's visualize the formulas for $x_2, y_2, z_2$."},{"metadata":{"trusted":true,"collapsed":true,"_uuid":"3f81c22446715f3a4e5d18adaaba619ba11d7f29"},"cell_type":"code","source":"# Use time parameter s = omega*t\n\ndef x(s,R,phi_0):\n    return R*(np.cos(phi_0-s)-np.cos(phi_0))\n\ndef y(s,R,phi_0):\n    return R*(np.sin(phi_0-s)-np.sin(phi_0))\n\n# p_L = p_longitudinal, p_T = p_transversal\ndef z(s,R,p_T,p_L):\n    return R*(p_L/p_T)*s","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"collapsed":true,"_uuid":"883113e1d8fb14b198523803177ba1130c0db4de"},"cell_type":"code","source":"def r1(s,R,phi_0,p_T,p_L):\n    return np.sqrt(x(s,R,phi_0)**2+y(s,R,phi_0)**2+z(s,R,p_T,p_L)**2)\n\ndef r2(s,R,phi_0):\n    return np.sqrt(x(s,R,phi_0)**2+y(s,R,phi_0)**2)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"collapsed":true,"_uuid":"6b6d5d26f4789c9ab583ad6203799333accd9f18"},"cell_type":"code","source":"def x2(s,R,phi_0,p_T,p_L):\n    return x(s,R,phi_0)/r1(s,R,phi_0,p_T,p_L)\n\ndef y2(s,R,phi_0,p_T,p_L):\n    return y(s,R,phi_0)/r1(s,R,phi_0,p_T,p_L)\n\ndef z2(s,R,phi_0,p_T,p_L):\n    return z(s,R,p_T,p_L)/r2(s,R,phi_0)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"collapsed":true,"_uuid":"0fe7cf55d44f26a3e677da75e01934b142fed86b"},"cell_type":"code","source":"# Set some values for radius R and momenta p_L, p_T, used throughout the examples\nR = 1\np_L = 100\np_T = 10","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"a657978e15162e6520f13c356e52532c405a798c","collapsed":true},"cell_type":"code","source":"# Start with a short time interval [0.01, 0.5]\n# We do not start in s=0 to avoid dividing by 0 when calculating x2, y2, z2\nS = np.linspace(0.01, 0.5, 200)\n\n# Plot (x,y) for 10 different values for phi_0, corresponding to different initial velocity vectors\nfor phi_0 in np.linspace(0, np.pi/2, 10):\n      \n    X = x(S,R,phi_0)\n    Y = y(S,R,phi_0)\n    plt.axis(\"equal\")\n    plt.xlabel(\"x\")\n    plt.ylabel(\"y\")\n    \n    plt.plot(X,Y)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"5a0d75894ac21d374d921cb94296559cc1e59cca","collapsed":true},"cell_type":"code","source":"# Plot transformed coordinates (x2,y2) for the same angles\n# The transformed curves align on a circle segment\nfor phi_0 in np.linspace(0, np.pi/2, 10):\n    \n    X2 = x2(S,R,phi_0,p_T,p_L)\n    Y2 = y2(S,R,phi_0,p_T,p_L)\n    \n    plt.axis(\"equal\")\n    plt.xlabel(\"x2\")\n    plt.ylabel(\"y2\")\n    \n    plt.plot(X2,Y2)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"cdf37c702f8d7e2c6617189703d05588f38d8151","collapsed":true},"cell_type":"code","source":"# If we make the time interval very short, the different starting points, depending on phi_0, become obvious\nS = np.linspace(0.01, 0.03, 200)\n\nfor phi_0 in np.linspace(0, np.pi/2, 10):\n    \n    X2 = x2(S,R,phi_0,p_T,p_L)\n    Y2 = y2(S,R,phi_0,p_T,p_L)\n    \n    plt.axis(\"equal\")\n    plt.xlabel(\"x2\")\n    plt.ylabel(\"y2\")\n    \n    plt.plot(X2,Y2)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"640bf6429022771539cd18808cb71ccc74355a15","collapsed":true},"cell_type":"code","source":"S = np.linspace(0.01, 0.5, 200)\n\n# The circle becomes clearer if we let phi_0 run from 0 to 2*pi\nfor phi_0 in np.linspace(0, 2*np.pi, 10):\n    \n    X2 = x2(S,R,phi_0,p_T,p_L)\n    Y2 = y2(S,R,phi_0,p_T,p_L)\n    \n    plt.axis(\"equal\")\n    plt.xlabel(\"x2\")\n    plt.ylabel(\"y2\")\n    \n    plt.plot(X2,Y2)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"761a673387104333785667eb3bd5d31eb391aa1b","collapsed":true},"cell_type":"code","source":"# If the time interval becomes large, the segments do not align as well on the circle as before\n\nS = np.linspace(0.01, 2, 200)\n\nfor phi_0 in np.linspace(0, np.pi/2, 10):\n    \n    X2 = x2(S,R,phi_0,p_T,p_L)\n    Y2 = y2(S,R,phi_0,p_T,p_L)\n    \n    plt.axis(\"equal\")\n    plt.xlabel(\"x2\")\n    plt.ylabel(\"y2\")\n    \n    plt.plot(X2,Y2)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"ab2640e19a2e247cc35ad8db294f56107ef69135","collapsed":true},"cell_type":"code","source":"for phi_0 in np.linspace(0,2*np.pi,10):\n    \n    X2 = x2(S,R,phi_0,p_T,p_L)\n    Y2 = y2(S,R,phi_0,p_T,p_L)\n    \n    plt.axis(\"equal\")\n    plt.xlabel(\"x2\")\n    plt.ylabel(\"y2\")\n    \n    plt.plot(X2,Y2)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"57bb763d357bd51a68fb92d535f34cc481dc7d63","collapsed":true},"cell_type":"code","source":"# We now plot z2\n# Back to the short time interval\nS = np.linspace(0.01, 0.5, 200)\n\n# z2 is independent of phi_0, set phi_0 = 0\nphi_0 = 0\n\n# Plot the time dependency of z2\n# z2 is approximately constant and the quadratic time dependency is apparent\nZ2 = z2(S,R,phi_0,p_T,p_L)\nplt.xlabel(\"s\")\nplt.ylabel(\"z2\")\nplt.plot(S,Z2)\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"c4929ff2557b4976464b07d7958a8b6f1f5a0b42","collapsed":true},"cell_type":"code","source":"# With a larger time interval, higher order terms kick in\nS = np.linspace(0.01, 5, 200)\n\nphi_0 = 0\n\nZ2 = z2(S,R,phi_0,p_T,p_L)\nplt.xlabel(\"s\")\nplt.ylabel(\"z2\")\nplt.plot(S,Z2)\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"47529da09822a5fc2e0920ad278a77baeeb3778e","collapsed":true},"cell_type":"code","source":"# The denominator of z2 becomes singular at integer multiples of 2*pi\nS = np.linspace(0.01, 7, 200)\n\nphi_0 = 0\n\nZ2 = z2(S,R,phi_0,p_T,p_L)\nplt.xlabel(\"s\")\nplt.ylabel(\"z2\")\nplt.plot(S,Z2)\nplt.show()","execution_count":null,"outputs":[]}],"metadata":{"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"name":"python","version":"3.6.5","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"}},"nbformat":4,"nbformat_minor":1}