From 0c442ba7452fa09b14805533914e8a058cfd4b43 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Sun, 14 Jun 2020 21:49:19 +0300 Subject: [PATCH] WIP Ewald 1D in 3D Former-commit-id: e4a1632ef3b0185b67a24bd93fdff8b3de79267c --- ...scuffNotFixed_vs_scuffFixed2-spheres.ipynb | 719 ++++++++++++++++++ notes/ewald_1D_in_3D.lyx | 85 ++- 2 files changed, 798 insertions(+), 6 deletions(-) create mode 100644 misc/eval_TMatrices_generator_vs_scuffNotFixed_vs_scuffFixed2-spheres.ipynb diff --git a/misc/eval_TMatrices_generator_vs_scuffNotFixed_vs_scuffFixed2-spheres.ipynb b/misc/eval_TMatrices_generator_vs_scuffNotFixed_vs_scuffFixed2-spheres.ipynb new file mode 100644 index 0000000..4eab7d8 --- /dev/null +++ b/misc/eval_TMatrices_generator_vs_scuffNotFixed_vs_scuffFixed2-spheres.ipynb @@ -0,0 +1,719 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "import qpms\n", + "import warnings\n", + "from qpms.cybspec import BaseSpec\n", + "from qpms.cytmatrices import CTMatrix, TMatrixGenerator, TMatrixInterpolator\n", + "from qpms.qpms_c import Particle, pgsl_ignore_error\n", + "from qpms.cymaterials import EpsMu, EpsMuGenerator, LorentzDrudeModel, lorentz_drude\n", + "from qpms.cycommon import DebugFlags, dbgmsg_enable\n", + "from qpms import FinitePointGroup, ScatteringSystem, BesselType, eV, hbar\n", + "import scipy.constants as sci\n", + "eh = eV/hbar" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [], + "source": [ + "#TODO\n", + "period = 520e-9\n", + "a1 = np.array([0,period]) \n", + "a2 = np.array([period,0])" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [], + "source": [ + "#Particle positions\n", + "orig_x = [0]\n", + "orig_y = [0]\n", + "orig_xy = np.stack(np.meshgrid(orig_x,orig_y),axis=-1)" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [], + "source": [ + "period = 0.52\n", + "refractive_index = 1.52 # for background medium\n", + "height = 50e-9 # Particle height\n", + "radius = 50e-9 # Particle radius\n", + "medium = EpsMu(refractive_index**2) # non-lossy background medium with constant refr. index #OK\n", + "# global symmetry group of the system\n", + "#sym = FinitePointGroup(point_group_info['D4h'])\n", + "omega = 1.58*eh\n", + "metal = lorentz_drude['Ag']\n", + "kx_lim = np.array([-0.2, 0.2], dtype=float)\n", + "N=501" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "(7.59633723939313, 8.305328715069821, 7.94387469344152, 8.001475225494097)" + ] + }, + "execution_count": 5, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "omega_scuff_min=1.5*eh/3e14\n", + "omega_scuff_max=1.64*eh/3e14\n", + "omega_scuff=omega/3e14\n", + "omega_scuff_slr=(2*np.pi*sci.c/(1.52*0.52e-6))/3e14\n", + "omega_scuff_min, omega_scuff_max, omega_scuff_slr, omega_scuff" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "outputs": [], + "source": [ + "bspec = BaseSpec(lMax = 1)\n", + "\n", + "tmfile_scuffOld = '/home/javier/tmatrices/sphereAg_50nm_oldScuff.TMatrix'\n", + "tmfile_scuffNew = '/home/javier/tmatrices/sphere50nm_newScuff.TMatrix'\n", + "interp_old = TMatrixInterpolator(tmfile_scuffOld, bspec, atol=1e-8) \n", + "interp_new = TMatrixInterpolator(tmfile_scuffNew, bspec, atol=1e-8) \n", + "tmscuff_not_fixed = interp_old(omega)\n", + "tmscuff_bugfixed = interp_new(omega)" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [], + "source": [ + "tmgen = TMatrixGenerator.sphere(medium, metal, radius)" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": {}, + "outputs": [], + "source": [ + "tmgen_omega=tmgen(bspec,omega)" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "array([[-0.07824951+0.25215549j, 0. +0.j , 0. +0.j , 0. +0.j , 0. +0.j , 0. +0.j ],\n", + " [ 0. +0.j , -0.07824951+0.25215549j, 0. +0.j , 0. +0.j , 0. +0.j , 0. +0.j ],\n", + " [ 0. +0.j , 0. +0.j , -0.07824951+0.25215549j, 0. +0.j , 0. +0.j , 0. +0.j ],\n", + " [ 0. +0.j , 0. +0.j , 0. +0.j , -0.00083788-0.01420874j, 0. +0.j , 0. +0.j ],\n", + " [ 0. +0.j , 0. +0.j , 0. +0.j , 0. +0.j , -0.00083788-0.01420874j, 0. +0.j ],\n", + " [ 0. +0.j , 0. +0.j , 0. +0.j , 0. +0.j , 0. +0.j , -0.00083788-0.01420874j]])" + ] + }, + "execution_count": 9, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "tmgen_omega.as_ndarray() # T-Matrix from generator" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "array([[-4.70886671e-03+1.79440815e-02j, 1.60098226e-13+4.20092462e-13j, 2.26455777e-07+3.98932962e-07j, -2.13258227e-14-6.01930824e-14j, 7.66251504e-08-1.32539814e-08j,\n", + " -7.76488937e-14+1.08532726e-13j],\n", + " [-4.41531068e-13+9.81757019e-14j, -4.70582367e-03+1.79389513e-02j, -1.60702588e-13-4.19762769e-13j, 6.32372612e-08+4.71302619e-08j, -2.36039655e-14+1.00983241e-13j,\n", + " 7.76590838e-08-1.38092126e-08j],\n", + " [-4.58579591e-07-7.49053314e-09j, 4.41175187e-13-9.87457125e-14j, -4.70886672e-03+1.79440815e-02j, 2.13790755e-14+1.31812722e-13j, 6.25535562e-08+4.61875295e-08j,\n", + " 4.59276044e-14-4.47400043e-14j],\n", + " [ 4.93366461e-14-4.29235249e-14j, -6.94450011e-09-1.60278760e-08j, -7.79609985e-14+1.07845028e-13j, -3.03166095e-05-2.01474828e-03j, -2.39253108e-13-1.34944425e-13j,\n", + " 1.16328171e-08+6.36863105e-09j],\n", + " [ 8.21008790e-09-1.53591010e-08j, -2.36580203e-14+1.00878213e-13j, -5.41800099e-09-1.65573897e-08j, 2.26431323e-13-1.55248143e-13j, -3.03231670e-05-2.01493037e-03j,\n", + " 2.38800352e-13+1.35401679e-13j],\n", + " [ 2.18731933e-14+1.31428484e-13j, 7.13657440e-09-1.65412452e-08j, -2.48667583e-14-6.12974376e-14j, -1.11447700e-08+7.19733229e-09j, -2.26916383e-13+1.54810715e-13j,\n", + " -3.03166170e-05-2.01474828e-03j]])" + ] + }, + "execution_count": 10, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "tmscuff_not_fixed.as_ndarray() # T-Matrix of not fixed version of Scuff-EM" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "array([[-9.43033280e-02+3.59361755e-01j, 3.20658236e-12+8.41416904e-12j, 5.13892659e-06+8.57797510e-06j, -4.28949536e-13-1.20456867e-12j, 1.41239218e-06-2.46862181e-07j,\n", + " -1.55480479e-12+2.17400571e-12j],\n", + " [-8.84386127e-12+1.96169624e-12j, -9.42395846e-02+3.59248336e-01j, -3.22166474e-12-8.40342496e-12j, 1.16868119e-06+8.68484583e-07j, -4.72874254e-13+2.02296639e-12j,\n", + " 1.43308826e-06-2.58864855e-07j],\n", + " [-9.99885677e-06+6.63229466e-08j, 8.83605486e-12-1.97841571e-12j, -9.43033282e-02+3.59361755e-01j, 4.27888728e-13+2.63976219e-12j, 1.15544356e-06+8.48824647e-07j,\n", + " 9.21093604e-13-8.94225073e-13j],\n", + " [ 9.85267406e-13-8.59818272e-13j, -1.43866720e-06+2.57952214e-07j, -1.56233879e-12+2.15974773e-12j, -6.07143255e-04-4.03488631e-02j, -4.79139425e-12-2.70322093e-12j,\n", + " 1.52945653e-07+8.36194839e-08j],\n", + " [-1.15196465e-06-8.46913318e-07j, -4.73514628e-13+2.02574744e-12j, -1.40806670e-06+2.47097086e-07j, 4.53481568e-12-3.10854698e-12j, -6.07256494e-04-4.03513114e-02j,\n", + " 4.78231690e-12+2.71099900e-12j],\n", + " [ 4.39050064e-13+2.63269196e-12j, -1.17350730e-06-8.70822943e-07j, -4.95069921e-13-1.22744130e-12j, -1.44529174e-07+9.78280211e-08j, -4.54412580e-12+3.10101417e-12j,\n", + " -6.07143303e-04-4.03488631e-02j]])" + ] + }, + "execution_count": 11, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "tmscuff_bugfixed.as_ndarray() # T-Matrix of FIXED version of Scuff-EM" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "(matrix([[-9.43033280e-02+3.59361755e-01j, 3.20658236e-12+8.41416904e-12j],\n", + " [-8.84386127e-12+1.96169624e-12j, -9.42395846e-02+3.59248336e-01j]]),\n", + " matrix([[-9.43033280e-02-3.59361755e-01j, -8.84386127e-12-1.96169624e-12j],\n", + " [ 3.20658236e-12-8.41416904e-12j, -9.42395846e-02-3.59248336e-01j]]))" + ] + }, + "execution_count": 12, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "#I play around with the different operations so everything is correct\n", + "tmscuffnew = tmscuff_bugfixed.as_ndarray()\n", + "tmscuffnew_mat = np.asmatrix(tmscuffnew)\n", + "tmscuffnew_dag = tmscuffnew_mat.getH()\n", + "tmscuffnew_mat[0:2,0:2], tmscuffnew_dag[0:2,0:2]" + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "(matrix([[1, 0],\n", + " [9, 4]]), matrix([[4, 9],\n", + " [0, 1]]))" + ] + }, + "execution_count": 13, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "A = np.matrix([[0,1],[2,3]])\n", + "B = np.matrix([[3,2],[1,0]])\n", + "AB = np.dot(A,B)\n", + "BA = np.dot(B,A)\n", + "AB,BA" + ] + }, + { + "cell_type": "code", + "execution_count": 14, + "metadata": {}, + "outputs": [], + "source": [ + "sum1new = np.dot(tmscuffnew_dag,tmscuffnew) #is this the right order? Regarding the above, yes it is.\n", + "sum2new = 0.5*tmscuffnew.__add__(tmscuffnew_dag) " + ] + }, + { + "cell_type": "code", + "execution_count": 15, + "metadata": {}, + "outputs": [], + "source": [ + "powermatrix_scuffnew = sum1new.__add__(sum2new) #powermatrix for bugfixed scuff" + ] + }, + { + "cell_type": "code", + "execution_count": 16, + "metadata": {}, + "outputs": [], + "source": [ + "#Power matrix for NOT bugfixed scuff\n", + "tmscuffold = tmscuff_not_fixed.as_ndarray()\n", + "tmscuffold_mat = np.asmatrix(tmscuffold)\n", + "tmscuffold_dag = tmscuffold_mat.getH()\n", + "sum1old = np.dot(tmscuffold_dag,tmscuffold) #is this the right order? Regarding the above, yes it is.\n", + "sum2old = 0.5*tmscuffold.__add__(tmscuffold_dag)\n", + "\n", + "powermatrix_scuffold = sum1old.__add__(sum2old) " + ] + }, + { + "cell_type": "code", + "execution_count": 17, + "metadata": {}, + "outputs": [], + "source": [ + "#Power matrix for T matrix generator\n", + "tmscuffgen = tmgen_omega.as_ndarray()\n", + "tmscuffgen_mat = np.asmatrix(tmscuffgen)\n", + "tmscuffgen_dag = tmscuffgen_mat.getH()\n", + "sum1gen = np.dot(tmscuffgen_dag,tmscuffgen) #is this the right order? Regarding the above, yes it is.\n", + "sum2gen = 0.5*tmscuffgen.__add__(tmscuffgen_dag)\n", + "\n", + "powermatrix_scuffgen = sum1gen.__add__(sum2gen) " + ] + }, + { + "cell_type": "code", + "execution_count": 18, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "((-0.008544129580638216+0j),\n", + " (-0.008544129580638216+0j),\n", + " (-0.008544129580638216+0j),\n", + " (-0.0006352854997266949+0j),\n", + " (-0.0006352854997266949+0j),\n", + " (-0.0006352854997266949+0j))" + ] + }, + "execution_count": 18, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "#Power matrix of generator is diagonal and all its eigenvalues are negative:\n", + "powermatrix_scuffgen[0,0], powermatrix_scuffgen[1,1], powermatrix_scuffgen[2,2], powermatrix_scuffgen[3,3], powermatrix_scuffgen[4,4], powermatrix_scuffgen[5,5]" + ] + }, + { + "cell_type": "code", + "execution_count": 19, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "((0.04373066071852283+0j),\n", + " (0.04370088172427582+0j),\n", + " (0.043730660553373296+0j),\n", + " (0.001021256123757952+0j),\n", + " (0.0010213406041430792+0j),\n", + " (0.0010212560750220145+0j))" + ] + }, + "execution_count": 19, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "powermatrix_scuffnew[0,0], powermatrix_scuffnew[1,1], powermatrix_scuffnew[2,2], powermatrix_scuffnew[3,3], powermatrix_scuffnew[4,4], powermatrix_scuffnew[5,5]" + ] + }, + { + "cell_type": "code", + "execution_count": 20, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "((-0.004364703222422517+0j),\n", + " (-0.004361872918812184+0j),\n", + " (-0.004364703232328262+0j),\n", + " (-2.6256479824067564e-05+0j),\n", + " (-2.6262303053141142e-05+0j),\n", + " (-2.6256487319428013e-05+0j))" + ] + }, + "execution_count": 20, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "powermatrix_scuffold[0,0], powermatrix_scuffold[1,1], powermatrix_scuffold[2,2], powermatrix_scuffold[3,3], powermatrix_scuffold[4,4], powermatrix_scuffold[5,5]" + ] + }, + { + "cell_type": "code", + "execution_count": 22, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "((-10.019160178833705-0j),\n", + " (-10.018834234211564-0j),\n", + " (-10.019160118257586-0j),\n", + " (-38.895393845668316-0j),\n", + " (-38.88998623145963-0j),\n", + " (-38.89538088616882-0j))" + ] + }, + "execution_count": 22, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "#So! There is a fixed factor between the N/M elements of Scuff new and old:\n", + "powermatrix_scuffnew[0,0]/powermatrix_scuffold[0,0], powermatrix_scuffnew[1,1]/powermatrix_scuffold[1,1], powermatrix_scuffnew[2,2]/powermatrix_scuffold[2,2], powermatrix_scuffnew[3,3]/powermatrix_scuffold[3,3], powermatrix_scuffnew[4,4]/powermatrix_scuffold[4,4], powermatrix_scuffnew[5,5]/powermatrix_scuffold[5,5]" + ] + }, + { + "cell_type": "code", + "execution_count": 23, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "((-5.118211317583541-0j),\n", + " (-5.11472600126595-0j),\n", + " (-5.118211298254535-0j),\n", + " (-1.6075545942687262-0j),\n", + " (-1.6076875744566317-0j),\n", + " (-1.60755451755371-0j))" + ] + }, + "execution_count": 23, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "#And there is a fixed factor for new / generated as well! (hence also for old / generated)\n", + "powermatrix_scuffnew[0,0]/powermatrix_scuffgen[0,0], powermatrix_scuffnew[1,1]/powermatrix_scuffgen[1,1], powermatrix_scuffnew[2,2]/powermatrix_scuffgen[2,2], powermatrix_scuffnew[3,3]/powermatrix_scuffgen[3,3], powermatrix_scuffnew[4,4]/powermatrix_scuffgen[4,4], powermatrix_scuffnew[5,5]/powermatrix_scuffgen[5,5]" + ] + }, + { + "cell_type": "code", + "execution_count": 21, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "((-1.60755451755371-0j), (-2.557363000632975+0j))" + ] + }, + "execution_count": 21, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "powermatrix_scuffnew[5,5]/powermatrix_scuffgen[5,5], powermatrix_scuffnew[1,1]/powermatrix_scuffgen[1,1]/2" + ] + }, + { + "cell_type": "code", + "execution_count": 22, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "((1.599227713011597e-16-0j),\n", + " (8.902299490051238e-14+0j),\n", + " (1.504486797691096e-21-0j))" + ] + }, + "execution_count": 22, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "#Let's also calculate the determinants\n", + "detnew = powermatrix_scuffnew[0,0]*powermatrix_scuffnew[1,1]*powermatrix_scuffnew[2,2]*powermatrix_scuffnew[3,3]*powermatrix_scuffnew[4,4]*powermatrix_scuffnew[5,5]\n", + "detold = powermatrix_scuffold[0,0]*powermatrix_scuffold[1,1]*powermatrix_scuffold[2,2]*powermatrix_scuffold[3,3]*powermatrix_scuffold[4,4]*powermatrix_scuffold[5,5]\n", + "detgen = powermatrix_scuffgen[0,0]*powermatrix_scuffgen[1,1]*powermatrix_scuffgen[2,2]*powermatrix_scuffgen[3,3]*powermatrix_scuffgen[4,4]*powermatrix_scuffgen[5,5]\n", + "detgen, detnew, detold" + ] + }, + { + "cell_type": "code", + "execution_count": 23, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "((1.5992277130116025e-16+0j),\n", + " (8.902299447416401e-14+2.6128569373609595e-39j),\n", + " (1.504486701753944e-21-2.892571608172536e-45j))" + ] + }, + "execution_count": 23, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "np.linalg.det(powermatrix_scuffgen), np.linalg.det(powermatrix_scuffnew), np.linalg.det(powermatrix_scuffold)" + ] + }, + { + "cell_type": "code", + "execution_count": 24, + "metadata": {}, + "outputs": [], + "source": [ + "#We try to normalize the power matrix elements for each case with the corresponding determinants." + ] + }, + { + "cell_type": "code", + "execution_count": 25, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "((-53426597795433.88+0j),\n", + " (-53426597795433.88+0j),\n", + " (-53426597795433.88+0j),\n", + " (-3972451793812.1055+0j),\n", + " (-3972451793812.1055+0j),\n", + " (-3972451793812.1055+0j))" + ] + }, + "execution_count": 25, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "#T-matrix generator: \n", + "powermatrix_scuffgen[0,0]/detgen, powermatrix_scuffgen[1,1]/detgen, powermatrix_scuffgen[2,2]/detgen, powermatrix_scuffgen[3,3]/detgen, powermatrix_scuffgen[4,4]/detgen, powermatrix_scuffgen[5,5]/detgen " + ] + }, + { + "cell_type": "code", + "execution_count": 26, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "((-2.901124309712079e+18+0j),\n", + " (-2.899243067806416e+18+0j),\n", + " (-2.901124316296214e+18+0j),\n", + " (-1.7452117136795634e+16+0j),\n", + " (-1.7455987711853198e+16+0j),\n", + " (-1.7452122118800434e+16+0j))" + ] + }, + "execution_count": 26, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "#Not bugfixed Scuff-EM: \n", + "powermatrix_scuffold[0,0]/detold, powermatrix_scuffold[1,1]/detold, powermatrix_scuffold[2,2]/detold, powermatrix_scuffold[3,3]/detold, powermatrix_scuffold[4,4]/detold, powermatrix_scuffold[5,5]/detold " + ] + }, + { + "cell_type": "code", + "execution_count": 27, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "((491228819782.9563+0j),\n", + " (490894310768.96173+0j),\n", + " (491228817927.8229+0j),\n", + " (11471823936.043226+0j),\n", + " (11472772908.667904+0j),\n", + " (11471823388.58987+0j))" + ] + }, + "execution_count": 27, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "#New Scuff-EM: \n", + "powermatrix_scuffnew[0,0]/detnew, powermatrix_scuffnew[1,1]/detnew, powermatrix_scuffnew[2,2]/detnew, powermatrix_scuffnew[3,3]/detnew, powermatrix_scuffnew[4,4]/detnew, powermatrix_scuffnew[5,5]/detnew " + ] + }, + { + "cell_type": "code", + "execution_count": 28, + "metadata": {}, + "outputs": [], + "source": [ + "#It might make more sense to renormalize the electric and magnetic parts of the power matrices separately:" + ] + }, + { + "cell_type": "code", + "execution_count": 29, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "((13698.226641508763-0j),\n", + " (13698.226641508763-0j),\n", + " (13698.226641508763-0j),\n", + " (2477776.413504498-0j),\n", + " (2477776.413504498-0j),\n", + " (2477776.413504498-0j))" + ] + }, + "execution_count": 29, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "#T-matrix generator:\n", + "detgen_el = powermatrix_scuffgen[0,0]*powermatrix_scuffgen[1,1]*powermatrix_scuffgen[2,2]\n", + "detgen_mag = powermatrix_scuffgen[3,3]*powermatrix_scuffgen[4,4]*powermatrix_scuffgen[5,5]\n", + "powermatrix_scuffgen[0,0]/detgen_el, powermatrix_scuffgen[1,1]/detgen_el, powermatrix_scuffgen[2,2]/detgen_el, powermatrix_scuffgen[3,3]/detgen_mag, powermatrix_scuffgen[4,4]/detgen_mag, powermatrix_scuffgen[5,5]/detgen_mag " + ] + }, + { + "cell_type": "code", + "execution_count": 30, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "((52525.751093170744-0j),\n", + " (52491.69062780147-0j),\n", + " (52525.75121237855-0j),\n", + " (1450208903.6280527-0j),\n", + " (1450530534.6580672-0j),\n", + " (1450209317.614944-0j))" + ] + }, + "execution_count": 30, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "#NOT bugfixed Scuff-EM:\n", + "detold_el = powermatrix_scuffold[0,0]*powermatrix_scuffold[1,1]*powermatrix_scuffold[2,2]\n", + "detold_mag = powermatrix_scuffold[3,3]*powermatrix_scuffold[4,4]*powermatrix_scuffold[5,5]\n", + "powermatrix_scuffold[0,0]/detold_el, powermatrix_scuffold[1,1]/detold_el, powermatrix_scuffold[2,2]/detold_el, powermatrix_scuffold[3,3]/detold_mag, powermatrix_scuffold[4,4]/detold_mag, powermatrix_scuffold[5,5]/detold_mag " + ] + }, + { + "cell_type": "code", + "execution_count": 31, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "((523.2675016470721+0j),\n", + " (522.9111754524715+0j),\n", + " (523.2674996709441+0j),\n", + " (958726.5381294343+0j),\n", + " (958805.8459399715+0j),\n", + " (958726.4923775065+0j))" + ] + }, + "execution_count": 31, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "#Bugfixed Scuff-EM:\n", + "detnew_el = powermatrix_scuffnew[0,0]*powermatrix_scuffnew[1,1]*powermatrix_scuffnew[2,2]\n", + "detnew_mag = powermatrix_scuffnew[3,3]*powermatrix_scuffnew[4,4]*powermatrix_scuffnew[5,5]\n", + "powermatrix_scuffnew[0,0]/detnew_el, powermatrix_scuffnew[1,1]/detnew_el, powermatrix_scuffnew[2,2]/detnew_el, powermatrix_scuffnew[3,3]/detnew_mag, powermatrix_scuffnew[4,4]/detnew_mag, powermatrix_scuffnew[5,5]/detnew_mag " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "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.6.6" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/notes/ewald_1D_in_3D.lyx b/notes/ewald_1D_in_3D.lyx index 796bf3f..5db7940 100644 --- a/notes/ewald_1D_in_3D.lyx +++ b/notes/ewald_1D_in_3D.lyx @@ -411,7 +411,7 @@ The angular integral (assuming it can be separated from the rest like this) is \begin_inset Formula \[ -I_{l}^{m}\equiv\int\ud\Omega_{\vect r}\,\ushD lm\left(\uvec r\right)e^{i\vect K\cdot\vect r_{\parallel}}e^{-2s_{\bot}r_{\bot}\cos\left(\phi-\Phi\right)} +I_{l}^{m}\equiv\int\ud\Omega_{\vect r}\,\ushD lm\left(\uvec r\right)e^{-\left(r_{\bot}^{2}+2s_{\bot}r_{\bot}\cos\left(\phi-\Phi\right)\right)^{2}\kappa^{2}\gamma_{\vect K}^{2}/4\tau} \] \end_inset @@ -420,23 +420,96 @@ I_{l}^{m}\equiv\int\ud\Omega_{\vect r}\,\ushD lm\left(\uvec r\right)e^{i\vect K\ \end_layout \begin_layout Standard -\begin_inset Note Note -status open +Let's further extract the azimuthal part +\begin_inset Formula $\left(w\equiv2r_{\bot}s_{\bot}\kappa^{2}\gamma_{\vect K}^{2}/4\tau\right)$ +\end_inset + -\begin_layout Plain Layout -which can be separated even more into two integrals \begin_inset Formula \[ -I_{l}^{m}=\lambda'_{lm}\left(\int_{0}^{2\pi}e^{-im\phi}e^{-2s_{\bot}r_{\bot}\cos\left(\phi-\Phi\right)}\ud\phi\right)\left(\int_{0}^{\pi}P_{l}^{-m}\left(\cos\theta\right)e^{i\left|\vect K\right|\left|\vect r\right|\cos\theta}\sin\theta\,\ud\theta\right) +e^{-im\Phi}A_{l}^{m}\equiv\int_{0}^{2\pi}e^{-im\phi}e^{-w\cos\left(\phi-\Phi\right)}\ud\phi=e^{-im\Phi}\int_{0}^{2\pi}e^{-im\varphi}e^{-w\cos\varphi}\ud\varphi \] \end_inset +\begin_inset Formula +\begin{align*} +A_{l}^{m} & =\int_{0}^{2\pi}e^{-im\varphi}\sum_{n=0}^{\infty}\frac{\left(-w\cos\varphi\right)^{n}}{n!}\ud\varphi\\ + & =\int_{0}^{2\pi}e^{-im\varphi}\sum_{n=0}^{\infty}\frac{\left(-w\right)^{n}}{2^{n}n!}\left(e^{i\varphi}+e^{-i\varphi}\right)^{n}\ud\varphi\\ + & =\sum_{n=0}^{\infty}\frac{\left(-w\right)^{n}}{2^{n}n!}\int_{0}^{2\pi}e^{-im\varphi}\sum_{k=0}^{n}\binom{n}{k}e^{ik\varphi}e^{-i\left(n-k\right)\varphi}\ud\varphi\\ + & =\sum_{n=0}^{\infty}\frac{\left(-w\right)^{n}}{2^{n}n!}\sum_{k=0}^{n}\binom{n}{k}\int_{0}^{2\pi}e^{i\left(2k-n-m\right)\varphi}\ud\varphi\\ + & =2\pi\sum_{n=0}^{\infty}\frac{\left(-w\right)^{n}}{2^{n}n!}\sum_{k=0}^{n}\binom{n}{k}\delta_{2k-n-m=0}\\ + & =2\pi\sum_{n=0}^{\infty}\frac{\left(-w\right)^{n}}{2^{n}n!}\sum_{k=0}^{n}\binom{n}{k}\delta_{2k-n-m=0}\\ + & =2\pi\sum_{k=0}^{\infty}\sum_{n=k}^{\infty}\frac{\left(-w\right)^{n}}{2^{n}n!}\binom{n}{k}\delta_{2k-n-m=0}\\ + & =2\pi\sum_{k=0}^{\infty}\frac{\left(-w\right)^{2k-m}}{2^{2k-m}\left(2k-m\right)!}\binom{2k-m}{k}\delta_{2k-m\ge k}\\ + & =2\pi\sum_{k=0}^{\infty}\frac{\left(-w\right)^{2k-m}}{2^{2k-m}}\frac{1}{k!\left(k-m\right)!}\delta_{k-m\ge0}\\ + & =2\pi\sum_{k=\max\left(m,0\right)}^{\infty}\left(-\frac{w}{2}\right)^{2k-m}\frac{1}{k!\left(k-m\right)!} +\end{align*} + +\end_inset + + +\begin_inset Note Note +status collapsed + +\begin_layout Plain Layout +\begin_inset Formula +\begin{align*} +A_{l}^{m} & =\int_{0}^{2\pi}e^{-im\varphi}\sum_{n=0}^{\infty}\frac{\left(-w\cos\varphi\right)^{n}}{n!}\ud\varphi\\ + & =\int_{0}^{2\pi}e^{-im\varphi}\sum_{n=0}^{\infty}\frac{\left(-w\right)^{n}}{2^{n}n!}\left(e^{i\varphi}+e^{-i\varphi}\right)^{n}\ud\varphi\\ + & =\sum_{n=0}^{\infty}\frac{\left(-w\right)^{n}}{2^{n}n!}\int_{0}^{2\pi}e^{-im\varphi}\sum_{k=0}^{n}\binom{n}{k}e^{i\left(n-k\right)\varphi}e^{-ik\varphi}\ud\varphi\\ + & =\sum_{n=0}^{\infty}\frac{\left(-w\right)^{n}}{2^{n}n!}\sum_{k=0}^{n}\binom{n}{k}\int_{0}^{2\pi}e^{i\left(-2k+n-m\right)\varphi}\ud\varphi\\ + & =2\pi\sum_{n=0}^{\infty}\frac{\left(-w\right)^{n}}{2^{n}n!}\sum_{k=0}^{n}\binom{n}{k}\delta_{-2k+n-m=0}\\ + & =2\pi\sum_{n=0}^{\infty}\frac{\left(-w\right)^{n}}{2^{n}n!}\sum_{k=0}^{n}\binom{n}{k}\delta_{-2k+n-m=0}\\ + & =2\pi\sum_{k=0}^{\infty}\sum_{n=k}^{\infty}\frac{\left(-w\right)^{n}}{2^{n}n!}\binom{n}{k}\delta_{-2k+n-m=0}\\ + & =2\pi\sum_{k=0}^{\infty}\frac{\left(-w\right)^{2k+m}}{2^{2k+m}\left(2k+m\right)!}\binom{2k+m}{k}\delta_{2k+m\ge k}\\ + & =2\pi\sum_{k=0}^{\infty}\frac{\left(-w\right)^{2k+m}}{2^{2k+m}}\frac{1}{k!\left(k+m\right)!}\delta_{k+m\ge0}\\ + & =2\pi\sum_{k=\max\left(-m,0\right)}^{\infty}\frac{\left(-w\right)^{2k+m}}{2^{2k+m}}\frac{1}{k!\left(k+m\right)!} +\end{align*} + +\end_inset + + \end_layout \end_inset +Althought it's not superobvious, this sum is symmetric w.r.t. + sign change in +\begin_inset Formula $m$ +\end_inset + +. +\end_layout + +\begin_layout Standard +Let's do the polar integration next: +\begin_inset Formula $r_{\bot}=r\sin\theta$ +\end_inset + + +\begin_inset Formula +\[ +B_{l}^{m}\equiv\int_{0}^{\pi}\sin\theta\ud\theta\,P_{l}^{-m}\left(\cos\theta\right)P_{l}^{0}\left(\cos\theta\right)e^{-\left(\sin\theta\right)^{2}r^{2}\kappa^{2}\gamma_{\vect K}^{2}/4\tau}\left(-\sin\theta\,rs_{\bot}\kappa^{2}\gamma_{\vect K}^{2}/4\tau\right)^{2k-m} +\] + +\end_inset + +Label +\begin_inset Formula $u\equiv r^{2}\kappa^{2}\gamma_{\vect K}^{2}/4\tau,v\equiv rs_{\bot}\kappa^{2}\gamma_{\vect K}^{2}/4\tau$ +\end_inset + +; then +\begin_inset Formula +\begin{align*} +B_{l}^{m} & =\int_{0}^{\pi}\sin\theta\ud\theta\,P_{l}^{-m}\left(\cos\theta\right)P_{l}^{0}\left(\cos\theta\right)e^{-u\left(\sin\theta\right)^{2}}\left(-v\sin\theta\right)^{2k-m}\\ + & =\int_{0}^{\pi}\sin\theta\ud\theta\,P_{l}^{-m}\left(\cos\theta\right)P_{l}^{0}\left(\cos\theta\right)\left(-v\sin\theta\right)^{2k-m}\sum_{a=0}^{\infty}\frac{\left(-u\right)^{a}}{a!}\left(\sin\theta\right)^{2a}\\ + & =\left(-v\right)^{2k-m}\sum_{a=0}^{\infty}\frac{\left(-u\right)^{a}}{a!}\int_{0}^{\pi}\sin\theta\ud\theta\,P_{l}^{-m}\left(\cos\theta\right)P_{l}^{0}\left(\cos\theta\right)\left(\sin\theta\right)^{2a+2k-m} +\end{align*} + +\end_inset + \end_layout