diff --git a/cal.py b/cal.py index 89b1d41..e7cec56 100644 --- a/cal.py +++ b/cal.py @@ -1,12 +1,12 @@ import numpy as np -from scipy.stats import multivariate_normal import matplotlib.pyplot as plt -import scipy import scipy.fftpack as sfft +import matplotlib.patches as patches +import matplotlib def deg_2_rad(winkel): - return winkel/180.0 * np.pi + return winkel / 180.0 * np.pi # all units in angstrom @@ -28,7 +28,6 @@ def mono_2_rutile(c_m, a_m): 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 @@ -41,8 +40,7 @@ class Lattice: 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 = 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 @@ -52,15 +50,15 @@ class Lattice: 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 + 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) + X, Y = self._generate_vec(x_len * 2, y_len * 2) self._get_mono(X, Y) self._get_rutile(X, Y) @@ -102,6 +100,7 @@ def test_lattice(): RESOLUTION = 0.1 +CMAP = "Greys" def image_from_pos(pos_x, pos_y): @@ -141,40 +140,175 @@ def gaussian(img): x = np.linspace(-ratio, ratio, img.shape[0]) y = np.linspace(-1, 1, img.shape[1]) X, Y = np.meshgrid(x, y) - sigma = 0.3 - z = (1/(2*np.pi*sigma*sigma) * np.exp(-(X**2/(2*sigma**2) - + Y**2/(2*sigma**2)))) + 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) - fig, [axs,axs2] = plt.subplots(2, 3) - axs[0].imshow(img) + 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) - plt.pause(.1) + 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))) - + 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) + axs2[0].imshow(img, **IMSHOW_WARGS) img = gaussian(img) - axs2[1].imshow(img) - plt.pause(.1) + 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))) + 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()