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

Handle "empty" irreps in finiterectlat-constant-driving.py

This commit is contained in:
Marek Nečada 2022-06-23 18:35:33 +03:00
parent 6b89feed70
commit bb15a2b035
1 changed files with 27 additions and 10 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

37
misc/finiterectlat-constant-driving.py
View File

@ -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)