9 "import numpy as np\n",
12 "from qpms.cybspec import BaseSpec\n",
13 "from qpms.cytmatrices import CTMatrix, TMatrixGenerator, TMatrixInterpolator\n",
14 "from qpms.qpms_c import Particle, pgsl_ignore_error\n",
15 "from qpms.cymaterials import EpsMu, EpsMuGenerator, LorentzDrudeModel, lorentz_drude\n",
16 "from qpms.cycommon import DebugFlags, dbgmsg_enable\n",
17 "from qpms import FinitePointGroup, ScatteringSystem, BesselType, eV, hbar\n",
18 "import scipy.constants as sci\n",
30 "a1 = np.array([0,period]) \n",
31 "a2 = np.array([period,0])"
40 "#Particle positions\n",
43 "orig_xy = np.stack(np.meshgrid(orig_x,orig_y),axis=-1)"
53 "refractive_index = 1.52 # for background medium\n",
54 "height = 50e-9 # Particle height\n",
55 "radius = 50e-9 # Particle radius\n",
56 "medium = EpsMu(refractive_index**2) # non-lossy background medium with constant refr. index #OK\n",
57 "# global symmetry group of the system\n",
58 "#sym = FinitePointGroup(point_group_info['D4h'])\n",
60 "metal = lorentz_drude['Ag']\n",
61 "kx_lim = np.array([-0.2, 0.2], dtype=float)\n",
73 "(7.59633723939313, 8.305328715069821, 7.94387469344152, 8.001475225494097)"
78 "output_type": "execute_result"
82 "omega_scuff_min=1.5*eh/3e14\n",
83 "omega_scuff_max=1.64*eh/3e14\n",
84 "omega_scuff=omega/3e14\n",
85 "omega_scuff_slr=(2*np.pi*sci.c/(1.52*0.52e-6))/3e14\n",
86 "omega_scuff_min, omega_scuff_max, omega_scuff_slr, omega_scuff"
95 "bspec = BaseSpec(lMax = 1)\n",
97 "tmfile_scuffOld = '/home/javier/tmatrices/sphereAg_50nm_oldScuff.TMatrix'\n",
98 "tmfile_scuffNew = '/home/javier/tmatrices/sphere50nm_newScuff.TMatrix'\n",
99 "interp_old = TMatrixInterpolator(tmfile_scuffOld, bspec, atol=1e-8) \n",
100 "interp_new = TMatrixInterpolator(tmfile_scuffNew, bspec, atol=1e-8) \n",
101 "tmscuff_not_fixed = interp_old(omega)\n",
102 "tmscuff_bugfixed = interp_new(omega)"
107 "execution_count": 7,
111 "tmgen = TMatrixGenerator.sphere(medium, metal, radius)"
116 "execution_count": 8,
120 "tmgen_omega=tmgen(bspec,omega)"
125 "execution_count": 9,
131 "array([[-0.07824951+0.25215549j, 0. +0.j , 0. +0.j , 0. +0.j , 0. +0.j , 0. +0.j ],\n",
132 " [ 0. +0.j , -0.07824951+0.25215549j, 0. +0.j , 0. +0.j , 0. +0.j , 0. +0.j ],\n",
133 " [ 0. +0.j , 0. +0.j , -0.07824951+0.25215549j, 0. +0.j , 0. +0.j , 0. +0.j ],\n",
134 " [ 0. +0.j , 0. +0.j , 0. +0.j , -0.00083788-0.01420874j, 0. +0.j , 0. +0.j ],\n",
135 " [ 0. +0.j , 0. +0.j , 0. +0.j , 0. +0.j , -0.00083788-0.01420874j, 0. +0.j ],\n",
136 " [ 0. +0.j , 0. +0.j , 0. +0.j , 0. +0.j , 0. +0.j , -0.00083788-0.01420874j]])"
139 "execution_count": 9,
141 "output_type": "execute_result"
145 "tmgen_omega.as_ndarray() # T-Matrix from generator"
150 "execution_count": 10,
156 "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",
157 " -7.76488937e-14+1.08532726e-13j],\n",
158 " [-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",
159 " 7.76590838e-08-1.38092126e-08j],\n",
160 " [-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",
161 " 4.59276044e-14-4.47400043e-14j],\n",
162 " [ 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",
163 " 1.16328171e-08+6.36863105e-09j],\n",
164 " [ 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",
165 " 2.38800352e-13+1.35401679e-13j],\n",
166 " [ 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",
167 " -3.03166170e-05-2.01474828e-03j]])"
170 "execution_count": 10,
172 "output_type": "execute_result"
176 "tmscuff_not_fixed.as_ndarray() # T-Matrix of not fixed version of Scuff-EM"
181 "execution_count": 11,
187 "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",
188 " -1.55480479e-12+2.17400571e-12j],\n",
189 " [-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",
190 " 1.43308826e-06-2.58864855e-07j],\n",
191 " [-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",
192 " 9.21093604e-13-8.94225073e-13j],\n",
193 " [ 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",
194 " 1.52945653e-07+8.36194839e-08j],\n",
195 " [-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",
196 " 4.78231690e-12+2.71099900e-12j],\n",
197 " [ 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",
198 " -6.07143303e-04-4.03488631e-02j]])"
201 "execution_count": 11,
203 "output_type": "execute_result"
207 "tmscuff_bugfixed.as_ndarray() # T-Matrix of FIXED version of Scuff-EM"
212 "execution_count": 12,
218 "(matrix([[-9.43033280e-02+3.59361755e-01j, 3.20658236e-12+8.41416904e-12j],\n",
219 " [-8.84386127e-12+1.96169624e-12j, -9.42395846e-02+3.59248336e-01j]]),\n",
220 " matrix([[-9.43033280e-02-3.59361755e-01j, -8.84386127e-12-1.96169624e-12j],\n",
221 " [ 3.20658236e-12-8.41416904e-12j, -9.42395846e-02-3.59248336e-01j]]))"
224 "execution_count": 12,
226 "output_type": "execute_result"
230 "#I play around with the different operations so everything is correct\n",
231 "tmscuffnew = tmscuff_bugfixed.as_ndarray()\n",
232 "tmscuffnew_mat = np.asmatrix(tmscuffnew)\n",
233 "tmscuffnew_dag = tmscuffnew_mat.getH()\n",
234 "tmscuffnew_mat[0:2,0:2], tmscuffnew_dag[0:2,0:2]"
239 "execution_count": 13,
245 "(matrix([[1, 0],\n",
246 " [9, 4]]), matrix([[4, 9],\n",
250 "execution_count": 13,
252 "output_type": "execute_result"
256 "A = np.matrix([[0,1],[2,3]])\n",
257 "B = np.matrix([[3,2],[1,0]])\n",
258 "AB = np.dot(A,B)\n",
259 "BA = np.dot(B,A)\n",
265 "execution_count": 14,
269 "sum1new = np.dot(tmscuffnew_dag,tmscuffnew) #is this the right order? Regarding the above, yes it is.\n",
270 "sum2new = 0.5*tmscuffnew.__add__(tmscuffnew_dag) "
275 "execution_count": 15,
279 "powermatrix_scuffnew = sum1new.__add__(sum2new) #powermatrix for bugfixed scuff"
284 "execution_count": 16,
288 "#Power matrix for NOT bugfixed scuff\n",
289 "tmscuffold = tmscuff_not_fixed.as_ndarray()\n",
290 "tmscuffold_mat = np.asmatrix(tmscuffold)\n",
291 "tmscuffold_dag = tmscuffold_mat.getH()\n",
292 "sum1old = np.dot(tmscuffold_dag,tmscuffold) #is this the right order? Regarding the above, yes it is.\n",
293 "sum2old = 0.5*tmscuffold.__add__(tmscuffold_dag)\n",
295 "powermatrix_scuffold = sum1old.__add__(sum2old) "
300 "execution_count": 17,
304 "#Power matrix for T matrix generator\n",
305 "tmscuffgen = tmgen_omega.as_ndarray()\n",
306 "tmscuffgen_mat = np.asmatrix(tmscuffgen)\n",
307 "tmscuffgen_dag = tmscuffgen_mat.getH()\n",
308 "sum1gen = np.dot(tmscuffgen_dag,tmscuffgen) #is this the right order? Regarding the above, yes it is.\n",
309 "sum2gen = 0.5*tmscuffgen.__add__(tmscuffgen_dag)\n",
311 "powermatrix_scuffgen = sum1gen.__add__(sum2gen) "
316 "execution_count": 18,
322 "((-0.008544129580638216+0j),\n",
323 " (-0.008544129580638216+0j),\n",
324 " (-0.008544129580638216+0j),\n",
325 " (-0.0006352854997266949+0j),\n",
326 " (-0.0006352854997266949+0j),\n",
327 " (-0.0006352854997266949+0j))"
330 "execution_count": 18,
332 "output_type": "execute_result"
336 "#Power matrix of generator is diagonal and all its eigenvalues are negative:\n",
337 "powermatrix_scuffgen[0,0], powermatrix_scuffgen[1,1], powermatrix_scuffgen[2,2], powermatrix_scuffgen[3,3], powermatrix_scuffgen[4,4], powermatrix_scuffgen[5,5]"
342 "execution_count": 19,
348 "((0.04373066071852283+0j),\n",
349 " (0.04370088172427582+0j),\n",
350 " (0.043730660553373296+0j),\n",
351 " (0.001021256123757952+0j),\n",
352 " (0.0010213406041430792+0j),\n",
353 " (0.0010212560750220145+0j))"
356 "execution_count": 19,
358 "output_type": "execute_result"
362 "powermatrix_scuffnew[0,0], powermatrix_scuffnew[1,1], powermatrix_scuffnew[2,2], powermatrix_scuffnew[3,3], powermatrix_scuffnew[4,4], powermatrix_scuffnew[5,5]"
367 "execution_count": 20,
373 "((-0.004364703222422517+0j),\n",
374 " (-0.004361872918812184+0j),\n",
375 " (-0.004364703232328262+0j),\n",
376 " (-2.6256479824067564e-05+0j),\n",
377 " (-2.6262303053141142e-05+0j),\n",
378 " (-2.6256487319428013e-05+0j))"
381 "execution_count": 20,
383 "output_type": "execute_result"
387 "powermatrix_scuffold[0,0], powermatrix_scuffold[1,1], powermatrix_scuffold[2,2], powermatrix_scuffold[3,3], powermatrix_scuffold[4,4], powermatrix_scuffold[5,5]"
392 "execution_count": 22,
398 "((-10.019160178833705-0j),\n",
399 " (-10.018834234211564-0j),\n",
400 " (-10.019160118257586-0j),\n",
401 " (-38.895393845668316-0j),\n",
402 " (-38.88998623145963-0j),\n",
403 " (-38.89538088616882-0j))"
406 "execution_count": 22,
408 "output_type": "execute_result"
412 "#So! There is a fixed factor between the N/M elements of Scuff new and old:\n",
413 "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]"
418 "execution_count": 23,
424 "((-5.118211317583541-0j),\n",
425 " (-5.11472600126595-0j),\n",
426 " (-5.118211298254535-0j),\n",
427 " (-1.6075545942687262-0j),\n",
428 " (-1.6076875744566317-0j),\n",
429 " (-1.60755451755371-0j))"
432 "execution_count": 23,
434 "output_type": "execute_result"
438 "#And there is a fixed factor for new / generated as well! (hence also for old / generated)\n",
439 "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]"
444 "execution_count": 21,
450 "((-1.60755451755371-0j), (-2.557363000632975+0j))"
453 "execution_count": 21,
455 "output_type": "execute_result"
459 "powermatrix_scuffnew[5,5]/powermatrix_scuffgen[5,5], powermatrix_scuffnew[1,1]/powermatrix_scuffgen[1,1]/2"
464 "execution_count": 22,
470 "((1.599227713011597e-16-0j),\n",
471 " (8.902299490051238e-14+0j),\n",
472 " (1.504486797691096e-21-0j))"
475 "execution_count": 22,
477 "output_type": "execute_result"
481 "#Let's also calculate the determinants\n",
482 "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",
483 "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",
484 "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",
485 "detgen, detnew, detold"
490 "execution_count": 23,
496 "((1.5992277130116025e-16+0j),\n",
497 " (8.902299447416401e-14+2.6128569373609595e-39j),\n",
498 " (1.504486701753944e-21-2.892571608172536e-45j))"
501 "execution_count": 23,
503 "output_type": "execute_result"
507 "np.linalg.det(powermatrix_scuffgen), np.linalg.det(powermatrix_scuffnew), np.linalg.det(powermatrix_scuffold)"
512 "execution_count": 24,
516 "#We try to normalize the power matrix elements for each case with the corresponding determinants."
521 "execution_count": 25,
527 "((-53426597795433.88+0j),\n",
528 " (-53426597795433.88+0j),\n",
529 " (-53426597795433.88+0j),\n",
530 " (-3972451793812.1055+0j),\n",
531 " (-3972451793812.1055+0j),\n",
532 " (-3972451793812.1055+0j))"
535 "execution_count": 25,
537 "output_type": "execute_result"
541 "#T-matrix generator: \n",
542 "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 "
547 "execution_count": 26,
553 "((-2.901124309712079e+18+0j),\n",
554 " (-2.899243067806416e+18+0j),\n",
555 " (-2.901124316296214e+18+0j),\n",
556 " (-1.7452117136795634e+16+0j),\n",
557 " (-1.7455987711853198e+16+0j),\n",
558 " (-1.7452122118800434e+16+0j))"
561 "execution_count": 26,
563 "output_type": "execute_result"
567 "#Not bugfixed Scuff-EM: \n",
568 "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 "
573 "execution_count": 27,
579 "((491228819782.9563+0j),\n",
580 " (490894310768.96173+0j),\n",
581 " (491228817927.8229+0j),\n",
582 " (11471823936.043226+0j),\n",
583 " (11472772908.667904+0j),\n",
584 " (11471823388.58987+0j))"
587 "execution_count": 27,
589 "output_type": "execute_result"
594 "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 "
599 "execution_count": 28,
603 "#It might make more sense to renormalize the electric and magnetic parts of the power matrices separately:"
608 "execution_count": 29,
614 "((13698.226641508763-0j),\n",
615 " (13698.226641508763-0j),\n",
616 " (13698.226641508763-0j),\n",
617 " (2477776.413504498-0j),\n",
618 " (2477776.413504498-0j),\n",
619 " (2477776.413504498-0j))"
622 "execution_count": 29,
624 "output_type": "execute_result"
628 "#T-matrix generator:\n",
629 "detgen_el = powermatrix_scuffgen[0,0]*powermatrix_scuffgen[1,1]*powermatrix_scuffgen[2,2]\n",
630 "detgen_mag = powermatrix_scuffgen[3,3]*powermatrix_scuffgen[4,4]*powermatrix_scuffgen[5,5]\n",
631 "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 "
636 "execution_count": 30,
642 "((52525.751093170744-0j),\n",
643 " (52491.69062780147-0j),\n",
644 " (52525.75121237855-0j),\n",
645 " (1450208903.6280527-0j),\n",
646 " (1450530534.6580672-0j),\n",
647 " (1450209317.614944-0j))"
650 "execution_count": 30,
652 "output_type": "execute_result"
656 "#NOT bugfixed Scuff-EM:\n",
657 "detold_el = powermatrix_scuffold[0,0]*powermatrix_scuffold[1,1]*powermatrix_scuffold[2,2]\n",
658 "detold_mag = powermatrix_scuffold[3,3]*powermatrix_scuffold[4,4]*powermatrix_scuffold[5,5]\n",
659 "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 "
664 "execution_count": 31,
670 "((523.2675016470721+0j),\n",
671 " (522.9111754524715+0j),\n",
672 " (523.2674996709441+0j),\n",
673 " (958726.5381294343+0j),\n",
674 " (958805.8459399715+0j),\n",
675 " (958726.4923775065+0j))"
678 "execution_count": 31,
680 "output_type": "execute_result"
684 "#Bugfixed Scuff-EM:\n",
685 "detnew_el = powermatrix_scuffnew[0,0]*powermatrix_scuffnew[1,1]*powermatrix_scuffnew[2,2]\n",
686 "detnew_mag = powermatrix_scuffnew[3,3]*powermatrix_scuffnew[4,4]*powermatrix_scuffnew[5,5]\n",
687 "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 "
692 "execution_count": null,
700 "display_name": "Python 3",
701 "language": "python",
709 "file_extension": ".py",
710 "mimetype": "text/x-python",
712 "nbconvert_exporter": "python",
713 "pygments_lexer": "ipython3",