Analytic potential.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Using matplotlib backend: TkAgg\n"
     ]
    }
   ],
   "source": [
    "# %pylab\n",
    "\n",
    "import numpy\n",
    "import matplotlib\n",
    "from matplotlib import pylab, mlab, pyplot\n",
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "np = numpy\n",
    "plt = pyplot\n",
    "\n",
    "from pylab import *\n",
    "from numpy import *\n",
    "\n",
    "%matplotlib"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [
    {
     "ename": "Exception",
     "evalue": "File `'/home/adrian/Documents/Programmes/Python/GiantMol/Functions/turbo_colormap.py'` not found.",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mOSError\u001b[0m                                   Traceback (most recent call last)",
      "\u001b[0;32m/usr/lib/python3/dist-packages/IPython/core/magics/execution.py\u001b[0m in \u001b[0;36mrun\u001b[0;34m(self, parameter_s, runner, file_finder)\u001b[0m\n\u001b[1;32m    713\u001b[0m             \u001b[0mfpath\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0marg_lst\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 714\u001b[0;31m             \u001b[0mfilename\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mfile_finder\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfpath\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    715\u001b[0m         \u001b[0;32mexcept\u001b[0m \u001b[0mIndexError\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/usr/lib/python3/dist-packages/IPython/utils/path.py\u001b[0m in \u001b[0;36mget_py_filename\u001b[0;34m(name, force_win32)\u001b[0m\n\u001b[1;32m    108\u001b[0m     \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 109\u001b[0;31m         \u001b[0;32mraise\u001b[0m \u001b[0mIOError\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'File `%r` not found.'\u001b[0m \u001b[0;34m%\u001b[0m \u001b[0mname\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    110\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mOSError\u001b[0m: File `'/home/adrian/Documents/Programmes/Python/GiantMol/Functions/turbo_colormap.py'` not found.",
      "\nDuring handling of the above exception, another exception occurred:\n",
      "\u001b[0;31mException\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[0;32m/tmp/ipykernel_217862/305484434.py\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mget_ipython\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mrun_line_magic\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'run'\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m'/home/adrian/Documents/Programmes/Python/GiantMol/Functions/turbo_colormap.py'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;32m/usr/lib/python3/dist-packages/IPython/core/interactiveshell.py\u001b[0m in \u001b[0;36mrun_line_magic\u001b[0;34m(self, magic_name, line, _stack_depth)\u001b[0m\n\u001b[1;32m   2362\u001b[0m                 \u001b[0mkwargs\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'local_ns'\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mget_local_scope\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mstack_depth\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   2363\u001b[0m             \u001b[0;32mwith\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mbuiltin_trap\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 2364\u001b[0;31m                 \u001b[0mresult\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mfn\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   2365\u001b[0m             \u001b[0;32mreturn\u001b[0m \u001b[0mresult\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   2366\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m<decorator-gen-52>\u001b[0m in \u001b[0;36mrun\u001b[0;34m(self, parameter_s, runner, file_finder)\u001b[0m\n",
      "\u001b[0;32m/usr/lib/python3/dist-packages/IPython/core/magic.py\u001b[0m in \u001b[0;36m<lambda>\u001b[0;34m(f, *a, **k)\u001b[0m\n\u001b[1;32m    185\u001b[0m     \u001b[0;31m# but it's overkill for just that one bit of state.\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    186\u001b[0m     \u001b[0;32mdef\u001b[0m \u001b[0mmagic_deco\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0marg\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 187\u001b[0;31m         \u001b[0mcall\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;32mlambda\u001b[0m \u001b[0mf\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m*\u001b[0m\u001b[0ma\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mk\u001b[0m\u001b[0;34m:\u001b[0m \u001b[0mf\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0ma\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mk\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    188\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    189\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0mcallable\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0marg\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/usr/lib/python3/dist-packages/IPython/core/magics/execution.py\u001b[0m in \u001b[0;36mrun\u001b[0;34m(self, parameter_s, runner, file_finder)\u001b[0m\n\u001b[1;32m    723\u001b[0m             \u001b[0;32mif\u001b[0m \u001b[0mos\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mname\u001b[0m \u001b[0;34m==\u001b[0m \u001b[0;34m'nt'\u001b[0m \u001b[0;32mand\u001b[0m \u001b[0mre\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mmatch\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34mr\"^'.*'$\"\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mfpath\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    724\u001b[0m                 \u001b[0mwarn\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'For Windows, use double quotes to wrap a filename: %run \"mypath\\\\myfile.py\"'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 725\u001b[0;31m             \u001b[0;32mraise\u001b[0m \u001b[0mException\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mmsg\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    726\u001b[0m         \u001b[0;32mexcept\u001b[0m \u001b[0mTypeError\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    727\u001b[0m             \u001b[0;32mif\u001b[0m \u001b[0mfpath\u001b[0m \u001b[0;32min\u001b[0m \u001b[0msys\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mmeta_path\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mException\u001b[0m: File `'/home/adrian/Documents/Programmes/Python/GiantMol/Functions/turbo_colormap.py'` not found."
     ]
    }
   ],
   "source": [
    "%run /home/adrian/Documents/Programmes/Python/GiantMol/Functions/turbo_colormap.py"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "$\\kappa_{simion}$ = 0.270471\n"
     ]
    }
   ],
   "source": [
    "# Constantes de la physique\n",
    "# ABSOLUMENT RECHARGER APRÈS AVOIR EXECUTÉ LES CASES D'IMPORT AU DESSUS\n",
    "\n",
    "C_e = 1.602e-19        # Coulomb\n",
    "kb = 1.38064852*1e-23  # Boltzman\n",
    "m_Ca = 40.078*1.66054e-27 # masse Ca 40.078\n",
    "m_GM = 1e6*1.66054e-27 # la masse de la GMol\n",
    "eps0 = 8.854187*1e-12  # permittivité électrique du vide\n",
    "\n",
    "\n",
    "r0 = 2.5e-3 # 2.5e-3   # rayon piège Ca+\n",
    "d0 = 4e-3/2            # longueur piège Ca+\n",
    "Omega = 2.045e6*2*pi # 2.047e6\n",
    "bk = 4 # nombre de barreaux par groupe (2 -> 4-pole , 4 -> 8-pole ...)\n",
    "\n",
    "mkappa = 0.23          # écrantage piège réel GiantMol\n",
    "wzLC = (2*pi*90806.9982303)**2\n",
    "kappa_simion = m_Ca*d0**2*wzLC/(2*C_e)\n",
    "print('%s = %f' % ('$\\kappa_{simion}$',kappa_simion) )\n",
    "\n",
    "zeta = kappa_simion*r0**2/d0**2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.25 mm\n"
     ]
    }
   ],
   "source": [
    "Mk = 40\n",
    "Ml = Mk\n",
    "\n",
    "DL = 1e-2/Mk\n",
    "DLz = 2*DL\n",
    "Mp = int(5e-3/DL//2)\n",
    "\n",
    "print(DL*1e3,'mm')\n",
    "\n",
    "# Nz = 50\n",
    "\n",
    "stgart = 1\n",
    "end = 20\n",
    "delta_reel = end - stgart\n",
    "Ec = zeros((2*Mk,2*Ml,delta_reel,3))\n",
    "Ed = zeros((2*Mk,2*Ml,delta_reel,3))\n",
    "\n",
    "rho = 0.001"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "start\n",
      "0 -40 1\n",
      "0 -40 2\n",
      "0 -40 3\n",
      "0 -40 4\n",
      "0 -40 5\n",
      "0 -40 6\n",
      "0 -40 7\n",
      "0 -40 8\n",
      "0 -40 9\n",
      "0 -40 10\n",
      "0 -40 11\n",
      "0 -40 12\n",
      "0 -40 13\n",
      "0 -40 14\n",
      "0 -40 15\n",
      "0 -40 16\n",
      "0 -40 17\n",
      "0 -40 18\n",
      "0 -40 19\n",
      "0 -39 1\n",
      "0 -39 2\n",
      "0 -39 3\n",
      "0 -39 4\n",
      "0 -39 5\n",
      "0 -39 6\n",
      "0 -39 7\n",
      "0 -39 8\n",
      "0 -39 9\n",
      "0 -39 10\n",
      "0 -39 11\n",
      "0 -39 12\n",
      "0 -39 13\n",
      "0 -39 14\n",
      "0 -39 15\n",
      "0 -39 16\n",
      "0 -39 17\n",
      "0 -39 18\n",
      "0 -39 19\n",
      "0 -38 1\n",
      "0 -38 2\n",
      "0 -38 3\n",
      "0 -38 4\n",
      "0 -38 5\n",
      "0 -38 6\n",
      "0 -38 7\n",
      "0 -38 8\n",
      "0 -38 9\n",
      "0 -38 10\n",
      "0 -38 11\n",
      "0 -38 12\n",
      "0 -38 13\n",
      "0 -38 14\n",
      "0 -38 15\n",
      "0 -38 16\n",
      "0 -38 17\n",
      "0 -38 18\n",
      "0 -38 19\n",
      "0 -37 1\n",
      "0 -37 2\n",
      "0 -37 3\n",
      "0 -37 4\n",
      "0 -37 5\n",
      "0 -37 6\n",
      "0 -37 7\n",
      "0 -37 8\n",
      "0 -37 9\n",
      "0 -37 10\n",
      "0 -37 11\n",
      "0 -37 12\n",
      "0 -37 13\n",
      "0 -37 14\n",
      "0 -37 15\n",
      "0 -37 16\n",
      "0 -37 17\n",
      "0 -37 18\n",
      "0 -37 19\n",
      "0 -36 1\n",
      "0 -36 2\n",
      "0 -36 3\n",
      "0 -36 4\n",
      "0 -36 5\n",
      "0 -36 6\n",
      "0 -36 7\n",
      "0 -36 8\n",
      "0 -36 9\n",
      "0 -36 10\n",
      "0 -36 11\n",
      "0 -36 12\n",
      "0 -36 13\n",
      "0 -36 14\n",
      "0 -36 15\n",
      "0 -36 16\n",
      "0 -36 17\n",
      "0 -36 18\n",
      "0 -36 19\n",
      "0 -35 1\n",
      "0 -35 2\n",
      "0 -35 3\n",
      "0 -35 4\n",
      "0 -35 5\n",
      "0 -35 6\n",
      "0 -35 7\n",
      "0 -35 8\n",
      "0 -35 9\n",
      "0 -35 10\n",
      "0 -35 11\n",
      "0 -35 12\n",
      "0 -35 13\n",
      "0 -35 14\n",
      "0 -35 15\n",
      "0 -35 16\n",
      "0 -35 17\n",
      "0 -35 18\n",
      "0 -35 19\n",
      "0 -34 1\n",
      "0 -34 2\n",
      "0 -34 3\n",
      "0 -34 4\n",
      "0 -34 5\n",
      "0 -34 6\n",
      "0 -34 7\n",
      "0 -34 8\n",
      "0 -34 9\n",
      "0 -34 10\n",
      "0 -34 11\n",
      "0 -34 12\n",
      "0 -34 13\n",
      "0 -34 14\n",
      "0 -34 15\n",
      "0 -34 16\n",
      "0 -34 17\n",
      "0 -34 18\n",
      "0 -34 19\n",
      "0 -33 1\n",
      "0 -33 2\n",
      "0 -33 3\n",
      "0 -33 4\n",
      "0 -33 5\n",
      "0 -33 6\n",
      "0 -33 7\n",
      "0 -33 8\n",
      "0 -33 9\n",
      "0 -33 10\n",
      "0 -33 11\n",
      "0 -33 12\n",
      "0 -33 13\n",
      "0 -33 14\n",
      "0 -33 15\n",
      "0 -33 16\n",
      "0 -33 17\n",
      "0 -33 18\n",
      "0 -33 19\n",
      "0 -32 1\n",
      "0 -32 2\n",
      "0 -32 3\n",
      "0 -32 4\n",
      "0 -32 5\n",
      "0 -32 6\n",
      "0 -32 7\n",
      "0 -32 8\n",
      "0 -32 9\n",
      "0 -32 10\n",
      "0 -32 11\n",
      "0 -32 12\n",
      "0 -32 13\n",
      "0 -32 14\n",
      "0 -32 15\n",
      "0 -32 16\n",
      "0 -32 17\n",
      "0 -32 18\n",
      "0 -32 19\n",
      "0 -31 1\n",
      "0 -31 2\n",
      "0 -31 3\n",
      "0 -31 4\n",
      "0 -31 5\n",
      "0 -31 6\n",
      "0 -31 7\n",
      "0 -31 8\n",
      "0 -31 9\n",
      "0 -31 10\n",
      "0 -31 11\n",
      "0 -31 12\n",
      "0 -31 13\n",
      "0 -31 14\n",
      "0 -31 15\n",
      "0 -31 16\n",
      "0 -31 17\n",
      "0 -31 18\n",
      "0 -31 19\n",
      "0 -30 1\n",
      "0 -30 2\n",
      "0 -30 3\n",
      "0 -30 4\n",
      "0 -30 5\n",
      "0 -30 6\n",
      "0 -30 7\n",
      "0 -30 8\n",
      "0 -30 9\n",
      "0 -30 10\n",
      "0 -30 11\n",
      "0 -30 12\n",
      "0 -30 13\n",
      "0 -30 14\n",
      "0 -30 15\n",
      "0 -30 16\n",
      "0 -30 17\n",
      "0 -30 18\n",
      "0 -30 19\n",
      "0 -29 1\n",
      "0 -29 2\n",
      "0 -29 3\n",
      "0 -29 4\n",
      "0 -29 5\n",
      "0 -29 6\n",
      "0 -29 7\n",
      "0 -29 8\n",
      "0 -29 9\n",
      "0 -29 10\n",
      "0 -29 11\n",
      "0 -29 12\n",
      "0 -29 13\n",
      "0 -29 14\n",
      "0 -29 15\n",
      "0 -29 16\n",
      "0 -29 17\n",
      "0 -29 18\n",
      "0 -29 19\n",
      "0 -28 1\n",
      "0 -28 2\n",
      "0 -28 3\n",
      "0 -28 4\n",
      "0 -28 5\n",
      "0 -28 6\n",
      "0 -28 7\n",
      "0 -28 8\n",
      "0 -28 9\n",
      "0 -28 10\n",
      "0 -28 11\n",
      "0 -28 12\n",
      "0 -28 13\n",
      "0 -28 14\n",
      "0 -28 15\n",
      "0 -28 16\n",
      "0 -28 17\n",
      "0 -28 18\n",
      "0 -28 19\n",
      "0 -27 1\n",
      "0 -27 2\n",
      "0 -27 3\n",
      "0 -27 4\n",
      "0 -27 5\n",
      "0 -27 6\n",
      "0 -27 7\n",
      "0 -27 8\n",
      "0 -27 9\n",
      "0 -27 10\n",
      "0 -27 11\n",
      "0 -27 12\n",
      "0 -27 13\n",
      "0 -27 14\n",
      "0 -27 15\n",
      "0 -27 16\n",
      "0 -27 17\n",
      "0 -27 18\n",
      "0 -27 19\n",
      "0 -26 1\n",
      "0 -26 2\n",
      "0 -26 3\n",
      "0 -26 4\n",
      "0 -26 5\n",
      "0 -26 6\n",
      "0 -26 7\n",
      "0 -26 8\n",
      "0 -26 9\n",
      "0 -26 10\n",
      "0 -26 11\n",
      "0 -26 12\n",
      "0 -26 13\n",
      "0 -26 14\n",
      "0 -26 15\n",
      "0 -26 16\n",
      "0 -26 17\n",
      "0 -26 18\n",
      "0 -26 19\n",
      "0 -25 1\n",
      "0 -25 2\n",
      "0 -25 3\n",
      "0 -25 4\n",
      "0 -25 5\n",
      "0 -25 6\n",
      "0 -25 7\n",
      "0 -25 8\n",
      "0 -25 9\n",
      "0 -25 10\n",
      "0 -25 11\n",
      "0 -25 12\n",
      "0 -25 13\n",
      "0 -25 14\n",
      "0 -25 15\n",
      "0 -25 16\n",
      "0 -25 17\n",
      "0 -25 18\n",
      "0 -25 19\n",
      "0 -24 1\n",
      "0 -24 2\n",
      "0 -24 3\n",
      "0 -24 4\n",
      "0 -24 5\n",
      "0 -24 6\n",
      "0 -24 7\n",
      "0 -24 8\n",
      "0 -24 9\n",
      "0 -24 10\n",
      "0 -24 11\n",
      "0 -24 12\n",
      "0 -24 13\n",
      "0 -24 14\n",
      "0 -24 15\n",
      "0 -24 16\n",
      "0 -24 17\n",
      "0 -24 18\n",
      "0 -24 19\n",
      "0 -23 1\n",
      "0 -23 2\n",
      "0 -23 3\n",
      "0 -23 4\n",
      "0 -23 5\n",
      "0 -23 6\n",
      "0 -23 7\n",
      "0 -23 8\n",
      "0 -23 9\n",
      "0 -23 10\n",
      "0 -23 11\n",
      "0 -23 12\n",
      "0 -23 13\n",
      "0 -23 14\n",
      "0 -23 15\n",
      "0 -23 16\n",
      "0 -23 17\n",
      "0 -23 18\n",
      "0 -23 19\n",
      "0 -22 1\n",
      "0 -22 2\n",
      "0 -22 3\n",
      "0 -22 4\n",
      "0 -22 5\n",
      "0 -22 6\n",
      "0 -22 7\n",
      "0 -22 8\n",
      "0 -22 9\n",
      "0 -22 10\n",
      "0 -22 11\n",
      "0 -22 12\n",
      "0 -22 13\n",
      "0 -22 14\n",
      "0 -22 15\n",
      "0 -22 16\n",
      "0 -22 17\n",
      "0 -22 18\n",
      "0 -22 19\n",
      "0 -21 1\n",
      "0 -21 2\n",
      "0 -21 3\n",
      "0 -21 4\n",
      "0 -21 5\n",
      "0 -21 6\n",
      "0 -21 7\n",
      "0 -21 8\n",
      "0 -21 9\n",
      "0 -21 10\n",
      "0 -21 11\n",
      "0 -21 12\n",
      "0 -21 13\n",
      "0 -21 14\n",
      "0 -21 15\n",
      "0 -21 16\n",
      "0 -21 17\n",
      "0 -21 18\n",
      "0 -21 19\n",
      "0 -20 1\n",
      "0 -20 2\n",
      "0 -20 3\n",
      "0 -20 4\n",
      "0 -20 5\n",
      "0 -20 6\n",
      "0 -20 7\n",
      "0 -20 8\n",
      "0 -20 9\n",
      "0 -20 10\n",
      "0 -20 11\n",
      "0 -20 12\n",
      "0 -20 13\n",
      "0 -20 14\n",
      "0 -20 15\n",
      "0 -20 16\n",
      "0 -20 17\n",
      "0 -20 18\n",
      "0 -20 19\n",
      "0 -19 1\n",
      "0 -19 2\n",
      "0 -19 3\n",
      "0 -19 4\n",
      "0 -19 5\n",
      "0 -19 6\n",
      "0 -19 7\n",
      "0 -19 8\n",
      "0 -19 9\n",
      "0 -19 10\n",
      "0 -19 11\n",
      "0 -19 12\n",
      "0 -19 13\n",
      "0 -19 14\n",
      "0 -19 15\n",
      "0 -19 16\n",
      "0 -19 17\n",
      "0 -19 18\n",
      "0 -19 19\n",
      "0 -18 1\n",
      "0 -18 2\n",
      "0 -18 3\n",
      "0 -18 4\n",
      "0 -18 5\n",
      "0 -18 6\n",
      "0 -18 7\n",
      "0 -18 8\n",
      "0 -18 9\n",
      "0 -18 10\n",
      "0 -18 11\n",
      "0 -18 12\n",
      "0 -18 13\n",
      "0 -18 14\n",
      "0 -18 15\n",
      "0 -18 16\n",
      "0 -18 17\n",
      "0 -18 18\n",
      "0 -18 19\n",
      "0 -17 1\n",
      "0 -17 2\n",
      "0 -17 3\n",
      "0 -17 4\n",
      "0 -17 5\n",
      "0 -17 6\n",
      "0 -17 7\n",
      "0 -17 8\n",
      "0 -17 9\n",
      "0 -17 10\n",
      "0 -17 11\n",
      "0 -17 12\n",
      "0 -17 13\n",
      "0 -17 14\n",
      "0 -17 15\n",
      "0 -17 16\n",
      "0 -17 17\n",
      "0 -17 18\n",
      "0 -17 19\n",
      "0 -16 1\n",
      "0 -16 2\n",
      "0 -16 3\n",
      "0 -16 4\n",
      "0 -16 5\n",
      "0 -16 6\n",
      "0 -16 7\n",
      "0 -16 8\n",
      "0 -16 9\n",
      "0 -16 10\n",
      "0 -16 11\n",
      "0 -16 12\n",
      "0 -16 13\n",
      "0 -16 14\n",
      "0 -16 15\n",
      "0 -16 16\n",
      "0 -16 17\n",
      "0 -16 18\n"
     ]
    },
    {
     "ename": "KeyboardInterrupt",
     "evalue": "",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mKeyboardInterrupt\u001b[0m                         Traceback (most recent call last)",
      "\u001b[0;32m/tmp/ipykernel_217862/3226942554.py\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m     13\u001b[0m                     \u001b[0mEc\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mx\u001b[0m\u001b[0;34m+\u001b[0m\u001b[0mMk\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0my\u001b[0m\u001b[0;34m+\u001b[0m\u001b[0mMl\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mz\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0mstgart\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m+=\u001b[0m \u001b[0mrho\u001b[0m\u001b[0;34m/\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m4\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0mpi\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0meps0\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m*\u001b[0m \u001b[0mDL\u001b[0m\u001b[0;34m**\u001b[0m\u001b[0;36m2\u001b[0m \u001b[0;34m*\u001b[0m \u001b[0mu_sur_l\u001b[0m \u001b[0;34m*\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0mx\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0mDL\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0mk\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0mDL\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     14\u001b[0m                     \u001b[0mEc\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mx\u001b[0m\u001b[0;34m+\u001b[0m\u001b[0mMk\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0my\u001b[0m\u001b[0;34m+\u001b[0m\u001b[0mMl\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mz\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0mstgart\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m+=\u001b[0m \u001b[0mrho\u001b[0m\u001b[0;34m/\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m4\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0mpi\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0meps0\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m*\u001b[0m \u001b[0mDL\u001b[0m\u001b[0;34m**\u001b[0m\u001b[0;36m2\u001b[0m \u001b[0;34m*\u001b[0m \u001b[0mu_sur_l\u001b[0m \u001b[0;34m*\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0my\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0mDL\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0ml\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0mDL\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 15\u001b[0;31m                     \u001b[0mEc\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mx\u001b[0m\u001b[0;34m+\u001b[0m\u001b[0mMk\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0my\u001b[0m\u001b[0;34m+\u001b[0m\u001b[0mMl\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mz\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0mstgart\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m2\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m+=\u001b[0m \u001b[0mrho\u001b[0m\u001b[0;34m/\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m4\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0mpi\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0meps0\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m*\u001b[0m \u001b[0mDL\u001b[0m\u001b[0;34m**\u001b[0m\u001b[0;36m2\u001b[0m \u001b[0;34m*\u001b[0m \u001b[0mu_sur_l\u001b[0m \u001b[0;34m*\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0mz\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0mDLz\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     16\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     17\u001b[0m             \u001b[0;32mfor\u001b[0m \u001b[0mk\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mrange\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0mMp\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m+\u001b[0m\u001b[0mMp\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mKeyboardInterrupt\u001b[0m: "
     ]
    }
   ],
   "source": [
    "print('start')\n",
    "\n",
    "u_sur_l = 0\n",
    "\n",
    "for x in [0]: # range(0)\n",
    "    for y in range(-Ml,Ml):\n",
    "        for z in range(stgart,end):\n",
    "            print(x,y,z)\n",
    "            for k in range(-Mk,+Mk):\n",
    "                for l in range(-Ml,Ml):\n",
    "#                     print(k,l)\n",
    "                    u_sur_l = ((x*DL-k*DL)**2 + (y*DL-l*DL)**2 + (z*DLz)**2)**(-3/2)\n",
    "                    Ec[x+Mk,y+Ml,z-stgart,0] += rho/(4*pi*eps0) * DL**2 * u_sur_l * (x*DL-k*DL)\n",
    "                    Ec[x+Mk,y+Ml,z-stgart,1] += rho/(4*pi*eps0) * DL**2 * u_sur_l * (y*DL-l*DL)\n",
    "                    Ec[x+Mk,y+Ml,z-stgart,2] += rho/(4*pi*eps0) * DL**2 * u_sur_l * (z*DLz)\n",
    "                \n",
    "            for k in range(-Mp,+Mp):\n",
    "                Mq = int(sqrt(Mp**2-k**2))\n",
    "                for l in range(-Mq,Mq):\n",
    "                    u_sur_l = ((x*DL-k*DL)**2 + (y*DL-l*DL)**2 + (z*DLz)**2)**(-3/2)\n",
    "                    Ed[x+Mk,y+Ml,z-stgart,0] += rho/(4*pi*eps0) * DL**2 * u_sur_l * (x*DL-k*DL)\n",
    "                    Ed[x+Mk,y+Ml,z-stgart,1] += rho/(4*pi*eps0) * DL**2 * u_sur_l * (y*DL-l*DL)\n",
    "                    Ed[x+Mk,y+Ml,z-stgart,2] += rho/(4*pi*eps0) * DL**2 * u_sur_l * (z*DLz)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# chjeck cercle\n",
    "fig = plt.figure('chjeck cercle',clear='True')\n",
    "ax = fig.add_subplot(111)\n",
    "\n",
    "for k in range(-Mp,+Mp):\n",
    "    Mq = sqrt(Mp**2-k**2)\n",
    "    ax.plot(k,Mq,'+')\n",
    "    print(k**2+Mq**2)\n",
    "ax.set_aspect('equal', adjustable='box')\n",
    "ax.grid()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(shape(Ec))\n",
    "E = Ed - Ec\n",
    "print(shape(E))\n",
    "E = E**2\n",
    "E = sum(E,axis=3)\n",
    "E = sqrt(E)\n",
    "print(shape(E))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X = np.arange(stgart*DL,end*DL,DL)\n",
    "Y = np.arange(-Ml*DL,Ml*DL,DL)\n",
    "print(len(X),len(Y))\n",
    "print(shape(Z),shape(E))\n",
    "Z = E[40,:,:]*1e-6\n",
    "\n",
    "fig, ax = plt.subplots(num='coucou',clear='True')\n",
    "im = ax.contourf(X,Y,Z,20,cmap='turbo', interpolation='bilinear') # ,cmap=cm.gray\n",
    "CS = ax.contour(X, Y, Z,20)\n",
    "ax.vlines(0,-Mk*DL,-Mp*DL)\n",
    "ax.vlines(0,+Mk*DL,+Mp*DL)\n",
    "ax.clabel(CS, inline=True, fontsize=10)\n",
    "ax.set_title('Simplest default with labels')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "E = Ed\n",
    "print(shape(E))\n",
    "E = E**2\n",
    "E = sum(E,axis=3)\n",
    "E = sqrt(E)\n",
    "print(shape(E))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X = np.arange(stgart*DL,end*DL,DL)\n",
    "Y = np.arange(-Ml*DL,Ml*DL,DL)\n",
    "print(len(X),len(Y))\n",
    "print(shape(Z),shape(E))\n",
    "Z = E[40,:,:]\n",
    "\n",
    "fig, ax = plt.subplots(num='coucou 2',clear='True')\n",
    "im = ax.contourf(X,Y,Z,20,cmap='turbo', interpolation='bilinear') # ,cmap=cm.gray\n",
    "CS = ax.contour(X, Y, Z,20)\n",
    "ax.vlines(0,-Mk*DL,+Mk*DL)\n",
    "ax.clabel(CS, inline=True, fontsize=10)\n",
    "ax.set_title('Simplest default with labels')"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.10.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}