{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# Body-oriented mill\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "This model is introduced in \n\nP. Degond, A. Diez, M. Na, Bulk topological states in a new collective dynamics model,  arXiv:2101.10864, 2021 \n\n\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "First of all, some standard imports. \n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import os\nimport sys\nimport time\nimport math\nimport torch\nimport numpy as np \nfrom matplotlib import pyplot as plt\nimport sisyphe.models as models\nfrom sisyphe.display import save\n\nuse_cuda = torch.cuda.is_available()\ndtype = torch.cuda.FloatTensor if use_cuda else torch.FloatTensor"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Body-oriented particles with initial perpendicular twist\n\nThe system is composed of body-oriented particles which are initially uniformly scattered in a periodic box but their body-orientations are \"twisted\". The body orientation of a particle at position $(x,y,z)$ is initially: \n\n\\begin{align}\\left(\\begin{array}{ccc} 1 & 0 & 0 \\\\ 0 & \\cos(2\\pi z) & -\\sin(2\\pi z) \\\\ 0 & \\sin(2\\pi z) & \\cos(2\\pi z)\\end{array}\\right).\\end{align}\n\n\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "from sisyphe.initial import cyclotron_twist_z\n\nN = 1500000\nL = 1\nR = .025\nnu = 40\nc = 1\nkappa = 10\n\npos, bo = cyclotron_twist_z(N,L,1,kappa,dtype)\n\nsimu = models.BOAsynchronousVicsek(pos=pos,bo=bo,\n                 v=c,\n                 jump_rate=nu,kappa=kappa,\n                 interaction_radius=R,\n                 box_size=L,\n                 boundary_conditions='periodic',\n                 variant = {\"name\" : \"normalised\", \"parameters\" : {}},\n                 options = {},\n                 sampling_method='vonmises',\n                 block_sparse_reduction=True,\n                 number_of_cells=15**3)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Run the simulation over 5 units of time and save the azimuthal angle of the mean direction of motion defined by: \n\n\\begin{align}\\varphi = \\mathrm{arg}(\\Omega^1+i\\Omega^2) \\in [0,2\\pi],\\end{align}\n\nwhere $\\Omega = (\\Omega^1,\\Omega^2,\\Omega^3)$ is the mean direction of motion of the particles with velocities $(\\Omega_i)_{1\\leq i\\leq N}$ : \n\n\\begin{align}\\Omega := \\frac{\\sum_{i=1}^N \\Omega_i}{|\\sum_{i=1}^N \\Omega_i|}\\end{align}\n\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "frames = [5.]\n\ns = time.time()\ndata = save(simu,frames,[],[\"phi\"],save_file=False)\ne = time.time()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Print the total simulation time and the average time per iteration. \n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "print('Total time: '+str(e-s)+' seconds')\nprint('Average time per iteration: '+str((e-s)/simu.iteration)+' seconds')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Plot the azimuthal angle $\\varphi$. \n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "plt.plot(data[\"time\"],data[\"phi\"])\nplt.xlabel(\"time\")\nplt.ylabel(\"Azimuthal angle\")"
      ]
    }
  ],
  "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.8.5"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}