diff --git a/array_design/array_gen.py b/array_design/array_gen.py new file mode 100644 index 0000000..62300bf --- /dev/null +++ b/array_design/array_gen.py @@ -0,0 +1,228 @@ +import numpy as np +from typing import Tuple, Optional +from .base_array import * + +def create_ula( + elem_num: int, + spacing: float, + position_errors: Optional[np.ndarray] = None, + gain_errors: Optional[np.ndarray] = None, + phase_errors: Optional[np.ndarray] = None +) -> Array1D: + """ + Create a Uniform Linear Array (ULA) along the x-axis with specified element count and spacing. + + Args: + elem_num (int): Number of elements M (must be a positive integer). + spacing (float): Distance between adjacent elements (in meters or wavelength units, must be > 0). + position_errors (Optional[np.ndarray]): Optional perturbation offsets for each element, (M,), (2, M), or (3, M). + gain_errors (Optional[np.ndarray]): Optional gain error factors (amplitude scale), shape (M,). + phase_errors (Optional[np.ndarray]): Optional phase errors (in radians), shape (M,). + + Returns: + Array1D: Initialized ULA object with 1D position array. + + Raises: + ValueError: If elem_num is not a positive integer or spacing is not positive. + + Notes: + - Element positions are on the x-axis only. + - `positions` is a 1D NumPy array of shape (M,), i.e., scalar x-coordinates. + """ + if not isinstance(elem_num, int) or elem_num <= 0: + raise ValueError(f"elem_num must be a positive integer, got {elem_num}.") + if not isinstance(spacing, (float, int)) or spacing <= 0: + raise ValueError(f"spacing must be a positive number, got {spacing}.") + + # Uniform 1D positions: [0, spacing, ..., (M-1)*spacing] + indices = np.arange(elem_num) + positions = spacing * indices # shape: (M,) + + return Array1D( + positions=positions, + type=ArrayType.ULA, + spacing=spacing, + position_errors=position_errors, + gain_errors=gain_errors, + phase_errors=phase_errors + ) + +def create_ura( + array_size: Tuple[int, int], + spacing: float, + position_errors: Optional[np.ndarray] = None, + gain_errors: Optional[np.ndarray] = None, + phase_errors: Optional[np.ndarray] = None +) -> Array2D: + """ + Create a Uniform Rectangular Array (URA) in the xy-plane. + + The array consists of elements placed on a 2D rectangular grid, with uniform + spacing in both x and y directions. The array lies entirely in the z=0 plane. + + Args: + array_size (Tuple[int, int]): Number of elements along (cols, rows), i.e., (x-axis, y-axis). + The total number of elements is M = cols * rows. + spacing (float): Distance between adjacent elements (in meters or wavelength units). + position_errors (Optional[np.ndarray]): Optional position perturbations of shape (2, M) or (3, M). + gain_errors (Optional[np.ndarray]): Optional gain errors of shape (M,). + phase_errors (Optional[np.ndarray]): Optional phase errors (in radians), shape (M,). + + Returns: + Array2D: Initialized uniform rectangular array object. + + Raises: + ValueError: If array dimensions are not positive integers or spacing is not positive. + + Notes: + - Output positions are of shape (2, M), each column is [x; y] of an element. + - Elements are ordered row-first: for each row, elements increase left to right (x). + - The array is centered at the origin only if manually shifted after creation. + """ + num_cols, num_rows = array_size + if not isinstance(num_cols, int) or num_cols <= 0 or not isinstance(num_cols, int) or num_cols <= 0: + raise ValueError(f"array_size must be two positive integers, got ({num_cols}, {num_cols}).") + if not isinstance(spacing, (float, int)) or spacing <= 0: + raise ValueError(f"spacing must be positive, got {spacing}.") + + total_elements = num_cols * num_rows + + # Generate 2D grid coordinates in xy-plane + x_coords = spacing * np.arange(num_cols) # x: fast-changing axis + y_coords = spacing * np.arange(num_rows) # y: slow-changing axis + xv, yv = np.meshgrid(x_coords, y_coords) # shape (cols, rows) + + # Flatten and stack into (2, M) position matrix + positions = np.zeros((2, total_elements), dtype=float) + positions[0, :] = xv.flatten() # x-axis + positions[1, :] = yv.flatten() # y-axis + + return Array2D( + positions=positions, + type=ArrayType.URA, + spacing=spacing, + position_errors=position_errors, + gain_errors=gain_errors, + phase_errors=phase_errors + ) + +def create_uca( + elem_num: int, + radius: float, + position_errors: Optional[np.ndarray] = None, + gain_errors: Optional[np.ndarray] = None, + phase_errors: Optional[np.ndarray] = None +) -> Array2D: + """ + Create a Uniform Circular Array (UCA) centered at the origin in the xy-plane. + + The array elements are evenly spaced along the circumference of a circle, + with the first element at angle 0 (x-axis), and proceeding counterclockwise. + + Args: + elem_num (int): Number of elements M (must be a positive integer). + radius (float): Radius of the circular array (in meters or wavelength units, must be > 0). + position_errors (Optional[np.ndarray]): Optional position perturbations, shape (2, M) or (3, M). + gain_errors (Optional[np.ndarray]): Optional gain error factors (amplitude scale), shape (M,). + phase_errors (Optional[np.ndarray]): Optional phase errors (in radians), shape (M,). + + Returns: + Array2D: Initialized array object with circular geometry. + + Raises: + ValueError: + - If elem_num is not a positive integer. + - If radius is not a positive number. + + Notes: + - The output positions are of shape (2, M), representing [x; y] coordinates. + - The nominal element spacing (along the arc) is returned as the 'spacing' field. + - The array lies entirely in the xy-plane (z = 0 is implied). + """ + if not isinstance(elem_num, int) or elem_num <= 0: + raise ValueError(f"elem_num ({elem_num}) must be a positive integer.") + if not isinstance(radius, (float, int)) or radius <= 0: + raise ValueError(f"radius ({radius}) must be a positive number.") + + # Compute angular positions (uniform spacing) + theta = np.linspace(0, 2 * np.pi, elem_num, endpoint=False) + + # Compute Cartesian coordinates in 2D + x = radius * np.cos(theta) + y = radius * np.sin(theta) + + # Store as 2×N position matrix: positions[0, :] = x, positions[1, :] = y + positions = np.zeros((2, elem_num), dtype=float) + positions[0, :] = x + positions[1, :] = y + + # Compute approximate arc-length spacing between adjacent elements + step = 2 * np.pi / elem_num + spacing = 2 * radius * np.sin(step / 2) # chord length ≈ arc length for small step + + return Array2D( + positions=positions, + type=ArrayType.UCA, + spacing=spacing, + position_errors=position_errors, + gain_errors=gain_errors, + phase_errors=phase_errors + ) + +def create_cyl_array_with_uca( + elem_per_layer: int, + radius: float, + num_layers: int, + layer_spacing: float, + gain_errors: Optional[np.ndarray] = None, + phase_errors: Optional[np.ndarray] = None, + position_errors: Optional[np.ndarray] = None +) -> Array3D: + """ + Create a cylindrical array composed of stacked UCAs along the z-axis. + + Args: + elem_per_layer (int): Number of elements per circular UCA layer (must be > 0). + radius (float): Radius of each UCA layer (must be > 0). + num_layers (int): Number of stacked layers along z-axis (must be > 0). + layer_spacing (float): Distance between adjacent layers along z-axis (must be > 0). + gain_errors (Optional[np.ndarray]): Optional gain errors, shape (M,). + phase_errors (Optional[np.ndarray]): Optional phase errors, shape (M,). + position_errors (Optional[np.ndarray]): Optional position errors, shape (3, M). + + Returns: + Array3D: Array object with cylindrical geometry. + + Raises: + ValueError: If any input is invalid (non-positive integer or float). + """ + if elem_per_layer <= 0 or num_layers <= 0: + raise ValueError("elem_per_layer and num_layers must be positive integers.") + if radius <= 0: + raise ValueError("radius must be a positive number.") + if layer_spacing <= 0: + raise ValueError("layer_spacing must be a positive number.") + + total_elements = elem_per_layer * num_layers + positions = np.zeros((3, total_elements), dtype=float) + + theta = np.linspace(0, 2 * np.pi, elem_per_layer, endpoint=False) + x = radius * np.cos(theta) + y = radius * np.sin(theta) + + for i in range(num_layers): + z = layer_spacing * i + idx_start = i * elem_per_layer + idx_end = (i + 1) * elem_per_layer + positions[0, idx_start:idx_end] = x + positions[1, idx_start:idx_end] = y + positions[2, idx_start:idx_end] = z + + return Array3D( + positions=positions, + type=ArrayType.CYLINDRICAL, + spacing=layer_spacing, + position_errors=position_errors, + gain_errors=gain_errors, + phase_errors=phase_errors + ) diff --git a/array_design/base_array.py b/array_design/base_array.py index 5694f03..86dda9c 100644 --- a/array_design/base_array.py +++ b/array_design/base_array.py @@ -11,6 +11,7 @@ class ArrayType(Enum): NESTED = 'nested' COPRIME = 'coprime' MRA = 'mra' + CYLINDRICAL = "cylindrical" CUSTOM = 'custom' class BaseArray(ABC): @@ -57,12 +58,12 @@ class BaseArray(ABC): ValueError: If input arrays have mismatched shapes or invalid dimensions. """ - # position-related - positions = np.array(positions) - if positions.ndim != 2: - raise ValueError(f"`positions` must be a 2D array, got ndim={positions.ndim}.") + # element positions + positions = np.array(positions, dtype=float) + if positions.ndim == 1: + positions = positions.reshape(1, -1) # reshape to row vector if positions.shape[0] not in (1, 2, 3): - raise ValueError(f"`positions` must be in 1D, 2D, or 3D space. Got dim={positions.shape[0]}.") + raise ValueError(f"positions must be in 1D, 2D, or 3D space. Got dim={positions.shape[0]}.") self.positions = positions self.dim = positions.shape[0] self.count = positions.shape[1] @@ -75,29 +76,31 @@ class BaseArray(ABC): self.indices = indices # errors and perturbations + # TODO: mutual coupling needs to be modeled + # TODO: if position error is provided, the dimension of the array should change w.s.t. the dimension of the error if position_errors is not None: - position_errors = np.array(position_errors) + position_errors = np.array(position_errors, dtype=float) if position_errors.shape != positions.shape: raise ValueError(f"Shape mismatch: position_errors {position_errors.shape} vs positions {positions.shape}") self.position_errors = position_errors else: - self.position_errors = np.zeros_like(positions) + self.position_errors = position_errors # None if gain_errors is not None: - gain_errors = np.array(gain_errors) + gain_errors = np.array(gain_errors, dtype=float) if gain_errors.shape != positions.shape: raise ValueError(f"Shape mismatch: gain_errors {gain_errors.shape} vs positions {positions.shape}") self.gain_errors = gain_errors else: - self.gain_errors = np.zeros_like(positions) + self.gain_errors = gain_errors # None if phase_errors is not None: - phase_errors = np.array(phase_errors) + phase_errors = np.array(phase_errors, dtype=float) if phase_errors.shape != positions.shape: raise ValueError(f"Shape mismatch: phase_errors {phase_errors.shape} vs positions {positions.shape}") self.phase_errors = phase_errors else: - self.phase_errors = np.zeros_like(positions) + self.phase_errors = phase_errors # None # other attributes self.spacing = spacing @@ -223,7 +226,7 @@ class Array1D(BaseArray): ValueError: If dimension of directions is not 1 or 2. """ directions = np.array(directions) - if directions.ndim == 1: + if directions.ndim in (0, 1): directions = np.reshape(directions, (1, -1)) elif directions.ndim == 2: if directions.shape[0] != 2: @@ -251,13 +254,16 @@ class Array1D(BaseArray): DA = 1j * 2 * np.pi / wavelength * np.outer(positions, np.cos(directions)) * A # apply gain/phase errors - A *= np.exp(1j * self.phase_errors.flatten()) - A *= self.gain_errors.flatten() - if DA is not None: - A_phase = np.exp(1j * self.phase_errors.flatten())[:, None] + if self.gain_errors is not None: A_gain = self.gain_errors.flatten()[:, None] - DA *= A_phase - DA *= A_gain + A *= A_gain + if DA is not None: + DA *= A_gain + if self.phase_errors is not None: + A_phase = np.exp(1j * self.phase_errors.flatten())[:, None] + A *= A_phase + if DA is not None: + DA *= A_phase return A, DA @@ -297,31 +303,36 @@ class Array1D(BaseArray): ValueError: If `center` is provided but its dimension does not match `target_normal`. """ - target_normal = np.reshape(np.array(target_normal, dtype=float), (-1, 1)) - dim = target_normal.shape[0] - if dim not in (2, 3): + target_normal = (np.array(target_normal, dtype=float)).flatten() + target_dim = target_normal.shape[0] + if target_dim not in (2, 3): raise ValueError(f"Target normal must be a 2D or 3D vector, got shape {target_normal.shape}.") positions = self.positions - elem_num = positions.shape[1] + elem_num = self.count # Determine rotation center if center is None: center = np.mean(positions, axis=1) center = np.reshape(np.array(center, dtype=float), (-1, 1)) - if center.shape[0] != dim: + if center.shape[0] != positions.shape[0]: raise ValueError( - f"Center and target_normal must have the same dimension. " - f"Got center dim = {center.shape[0]}, target_normal dim = {dim}." + f"Array center and element positions must have the same dimension. " + f"Got center dim = {center.shape[0]}, positions dim = {positions.shape[0]}." ) # Centered positions aligned along y-axis - if dim == 2: + if target_dim == 2: # Arrange elements along Y-axis and center at origin centered_pos = np.vstack((np.zeros((1, elem_num)), positions - center)) # Get normal vector if current_normal is None: current_normal = np.array([1.0, 0.0]) # assume + else: + current_normal = (np.array(current_normal, dtype=float)).flatten() + current_dim = current_normal.shape[0] + if current_dim != target_dim: + raise ValueError(f"Dimension of target normal ({target_dim}) and current normal ({current_dim}) must match.") current_normal = current_normal / np.linalg.norm(current_normal) target_normal = target_normal / np.linalg.norm(target_normal) # Get the angle between current and target normal vector @@ -330,7 +341,13 @@ class Array1D(BaseArray): [np.cos(angle), -np.sin(angle)], [np.sin(angle), np.cos(angle)] ]) - positions_rot = rot_mat @ centered_pos + center[:, np.newaxis] + # expand the dimension of rotation center + center_exp = np.vstack([ + np.zeros((1, 1)), + center.reshape(1, 1) + ]) + # perform rotation and translation + positions_rot = rot_mat @ centered_pos + center_exp return Array2D( positions = positions_rot, type = self.type, @@ -340,8 +357,8 @@ class Array1D(BaseArray): gain_errors = self.gain_errors, phase_errors = self.phase_errors ), rot_mat - elif dim == 3: - # Arrange elements along xOy plane and center at origin + elif target_dim == 3: + # Arrange elements along y-axis and center at origin centered_pos = np.vstack(( np.zeros((1, elem_num)), positions - center, @@ -350,6 +367,11 @@ class Array1D(BaseArray): # Get normal vector if current_normal is None: current_normal = np.array([0, 0, 1]) + else: + current_normal = (np.array(current_normal, dtype=float)).flatten() + current_dim = current_normal.shape[0] + if current_dim != target_dim: + raise ValueError(f"Dimension of target normal ({target_dim}) and current normal ({current_dim}) must match.") current_normal = current_normal / np.linalg.norm(current_normal) target_normal = target_normal / np.linalg.norm(target_normal) # Rotate to a 3D direction using axis-angle method @@ -370,7 +392,15 @@ class Array1D(BaseArray): [-axis[1], axis[0], 0] ]) rot_mat = np.eye(3) + np.sin(angle) * K + (1 - np.cos(angle)) * K @ K - positions_rot = rot_mat @ centered_pos + center[:, np.newaxis] + + # expand the dimension of rotation center + center_exp = np.vstack([ + np.zeros((1, 1)), + center.reshape(1, 1), + np.zeros((1, 1)) + ]) + # perform rotation and translation + positions_rot = rot_mat @ centered_pos + center return Array3D( positions = positions_rot, type = self.type, @@ -403,14 +433,13 @@ class Array2D(BaseArray): Args: wavelength (float): Signal wavelength. directions (np.ndarray): Input directions, must be in one of the following formats: - - shape (D,) or (1, D): D broadside angles (in radians), valid for 1D array - shape (2, D): D pairs of azimuth and elevation angles (in radians), used for projection to broadside Returns: A (np.ndarray): Steering matrix of shape (M, D) Raises: - ValueError: If dimension of directions is not 1 or 2. + ValueError: If dimension of directions is 2. """ directions = np.array(directions) @@ -422,12 +451,12 @@ class Array2D(BaseArray): ) az, el = directions[0, :], directions[1, :] cos_el = np.cos(el) - x_proj = cos_el * np.sin(az) - y_proj = cos_el * np.cos(az) + x_proj = cos_el * np.cos(az) + y_proj = cos_el * np.sin(az) else: raise ValueError( f"Invalid shape for directions ({directions.ndim}). " - "Expected shape (D,) or (1, D) for 1D broadside angles or shape (2, D) for azimuth/elevation." + "Expected shape (2, D) for azimuth/elevation." ) positions = self.positions @@ -438,15 +467,115 @@ class Array2D(BaseArray): A = np.exp(1j * phase) # apply gain/phase errors - A *= np.exp(1j * self.phase_errors.flatten()) - A *= self.gain_errors.flatten() + if self.gain_errors is not None: + A *= self.gain_errors.flatten() + if self.phase_errors is not None: + A *= np.exp(1j * self.phase_errors.flatten()) return A - def rotate_to(self, target_dim: int, rotation_matrix: np.ndarray) -> BaseArray: - assert target_dim == 3 and rotation_matrix.shape == (3, 2) - new_positions = self.element_positions @ rotation_matrix.T - return Array3D(...) + def rotate( + self, + target_normal: np.ndarray, + current_normal: Optional[np.ndarray] = None, + center: Optional[np.ndarray] = None + ) -> BaseArray: + """ + Rotate a 2D array to align its normal vector with a target direction. + + The array is assumed to be initially aligned along xOy plane (in 3D). The method performs a rigid-body rotation to align the + array's normal with the specified `target_normal`. + + The rotation is performed in three steps: + 1. Translate the centroid of the array to the origin. + 2. Apply the rotation based on the angle between current and target normal vectors. + 3. Translate the array back to the specified or original centroid. + + Args: + target_normal (np.ndarray): Target direction to align the array's normal to. + Must be a 3D vector of shape (3,). + current_normal (Optional[np.ndarray]): Original normal direction of the array. + Defaults to [0, 0, 1] in 3D. + center (Optional[np.ndarray]): Rotation center, defaults to the array's centroid. + Must be the same dimension as `target_normal`. + + Returns: + Tuple[BaseArray, np.ndarray]: + - A new rotated array object (Array3D) with updated element positions. + - The rotation matrix used in the transformation (3x3). + + Raises: + ValueError: If `target_normal` is not 3D. + ValueError: If `center` is provided but its dimension does not match `target_normal`. + """ + + target_normal = (np.array(target_normal, dtype=float)).flatten() + target_dim = target_normal.shape[0] + if target_dim != 3: + raise ValueError(f"Target normal must be a 3D vector for 2D array, got shape {target_normal.shape}.") + + positions = self.positions + elem_num = self.count + + # Determine rotation center + if center is None: + center = np.mean(positions, axis=1) + center = np.reshape(np.array(center, dtype=float), (-1, 1)) + if center.shape[0] != positions.shape[0]: + raise ValueError( + f"Array center and element positions must have the same dimension. " + f"Got center dim = {center.shape[0]}, positions dim = {positions.shape[0]}." + ) + + # Arrange elements along xOy plane and center at origin + centered_pos = np.vstack(( + positions - center, + np.zeros((1, elem_num)) + )) + # Get normal vector + if current_normal is None: + current_normal = np.array([0, 0, 1]) + else: + current_normal = (np.array(current_normal, dtype=float)).flatten() + current_dim = current_normal.shape[0] + if current_dim != target_dim: + raise ValueError(f"Dimension of target normal ({target_dim}) and current normal ({current_dim}) must match.") + current_normal = current_normal / np.linalg.norm(current_normal) + target_normal = target_normal / np.linalg.norm(target_normal) + # Rotate to a 3D direction using axis-angle method + axis = np.cross(current_normal, target_normal) + norm_axis = np.linalg.norm(axis) + if norm_axis < 1e-8: + # If vectors are parallel or anti-parallel + if np.dot(current_normal, target_normal) > 0: + rot_mat = np.eye(3) # No rotation + else: + rot_mat = np.diag([-1, -1, 1]) # 180-degree rotation + else: + axis = axis / norm_axis + angle = np.arccos(np.dot(current_normal, target_normal)) + K = np.array([ + [0, -axis[2], axis[1]], + [axis[2], 0, -axis[0]], + [-axis[1], axis[0], 0] + ]) + rot_mat = np.eye(3) + np.sin(angle) * K + (1 - np.cos(angle)) * K @ K + # expand the dimension of rotation center + center_exp = np.vstack([ + center.reshape(2, 1), + np.zeros((1, 1)) + ]) + # perform rotation and translation + positions_rot = rot_mat @ centered_pos + center_exp + return Array3D( + positions = positions_rot, + type = self.type, + spacing = self.spacing, + indices = self.indices, + position_errors = self.position_errors, + gain_errors = self.gain_errors, + phase_errors = self.phase_errors + ), rot_mat class Array3D(BaseArray): """ @@ -457,16 +586,151 @@ class Array3D(BaseArray): super().__init__(*args, **kwargs) assert self.dim == 3, "Array3D requires 3D element positions." - def steering_matrix(self, wavelength: float, direction: np.ndarray) -> np.ndarray: - """ - Compute 3D steering vector. + def steering_matrix( + self, + wavelength: float, + directions: np.ndarray, + ) -> np.ndarray: + """Compute the steering matrix of 2D array. Args: - wavelength (float): Wavelength. - direction (np.ndarray): 3D unit vector. + wavelength (float): Signal wavelength. + directions (np.ndarray): Input directions, must be in one of the following formats: + - shape (2, D): D pairs of azimuth and elevation angles (in radians), used for projection to broadside Returns: - np.ndarray: Complex steering vector. + A (np.ndarray): Steering matrix of shape (M, D) + + Raises: + ValueError: If dimension of directions is 2. """ + directions = np.array(directions) + + if directions.ndim == 2: + if directions.shape[0] != 2: + raise ValueError( + f"If directions is 2D, its shape must be (2, D), " + f"but got shape {directions.shape}." + ) + az, el = directions[0, :], directions[1, :] + cos_el = np.cos(el) + x_proj = cos_el * np.cos(az) + y_proj = cos_el * np.sin(az) + z_proj = np.sin(el) + else: + raise ValueError( + f"Invalid shape for directions ({directions.ndim}). " + "Expected shape (2, D) for azimuth/elevation." + ) + + positions = self.positions + + # Compute phase shifts k = 2 * np.pi / wavelength - return np.exp(1j * k * self.element_positions @ direction) + phase = k * (positions[0, :, None] * x_proj + positions[1, :, None] * y_proj + + positions[2, :, None] * z_proj) + A = np.exp(1j * phase) + + # apply gain/phase errors + if self.gain_errors is not None: + A *= self.gain_errors.flatten() + if self.phase_errors is not None: + A *= np.exp(1j * self.phase_errors.flatten()) + + return A + + def rotate( + self, + target_normal: np.ndarray, + current_normal: Optional[np.ndarray] = None, + center: Optional[np.ndarray] = None + ) -> BaseArray: + """ + Rotate a 3D array to align its normal vector with a target direction. + + The method performs a rigid-body rotation to align the array's normal with the specified `target_normal`. + + The rotation is performed in three steps: + 1. Translate the centroid of the array to the origin. + 2. Apply the rotation based on the angle between current and target normal vectors. + 3. Translate the array back to the specified or original centroid. + + Args: + target_normal (np.ndarray): Target direction to align the array's normal to. + Must be a 3D vector of shape (3,). + current_normal (Optional[np.ndarray]): Original normal direction of the array. + Defaults to [0, 0, 1] in 3D. + center (Optional[np.ndarray]): Rotation center, defaults to the array's centroid. + Must be the same dimension as `target_normal`. + + Returns: + Tuple[BaseArray, np.ndarray]: + - A new rotated array object (Array3D) with updated element positions. + - The rotation matrix used in the transformation (3x3). + + Raises: + ValueError: If `target_normal` is not 3D. + ValueError: If `center` is provided but its dimension does not match element `positions`. + """ + + target_normal = (np.array(target_normal, dtype=float)).flatten() + target_dim = target_normal.shape[0] + if target_dim != 3: + raise ValueError(f"Target normal must be a 3D vector for 2D array, got shape {target_normal.shape}.") + + positions = self.positions + elem_num = self.count + + # Determine rotation center + if center is None: + center = np.mean(positions, axis=1) + center = np.reshape(np.array(center, dtype=float), (-1, 1)) + if center.shape[0] != positions.shape[0]: + raise ValueError( + f"Array center and element positions must have the same dimension. " + f"Got center dim = {center.shape[0]}, positions dim = {positions.shape[0]}." + ) + + # Arrange elements along xOy plane and center at origin + centered_pos = positions - center + # Get normal vector + if current_normal is None: + current_normal = np.array([0, 0, 1]) + else: + current_normal = (np.array(current_normal, dtype=float)).flatten() + current_dim = current_normal.shape[0] + if current_dim != target_dim: + raise ValueError(f"Dimension of target normal ({target_dim}) and current normal ({current_dim}) must match.") + current_normal = current_normal / np.linalg.norm(current_normal) + target_normal = target_normal / np.linalg.norm(target_normal) + # Rotate to a 3D direction using axis-angle method + axis = np.cross(current_normal, target_normal) + norm_axis = np.linalg.norm(axis) + if norm_axis < 1e-8: + # If vectors are parallel or anti-parallel + if np.dot(current_normal, target_normal) > 0: + rot_mat = np.eye(3) # No rotation + else: + rot_mat = np.diag([-1, -1, 1]) # 180-degree rotation + else: + axis = axis / norm_axis + angle = np.arccos(np.dot(current_normal, target_normal)) + K = np.array([ + [0, -axis[2], axis[1]], + [axis[2], 0, -axis[0]], + [-axis[1], axis[0], 0] + ]) + rot_mat = np.eye(3) + np.sin(angle) * K + (1 - np.cos(angle)) * K @ K + # expand the dimension of rotation center + center_exp = center.reshape(3, 1) + # perform rotation and translation + positions_rot = rot_mat @ centered_pos + center_exp + return Array3D( + positions = positions_rot, + type = self.type, + spacing = self.spacing, + indices = self.indices, + position_errors = self.position_errors, + gain_errors = self.gain_errors, + phase_errors = self.phase_errors + ), rot_mat