misc/201903_finiterectlat_AaroBEC_fatMie.py

   1 #!/usr/bin/env python3
   2 # coding: utf-8
   3 from qpms import Particle, CTMatrix, BaseSpec, FinitePointGroup, ScatteringSystem, TMatrixInterpolator, eV, hbar, c, MaterialInterpolator, scatsystem_set_nthreads
   4 from qpms.symmetries import point_group_info
   5 import numpy as np
   6 import os
   7 import sys
   8 from pathlib import Path
   9
  10 if 'SLURM_CPUS_PER_TASK' in os.environ:
  11     scatsystem_set_nthreads(int(os.environ['SLURM_CPUS_PER_TASK']))
  12
  13 nm = 1e-9
  14
  15 rewrite_output = '--rewrite-output' in sys.argv
  16
  17 radiusfactor = float(sys.argv[5])
  18
  19 cyr_part_height = 50*nm
  20 cyr_part_radius = 50*nm
  21 cyr_part_volume = cyr_part_height * np.pi * cyr_part_radius**2
  22 eqv_sph_radius = (3/4/np.pi*cyr_part_volume)**(1/3) * radiusfactor
  23
  24 sym = FinitePointGroup(point_group_info['D2h'])
  25 bspec = BaseSpec(lMax = 2)
  26 #tmfile = '/m/phys/project/qd/Marek/tmatrix-experiments/Cylinder/AaroBEC/cylinder_50nm_lMax4_cleaned.TMatrix'
  27 materialfile = '/home/necadam1/wrkdir/repo/refractiveindex.info-database/database/data/main/Au/Johnson.yml'
  28
  29 #outputdatadir = '/home/necadam1/wrkdir/AaroBECfinite_new'
  30 #outputdatadir = '/u/46/necadam1/unix/project/AaroBECfinite_sph'
  31 outputdatadir = '/home/necadam1/wrkdir/AaroBECfinite_fatsph'
  32 os.makedirs(outputdatadir, exist_ok = True)
  33 mi = MaterialInterpolator(materialfile)
  34 #interp = TMatrixInterpolator(tmfile, bspec, symmetrise = sym, atol = 1e-8)
  35 # There is only one t-matrix in the system for each frequency. We initialize the matrix with the lowest frequency data.
  36 # Later, we can replace it using the tmatrix[...] = interp(freq) and s.update_tmatrices NOT YET; TODO
  37
  38 omega = float(sys.argv[3]) * eV/hbar
  39 sv_threshold = float(sys.argv[4])
  40
  41
  42 # Now place the particles and set background index.
  43 px = 571*nm; py = 621*nm
  44 n = 1.52
  45 Nx = int(sys.argv[1])
  46 Ny = int(sys.argv[2])
  47
  48 orig_x = (np.arange(Nx/2) + (0 if (Nx % 2) else .5)) * px
  49 orig_y = (np.arange(Ny/2) + (0 if (Ny % 2) else .5)) * py
  50
  51 orig_xy = np.stack(np.meshgrid(orig_x, orig_y), axis = -1)
  52
  53 #tmatrix = interp(omega)
  54 tmatrix = CTMatrix.spherical_perm(bspec, eqv_sph_radius, omega, mi(omega), n**2)
  55 particles = [Particle(orig_xy[i], tmatrix) for i in np.ndindex(orig_xy.shape[:-1])]
  56
  57
  58 ss = ScatteringSystem(particles, sym)
  59
  60
  61 k = n * omega / c
  62
  63
  64 for iri in range(ss.nirreps):
  65     destpath = os.path.join(outputdatadir, 'Nx%d_Ny%d_%geV_ir%d_r%gnm.npz'%(Nx, Ny, omega/eV*hbar, iri, eqv_sph_radius/nm))
  66     touchpath = os.path.join(outputdatadir, 'Nx%d_Ny%d_%geV_ir%d_r%gnm.done'%(Nx, Ny, omega/eV*hbar, iri, eqv_sph_radius/nm))
  67     if (os.path.isfile(destpath) or os.path.isfile(touchpath)) and not rewrite_output:
  68       print(destpath, 'already exists, skipping')
  69       continue
  70     mm_iri = ss.modeproblem_matrix_packed(k, iri)
  71     U, S, Vh = np.linalg.svd(mm_iri)
  72     del U
  73     print(iri, ss.irrep_names[iri], S[-1])
  74     starti = max(0,len(S) - np.searchsorted(S[::-1], sv_threshold, side='left')-1)
  75     np.savez(destpath,
  76         S=S[starti:], omega=omega, Vh = Vh[starti:], iri=iri, Nx = Nx, Ny= Ny )
  77     del S
  78     del Vh
  79     Path(touchpath).touch()
  80     # Don't forget to conjugate Vh before transforming it to the full vector!
  81