diff --git a/qpms/hexpoints.py b/qpms/hexpoints.py deleted file mode 100644 index c5871e2..0000000 --- a/qpms/hexpoints.py +++ /dev/null @@ -1,588 +0,0 @@ -import math -import numpy as np -nx = None - -_s3 = math.sqrt(3) - - -def generate_trianglepoints(maxlayer, include_origin = False, v3d = True, circular = True, sixthindices = False, mirrorindices = False): - _e6 = np.array([[math.cos(2*math.pi*i/6),math.sin(2*math.pi*i/6),0] if v3d else [math.cos(2*math.pi*i/6),math.sin(2*math.pi*i/6)] for i in range(6)]) - points = np.empty((3*maxlayer*(maxlayer+1)+(1 if include_origin else 0), 3 if v3d else 2)) - point_i = 0 - if (include_origin): - points[0] = np.array((0,0,0) if v3d else (0,0)) - point_i = 1 - if sixthindices: - si = np.empty((6,(maxlayer*(maxlayer+1))//2), dtype=int) - sii = [0,0,0,0,0,0] - if mirrorindices: - if(maxlayer < 3): - mi = np.empty((2,0,), dtype=int) - else: - #layer indices start from one! - ilayer = np.arange(1, maxlayer+1) # We need first to "copy" layer indices to correspond to the Muster count - mustercount = (ilayer - 1)//2 - mustercum = np.cumsum(mustercount) - layerstart = np.zeros((mustercum[maxlayer - 1]), dtype=int) - layerstart[mustercum[:(maxlayer-1)]] = 1 - layer = np.cumsum(layerstart) + 2 # That's it - lb = 3*layer*(layer-1) # layer base (lowest) index - li = np.arange(len(layer)) - mustercum[layer-2] # muster indices for each layers - mi = np.empty((2, len(layer)), dtype=int) - mi[0] = lb + 1 + li + include_origin - mi[1] = lb + layer - (1 + li) + include_origin - # there are two non-musters in each even layer, one non-muster in each odd - layer = (2*np.arange(((3*maxlayer)//2))+1)//3 + 1 - nmi = 3*layer*(layer-1) - nmi[2::3] += layer[2::3] // 2 # second non-musters in even layers - nmi += include_origin - - for layer in range(1,maxlayer+1): - for i in range(6): - base = _e6[i]*layer - shift = _e6[(i+2)%6] - ar = np.arange(layer) - points[point_i:(point_i+layer)] = base[nx,:] + ar[:,nx] * shift[nx,:] - if sixthindices: - si[i, sii[i]:sii[i]+layer] = point_i + ar - sii[i] += layer - point_i += layer - if (circular): - mask = (np.sum(points * points, axis = -1) <= maxlayer * maxlayer * 3/ 4 + 0.1) # UGLY FIX OF ASYMMETRY BECAUSE OF ROUNDING ERROR - points = points[mask] - if sixthindices: - cum = np.cumsum(mask) - 1 - mask0 = mask[si[0]] - si_ = si[:,mask0] - si = cum[si_] - if mirrorindices: - cum = np.cumsum(mask) - 1 - mask0 = mask[mi[0]] - mi_ = mi[:,mask0] - mi = cum[mi_] - mask0 = mask[nmi] - nmi_ = nmi[mask0] - nmi = cum[nmi_] - if not (mirrorindices or sixthindices): - return points - else: - return {'points': points, - 'si' : si if sixthindices else None, - 'mi' : mi if mirrorindices else None, - 'nmi' : nmi if mirrorindices else None} - -def generate_trianglepoints_hexcomplement(maxlayer, v3d = True, circular = True, thirdindices = False, mirrorindices=False): - _e6 = np.array([[math.cos(2*math.pi*i/6),math.sin(2*math.pi*i/6),0] if v3d else [math.cos(2*math.pi*i/6),math.sin(2*math.pi*i/6)] for i in range(6)]) - _f6 = np.array([[-math.sin(2*math.pi*i/6),math.cos(2*math.pi*i/6),0] if v3d else [math.sin(2*math.pi*i/6),-math.cos(2*math.pi*i/6)] for i in range(6)]) - points = np.empty((3*maxlayer*maxlayer, 3 if v3d else 2)) - point_i = 0 - # 3 * layer ** 2 is the basis index for a layer, a layer contains 3 * (2*layer + 1) points - if thirdindices: - ii = np.arange(maxlayer**2) - layer = np.empty((maxlayer**2), dtype=int) - layer = np.sqrt(ii, out=layer, casting='unsafe') - #ti0 = 2*layer**2 + ii - ti = np.arange(3)[:, nx] * (2*layer + 1)[nx, :] + (2*layer**2 + ii)[nx,:] - if mirrorindices: - ii = np.arange(((maxlayer-1)*maxlayer)//2) - layer = np.empty((((maxlayer-1)*maxlayer)//2), dtype=int) - layer = (np.sqrt(1+8*ii, out=layer, casting='unsafe')+1)//2 - li = ii - ((layer ) * (layer-1))//2# numbers indices in a each layer - lb = 3*layer **2 # base index of a layer - mi = np.empty((2,len(layer)), dtype=int) - mi[0] = lb + li + layer % 2 - mi[1] = lb + 2*layer - li - # indices of non-mirrored/self-mirrored - layer = np.arange(maxlayer) - lb = 3 * layer**2 - nmi = lb + ((layer + 1) % 2) * layer - for layer in range(0,maxlayer): - if (layer % 2): # odd layer - for i in range(3): - base = _f6[(2*i-1)%6] * ((0.5 + 1.5 * layer) / _s3) - shift = _e6[(2*i+2)%6] - count = (layer + 1) // 2 - ar = np.arange(count) - points[point_i:point_i+count,:] = base + ar[:,nx]*shift[nx,:] - point_i += count - - base = _e6[(2*i+1)%6]*layer + _f6[(2*i)%6] / _s3 - shift = _e6[(2*i+3)%6] - count = layer - ar = np.arange(count) - points[point_i:point_i+count,:] = base + ar[:,nx]*shift[nx,:] - point_i += count - - base = _e6[(2*i+2)%6]*layer + _f6[(2*i)%6] / _s3 - shift = _e6[(2*i+4)%6] - count = (layer + 1) // 2 - ar = np.arange(count) - points[point_i:point_i+count,:] = base + ar[:,nx]*shift[nx,:] - point_i += count - else: # even layer: - for i in range(3): - shift = _e6[(2*i+2)%6] - base = _f6[(2*i-1)%6] * ((0.5 + 1.5 * layer) / _s3) + shift / 2 - count = layer // 2 - ar = np.arange(count) - points[point_i:point_i+count,:] = base + ar[:,nx]*shift[nx,:] - point_i += count - - base = _e6[(2*i+1)%6]*layer + _f6[(2*i)%6] / _s3 - shift = _e6[(2*i+3)%6] - count = layer - ar = np.arange(count) - points[point_i:point_i+count,:] = base + ar[:,nx]*shift[nx,:] - point_i += count - - base = _e6[(2*i+2)%6]*layer + _f6[(2*i)%6] / _s3 - shift = _e6[(2*i+4)%6] - count = (layer + 2) // 2 - ar = np.arange(count) - points[point_i:point_i+count,:] = base + ar[:,nx]*shift[nx,:] - point_i += count - #if (mirrorindices): - - if (circular): - mask = (np.sum(points * points, axis = -1) <= maxlayer * maxlayer * 3/ 4 + 0.01) # UGLY FIX OF ASYMMETRY BECAUSE OF ROUNDING ERROR - points = points[mask] - if thirdindices: - cum = np.cumsum(mask) - 1 - mask0 = mask[ti[0]] - ti_ = ti[:,mask0] - ti = cum[ti_] - if mirrorindices: - cum = np.cumsum(mask) - 1 - mask0 = mask[mi[0]] - mi_ = mi[:,mask0] - mi = cum[mi_] - mask0 = mask[nmi] - nmi_ = nmi[mask0] - nmi = cum[nmi_] - if not (mirrorindices or thirdindices): - return points - else: - return {'points': points, - 'ti' : ti if thirdindices else None, - 'mi' : mi if mirrorindices else None, - 'nmi' : nmi if mirrorindices else None - } - - -from .cycommon import get_mn_y -from .cytranslations import trans_calculator -from .qpms_p import cart2sph - -def hexlattice_precalc_AB_save(file, lMax, k_hexside, maxlayer, circular=True, savepointinfo = False, J_scat=3): - params = { - 'lMax' : lMax, - 'k_hexside' : k_hexside, - 'maxlayer' : maxlayer, - 'circular' : circular, - 'savepointinfo' : savepointinfo, - 'J_scat' : J_scat - } - tpdict = generate_trianglepoints(maxlayer, v3d=True, circular=circular, sixthindices=True, mirrorindices=True) - tphcdict = generate_trianglepoints_hexcomplement(maxlayer, v3d=True, circular=circular, thirdindices=True, mirrorindices=True) - my, ny = get_mn_y(lMax) - nelem = len(my) - a_self_nm = np.empty((tpdict['nmi'].shape[0],nelem,nelem), dtype=complex) - b_self_nm = np.empty((tpdict['nmi'].shape[0],nelem,nelem), dtype=complex) - a_self_m0 = np.empty((tpdict['mi'].shape[1],nelem,nelem), dtype=complex) - b_self_m0 = np.empty((tpdict['mi'].shape[1],nelem,nelem), dtype=complex) - a_d2u_nm = np.empty((tphcdict['nmi'].shape[0],nelem,nelem), dtype=complex) - b_d2u_nm = np.empty((tphcdict['nmi'].shape[0],nelem,nelem), dtype=complex) - a_d2u_m0 = np.empty((tphcdict['mi'].shape[1],nelem,nelem), dtype=complex) - b_d2u_m0 = np.empty((tphcdict['mi'].shape[1],nelem,nelem), dtype=complex) - - k_0 = k_hexside*_s3 # not really a wave vector here because of the normalisation! - tc = trans_calculator(lMax) - - y = np.arange(nelem) - - points = tpdict['points'][tpdict['nmi']] - d_i2j = cart2sph(points) - a_self_nm, b_self_nm = tc.get_AB_arrays(k_0*d_i2j[:,0],d_i2j[:,1],d_i2j[:,2],np.array([False]),J_scat) - - points = tpdict['points'][tpdict['mi'][0]] - d_i2j = cart2sph(points) - a_self_m0, b_self_m0 = tc.get_AB_arrays(k_0*d_i2j[:,0],d_i2j[:,1],d_i2j[:,2],np.array([False]),J_scat) - - points = tphcdict['points'][tphcdict['nmi']] - d_i2j = cart2sph(points) - a_d2u_nm, b_d2u_nm = tc.get_AB_arrays(k_0*d_i2j[:,0],d_i2j[:,1],d_i2j[:,2],np.array([False]),J_scat) - - points = tphcdict['points'][tphcdict['mi'][0]] - d_i2j = cart2sph(points) - a_d2u_m0, b_d2u_m0 = tc.get_AB_arrays(k_0*d_i2j[:,0],d_i2j[:,1],d_i2j[:,2],np.array([False]),J_scat) - - tosave = { - 'a_self_nm' : a_self_nm, - 'a_self_m0' : a_self_m0, - 'b_self_nm' : b_self_nm, - 'b_self_m0' : b_self_m0, - 'a_d2u_nm' : a_d2u_nm, - 'a_d2u_m0' : a_d2u_m0, - 'b_d2u_nm' : b_d2u_nm, - 'b_d2u_m0' : b_d2u_m0, - 'precalc_params' : params - } - if savepointinfo: - tosave['tp_points'] = tpdict['points'], - tosave['tp_si'] = tpdict['si'], - tosave['tp_mi'] = tpdict['mi'], - tosave['tp_nmi'] = tpdict['nmi'] - tosave['tphc_points'] = tphcdict['points'], - tosave['tphc_ti'] = tphcdict['ti'], - tosave['tphc_mi'] = tphcdict['mi'], - tosave['tphc_nmi'] = tphcdict['nmi'] - np.savez(file, **tosave) - -def hexlattice_precalc_AB_loadunwrap(file, tpdict = None, tphcdict = None, return_points = False): - npz = np.load(file) - precalc_params = npz['precalc_params'][()] - my, ny = get_mn_y(precalc_params['lMax']) - nelem = len(my) - # this I should have made more universal... - if precalc_params['savepointinfo']: - if not tpdict: - tpdict = { - 'points' : npz['tp_points'], - 'si' : npz['tp_si'], - 'mi' : npz['tp_mi'], - 'nmi' : npz['tp_nmi'], - } - if not tphcdict: - tphcdict = { - 'points' : npz['tphc_points'], - 'ti' : npz['tphc_ti'], - 'mi' : npz['tphc_mi'], - 'nmi' : npz['tphc_nmi'] - } - else: - if not tpdict: - tpdict = generate_trianglepoints(maxlayer = precalc_params['maxlayer'], v3d=True, - circular=precalc_params['circular'], sixthindices=True, mirrorindices=True) - if not tphcdict: - tphcdict = generate_trianglepoints_hexcomplement(maxlayer=precalc_params['maxlayer'], v3d=True, - circular=precalc_params['circular'], thirdindices=True, mirrorindices=True) - - # For some obscure reason, I keep getting trailing single-dimension in the beginning for these arrays - for a in (tpdict['points'], tphcdict['points'], tpdict['si'], tpdict['mi'], - tphcdict['ti'], tphcdict['mi']): - if len(a.shape) > 2: - a.shape = a.shape[1::] - - self_tr = tpdict['points'] - d2u_tr = tphcdict['points'] - if len(self_tr.shape)>2: - self_tr = np.reshape(self_tr, self_tr.shape[1::]) - if len(d2u_tr.shape)>2: - d2u_tr = np.reshape(d2u_tr, d2u_tr.shape[1::]) - u2d_tr = -d2u_tr - a_self = np.empty((self_tr.shape[0],nelem,nelem), dtype=complex) - b_self = np.empty((self_tr.shape[0],nelem,nelem), dtype=complex) - a_d2u = np.empty(( d2u_tr.shape[0],nelem,nelem), dtype=complex) - b_d2u = np.empty(( d2u_tr.shape[0],nelem,nelem), dtype=complex) - a_self[tpdict['nmi']]=npz['a_self_nm'] - a_self[tpdict['mi'][0]]=npz['a_self_m0'] - b_self[tpdict['nmi']]=npz['b_self_nm'] - b_self[tpdict['mi'][0]]=npz['b_self_m0'] - mirrorangles = cart2sph(self_tr[tpdict['mi'][1]])[:,2] - cart2sph(self_tr[tpdict['mi'][0]])[:,2] - a_self[tpdict['mi'][1],:,:] = a_self[tpdict['mi'][0],:,:] * np.exp(1j*mirrorangles[:,nx,nx]*(my[nx,nx,:]-my[nx,:,nx])) - b_self[tpdict['mi'][1],:,:] = b_self[tpdict['mi'][0],:,:] * np.exp(1j*mirrorangles[:,nx,nx]*(my[nx,nx,:]-my[nx,:,nx])) - for i in range(1,6): - a_self[tpdict['si'][i],:,:] = a_self[tpdict['si'][0],:,:] * np.exp(1j*math.pi/3*i*(my[nx,:]-my[:,nx])) - b_self[tpdict['si'][i],:,:] = b_self[tpdict['si'][0],:,:] * np.exp(1j*math.pi/3*i*(my[nx,:]-my[:,nx])) - a_d2u[tphcdict['nmi']]=npz['a_d2u_nm'] - a_d2u[tphcdict['mi'][0]]=npz['a_d2u_m0'] - b_d2u[tphcdict['nmi']]=npz['b_d2u_nm'] - b_d2u[tphcdict['mi'][0]]=npz['b_d2u_m0'] - mirrorangles = cart2sph(self_tr[tphcdict['mi'][1]])[:,2] - cart2sph(self_tr[tphcdict['mi'][0]])[:,2] - a_d2u[tphcdict['mi'][1],:,:] = a_d2u[tphcdict['mi'][0],:,:] * np.exp(1j*mirrorangles[:,nx,nx]*(my[nx,nx,:]-my[nx,:,nx])) - b_d2u[tphcdict['mi'][1],:,:] = b_d2u[tphcdict['mi'][0],:,:] * np.exp(1j*mirrorangles[:,nx,nx]*(my[nx,nx,:]-my[nx,:,nx])) - for i in (1,-1): - a_d2u[tphcdict['ti'][i],:,:] = a_d2u[tphcdict['ti'][0],:,:] * np.exp(i*2j*math.pi/3*(my[nx,:]-my[:,nx])) - b_d2u[tphcdict['ti'][i],:,:] = b_d2u[tphcdict['ti'][0],:,:] * np.exp(i*2j*math.pi/3*(my[nx,:]-my[:,nx])) - a_u2d = a_d2u * (-1)**(my[nx,:]-my[:,nx]) - b_u2d = b_d2u * (-1)**(my[nx,:]-my[:,nx]) - d = { - 'a_self' : a_self, - 'b_self' : b_self, - 'a_d2u' : a_d2u, - 'b_d2u' : b_d2u, - 'a_u2d' : a_u2d, - 'b_u2d' : b_u2d, - } - for k in precalc_params.keys(): - d[k] = precalc_params[k] - if return_points: - d['d2u_tr'] = tphcdict['points'] - d['u2d_tr'] = -tphcdict['points'] - d['self_tr'] = tpdict['points'] - return d - -def hexlattice_get_AB(lMax, k_hexside, maxlayer, circular=True, return_points = True, J_scat=3): - params = { - 'lMax' : lMax, - 'k_hexside' : k_hexside, - 'maxlayer' : maxlayer, - 'circular' : circular, - 'savepointinfo' : return_points, # should I delete this key? - 'J_scat' : J_scat - } - tpdict = generate_trianglepoints(maxlayer, v3d=True, circular=circular, sixthindices=True, mirrorindices=True) - tphcdict = generate_trianglepoints_hexcomplement(maxlayer, v3d=True, circular=circular, thirdindices=True, mirrorindices=True) - my, ny = get_mn_y(lMax) - nelem = len(my) - a_self_nm = np.empty((tpdict['nmi'].shape[0],nelem,nelem), dtype=complex) - b_self_nm = np.empty((tpdict['nmi'].shape[0],nelem,nelem), dtype=complex) - a_self_m0 = np.empty((tpdict['mi'].shape[1],nelem,nelem), dtype=complex) - b_self_m0 = np.empty((tpdict['mi'].shape[1],nelem,nelem), dtype=complex) - a_d2u_nm = np.empty((tphcdict['nmi'].shape[0],nelem,nelem), dtype=complex) - b_d2u_nm = np.empty((tphcdict['nmi'].shape[0],nelem,nelem), dtype=complex) - a_d2u_m0 = np.empty((tphcdict['mi'].shape[1],nelem,nelem), dtype=complex) - b_d2u_m0 = np.empty((tphcdict['mi'].shape[1],nelem,nelem), dtype=complex) - - k_0 = k_hexside*_s3 # not really a wave vector here because of the normalisation! - tc = trans_calculator(lMax) - - y = np.arange(nelem) - - points = tpdict['points'][tpdict['nmi']] - d_i2j = cart2sph(points) - a_self_nm, b_self_nm = tc.get_AB_arrays(k_0*d_i2j[:,0],d_i2j[:,1],d_i2j[:,2],np.array([False]),J_scat) - - points = tpdict['points'][tpdict['mi'][0]] - d_i2j = cart2sph(points) - a_self_m0, b_self_m0 = tc.get_AB_arrays(k_0*d_i2j[:,0],d_i2j[:,1],d_i2j[:,2],np.array([False]),J_scat) - - points = tphcdict['points'][tphcdict['nmi']] - d_i2j = cart2sph(points) - a_d2u_nm, b_d2u_nm = tc.get_AB_arrays(k_0*d_i2j[:,0],d_i2j[:,1],d_i2j[:,2],np.array([False]),J_scat) - - points = tphcdict['points'][tphcdict['mi'][0]] - d_i2j = cart2sph(points) - a_d2u_m0, b_d2u_m0 = tc.get_AB_arrays(k_0*d_i2j[:,0],d_i2j[:,1],d_i2j[:,2],np.array([False]),J_scat) - ''' - tosave = { - 'a_self_nm' : a_self_nm, - 'a_self_m0' : a_self_m0, - 'b_self_nm' : b_self_nm, - 'b_self_m0' : b_self_m0, - 'a_d2u_nm' : a_d2u_nm, - 'a_d2u_m0' : a_d2u_m0, - 'b_d2u_nm' : b_d2u_nm, - 'b_d2u_m0' : b_d2u_m0, - 'precalc_params' : params - } - if savepointinfo: - tosave['tp_points'] = tpdict['points'], - tosave['tp_si'] = tpdict['si'], - tosave['tp_mi'] = tpdict['mi'], - tosave['tp_nmi'] = tpdict['nmi'] - tosave['tphc_points'] = tphcdict['points'], - tosave['tphc_ti'] = tphcdict['ti'], - tosave['tphc_mi'] = tphcdict['mi'], - tosave['tphc_nmi'] = tphcdict['nmi'] - np.savez(file, **tosave) - ''' - self_tr = tpdict['points'] - d2u_tr = tphcdict['points'] - if len(self_tr.shape)>2: - self_tr = np.reshape(self_tr, self_tr.shape[1::]) - if len(d2u_tr.shape)>2: - d2u_tr = np.reshape(d2u_tr, d2u_tr.shape[1::]) - u2d_tr = -d2u_tr - a_self = np.empty((self_tr.shape[0],nelem,nelem), dtype=complex) - b_self = np.empty((self_tr.shape[0],nelem,nelem), dtype=complex) - a_d2u = np.empty(( d2u_tr.shape[0],nelem,nelem), dtype=complex) - b_d2u = np.empty(( d2u_tr.shape[0],nelem,nelem), dtype=complex) - a_self[tpdict['nmi']]=a_self_nm - a_self[tpdict['mi'][0]]=a_self_m0 - b_self[tpdict['nmi']]=b_self_nm - b_self[tpdict['mi'][0]]=b_self_m0 - mirrorangles = cart2sph(self_tr[tpdict['mi'][1]])[:,2] - cart2sph(self_tr[tpdict['mi'][0]])[:,2] - a_self[tpdict['mi'][1],:,:] = a_self[tpdict['mi'][0],:,:] * np.exp(1j*mirrorangles[:,nx,nx]*(my[nx,nx,:]-my[nx,:,nx])) - b_self[tpdict['mi'][1],:,:] = b_self[tpdict['mi'][0],:,:] * np.exp(1j*mirrorangles[:,nx,nx]*(my[nx,nx,:]-my[nx,:,nx])) - for i in range(1,6): - a_self[tpdict['si'][i],:,:] = a_self[tpdict['si'][0],:,:] * np.exp(1j*math.pi/3*i*(my[nx,:]-my[:,nx])) - b_self[tpdict['si'][i],:,:] = b_self[tpdict['si'][0],:,:] * np.exp(1j*math.pi/3*i*(my[nx,:]-my[:,nx])) - a_d2u[tphcdict['nmi']]=a_d2u_nm - a_d2u[tphcdict['mi'][0]]=a_d2u_m0 - b_d2u[tphcdict['nmi']]=b_d2u_nm - b_d2u[tphcdict['mi'][0]]=b_d2u_m0 - mirrorangles = cart2sph(self_tr[tphcdict['mi'][1]])[:,2] - cart2sph(self_tr[tphcdict['mi'][0]])[:,2] - a_d2u[tphcdict['mi'][1],:,:] = a_d2u[tphcdict['mi'][0],:,:] * np.exp(1j*mirrorangles[:,nx,nx]*(my[nx,nx,:]-my[nx,:,nx])) - b_d2u[tphcdict['mi'][1],:,:] = b_d2u[tphcdict['mi'][0],:,:] * np.exp(1j*mirrorangles[:,nx,nx]*(my[nx,nx,:]-my[nx,:,nx])) - for i in (1,-1): - a_d2u[tphcdict['ti'][i],:,:] = a_d2u[tphcdict['ti'][0],:,:] * np.exp(i*2j*math.pi/3*(my[nx,:]-my[:,nx])) - b_d2u[tphcdict['ti'][i],:,:] = b_d2u[tphcdict['ti'][0],:,:] * np.exp(i*2j*math.pi/3*(my[nx,:]-my[:,nx])) - a_u2d = a_d2u * (-1)**(my[nx,:]-my[:,nx]) - b_u2d = b_d2u * (-1)**(my[nx,:]-my[:,nx]) - d = { - 'a_self' : a_self, - 'b_self' : b_self, - 'a_d2u' : a_d2u, - 'b_d2u' : b_d2u, - 'a_u2d' : a_u2d, - 'b_u2d' : b_u2d, - } - for k in params.keys(): - d[k] = params[k] - if return_points: - d['d2u_tr'] = tphcdict['points'] - d['u2d_tr'] = -tphcdict['points'] - d['self_tr'] = tpdict['points'] - return d - -from scipy.constants import c -from .timetrack import _time_b, _time_e -from .tmatrices import symz_indexarrays - -def hexlattice_zsym_getSVD(lMax, TMatrices_om, epsilon_b, hexside, maxlayer, omega, klist, gaussianSigma=False, onlyNmin=0, verbose=False): - btime = _time_b(verbose) - nelem = lMax * (lMax + 2) - n2id = np.identity(2*nelem) - n2id.shape = (2,nelem,2,nelem) - nan = float('nan') - k_0 = omega * math.sqrt(epsilon_b) / c - tdic = hexlattice_get_AB(lMax,k_0*hexside,maxlayer) - a_self = tdic['a_self'][:,:nelem,:nelem] - b_self = tdic['b_self'][:,:nelem,:nelem] - a_u2d = tdic['a_u2d'][:,:nelem,:nelem] - b_u2d = tdic['b_u2d'][:,:nelem,:nelem] - a_d2u = tdic['a_d2u'][:,:nelem,:nelem] - b_d2u = tdic['b_d2u'][:,:nelem,:nelem] - unitcell_translations = tdic['self_tr']*hexside*_s3 - u2d_translations = tdic['u2d_tr']*hexside*_s3 - d2u_translations = tdic['d2u_tr']*hexside*_s3 - - if gaussianSigma: - sbtime = _time_b(verbose, step='Calculating gaussian envelope') - unitcell_envelope = np.exp(-np.sum(tdic['self_tr']**2,axis=-1)/(2*gaussianSigma**2)) - u2d_envelope = np.exp(-np.sum(tdic['u2d_tr']**2,axis=-1)/(2*gaussianSigma**2)) - d2u_envelope = np.exp(-np.sum(tdic['d2u_tr']**2,axis=-1)/(2*gaussianSigma**2)) - _time_e(sbtime, verbose, step='Calculating gaussian envelope') - - #TMatrices_om = TMatrices_interp(omega) - if(not onlyNmin): - svUfullTElist = np.full((klist.shape[0], 2*nelem, 2*nelem), np.nan, dtype=complex) - svVfullTElist = np.full((klist.shape[0], 2*nelem, 2*nelem), np.nan, dtype=complex) - svSfullTElist = np.full((klist.shape[0], 2*nelem), np.nan, dtype=complex) - svUfullTMlist = np.full((klist.shape[0], 2*nelem, 2*nelem), np.nan, dtype=complex) - svVfullTMlist = np.full((klist.shape[0], 2*nelem, 2*nelem), np.nan, dtype=complex) - svSfullTMlist = np.full((klist.shape[0], 2*nelem), np.nan, dtype=complex) - else: - minsvTElist = np.full((klist.shape[0], onlyNmin),np.nan) - minsvTMlist = np.full((klist.shape[0], onlyNmin),np.nan) - - leftmatrixlist = np.full((klist.shape[0],2,2,nelem,2,2,nelem),np.nan,dtype=complex) - #isNaNlist = np.zeros((klist.shape[0]), dtype=bool) - isNaNlist = (k_0*k_0 - klist[:,0]**2 - klist[:,1]**2 < 0) - nnlist = np.logical_not(isNaNlist) - - sbtime = _time_b(verbose, step='Initialization of matrices for SVD for a given list of k\'s') - - - #ki = np.arange(klist.shape[0])[k_0*k_0 - klist[:,0]**2 - klist[:,1]**2 >= 0] - k = klist[nnlist] - - phases_self = np.exp(1j*np.tensordot(k,unitcell_translations,axes=(-1,-1))) - phases_u2d = np.exp(1j*np.tensordot(k,u2d_translations,axes=(-1,-1))) - phases_d2u = np.exp(1j*np.tensordot(k,d2u_translations,axes=(-1,-1))) - if gaussianSigma: - phases_self *= unitcell_envelope - phases_u2d *= u2d_envelope - phases_d2u *= d2u_envelope - leftmatrix = np.zeros((k.shape[0],2,2,nelem, 2,2,nelem), dtype=complex) - # 0:[u,E<--u,E ] - # 1:[d,M<--d,M ] - leftmatrix[:,0,0,:,0,0,:] = np.tensordot(phases_self,a_self, axes=(-1,0)) # u2u, E2E - leftmatrix[:,1,0,:,1,0,:] = leftmatrix[:,0,0,:,0,0,:] # d2d, E2E - leftmatrix[:,0,1,:,0,1,:] = leftmatrix[:,0,0,:,0,0,:] # u2u, M2M - leftmatrix[:,1,1,:,1,1,:] = leftmatrix[:,0,0,:,0,0,:] # d2d, M2M - leftmatrix[:,0,0,:,0,1,:] = np.tensordot(phases_self,b_self, axes=(-1,0)) # u2u, M2E - leftmatrix[:,0,1,:,0,0,:] = leftmatrix[:,0,0,:,0,1,:] # u2u, E2M - leftmatrix[:,1,1,:,1,0,:] = leftmatrix[:,0,0,:,0,1,:] # d2d, E2M - leftmatrix[:,1,0,:,1,1,:] = leftmatrix[:,0,0,:,0,1,:] # d2d, M2E - leftmatrix[:,0,0,:,1,0,:] = np.tensordot(phases_d2u, a_d2u,axes=(-1,0)) #d2u,E2E - leftmatrix[:,0,1,:,1,1,:] = leftmatrix[:,0,0,:,1,0,:] #d2u, M2M - leftmatrix[:,1,0,:,0,0,:] = np.tensordot(phases_u2d, a_u2d,axes=(-1,0)) #u2d,E2E - leftmatrix[:,1,1,:,0,1,:] = leftmatrix[:,1,0,:,0,0,:] #u2d, M2M - leftmatrix[:,0,0,:,1,1,:] = np.tensordot(phases_d2u, b_d2u,axes=(-1,0)) #d2u,M2E - leftmatrix[:,0,1,:,1,0,:] = leftmatrix[:,0,0,:,1,1,:] #d2u, E2M - leftmatrix[:,1,0,:,0,1,:] = np.tensordot(phases_u2d, b_u2d,axes=(-1,0)) #u2d,M2E - leftmatrix[:,1,1,:,0,0,:] = leftmatrix[:,1,0,:,0,1,:] #u2d, E2M - #leftmatrix is now the translation matrix T - for j in range(2): - leftmatrix[:,j] = np.rollaxis(-np.tensordot(TMatrices_om[j], leftmatrix[:,j], axes=([-2,-1],[1,2])),2) - # at this point, jth row of leftmatrix is that of -MT - leftmatrix[:,j,:,:,j,:,:] += n2id - #now we are done, 1-MT - - leftmatrixlist[nnlist] = leftmatrix - - ''' - # sem nějaká rozumná smyčka - for ki in range(klist.shape[0]): - k = klist[ki] - if (k_0*k_0 - k[0]*k[0] - k[1]*k[1] < 0): - isNaNlist[ki] = True - continue - - phases_self = np.exp(1j*np.tensordot(k,unitcell_translations,axes=(0,-1))) - phases_u2d = np.exp(1j*np.tensordot(k,u2d_translations,axes=(0,-1))) - phases_d2u = np.exp(1j*np.tensordot(k,d2u_translations,axes=(0,-1))) - if gaussianSigma: - phases_self *= unitcell_envelope - phases_u2d *= u2d_envelope - phases_d2u *= d2u_envelope - leftmatrix = np.zeros((2,2,nelem, 2,2,nelem), dtype=complex) - # 0:[u,E<--u,E ] - # 1:[d,M<--d,M ] - leftmatrix[0,0,:,0,0,:] = np.tensordot(a_self,phases_self, axes=(0,-1)) # u2u, E2E - leftmatrix[1,0,:,1,0,:] = leftmatrix[0,0,:,0,0,:] # d2d, E2E - leftmatrix[0,1,:,0,1,:] = leftmatrix[0,0,:,0,0,:] # u2u, M2M - leftmatrix[1,1,:,1,1,:] = leftmatrix[0,0,:,0,0,:] # d2d, M2M - leftmatrix[0,0,:,0,1,:] = np.tensordot(b_self,phases_self, axes=(0,-1)) # u2u, M2E - leftmatrix[0,1,:,0,0,:] = leftmatrix[0,0,:,0,1,:] # u2u, E2M - leftmatrix[1,1,:,1,0,:] = leftmatrix[0,0,:,0,1,:] # d2d, E2M - leftmatrix[1,0,:,1,1,:] = leftmatrix[0,0,:,0,1,:] # d2d, M2E - leftmatrix[0,0,:,1,0,:] = np.tensordot(a_d2u, phases_d2u,axes=(0,-1)) #d2u,E2E - leftmatrix[0,1,:,1,1,:] = leftmatrix[0,0,:,1,0,:] #d2u, M2M - leftmatrix[1,0,:,0,0,:] = np.tensordot(a_u2d, phases_u2d,axes=(0,-1)) #u2d,E2E - leftmatrix[1,1,:,0,1,:] = leftmatrix[1,0,:,0,0,:] #u2d, M2M - leftmatrix[0,0,:,1,1,:] = np.tensordot(b_d2u, phases_d2u,axes=(0,-1)) #d2u,M2E - leftmatrix[0,1,:,1,0,:] = leftmatrix[0,0,:,1,1,:] #d2u, E2M - leftmatrix[1,0,:,0,1,:] = np.tensordot(b_u2d, phases_u2d,axes=(0,-1)) #u2d,M2E - leftmatrix[1,1,:,0,0,:] = leftmatrix[1,0,:,0,1,:] #u2d, E2M - #leftmatrix is now the translation matrix T - for j in range(2): - leftmatrix[j] = -np.tensordot(TMatrices_om[j], leftmatrix[j], axes=([-2,-1],[0,1])) - # at this point, jth row of leftmatrix is that of -MT - leftmatrix[j,:,:,j,:,:] += n2id - #now we are done, 1-MT - - leftmatrixlist[ki] = leftmatrix - ''' - - leftmatrixlist_s = np.reshape(leftmatrixlist,(klist.shape[0], 2*2*nelem,2*2*nelem))[nnlist] - TEč, TMč = symz_indexarrays(lMax, 2) - leftmatrixlist_TE = leftmatrixlist_s[np.ix_(np.arange(leftmatrixlist_s.shape[0]),TEč,TEč)] - leftmatrixlist_TM = leftmatrixlist_s[np.ix_(np.arange(leftmatrixlist_s.shape[0]),TMč,TMč)] - _time_e(sbtime, verbose, step='Initializing matrices for SVD for a given list of k\'s') - - sbtime = _time_b(verbose, step='Calculating SVDs for a given list of k\'s.') - if(not onlyNmin): - svUfullTElist[nnlist], svSfullTElist[nnlist], svVfullTElist[nnlist] = np.linalg.svd(leftmatrixlist_TE, compute_uv=True) - svUfullTMlist[nnlist], svSfullTMlist[nnlist], svVfullTMlist[nnlist] = np.linalg.svd(leftmatrixlist_TM, compute_uv=True) - _time_e(sbtime, verbose, step='Calculating SVDs for a given list of k\'s.') - _time_e(btime, verbose) - return ((svUfullTElist, svSfullTElist, svVfullTElist), (svUfullTMlist, svSfullTMlist, svVfullTMlist)) - else: - minsvTElist[nnlist] = np.linalg.svd(leftmatrixlist_TE, compute_uv=False)[...,-onlyNmin:] - minsvTMlist[nnlist] = np.linalg.svd(leftmatrixlist_TM, compute_uv=False)[...,-onlyNmin:] - _time_e(sbtime, verbose, step='Calculating SVDs for a given list of k\'s.') - _time_e(btime, verbose) - return (minsvTElist, minsvTMlist) diff --git a/qpms/hexpoints_c.pyx b/qpms/hexpoints_c.pyx deleted file mode 100644 index ebcf0fb..0000000 --- a/qpms/hexpoints_c.pyx +++ /dev/null @@ -1,130 +0,0 @@ -import math -import numpy as np -cimport numpy as np -nx = None - -cdef double _s3 = math.sqrt(3) - -from scipy.constants import c -from .timetrack import _time_b, _time_e -from .tmatrices import symz_indexarrays -from .hexpoints import hexlattice_get_AB - -cpdef hexlattice_zsym_getSVD(int lMax, TMatrices_om, double epsilon_b, double hexside, size_t maxlayer, double omega, klist, gaussianSigma=False, int onlyNmin=0, verbose=False): - cdef np.ndarray[np.npy_double, ndim = 2] klist_c = klist - btime = _time_b(verbose) - cdef size_t nelem = lMax * (lMax + 2) - _n2id = np.identity(2*nelem) - _n2id.shape = (2,nelem,2,nelem) - cdef np.ndarray[np.npy_double, ndim = 4] n2id = _n2id - cdef double nan = float('nan') - k_0 = omega * math.sqrt(epsilon_b) / c - tdic = hexlattice_get_AB(lMax,k_0*hexside,maxlayer) - cdef np.ndarray[np.npy_cdouble, ndim = 3] a_self = tdic['a_self'][:,:nelem,:nelem] - cdef np.ndarray[np.npy_cdouble, ndim = 3] b_self = tdic['b_self'][:,:nelem,:nelem] - cdef np.ndarray[np.npy_cdouble, ndim = 3] a_u2d = tdic['a_u2d'][:,:nelem,:nelem] - cdef np.ndarray[np.npy_cdouble, ndim = 3] b_u2d = tdic['b_u2d'][:,:nelem,:nelem] - cdef np.ndarray[np.npy_cdouble, ndim = 3] a_d2u = tdic['a_d2u'][:,:nelem,:nelem] - cdef np.ndarray[np.npy_cdouble, ndim = 3] b_d2u = tdic['b_d2u'][:,:nelem,:nelem] - cdef np.ndarray[np.npy_double, ndim = 2] unitcell_translations = tdic['self_tr']*hexside*_s3 - cdef np.ndarray[np.npy_double, ndim = 2] u2d_translations = tdic['u2d_tr']*hexside*_s3 - cdef np.ndarray[np.npy_double, ndim = 2] d2u_translations = tdic['d2u_tr']*hexside*_s3 - - cdef np.ndarray[np.npy_double, ndim = 1] unitcell_envelope, u2d_envelope, d2u_envelope - if gaussianSigma: - sbtime = _time_b(verbose, step='Calculating gaussian envelope') - unitcell_envelope = np.exp(-np.sum(tdic['self_tr']**2,axis=-1)/(2*gaussianSigma**2)) - u2d_envelope = np.exp(-np.sum(tdic['u2d_tr']**2,axis=-1)/(2*gaussianSigma**2)) - d2u_envelope = np.exp(-np.sum(tdic['d2u_tr']**2,axis=-1)/(2*gaussianSigma**2)) - _time_e(sbtime, verbose, step='Calculating gaussian envelope') - - cdef np.ndarray[np.npy_cdouble, ndim = 3] svUfullTElist, svVfullTElist, svUfullTMlist, svVfullTMlist - cdef np.ndarray[np.npy_cdouble, ndim = 2] svSfullTElist, svSfullTMlist - cdef np.ndarray[np.npy_double, ndim = 2] minsvTElist, minsvTMlist - - #TMatrices_om = TMatrices_interp(omega) - if(not onlyNmin): - svUfullTElist = np.full((klist_c.shape[0], 2*nelem, 2*nelem), np.nan, dtype=complex) - svVfullTElist = np.full((klist_c.shape[0], 2*nelem, 2*nelem), np.nan, dtype=complex) - svSfullTElist = np.full((klist_c.shape[0], 2*nelem), np.nan, dtype=complex) - svUfullTMlist = np.full((klist_c.shape[0], 2*nelem, 2*nelem), np.nan, dtype=complex) - svVfullTMlist = np.full((klist_c.shape[0], 2*nelem, 2*nelem), np.nan, dtype=complex) - svSfullTMlist = np.full((klist_c.shape[0], 2*nelem), np.nan, dtype=complex) - else: - minsvTElist = np.full((klist_c.shape[0], onlyNmin),np.nan) - minsvTMlist = np.full((klist_c.shape[0], onlyNmin),np.nan) - - cdef np.ndarray[np.npy_cdouble] leftmatrixlist = np.full((klist_c.shape[0],2,2,nelem,2,2,nelem),np.nan,dtype=complex) - cdef np.ndarray[np.npy_bool, ndim=1] isNaNlist = np.zeros((klist_c.shape[0]), dtype=bool) - - sbtime = _time_b(verbose, step='Initialization of matrices for SVD for a given list of k\'s') - # sem nějaká rozumná smyčka - - # declarations for the ki loop: - cdef size_t ki - cdef np.ndarray[np.npy_cdouble, ndim = 1] phases_self - cdef np.ndarray[np.npy_cdouble, ndim = 1] phases_u2d - cdef np.ndarray[np.npy_cdouble, ndim = 1] phases_d2u - cdef np.ndarray[np.npy_cdouble, ndim=6] leftmatrix - cdef np.ndarray[np.npy_double, ndim=1] k - cdef int j - - for ki in range(klist_c.shape[0]): - k = klist_c[ki] - if (k_0*k_0 - k[0]*k[0] - k[1]*k[1] < 0): - isNaNlist[ki] = True - continue - - phases_self = np.exp(1j*np.tensordot(k,unitcell_translations,axes=(0,-1))) - phases_u2d = np.exp(1j*np.tensordot(k,u2d_translations,axes=(0,-1))) - phases_d2u = np.exp(1j*np.tensordot(k,d2u_translations,axes=(0,-1))) - if gaussianSigma: - phases_self *= unitcell_envelope - phases_u2d *= u2d_envelope - phases_d2u *= d2u_envelope - leftmatrix = np.zeros((2,2,nelem, 2,2,nelem), dtype=complex) - # 0:[u,E<--u,E ] - # 1:[d,M<--d,M ] - leftmatrix[0,0,:,0,0,:] = np.tensordot(a_self,phases_self, axes=(0,-1)) # u2u, E2E - leftmatrix[1,0,:,1,0,:] = leftmatrix[0,0,:,0,0,:] # d2d, E2E - leftmatrix[0,1,:,0,1,:] = leftmatrix[0,0,:,0,0,:] # u2u, M2M - leftmatrix[1,1,:,1,1,:] = leftmatrix[0,0,:,0,0,:] # d2d, M2M - leftmatrix[0,0,:,0,1,:] = np.tensordot(b_self,phases_self, axes=(0,-1)) # u2u, M2E - leftmatrix[0,1,:,0,0,:] = leftmatrix[0,0,:,0,1,:] # u2u, E2M - leftmatrix[1,1,:,1,0,:] = leftmatrix[0,0,:,0,1,:] # d2d, E2M - leftmatrix[1,0,:,1,1,:] = leftmatrix[0,0,:,0,1,:] # d2d, M2E - leftmatrix[0,0,:,1,0,:] = np.tensordot(a_d2u, phases_d2u,axes=(0,-1)) #d2u,E2E - leftmatrix[0,1,:,1,1,:] = leftmatrix[0,0,:,1,0,:] #d2u, M2M - leftmatrix[1,0,:,0,0,:] = np.tensordot(a_u2d, phases_u2d,axes=(0,-1)) #u2d,E2E - leftmatrix[1,1,:,0,1,:] = leftmatrix[1,0,:,0,0,:] #u2d, M2M - leftmatrix[0,0,:,1,1,:] = np.tensordot(b_d2u, phases_d2u,axes=(0,-1)) #d2u,M2E - leftmatrix[0,1,:,1,0,:] = leftmatrix[0,0,:,1,1,:] #d2u, E2M - leftmatrix[1,0,:,0,1,:] = np.tensordot(b_u2d, phases_u2d,axes=(0,-1)) #u2d,M2E - leftmatrix[1,1,:,0,0,:] = leftmatrix[1,0,:,0,1,:] #u2d, E2M - #leftmatrix is now the translation matrix T - for j in range(2): - leftmatrix[j] = -np.tensordot(TMatrices_om[j], leftmatrix[j], axes=([-2,-1],[0,1])) - # at this point, jth row of leftmatrix is that of -MT - leftmatrix[j,:,:,j,:,:] += n2id - #now we are done, 1-MT - - leftmatrixlist[ki] = leftmatrix - - nnlist = np.logical_not(isNaNlist) - leftmatrixlist_s = np.reshape(leftmatrixlist,(klist_c.shape[0], 2*2*nelem,2*2*nelem))[nnlist] - TEc, TMc = symz_indexarrays(lMax, 2) - leftmatrixlist_TE = leftmatrixlist_s[np.ix_(np.arange(leftmatrixlist_s.shape[0]),TEc,TEc)] - leftmatrixlist_TM = leftmatrixlist_s[np.ix_(np.arange(leftmatrixlist_s.shape[0]),TMc,TMc)] - _time_e(sbtime, verbose, step='Initializing matrices for SVD for a given list of k\'s') - - sbtime = _time_b(verbose, step='Calculating SVDs for a given list of k\'s.') - if(not onlyNmin): - svUfullTElist[nnlist], svSfullTElist[nnlist], svVfullTElist[nnlist] = np.linalg.svd(leftmatrixlist_TE, compute_uv=True) - svUfullTMlist[nnlist], svSfullTMlist[nnlist], svVfullTMlist[nnlist] = np.linalg.svd(leftmatrixlist_TM, compute_uv=True) - _time_e(sbtime, verbose, step='Calculating SVDs for a given list of k\'s.') - return ((svUfullTElist, svSfullTElist, svVfullTElist), (svUfullTMlist, svSfullTMlist, svVfullTMlist)) - else: - minsvTElist[nnlist] = np.linalg.svd(leftmatrixlist_TE, compute_uv=False)[...,-onlyNmin:] - minsvTMlist[nnlist] = np.linalg.svd(leftmatrixlist_TM, compute_uv=False)[...,-onlyNmin:] - _time_e(sbtime, verbose, step='Calculating SVDs for a given list of k\'s.') - return (minsvTElist, minsvTMlist) diff --git a/qpms/lattices2d.py b/qpms/lattices2d.py index eb5f68b..b5470f6 100644 --- a/qpms/lattices2d.py +++ b/qpms/lattices2d.py @@ -1,3 +1,7 @@ +""" +These functions are mostly deprecated by the C counterparts from lattices.h. +This file is still kept for reference, but might be removed in the future. +""" import numpy as np from enum import Enum from math import floor diff --git a/qpms/legacy.py b/qpms/legacy.py deleted file mode 100644 index 3fb2d71..0000000 --- a/qpms/legacy.py +++ /dev/null @@ -1,1058 +0,0 @@ -import math -import numpy as np -nx = None - -_s3 = math.sqrt(3) - -from qpms_c import get_mn_y, trans_calculator -from .qpms_p import cart2sph - - - - -# Functions copied from scattering_xu, additionaly normalized -from py_gmm.gmm import vec_trans as vc - -#@ujit -def q_max(m,n,μ,ν): - return min(n,ν,(n+ν-abs(m+μ))/2) - -q_max_v = np.vectorize(q_max) - -# returns array with indices corresponding to q -# argument q does nothing for now -#@ujit -def a_q(m,n,μ,ν,q = None): - qm=q_max(m,n,μ,ν) - res, err= vc.gaunt_xu(m,n,μ,ν,qm) - if(err): - print("m,n,μ,ν,qm = ",m,n,μ,ν,qm) - raise ValueError('Something bad in the fortran subroutine gaunt_xu happened') - return res - -a_q_v = np.vectorize(a_q) - - -# All arguments are single numbers (for now) -# ZDE VYCHÁZEJÍ DIVNÁ ZNAMÉNKA -#@ujit -def Ã(m,n,μ,ν,kdlj,θlj,φlj,r_ge_d,J): - """ - The à translation coefficient for spherical vector waves. - - Parameters - ---------- - m, n: int - The indices (degree and order) of the destination basis. - μ, ν: int - The indices of the source basis wave. - kdlj, θlj, φlj: float - The spherical coordinates of the relative position of - the new center vs. the old one (R_new - R_old); - the distance has to be already multiplied by the wavenumber! - r_ge_d: TODO - J: 1, 2, 3 or 4 - Type of the wave in the old center. - - Returns - ------- - TODO - - Bugs - ---- - gevero's gaunt coefficient implementation fails for large m, n (the unsafe territory - is somewhere around -72, 80) - - """ - exponent=(math.lgamma(2*n+1)-math.lgamma(n+2)+math.lgamma(2*ν+3)-math.lgamma(ν+2) - +math.lgamma(n+ν+m-μ+1)-math.lgamma(n-m+1)-math.lgamma(ν+μ+1) - +math.lgamma(n+ν+1) - math.lgamma(2*(n+ν)+1)) - presum = math.exp(exponent) - presum = presum * np.exp(1j*(μ-m)*φlj) * (-1)**m * 1j**(ν+n) / (4*n) - qmax = math.floor(q_max(-m,n,μ,ν)) #nemá tu být +m? - q = np.arange(qmax+1, dtype=int) - # N.B. -m !!!!!! - a1q = a_q(-m,n,μ,ν) # there is redundant calc. of qmax - ã1q = a1q / a1q[0] - p = n+ν-2*q - if(r_ge_d): - J = 1 - zp = zJn(n+ν,kdlj,J)[0][p] - Pp = lpmv(μ-m,p,math.cos(θlj)) - summandq = (n*(n+1) + ν*(ν+1) - p*(p+1)) * (-1)**q * ã1q * zp * Pp - - # Taylor normalisation v2, proven to be equivalent (NS which is better) - prenormratio = 1j**(ν-n) * math.sqrt(((2*ν+1)/(2*n+1))* math.exp( - math.lgamma(n+m+1)-math.lgamma(n-m+1)+math.lgamma(ν-μ+1)-math.lgamma(ν+μ+1))) - presum = presum / prenormratio - - # Taylor normalisation - #prenormmn = math.sqrt((2*n + 1)*math.factorial(n-m)/(4*π*factorial(n+m))) - #prenormμν = math.sqrt((2*ν + 1)*math.factorial(ν-μ)/(4*π*factorial(ν+μ))) - #presum = presum * prenormμν / prenormmn - - return presum * np.sum(summandq) - -# ZDE OPĚT JINAK ZNAMÉNKA než v Xu (J. comp. phys 127, 285) -#@ujit -def B̃(m,n,μ,ν,kdlj,θlj,φlj,r_ge_d,J): - """ - The B̃ translation coefficient for spherical vector waves. - - Parameters - ---------- - m, n: int - The indices (degree and order) of the destination basis. - μ, ν: int - The indices of the source basis wave. - kdlj, θlj, φlj: float - The spherical coordinates of the relative position of - the new center vs. the old one (R_new - R_old); - the distance has to be already multiplied by the wavenumber! - r_ge_d: TODO - J: 1, 2, 3 or 4 - Type of the wave in the old center. - - Returns: - -------- - TODO - """ - exponent=(math.lgamma(2*n+3)-math.lgamma(n+2)+math.lgamma(2*ν+3)-math.lgamma(ν+2) - +math.lgamma(n+ν+m-μ+2)-math.lgamma(n-m+1)-math.lgamma(ν+μ+1) - +math.lgamma(n+ν+2) - math.lgamma(2*(n+ν)+3)) - presum = math.exp(exponent) - presum = presum * np.exp(1j*(μ-m)*φlj) * (-1)**m * 1j**(ν+n+1) / ( - (4*n)*(n+1)*(n+m+1)) - Qmax = math.floor(q_max(-m,n+1,μ,ν)) - q = np.arange(Qmax+1, dtype=int) - if (μ == ν): # it would disappear in the sum because of the factor (ν-μ) anyway - ã2q = 0 - else: - a2q = a_q(-m-1,n+1,μ+1,ν) - ã2q = a2q / a2q[0] - a3q = a_q(-m,n+1,μ,ν) - ã3q = a3q / a3q[0] - #print(len(a2q),len(a3q)) - p = n+ν-2*q - if(r_ge_d): - J = 1 - zp_ = zJn(n+1+ν,kdlj,J)[0][p+1] # je ta +1 správně? - Pp_ = lpmv(μ-m,p+1,math.cos(θlj)) - summandq = ((2*(n+1)*(ν-μ)*ã2q - -(-ν*(ν+1) - n*(n+3) - 2*μ*(n+1)+p*(p+3))* ã3q) - *(-1)**q * zp_ * Pp_) - - # Taylor normalisation v2, proven to be equivalent - prenormratio = 1j**(ν-n) * math.sqrt(((2*ν+1)/(2*n+1))* math.exp( - math.lgamma(n+m+1)-math.lgamma(n-m+1)+math.lgamma(ν-μ+1)-math.lgamma(ν+μ+1))) - presum = presum / prenormratio - - ## Taylor normalisation - #prenormmn = math.sqrt((2*n + 1)*math.factorial(n-m)/(4*π*factorial(n+m))) - #prenormμν = math.sqrt((2*ν + 1)*math.factorial(ν-μ)/(4*π*factorial(ν+μ))) - #presum = presum * prenormμν / prenormmn - - return presum * np.sum(summandq) - - - - - - - - -#@ujit -def G_Mie_scat_precalc_cart(source_cart, dest_cart, RH, RV, a, nmax, k_i, k_e, μ_i=1, μ_e=1, J_ext=1, J_scat=3): - """ - r1_cart (destination), r2_cart (source) and the result are in cartesian coordinates - the result indices are in the source-destination order - TODO - """ - my, ny = get_mn_y(nmax) - nelem = len(my) - #source to origin - so_sph = cart2sph(-source_cart) - kd_so = k_e * so_sph[0] - θ_so = so_sph[1] - φ_so = so_sph[2] - # Decomposition of the source N_0,1, N_-1,1, and N_1,1 in the nanoparticle center - p_0 = np.empty((nelem), dtype=np.complex_) - q_0 = np.empty((nelem), dtype=np.complex_) - p_minus = np.empty((nelem), dtype=np.complex_) - q_minus = np.empty((nelem), dtype=np.complex_) - p_plus = np.empty((nelem), dtype=np.complex_) - q_plus = np.empty((nelem), dtype=np.complex_) - for y in range(nelem): - m = my[y] - n = ny[y] - p_0[y] = Ã(m,n, 0,1,kd_so,θ_so,φ_so,False,J=J_scat) - q_0[y] = B̃(m,n, 0,1,kd_so,θ_so,φ_so,False,J=J_scat) - p_minus[y] = Ã(m,n,-1,1,kd_so,θ_so,φ_so,False,J=J_scat) - q_minus[y] = B̃(m,n,-1,1,kd_so,θ_so,φ_so,False,J=J_scat) - p_plus[y] = Ã(m,n, 1,1,kd_so,θ_so,φ_so,False,J=J_scat) - q_plus[y] = B̃(m,n, 1,1,kd_so,θ_so,φ_so,False,J=J_scat) - a_0 = RV[ny] * p_0 - b_0 = RH[ny] * q_0 - a_plus = RV[ny] * p_plus - b_plus = RH[ny] * q_plus - a_minus = RV[ny] * p_minus - b_minus = RH[ny] * q_minus - orig2dest_sph = cart2sph(dest_cart) - orig2dest_sph[0] = k_e*orig2dest_sph[0] - M_dest_y, N_dest_y = vswf_yr1(orig2dest_sph,nmax,J=J_scat) - # N.B. these are in the local cartesian coordinates (r̂,θ̂,φ̂) - N_dest_0 = np.sum(a_0[:,ň] * N_dest_y, axis=-2) - M_dest_0 = np.sum(b_0[:,ň] * M_dest_y, axis=-2) - N_dest_plus = np.sum(a_plus[:,ň] * N_dest_y, axis=-2) - M_dest_plus = np.sum(b_plus[:,ň] * M_dest_y, axis=-2) - N_dest_minus = np.sum(a_minus[:,ň]* N_dest_y, axis=-2) - M_dest_minus = np.sum(b_minus[:,ň]* M_dest_y, axis=-2) - prefac = math.sqrt(1/(4*3*π))#/ε_0 - G_sourcez_dest = prefac * (N_dest_0+M_dest_0) - G_sourcex_dest = prefac * (N_dest_minus+M_dest_minus-N_dest_plus-M_dest_plus)/math.sqrt(2) - G_sourcey_dest = prefac * (N_dest_minus+M_dest_minus+N_dest_plus+M_dest_plus)/(1j*math.sqrt(2)) - G_source_dest = np.array([G_sourcex_dest, G_sourcey_dest, G_sourcez_dest]) - # To global cartesian coordinates: - G_source_dest = sph_loccart2cart(G_source_dest, sph=orig2dest_sph, axis=-1) - return G_source_dest - -#@ujit -def G_Mie_scat_cart(source_cart, dest_cart, a, nmax, k_i, k_e, μ_i=1, μ_e=1, J_ext=1, J_scat=3): - """ - TODO - """ - RH, RV, TH, TV = mie_coefficients(a=a, nmax=nmax, k_i=k_i, k_e=k_e, μ_i=μ_i, μ_e=μ_e, J_ext=J_ext, J_scat=J_scat) - return G_Mie_scat_precalc_cart_new(source_cart, dest_cart, RH, RV, a, nmax, k_i, k_e, μ_i, μ_e, J_ext, J_scat) - - -#TODO -def cross_section_Mie_precalc(): - pass - -def cross_section_Mie(a, nmax, k_i, k_e, μ_i, μ_e,): - pass - - - - -#################### -# Array simulations -#################### - - - -def scatter_plane_wave(omega, epsilon_b, positions, Tmatrices, k_dirs, E_0s, #saveto = None - ): - """ - Solves the plane wave linear scattering problem for a structure of "non-touching" particles - for one frequency and arbitrary number K of incoming plane waves. - - Parameters - ---------- - omega : positive number - The frequency of the field. - epsilon_b : complex number - Permittivity of the background medium (which has to be isotropic). - positions : (N,3)-shaped real array - Cartesian positions of the particles. - TMatrices : (N,2,nelem,2,nelem) or compatible - The T-matrices in the "Taylor convention" describing the scattering on a single nanoparticle. - If all the particles are identical and equally oriented, only one T-matrix can be given. - nelems = (lMax + 2) * lMax, where lMax is the highest multipole order to which the scattering - is calculated. - k_dirs : (K,3)-shaped real array or compatible - The direction of the incident field wave vector, normalized to one. - E_0s : (K,3)-shaped complex array or compatible - The electric intensity amplitude of the incident field. - - Returns - ------- - ab : (K, N, 2, nelem)-shaped complex array - The a (electric wave), b (magnetic wave) coefficients of the outgoing field for each particle - # Fuck this, it will be wiser to make separate function to calculate those from ab: - # sigma_xxx : TODO (K, 2, nelem) - # TODO partial (TODO which?) cross-section for each type of outgoing waves, summed over all - # nanoparticles (total cross section is given by the sum of this.) - """ - nelem = TMatrices.shape[-1] - if ((nelem != TMatrices.shape[-3]) or (2 != TMatrices.shape[-2]) or (2 != TMatrices.shape[-4])): - raise ValueError('The T-matrices must be of shape (N, 2, nelem, 2, nelem) but are of shape %s' % (str(TMatrices.shape),)) - lMax = nelem2lMax(nelem) - if not lMax: - raise ValueError('The "nelem" dimension of T-matrix has invalid value (%d).' % nelem) - # TODO perhaps more checks. - raise Error('Not implemented.') - pass - - - - -import warnings -#@ujit -def scatter_plane_wave_rectarray(omega, epsilon_b, xN, yN, xd, yd, TMatrices, k_dirs, E_0s, - return_pq_0 = False, return_pq= False, return_xy = False, watch_time = False): - """ - Solves the plane wave linear scattering problem for a rectangular array of particles - for one frequency and arbitrary number K of incoming plane waves. - - Parameters - ---------- - omega : positive number - The frequency of the field. - epsilon_b : complex number - Permittivity of the background medium (which has to be isotropic). - xN, yN : positive integers - Particle numbers in the x and y dimensions - xd, yd : positive numbers - Periodicities in the x and y direction - TMatrices : (xN, yN,2,nelem,2,nelem) or compatible or (2,nelem,2,nelem) - The T-matrices in the "Taylor convention" describing the scattering on a single nanoparticle. - If all the particles are identical and equally oriented, only one T-matrix can be given. - nelems = (lMax + 2) * lMax, where lMax is the highest multipole order to which the scattering - is calculated. - Electric wave index is 0, magnetic wave index is 1. - k_dirs : (K,3)-shaped real array or compatible - The direction of the incident field wave vector, normalized to one. - E_0s : (K,3)-shaped complex array or compatible - The electric intensity amplitude of the incident field. - return_pq_0 : bool - Return also the multipole decomposition coefficients of the incoming plane wave. - return_pq : bool NOT IMPLEMENTED - Return also the multipole decomposition coefficients of the field incoming to each - particle (inc. the field scattered from other particles. - return_xy : bool - Return also the cartesian x, y positions of the particles. - watch_time : bool - Inform about the progress on stderr - - Returns - ------- - ab : (K, xN, yN, 2, nelem)-shaped complex array - The a (electric wave), b (magnetic wave) coefficients of the outgoing field for each particle. - If none of return_pq or return_xy is set, the array is not enclosed in a tuple. - pq_0 : (K, xN, yn, 2, nelem)-shaped complex array - The p_0 (electric wave), b_0 (magnetic wave) coefficients of the incoming plane wave for each particle. - pq : (K, xN, yN, 2, nelem)-shaped complex array NOT IMPLEMENTED - The p (electric wave), q (magnetic wave) coefficients of the total exciting field - for each particle (including the field scattered from other particles) - x, y : (xN, yN)-shaped real array - The x,y positions of the nanoparticles. - """ - if (watch_time): - timec = time.time() - print('%.4f: running scatter_plane_wave_rectarray' % timec, file = sys.stderr) - sys.stderr.flush() - nelem = TMatrices.shape[-1] - if ((nelem != TMatrices.shape[-3]) or (2 != TMatrices.shape[-2]) or (2 != TMatrices.shape[-4])): - raise ValueError('The T-matrices must be of shape (N, 2, nelem, 2, nelem) but are of shape %s' % (str(TMatrices.shape),)) - lMax = nelem2lMax(nelem) - if not lMax: - raise ValueError('The "nelem" dimension of T-matrix has invalid value (%d).' % nelem) - if (watch_time): - print('xN = %d, yN = %d, lMax = %d' % (xN, yN, lMax), file = sys.stderr) - sys.stderr.flush() - # TODO perhaps more checks. - k_out = omega * math.sqrt(epsilon_b) / c # wave number - my, ny = get_mn_y(lMax) - N = yN * xN - - J_scat=3 - J_ext=1 - - # Do something with this ugly indexing crap - xind, yind = np.meshgrid(np.arange(xN),np.arange(yN), indexing='ij') - xind = xind.flatten() - yind = yind.flatten() - xyind = np.stack((xind, yind, np.zeros((xind.shape),dtype=int)),axis=-1) - cart_lattice=xyind * np.array([xd, yd, 0]) - x=cart_lattice[:,0] - y=cart_lattice[:,1] - xyind = xyind[:,0:2] - - # Lattice speedup - if (watch_time): - timec = time.time() - print('%.4f: calculating the %d translation matrix elements' % (timec, 8*nelem*nelem*xN*yN), file = sys.stderr) - sys.stderr.flush() - Agrid = np.zeros((nelem, 2*xN, 2*yN, nelem),dtype=np.complex_) - Bgrid = np.zeros((nelem, 2*xN, 2*yN, nelem),dtype=np.complex_) - for yl in range(nelem): # source - for xij in range(2*xN): - for yij in range(2*yN): - for yj in range(nelem): #dest - if((yij != yN) or (xij != xN)): - d_l2j = cart2sph(np.array([(xij-xN)*xd, (yij-yN)*yd, 0])) - Agrid[yj, xij, yij, yl] = Ã(my[yj],ny[yj],my[yl],ny[yl],kdlj=d_l2j[0]*k_out,θlj=d_l2j[1],φlj=d_l2j[2],r_ge_d=False,J=J_scat) - Bgrid[yj, xij, yij, yl] = B̃(my[yj],ny[yj],my[yl],ny[yl],kdlj=d_l2j[0]*k_out,θlj=d_l2j[1],φlj=d_l2j[2],r_ge_d=False,J=J_scat) - - # Translation coefficient matrix T - if (watch_time): - timecold = timec - timec = time.time() - print('%4f: translation matrix elements calculated (elapsed %.2f s), filling the matrix' - % (timec, timec-timecold), file = sys.stderr) - sys.stderr.flush() - transmat = np.zeros((xN* yN, 2, nelem, xN* yN, 2, nelem),dtype=np.complex_) - for l in range(N): - xil, yil = xyind[l] - for j in range(N): - xij, yij = xyind[j] - if (l!=j): - transmat[j,0,:,l,0,:] = Agrid[:, xij - xil + xN, yij - yil + yN, :] - transmat[j,0,:,l,1,:] = Bgrid[:, xij - xil + xN, yij - yil + yN, :] - transmat[j,1,:,l,0,:] = Bgrid[:, xij - xil + xN, yij - yil + yN, :] - transmat[j,1,:,l,1,:] = Agrid[:, xij - xil + xN, yij - yil + yN, :] - Agrid = None - Bgrid = None - if (watch_time): - timecold = timec - timec = time.time() - print('%4f: translation matrix filled (elapsed %.2f s), building the interaction matrix' - % (timec, timec-timecold), file=sys.stderr) - sys.stderr.flush() - - # Now we solve a linear problem (1 - M T) A = M P_0 where M is the T-matrix :-) - MT = np.empty((N,2,nelem,N,2,nelem),dtype=np.complex_) - - TMatrices = np.broadcast_to(TMatrices, (xN, yN, 2, nelem, 2, nelem)) - for j in range(N): # I wonder how this can be done without this loop... - xij, yij = xyind[j] - MT[j] = np.tensordot(TMatrices[xij, yij],transmat[j],axes=([-2,-1],[0,1])) - MT.shape = (N*2*nelem, N*2*nelem) - leftmatrix = np.identity(N*2*nelem) - MT - MT = None - if (watch_time): - timecold = timec - timec = time.time() - print('%.4f: interaction matrix complete (elapsed %.2f s)' % (timec, timec-timecold), - file=sys.stderr) - sys.stderr.flush() - - if ((1 == k_dirs.ndim) and (1 == E_0s.ndim)): - k_cart = k_dirs * k_out # wave vector of the incident plane wave - pq_0 = np.zeros((N,2,nelem), dtype=np.complex_) - p_y0, q_y0 = plane_pq_y(lMax, k_cart, E_0s) - pq_0[:,0,:] = np.exp(1j*np.sum(k_cart[ň,:]*cart_lattice,axis=-1))[:, ň] * p_y0[ň, :] - pq_0[:,1,:] = np.exp(1j*np.sum(k_cart[ň,:]*cart_lattice,axis=-1))[:, ň] * q_y0[ň, :] - if (return_pq_0): - pq_0_arr = pq_0 - MP_0 = np.empty((N,2,nelem),dtype=np.complex_) - #if (watch_time): - # print('%4f: building the interaction matrix' % time.time(), file=sys.stderr) - - for j in range(N): # I wonder how this can be done without this loop... - MP_0[j] = np.tensordot(TMatrices[xij, yij],pq_0[j],axes=([-2,-1],[-2,-1])) - MP_0.shape = (N*2*nelem,) - - if (watch_time): - timecold = time.time() - print('%4f: solving the scattering problem for single incoming wave' % timecold, - file = sys.stderr) - sys.stderr.flush() - ab = np.linalg.solve(leftmatrix, MP_0) - if watch_time: - timec = time.time() - print('%4f: solved (elapsed %.2f s)' % (timec, timec-timecold), file=sys.stderr) - sys.stderr.flush() - - ab.shape = (xN, yN, 2, nelem) - else: - # handle "broadcasting" for k, E - if 1 == k_dirs.ndim: - k_dirs = k_dirs[ň,:] - if 1 == E_0s.ndim: - E_0s = E_0s[ň,:] - K = max(E_0s.shape[-2], k_dirs.shape[-2]) - k_dirs = np.broadcast_to(k_dirs,(K,3)) - E_0s = np.broadcast_to(E_0s, (K,3)) - - # А ну, чики-брики и в дамки! - if watch_time: - timecold = time.time() - print('%.4f: factorizing the interaction matrix' % timecold, file=sys.stderr) - sys.stderr.flush() - lupiv = scipy.linalg.lu_factor(leftmatrix, overwrite_a=True) - leftmatrix = None - if watch_time: - timec = time.time() - print('%.4f: factorization complete (elapsed %.2f s)' % (timec, timec-timecold), - file = sys.stderr) - print('%.4f: solving the scattering problem for %d incoming waves' % (timec, K), - file=sys.stderr) - sys.stderr.flush() - timecold = timec - - if (return_pq_0): - pq_0_arr = np.zeros((K,N,2,nelem), dtype=np.complex_) - ab = np.empty((K,N*2*nelem), dtype=complex) - for ki in range(K): - k_cart = k_dirs[ki] * k_out - pq_0 = np.zeros((N,2,nelem), dtype=np.complex_) - p_y0, q_y0 = plane_pq_y(lMax, k_cart, E_0s[ki]) - pq_0[:,0,:] = np.exp(1j*np.sum(k_cart[ň,:]*cart_lattice,axis=-1))[:, ň] * p_y0[ň, :] - pq_0[:,1,:] = np.exp(1j*np.sum(k_cart[ň,:]*cart_lattice,axis=-1))[:, ň] * q_y0[ň, :] - if (return_pq_0): - pq_0_arr[ki] = pq_0 - MP_0 = np.empty((N,2,nelem),dtype=np.complex_) - for j in range(N): # I wonder how this can be done without this loop... - MP_0[j] = np.tensordot(TMatrices[xij, yij],pq_0[j],axes=([-2,-1],[-2,-1])) - MP_0.shape = (N*2*nelem,) - - ab[ki] = scipy.linalg.lu_solve(lupiv, MP_0) - ab.shape = (K, xN, yN, 2, nelem) - if watch_time: - timec = time.time() - print('%.4f: done (elapsed %.2f s)' % (timec, timec-timecold),file = sys.stderr) - sys.stderr.flush() - if not (return_pq_0 + return_pq + return_xy): - return ab - returnlist = [ab] - if (return_pq_0): - pq_0_arr.shape = ab.shape - returnlist.append(pq_0_arr) - if (return_pq): - warnings.warn("return_pq not implemented, ignoring") - # returnlist.append(pq_arr) - if (return_xy): - returnlist.append(x) - returnlist.append(y) - return tuple(returnlist) - - -import warnings -#@ujit -def scatter_constmultipole_rectarray(omega, epsilon_b, xN, yN, xd, yd, TMatrices, pq_0_c = 1, - return_pq= False, return_xy = False, watch_time = False): - """ - Solves the plane wave linear scattering problem for a rectangular array of particles - for one frequency and constant exciting spherical waves throughout the array. - - Parameters - ---------- - omega : positive number - The frequency of the field. - epsilon_b : complex number - Permittivity of the background medium (which has to be isotropic). - xN, yN : positive integers - Particle numbers in the x and y dimensions - xd, yd : positive numbers - Periodicities in the x and y direction - TMatrices : (xN, yN,2,nelem,2,nelem) or compatible or (2,nelem,2,nelem) - The T-matrices in the "Taylor convention" describing the scattering on a single nanoparticle. - If all the particles are identical and equally oriented, only one T-matrix can be given. - nelems = (lMax + 2) * lMax, where lMax is the highest multipole order to which the scattering - is calculated. - Electric wave index is 0, magnetic wave index is 1. - pq_0_c : (nelem)-shaped complex array or compatible - The initial excitation coefficients for the ("complex") multipole waves, in Taylor's convention. - return_pq : bool NOT IMPLEMENTED - Return also the multipole decomposition coefficients of the field incoming to each - particle (inc. the field scattered from other particles. - return_xy : bool - Return also the cartesian x, y positions of the particles. - watch_time : bool - Inform about the progress on stderr - - Returns - ------- - ab : (nelem, xN, yN, 2, nelem)-shaped complex array - The a (electric wave), b (magnetic wave) coefficients of the outgoing field for each particle. - If none of return_pq or return_xy is set, the array is not enclosed in a tuple. - pq : (nelem, xN, yN, 2, nelem)-shaped complex array NOT IMPLEMENTED - The p (electric wave), q (magnetic wave) coefficients of the total exciting field - for each particle (including the field scattered from other particles) - x, y : (xN, yN)-shaped real array - The x,y positions of the nanoparticles. - """ - if (watch_time): - timec = time.time() - print('%.4f: running scatter_plane_wave_rectarray' % timec, file = sys.stderr) - sys.stderr.flush() - nelem = TMatrices.shape[-1] - if ((nelem != TMatrices.shape[-3]) or (2 != TMatrices.shape[-2]) or (2 != TMatrices.shape[-4])): - raise ValueError('The T-matrices must be of shape (N, 2, nelem, 2, nelem) but are of shape %s' % (str(TMatrices.shape),)) - lMax = nelem2lMax(nelem) - if not lMax: - raise ValueError('The "nelem" dimension of T-matrix has invalid value (%d).' % nelem) - if (watch_time): - print('xN = %d, yN = %d, lMax = %d' % (xN, yN, lMax), file = sys.stderr) - sys.stderr.flush() - # TODO perhaps more checks. - k_out = omega * math.sqrt(epsilon_b) / c # wave number - my, ny = get_mn_y(lMax) - N = yN * xN - - J_scat=3 - J_ext=1 - - # Do something with this ugly indexing crap - xind, yind = np.meshgrid(np.arange(xN),np.arange(yN), indexing='ij') - xind = xind.flatten() - yind = yind.flatten() - xyind = np.stack((xind, yind, np.zeros((xind.shape),dtype=int)),axis=-1) - cart_lattice=xyind * np.array([xd, yd, 0]) - x=cart_lattice[:,0] - y=cart_lattice[:,1] - xyind = xyind[:,0:2] - - # Lattice speedup - if (watch_time): - timec = time.time() - print('%.4f: calculating the %d translation matrix elements' % (timec, 8*nelem*nelem*xN*yN), file = sys.stderr) - sys.stderr.flush() - Agrid = np.zeros((nelem, 2*xN, 2*yN, nelem),dtype=np.complex_) - Bgrid = np.zeros((nelem, 2*xN, 2*yN, nelem),dtype=np.complex_) - for yl in range(nelem): # source - for xij in range(2*xN): - for yij in range(2*yN): - for yj in range(nelem): #dest - if((yij != yN) or (xij != xN)): - d_l2j = cart2sph(np.array([(xij-xN)*xd, (yij-yN)*yd, 0])) - Agrid[yj, xij, yij, yl] = Ã(my[yj],ny[yj],my[yl],ny[yl],kdlj=d_l2j[0]*k_out,θlj=d_l2j[1],φlj=d_l2j[2],r_ge_d=False,J=J_scat) - Bgrid[yj, xij, yij, yl] = B̃(my[yj],ny[yj],my[yl],ny[yl],kdlj=d_l2j[0]*k_out,θlj=d_l2j[1],φlj=d_l2j[2],r_ge_d=False,J=J_scat) - - # Translation coefficient matrix T - if (watch_time): - timecold = timec - timec = time.time() - print('%4f: translation matrix elements calculated (elapsed %.2f s), filling the matrix' - % (timec, timec-timecold), file = sys.stderr) - sys.stderr.flush() - transmat = np.zeros((xN* yN, 2, nelem, xN* yN, 2, nelem),dtype=np.complex_) - for l in range(N): - xil, yil = xyind[l] - for j in range(N): - xij, yij = xyind[j] - if (l!=j): - transmat[j,0,:,l,0,:] = Agrid[:, xij - xil + xN, yij - yil + yN, :] - transmat[j,0,:,l,1,:] = Bgrid[:, xij - xil + xN, yij - yil + yN, :] - transmat[j,1,:,l,0,:] = Bgrid[:, xij - xil + xN, yij - yil + yN, :] - transmat[j,1,:,l,1,:] = Agrid[:, xij - xil + xN, yij - yil + yN, :] - Agrid = None - Bgrid = None - if (watch_time): - timecold = timec - timec = time.time() - print('%4f: translation matrix filled (elapsed %.2f s), building the interaction matrix' - % (timec, timec-timecold), file=sys.stderr) - sys.stderr.flush() - - # Now we solve a linear problem (1 - M T) A = M P_0 where M is the T-matrix :-) - MT = np.empty((N,2,nelem,N,2,nelem),dtype=np.complex_) - - TMatrices = np.broadcast_to(TMatrices, (xN, yN, 2, nelem, 2, nelem)) - for j in range(N): # I wonder how this can be done without this loop... - xij, yij = xyind[j] - MT[j] = np.tensordot(TMatrices[xij, yij],transmat[j],axes=([-2,-1],[0,1])) - MT.shape = (N*2*nelem, N*2*nelem) - leftmatrix = np.identity(N*2*nelem) - MT - MT = None - if (watch_time): - timecold = timec - timec = time.time() - print('%.4f: interaction matrix complete (elapsed %.2f s)' % (timec, timec-timecold), - file=sys.stderr) - sys.stderr.flush() - - # А ну, чики-брики и в дамки! - if watch_time: - timecold = time.time() - print('%.4f: factorizing the interaction matrix' % timecold, file=sys.stderr) - sys.stderr.flush() - lupiv = scipy.linalg.lu_factor(leftmatrix, overwrite_a=True) - leftmatrix = None - if watch_time: - timec = time.time() - print('%.4f: factorization complete (elapsed %.2f s)' % (timec, timec-timecold), - file = sys.stderr) - print('%.4f: solving the scattering problem for %d incoming multipoles' % (timec, nelem*2), - file=sys.stderr) - sys.stderr.flush() - timecold = timec - - if(pq_0_c == 1): - pq_0_c = np.full((2,nelem),1) - ab = np.empty((2,nelem,N*2*nelem), dtype=complex) - for N_or_M in range(2): - for yy in range(nelem): - pq_0 = np.zeros((2,nelem), dtype=np.complex_) - pq_0[N_or_M,yy] = pq_0_c[N_or_M,yy] - pq_0 = np.broadcast_to(pq_0, (N, 2, nelem)) - MP_0 = np.empty((N,2,nelem),dtype=np.complex_) - for j in range(N): # I wonder how this can be done without this loop... - xij, yij = xyind[j] - MP_0[j] = np.tensordot(TMatrices[xij, yij],pq_0[j],axes=([-2,-1],[-2,-1])) - MP_0.shape = (N*2*nelem,) - - ab[N_or_M, yy] = scipy.linalg.lu_solve(lupiv, MP_0) - ab.shape = (2,nelem, xN, yN, 2, nelem) - if watch_time: - timec = time.time() - print('%.4f: done (elapsed %.2f s)' % (timec, timec-timecold),file = sys.stderr) - sys.stderr.flush() - if not (return_pq + return_xy): - return ab - returnlist = [ab] - if (return_pq): - warnings.warn("return_pq not implemented, ignoring") - # returnlist.append(pq_arr) - if (return_xy): - returnlist.append(x) - returnlist.append(y) - return tuple(returnlist) - - - - - - - - - - - - - -# -------------- hexagonal lattice translation coefficients --------------------------- -# Implementation using (C) one-by-one AB-coefficient calculation ufunc -def hexlattice_precalc_AB_save2(file, lMax, k_hexside, maxlayer, circular=True, savepointinfo = False, J_scat=3): - params = { - 'lMax' : lMax, - 'k_hexside' : k_hexside, - 'maxlayer' : maxlayer, - 'circular' : circular, - 'savepointinfo' : savepointinfo, - 'J_scat' : J_scat - } - tpdict = generate_trianglepoints(maxlayer, v3d=True, circular=circular, sixthindices=True, mirrorindices=True) - tphcdict = generate_trianglepoints_hexcomplement(maxlayer, v3d=True, circular=circular, thirdindices=True, mirrorindices=True) - my, ny = get_mn_y(lMax) - nelem = len(my) - a_self_nm = np.empty((tpdict['nmi'].shape[0],nelem,nelem), dtype=complex) - b_self_nm = np.empty((tpdict['nmi'].shape[0],nelem,nelem), dtype=complex) - a_self_m0 = np.empty((tpdict['mi'].shape[1],nelem,nelem), dtype=complex) - b_self_m0 = np.empty((tpdict['mi'].shape[1],nelem,nelem), dtype=complex) - a_d2u_nm = np.empty((tphcdict['nmi'].shape[0],nelem,nelem), dtype=complex) - b_d2u_nm = np.empty((tphcdict['nmi'].shape[0],nelem,nelem), dtype=complex) - a_d2u_m0 = np.empty((tphcdict['mi'].shape[1],nelem,nelem), dtype=complex) - b_d2u_m0 = np.empty((tphcdict['mi'].shape[1],nelem,nelem), dtype=complex) - - k_0 = k_hexside*_s3 # not really a wave vector here because of the normalisation! - tc = trans_calculator(lMax) - - y = np.arange(nelem) - - points = tpdict['points'][tpdict['nmi']] - d_i2j = cart2sph(points) - a_self_nm, b_self_nm = tc.get_AB(my[nx,:,nx],ny[nx,:,nx],my[nx,nx,:],ny[nx,nx,:],k_0*d_i2j[:,nx,nx,0],d_i2j[:,nx,nx,1],d_i2j[:,nx,nx,2],False,J_scat) - - points = tpdict['points'][tpdict['mi'][0]] - d_i2j = cart2sph(points) - a_self_m0, b_self_m0 = tc.get_AB(my[nx,:,nx],ny[nx,:,nx],my[nx,nx,:],ny[nx,nx,:],k_0*d_i2j[:,nx,nx,0],d_i2j[:,nx,nx,1],d_i2j[:,nx,nx,2],False,J_scat) - - points = tphcdict['points'][tphcdict['nmi']] - d_i2j = cart2sph(points) - a_d2u_nm, b_d2u_nm = tc.get_AB(my[nx,:,nx],ny[nx,:,nx],my[nx,nx,:],ny[nx,nx,:],k_0*d_i2j[:,nx,nx,0],d_i2j[:,nx,nx,1],d_i2j[:,nx,nx,2],False,J_scat) - - points = tphcdict['points'][tphcdict['mi'][0]] - d_i2j = cart2sph(points) - a_d2u_m0, b_d2u_m0 = tc.get_AB(my[nx,:,nx],ny[nx,:,nx],my[nx,nx,:],ny[nx,nx,:],k_0*d_i2j[:,nx,nx,0],d_i2j[:,nx,nx,1],d_i2j[:,nx,nx,2],False,J_scat) - - tosave = { - 'a_self_nm' : a_self_nm, - 'a_self_m0' : a_self_m0, - 'b_self_nm' : b_self_nm, - 'b_self_m0' : b_self_m0, - 'a_d2u_nm' : a_d2u_nm, - 'a_d2u_m0' : a_d2u_m0, - 'b_d2u_nm' : b_d2u_nm, - 'b_d2u_m0' : b_d2u_m0, - 'precalc_params' : params - } - if savepointinfo: - tosave['tp_points'] = tpdict['points'], - tosave['tp_si'] = tpdict['si'], - tosave['tp_mi'] = tpdict['mi'], - tosave['tp_nmi'] = tpdict['nmi'] - tosave['tphc_points'] = tphcdict['points'], - tosave['tphc_ti'] = tphcdict['ti'], - tosave['tphc_mi'] = tphcdict['mi'], - tosave['tphc_nmi'] = tphcdict['nmi'] - np.savez(file, **tosave) - - - -# The oldest implementation, using the super-inefficient pure python translation coefficients -def hexlattice_precalc_AB_save_purepy(file, lMax, k_hexside, maxlayer, circular=True, savepointinfo = False, J_scat=3): - params = { - 'lMax' : lMax, - 'k_hexside' : k_hexside, - 'maxlayer' : maxlayer, - 'circular' : circular, - 'savepointinfo' : savepointinfo, - 'J_scat' : J_scat - } - tpdict = generate_trianglepoints(maxlayer, v3d=True, circular=circular, sixthindices=True, mirrorindices=True) - tphcdict = generate_trianglepoints_hexcomplement(maxlayer, v3d=True, circular=circular, thirdindices=True, mirrorindices=True) - my, ny = get_mn_y(lMax) - nelem = len(my) - a_self_nm = np.empty((tpdict['nmi'].shape[0],nelem,nelem), dtype=complex) - b_self_nm = np.empty((tpdict['nmi'].shape[0],nelem,nelem), dtype=complex) - a_self_m0 = np.empty((tpdict['mi'].shape[1],nelem,nelem), dtype=complex) - b_self_m0 = np.empty((tpdict['mi'].shape[1],nelem,nelem), dtype=complex) - a_d2u_nm = np.empty((tphcdict['nmi'].shape[0],nelem,nelem), dtype=complex) - b_d2u_nm = np.empty((tphcdict['nmi'].shape[0],nelem,nelem), dtype=complex) - a_d2u_m0 = np.empty((tphcdict['mi'].shape[1],nelem,nelem), dtype=complex) - b_d2u_m0 = np.empty((tphcdict['mi'].shape[1],nelem,nelem), dtype=complex) - - k_0 = k_hexside*_s3 # not really a wave vector here because of the normalisation! - - points = tpdict['points'][tpdict['nmi']] - for j in range(points.shape[0]): - d_i2j = cart2sph(points[j]) - for yi in range(nelem): - for yj in range(nelem): - a_self_nm[j, yj, yi] = Ã(my[yj],ny[yj],my[yi],ny[yi],kdlj=d_i2j[0]*k_0,θlj=d_i2j[1],φlj=d_i2j[2],r_ge_d=False,J=J_scat) - b_self_nm[j, yj, yi] = B̃(my[yj],ny[yj],my[yi],ny[yi],kdlj=d_i2j[0]*k_0,θlj=d_i2j[1],φlj=d_i2j[2],r_ge_d=False,J=J_scat) - points = tpdict['points'][tpdict['mi'][0]] - for j in range(points.shape[0]): - d_i2j = cart2sph(points[j]) - for yi in range(nelem): - for yj in range(nelem): - a_self_m0[j, yj, yi] = Ã(my[yj],ny[yj],my[yi],ny[yi],kdlj=d_i2j[0]*k_0,θlj=d_i2j[1],φlj=d_i2j[2],r_ge_d=False,J=J_scat) - b_self_m0[j, yj, yi] = B̃(my[yj],ny[yj],my[yi],ny[yi],kdlj=d_i2j[0]*k_0,θlj=d_i2j[1],φlj=d_i2j[2],r_ge_d=False,J=J_scat) - points = tphcdict['points'][tphcdict['nmi']] - for j in range(points.shape[0]): - d_i2j = cart2sph(points[j]) - for yi in range(nelem): - for yj in range(nelem): - a_d2u_nm[j, yj, yi] = Ã(my[yj],ny[yj],my[yi],ny[yi],kdlj=d_i2j[0]*k_0,θlj=d_i2j[1],φlj=d_i2j[2],r_ge_d=False,J=J_scat) - b_d2u_nm[j, yj, yi] = B̃(my[yj],ny[yj],my[yi],ny[yi],kdlj=d_i2j[0]*k_0,θlj=d_i2j[1],φlj=d_i2j[2],r_ge_d=False,J=J_scat) - points = tphcdict['points'][tphcdict['mi'][0]] - for j in range(points.shape[0]): - d_i2j = cart2sph(points[j]) - for yi in range(nelem): - for yj in range(nelem): - a_d2u_m0[j, yj, yi] = Ã(my[yj],ny[yj],my[yi],ny[yi],kdlj=d_i2j[0]*k_0,θlj=d_i2j[1],φlj=d_i2j[2],r_ge_d=False,J=J_scat) - b_d2u_m0[j, yj, yi] = B̃(my[yj],ny[yj],my[yi],ny[yi],kdlj=d_i2j[0]*k_0,θlj=d_i2j[1],φlj=d_i2j[2],r_ge_d=False,J=J_scat) - tosave = { - 'a_self_nm' : a_self_nm, - 'a_self_m0' : a_self_m0, - 'b_self_nm' : b_self_nm, - 'b_self_m0' : b_self_m0, - 'a_d2u_nm' : a_d2u_nm, - 'a_d2u_m0' : a_d2u_m0, - 'b_d2u_nm' : b_d2u_nm, - 'b_d2u_m0' : b_d2u_m0, - 'precalc_params' : params - } - if savepointinfo: - tosave['tp_points'] = tpdict['points'], - tosave['tp_si'] = tpdict['si'], - tosave['tp_mi'] = tpdict['mi'], - tosave['tp_nmi'] = tpdict['nmi'] - tosave['tphc_points'] = tphcdict['points'], - tosave['tphc_ti'] = tphcdict['ti'], - tosave['tphc_mi'] = tphcdict['mi'], - tosave['tphc_nmi'] = tphcdict['nmi'] - np.savez(file, **tosave) - -def hexlattice_precalc_AB_loadunwrap(file, tpdict = None, tphcdict = None, return_points = False): - npz = np.load(file) - precalc_params = npz['precalc_params'][()] - my, ny = get_mn_y(precalc_params['lMax']) - nelem = len(my) - # this I should have made more universal... - if precalc_params['savepointinfo']: - if not tpdict: - tpdict = { - 'points' : npz['tp_points'], - 'si' : npz['tp_si'], - 'mi' : npz['tp_mi'], - 'nmi' : npz['tp_nmi'], - } - if not tphcdict: - tphcdict = { - 'points' : npz['tphc_points'], - 'ti' : npz['tphc_ti'], - 'mi' : npz['tphc_mi'], - 'nmi' : npz['tphc_nmi'] - } - else: - if not tpdict: - tpdict = generate_trianglepoints(maxlayer = precalc_params['maxlayer'], v3d=True, - circular=precalc_params['circular'], sixthindices=True, mirrorindices=True) - if not tphcdict: - tphcdict = generate_trianglepoints_hexcomplement(maxlayer=precalc_params['maxlayer'], v3d=True, - circular=precalc_params['circular'], thirdindices=True, mirrorindices=True) - - # For some obscure reason, I keep getting trailing single-dimension in the beginning for these arrays - for a in (tpdict['points'], tphcdict['points'], tpdict['si'], tpdict['mi'], - tphcdict['ti'], tphcdict['mi']): - if len(a.shape) > 2: - a.shape = a.shape[1::] - - self_tr = tpdict['points'] - d2u_tr = tphcdict['points'] - if len(self_tr.shape)>2: - self_tr = np.reshape(self_tr, self_tr.shape[1::]) - if len(d2u_tr.shape)>2: - d2u_tr = np.reshape(d2u_tr, d2u_tr.shape[1::]) - u2d_tr = -d2u_tr - a_self = np.empty((self_tr.shape[0],nelem,nelem), dtype=complex) - b_self = np.empty((self_tr.shape[0],nelem,nelem), dtype=complex) - a_d2u = np.empty(( d2u_tr.shape[0],nelem,nelem), dtype=complex) - b_d2u = np.empty(( d2u_tr.shape[0],nelem,nelem), dtype=complex) - a_self[tpdict['nmi']]=npz['a_self_nm'] - a_self[tpdict['mi'][0]]=npz['a_self_m0'] - b_self[tpdict['nmi']]=npz['b_self_nm'] - b_self[tpdict['mi'][0]]=npz['b_self_m0'] - mirrorangles = cart2sph(self_tr[tpdict['mi'][1]])[:,2] - cart2sph(self_tr[tpdict['mi'][0]])[:,2] - a_self[tpdict['mi'][1],:,:] = a_self[tpdict['mi'][0],:,:] * np.exp(1j*mirrorangles[:,nx,nx]*(my[nx,nx,:]-my[nx,:,nx])) - b_self[tpdict['mi'][1],:,:] = b_self[tpdict['mi'][0],:,:] * np.exp(1j*mirrorangles[:,nx,nx]*(my[nx,nx,:]-my[nx,:,nx])) - for i in range(1,6): - a_self[tpdict['si'][i],:,:] = a_self[tpdict['si'][0],:,:] * np.exp(1j*math.pi/3*i*(my[nx,:]-my[:,nx])) - b_self[tpdict['si'][i],:,:] = b_self[tpdict['si'][0],:,:] * np.exp(1j*math.pi/3*i*(my[nx,:]-my[:,nx])) - a_d2u[tphcdict['nmi']]=npz['a_d2u_nm'] - a_d2u[tphcdict['mi'][0]]=npz['a_d2u_m0'] - b_d2u[tphcdict['nmi']]=npz['b_d2u_nm'] - b_d2u[tphcdict['mi'][0]]=npz['b_d2u_m0'] - mirrorangles = cart2sph(self_tr[tphcdict['mi'][1]])[:,2] - cart2sph(self_tr[tphcdict['mi'][0]])[:,2] - a_d2u[tphcdict['mi'][1],:,:] = a_d2u[tphcdict['mi'][0],:,:] * np.exp(1j*mirrorangles[:,nx,nx]*(my[nx,nx,:]-my[nx,:,nx])) - b_d2u[tphcdict['mi'][1],:,:] = b_d2u[tphcdict['mi'][0],:,:] * np.exp(1j*mirrorangles[:,nx,nx]*(my[nx,nx,:]-my[nx,:,nx])) - for i in (1,-1): - a_d2u[tphcdict['ti'][i],:,:] = a_d2u[tphcdict['ti'][0],:,:] * np.exp(i*2j*math.pi/3*(my[nx,:]-my[:,nx])) - b_d2u[tphcdict['ti'][i],:,:] = b_d2u[tphcdict['ti'][0],:,:] * np.exp(i*2j*math.pi/3*(my[nx,:]-my[:,nx])) - a_u2d = a_d2u * (-1)**(my[nx,:]-my[:,nx]) - b_u2d = b_d2u * (-1)**(my[nx,:]-my[:,nx]) - d = { - 'a_self' : a_self, - 'b_self' : b_self, - 'a_d2u' : a_d2u, - 'b_d2u' : b_d2u, - 'a_u2d' : a_u2d, - 'b_u2d' : b_u2d, - } - for k in precalc_params.keys(): - d[k] = precalc_params[k] - if return_points: - d['d2u_tr'] = tphcdict['points'] - d['u2d_tr'] = -tphcdict['points'] - d['self_tr'] = tpdict['points'] - return d - -def hexlattice_get_AB(lMax, k_hexside, maxlayer, circular=True, return_points = True, J_scat=3): - params = { - 'lMax' : lMax, - 'k_hexside' : k_hexside, - 'maxlayer' : maxlayer, - 'circular' : circular, - 'savepointinfo' : return_points, # should I delete this key? - 'J_scat' : J_scat - } - tpdict = generate_trianglepoints(maxlayer, v3d=True, circular=circular, sixthindices=True, mirrorindices=True) - tphcdict = generate_trianglepoints_hexcomplement(maxlayer, v3d=True, circular=circular, thirdindices=True, mirrorindices=True) - my, ny = get_mn_y(lMax) - nelem = len(my) - a_self_nm = np.empty((tpdict['nmi'].shape[0],nelem,nelem), dtype=complex) - b_self_nm = np.empty((tpdict['nmi'].shape[0],nelem,nelem), dtype=complex) - a_self_m0 = np.empty((tpdict['mi'].shape[1],nelem,nelem), dtype=complex) - b_self_m0 = np.empty((tpdict['mi'].shape[1],nelem,nelem), dtype=complex) - a_d2u_nm = np.empty((tphcdict['nmi'].shape[0],nelem,nelem), dtype=complex) - b_d2u_nm = np.empty((tphcdict['nmi'].shape[0],nelem,nelem), dtype=complex) - a_d2u_m0 = np.empty((tphcdict['mi'].shape[1],nelem,nelem), dtype=complex) - b_d2u_m0 = np.empty((tphcdict['mi'].shape[1],nelem,nelem), dtype=complex) - - k_0 = k_hexside*_s3 # not really a wave vector here because of the normalisation! - tc = trans_calculator(lMax) - - y = np.arange(nelem) - - points = tpdict['points'][tpdict['nmi']] - d_i2j = cart2sph(points) - a_self_nm, b_self_nm = tc.get_AB_arrays(k_0*d_i2j[:,0],d_i2j[:,1],d_i2j[:,2],np.array([False]),J_scat) - - points = tpdict['points'][tpdict['mi'][0]] - d_i2j = cart2sph(points) - a_self_m0, b_self_m0 = tc.get_AB_arrays(k_0*d_i2j[:,0],d_i2j[:,1],d_i2j[:,2],np.array([False]),J_scat) - - points = tphcdict['points'][tphcdict['nmi']] - d_i2j = cart2sph(points) - a_d2u_nm, b_d2u_nm = tc.get_AB_arrays(k_0*d_i2j[:,0],d_i2j[:,1],d_i2j[:,2],np.array([False]),J_scat) - - points = tphcdict['points'][tphcdict['mi'][0]] - d_i2j = cart2sph(points) - a_d2u_m0, b_d2u_m0 = tc.get_AB_arrays(k_0*d_i2j[:,0],d_i2j[:,1],d_i2j[:,2],np.array([False]),J_scat) - ''' - tosave = { - 'a_self_nm' : a_self_nm, - 'a_self_m0' : a_self_m0, - 'b_self_nm' : b_self_nm, - 'b_self_m0' : b_self_m0, - 'a_d2u_nm' : a_d2u_nm, - 'a_d2u_m0' : a_d2u_m0, - 'b_d2u_nm' : b_d2u_nm, - 'b_d2u_m0' : b_d2u_m0, - 'precalc_params' : params - } - if savepointinfo: - tosave['tp_points'] = tpdict['points'], - tosave['tp_si'] = tpdict['si'], - tosave['tp_mi'] = tpdict['mi'], - tosave['tp_nmi'] = tpdict['nmi'] - tosave['tphc_points'] = tphcdict['points'], - tosave['tphc_ti'] = tphcdict['ti'], - tosave['tphc_mi'] = tphcdict['mi'], - tosave['tphc_nmi'] = tphcdict['nmi'] - np.savez(file, **tosave) - ''' - self_tr = tpdict['points'] - d2u_tr = tphcdict['points'] - if len(self_tr.shape)>2: - self_tr = np.reshape(self_tr, self_tr.shape[1::]) - if len(d2u_tr.shape)>2: - d2u_tr = np.reshape(d2u_tr, d2u_tr.shape[1::]) - u2d_tr = -d2u_tr - a_self = np.empty((self_tr.shape[0],nelem,nelem), dtype=complex) - b_self = np.empty((self_tr.shape[0],nelem,nelem), dtype=complex) - a_d2u = np.empty(( d2u_tr.shape[0],nelem,nelem), dtype=complex) - b_d2u = np.empty(( d2u_tr.shape[0],nelem,nelem), dtype=complex) - a_self[tpdict['nmi']]=a_self_nm - a_self[tpdict['mi'][0]]=a_self_m0 - b_self[tpdict['nmi']]=b_self_nm - b_self[tpdict['mi'][0]]=b_self_m0 - mirrorangles = cart2sph(self_tr[tpdict['mi'][1]])[:,2] - cart2sph(self_tr[tpdict['mi'][0]])[:,2] - a_self[tpdict['mi'][1],:,:] = a_self[tpdict['mi'][0],:,:] * np.exp(1j*mirrorangles[:,nx,nx]*(my[nx,nx,:]-my[nx,:,nx])) - b_self[tpdict['mi'][1],:,:] = b_self[tpdict['mi'][0],:,:] * np.exp(1j*mirrorangles[:,nx,nx]*(my[nx,nx,:]-my[nx,:,nx])) - for i in range(1,6): - a_self[tpdict['si'][i],:,:] = a_self[tpdict['si'][0],:,:] * np.exp(1j*math.pi/3*i*(my[nx,:]-my[:,nx])) - b_self[tpdict['si'][i],:,:] = b_self[tpdict['si'][0],:,:] * np.exp(1j*math.pi/3*i*(my[nx,:]-my[:,nx])) - a_d2u[tphcdict['nmi']]=a_d2u_nm - a_d2u[tphcdict['mi'][0]]=a_d2u_m0 - b_d2u[tphcdict['nmi']]=b_d2u_nm - b_d2u[tphcdict['mi'][0]]=b_d2u_m0 - mirrorangles = cart2sph(self_tr[tphcdict['mi'][1]])[:,2] - cart2sph(self_tr[tphcdict['mi'][0]])[:,2] - a_d2u[tphcdict['mi'][1],:,:] = a_d2u[tphcdict['mi'][0],:,:] * np.exp(1j*mirrorangles[:,nx,nx]*(my[nx,nx,:]-my[nx,:,nx])) - b_d2u[tphcdict['mi'][1],:,:] = b_d2u[tphcdict['mi'][0],:,:] * np.exp(1j*mirrorangles[:,nx,nx]*(my[nx,nx,:]-my[nx,:,nx])) - for i in (1,-1): - a_d2u[tphcdict['ti'][i],:,:] = a_d2u[tphcdict['ti'][0],:,:] * np.exp(i*2j*math.pi/3*(my[nx,:]-my[:,nx])) - b_d2u[tphcdict['ti'][i],:,:] = b_d2u[tphcdict['ti'][0],:,:] * np.exp(i*2j*math.pi/3*(my[nx,:]-my[:,nx])) - a_u2d = a_d2u * (-1)**(my[nx,:]-my[:,nx]) - b_u2d = b_d2u * (-1)**(my[nx,:]-my[:,nx]) - d = { - 'a_self' : a_self, - 'b_self' : b_self, - 'a_d2u' : a_d2u, - 'b_d2u' : b_d2u, - 'a_u2d' : a_u2d, - 'b_u2d' : b_u2d, - } - for k in params.keys(): - d[k] = params[k] - if return_points: - d['d2u_tr'] = tphcdict['points'] - d['u2d_tr'] = -tphcdict['points'] - d['self_tr'] = tpdict['points'] - return d - - diff --git a/qpms/přehled.md b/qpms/přehled.md index b6db483..344ced9 100644 --- a/qpms/přehled.md +++ b/qpms/přehled.md @@ -1,7 +1,5 @@ # Kde co je ## Staré věci -hexpoints.py -legacy.py qpms_p.py (až na změny souřadnic???) ## Nové věci @@ -12,31 +10,8 @@ tmatrices.py types.py svwf.c -## Smíšené / v přepisu -scattering.py qpms_c.pyx -## ??? -hexpoints_c.pyx - -# hexpoints.py -Asi hlavně starý kód pro vytváření trojúhelníkových a hexagonálních mřížek -a počítání (a ukládání) interakčních matic - -## funkce -generate_trianglepoints -generate_trianglepoints_hexcomplement -hexlattice_precalc_AB_save -hexlattice_precalc_AB_loadunwrap -hexlattice_get_AB -hexlattice_zsym_getSVD - -# hexpoints_c.pyx -Obsahuje pouze jedinou funkci (která je i v hexpoints.py). -Používá se tohle vůbec někde? -## funkce -hexlattice_zsym_getSVD - # lattices2d.py Nový kód, manipulace s basemi, vytváření mřížek atd. @@ -58,23 +33,6 @@ filledWS filledWS2 change_basis -# legacy.py -Stařičký kód - -## funkce -q_max -a_q -à -B̃ -G_Mie_scat_precalc_cart -G_Mie_scat_cart -scatter_plane_wave -scatter_plane_wave_rectarray -scatter_constmultipole_rectarray -hexlattice_precalc_AB_save2 -hexlattice_precalc_AB_save_purepy -hexlattice_precalc_AB_loadunwrap -hexlattice_get_AB # qpms_p.py ## funkce @@ -111,12 +69,6 @@ G0L_analytical G0T_analytical G0_sum_1_slow -# scattering.py -## třídy -Scattering -LatticeScattering (neimplementováno nic, asi zrovna rozepsáno) -Scattering_2D_zsym - # scripts_common.py ## funkce make_action_sharedlist diff --git a/qpms/qpms_p.py b/qpms/qpms_p.py index 43681a0..f1a0202 100644 --- a/qpms/qpms_p.py +++ b/qpms/qpms_p.py @@ -1,3 +1,7 @@ +""" +This file contains mostly legacy code, but is still kept for reference. Avoid using this. +""" + import numpy as np from .qpms_c import * ň = np.newaxis diff --git a/qpms/scattering.py b/qpms/scattering.py deleted file mode 100644 index 2344253..0000000 --- a/qpms/scattering.py +++ /dev/null @@ -1,386 +0,0 @@ -''' -Object oriented approach for the classical multiple scattering problem. -''' - -__TODO__ = ''' -- Implement per-scatterer lMax - - This means that Scattering.TMatrices either can not be a single array with a fixed - (N, 2, nelem, 2, nelem) shape but rather list of (2, nelem, 2, nelem) with nelem varying - per particle or some of its elements have to be unused. Anyways, there has to be some kind of - list with the lMaxes. -''' - -import numpy as np -nx = np.newaxis -import time -import scipy -import sys -import warnings -import math -from qpms_c import get_mn_y, trans_calculator # TODO be explicit about what is imported -from .qpms_p import cart2sph, nelem2lMax # TODO be explicit about what is imported -from .timetrack import _time_b, _time_e - -class Scattering(object): - ''' - - This is the most general class for a system of scatterers - in a non-lossy homogeneous background - to be solved with the multiple_scattering method. The scatterers, - as long as they comply with the disjoint circumscribed sphere - hypothesis, can each have any position in the 3D space and - any T-matrix. - - Note that this object describes the scattering problem only for - a single given frequency, as the T-matrices and wavelenght - otherwise differ and all the computationally demanding - parts have to be done for each frequency. However, - the object can be recycled for many incident field shapes - at the given frequency. - - Attributes should be perhaps later redefined to be read-only - (or make descriptors for them). - - Args: - positions: (N,3)-shaped real array - TMatrices: (N,2,nelem,2,nelem)-shaped array - k_0 (float): Wave number for the space between scatterers. - - Attributes: - positions: - TMatrices: - k_0 (float): Wave number for the space between scatterers. - lMax (int): Absolute maximum l for all scatterers. Depending on implementation, - lMax can be smaller for some individual scatterers in certain subclasses. - FIXME: here it is still implemented as constant lMax for all sites, see #! - prepared (bool): Keeps information whether the interaction matrix has - already been built and factorized. - - - ''' - - def __init__(self, positions, TMatrices, k_0, lMax = None, verbose=False, J_scat=3): - self.J_scat = J_scat - self.positions = positions.reshape((-1, positions.shape[-1])) - self.interaction_matrix = None - self.N = self.positions.shape[0] - self.k_0 = k_0 - self.lMax = lMax if lMax else nelem2lMax(TMatrices.shape[-1]) - self.tc = trans_calculator(self.lMax) - nelem = self.lMax * (self.lMax + 2) #! - self.nelem = nelem #! - self.prepared = False - self.TMatrices = np.broadcast_to(TMatrices, (self.N,2,nelem,2,nelem)) - if np.isnan(np.min(TMatrices)): - warnings.warn("Some TMatrices contain NaNs. Expect invalid results") - if np.isnan(np.min(positions)): - warnings.warn("positions contain NaNs. Expect invalid results") - if math.isnan(k_0): - warnings.warn("k_0 is NaN. Expect invalid results") - - - - def prepare(self, keep_interaction_matrix = False, verbose=False): - btime = _time_b(verbose) - if not self.prepared: - if self.interaction_matrix is None: - self.build_interaction_matrix(verbose=verbose) - self.lupiv = scipy.linalg.lu_factor(self.interaction_matrix,overwrite_a = not keep_interaction_matrix) - if not keep_interaction_matrix: - self.interaction_matrix = None - self.prepared = True - _time_e(btime, verbose) - - def build_interaction_matrix(self,verbose = False): - btime = _time_b(verbose) - N = self.N - my, ny = get_mn_y(self.lMax) - nelem = len(my) - leftmatrix = np.zeros((N,2,nelem,N,2,nelem), dtype=complex) - sbtime = _time_b(verbose, step = 'Calculating interparticle translation coefficients') - """ - for i in range(N): - for j in range(N): - for yi in range(nelem): - for yj in range(nelem): - if(i != j): - d_i2j = cart2sph(self.positions[j]-self.positions[i]) - a = Ã(my[yj],ny[yj],my[yi],ny[yi],kdlj=d_i2j[0]*self.k_0,θlj=d_i2j[1],φlj=d_i2j[2],r_ge_d=False,J=self.J_scat) - b = B̃(my[yj],ny[yj],my[yi],ny[yi],kdlj=d_i2j[0]*self.k_0,θlj=d_i2j[1],φlj=d_i2j[2],r_ge_d=False,J=self.J_scat) - leftmatrix[j,0,yj,i,0,yi] = a - leftmatrix[j,1,yj,i,1,yi] = a - leftmatrix[j,0,yj,i,1,yi] = b - leftmatrix[j,1,yj,i,0,yi] = b - """ - kdji = cart2sph(self.positions[:,nx,:] - self.positions[nx,:,:]) - kdji[:,:,0] *= self.k_0 - # get_AB array structure: [j,yj,i,yi] - a, b = self.tc.get_AB(my[nx,:,nx,nx],ny[nx,:,nx,nx],my[nx,nx,nx,:],ny[nx,nx,nx,:], - (kdji[:,:,0])[:,nx,:,nx], (kdji[:,:,1])[:,nx,:,nx], (kdji[:,:,2])[:,nx,:,nx], - False,self.J_scat) - mask = np.broadcast_to(np.eye(N,dtype=bool)[:,nx,:,nx],(N,nelem,N,nelem)) - a[mask] = 0 # no self-translations - b[mask] = 0 - leftmatrix[:,0,:,:,0,:] = a - leftmatrix[:,1,:,:,1,:] = a - leftmatrix[:,0,:,:,1,:] = b - leftmatrix[:,1,:,:,0,:] = b - _time_e(sbtime, verbose, step = 'Calculating interparticle translation coefficients') - # at this point, leftmatrix is the translation matrix - n2id = np.identity(2*nelem) - n2id.shape = (2,nelem,2,nelem) - for j in range(N): - leftmatrix[j] = - np.tensordot(self.TMatrices[j],leftmatrix[j],axes=([-2,-1],[0,1])) - # at this point, jth row of leftmatrix is that of -MT - leftmatrix[j,:,:,j,:,:] += n2id - # now we are done, 1-MT - leftmatrix.shape=(N*2*nelem,N*2*nelem) - self.interaction_matrix = leftmatrix - _time_e(btime, verbose) - - def scatter(self, pq_0, verbose = False): - ''' - pq_0 is (N, nelem, 2)-shaped array - ''' - btime = _time_b(verbose) - if math.isnan(np.min(pq_0)): - warnings.warn("The incident field expansion contains NaNs. Expect invalid results.") - self.prepare(verbose=verbose) - pq_0 = np.broadcast_to(pq_0, (self.N,2,self.nelem)) - MP_0 = np.empty((self.N,2,self.nelem),dtype=np.complex_) - for j in range(self.N): - MP_0[j] = np.tensordot(self.TMatrices[j], pq_0[j],axes=([-2,-1],[-2,-1])) - MP_0.shape = (self.N*2*self.nelem,) - solvebtime = _time_b(verbose,step='Solving the linear equation') - ab = scipy.linalg.lu_solve(self.lupiv, MP_0) - if math.isnan(np.min(ab)): - warnings.warn("Got NaN in the scattering result. Damn.") - raise - _time_e(solvebtime, verbose, step='Solving the linear equation') - ab.shape = (self.N,2,self.nelem) - _time_e(btime, verbose) - return ab - - def scatter_constmultipole(self, pq_0_c, verbose = False): - btime = _time_b(verbose) - N = self.N - self.prepare(verbose=verbose) - nelem = self.nelem - if(pq_0_c ==1): - pq_0_c = np.full((2,nelem),1) - ab = np.empty((2,nelem,N*2*nelem), dtype=complex) - for N_or_M in range(2): - for yy in range(nelem): - pq_0 = np.zeros((2,nelem),dtype=np.complex_) - pq_0[N_or_M,yy] = pq_0_c[N_or_M,yy] - pq_0 = np.broadcast_to(pq_0, (N,2,nelem)) - MP_0 = np.empty((N,2,nelem),dtype=np.complex_) - for j in range(N): - MP_0[j] = np.tensordot(self.TMatrices[j], pq_0[j],axes=([-2,-1],[-2,-1])) - MP_0.shape = (N*2*nelem,) - ab[N_or_M,yy] = scipy.linalg.lu_solve(self.lupiv,MP_0) - ab.shape = (2,nelem,N,2,nelem) - _time_e(btime, verbose) - return ab - -class LatticeScattering(Scattering): - def __init__(self, lattice_spec, k_0, zSym = False): - pass - - -""" -class Scattering_2D_lattice_rectcells(Scattering): - def __init__(self, rectcell_dims, rectcell_elem_positions, cellspec, k_0, rectcell_TMatrices = None, TMatrices = None, lMax = None, verbose=False, J_scat=3): - ''' - cellspec: dvojice ve tvaru (seznam_zaplněnosti, seznam_pozic) - ''' - if (rectcell_TMatrices is None) == (TMatrices is None): - raise ValueError('Either rectcell_TMatrices or TMatrices has to be given') - ###self.positions = ZDE JSEM SKONČIL - self.J_scat = J_scat - self.positions = positions - self.interaction_matrix = None - self.N = positions.shape[0] - self.k_0 = k_0 - self.lMax = lMax if lMax else nelem2lMax(TMatrices.shape[-1]) - nelem = lMax * (lMax + 2) #! - self.nelem = nelem #! - self.prepared = False - self.TMatrices = np.broadcast_to(TMatrices, (self.N,2,nelem,2,nelem)) -""" - -class Scattering_2D_zsym(Scattering): - def __init__(self, positions, TMatrices, k_0, lMax = None, verbose=False, J_scat=3): - Scattering.__init__(self, positions, TMatrices, k_0, lMax, verbose, J_scat) - #TODO some checks on TMatrices symmetry - self.TE_yz = np.arange(self.nelem) - self.TM_yz = self.TE_yz - self.my, self.ny = get_mn_y(self.lMax) - self.TE_NMz = (self.my + self.ny) % 2 - self.TM_NMz = 1 - self.TE_NMz - self.tc = trans_calculator(self.lMax) - # TODO možnost zadávat T-matice rovnou ve zhuštěné podobě - TMatrices_TE = TMatrices[...,self.TE_NMz[:,nx],self.TE_yz[:,nx],self.TE_NMz[nx,:],self.TE_yz[nx,:]] - TMatrices_TM = TMatrices[...,self.TM_NMz[:,nx],self.TM_yz[:,nx],self.TM_NMz[nx,:],self.TM_yz[nx,:]] - self.TMatrices_TE = np.broadcast_to(TMatrices_TE, (self.N, self.nelem, self.nelem)) - self.TMatrices_TM = np.broadcast_to(TMatrices_TM, (self.N, self.nelem, self.nelem)) - self.prepared_TE = False - self.prepared_TM = False - self.interaction_matrix_TE = None - self.interaction_matrix_TM= None - - def prepare_partial(self, TE_or_TM, keep_interaction_matrix = False, verbose=False): - ''' - TE is 0, TM is 1. - ''' - btime = _time_b(verbose) - if (TE_or_TM == 0): #TE - if not self.prepared_TE: - if self.interaction_matrix_TE is None: - self.build_interaction_matrix(0, verbose) - sbtime = _time_b(verbose, step = 'Calculating LU decomposition of the interaction matrix, TE part') - self.lupiv_TE = scipy.linalg.lu_factor(self.interaction_matrix_TE, overwrite_a = not keep_interaction_matrix) - _time_e(sbtime, verbose, step = 'Calculating LU decomposition of the interaction matrix, TE part') - if(np.isnan(np.min(self.lupiv_TE[0])) or np.isnan(np.min(self.lupiv_TE[1]))): - warnings.warn("LU factorisation of interaction matrix contains NaNs. Expect invalid results.") - self.prepared_TE = True - if (TE_or_TM == 1): #TM - if not self.prepared_TM: - if self.interaction_matrix_TM is None: - self.build_interaction_matrix(1, verbose) - sbtime = _time_b(verbose, step = 'Calculating LU decomposition of the interaction matrix, TM part') - self.lupiv_TM = scipy.linalg.lu_factor(self.interaction_matrix_TM, overwrite_a = not keep_interaction_matrix) - _time_e(sbtime, verbose, step = 'Calculating LU decomposition of the interaction matrix, TM part') - if(np.isnan(np.min(self.lupiv_TM[0])) or np.isnan(np.min(self.lupiv_TM[1]))): - warnings.warn("LU factorisation of interaction matrix contains NaNs. Expect invalid results.") - self.prepared_TM = True - _time_e(btime, verbose) - - def prepare(self, keep_interaction_matrix = False, verbose=False): - btime = _time_b(verbose) - if not self.prepared: - self.prepare_partial(0, keep_interaction_matrix, verbose) - self.prepare_partial(1, keep_interaction_matrix, verbose) - self.prepared = True - _time_e(btime, verbose) - - def build_interaction_matrix(self,TE_or_TM = None, verbose = False): - #None means both - btime = _time_b(verbose) - N = self.N - my, ny = get_mn_y(self.lMax) - nelem = len(my) - idm = np.identity(nelem) - if (TE_or_TM == 0): - EoMl = (0,) - elif (TE_or_TM == 1): - EoMl = (1,) - elif (TE_or_TM is None): - EoMl = (0,1) - sbtime = _time_b(verbose, step = 'Calculating interparticle translation coefficients') - kdji = cart2sph(self.positions[:,nx,:] - self.positions[nx,:,:], allow2d=True) - kdji[:,:,0] *= self.k_0 - # get_AB array structure: [j,yj,i,yi] - # FIXME I could save some memory by calculating only half of these coefficients - a, b = self.tc.get_AB(my[nx,:,nx,nx],ny[nx,:,nx,nx],my[nx,nx,nx,:],ny[nx,nx,nx,:], - (kdji[:,:,0])[:,nx,:,nx], (kdji[:,:,1])[:,nx,:,nx], (kdji[:,:,2])[:,nx,:,nx], - False,self.J_scat) - mask = np.broadcast_to(np.eye(N,dtype=bool)[:,nx,:,nx],(N,nelem,N,nelem)) - a[mask] = 0 # no self-translations - b[mask] = 0 - if np.isnan(np.min(a)) or np.isnan(np.min(b)): - warnings.warn("Some of the translation coefficients is a NaN. Expect invalid results.") - _time_e(sbtime, verbose, step = 'Calculating interparticle translation coefficients') - for EoM in EoMl: - leftmatrix = np.zeros((N,nelem,N,nelem), dtype=complex) - y = np.arange(nelem) - yi = y[nx,nx,nx,:] - yj = y[nx,:,nx,nx] - mask = np.broadcast_to((((yi - yj) % 2) == 0),(N,nelem,N,nelem)) - leftmatrix[mask] = a[mask] - mask = np.broadcast_to((((yi - yj) % 2) != 0),(N,nelem,N,nelem)) - leftmatrix[mask] = b[mask] - """ # we use to calculate the AB coefficients here - for i in range(N): - for j in range(i): - for yi in range(nelem): - for yj in range(nelem): - d_i2j = cart2sph(self.positions[j]-self.positions[i]) - if ((yi - yj) % 2) == 0: - tr = Ã(my[yj],ny[yj],my[yi],ny[yi],kdlj=d_i2j[0]*self.k_0,θlj=d_i2j[1],φlj=d_i2j[2],r_ge_d=False,J=self.J_scat) - else: - tr = B̃(my[yj],ny[yj],my[yi],ny[yi],kdlj=d_i2j[0]*self.k_0,θlj=d_i2j[1],φlj=d_i2j[2],r_ge_d=False,J=self.J_scat) - leftmatrix[j,yj,i,yi] = tr - leftmatrix[i,yi,j,yj] = tr if (0 == (my[yj]+my[yi]) % 2) else -tr - _time_e(sbtime, verbose, step = 'Calculating interparticle translation coefficients, T%s part' % ('M' if EoM else 'E')) - """ - for j in range(N): - leftmatrix[j] = - np.tensordot(self.TMatrices_TM[j] if EoM else self.TMatrices_TE[j],leftmatrix[j], - axes = ([-1],[0])) - leftmatrix[j,:,j,:] += idm - leftmatrix.shape = (self.N*self.nelem, self.N*self.nelem) - if np.isnan(np.min(leftmatrix)): - warnings.warn("Interaction matrix contains some NaNs. Expect invalid results.") - if EoM == 0: - self.interaction_matrix_TE = leftmatrix - if EoM == 1: - self.interaction_matrix_TM = leftmatrix - a = None - b = None - _time_e(btime, verbose) - - def scatter_partial(self, TE_or_TM, pq_0, verbose = False): - ''' - pq_0 is (N, nelem)-shaped array - ''' - btime = _time_b(verbose) - self.prepare_partial(TE_or_TM, verbose = verbose) - pq_0 = np.broadcast_to(pq_0, (self.N, self.nelem)) - MP_0 = np.empty((self.N,self.nelem),dtype=np.complex_) - for j in range(self.N): - if TE_or_TM: #TM - MP_0[j] = np.tensordot(self.TMatrices_TM[j], pq_0[j], axes=([-1],[-1])) - else: #TE - MP_0[j] = np.tensordot(self.TMatrices_TE[j], pq_0[j], axes=([-1],[-1])) - MP_0.shape = (self.N*self.nelem,) - solvebtime = _time_b(verbose,step='Solving the linear equation') - ab = scipy.linalg.lu_solve(self.lupiv_TM if TE_or_TM else self.lupiv_TE, MP_0) - _time_e(solvebtime, verbose, step='Solving the linear equation') - ab.shape = (self.N, self.nelem) - _time_e(btime,verbose) - return ab - - def scatter(self, pq_0, verbose = False): - ''' - FI7ME - pq_0 is (N, nelem, 2)-shaped array - ''' - btime = _time_b(verbose) - raise Exception('Not yet implemented') - self.prepare(verbose=verbose) - pq_0 = np.broadcast_to(pq_0, (self.N,2,self.nelem)) - MP_0 = np.empty((N,2,nelem),dtype=np.complex_) - for j in range(self.N): - MP_0[j] = np.tensordot(self.TMatrices[j], pq_0[j],axes=([-2,-1],[-2,-1])) - MP_0.shape = (N*2*self.nelem,) - solvebtime = _time_b(verbose,step='Solving the linear equation') - ab = scipy.linalg.lu_solve(self.lupiv, MP_0) - _time_e(solvebtime, verbose, step='Solving the linear equation') - ab.shape = (N,2,nelem) - _time_e(btime, verbose) - return ab - - def forget_matrices(self): - ''' - Free interaction matrices and set the respective flags - (useful when memory is a bottleneck). - ''' - self.interaction_matrix_TE = None - self.interaction_matrix_TM = None - self.lupiv_TE = None - self.lupiv_TM = None - self.prepared_TE = False - self.prepared_TM = False - self.prepared = False - - diff --git a/qpms/timetrack.py b/qpms/timetrack.py index 358c2b8..89364bd 100644 --- a/qpms/timetrack.py +++ b/qpms/timetrack.py @@ -1,3 +1,6 @@ +""" +Legacy code, used only in scripts_common.py. To be removed in future versions. +""" import time import sys diff --git a/qpms/tmatrices.py b/qpms/tmatrices.py index 347d4d8..1c9064c 100644 --- a/qpms/tmatrices.py +++ b/qpms/tmatrices.py @@ -1,3 +1,8 @@ +""" +Mostly legacy code, but still kept for reference. +Might be removed in future versions. +""" + import numpy as np use_moble_quaternion = False try: