QPMS
/
qpms
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

Fix the alternative order modematrix implementation 

Former-commit-id: 7daaf636ca803565e2a72b03e1ce184ece24e4c0
This commit is contained in:
Marek Nečada 2019-03-17 23:32:51 +00:00
parent 7db7123c91
commit 86b87c955c
2 changed files with 46 additions and 2 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
qpms/scatsystem.c
View File

@ -1472,7 +1472,7 @@ complex double *qpms_scatsys_build_modeproblem_matrix_irrep_packed_orbitorderR(
const complex double one = 1, zero = 0; const complex double one = 1, zero = 0;
for(qpms_ss_pi_t opistartR = 0; opistartR < ss->p_count; for(qpms_ss_pi_t opistartR = 0; opistartR < ss->p_count;
opistartR += ss->orbit_types[ss->p_orbitinfo[opistartR].t].size //orbit_p_countR; might write a while() instead opistartR += ss->orbit_types[ss->p_orbitinfo[ss->p_by_orbit[opistartR]].t].size //orbit_p_countR; might write a while() instead
) { ) {
const qpms_ss_pi_t orbitstartpiR = ss->p_by_orbit[opistartR]; const qpms_ss_pi_t orbitstartpiR = ss->p_by_orbit[opistartR];
const qpms_ss_oti_t otiR = ss->p_orbitinfo[orbitstartpiR].t; const qpms_ss_oti_t otiR = ss->p_orbitinfo[orbitstartpiR].t;
@ -1495,7 +1495,7 @@ complex double *qpms_scatsys_build_modeproblem_matrix_irrep_packed_orbitorderR(
for(qpms_ss_orbit_pi_t opiR = 0; opiR < orbit_p_countR; ++opiR) { for(qpms_ss_orbit_pi_t opiR = 0; opiR < orbit_p_countR; ++opiR) {
qpms_ss_pi_t piR = ss->p_by_orbit[opistartR + opiR]; qpms_ss_pi_t piR = ss->p_by_orbit[opistartR + opiR];
assert(opiR == ss->p_orbitinfo[piR].p); assert(opiR == ss->p_orbitinfo[piR].p);
const qpms_ss_oti_t otiR = ss->p_orbitinfo[piR].t; assert(otiR == ss->p_orbitinfo[piR].t);
assert(ss->p_orbitinfo[piR].osn == osnR); assert(ss->p_orbitinfo[piR].osn == osnR);
const cart3_t posR = ss->p[piR].pos; const cart3_t posR = ss->p[piR].pos;
// dest particle T-matrix // dest particle T-matrix

44
tests/modeproblem_implementations.py Normal file
View File

@ -0,0 +1,44 @@
from qpms import Particle, CTMatrix, BaseSpec, FinitePointGroup, ScatteringSystem, TMatrixInterpolator, eV, hbar, c
from qpms.symmetries import point_group_info
import numpy as np
import os
import sys
nm = 1e-9
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'
tmfile = '/home/mmn/repo/tmatrix-experiments/Cylinder/AaroBEC/cylinder_50nm_lMax4_cleaned.TMatrix'
#outputdatadir = '/home/necadam1/wrkdir/AaroBECfinite_new'
outputdatadir = '/tmp/u/46/necadam1/unix/project/AaroBECfinite_new'
os.makedirs(outputdatadir, exist_ok = True)
interp = TMatrixInterpolator(tmfile, bspec, symmetrise = sym, atol = 1e-8)
# There is only one t-matrix in the system for each frequency. We initialize the matrix with the lowest frequency data.
# Later, we can replace it using the tmatrix[...] = interp(freq) and s.update_tmatrices NOT YET; TODO
omega = 1.475 * eV/hbar
sv_threshold = 0.5
# Now place the particles and set background index.
px = 571*nm; py = 621*nm
n = 1.51
Nx = 5
Ny = 7
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)
k = n * omega / c
for iri in range(ss.nirreps):
mm_iri_orig = ss.modeproblem_matrix_packed(k, iri)
mm_iri_alt = ss.modeproblem_matrix_packed(k, iri, version='R')
print(np.amax(abs(mm_iri_orig-mm_iri_alt)))