diff --git a/rust/Cargo.toml b/rust/Cargo.toml new file mode 100644 index 0000000..97daa1e --- /dev/null +++ b/rust/Cargo.toml @@ -0,0 +1,11 @@ +[package] +name = "rust" +version = "0.1.0" +edition = "2021" + +# See more keys and their definitions at https://doc.rust-lang.org/cargo/reference/manifest.html + +[dependencies] +ndarray = "0.15.6" +ndarray-npy = "0.8.1" +plotters = "0.3.4" diff --git a/rust/plot.py b/rust/plot.py new file mode 100644 index 0000000..ee88bbf --- /dev/null +++ b/rust/plot.py @@ -0,0 +1,9 @@ +import numpy as np + +import matplotlib.pyplot as plt + +file = np.load("./rutile.npz") +plt.scatter(file["x"], file["y"]) +file = np.load("./mono.npz") +plt.scatter(file["x"], file["y"]) +plt.show() diff --git a/rust/src/main.rs b/rust/src/main.rs new file mode 100644 index 0000000..1b4cbd1 --- /dev/null +++ b/rust/src/main.rs @@ -0,0 +1,19 @@ +mod vo2; +use vo2::get_mono; + +use crate::vo2::get_rutile; +use ndarray_npy::NpzWriter; +use std::fs::File; +fn main() { + let (mono_x, mono_y) = get_mono(20, 20); + let (rutile_x, rutile_y) = get_rutile(20, 20); + let mut npz = NpzWriter::new(File::create("mono.npz").unwrap()); + npz.add_array("x", &mono_x); + npz.add_array("y", &mono_y); + npz.finish().unwrap(); + + let mut npz = NpzWriter::new(File::create("rutile.npz").unwrap()); + npz.add_array("x", &rutile_x); + npz.add_array("y", &rutile_y); + npz.finish().unwrap(); +} diff --git a/rust/src/vo2.rs b/rust/src/vo2.rs new file mode 100644 index 0000000..743418d --- /dev/null +++ b/rust/src/vo2.rs @@ -0,0 +1,143 @@ +use ndarray::{Array, Array2}; + +const BASE_A_M: f64 = 5.75; +const BASE_B_M: f64 = 4.5; +const BASE_C_M: f64 = 5.38; + +const BASE_C_R: f64 = 2.856; +const BASE_B_R: f64 = 4.554; +const BASE_A_R: f64 = BASE_B_R; + +const ALPHA_M: f64 = 122.64; // degree + +fn deg_2_rad(angle_degree: f64) -> f64 { + return std::f64::consts::PI / 180. * angle_degree; +} + +fn mono_2_rutile(c_m: f64, a_m: f64) -> (f64, f64) { + let a_r = deg_2_rad(ALPHA_M - 90.).cos() * c_m * BASE_C_M; + let c_r = (a_m) * BASE_A_M + deg_2_rad(ALPHA_M - 90.).sin() * c_m * BASE_C_M; + return (a_r, c_r); +} +pub fn get_rutile(len_x: usize, len_y: usize) -> (Array2, Array2) { + let mut rutile_x_pos: Array2 = Array::zeros((len_x * 2, len_y * 2)); + let mut rutile_y_pos: Array2 = Array::zeros((len_x * 2, len_y * 2)); + for x in 0..2 * len_x { + for y in 0..2 * len_y { + let tmp_x: f64 = if y % 4 >= 2 { x as f64 - 1. } else { x as f64 }; + rutile_x_pos[[x, y]] = tmp_x as f64 * BASE_C_R + (y % 4) as f64 * 0.5 * BASE_C_R; + rutile_y_pos[[x, y]] = y as f64 * 0.5 * BASE_A_R; + } + } + return (rutile_x_pos, rutile_y_pos); +} + +pub fn get_mono(len_x: usize, len_y: usize) -> (Array2, Array2) { + const OFFSET_A_M: f64 = 0.25 - 0.23947; + const OFFSET_C_M: f64 = 0.02646; + let (offset_a_r, offset_c_r) = mono_2_rutile(OFFSET_C_M, OFFSET_A_M); + let mut mono_x_pos: Array2 = Array::zeros((len_x * 2, len_y * 2)); + let mut mono_y_pos: Array2 = Array::zeros((len_x * 2, len_y * 2)); + for x in 0..2 * len_x { + for y in 0..2 * len_y { + let tmp_x: i64 = if y % 4 >= 2 { x as i64 - 1 } else { x as i64 }; + let offset_sign = if tmp_x % 2 == 0 { -1. } else { 1. }; + mono_x_pos[[x, y]] = offset_sign * 2. * offset_a_r + + tmp_x as f64 * BASE_C_R + + (y % 4) as f64 * 0.5 * BASE_C_R; + mono_y_pos[[x, y]] = offset_sign * 2. * offset_c_r + 0.5 * y as f64 * BASE_A_R; + } + } + return (mono_x_pos, mono_y_pos); +} /* + def _mono_2_rutile(self, c_m, a_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_ndc_m = 0.02646 + + offset_a_r, offset_c_r = self._mono_2_rutile(offset_c_m, offset_a_m) + + res = 0.05 + offset_a_r = res * int(offset_a_r/res) + offset_c_r = res * int(offset_c_r/res) + + print(offset_a_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 + + @timeit + 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 + + def reci_rutile(self): + num = 20 + #num = 2 + x = np.arange(-num, num + 1) + y = np.arange(-num, num + 1) + X, Y = np.meshgrid(x, y) + return (X * 0.22 + Y * 0.44).flatten(), (X * 0.349).flatten() + + def reci_mono(self): + x, y = self.reci_rutile() + return x + 0.1083, y + 0.1719 + + def reci_mono_2(self): + x, y = self.reci_rutile() + return x - 0.1083, y + 0.1719 + + def reci(self): + cutoff = 5. + x, y = self.reci_rutile() + mask = np.logical_and(np.abs(x) < cutoff, np.abs(y) < cutoff) + p1 = (x[mask], y[mask]) + x, y = self.reci_mono() + mask = np.logical_and(np.abs(x) < cutoff, np.abs(y) < cutoff) + p2 = (x[mask], y[mask]) + x, y = self.reci_mono_2() + mask = np.logical_and(np.abs(x) < cutoff, np.abs(y) < cutoff) + p3 = (x[mask], y[mask]) + return [p1, p2, p3] + + */