3 from qpms
.argproc
import ArgParser
, annotate_pdf_metadata
8 ap
= ArgParser(['rectlattice2d_finite', 'single_particle', 'single_lMax', 'omega_seq_real_ng', 'planewave'])
9 ap
.add_argument("-o", "--output", type=str, required
=False, help='output path (if not provided, will be generated automatically)')
10 ap
.add_argument("-O", "--plot-out", type=str, required
=False, help="path to plot output (optional)")
11 ap
.add_argument("-P", "--plot", action
='store_true', help="if -p not given, plot to a default path")
12 ap
.add_argument("-g", "--save-gradually", action
='store_true', help="saves the partial result after computing each irrep")
18 logging
.basicConfig(format
='%(asctime)s %(message)s', level
=logging
.INFO
)
22 from qpms
.qpms_p
import cart2sph
, sph2cart
, sph_loccart2cart
, sph_loccart_basis
23 from qpms
.cybspec
import BaseSpec
24 from qpms
.cytmatrices
import CTMatrix
, TMatrixGenerator
25 from qpms
.qpms_c
import Particle
26 from qpms
.cymaterials
import EpsMu
, EpsMuGenerator
, LorentzDrudeModel
, lorentz_drude
27 from qpms
.cycommon
import DebugFlags
, dbgmsg_enable
28 from qpms
import FinitePointGroup
, ScatteringSystem
, BesselType
, eV
, hbar
29 from qpms
.symmetries
import point_group_info
32 dbgmsg_enable(DebugFlags
.INTEGRATION
)
37 particlestr
= ("sph" if a
.height
is None else "cyl") + ("_r%gnm" % (a
.radius
*1e9
))
38 if a
.height
is not None: particlestr
+= "_h%gnm" % (a
.height
* 1e9
)
39 defaultprefix
= "%s_p%gnmx%gnm_%dx%d_m%s_bg%s_φ%gπ_θ(%g_%g)π_ψ%gπ_χ%gπ_f%s_L%d" % (
40 particlestr
, px
*1e9
, py
*1e9
, Nx
, Ny
, str(a
.material
), str(a
.background
), a
.phi
/pi
, np
.amin(a
.theta
)/pi
, np
.amax(a
.theta
)/pi
, a
.psi
/pi
, a
.chi
/pi
, ap
.omega_descr
, a
.lMax
, )
41 logging
.info("Default file prefix: %s" % defaultprefix
)
45 orig_x
= (np
.arange(Nx
/2) + (0 if (Nx
% 2) else .5)) * px
46 orig_y
= (np
.arange(Ny
/2) + (0 if (Ny
% 2) else .5)) * py
48 orig_xy
= np
.stack(np
.meshgrid(orig_x
, orig_y
), axis
= -1)
50 bspec
= BaseSpec(lMax
= a
.lMax
)
51 particles
= [Particle(orig_xy
[i
], ap
.tmgen
, bspec
=bspec
) for i
in np
.ndindex(orig_xy
.shape
[:-1])]
53 sym
= FinitePointGroup(point_group_info
['D2h'])
54 ss
, ssw
= ScatteringSystem
.create(particles
, ap
.background_emg
, ap
.allomegas
[0], sym
=sym
)
58 a
.theta
= np
.atleast_1d(np
.array(a
.theta
))
59 dir_sph_list
= np
.stack((np
.broadcast_to(1, a
.theta
.shape
), a
.theta
, np
.broadcast_to(a
.phi
, a
.theta
.shape
)), axis
=-1)
60 sψ
, cψ
= math
.sin(a
.psi
), math
.cos(a
.psi
)
61 sχ
, cχ
= math
.sin(a
.chi
), math
.cos(a
.chi
)
62 E_sph
= (0., cψ
*cχ
+ 1j
*sψ
*sχ
, sψ
*cχ
+ 1j
*cψ
*sχ
)
64 dir_cart_list
= sph2cart(dir_sph_list
)
65 E_cart_list
= sph_loccart2cart(E_sph
, dir_sph_list
)
67 nfreq
= len(ap
.allomegas
)
68 ndir
= a
.theta
.shape
[0]
70 k_cart_arr
= np
.empty((nfreq
, ndir
, 3), dtype
=float)
71 wavenumbers
= np
.empty((nfreq
,), dtype
=float)
73 σ_ext_arr_ir
= np
.empty((nfreq
, ndir
, ss
.nirreps
), dtype
=float)
74 σ_scat_arr_ir
= np
.empty((nfreq
, ndir
, ss
.nirreps
), dtype
=float)
76 outfile_tmp
= defaultprefix
+ ".tmp" if a
.output
is None else a
.output
+ ".tmp"
78 for i
, omega
in enumerate(ap
.allomegas
):
79 logging
.info("Processing frequency %g eV" % (omega
/ eV
,))
82 if ssw
.wavenumber
.imag
!= 0:
83 warnings
.warn("The background medium wavenumber has non-zero imaginary part. Don't expect emaningful results for cross sections.")
84 wavenumber
= ssw
.wavenumber
.real
85 wavenumbers
[i
] = wavenumber
87 k_sph_list
= np
.array(dir_sph_list
, copy
=True)
88 k_sph_list
[:,0] = wavenumber
90 k_cart_arr
[i
] = sph2cart(k_sph_list
)
92 for iri
in range(ss
.nirreps
):
93 logging
.info("processing irrep %d/%d" % (iri
, ss
.nirreps
))
94 LU
= None # to trigger garbage collection before the next call
95 translation_matrix
= None
96 LU
= ssw
.scatter_solver(iri
)
97 logging
.info("LU solver created")
98 translation_matrix
= ssw
.translation_matrix_packed(iri
, BesselType
.REGULAR
) + np
.eye(ss
.saecv_sizes
[iri
])
99 logging
.info("auxillary translation matrix created")
101 for j
in range(ndir
):
102 k_cart
= k_cart_arr
[i
,j
]
103 # the following two could be calculated only once, but probably not a big deal
104 ã
= ss
.planewave_full(k_cart
=k_cart_arr
[i
,j
], E_cart
=E_cart_list
[j
])
105 Tã
= ssw
.apply_Tmatrices_full(ã
)
107 Tãi
= ss
.pack_vector(Tã
, iri
)
108 ãi
= ss
.pack_vector(ã
, iri
)
110 σ_ext_arr_ir
[i
, j
, iri
] = -np
.vdot(ãi
, fi
).real
/wavenumber
**2
111 σ_scat_arr_ir
[i
, j
, iri
] = np
.vdot(fi
,np
.dot(translation_matrix
, fi
)).real
/wavenumber
**2
113 iriout
= outfile_tmp
+ ".%d.%d" % (i
, iri
)
114 np
.savez(iriout
, omegai
=i
, iri
=iri
, meta
={**vars(a
), 'qpms_version' : qpms
.__version
__()}, omega
=omega
, k_sph
=k_sph_list
, k_cart
= k_cart_arr
, E_cart
=E_cart_list
, E_sph
=np
.array(E_sph
),
115 wavenumber
=wavenumber
, σ_ext_list_ir
=σ_ext_arr_ir
[i
,:,iri
], σ_scat_list_ir
=σ_scat_list_ir
[i
,:,iri
])
116 logging
.info("partial results saved to %s"%iriout
)
118 σ_abs_arr_ir
= σ_ext_arr_ir
- σ_scat_arr_ir
119 σ_abs_arr
= np
.sum(σ_abs_arr_ir
, axis
=-1)
120 σ_scat_arr
= np
.sum(σ_scat_arr_ir
, axis
=-1)
121 σ_ext_arr
= np
.sum(σ_ext_arr_ir
, axis
=-1)
124 outfile
= defaultprefix
+ ".npz" if a
.output
is None else a
.output
125 np
.savez(outfile
, meta
={**vars(a
), 'qpms_version' : qpms
.__version
__()},
126 k_sph
=k_sph_list
, k_cart
= k_cart_arr
, E_cart
=E_cart_list
, E_sph
=np
.array(E_sph
),
127 σ_ext
=σ_ext_arr
,σ_abs
=σ_abs_arr
,σ_scat
=σ_scat_arr
,
128 σ_ext_ir
=σ_ext_arr_ir
,σ_abs_ir
=σ_abs_arr_ir
,σ_scat_ir
=σ_scat_arr_ir
, omega
=ap
.allomegas
, wavenumbers
=wavenumbers
130 logging
.info("Saved to %s" % outfile
)
132 if a
.plot
or (a
.plot_out
is not None):
134 from matplotlib
.backends
.backend_pdf
import PdfPages
135 matplotlib
.use('pdf')
136 from matplotlib
import pyplot
as plt
137 from scipy
.interpolate
import griddata
139 plotfile
= defaultprefix
+ ".pdf" if a
.plot_out
is None else a
.plot_out
140 with
PdfPages(plotfile
) as pdf
:
142 sintheta
= np
.sin(a
.theta
)
143 if False: #len(ap.omega_ranges) != 0:
144 # angle plot ---------------------------------
145 fig
= plt
.figure(figsize
=(210/25.4, 297/25.4))
146 vmax
= max(np
.amax(σ_ext_arr
), np
.amax(σ_scat_arr
), np
.amax(σ_abs_arr
))
147 vmin
= min(np
.amin(σ_ext_arr
), np
.amin(σ_scat_arr
), np
.amin(σ_abs_arr
))
149 ax
= fig
.add_subplot(311)
150 ax
.pcolormesh(a
.theta
, ap
.allomegas
/eh
, σ_ext_arr
, vmin
=vmin
, vmax
=vmax
)
151 ax
.set_xlabel('$\\theta$')
152 ax
.set_ylabel('$\\hbar\\omega / \\mathrm{eV}$')
153 ax
.set_title('$\\sigma_\\mathrm{ext}$')
155 ax
= fig
.add_subplot(312)
156 ax
.pcolormesh(a
.theta
, ap
.allomegas
/eh
, σ_scat_arr
, vmin
=vmin
, vmax
=vmax
)
157 ax
.set_xlabel('$\\theta$')
158 ax
.set_ylabel('$\\hbar\\omega / \\mathrm{eV}$')
159 ax
.set_title('$\\sigma_\\mathrm{scat}$')
161 ax
= fig
.add_subplot(313)
162 im
= ax
.pcolormesh(a
.theta
, ap
.allomegas
/eh
, σ_abs_arr
, vmin
=vmin
, vmax
=vmax
)
163 ax
.set_xlabel('$\\theta$')
164 ax
.set_ylabel('$\\hbar\\omega / \\mathrm{eV}$')
165 ax
.set_title('$\\sigma_\\mathrm{abs}$')
168 fig
.subplots_adjust(right
=0.8)
169 fig
.colorbar(im
, cax
= fig
.add_axes([0.85, 0.15, 0.05, 0.7]))
174 if len(ap
.omega_ranges
) != 0:
175 # "k-space" plot -----------------------------
176 domega
= np
.amin(np
.diff(ap
.allomegas
))
177 dsintheta
= np
.amin(abs(np
.diff(sintheta
)))
178 dk
= dsintheta
* wavenumbers
[0]
181 grid_y
, grid_x
= np
.mgrid
[ap
.allomegas
[0] : ap
.allomegas
[-1] : domega
, np
.amin(sintheta
) * wavenumbers
[-1] : np
.amax(sintheta
) * wavenumbers
[-1] : dk
]
182 imextent
= (np
.amin(sintheta
) * wavenumbers
[-1] / 1e6
, np
.amax(sintheta
) * wavenumbers
[-1] / 1e6
, ap
.allomegas
[0] / eh
, ap
.allomegas
[-1] / eh
)
184 # source coordinates for griddata
185 ktheta
= sintheta
[None, :] * wavenumbers
[:, None]
186 omegapoints
= np
.broadcast_to(ap
.allomegas
[:, None], ktheta
.shape
)
187 points
= np
.stack( (ktheta
.flatten(), omegapoints
.flatten()), axis
= -1)
189 fig
= plt
.figure(figsize
=(210/25.4, 297/25.4))
190 vmax
= np
.amax(σ_ext_arr
)
192 ax
= fig
.add_subplot(311)
193 grid_z
= griddata(points
, σ_ext_arr
.flatten(), (grid_x
, grid_y
), method
= ipm
)
194 ax
.imshow(grid_z
, extent
= imextent
, origin
= 'lower', vmin
= 0, vmax
= vmax
, aspect
= 'auto', interpolation
='none')
195 ax
.set_xlabel('$k_\\theta / \\mathrm{\\mu m^{-1}}$')
196 ax
.set_ylabel('$\\hbar\\omega / \\mathrm{eV}$')
197 ax
.set_title('$\\sigma_\\mathrm{ext}$')
199 ax
= fig
.add_subplot(312)
200 grid_z
= griddata(points
, σ_scat_arr
.flatten(), (grid_x
, grid_y
), method
= ipm
)
201 ax
.imshow(grid_z
, extent
= imextent
, origin
= 'lower', vmin
= 0, vmax
= vmax
, aspect
= 'auto', interpolation
='none')
202 ax
.set_xlabel('$k_\\theta / \\mathrm{\\mu m^{-1}}$')
203 ax
.set_ylabel('$\\hbar\\omega / \\mathrm{eV}$')
204 ax
.set_title('$\\sigma_\\mathrm{scat}$')
206 ax
= fig
.add_subplot(313)
207 grid_z
= griddata(points
, σ_abs_arr
.flatten(), (grid_x
, grid_y
), method
= ipm
)
208 im
= ax
.imshow(grid_z
, extent
= imextent
, origin
= 'lower', vmin
= 0, vmax
= vmax
, aspect
= 'auto', interpolation
='none')
209 ax
.set_xlabel('$k_\\theta / \\mathrm{\\mu m^{-1}}$')
210 ax
.set_ylabel('$\\hbar\\omega / \\mathrm{eV}$')
211 ax
.set_title('$\\sigma_\\mathrm{abs}$')
213 fig
.subplots_adjust(right
=0.8)
214 fig
.colorbar(im
, cax
= fig
.add_axes([0.85, 0.15, 0.05, 0.7]))
219 for omega
in ap
.omega_singles
:
220 i
= np
.searchsorted(ap
.allomegas
, omega
)
223 fig
.suptitle("%g eV" % (omega
/ eh
))
224 ax
= fig
.add_subplot(111)
225 sintheta
= np
.sin(a
.theta
)
226 ax
.plot(sintheta
, σ_ext_arr
[i
]*1e12
,label
='$\sigma_\mathrm{ext}$')
227 ax
.plot(sintheta
, σ_scat_arr
[i
]*1e12
, label
='$\sigma_\mathrm{scat}$')
228 ax
.plot(sintheta
, σ_abs_arr
[i
]*1e12
, label
='$\sigma_\mathrm{abs}$')
230 ax
.set_xlabel('$\sin\\theta$')
231 ax
.set_ylabel('$\sigma/\mathrm{\mu m^2}$')
235 annotate_pdf_metadata(pdf
, scriptname
="finiterectlat-scatter.py")