From 3549ceadcca0b62d2b1219ad65c7b75b6d70e412 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Thu, 13 Jun 2019 09:59:32 +0300 Subject: [PATCH] Finite systems tutorial Former-commit-id: 46d65d333b2f0c4dc85359f49d4253cefb7df64b --- finite_systems.md | 105 +++++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 104 insertions(+), 1 deletion(-) diff --git a/finite_systems.md b/finite_systems.md index c623c51..bbe8135 100644 --- a/finite_systems.md +++ b/finite_systems.md @@ -9,7 +9,8 @@ which holds information about particle positions and their T-matrices keeps track about the symmetry group and how the particles transform under the symmetry operations. - +Let's have look how thinks are done on a small python script. +The following script is located in `misc/201903_finiterectlat_AaroBEC.py`. ``` #!/usr/bin/env python @@ -64,5 +65,107 @@ for iri in range(ss.nirreps): # Don't forget to conjugate Vh before transforming it to the full vector! ``` +Let's have a look at the imports. + +``` +from qpms import Particle, CTMatrix, BaseSpec, FinitePointGroup, ScatteringSystem, TMatrixInterpolator, eV, hbar, c +from qpms.symmetries import point_group_info +``` + + * `Particle` is a wrapper over the C structure `qpms_particle_t`, + containing information about particle position and T-matrix. + * `CTMatrix` is a wrapper over the C structure `qpms_tmatrix_t`, + containing a T-matrix. + * `BaseSpec` is a wrapper over the C structure `qpms_vswf_set_spec_t`, + defining with which subset of VSWFs we are working with and how their + respective coefficients are ordered in memory. Typically, this + just means having all electric and magnetic VSWFs up to a given multipole + order `lMax` in the "standard" ordering, but other ways are possible. + Note that different `Particle`s (or, more specifically, `CTMatrix`es) + can have different `BaseSpec`s and happily coexist in the same + `ScatteringSystem`. + This makes sense if the system contains particles with different sizes, + where the larger particles need cutoff at higher multipole orders. + * `FinitePointGroup` is a wrapper over the C structure `qpms_finite_group_t` + containing info about a 3D point group and its representation. Its contents + are currently *not* generated using C code. Rather, it is populated using a + `SVWFPointGroupInfo` instance from the + `point_group_info` python dictionary, which uses sympy to generate the group + and its representation from generators and some metadata. + * `ScatteringSystem` is a wrapper over the C structure `qpms_scatsys_t`, + mentioned earlier, containing info about the whole structure. + * `TMatrixInterpolator` in a wrapper over the C structure + `qpms_tmatrix_interpolator_t` which contains tabulated T-matrices + (calculated e.g. using `scuff-tmatrix`) + and generates frequency-interpolated T-matrices base on these. + * `eV`, `hbar`, `c` are numerical constants with rather obvious meanings. + +Let's go on: +``` +sym = FinitePointGroup(point_group_info['D2h']) +bspec = BaseSpec(lMax = 2) +tmfile = '/m/phys/project/qd/Marek/tmatrix-experiments/Cylinder/AaroBEC/cylinder_50nm_lMax4_cleaned.TMatrix' +... +interp = TMatrixInterpolator(tmfile, bspec, symmetrise = sym, atol = 1e-8) +``` + +The `D2h` choice indicates that our system will have mirror symmetries along +the *xy*, *xz* and *yz* axes. Using the `BaseSpec` with the standard +constructor with `lMax = 2` we declare that we include all the VSWFs up to +quadrupole order. Next, we create a `TMatrixInterpolator` based on a file +created by `scuff-tmatrix`. We force the symmetrisation of the T-matrices with +the same point group as the overall system symmetry in order to eliminate the +possible asymmetries caused by the used mesh. The `atol` parameter just says +that if the absolute value of a given T-matrix element is smaller than the +`atol` value, it is set to zero. + +``` +omega = float(sys.argv[3]) * eV/hbar +sv_threshold = float(sys.argv[4]) + +# Now place the particles and set background index. +px = 571*nm; py = 621*nm +n = 1.51 +Nx = int(sys.argv[1]) +Ny = int(sys.argv[2]) + +orig_x = (np.arange(Nx/2) + (0 if (Nx % 2) else .5)) * px +orig_y = (np.arange(Ny/2) + (0 if (Ny % 2) else .5)) * py + +orig_xy = np.stack(np.meshgrid(orig_x, orig_y), axis = -1) + +tmatrix = interp(omega) +particles = [Particle(orig_xy[i], tmatrix) for i in np.ndindex(orig_xy.shape[:-1])] + +ss = ScatteringSystem(particles, sym) +``` +This chunk sets the light frequency and array size based on a command +line argument. Then it generates a list of particles covering a quarter +of a rectangular array. Finally, these particles are used to generate +the final scattering system – the rest of the particles is generated +automatically to satisfy the specified system symmetry. + +``` +for iri in range(ss.nirreps): + mm_iri = ss.modeproblem_matrix_packed(k, iri) + U, S, Vh = np.linalg.svd(mm_iri) + print(iri, ss.irrep_names[iri], S[-1]) + starti = max(0,len(S) - np.searchsorted(S[::-1], sv_threshold, side='left')-1) + np.savez(os.path.join(outputdatadir, 'Nx%d_Ny%d_%geV_ir%d.npz'%(Nx, Ny, omega/eV*hbar, iri)), + S=S[starti:], omega=omega, Vh = Vh[starti:], iri=iri, Nx = Nx, Ny= Ny ) +``` +The last part iterates over the irreducible representations of the systems. +It generates scattering problem LHS (TODO ref) matrix reduced (projected) +onto each irrep, and performs SVD on that reduced matrix, +and saving the lowest singular values (or all singular values smaller than +`sv_threshold`) together with their respective singular vectors to files. + +The singular vectors corresponding to zero singular values represent the +"modes" of the finite array. + +*TODO analyzing the resulting files.* + + + [scatsystem.h]: @ref scatsystem.h [qpms_scatsys_t]: @ref qpms_scatsys_t