Clean Up
This commit is contained in:
parent
3ca11080d7
commit
a6e526264e
GPG Key ID:
2194FC747048A7FD
97
2d_fourie/lattices.py
Normal file
97
2d_fourie/lattices.py
Normal file
@ -0,0 +1,97 @@
|
||||
import numpy as np
|
||||
|
||||
|
||||
def deg_2_rad(winkel):
|
||||
return winkel / 180.0 * np.pi
|
||||
|
||||
|
||||
# all units in angstrom
|
||||
|
||||
|
||||
class Lattice:
|
||||
def __init__(self, x_len, y_len):
|
||||
pass
|
||||
|
||||
def get_from_mask(self, maske):
|
||||
pass
|
||||
|
||||
|
||||
class SCC_Lattice(Lattice):
|
||||
def __init__(self, x_len, y_len):
|
||||
x = np.arange(x_len) * 5
|
||||
y = np.arange(x_len) * 5
|
||||
self.X, self.Y = np.meshgrid(x, y)
|
||||
|
||||
def get_from_mask(self, maske):
|
||||
return self.X, self.Y
|
||||
|
||||
|
||||
class VO2_Lattice(Lattice):
|
||||
base_a_m = 5.75
|
||||
base_b_m = 4.5
|
||||
base_c_m = 5.38
|
||||
|
||||
base_c_r = 2.856
|
||||
base_b_r = 4.554
|
||||
base_a_r = base_b_r
|
||||
|
||||
alpha_m = 122.64 # degree
|
||||
|
||||
def _mono_2_rutile(self, c_m, a_m):
|
||||
a_r = np.cos(deg_2_rad(self.alpha_m - 90)) * c_m * self.base_c_m
|
||||
c_r = (a_m) * self.base_a_m + \
|
||||
np.sin(deg_2_rad(self.alpha_m - 90)) * c_m * self.base_c_m
|
||||
return a_r, c_r
|
||||
|
||||
def _get_rutile(self, X, Y):
|
||||
self.atom_x_rut = X * self.base_c_r + \
|
||||
np.mod(Y, 4) * 0.5 * self.base_c_r
|
||||
self.atom_y_rut = Y * 0.5 * self.base_a_r
|
||||
|
||||
def _get_mono(self, X, Y):
|
||||
offset_a_m = 0.25 - 0.23947
|
||||
offset_c_m = 0.02646
|
||||
|
||||
offset_a_r, offset_c_r = self.mono_2_rutile(offset_c_m, offset_a_m)
|
||||
|
||||
print("A_r: ", offset_a_r, "C_r: ", offset_c_r)
|
||||
|
||||
self.atom_x_mono = offset_a_r + X * \
|
||||
self.base_c_r + np.mod(Y, 4) * 0.5 * self.base_c_r
|
||||
self.atom_x_mono[np.mod(X, 2) == 0] -= 2 * offset_a_r
|
||||
|
||||
self.atom_y_mono = offset_c_r + 0.5 * Y * self.base_a_r
|
||||
self.atom_y_mono[np.mod(X, 2) == 0] -= 2 * offset_c_r
|
||||
|
||||
def _generate_vec(self, x_len: int, y_len: int):
|
||||
x = np.arange(x_len)
|
||||
y = np.arange(y_len)
|
||||
X, Y = np.meshgrid(x, y)
|
||||
X[np.mod(Y, 4) == 3] = X[np.mod(Y, 4) == 3] - 1
|
||||
X[np.mod(Y, 4) == 2] = X[np.mod(Y, 4) == 2] - 1
|
||||
assert np.mod(x.size, 2) == 0
|
||||
assert np.mod(y.size, 2) == 0
|
||||
|
||||
return X, Y
|
||||
|
||||
def __init__(self, x_len: int, y_len: int):
|
||||
X, Y = self._generate_vec(x_len * 2, y_len * 2)
|
||||
self._get_mono(X, Y)
|
||||
self._get_rutile(X, Y)
|
||||
|
||||
def get_from_mask(self, maske: np.array):
|
||||
inplace_pos_x = np.zeros_like(self.atom_x_mono)
|
||||
inplace_pos_y = np.zeros_like(self.atom_x_mono)
|
||||
|
||||
mask = np.empty_like(self.atom_x_mono, dtype=bool)
|
||||
mask[0::2, 0::2] = maske
|
||||
mask[1::2, 0::2] = maske
|
||||
mask[0::2, 1::2] = maske
|
||||
mask[1::2, 1::2] = maske
|
||||
|
||||
inplace_pos_x[mask] = self.atom_x_rut[mask]
|
||||
inplace_pos_y[mask] = self.atom_y_rut[mask]
|
||||
mask = np.invert(mask)
|
||||
inplace_pos_x[mask] = self.atom_x_mono[mask]
|
||||
inplace_pos_y[mask] = self.atom_y_mono[mask]
|
||||
return inplace_pos_x, inplace_pos_y
|
||||
394
2d_fourie/main.py
Normal file
394
2d_fourie/main.py
Normal file
@ -0,0 +1,394 @@
|
||||
from lattices import SCC_Lattice
|
||||
|
||||
|
||||
import numpy as np
|
||||
import matplotlib.pyplot as plt
|
||||
import scipy.fftpack as sfft
|
||||
import matplotlib.patches as patches
|
||||
import matplotlib
|
||||
import scipy
|
||||
import tqdm
|
||||
|
||||
|
||||
class SpinImage:
|
||||
resolution = 0.1
|
||||
|
||||
def __init__(self, x_pos, y_pos):
|
||||
self.length_x = np.max(x_pos) + self.resolution
|
||||
self.length_y = np.max(y_pos) + self.resolution
|
||||
self.img = self.image_from_pos(x_pos, y_pos)
|
||||
|
||||
def image_from_pos(self, pos_x, pos_y):
|
||||
x_ind = np.arange(0, self.length_x, self.resolution) # angstrom
|
||||
y_ind = np.arange(0, self.length_y, self.resolution) # angstrom
|
||||
img = np.zeros((x_ind.size, y_ind.size))
|
||||
xind = np.searchsorted(x_ind, pos_x)
|
||||
yind = np.searchsorted(y_ind, pos_y)
|
||||
img[xind, yind] = 1
|
||||
return img
|
||||
|
||||
def fft(self):
|
||||
Z_fft = sfft.fft2(self.img)
|
||||
Z_shift = sfft.fftshift(Z_fft)
|
||||
fft_freqx = sfft.fftfreq(self.img.shape[0], self.resolution)
|
||||
fft_freqy = sfft.fftfreq(self.img.shape[1], self.resolution)
|
||||
fft_freqx_clean = sfft.fftshift(fft_freqx)
|
||||
fft_freqy_clean = sfft.fftshift(fft_freqy)
|
||||
return fft_freqx_clean, fft_freqy_clean, np.abs(Z_shift) ** 2
|
||||
|
||||
def pad_it_square(self, additional_pad=0):
|
||||
h = self.img.shape[0]
|
||||
w = self.img.shape[1]
|
||||
print(h, w)
|
||||
xx = np.maximum(h, w) + 2 * additional_pad
|
||||
yy = xx
|
||||
self.length_x = xx * self.resolution
|
||||
self.length_y = yy * self.resolution
|
||||
print("Pad to: ", xx, yy)
|
||||
|
||||
a = (xx - h) // 2
|
||||
aa = xx - a - h
|
||||
|
||||
b = (yy - w) // 2
|
||||
bb = yy - b - w
|
||||
|
||||
self.img = np.pad(self.img, pad_width=(
|
||||
(a, aa), (b, bb)), mode="constant")
|
||||
|
||||
def gaussian(self, sigma):
|
||||
x = np.arange(-self.length_x/2,
|
||||
self.length_x/2, self.resolution)
|
||||
y = np.arange(-self.length_y/2,
|
||||
self.length_y/2, self.resolution)
|
||||
X, Y = np.meshgrid(x, y)
|
||||
z = (
|
||||
1 / (2 * np.pi * sigma * sigma)
|
||||
* np.exp(-(X**2 / (2 * sigma**2) + Y**2 / (2 * sigma**2)))
|
||||
)
|
||||
self.img = np.multiply(self.img, z.T)
|
||||
|
||||
def plot(self, ax, scale=None):
|
||||
if scale is None:
|
||||
ax.imshow(self.img)
|
||||
else:
|
||||
quad = np.ones((int(scale/self.resolution),
|
||||
int(scale/self.resolution)))
|
||||
img = scipy.signal.convolve2d(self.img, quad)
|
||||
ax.imshow(img)
|
||||
|
||||
def blur(self, sigma):
|
||||
self.img = scipy.ndimage.gaussian_filter(self.img, sigma)
|
||||
|
||||
|
||||
def plot(freqx, freqy, intens, ax_log=None, ax_lin=None):
|
||||
#point_x, point_y = reci_rutile()
|
||||
# for px, py in zip(point_x, point_y):
|
||||
# rect = rect_at_point(px, py, "r")
|
||||
# ax.add_patch(rect)
|
||||
# ax.text(
|
||||
# px, py, f"{reduce(extract_rect(intens, px, py, freqx, freqy)):2.2}", clip_on=True
|
||||
# )
|
||||
|
||||
#point_x, point_y = reci_mono()
|
||||
# for px, py in zip(point_x, point_y):
|
||||
# rect = rect_at_point(px, py, "b")
|
||||
# ax.add_patch(rect)
|
||||
# ax.text(
|
||||
# px, py, f"{reduce(extract_rect(intens, px, py, freqx, freqy)):2.2}", clip_on=True
|
||||
# )
|
||||
if ax_log:
|
||||
t = ax_log.imshow(
|
||||
intens,
|
||||
extent=(np.min(freqx), np.max(freqx),
|
||||
np.min(freqy), np.max(freqy)),
|
||||
norm=matplotlib.colors.LogNorm(),
|
||||
cmap="viridis"
|
||||
)
|
||||
plt.colorbar(t)
|
||||
if ax_lin:
|
||||
t = ax_lin.imshow(
|
||||
intens,
|
||||
extent=(np.min(freqx), np.max(freqx),
|
||||
np.min(freqy), np.max(freqy)),
|
||||
cmap="viridis"
|
||||
)
|
||||
plt.colorbar(t)
|
||||
|
||||
|
||||
def test_square():
|
||||
lat = SCC_Lattice(300, 300)
|
||||
si = SpinImage(*lat.get_from_mask(None))
|
||||
fig, axs = plt.subplots(2, 2)
|
||||
si.pad_it_square(10)
|
||||
si.plot(axs[0, 0], 2)
|
||||
si.gaussian(300)
|
||||
# si.blur(3)
|
||||
si.plot(axs[0, 1], 2)
|
||||
|
||||
plt.pause(0.1)
|
||||
fx, fy, intens = si.fft()
|
||||
plot(fx, fy, intens, axs[1, 0], axs[1, 1])
|
||||
print("Done")
|
||||
plt.show()
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
test_square()
|
||||
# def test_lattice():
|
||||
# lat = VO2_Lattice(10, 10)
|
||||
# maske = np.zeros((10, 10), dtype=bool)
|
||||
# x, y = lat.get_from_mask(maske)
|
||||
#
|
||||
# plt.scatter(x, y)
|
||||
# maske = np.invert(maske)
|
||||
# x, y = lat.get_from_mask(maske)
|
||||
# plt.scatter(x, y)
|
||||
#
|
||||
# maske[:3, :5] = False
|
||||
# x, y = lat.get_from_mask(maske)
|
||||
# plt.scatter(x, y)
|
||||
# plt.show()
|
||||
#
|
||||
#
|
||||
# self.resolution = 0.1
|
||||
# CMAP = "Greys"
|
||||
#
|
||||
#
|
||||
# def test_img():
|
||||
# lat = VO2_Lattice(10, 10)
|
||||
# maske = np.ones((10, 10), dtype=bool)
|
||||
# x, y = lat.get_from_mask(maske)
|
||||
# img = image_from_pos(x, y)
|
||||
# plt.imshow(img.T, origin="lower", extent=(0, np.max(x), 0, np.max(y)))
|
||||
# plt.scatter(x, y)
|
||||
# plt.show()
|
||||
#
|
||||
#
|
||||
# def gaussian(img):
|
||||
# x = np.arange(-self.resolution * img.shape[0]/2,
|
||||
# self.resolution * img.shape[0]/2, self.resolution)
|
||||
# y = np.arange(-self.resolution * img.shape[1]/2,
|
||||
# self.resolution * img.shape[1]/2, self.resolution)
|
||||
# X, Y = np.meshgrid(x, y)
|
||||
# sigma = self.resolution * img.shape[0] / 10
|
||||
# print("Sigma: ", sigma)
|
||||
# z = (
|
||||
# 1 / (2 * np.pi * sigma * sigma)
|
||||
# * np.exp(-(X**2 / (2 * sigma**2) + Y**2 / (2 * sigma**2)))
|
||||
# )
|
||||
# return np.multiply(img, z.T)
|
||||
#
|
||||
#
|
||||
# def rect_at_point(x, y, color):
|
||||
# length_2 = 0.08
|
||||
# rect = patches.Rectangle(
|
||||
# (x - length_2, y - length_2),
|
||||
# 2 * length_2,
|
||||
# 2 * length_2,
|
||||
# linewidth=1,
|
||||
# edgecolor=color,
|
||||
# facecolor="none",
|
||||
# )
|
||||
# return rect
|
||||
#
|
||||
#
|
||||
# def reci_rutile():
|
||||
# x = np.arange(-2, 3)
|
||||
# y = np.arange(-2, 3)
|
||||
# X, Y = np.meshgrid(x, y)
|
||||
# return (X * 0.22 + Y * 0.44).flatten(), (X * 0.349).flatten()
|
||||
#
|
||||
#
|
||||
# def reci_mono():
|
||||
# x, y = reci_rutile()
|
||||
# return x + 0.1083, y + 0.1719
|
||||
#
|
||||
#
|
||||
# def draw_big_val_rect(img, x, y, x_index, y_index):
|
||||
# length_2 = 0.08
|
||||
# pos_x_lower = x - length_2
|
||||
# pos_x_upper = x + length_2
|
||||
#
|
||||
# pos_y_lower = y - length_2
|
||||
# pos_y_upper = y + length_2
|
||||
# x_lower = np.searchsorted(x_index, pos_x_lower)
|
||||
# x_upper = np.searchsorted(x_index, pos_x_upper)
|
||||
# y_lower = np.searchsorted(y_index, pos_y_lower)
|
||||
# y_upper = np.searchsorted(y_index, pos_y_upper)
|
||||
#
|
||||
# img[y_lower:y_upper, x_lower:x_upper] = 1e4
|
||||
# return img
|
||||
#
|
||||
#
|
||||
# def extract_rect(img, x, y, x_index, y_index):
|
||||
# length_2 = 0.08
|
||||
#
|
||||
# pos_x_lower = x - length_2
|
||||
# pos_x_upper = x + length_2
|
||||
#
|
||||
# pos_y_lower = y - length_2
|
||||
# pos_y_upper = y + length_2
|
||||
#
|
||||
# x_lower = np.searchsorted(x_index, pos_x_lower)
|
||||
# x_upper = np.searchsorted(x_index, pos_x_upper)
|
||||
#
|
||||
# y_lower = np.searchsorted(y_index, pos_y_lower)
|
||||
# y_upper = np.searchsorted(y_index, pos_y_upper)
|
||||
#
|
||||
# # fix different number of spins possible
|
||||
# if x_upper - x_lower < 10:
|
||||
# x_upper += 1
|
||||
# if y_upper - y_lower < 10:
|
||||
# y_upper += 1
|
||||
#
|
||||
# return img[y_lower:y_upper, x_lower:x_upper]
|
||||
#
|
||||
#
|
||||
# def extract_peaks(freqx, freqy, intens):
|
||||
# rutile = []
|
||||
# point_x, point_y = reci_rutile()
|
||||
# for px, py in zip(point_x, point_y):
|
||||
# rutile.append(reduce(extract_rect(intens, px, py, freqx, freqy)))
|
||||
#
|
||||
# mono = []
|
||||
# point_x, point_y = reci_mono()
|
||||
# for px, py in zip(point_x, point_y):
|
||||
# mono.append(reduce(extract_rect(intens, px, py, freqx, freqy)))
|
||||
# return rutile, mono
|
||||
#
|
||||
#
|
||||
# def plot(ax, freqx, freqy, intens):
|
||||
# point_x, point_y = reci_rutile()
|
||||
# for px, py in zip(point_x, point_y):
|
||||
# rect = rect_at_point(px, py, "r")
|
||||
# ax.add_patch(rect)
|
||||
# ax.text(
|
||||
# px, py, f"{reduce(extract_rect(intens, px, py, freqx, freqy)):2.2}", clip_on=True
|
||||
# )
|
||||
#
|
||||
# point_x, point_y = reci_mono()
|
||||
# for px, py in zip(point_x, point_y):
|
||||
# rect = rect_at_point(px, py, "b")
|
||||
# ax.add_patch(rect)
|
||||
# ax.text(
|
||||
# px, py, f"{reduce(extract_rect(intens, px, py, freqx, freqy)):2.2}", clip_on=True
|
||||
# )
|
||||
# ax.imshow(
|
||||
# intens,
|
||||
# extent=(np.min(freqx), np.max(freqx), np.min(freqy), np.max(freqy)),
|
||||
# norm=matplotlib.colors.LogNorm(),
|
||||
# cmap="Greys"
|
||||
# )
|
||||
#
|
||||
#
|
||||
# def test_all():
|
||||
# LEN = 100
|
||||
# SIZE = 60 * LEN + 1
|
||||
# quad = np.ones((3, 3))
|
||||
#
|
||||
# fig, ax = plt.subplots(1, 3)
|
||||
#
|
||||
# lat = VO2_Lattice(LEN, LEN)
|
||||
# maske = np.ones((LEN, LEN), dtype=bool)
|
||||
# x, y = lat.get_from_mask(maske)
|
||||
# img = image_from_pos(x, y)
|
||||
# img = padding(img, SIZE, SIZE)
|
||||
# #img = scipy.signal.convolve2d(img, quad)
|
||||
# img = gaussian(img)
|
||||
# freqx, freqy, intens_rutile = fft(img)
|
||||
#
|
||||
# img = scipy.signal.convolve2d(img, quad)
|
||||
# ax[0].imshow(img)
|
||||
#
|
||||
# maske = np.zeros((LEN, LEN), dtype=bool)
|
||||
# x, y = lat.get_from_mask(maske)
|
||||
# img = image_from_pos(x, y)
|
||||
# img = padding(img, SIZE, SIZE)
|
||||
# img = gaussian(img)
|
||||
# freqx, freqy, intens_mono = fft(img)
|
||||
#
|
||||
# img = scipy.signal.convolve2d(img, quad)
|
||||
# ax[2].imshow(img)
|
||||
#
|
||||
# maske = np.zeros((LEN, LEN), dtype=bool)
|
||||
# ind = np.arange(LEN*LEN)
|
||||
# np.random.shuffle(ind)
|
||||
# ind = np.unravel_index(ind[:int(LEN*LEN/2)], (LEN, LEN))
|
||||
# maske[ind] = True
|
||||
# x, y = lat.get_from_mask(maske)
|
||||
# img = image_from_pos(x, y)
|
||||
# img = padding(img, SIZE, SIZE)
|
||||
# img = gaussian(img)
|
||||
# freqx, freqy, intens_mono = fft(img)
|
||||
#
|
||||
# img = scipy.signal.convolve2d(img, quad)
|
||||
# ax[1].imshow(img)
|
||||
#
|
||||
# print(np.mean(maske))
|
||||
# x, y = lat.get_from_mask(maske)
|
||||
# img = image_from_pos(x, y)
|
||||
# img = padding(img, SIZE, SIZE)
|
||||
# img = gaussian(img)
|
||||
# freqx, freqy, intens_50 = fft(img)
|
||||
#
|
||||
# fig, axs = plt.subplots(1, 3)
|
||||
# plot(axs[0], freqx=freqx, freqy=freqy, intens=intens_rutile)
|
||||
# plot(axs[2], freqx=freqx, freqy=freqy, intens=intens_mono)
|
||||
# plot(axs[1], freqx=freqx, freqy=freqy, intens=intens_50)
|
||||
# axs[0].set_title("Rutile")
|
||||
# axs[2].set_title("Mono")
|
||||
# axs[1].set_title("50/50")
|
||||
#
|
||||
# for ax in axs:
|
||||
# ax.set_xlim(-1.0, 1.0)
|
||||
# ax.set_ylim(-1.0, 1.0)
|
||||
#
|
||||
#
|
||||
# def eval(maske, lat, LEN):
|
||||
# x, y = lat.get_from_mask(maske)
|
||||
# SIZE = 60 * LEN + 1
|
||||
# img = image_from_pos(x, y)
|
||||
# img = padding(img, SIZE, SIZE)
|
||||
# img = gaussian(img)
|
||||
# freqx, freqy, intens = fft(img)
|
||||
# return extract_peaks(freqx, freqy, intens)
|
||||
#
|
||||
#
|
||||
# def reduce(arr):
|
||||
# arr = np.array(arr)
|
||||
# arr = arr.flatten()
|
||||
# return np.sum(arr[np.argpartition(arr, -8)[-8:]])
|
||||
#
|
||||
#
|
||||
# def main():
|
||||
# LEN = 80
|
||||
# lat = VO2_Lattice(LEN, LEN)
|
||||
# maske = np.zeros((LEN, LEN), dtype=bool)
|
||||
# ind = np.arange(LEN*LEN)
|
||||
# np.random.shuffle(ind)
|
||||
# percentage = []
|
||||
# rutile = []
|
||||
# monoclinic = []
|
||||
# counter = 0
|
||||
# for i in tqdm.tqdm(ind):
|
||||
# i_unravel = np.unravel_index(i, (LEN, LEN))
|
||||
# maske[i_unravel] = True
|
||||
# if np.mod(counter, 300) == 0:
|
||||
# rut, mono = eval(maske, lat, LEN)
|
||||
# percentage.append(np.mean(maske))
|
||||
# rutile.append(reduce(rut))
|
||||
# monoclinic.append(reduce(mono))
|
||||
# counter += 1
|
||||
#
|
||||
# print(len(percentage), len(mono), len(rutile))
|
||||
# print(mono)
|
||||
# plt.figure()
|
||||
# plt.scatter(percentage, np.array(monoclinic)/monoclinic[0], label="mono")
|
||||
# plt.scatter(percentage, np.array(rutile)/rutile[0], label="rut")
|
||||
# plt.legend()
|
||||
#
|
||||
#
|
||||
# if __name__ == "__main__":
|
||||
# test_all()
|
||||
# # main()
|
||||
# plt.show()
|
||||
316
cal.py
316
cal.py
@ -1,316 +0,0 @@
|
||||
import numpy as np
|
||||
import matplotlib.pyplot as plt
|
||||
import scipy.fftpack as sfft
|
||||
import matplotlib.patches as patches
|
||||
import matplotlib
|
||||
|
||||
|
||||
def deg_2_rad(winkel):
|
||||
return winkel / 180.0 * np.pi
|
||||
|
||||
|
||||
# all units in angstrom
|
||||
base_a_m = 5.75
|
||||
base_b_m = 4.5
|
||||
base_c_m = 5.38
|
||||
|
||||
base_c_r = 2.856
|
||||
base_b_r = 4.554
|
||||
base_a_r = base_b_r
|
||||
|
||||
alpha_m = 122.64 # degree
|
||||
|
||||
|
||||
def mono_2_rutile(c_m, a_m):
|
||||
a_r = np.cos(deg_2_rad(alpha_m - 90)) * c_m * base_c_m
|
||||
c_r = (a_m) * base_a_m + np.sin(deg_2_rad(alpha_m - 90)) * c_m * base_c_m
|
||||
return a_r, c_r
|
||||
|
||||
|
||||
class Lattice:
|
||||
def _get_rutile(self, X, Y):
|
||||
self.atom_x_rut = X * base_c_r + np.mod(Y, 4) * 0.5 * base_c_r
|
||||
self.atom_y_rut = Y * 0.5 * base_a_r
|
||||
|
||||
def _get_mono(self, X, Y):
|
||||
offset_a_m = 0.25 - 0.23947
|
||||
offset_c_m = 0.02646
|
||||
|
||||
offset_a_r, offset_c_r = mono_2_rutile(offset_c_m, offset_a_m)
|
||||
|
||||
print("A_r: ", offset_a_r, "C_r: ", offset_c_r)
|
||||
|
||||
self.atom_x_mono = offset_a_r + X * base_c_r + np.mod(Y, 4) * 0.5 * base_c_r
|
||||
self.atom_x_mono[np.mod(X, 2) == 0] -= 2 * offset_a_r
|
||||
|
||||
self.atom_y_mono = offset_c_r + 0.5 * Y * base_a_r
|
||||
self.atom_y_mono[np.mod(X, 2) == 0] -= 2 * offset_c_r
|
||||
|
||||
def _generate_vec(self, x_len: int, y_len: int):
|
||||
x = np.arange(x_len)
|
||||
y = np.arange(y_len)
|
||||
X, Y = np.meshgrid(x, y)
|
||||
X[np.mod(Y, 4) == 3] = X[np.mod(Y, 4) == 3] - 1
|
||||
X[np.mod(Y, 4) == 2] = X[np.mod(Y, 4) == 2] - 1
|
||||
assert np.mod(x.size, 2) == 0
|
||||
assert np.mod(y.size, 2) == 0
|
||||
|
||||
return X, Y
|
||||
|
||||
def __init__(self, x_len: int, y_len: int):
|
||||
X, Y = self._generate_vec(x_len * 2, y_len * 2)
|
||||
self._get_mono(X, Y)
|
||||
self._get_rutile(X, Y)
|
||||
|
||||
def get_from_mask(self, maske: np.array, inplace_pos_x=None, inplace_pos_y=None):
|
||||
if inplace_pos_x is None:
|
||||
inplace_pos_x = np.zeros_like(self.atom_x_mono)
|
||||
if inplace_pos_y is None:
|
||||
inplace_pos_y = np.zeros_like(self.atom_x_mono)
|
||||
|
||||
mask = np.empty_like(self.atom_x_mono, dtype=bool)
|
||||
print(mask.shape, maske.shape)
|
||||
mask[0::2, 0::2] = maske
|
||||
mask[1::2, 0::2] = maske
|
||||
mask[0::2, 1::2] = maske
|
||||
mask[1::2, 1::2] = maske
|
||||
|
||||
inplace_pos_x[mask] = self.atom_x_rut[mask]
|
||||
inplace_pos_y[mask] = self.atom_y_rut[mask]
|
||||
mask = np.invert(mask)
|
||||
inplace_pos_x[mask] = self.atom_x_mono[mask]
|
||||
inplace_pos_y[mask] = self.atom_y_mono[mask]
|
||||
return inplace_pos_x, inplace_pos_y
|
||||
|
||||
|
||||
def test_lattice():
|
||||
lat = Lattice(10, 10)
|
||||
maske = np.zeros((10, 10), dtype=bool)
|
||||
x, y = lat.get_from_mask(maske)
|
||||
|
||||
plt.scatter(x, y)
|
||||
maske = np.invert(maske)
|
||||
x, y = lat.get_from_mask(maske)
|
||||
plt.scatter(x, y)
|
||||
|
||||
maske[:3, :5] = False
|
||||
x, y = lat.get_from_mask(maske)
|
||||
plt.scatter(x, y)
|
||||
plt.show()
|
||||
|
||||
|
||||
RESOLUTION = 0.1
|
||||
CMAP = "Greys"
|
||||
|
||||
|
||||
def image_from_pos(pos_x, pos_y):
|
||||
length_x = np.max(pos_x) + RESOLUTION
|
||||
length_y = np.max(pos_y) + RESOLUTION
|
||||
x_ind = np.arange(0, length_x, RESOLUTION) # angstrom
|
||||
y_ind = np.arange(0, length_y, RESOLUTION) # angstrom
|
||||
img = np.zeros((x_ind.size, y_ind.size))
|
||||
xind = np.searchsorted(x_ind, pos_x)
|
||||
yind = np.searchsorted(y_ind, pos_y)
|
||||
img[xind, yind] = 1
|
||||
return img
|
||||
|
||||
|
||||
def test_img():
|
||||
lat = Lattice(10, 10)
|
||||
maske = np.ones((10, 10), dtype=bool)
|
||||
x, y = lat.get_from_mask(maske)
|
||||
img = image_from_pos(x, y)
|
||||
plt.imshow(img.T, origin="lower", extent=(0, np.max(x), 0, np.max(y)))
|
||||
plt.scatter(x, y)
|
||||
plt.show()
|
||||
|
||||
|
||||
def fft(img):
|
||||
Z_fft = sfft.fft2(img)
|
||||
Z_shift = sfft.fftshift(Z_fft)
|
||||
fft_freqx = sfft.fftfreq(img.shape[0], RESOLUTION)
|
||||
fft_freqy = sfft.fftfreq(img.shape[1], RESOLUTION)
|
||||
fft_freqx_clean = sfft.fftshift(fft_freqx)
|
||||
fft_freqy_clean = sfft.fftshift(fft_freqy)
|
||||
return fft_freqx_clean, fft_freqy_clean, np.abs(Z_shift) ** 2
|
||||
|
||||
|
||||
def gaussian(img):
|
||||
ratio = img.shape[0] / img.shape[1]
|
||||
x = np.linspace(-ratio, ratio, img.shape[0])
|
||||
y = np.linspace(-1, 1, img.shape[1])
|
||||
X, Y = np.meshgrid(x, y)
|
||||
sigma = 0.5
|
||||
z = (
|
||||
1
|
||||
/ (2 * np.pi * sigma * sigma)
|
||||
* np.exp(-(X**2 / (2 * sigma**2) + Y**2 / (2 * sigma**2)))
|
||||
)
|
||||
return np.multiply(img, z.T)
|
||||
|
||||
|
||||
def padding(array, xx, yy):
|
||||
"""
|
||||
:param array: numpy array
|
||||
:param xx: desired height
|
||||
:param yy: desirex width
|
||||
:return: padded array
|
||||
"""
|
||||
|
||||
h = array.shape[0]
|
||||
w = array.shape[1]
|
||||
|
||||
a = (xx - h) // 2
|
||||
aa = xx - a - h
|
||||
|
||||
b = (yy - w) // 2
|
||||
bb = yy - b - w
|
||||
|
||||
return np.pad(array, pad_width=((a, aa), (b, bb)), mode="constant")
|
||||
|
||||
|
||||
def rect_at_point(x, y, color):
|
||||
length_2 = 0.08
|
||||
rect = patches.Rectangle(
|
||||
(x - length_2, y - length_2),
|
||||
2 * length_2,
|
||||
2 * length_2,
|
||||
linewidth=1,
|
||||
edgecolor=color,
|
||||
facecolor="none",
|
||||
)
|
||||
return rect
|
||||
|
||||
|
||||
def reci_rutile():
|
||||
x = np.arange(-2, 3)
|
||||
y = np.arange(-2, 3)
|
||||
X, Y = np.meshgrid(x, y)
|
||||
return (X * 0.22 + Y * 0.44).flatten(), (X * 0.349).flatten()
|
||||
|
||||
|
||||
def reci_mono():
|
||||
x, y = reci_rutile()
|
||||
return x + 0.1083, y + 0.1719
|
||||
|
||||
|
||||
def draw_big_val_rect(img, x, y, x_index, y_index):
|
||||
length_2 = 0.08
|
||||
pos_x_lower = x - length_2
|
||||
pos_x_upper = x + length_2
|
||||
|
||||
pos_y_lower = y - length_2
|
||||
pos_y_upper = y + length_2
|
||||
x_lower = np.searchsorted(x_index, pos_x_lower)
|
||||
x_upper = np.searchsorted(x_index, pos_x_upper)
|
||||
y_lower = np.searchsorted(y_index, pos_y_lower)
|
||||
y_upper = np.searchsorted(y_index, pos_y_upper)
|
||||
|
||||
img[y_lower:y_upper, x_lower:x_upper] = 1e4
|
||||
return img
|
||||
|
||||
|
||||
def extract_rect(img, x, y, x_index, y_index):
|
||||
length_2 = 0.08
|
||||
|
||||
pos_x_lower = x - length_2
|
||||
pos_x_upper = x + length_2
|
||||
|
||||
pos_y_lower = y - length_2
|
||||
pos_y_upper = y + length_2
|
||||
|
||||
x_lower = np.searchsorted(x_index, pos_x_lower)
|
||||
x_upper = np.searchsorted(x_index, pos_x_upper)
|
||||
|
||||
y_lower = np.searchsorted(y_index, pos_y_lower)
|
||||
y_upper = np.searchsorted(y_index, pos_y_upper)
|
||||
|
||||
return img[y_lower:y_upper, x_lower:x_upper]
|
||||
|
||||
|
||||
def main():
|
||||
FFT_KWARGS = {"norm": matplotlib.colors.LogNorm(vmin=1), "cmap": "Greys"}
|
||||
IMSHOW_WARGS = {"cmap": "Greys"}
|
||||
SIZE = 601
|
||||
|
||||
lat = Lattice(10, 10)
|
||||
maske = np.ones((10, 10), dtype=bool)
|
||||
x, y = lat.get_from_mask(maske)
|
||||
img = image_from_pos(x, y)
|
||||
|
||||
img = padding(img, SIZE, SIZE)
|
||||
|
||||
fig, [axs, axs2] = plt.subplots(2, 3)
|
||||
axs[0].imshow(img, **IMSHOW_WARGS)
|
||||
img = gaussian(img)
|
||||
axs[1].imshow(img, **IMSHOW_WARGS)
|
||||
plt.pause(0.1)
|
||||
|
||||
freqx, freqy, intens = fft(img)
|
||||
axs[2].imshow(
|
||||
intens,
|
||||
extent=(np.min(freqx), np.max(freqx), np.min(freqy), np.max(freqy)),
|
||||
**FFT_KWARGS,
|
||||
)
|
||||
|
||||
intens_rut = intens
|
||||
maske = np.zeros((10, 10), dtype=bool)
|
||||
x, y = lat.get_from_mask(maske)
|
||||
img = image_from_pos(x, y)
|
||||
img = padding(img, SIZE, SIZE)
|
||||
|
||||
axs2[0].imshow(img, **IMSHOW_WARGS)
|
||||
img = gaussian(img)
|
||||
axs2[1].imshow(img, **IMSHOW_WARGS)
|
||||
plt.pause(0.1)
|
||||
|
||||
freqx, freqy, intens = fft(img)
|
||||
axs2[2].imshow(
|
||||
intens,
|
||||
extent=(np.min(freqx), np.max(freqx), np.min(freqy), np.max(freqy)),
|
||||
**FFT_KWARGS,
|
||||
)
|
||||
|
||||
# Create a Rectangle patch
|
||||
# Add the patch to the Axes
|
||||
point_x, point_y = reci_rutile()
|
||||
for px, py in zip(point_x, point_y):
|
||||
rect = rect_at_point(px, py, "r")
|
||||
axs[2].add_patch(rect)
|
||||
axs[2].text(
|
||||
px, py, f"{np.sum(extract_rect(intens_rut, px, py, freqx, freqy)):0.2}"
|
||||
)
|
||||
rect = rect_at_point(px, py, "r")
|
||||
axs2[2].add_patch(rect)
|
||||
axs2[2].text(
|
||||
px, py, f"{np.sum(extract_rect(intens, px, py, freqx, freqy)):0.2}"
|
||||
)
|
||||
|
||||
point_x, point_y = reci_mono()
|
||||
for px, py in zip(point_x, point_y):
|
||||
# rect = rect_at_point(px, py,"b")
|
||||
# axs[2].add_patch(rect)
|
||||
rect = rect_at_point(px, py, "b")
|
||||
axs2[2].add_patch(rect)
|
||||
axs2[2].text(
|
||||
px, py, f"{np.sum(extract_rect(intens, px, py, freqx, freqy)):0.2}"
|
||||
)
|
||||
|
||||
axs[2].set_xlim(-1.0, 1.0)
|
||||
axs[2].set_ylim(-1.0, 1.0)
|
||||
axs2[2].set_xlim(-1.0, 1.0)
|
||||
axs2[2].set_ylim(-1.0, 1.0)
|
||||
|
||||
plt.figure()
|
||||
diff = intens_rut - intens
|
||||
plt.imshow(
|
||||
diff, extent=(np.min(freqx), np.max(freqx), np.min(freqy), np.max(freqy))
|
||||
)
|
||||
|
||||
plt.xlim(-1.0, 1.0)
|
||||
plt.ylim(-1.0, 1.0)
|
||||
plt.show()
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
main()
|
||||
Loading…
Reference in New Issue
Block a user