{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Calculation of Matrix Elements"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We show how to compute matrix elements using the pairinteraction Python API. As an introductory example, we consider Rubidium and calculate the values of the radial matrix element $\\left|\\left\\langle ns_{1/2},m_j=1/2 \\right| r \\left|n'p_{1/2},m_j=1/2\\right\\rangle\\right|$ as a function of the principal quantum numbers $n$ and $n'$. This [Jupyter notebook](https://github.com/pairinteraction/pairinteraction/blob/master/doc/sphinx/examples_python/matrix_elements.ipynb) is available on GitHub."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As described in the [introduction](introduction.ipynb), we start our code with some preparations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "# Arrays\n",
    "import numpy as np\n",
    "\n",
    "# Plotting\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib.ticker import MaxNLocator\n",
    "\n",
    "# Operating system interfaces\n",
    "import os\n",
    "\n",
    "# pairinteraction :-)\n",
    "from pairinteraction import pireal as pi\n",
    "\n",
    "# Create cache for matrix elements\n",
    "if not os.path.exists(\"./cache\"):\n",
    "    os.makedirs(\"./cache\")\n",
    "cache = pi.MatrixElementCache(\"./cache\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We use pairinteraction's `StateOne` class to define the single-atom states $\\left|n,l,j,m_j\\right\\rangle$ for which the matrix elements should be calculated."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "array_n = range(51, 61)\n",
    "array_nprime = range(51, 61)\n",
    "array_state_final = [pi.StateOne(\"Rb\", n, 0, 0.5, 0.5) for n in array_n]\n",
    "array_state_initial = [pi.StateOne(\"Rb\", n, 1, 0.5, 0.5) for n in array_nprime]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The method `MatrixElementCache.getRadial(state_f, state_i, power)` returns the value of the radial matrix element of $r^p$ in units of $\\mu\\text{m}^p$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "matrixelements = np.empty((len(array_state_final), len(array_state_initial)))\n",
    "for idx_f, state_f in enumerate(array_state_final):\n",
    "    for idx_i, state_i in enumerate(array_state_initial):\n",
    "        matrixelements[idx_f, idx_i] = np.abs(cache.getRadial(state_f, state_i, 1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We visualize the calculated matrix elements with `matplotlib`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "tags": [
     "nbsphinx-thumbnail"
    ]
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQ0AAAEKCAYAAAAIFwCwAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAEiJJREFUeJzt3X+wXHV5x/H3E0NyY+SnITEtTIOWoDMqaYxY2xCMWgSqRqVELVV+2MYfSGvHyqA4Fss4o4J1qH+oAaEMRQZUUKZUJDqjju0AJpooSBSRUBIkEBHUKgnJffrHOXFuwk2y35vz3c1m3q+ZO3f37Nlnn9l77+eec/bss5GZSFKvJg26AUnDxdCQVMTQkFTE0JBUxNCQVMTQkFTE0JBUxNCQVMTQkFRk8qAb6MWUmJojTB90G/unqFW2UuFJlf7PHdD9n8Lsox/tvCbAM6L753bdA0+y6dFtPRUeitAYYToviVcMuo3eVfihEnX+WGJSnT/umFznVyumTatSlz+Y2XnJC266tvOaAH8+0v3vwnGveqDndd09kVTE0JBUxNCQVMTQkFTE0JBUxNCQVKRqaETEIRHxxYhYGxF3R8RLI+KwiFgREfe03w+t2YOkbtXe0rgUuCUznwscC9wNnA98IzOPBr7RXpc0JKqFRkQcDCwCPgeQmVsy8zFgCXBVu9pVwOtq9SCpezW3NI4CHgGujIjvR8TlETEdmJWZP2/XeQiYNd6dI2JZRKyMiJVPsrlim5JK1AyNycB84NOZ+SfA/7HTrkg2o9DHHYeemcszc0FmLjiAqRXblFSiZmisB9Zn5u3t9S/ShMjGiJgN0H5/uGIPkjpWLTQy8yHggYg4pl30CuBHwE3AGe2yM4Cv1OpBUvdqv8v1XOCaiJgC/Aw4iyaoro+ItwH3A0sr9yCpQ1VDIzNXAwvGuWmI3ucuaSzPCJVUxNCQVMTQkFTE0JBUxNCQVGQoBgtXU2MAMFQZAuwA4NazZlQpu+zLN3de85gDftd5TYBTXvbWzmves+6qPa/UcktDUhFDQ1IRQ0NSEUNDUhFDQ1IRQ0NSkaovuUbEOuDXwDZga2YuiIhjgc8AzwDWAadn5q9q9iGpO/3Y0licmfMyc/u7XS8Hzs/MFwA3Au/rQw+SOjKI3ZO5wLfbyyuAUwfQg6QJqh0aCdwaEasiYlm77C6aieQApwFHVu5BUodqh8bCzJwPnAycExGLgLOBd0XEKuBAYMt4d3QaubRvqhoambmh/f4wzfGL4zJzbWaemJkvAq4F7t3FfZ1GLu2Dan5Y0vSIOHD7ZeBE4M6ImNkumwR8kOaVFElDouaWxizgOxGxBrgDuDkzbwHeHBE/AdYCDwJXVuxBUseqnaeRmT+j+fzWnZdfSvMZr5KGkGeESipiaEgqYmhIKmJoSCpiaEgqYmhIKjI808hrTA6vMDW8KVuh16c9rfuaQIxUOtt25jOrlH3dDf9dpe7RBzzSec23vvwtndcEYGP3vbL1yZ5XdUtDUhFDQ1IRQ0NSEUNDUhFDQ1IRQ0NSkUFMI59HM0NjBNgKvCsz76jZh6Tu9OM8jcWZuWnM9Y8DH87Mr0bEKe31l/WhD0kdGMTuSQIHtZcPphnEI2lI1N7S2D6NPIHPZuZy4D3A1yLiEprQ+rPKPUjqUO3QWJiZG9q5oCsiYi3wV8A/ZuaXImIp8DnglTvfsf3Ig2UAIzy9cpuSetX3aeTAGcAN7SpfaJeNd1+nkUv7oL5PI6c5hnFCu9rLgXtq9SCpezV3T2YBN0bz7tTJwOcz85aI+A1waURMBp6g3QWRNBwGMY38O8CLaj2upLo8I1RSEUNDUhFDQ1IRQ0NSEUNDUpEhGizcfb5VGQAMVYYAT5o20nlNAGbOqFJ24Rd+WKXuc6ZsrFL3fa/86+6LPrxpz+tMwOjmzZ3XzNHseV23NCQVMTQkFTE0JBUxNCQVMTQkFTE0JBUxNCQVGcQ08uuAY9pVDgEey8x5NfuQ1J2+TyPPzDduvxwRnwAe70MPkjoysDNCo5nOs5RmepekIVH7mMb2aeSr2kHBYx0PbMzMccf9RcSyiFgZESufpPvTZiVNTN+nkWfmt9vb3gxcu6s7th93sBzgoDis9xPjJVU1iGnktPNB3wBcV/PxJXVvENPIofmck7WZub7W40uqo+/TyNvb3sRudk0k7bv6Po28ve3MWo8rqS7PCJVUxNCQVMTQkFTE0JBUxNCQVGQ4ppFHpcnhFaaGA0yaOrX7ooc/s/uawPOuva9K3WdPfbhK3X/9i9dUqZuPdD85PJ+o8/aH3LatQlGnkUuqxNCQVMTQkFTE0JBUxNCQVMTQkFTE0JBUpGpoRMS6iPhhRKyOiJVjlp8bEWsj4q6I+HjNHiR1q+/TyCNiMbAEODYzN7ejACUNiUHsnrwT+GhmbobfjwKUNCQGMY18LnB8RNweEd+KiBePd8cdppGn08ilfUXfp5G3j3kY8KfAi4HrI+LZmTue/L7DNPJJTiOX9hWDmEa+HrghG3cAo8CMmn1I6s4gppF/GVjcLp8LTAG6f4uhpCr6Po08IqYAV0TEncAW4Iydd00k7bv6Po08M7cAf1PrcSXV5RmhkooYGpKKGBqSihgakor0dCA0IqYCpwJzxt4nM/+lTluS9lW9vnryFeBxYBUwgHO6o8rk8CpTwwFmdj85fPbVdd6iM2fkF1XqXn3Soip1RzfV6bfG5PAqU8MBcrRO3R71GhpHZOZJVTuRNBR6PabxPxHxgqqdSBoKvW5pLATOjIj7aHZPAsjMfGG1ziTtk3oNjZOrdiFpaPQUGpl5f+1GJA0Hz9OQVKTvg4Uj4sKI2NAuWx0Rp9TsQVK3+j5YuPXJzLykD48tqWPunkgqMojBwgDvjogfRMQVEXFo5R4kdah2aCzMzPk0L9meExGLgE8DzwHmAT8HPjHeHXecRv5E5TYl9arvg4Uzc2NmbsvMUeAymmHD4913eWYuyMwFB8RIzTYlFej7YOGImD1mtdfTDBuWNCQGMVj46oiYR3O8Yx3w9oo9SOrYIAYLv6XWY0qqz5dcJRUxNCQVMTQkFTE0JBUxNCQVMTQkFenHu1z3WkQQU6Z0X3fGYZ3XBJh6xW86r3nEtF92XhPglr98yqvinRh9ZOc3NncjtzxZp26NyeG1poYP+PPS3dKQVMTQkFTE0JBUxNCQVMTQkFTE0JBUpO/TyMfc9t6IyIiYUbMHSd0ayDTyiDiSZijP//bh8SV1aFC7J58EzqMZxCNpiPR9GnlELAE2ZOaayo8tqYLauycLM3NDRMwEVkTEWuADNLsmu9WGzDKAkZhet0tJPev3NPITgKOANRGxDjgC+F5EPGuc+/5+GvkUp5FL+4x+TyP/bmbOzMw5mTkHWA/Mz8yHavUhqVt9n0Ze8fEk9UHfp5HvtM6cWo8vqQ7PCJVUxNCQVMTQkFTE0JBUxNCQVMTQkFRkKKaRM3kykypMDt9yWYUJ1MBR037Vec2Vr3125zUBRjf9okrdoZoaDnUmhw94angtbmlIKmJoSCpiaEgqYmhIKmJoSCpiaEgqUvUl13bQzq+BbcDWzFwQERcBS4BR4GHgzMx8sGYfkrrTjy2NxZk5LzMXtNcvzswXZuY84D+BD/WhB0kd6fvuSWaOPfNpOk4kl4ZK7TNCt08jT+CzmbkcICI+ArwVeBxYPN4ddxgsPPnAym1K6lXtLY2FmTkfOBk4JyIWAWTmBZl5JHAN8O7x7rjDYOFJT6/cpqRe9Xsa+XE7rXINcGrNHiR1q9/TyO+MiKPHrLYEWFurB0nd6/s08oj4UkQcQ/OS6/3AOyr2IKljfZ9GnpnujkhDzDNCJRUxNCQVMTQkFTE0JBUxNCQVGYrBwlv/KHj0U923etTIo53XBLjvDYd3XnP00Tq9OgB4e13fAtUrtzQkFTE0JBUxNCQVMTQkFTE0JBUxNCQVMTQkFRnENPKLgdcAW4B7gbMy87GafUjqziCmka8Anp+ZLwR+Ary/Dz1I6sggppHfmplb26u3AUf0uwdJE1c7NLZPI1/VThff2dnAV8e7Y0Qsi4iVEbFy6+O/rdqkpN7Vfu/JwszcEBEzgRURsTYzvw0QERcAW2mGCz9F+3EHywGmz53tGwOkfcRAppFHxJnAq4HTM32nkDRMBjGN/CTgPOC1mel+hzRkBjGN/KfAVJrdFYDbMtOJ5NKQGMQ08j+u9ZiS6vOMUElFDA1JRQwNSUUMDUlFDA1JRYZiGvnkGOXQkd91XvfxNz6985pQZ3J4btnSeU2AHK10bp1Tw/dbbmlIKmJoSCpiaEgqYmhIKmJoSCpiaEgqUjU0ImJdRPwwIlZHxMp22WkRcVdEjEbEgj3VkLRv6cd5Goszc9OY63cCbwA+24fHltSxvp/clZl3A7SzNCQNmUEPFpY0ZAY2WHhP2pBZBjAy68CaPUoqMJDBwj3ed3lmLsjMBVMOnlarRUmF+j5YuNbjSeqPmlsas4DvRMQa4A7g5naw8OsjYj3wUuDmiPhaxR4kdWwQg4VvpNlVkTSEPCNUUhFDQ1IRQ0NSEUNDUhFDQ1IRQ0NSkcghmO4cEY8A9/e4+gxg0x7XKlej7jD1Omx1h6nXfaHuH2Xm4b0UHIrQKBERKzOz8zkdNeoOU6/DVneYeh22uu6eSCpiaEgqsj+GxvIhqjtMvQ5b3WHqdajq7nfHNCTVtT9uaUiqaCg+AHpXImId8GtgG7A1MxdExGnAhcDzgOMyc2VHdS8GXgNsAe4FzsrMxzqoexGwBBgFHgbOzMwH97bumNveC1wCHL7TgOeJ9Hoh8HfAI+1qH8jM/+qi14g4FzinXX5zZp63t3Uj4jrgmHaVQ4DHMnNeB3XnAZ8BRoCtwLsy8469rHlsW/MZwDrg9Mz8VWGvhwCXA8+nGbV5NvBj4DpgTlt3aWb+sqTuU2Tm0H61T8KMnZY9j+YX5ZvAgg7rnghMbi9/DPhYR3UPGnP574HPdFG3XX4k8DWac1yecvsEer0Q+KcKP7PFwNeBqe31mV09B2Nu/wTwoY76vRU4ub18CvDNDmp+FzihvXw2cNEEer0K+Nv28hSaoPw4cH677PyJ/N7u/LXf7Z5k5t2Z+eMKdW/NzK3t1duAIzqqO/a/yXSa/xBd+SRwXsc1a3gn8NHM3Ay/Hw/ZmWhG3y8Fru2oZAIHtZcPBoq2DHdhLrB9fu4K4NSSO0fEwcAi4HMAmbklmy3hJTRhQvv9dXvb6LCHRq1p53uqezbw1a7qRsRHIuIB4HTgQ13UjYglwIbMXDOBervsFXh3RPwgIq6IiEM7qjsXOD4ibo+Ib0XEizvsF+B4YGNm3tNR3fcAF7c/s0uA93dQ8y6aP3CA02i2EkscRbPbeGVEfD8iLm/HbM7KzJ+36zxEM1Fv7+ztpsogv4A/bL/PBNYAi8bc9k0mvnuyu7oX0Eweiy7rtsvfD3y4i7rA7cDB7fJ1lO+ejFdzFvA0mn82HwGu6KjXO4FPAUEzfPq+0ud3Dz+zTwPv7ep3Afg34NR2+VLg6x3UfC7Nbs8q4J+BXxTWXEBzfOUl7fVLgYtojuOMXe+XE3kexn4N9ZZG7sW084nUjYgzgVfTHKQq3uTvod9rKNws3UXdE2j+86xpD7odAXwvIp61N71m5sbM3JaZo8Bl4/Q/obrAeuCGbNxBc1B4Rgd1iYjJNJ/od11pr7upewZwQ7vKFyh8Hnbx3K7NzBMz80U0u1H3Fra6Hlifmbe3178IzAc2RsRsgPb7Xu/6DW1o1Jp2vqu6EXESzfGB12bmbzuse/SY1ZYAazuo+93MnJmZczJzDs0v1PzMfGgve509ZrXXU/h87+Zn9mWag6FExFyag3glr/Ts7nfhlcDazFxf0use6j5IE8wALwd63u3ZzXM7s102CfggzSspPWt/tg9ExPZXi14B/Ai4iSbkaL9/paTueIb5JddZwI3txztOBj6f7bRzmk3dw2mmna/OzFd1UPenwFSaD30CuC0z39FB3S+1P+hRmlc5Smrusm5hjV57vbp9uTFpdnne3lHdKcAVEXEnzUvaZxRuye3uOXgTEz8Auqt+fwNc2m7FPEH7oV57WfMfIuKcdp0bgCsn0O+5wDXt8/kz4CyaDYPrI+JtNL9fSydQdweeESqpyNDunkgaDENDUhFDQ1IRQ0NSEUNDUhFDQ1VFxL9HxMsG3Ye6Y2hIKmJoqBMRMSci7o6IyyLiroi4NSKmAY/TnLCl/YQnd6kTETEH+CnNmwRXR8T1wE2Z+R8DbUydc0tDXbovM1e3l1fRTIvSfsbQUJc2j7m8jeF+b5N2wdCQVMTQkFTEA6GSirilIamIoSGpiKEhqYihIamIoSGpiKEhqYihIamIoSGpyP8DgK5E5lbqeQcAAAAASUVORK5CYII=",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10c1a8e48>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig = plt.figure()\n",
    "ax = fig.add_subplot(1, 1, 1)\n",
    "ax.imshow(\n",
    "    matrixelements,\n",
    "    extent=(\n",
    "        array_nprime[0] - 0.5,\n",
    "        array_nprime[-1] + 0.5,\n",
    "        array_n[0] - 0.5,\n",
    "        array_n[-1] + 0.5,\n",
    "    ),\n",
    "    origin=\"lower\",\n",
    ")\n",
    "ax.set_ylabel(r\"n\")\n",
    "ax.set_xlabel(r\"n'\")\n",
    "ax.yaxis.set_major_locator(MaxNLocator(integer=True))\n",
    "ax.xaxis.set_major_locator(MaxNLocator(integer=True));"
   ]
  }
 ],
 "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.10.6 (main, Mar 10 2023, 10:55:28) [GCC 11.3.0]"
  },
  "vscode": {
   "interpreter": {
    "hash": "e7370f93d1d0cde622a1f8e1c04877d8463912d04d973331ad4851f04de6915a"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}