qpms/lepaper/figs/ewald_branchcuts.ipynb

152 lines
43 KiB
Plaintext
Raw Normal View History

{
"cells": [
{
"cell_type": "code",
"execution_count": 123,
"metadata": {},
"outputs": [],
"source": [
"%matplotlib inline\n",
"import numpy as np\n",
"from matplotlib import pyplot as plt, gridspec\n",
"from qpms import EpsMuGenerator, lorentz_drude, unitcell, eV, hbar, BaseSpec, CTMatrix, Particle, pi\n",
"from matplotlib.colors import BoundaryNorm, Normalize\n",
"eh = eV/hbar\n",
"ccycle=plt.rcParams['axes.prop_cycle'].by_key()['color']"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [],
"source": [
"# I only need to build a valid unitcell to use its empty_lattice_modes method\n",
"px = 580e-9\n",
"py = px\n",
"\n",
"a1 = np.array([px, 0])\n",
"a2 = np.array([0, py])\n",
"\n",
"bspec = BaseSpec(lMax=3)\n",
"tm = CTMatrix(bspec, None)\n",
"emg = EpsMuGenerator(lorentz_drude['Au'])\n",
"pp = Particle((0,0), tm, r=50e-9, material=emg)\n",
"u = unitcell(a1, a2, [pp])\n"
]
},
{
"cell_type": "code",
"execution_count": 45,
"metadata": {},
"outputs": [],
"source": [
"fn = 10\n",
"kxs = np.linspace(0, pi/px)\n",
"emptyfreqs = np.empty(kxs.shape+(fn,))\n",
"for i, kx in enumerate(kxs):\n",
" emptyfreqs[i,:] = u.empty_lattice_modes((kx,0),3*eh)[:fn]\n",
"emptykappas = u.omega2k(emptyfreqs)\n",
"mykx = 0.2*pi/px"
]
},
{
"cell_type": "code",
"execution_count": 132,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAY4AAAEOCAYAAACetPCkAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzsnXd8VFX6h58zM+m9FxISSgiEntAJSrEgirJrpQkqAq7guurqT11ABHVV1oqgFEWKBSsqIiiK1AQIECAJkAAJCaSQRnqZzPn9cZMQMJA2k0m5z+dzDJl777nvxHPmO+e+531fIaVERUVFRUWlvmjMbYCKioqKSutCFQ4VFRUVlQahCoeKioqKSoNQhUNFRUVFpUGowqGioqKi0iBU4VBRUVFRaRDNJhxCCGshxH4hRLQQIkYIsbCWc6yEEF8KIRKEEJFCiMDmsk9FRUWdpyr1ozlXHKXAaCllX6AfMFYIMeSqcx4BcqSUXYG3gdeb0T4VFRV1nqrUg2YTDqlQUPmrRWW7OvrwLuDTyn9/DYwRQohmMlFFpd2jzlOV+qBrzpsJIbRAFNAV+EBKGXnVKR2AZAAppV4IcQlwAzKv6mcmMBPAzs4urHv37qY2vRmQUF4C5cWgLwZ9GVSUgUEPsgKkAa4b5S9AVDUNCC1oqpquslmAtvJn5b9LSw3kZBfi4maPlbVFk96BPrcUQ2H5X17X2Fmgc7ZqUt8qTScqKipTSulR13ktdZ5Wja9zGADoiAYpDSRfSqOoooSAgACcnJyadI+2jpSSc+fOkZmZCVCv8VAbzSocUsoKoJ8Qwhn4TgjRS0p5vBH9rABWAAwYMEAePHjQyJa2UAwGyDsPF09AVgLkJEHeBSjMgKJsKL0EZUWgL6kUHANgAMoq21+Rlf+RZKHRVgpMlbhoLTBoLCjDghIsKZJWlEod5cKCcrSUY4FeailHg15q8c96GB12f72JAL/XRpju76JSL4QQSfU5r6XO05Tnd4GEORQCsLRyrJVWlNHlzTFN6rs9kJiYyJQpU8jMzOS+++5j48aN9RoPtdGswlGFlDJXCPEHMBaoOSDPA/5AihBCBzgBWWYwsWWi0YCzv9KCbq77fIMBCtIhNxnyUqAgAwozKCvIJis7h9z8fIqLi8nPL0JIPZaATl+ODj06KtChx4IKrEQJluixpRxHKrCoPK4VV66AUni8djvUdGitkpY0Tw0GA9Igqe2JmJXW0pS3bhN89tlnPPbYYwBs2LCBSZMmsXHjxkb312zCIYTwAMorB6MNcDN/dar9AEwD9gH3AL9LNQtjgzEYDJxIz2dPfCZHknM5l11KZoETl4ptKSn3x1DbX1RKtAaBpU6LlYUGW0sd9lZaHK0tcLa1wMXOEg97KxytLbCx1GJtocFaK7DVGbDV6LHSgsfa82hr61t9+t1qaGnz1GAwsGTJEl5++WV+m/oJfk7eTOOqx57q+LomeXl5PP7446xfv55hw4axYcMGAgMDm9xvc644fIBPK5+faoCNUsqfhBAvAwellD8Aq4F1QogEIBt4oBnta3UYDAaiUy6xJyGT6ORcTl8sJD2/hKLSir98ybfQCmwtdXg6WOHpaI2/iw1dPO3p7u1IiJc9s8cuYcjI7vzfovsabU/2YANFEWl/ed12sHej+1RpdlrMPF2xYgXPPfccubm5aLVaosvO0AEvBl71saWOr9rZu3cvkydPJjk5mYULF/LCCy+g0xnnI7/ZhENKeRToX8vr82v8uwS4t7lsak0YDAYOJ19iy/FUDpzN5kxmIfkl+ivOEYC9tY6unvZ09bSnn78zw7u6E+LjgEZz/Q10ocO6cmD3KQwGQ53nXgvXCUEAFEWmKY+nhDKpq15Xafm0hHm6ceNG5s6dS0ZGBkII7r77bj7++GMcHR3J/j6e6IjzAAQJrTq+akGv1/Pqq6/y8ssv4+/vz65duxg6dKhR72EWH4dK3STnFPFj9AX2JmRyIq2A7MLSKx4x2VpqCfZyoLuPA/38nRnR1Z2uXg6Nvt/A8G7s3R5L7JFz9AoNbHQ/rhOC1Ims0ii2bNnCrFmzSE5ORgjBLbfcwqeffoq39+UVheuEID5MVzZvfTnLuB+GbYEqB/iePXuYMmUKH3zwAY6Ojka/jyocLYT0vBI2RCSxPS6DhIsFlOoN1ccstAI/Fxv6+DkzKtiDW3v6YG9t3P91Q0d256M3LPn56wNNEg4VlYayY8cOZsyYwenTpwEYPnw469ato1OnTma2rHVR0wG+fv16Jk+ebLJ7tX7haKW+84ISPV9FJbPlWCoxqXkUllZUH3O2sWBgoAvhXT0Y38+HDs62JrfH2c2esGFdidhxAilr372iomJMIiMjmT59OidOnAAgLCyMTz/9lJ49e5rZstbF1Q7w9evXm1x0W79wpB6BNzqDTz/ocSf0vges7M1tVa3sP5vFql1nOZCYTU7R5UA5W0stYQHOjO3lzX1h/jjZmmd74cAR3dizPZaYw0nqqkPFZBw5coQHH3yQY8eOAdCrVy8++eQTBgwYYGbLWh/79u1j8uTJJCUl8dJLL/Hiiy8azQF+PVq/cOisoTQfTm9X2k//BGsn8AyBrjdD3wfAqYNZTDMYDGyNTWfdviQOncuhpFx5/GShFXT3dmB0d08mDu6Iv4vpVxT1YfCNSmTv/l2nVOFQMToxMTFMmzaNqKgoAIKDg1m9ejXDhw83s2Wtj9oc4MOGDWu2+4vWHiZRHZGafgKi18OZHZCZoKTtqEJnDa6doOMw6Pk3CBiuBNOZAL3ewJdRyXx5IJnYC3noKz3adlZaBgW68kh4J8KDGhXl3yw8ft8HVFQY+PCbueY2RcXICCGipJRm+Vpvb28vCwuViO8uXbqwYsUKRo8e3ai+opKyAQgLcDWafa2Jmg7wyZMn88EHHzQq1UpTxkPrX3FU4dUdbll8+fe8VIj+HBJ+hYy4y+3gakCAvQd49oQuo6Dn35Vo7EZiMBj4ITqVj3ae5mRafvXuJxdbC27o5sGsGzoT4ts6cugMuiGYz1fsICsjDzdP4+/GUGmfFBYW4uPjw+rVq7ntttua1Fd7FQyAzz//nNmzZyOlNLkD/Hq0nRVHXejL4MRPcGIznI+CSylgqJGQT2sJDr7g3UcRk5A7wc79ul3Gpl7if1tPsTshs3oXlJeDFbf09Gb2yM7N4tQ2NnHRyfxr6kfMnXcnt987yNzmqBgRc644hFDy0/Tv35+PP/6Yfv36Nbqv9rjiyMvLY86cOaxbt46hQ4eyYcOGJjvAmzIe2o9w1EZOIhz7Bs7+qaxGijIrEwNWorUCRx9lZdJpBASPI9/Gl3e3x/P9ofNkFiqJA20ttYzp4ckztwQT4FZLkr9WhJSS+254lV6hgSx41zzfZlRMgzmFIzg4WBoMBhISEgAIDQ3lk08+oU+fPg3u6/6P9gHtJ44jIiKCSZMmkZSUxLx58/jPf/5jFAe4+qiqsbgEwg1PK62K84fhxI+QtFfxleSeg5xEtsee541NcEr6IxFoMNDfsYC5g10ZfeMo0LWNRGtCCLr19CMxPk3dlqtiNBwcHDh48CDbt29n9uzZHDp0iL59+9KvXz9Wr15NaGiouU1scZjbAX492rdw1EaH/kpDcXS//0c8n+45Q265AZD4aXKYqN3Oo2ITlmUG2IXSdNZg7w3uXcG3PwSOgI5DW6Wg9B4QSNTeeFKTs/Ht6GZuc1TaEGPGjCE+Pr5aQI4cOUJYWBi9e/dm9erVDBw40NwmtgiM5QA3Fapw1EJ6Xgkv/RDDr7Hp6A0SjYARQe4suCOkMq3HVKX+Rfw2SNwNqUfhUjJcOge5iZDwG+x8U+lMawW2buDcETx7QIcw5bGXS6AZ3+H1GXxDMGve+5WIP0/w96nqVkkV41MlIL///juzZ8/m2LFjDBo0iJ49e7Jy5Uqj51ZqTbQUB/j1UIWjBnsTMlm8OZbY1HxA8V1MHuDPc2ODsbW86k9l66rEiPS9KjFoThKc/gNS9kNGrOKEL8yA/AuQHAFRn1SeKMDCRhEVJz9w7QLevcA3DHz6gIW16d/wNej
"text/plain": [
"<Figure size 432x288 with 2 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"f = plt.figure()\n",
"gspec = f.add_gridspec(ncols=2, nrows=1, width_ratios=[2,1],\n",
" # height_ratios=heights\n",
" )\n",
"ax = f.add_subplot(gspec[0,1])\n",
"for i in range(fn):\n",
" ax.plot(kxs * px / pi, emptykappas[:,i] * px / pi, c='k')\n",
" ax.plot(mykx * px / pi, \n",
" u.omega2k(u.empty_lattice_modes((mykx,0),3*eh)[i]) * px / pi,\n",
" marker='o', c=ccycle[6])\n",
"ax.axvline(mykx * px / pi, ls='--')\n",
"ax.set_xlim([0,1])\n",
"ax.set_ylim([0,3])\n",
"ax.set_ylabel('$\\\\kappa p_x / \\\\pi$')\n",
"ax.set_xlabel('$k_x p_x / \\\\pi$')\n",
"\n",
"\n",
"# branch cut image\n",
"mybetas = u.omega2k(u.empty_lattice_modes((mykx, 0), 3*eh))[:fn]\n",
"angles = np.linspace(pi,2*pi)\n",
"# branch cuts due to the 'lilgamma' function\n",
"gammacuts[:,:] = mybetas[:,None]*0.5*(1 + np.exp(1j*angles[None,:]))\n",
"# branch cuts due to incomplete Γ\n",
"gci = np.linspace(0, -1) * pi / px\n",
"Gammacuts = 1j*gci[None,:] + (np.sqrt(mybetas[:, None]**2+gci[None,:]**2))\n",
"\n",
"gr = np.linspace(0, 3, num=500) * pi / px\n",
"gi = np.linspace(-1.0, 0,num=400) * pi / px\n",
"gr, gi = np.meshgrid(gr, gi)\n",
"g = gr + 1j*gi\n",
"eps = np.empty(g.shape, dtype=complex)\n",
"for i in np.ndindex(g.shape):\n",
" eps[i] = emg(u.k2omega(g[i])).n#em['eps']\n",
"\n",
"ax = f.add_subplot(gspec[0,0])\n",
"for i in range(fn):\n",
" ax.plot(gammacuts[i].imag * px / pi, gammacuts[i].real * px / pi, c=ccycle[0])\n",
" ax.plot(Gammacuts[i].imag * px / pi, Gammacuts[i].real * px / pi, c=ccycle[1])\n",
" ax.plot(0, \n",
" u.omega2k(u.empty_lattice_modes((mykx,0),3*eh)[i]) * px / pi,\n",
" marker='o', c=ccycle[6])\n",
"ax.contour(g.imag * px / pi, g.real * px / pi, eps.imag)\n",
"#cax = ax.pcolormesh(g.imag * px / pi, g.real * px / pi, eps.imag, cmap = plt.get_cmap('plasma'), norm = Normalize(vmin=-5, vmax=5, clip=True))\n",
"ax.set_xlim([-1,1])\n",
"ax.set_ylim([0,3])\n",
"ax.set_ylabel('$\\\\Re\\\\kappa p_x / \\\\pi$')\n",
"ax.set_xlabel('$\\\\Im\\\\kappa p_x / \\\\pi$')\n",
"f.savefig(\"ewald_branchcuts.pdf\")"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Python [default]",
"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.5.6"
}
},
"nbformat": 4,
"nbformat_minor": 2
}