diff --git a/.lefthook.yaml b/.lefthook.yaml index 2aeb41b..e37e0ae 100644 --- a/.lefthook.yaml +++ b/.lefthook.yaml @@ -31,18 +31,19 @@ pre-commit: - name: ruff glob: "*.{py,pyi}" - stage_fixed: true group: piped: true jobs: - name: ruff check run: pixi {run} ruff check --fix {staged_files} + stage_fixed: true - name: ruff format run: pixi {run} ruff format {staged_files} + stage_fixed: true - # - name: pyright + # - name: pyrefly # glob: "*.{py,pyi}" - # run: pixi {run} pyright {staged_files} + # run: pixi {run} pyrefly-check {staged_files} # - name: mypy # glob: "*.{py,pyi}" diff --git a/CHANGELOG.md b/CHANGELOG.md index 05b7371..d40eb88 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -5,6 +5,16 @@ All notable changes to this project will be documented in this file. The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.1.0/), and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html). +## [Unreleased] + +### Added + +- Implemented Bolometry Observer functionality + +### Changed + +- Add `pyrefly` package for type checking (still experimental) + ## [0.2.1] - 2025-11-18 ### Added diff --git a/pixi.toml b/pixi.toml index 4c0b016..3386284 100644 --- a/pixi.toml +++ b/pixi.toml @@ -105,7 +105,7 @@ lefthook = "*" ruff = "*" typos = "*" mypy = "*" -basedpyright = "*" +pyrefly = "*" actionlint = "*" shellcheck = "*" validate-pyproject = "*" @@ -118,7 +118,7 @@ lefthook = { cmd = "lefthook", description = "🔗 Run lefthook" } hooks = { cmd = "lefthook install", description = "🔗 Install pre-commit hooks" } pre-commit = { cmd = "lefthook run pre-commit", description = "🔗 Run pre-commit checks" } mypy = { cmd = "mypy", description = "Type check with mypy" } -pyright = { cmd = "basedpyright", description = "Type check with basedpyright" } +pyrefly-check = { cmd = "pyrefly check", description = "Type check with pyrefly" } ruff-check = { cmd = "ruff check --fix", description = "Lint with ruff" } ruff-format = { cmd = "ruff format", description = "Format with ruff" } dprint = { cmd = "dprint fmt", description = "Format with dprint" } diff --git a/pyproject.toml b/pyproject.toml index ccde31c..c197d79 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -164,23 +164,12 @@ max-line-length = 100 # === Type checking settings === # ------------------------------ [tool.mypy] -files = ["src", "tests"] +files = ["src"] warn_unused_configs = true strict = true enable_error_code = ["ignore-without-code", "truthy-bool"] disable_error_code = ["no-any-return"] plugins = ["numpy.typing.mypy_plugin"] -[tool.basedpyright] -include = ["src", "tests"] -pythonPlatform = "All" -typeCheckingMode = "all" -reportMissingTypeStubs = true - -# no raysect/cherab type stubs; pytest fixtures -reportUnknownMemberType = false -reportUnknownVariableType = false -reportUnknownParameterType = false - -# ruff handles this -reportUnusedParameter = false +[tool.pyrefly] +project-includes = ["src"] diff --git a/src/cherab/imas/__init__.py b/src/cherab/imas/__init__.py index 37746bf..5368303 100644 --- a/src/cherab/imas/__init__.py +++ b/src/cherab/imas/__init__.py @@ -17,8 +17,6 @@ # under the Licence. """The Cherab subpackage for the IMAS.""" -from __future__ import annotations - from importlib.metadata import version as _version __version__ = _version("cherab-imas") diff --git a/src/cherab/imas/ids/bolometer/__init__.py b/src/cherab/imas/ids/bolometer/__init__.py new file mode 100644 index 0000000..e00f925 --- /dev/null +++ b/src/cherab/imas/ids/bolometer/__init__.py @@ -0,0 +1,23 @@ +# Copyright 2023 Euratom +# Copyright 2023 United Kingdom Atomic Energy Authority +# Copyright 2023 Centro de Investigaciones Energéticas, Medioambientales y Tecnológicas +# +# Licensed under the EUPL, Version 1.1 or – as soon they will be approved by the +# European Commission - subsequent versions of the EUPL (the "Licence"); +# You may not use this work except in compliance with the Licence. +# You may obtain a copy of the Licence at: +# +# https://joinup.ec.europa.eu/software/page/eupl5 +# +# Unless required by applicable law or agreed to in writing, software distributed +# under the Licence is distributed on an "AS IS" basis, WITHOUT WARRANTIES OR +# CONDITIONS OF ANY KIND, either express or implied. +# +# See the Licence for the specific language governing permissions and limitations +# under the Licence. +"""Subpackage for loading bolometer data from IMAS IDS structures.""" + +from . import utility +from ._camera import BoloCamera, BoloChannel, Geometry, load_cameras + +__all__ = ["load_cameras", "utility", "BoloCamera", "BoloChannel", "Geometry"] diff --git a/src/cherab/imas/ids/bolometer/_camera.py b/src/cherab/imas/ids/bolometer/_camera.py new file mode 100644 index 0000000..926db5c --- /dev/null +++ b/src/cherab/imas/ids/bolometer/_camera.py @@ -0,0 +1,322 @@ +# Copyright 2023 Euratom +# Copyright 2023 United Kingdom Atomic Energy Authority +# Copyright 2023 Centro de Investigaciones Energéticas, Medioambientales y Tecnológicas +# +# Licensed under the EUPL, Version 1.1 or – as soon they will be approved by the +# European Commission - subsequent versions of the EUPL (the "Licence"); +# You may not use this work except in compliance with the Licence. +# You may obtain a copy of the Licence at: +# +# https://joinup.ec.europa.eu/software/page/eupl5 +# +# Unless required by applicable law or agreed to in writing, software distributed +# under the Licence is distributed on an "AS IS" basis, WITHOUT WARRANTIES OR +# CONDITIONS OF ANY KIND, either express or implied. +# +# See the Licence for the specific language governing permissions and limitations +# under the Licence. +"""Module for loading bolometer cameras from IMAS bolometer IDS.""" + +from __future__ import annotations + +from collections.abc import Iterator +from dataclasses import dataclass, field + +import numpy as np +from numpy.typing import NDArray +from raysect.core.math import Point3D, Vector3D, from_cylindrical + +from imas.ids_structure import IDSStructure +from imas.ids_toplevel import IDSToplevel + +from .utility import CameraType, GeometryType + + +@dataclass +class Geometry: + """Represent the geometric specification of a bolometer sensor head or slit aperture. + + The Geometry describes both simple rectangular / circular sensor faces and + arbitrary polygonal / polyline definitions via explicit coordinates. + It encapsulates a local orthonormal (or user supplied) basis, extents, and derived surface + properties. + """ + + centre: Point3D = field(default_factory=Point3D) + """Geometric centre (reference point) of the sensor/slit in global 3D coordinates.""" + type: GeometryType = GeometryType.RECTANGLE + """Enumerated shape type. + + Common values include `RECTANGLE`, `CIRCLE`, `OUTLINE`, etc. + Defaults to `.GeometryType.RECTANGLE`. + """ + basis_x: Vector3D = field(default_factory=lambda: Vector3D(1, 0, 0)) + """Local x-axis direction vector lying in the sensor/slit plane.""" + basis_y: Vector3D = field(default_factory=lambda: Vector3D(0, 1, 0)) + """Local y-axis direction vector lying in the sensor/slit plane.""" + basis_z: Vector3D = field(default_factory=lambda: Vector3D(0, 0, 1)) + """Local outward-facing normal vector perpendicular to the sensor/slit plane. + + This vector must be directed toward the radiation sources. + When None, it can be derived as the cross product of `basis_x` and `basis_y`. + """ + dx: float = 0.0 + """Width along `basis_x` for rectangular geometry. + + None when not applicable (e.g., circular or polygonal types). + """ + dy: float = 0.0 + """Width along `basis_y` for rectangular geometry. + + None when not applicable (e.g., circular or polygonal types). + """ + surface: float = 0.0 + """Precomputed surface area of the aperture face. + + If None, may be derived from `dx`/`dy` (rectangles), `radius` (circles), or `coords` (polygons) + during validation or runtime. + """ + radius: float = 0.0 + """Radius for circular geometry types. + + None if the geometry is not circular. + """ + coords: NDArray[np.float64] | None = None + """Coordinate array defining the outline in the `basis_x` and `basis_y` plane. + + This array is used to represent a complex geometric outline. + The array shape is ``(2, N)`` where ``N`` is the number of points. + """ + + +@dataclass +class BoloChannel: + """Represent a bolometer camera channel. + + Each channel contains one foil and associated several slits. + """ + + foil: Geometry + """Geometry of the foil used in the bolometer channel.""" + slits: list[Geometry] + """List of geometries representing the slits associated with the bolometer channel.""" + num_subcollimator: int + """Number of sub-collimators in the bolometer channel.""" + thickness_subcollimator: float + """Thickness of each sub-collimator in the bolometer channel [m].""" + + +@dataclass +class BoloCamera: + """Represent a bolometer camera and its associated data. + + This class encapsulates all the properties and data channels associated with + a bolometer camera used for plasma diagnostics. + """ + + name: str + """Unique identifier or label for the bolometer camera.""" + description: str + """Detailed description of the bolometer camera, including its purpose, location, or other relevant information.""" + type: CameraType + """Type of Bolometer camera: Pinhole/Collimator""" + channels: list[BoloChannel] + """List of individual bolometer channels belonging to this camera.""" + + def __repre__(self) -> str: + return f"" + + def __len__(self) -> int: + """Return the number of channels in the bolometer camera. + + Returns + ------- + int + Number of channels. + """ + return len(self.channels) + + def __iter__(self) -> Iterator[BoloChannel]: + """Return an iterator over the bolometer camera channels. + + Returns + ------- + `Iterator[BoloChannel]` + An iterator over the bolometer camera channels. + """ + return iter(self.channels) + + def __getitem__(self, index: int) -> BoloChannel: + """Return the bolometer channel at the specified index. + + Parameters + ---------- + index + Index of the channel to retrieve. + + Returns + ------- + `.BoloChannel` + The bolometer channel at the specified index. + """ + return self.channels[index] + + +def load_cameras(ids: IDSToplevel) -> list[BoloCamera]: + """Load bolometer cameras from the bolometer IDS. + + This function retrieves the camera information from the bolometer IDS and organizes it into a + structured format. + + Parameters + ---------- + ids + The bolometer IDS. + + Returns + ------- + `.BoloCamera` + A list of bolometer camera data structures. + + Raises + ------ + ValueError + If the provided IDS is not a bolometer IDS. + RuntimeError + If no cameras are found in the IDS. + """ + if not str(ids.metadata.name) == "bolometer": + raise ValueError(f"Invalid bolometer IDS ({ids.metadata.name}).") + + cameras = getattr(ids, "camera", []) + if not len(cameras): + raise RuntimeError("No camera found in IDS.") + + bolo_data: list[BoloCamera] = [] + + for camera in cameras: + name = str(camera.name) + description = str(camera.description) + camera_type = CameraType.from_value(camera.type.index.value) + + channels: list[BoloChannel] = [] + for channel in camera.channel: + # Sub-collimator + n_subcollimators = int(channel.subcollimators_n.value) + thickness_subcollimator = float(channel.subcollimators_separation.value) + + channels.append( + BoloChannel( + foil=load_geometry(channel.detector), + slits=[load_geometry(aperture) for aperture in channel.aperture], + num_subcollimator=n_subcollimators if n_subcollimators > 1 else 1, + thickness_subcollimator=thickness_subcollimator + if thickness_subcollimator > 0 + else 0.0, + ) + ) + + bolo_data.append( + BoloCamera( + name=name, + description=description, + type=camera_type, + channels=channels, + ) + ) + + return bolo_data + + +def load_geometry(sensor: IDSStructure) -> Geometry: + """Load the geometry of a sensor or aperture from the bolometer IDS. + + Parameters + ---------- + sensor + Detector or aperture structure object. + + Returns + ------- + `.Geometry` + Object with the following attributes: + + - ``'centre'``: Point3D with the coordinates of the sensor centre, + - ``'type'``: Geometry type of the sensor, + - ``'basis_x'``: Vector3D with the x-basis vector of the sensor, + - ``'basis_y'``: Vector3D with the y-basis vector of the sensor, + - ``'basis_z'``: Vector3D with the z-basis vector of the sensor, + - ``'dx'``: Width of the sensor in the x-direction [m], + - ``'dy'``: Width of the sensor in the y-direction [m], + - ``'surface'``: Surface area of the sensor [m²], + - ``'radius'``: Radius of the sensor [m] if the geometry type is circular, + - ``'coords'``: Outline coordinates for basis_x and basis_y if the geometry type is outline. + + Some keys may not be present depending on the geometry type. + + Raises + ------ + ValueError + If the geometry type is invalid. + """ + geometry = Geometry() + + geometry.centre = from_cylindrical( + sensor.centre.r, sensor.centre.z, np.rad2deg(sensor.centre.phi) + ) + + geometry_type = GeometryType.from_value(int(sensor.geometry_type.value)) + geometry.type = geometry_type + + match geometry_type: + case GeometryType.RECTANGLE: + # Unit vectors + geometry.basis_x = Vector3D( + sensor.x2_unit_vector.x, sensor.x2_unit_vector.y, sensor.x2_unit_vector.z + ) + geometry.basis_y = Vector3D( + sensor.x1_unit_vector.x, sensor.x1_unit_vector.y, sensor.x1_unit_vector.z + ) + geometry.basis_z = Vector3D( + sensor.x3_unit_vector.x, sensor.x3_unit_vector.y, sensor.x3_unit_vector.z + ) + + # Dimensions + geometry.dx = float(sensor.x2_width.value) + geometry.dy = float(sensor.x1_width.value) + + case GeometryType.CIRCULAR: + # Radius + radius = getattr(sensor, "radius", None) + if radius is None or radius <= 0: + raise ValueError(f"Invalid radius ({radius}).") + + geometry.radius = float(radius.value) + + # Unit vectors + geometry.basis_z = Vector3D( + sensor.x3_unit_vector.x, sensor.x3_unit_vector.y, sensor.x3_unit_vector.z + ) + geometry.basis_x = geometry.basis_z.orthogonal().normalise() + geometry.basis_y = geometry.basis_z.cross(geometry.basis_x).normalise() + + case GeometryType.OUTLINE: + # Outline coordinates for basis_x and basis_y + geometry.coords = np.vstack((sensor.outline.x2, sensor.outline.x1)) + + # Unit vectors + geometry.basis_x = Vector3D( + sensor.x2_unit_vector.x, sensor.x2_unit_vector.y, sensor.x2_unit_vector.z + ) + geometry.basis_y = Vector3D( + sensor.x1_unit_vector.x, sensor.x1_unit_vector.y, sensor.x1_unit_vector.z + ) + geometry.basis_z = Vector3D( + sensor.x3_unit_vector.x, sensor.x3_unit_vector.y, sensor.x3_unit_vector.z + ) + + # Surface area + surface = getattr(sensor, "surface", None) + geometry.surface = float(surface.value) if surface is not None else 0.0 + + return geometry diff --git a/src/cherab/imas/ids/bolometer/utility.py b/src/cherab/imas/ids/bolometer/utility.py new file mode 100644 index 0000000..3c4b915 --- /dev/null +++ b/src/cherab/imas/ids/bolometer/utility.py @@ -0,0 +1,101 @@ +# Copyright 2023 Euratom +# Copyright 2023 United Kingdom Atomic Energy Authority +# Copyright 2023 Centro de Investigaciones Energéticas, Medioambientales y Tecnológicas +# +# Licensed under the EUPL, Version 1.1 or – as soon they will be approved by the +# European Commission - subsequent versions of the EUPL (the "Licence"); +# You may not use this work except in compliance with the Licence. +# You may obtain a copy of the Licence at: +# +# https://joinup.ec.europa.eu/software/page/eupl5 +# +# Unless required by applicable law or agreed to in writing, software distributed +# under the Licence is distributed on an "AS IS" basis, WITHOUT WARRANTIES OR +# CONDITIONS OF ANY KIND, either express or implied. +# +# See the Licence for the specific language governing permissions and limitations +# under the Licence. +"""Module for bolometer utility functions.""" + +from enum import Enum +from typing import Self + +__all__ = ["CameraType", "GeometryType"] + + +class CameraType(Enum): + """Enum for camera type. + + The type of the bolometer camera. + + Attributes + ---------- + PINHOLE + The camera is a pinhole camera. + COLLIMATOR + The camera is a collimator camera. + OTHER + The camera is of other type. + """ + + PINHOLE = 1 + COLLIMATOR = 2 + OTHER = 0 + + @classmethod + def from_value(cls, value: int) -> Self: + """Get the camera type from a value. + + Parameters + ---------- + value + The integer value to convert to a camera type. + If the value is not a valid camera type, the default is `OTHER`. + + Returns + ------- + `.CameraType` + The corresponding `.CameraType` enum member. + """ + if value in cls._value2member_map_: + return cls(value) + return cls.OTHER # pyright: ignore[reportReturnType] + + +class GeometryType(Enum): + """Enum for geometry type. + + The geometry type of the bolometer foil or slit. + + Attributes + ---------- + OUTLINE + The geometry is defined by an outline. + CIRCULAR + The geometry is circular. + RECTANGLE + The geometry is rectangular. + """ + + OUTLINE = 1 + CIRCULAR = 2 + RECTANGLE = 3 + + @classmethod + def from_value(cls, value: int) -> Self: + """Get the geometry type from a value. + + Parameters + ---------- + value + The integer value to convert to a geometry type. + If the value is not a valid geometry type, the default is `RECTANGLE`. + + Returns + ------- + `.GeometryType` + The corresponding `.GeometryType` enum member. + """ + if value in cls._value2member_map_: + return cls(value) + return cls.RECTANGLE # pyright: ignore[reportReturnType] diff --git a/src/cherab/imas/ids/common/slice.py b/src/cherab/imas/ids/common/slice.py index 1711941..93ecd89 100644 --- a/src/cherab/imas/ids/common/slice.py +++ b/src/cherab/imas/ids/common/slice.py @@ -21,7 +21,7 @@ from numpy import inf -from imas import DBEntry +from imas.db_entry import DBEntry from imas.ids_defs import CLOSEST_INTERP from imas.ids_toplevel import IDSToplevel diff --git a/src/cherab/imas/ids/wall/load3d.py b/src/cherab/imas/ids/wall/load3d.py index a28ee2f..29131ee 100644 --- a/src/cherab/imas/ids/wall/load3d.py +++ b/src/cherab/imas/ids/wall/load3d.py @@ -46,8 +46,17 @@ def load_wall_3d( dict[str, dict[str, np.ndarray]] Dictionary of wall components defined by vertices and triangles. The dictionary keys for components are assigns as follows: - ``"{grid_name}.{subset_name}.{material_name}"`` - E.g.: ``"FullTokamak.full_main_chamber_wall.Be"``. + ``"{grid_name}.{subset_name}.{material_name}"``. + So, the returning dictionary looks like: + ```python + { + "FullTokamak.full_main_chamber_wall.Be": { + "vertices": np.ndarray of shape (N, 3), + "triangles": np.ndarray of shape (M, 3), + }, + ... + } + ``` Raises ------ diff --git a/src/cherab/imas/observer/__init__.py b/src/cherab/imas/observer/__init__.py new file mode 100644 index 0000000..a10d71c --- /dev/null +++ b/src/cherab/imas/observer/__init__.py @@ -0,0 +1,22 @@ +# Copyright 2023 Euratom +# Copyright 2023 United Kingdom Atomic Energy Authority +# Copyright 2023 Centro de Investigaciones Energéticas, Medioambientales y Tecnológicas +# +# Licensed under the EUPL, Version 1.1 or – as soon they will be approved by the +# European Commission - subsequent versions of the EUPL (the "Licence"); +# You may not use this work except in compliance with the Licence. +# You may obtain a copy of the Licence at: +# +# https://joinup.ec.europa.eu/software/page/eupl5 +# +# Unless required by applicable law or agreed to in writing, software distributed +# under the Licence is distributed on an "AS IS" basis, WITHOUT WARRANTIES OR +# CONDITIONS OF ANY KIND, either express or implied. +# +# See the Licence for the specific language governing permissions and limitations +# under the Licence. +"""Subpackage for creating observer objects from IMAS.""" + +from .bolometer import load_bolometers + +__all__ = ["load_bolometers"] diff --git a/src/cherab/imas/observer/bolometer.py b/src/cherab/imas/observer/bolometer.py new file mode 100644 index 0000000..4a22ab7 --- /dev/null +++ b/src/cherab/imas/observer/bolometer.py @@ -0,0 +1,410 @@ +# Copyright 2023 Euratom +# Copyright 2023 United Kingdom Atomic Energy Authority +# Copyright 2023 Centro de Investigaciones Energéticas, Medioambientales y Tecnológicas +# +# Licensed under the EUPL, Version 1.1 or – as soon they will be approved by the +# European Commission - subsequent versions of the EUPL (the "Licence"); +# You may not use this work except in compliance with the Licence. +# You may obtain a copy of the Licence at: +# +# https://joinup.ec.europa.eu/software/page/eupl5 +# +# Unless required by applicable law or agreed to in writing, software distributed +# under the Licence is distributed on an "AS IS" basis, WITHOUT WARRANTIES OR +# CONDITIONS OF ANY KIND, either express or implied. +# +# See the Licence for the specific language governing permissions and limitations +# under the Licence. +"""Module for loading bolometer cameras from IMAS bolometer IDS.""" + +from __future__ import annotations + +from raysect.core.constants import ORIGIN +from raysect.core.math import Point3D, Vector3D, rotate_basis, translate +from raysect.core.scenegraph._nodebase import _NodeBase +from raysect.optical.material import AbsorbingSurface +from raysect.primitive import Box, Subtract, Union +from raysect.primitive.csg import CSGPrimitive + +from cherab.tools.observers.bolometry import BolometerCamera, BolometerFoil, BolometerSlit +from imas.db_entry import DBEntry + +from ..ids.bolometer import load_cameras +from ..ids.bolometer._camera import BoloCamera, Geometry +from ..ids.bolometer.utility import CameraType, GeometryType +from ..ids.common import get_ids_time_slice + +__all__ = ["load_bolometers"] + +X_AXIS = Vector3D(1, 0, 0) +Y_AXIS = Vector3D(0, 1, 0) +Z_AXIS = Vector3D(0, 0, 1) +THICKNESS = 2.0e-3 # Thickness of camera box walls +THICKNESS_INNER_LAYER = 0.5e-3 # Thickness of inner aperture layers +EPS = 1e-4 # Small epsilon value to avoid numerical issues + + +def load_bolometers(*args, parent: _NodeBase | None = None, **kwargs) -> list[BolometerCamera]: + """Load bolometer cameras from IMAS bolometer IDS. + + .. note:: + This function requires the Data dictionary v4.1.0 or later. + + Parameters + ---------- + *args + Arguments passed to `~imas.db_entry.DBEntry`. + parent + The parent node of `~cherab.tools.observers.bolometry.BolometerCamera` in the Raysect + scene-graph, by default None. + **kwargs + Keyword arguments passed to `~imas.db_entry.DBEntry` constructor. + + Returns + ------- + `list[BolometerCamera]` + List of `~cherab.tools.observers.bolometry.BolometerCamera` objects. + + Examples + -------- + >>> from raysect.optical import World + >>> world = World() + + If you have a local IMAS database and store the "bolometer.h5" file there: + + >>> bolometers = load_bolometers("imas:hdf5?path=path/to/db/", "r", parent=world) + + If you want to load netCDF files directly: + + >>> bolometers = load_bolometers("path/to/bolometer_file.nc", "r", parent=world) + """ + # Load bolometer IDS + with DBEntry(*args, **kwargs) as entry: + # Get available time slices + ids = get_ids_time_slice(entry, "bolometer") + + # Extract bolometer data + bolo_data = load_cameras(ids) + + bolometers: list[BolometerCamera] = [] + + for data in bolo_data: + # Skip empty cameras + if len(data.channels) == 0: + continue + + # ------------------ + # === Camera Box === + # ------------------ + camera_box = _create_camera_box(data) + camera = BolometerCamera(camera_geometry=camera_box, name=data.name, parent=parent) + + match data.type: + case CameraType.PINHOLE: + # ---------------------- + # === Slit (Pinhole) === + # ---------------------- + # Pick up only first aperture and use it for all channels + slit_data = data.channels[0].slits[0] + slit = BolometerSlit( + f"slit-{data.name}", + slit_data.centre, + slit_data.basis_x, + slit_data.dx, + slit_data.basis_y, + slit_data.dy, + curvature_radius=(slit_data.radius or 0.0) + if slit_data.type == GeometryType.CIRCULAR + else 0.0, + parent=camera, + ) + case CameraType.COLLIMATOR: + slit = None # Defined per channel below + case _: + raise NotImplementedError(f"Camera type {data.type} not supported yet.") + + for i_channel, channel in enumerate(data.channels): + # ------------------------- + # === Slit (Collimator) === + # ------------------------- + # Concatenate top plate slits into one slit if multiple sub-collimators + # NOTE: Top plate slits are accumulated in the beginning of the slits list + # NOTE: All slit geometry types are assumed to be the RECTANGLE. + # NOTE: Sub-collimators are assumed to be aligned in a row along the basis_y direction. + if (num_subcol := channel.num_subcollimator) > 1: + _slit_data = channel.slits[:num_subcol] + centre = ( + ORIGIN + + sum([ORIGIN.vector_to(s.centre) for s in _slit_data], Vector3D(0, 0, 0)) + / num_subcol + ) + basis_x = _slit_data[0].basis_x + basis_y = _slit_data[0].basis_y + dx = _slit_data[0].dx + dy = ( + ORIGIN + + ORIGIN.vector_to(_slit_data[0].centre) + - _slit_data[0].dy * 0.5 * basis_y + ).distance_to( + ORIGIN + + ORIGIN.vector_to(_slit_data[-1].centre) + + _slit_data[-1].dy * 0.5 * basis_y + ) + + slit_data = Geometry( + type=GeometryType.RECTANGLE, + centre=centre, + basis_x=basis_x, + basis_y=basis_y, + dx=dx, + dy=dy, + ) + else: + slit_data = channel.slits[0] + + if data.type == CameraType.COLLIMATOR: + # Create slit object + match slit_data.type: + case GeometryType.CIRCULAR | GeometryType.RECTANGLE: + slit = BolometerSlit( + f"slit-{data.name}-ch{i_channel}", + slit_data.centre, + slit_data.basis_x, + slit_data.dx, + slit_data.basis_y, + slit_data.dy, + curvature_radius=(slit_data.radius or 0.0) + if slit_data.type == GeometryType.CIRCULAR + else 0.0, + parent=camera, + ) + case _: + raise NotImplementedError("Outline geometry not supported yet.") + + # ------------ + # === Foil === + # ------------ + match channel.foil.type: + case GeometryType.CIRCULAR | GeometryType.RECTANGLE: + foil = BolometerFoil( + f"foil-{data.name}-ch{i_channel}", + channel.foil.centre, + channel.foil.basis_x, + channel.foil.dx, + channel.foil.basis_y, + channel.foil.dy, + slit, + curvature_radius=(channel.foil.radius or 0.0) + if channel.foil.type == GeometryType.CIRCULAR + else 0.0, + parent=camera, + ) + case _: + raise NotImplementedError("Outline geometry not supported yet.") + + # Add foil to the camera + camera.add_foil_detector(foil) + + bolometers.append(camera) + + return bolometers + + +def _create_camera_box(bolo_data: BoloCamera) -> CSGPrimitive: + """Create a camera housing box geometry. + + This box is represented as a Subtract primitive, where a smaller inner box is subtracted + from a larger outer box to create a hollow structure. + Additionally, top plate apertures and other apertures (inner layers) are included as Unions. + + Returns + ------- + `CSGPrimitive` + The camera housing box geometry. + """ + # Extract slits data + slits_top: list[Geometry] + slits_inner: list[list[Geometry]] = [] + if bolo_data.type == CameraType.PINHOLE: + # Pick up only first aperture and use it for all channels + slits_top = [bolo_data[0].slits[0]] + else: + # Consider sub-collimators for collimator cameras + num_subcol = bolo_data[0].num_subcollimator + num_channel = len(bolo_data) + _s: list[Geometry] = [] + slits_top = sum( + [bolo_data[i_ch].slits[:num_subcol] for i_ch in range(num_channel)], + _s, + ) + slits_inner = [ + sum( + [ + bolo_data[i_ch].slits[i_slit : i_slit + num_subcol] + for i_ch in range(num_channel) + ], + _s, + ) + for i_slit in range(num_subcol, len(bolo_data[0].slits), num_subcol) + ] + + # ------------------------------- + # === Local coordinate system === + # ------------------------------- + # The local origin is placed at the average position of top plate slits of all foil detectors, + # and the basis vectors are defined based on the average normal and basis_x vectors of the slits. + origin = ORIGIN + sum( + [ORIGIN.vector_to(slit.centre) for slit in slits_top], + Vector3D(0, 0, 0), + ) / len(slits_top) + + basis_z: Vector3D = sum( + [slit.basis_z for slit in slits_top], + Vector3D(0, 0, 0), + ).normalise() + + basis_x: Vector3D = sum( + [slit.basis_x for slit in slits_top], + Vector3D(0, 0, 0), + ).normalise() + + basis_y = basis_z.cross(basis_x).normalise() + + # Transformation matrix from local to global coordinate system + to_root = translate(*origin) * rotate_basis(basis_z, basis_y) + + # ------------------------------ + # === Determine box geometry === + # ------------------------------ + EPS_WIDTH = 1e-3 + EPS_HEIGHT = 1e-3 + EPS_DEPTH = 1e-3 + + # Collect all corner points of foils and slits + pts: list[Point3D] = [] + for geometry in [bolo_data[i].foil for i in range(len(bolo_data))] + slits_top: + pts += _get_verts(geometry) + + # Camera box dimensions (in local coordinate system) + box_width_u = max([origin.vector_to(p).dot(basis_x) for p in pts]) + EPS_WIDTH + box_width_l = min([origin.vector_to(p).dot(basis_x) for p in pts]) - EPS_WIDTH + box_height_u = max([origin.vector_to(p).dot(basis_y) for p in pts]) + EPS_HEIGHT + box_height_l = min([origin.vector_to(p).dot(basis_y) for p in pts]) - EPS_HEIGHT + box_depth_u = max([origin.vector_to(p).dot(basis_z) for p in pts]) # slit top plate + box_depth_l = min([origin.vector_to(p).dot(basis_z) for p in pts]) - EPS_DEPTH + + # ------------------------- + # === Create Hollow Box === + # ------------------------- + camera_box = Box( + lower=ORIGIN + box_width_l * X_AXIS + box_height_l * Y_AXIS + box_depth_l * Z_AXIS, + upper=ORIGIN + box_width_u * X_AXIS + box_height_u * Y_AXIS + box_depth_u * Z_AXIS, + name=f"inner-box-{bolo_data.name}", + ) + outside_box = Box( + lower=camera_box.lower - Vector3D(THICKNESS, THICKNESS, THICKNESS), + upper=camera_box.upper + Vector3D(THICKNESS, THICKNESS, THICKNESS), + name=f"outside-box-{bolo_data.name}", + ) + camera_box = Subtract(outside_box, camera_box, name=f"hollow-box-{bolo_data.name}") + + # ------------------------------------------ + # === Clip out Top Plate Apertures Holes === + # ------------------------------------------ + for slit in slits_top: + slit_verts = _get_verts(slit) + aperture_box = Box( + lower=ORIGIN + + min([origin.vector_to(p).dot(basis_x) for p in slit_verts]) * X_AXIS + + min([origin.vector_to(p).dot(basis_y) for p in slit_verts]) * Y_AXIS + + (box_depth_u - EPS) * Z_AXIS, + upper=ORIGIN + + max([origin.vector_to(p).dot(basis_x) for p in slit_verts]) * X_AXIS + + max([origin.vector_to(p).dot(basis_y) for p in slit_verts]) * Y_AXIS + + (box_depth_u + THICKNESS + EPS) * Z_AXIS, + name=f"aperture-box-top-{bolo_data.name}", + ) + camera_box = Subtract(camera_box, aperture_box, name=f"aperture-top-{bolo_data.name}") + + # ---------------------------------- + # === Add Inner Apertures Plates === + # ---------------------------------- + # NOTO: Assume all inner apertures use the same local coordinate system as the top plate. + for slits in slits_inner: + # Create Inner Aperture Plate layer + layer_depth_z = basis_z.dot(origin.vector_to(slits[0].centre)) + layer = Box( + lower=Point3D( + outside_box.lower.x, + outside_box.lower.y, + layer_depth_z - THICKNESS_INNER_LAYER * 0.5, + ), + upper=Point3D( + outside_box.upper.x, + outside_box.upper.y, + layer_depth_z + THICKNESS_INNER_LAYER * 0.5, + ), + name=f"inner-plate-{bolo_data.name}", + ) + + for slit in slits: + slit_verts = _get_verts(slit) + aperture_box = Box( + lower=ORIGIN + + min([origin.vector_to(p).dot(basis_x) for p in slit_verts]) * X_AXIS + + min([origin.vector_to(p).dot(basis_y) for p in slit_verts]) * Y_AXIS + + (layer_depth_z - THICKNESS_INNER_LAYER * 0.5 - EPS) * Z_AXIS, + upper=ORIGIN + + max([origin.vector_to(p).dot(basis_x) for p in slit_verts]) * X_AXIS + + max([origin.vector_to(p).dot(basis_y) for p in slit_verts]) * Y_AXIS + + (layer_depth_z + THICKNESS_INNER_LAYER * 0.5 + EPS) * Z_AXIS, + name=f"inner-aperture-box-{bolo_data.name}", + ) + layer = Subtract(layer, aperture_box, name=f"inner-aperture-plate-{bolo_data.name}") + + camera_box = Union(camera_box, layer) + + # Transform to global coordinate system + camera_box.transform = to_root + + # Apply absorbing material + camera_box.material = AbsorbingSurface() + + # Name the camera box + camera_box.name = f"camera-box-{bolo_data.name}" + + return camera_box + + +def _get_verts(geometry: Geometry) -> list[Point3D]: + """Get the geometry vertices. + + For example, if the geometry is rectangular, return 4 corner vertices. + + Parameters + ---------- + geometry + Geometry structure object. + + Returns + ------- + `list[Point3D]` + List of geometry vertices. + """ + dx = geometry.dx + dy = geometry.dy + basis_x = geometry.basis_x + basis_y = geometry.basis_y + center = geometry.centre + match geometry.type: + case GeometryType.RECTANGLE: + verts = [ + center + 0.5 * dx * basis_x + 0.5 * dy * basis_y, + center - 0.5 * dx * basis_x + 0.5 * dy * basis_y, + center - 0.5 * dx * basis_x - 0.5 * dy * basis_y, + center + 0.5 * dx * basis_x - 0.5 * dy * basis_y, + ] + case _: + raise NotImplementedError(f"Geometry type {geometry.type} not implemented yet.") + + return verts