From bb15a2b03502cfe9afd1c284d2adecbd42e3f1ff Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Thu, 23 Jun 2022 18:35:33 +0300 Subject: [PATCH] Handle "empty" irreps in finiterectlat-constant-driving.py --- misc/finiterectlat-constant-driving.py | 37 +++++++++++++++++++------- 1 file changed, 27 insertions(+), 10 deletions(-) diff --git a/misc/finiterectlat-constant-driving.py b/misc/finiterectlat-constant-driving.py index 06c70be..795c6b5 100755 --- a/misc/finiterectlat-constant-driving.py +++ b/misc/finiterectlat-constant-driving.py @@ -71,13 +71,17 @@ def realdipfieldlabels(yp): if yp == 1: return 'y' if yp == 2: return 'z' raise ValueError -def realdipfields(vecgrid, yp): +def realdipfields(vecgrid, yp, bspec=BaseSpec(lMax=1)): + if yp == 0 or yp == 1: + im = np.where(bspec.ilist == 6)[0][0] + ip = np.where(bspec.ilist == 14)[0][0] if yp == 1: - return vecgrid[...,0] + vecgrid[...,2] + return vecgrid[...,ip] + vecgrid[...,ip] if yp == 0: - return -1j*(vecgrid[...,0] - vecgrid[...,2]) + return -1j*(vecgrid[...,im] - vecgrid[...,ip]) if yp == 2: - return vecgrid[...,1] + i0 = np.where(bspec.ilist == 10)[0][0] + return vecgrid[...,i0] raise ValueError def float_nicestr(x, tol=1e-5): @@ -199,6 +203,9 @@ scattered_ir = [None for iri in range(ss.nirreps)] ir_contained = np.ones((nsp, nelem, ss.nirreps), dtype=bool) for iri in range(ss.nirreps): + if ss.saecv_sizes[iri] == 0: + logging.info('irrep %d/%d has an empty VSWF set, skipping') + continue logging.info("processing irrep %d/%d" % (iri, ss.nirreps)) LU = None # to trigger garbage collection before the next call translation_matrix = None @@ -318,15 +325,25 @@ if a.plot or (a.plot_out is not None): vecgrid_ff = np.fft.fftshift(np.fft.fft2(vecgrid, axes=(0,1)),axes=(0,1)) lemax = np.amax(abs(vecgrid)) for yp in range(0,3): - if(np.amax(abs(realdipfields(vecgrid,yp))) > lemax*1e-5): - axes[y,yp*4].imshow(abs(realdipfields(vecgrid,yp)), vmin=0, interpolation='none') - axes[y,yp*4].text(0.5, 0.5, '%g' % np.amax(abs(realdipfields(vecgrid,yp))), horizontalalignment='center', verticalalignment='center', transform=axes[y,yp*4].transAxes) - axes[y,yp*4+1].imshow(np.angle(realdipfields(vecgrid,yp)), vmin=-np.pi, vmax=np.pi, cmap=phasecm, interpolation='none') - axes[y,yp*4+2].imshow(abs(realdipfields(vecgrid_ff,yp)), vmin=0, interpolation='none') - axes[y,yp*4+3].imshow(np.angle(realdipfields(vecgrid_ff,yp)), vmin=-np.pi, vmax=np.pi, cmap=phasecm, interpolation='none') + try: + if(np.amax(abs(realdipfields(vecgrid,yp,bspec))) > lemax*1e-5): + axes[y,yp*4].imshow(abs(realdipfields(vecgrid,yp,bspec)), vmin=0, interpolation='none') + axes[y,yp*4].text(0.5, 0.5, '%g' % np.amax(abs(realdipfields(vecgrid,yp,bspec))), horizontalalignment='center', verticalalignment='center', transform=axes[y,yp*4].transAxes) + axes[y,yp*4+1].imshow(np.angle(realdipfields(vecgrid,yp,bspec)), vmin=-np.pi, vmax=np.pi, cmap=phasecm, interpolation='none') + axes[y,yp*4+2].imshow(abs(realdipfields(vecgrid_ff,yp,bspec)), vmin=0, interpolation='none') + axes[y,yp*4+3].imshow(np.angle(realdipfields(vecgrid_ff,yp,bspec)), vmin=-np.pi, vmax=np.pi, cmap=phasecm, interpolation='none') else: for c in range(0,4): axes[y,yp*4+c].tick_params(bottom=False, left=False, labelbottom=False, labelleft=False) + except (KeyError, IndexError) as e: + logging.info("skipping dipole plot #%d for y=%d, spi=%d (dipole probably not included in the VSWF set) %s" % + (yp, y, spi, e)) + for c in range(0,4): + axx= axes[y,yp*4+c] + axx.tick_params(bottom=False, left=False, labelbottom=False, labelleft=False) + axx.text(0.5, 0.5, 'skipped', horizontalalignment='center',verticalalignment='center', transform=axx.transAxes) + + if not math.isnan(a.ccd_distance): fxye=(-ccd_size/2, ccd_size/2, -ccd_size/2, ccd_size/2) e2vmax = np.amax(np.linalg.norm(ccd_fields[spi,y], axis=-1)**2)