From 40863c7a330c91e61f0a8cf7c4816598c5ea4b98 Mon Sep 17 00:00:00 2001 From: Jack Lovell Date: Mon, 4 Mar 2024 10:58:54 +0000 Subject: [PATCH 1/5] Improve regularisation documentation and add ADMT demo * Add a Regularisation section to the tomography documentation, describing the functionality inside the admt_utils module. * Fixup the docstrings in admt_utils. * Add calculation of "skewed" second derivative operators, which operate along the diagonals of the grid. * Add a bolometer diagonstic system to Generomak. * Add a new example using the Generomak bolometers and the regularisation operators to perform isotropic and ADMT inversions. --- cherab/generomak/diagnostics/__init__.py | 1 + cherab/generomak/diagnostics/bolometers.py | 167 ++++++++++ cherab/tools/inversions/admt_utils.py | 98 ++++-- .../bolometry/admt_tomographic_inversion.py | 285 ++++++++++++++++++ docs/source/tools/tomography.rst | 33 ++ 5 files changed, 556 insertions(+), 28 deletions(-) create mode 100644 cherab/generomak/diagnostics/__init__.py create mode 100644 cherab/generomak/diagnostics/bolometers.py create mode 100644 demos/observers/bolometry/admt_tomographic_inversion.py diff --git a/cherab/generomak/diagnostics/__init__.py b/cherab/generomak/diagnostics/__init__.py new file mode 100644 index 00000000..ee4a0fc8 --- /dev/null +++ b/cherab/generomak/diagnostics/__init__.py @@ -0,0 +1 @@ +from .bolometers import load_bolometers diff --git a/cherab/generomak/diagnostics/bolometers.py b/cherab/generomak/diagnostics/bolometers.py new file mode 100644 index 00000000..e4693712 --- /dev/null +++ b/cherab/generomak/diagnostics/bolometers.py @@ -0,0 +1,167 @@ +""" +Some foil bolometers for measuring total radiated power. +""" +from raysect.core import (Node, Point3D, Vector3D, rotate_basis, + rotate_x, rotate_y, rotate_z, translate) +from raysect.core.math import rotate +from raysect.optical.material import AbsorbingSurface +from raysect.primitive import Box, Subtract + +from cherab.tools.observers import BolometerCamera, BolometerSlit, BolometerFoil + + +# Convenient constants +XAXIS = Vector3D(1, 0, 0) +YAXIS = Vector3D(0, 1, 0) +ZAXIS = Vector3D(0, 0, 1) +ORIGIN = Point3D(0, 0, 0) +# Bolometer geometry, independent of camera. +BOX_WIDTH = 0.05 +BOX_WIDTH = 0.2 +BOX_HEIGHT = 0.07 +BOX_DEPTH = 0.2 +SLIT_WIDTH = 0.004 +SLIT_HEIGHT = 0.005 +FOIL_WIDTH = 0.0013 +FOIL_HEIGHT = 0.0038 +FOIL_CORNER_CURVATURE = 0.0005 +FOIL_SEPARATION = 0.00508 # 0.2 inch between foils + + +def _make_bolometer_camera(slit_sensor_separation, sensor_angles): + """ + Build a single bolometer camera. + + The camera consists of a box with a rectangular slit and 4 sensors, + each of which has 4 foils. + In its local coordinate system, the camera's slit is located at the + origin and the sensors below the z=0 plane, looking up towards the slit. + """ + camera_box = Box(lower=Point3D(-BOX_WIDTH / 2, -BOX_HEIGHT / 2, -BOX_DEPTH), + upper=Point3D(BOX_WIDTH / 2, BOX_HEIGHT / 2, 0)) + # Hollow out the box + inside_box = Box(lower=camera_box.lower + Vector3D(1e-5, 1e-5, 1e-5), + upper=camera_box.upper - Vector3D(1e-5, 1e-5, 1e-5)) + camera_box = Subtract(camera_box, inside_box) + # The slit is a hole in the box + aperture = Box(lower=Point3D(-SLIT_WIDTH / 2, -SLIT_HEIGHT / 2, -1e-4), + upper=Point3D(SLIT_WIDTH / 2, SLIT_HEIGHT / 2, 1e-4)) + camera_box = Subtract(camera_box, aperture) + camera_box.material = AbsorbingSurface() + bolometer_camera = BolometerCamera(camera_geometry=camera_box) + # The bolometer slit in this instance just contains targeting information + # for the ray tracing, since we have already given our camera a geometry + # The slit is defined in the local coordinate system of the camera + slit = BolometerSlit(slit_id="Example slit", centre_point=ORIGIN, + basis_x=XAXIS, dx=SLIT_WIDTH, basis_y=YAXIS, dy=SLIT_HEIGHT, + parent=bolometer_camera) + for j, angle in enumerate(sensor_angles): + # 4 bolometer foils, spaced at equal intervals along the local X axis + sensor = Node(name="Bolometer sensor", parent=bolometer_camera, + transform=rotate_y(angle) * translate(0, 0, -slit_sensor_separation)) + for i, shift in enumerate([-1.5, -0.5, 0.5, 1.5]): + # Note that the foils will be parented to the camera rather than the + # sensor, so we need to define their transform relative to the camera. + foil_transform = sensor.transform * translate(shift * FOIL_SEPARATION, 0, 0) + foil = BolometerFoil(detector_id="Foil {} sensor {}".format(i + 1, j + 1), + centre_point=ORIGIN.transform(foil_transform), + basis_x=XAXIS.transform(foil_transform), dx=FOIL_WIDTH, + basis_y=YAXIS.transform(foil_transform), dy=FOIL_HEIGHT, + slit=slit, parent=bolometer_camera, units="Power", + accumulate=False, curvature_radius=FOIL_CORNER_CURVATURE) + bolometer_camera.add_foil_detector(foil) + return bolometer_camera + + +def load_bolometers(parent=None): + """ + Load the Generomak bolometers. + + The Generomak bolometer diagnostic consists of multiple 12-channel + cameras. Each camera has 3 4-channel sensors inside. + + * 2 cameras are located at the midplane with purely-poloidal, + horizontal views. + * 1 camera is located at the top of the machine with purely-poloidal, + vertical views. + * 2 cameras have purely tangential views at the midplane. + * 1 camera has combined poloidal+tangential views, which look like + curved lines of sight in the poloidal plane. It looks at the lower + divertor. + + :param parent: the scenegraph node the bolometers will belong to. + :return: a list of BolometerCamera instances, one for each of the + cameras described above. + """ + poloidal_camera_rotations = [ + 30, -30, # Horizontal poloidal, + -90, # Vertical poloidal, + 0, # Tangential, + 0, # Combined poloidal/tangential, + ] + toroidal_camera_rotations = [ + 0, 0, # Horizontal poloidal + 0, # Vertical poloidal + -40, # Tangential + 40, # Combined poloidal/tangential + ] + radial_camera_rotations = [ + 0, 0, # Horizontal poloidal + 0, # Vertical poloidal + 90, # Tangential + 90, # Combined poloidal/tangential + ] + camera_origins = [ + Point3D(2.5, 0.05, 0), Point3D(2.5, -0.05, 0), # Horizontal poloidal + Point3D(1.3, 0, 1.42), # Vertical poloidal + Point3D(2.5, 0, 0), # Midplane tangential horizontal + Point3D(2.5, 0, -0.2), # Combined poloidal/tangential + ] + slit_sensor_separations = [ + 0.08, 0.08, # Horizontal poloidal + 0.05, # Vertical poloidal + 0.1, # Tangential + 0.1, # Combined poloidal/tangential + ] + all_sensor_angles = [ + [-22.5, -7.5, 7.5, 22.5], [-22.5, -7.5, 7.5, 22.5], # Horizontal poloidal + [-36, -12, 12, 36], # Vertical poloidal + [-18, -6, 6, 18], # Tangential + [-18, -6, 6, 18], # Combined poloidal/tangential + ] + toroidal_angles = [ + 10, 10, # Horizontal poloidal, need to avoid LFS limiters. + 0, # Vertical poloidal, happy to hit LFS limiters. + -15, # Tangential, avoid LFS limiters. + 15, # Combined poloidal/tangential, avoid LFS limiters. + ] + names = [ + 'HozPol1', 'HozPol2', + 'VertPol', + 'TanMid1', + 'TanPol1', + ] + cameras = [] + # FIXME: this for loop definition is really ugly! + for ( + poloidal_rotation, toroidal_rotation, radial_rotation, camera_origin, + slit_sensor_separation, sensor_angles, toroidal_angle, name + ) in zip( + poloidal_camera_rotations, toroidal_camera_rotations, radial_camera_rotations, + camera_origins, + slit_sensor_separations, all_sensor_angles, toroidal_angles, names + ): + camera = _make_bolometer_camera(slit_sensor_separation, sensor_angles) + # FIXME: work out how to combine tangential and poloidal rotations. + camera.transform = ( + rotate_z(toroidal_angle) + * translate(camera_origin.x, camera_origin.y, camera_origin.z) + * rotate_z(toroidal_rotation) + * rotate_x(radial_rotation) + * rotate_y(poloidal_rotation + 90) + * rotate_basis(-ZAXIS, YAXIS) + ) + camera.parent = parent + camera.name = "{} at angle {}".format(name, poloidal_rotation) + cameras.append(camera) + return cameras diff --git a/cherab/tools/inversions/admt_utils.py b/cherab/tools/inversions/admt_utils.py index 549c4d42..74fac859 100644 --- a/cherab/tools/inversions/admt_utils.py +++ b/cherab/tools/inversions/admt_utils.py @@ -36,15 +36,15 @@ def generate_derivative_operators(voxel_vertices, grid_index_1d_to_2d_map, Generate the first and second derivative operators for a regular grid. :param ndarray voxel_vertices: an Nx4x2 array of coordinates of the - vertices of each voxel, (R, Z) + vertices of each voxel, (R, Z) :param dict grid_1d_to_2d_map: a mapping from the 1D array of - voxels in the grid to a 2D array of voxels if they were arranged - spatially. + voxels in the grid to a 2D array of voxels if they were arranged + spatially. :param dict grid_2d_to_1d_map: the inverse mapping from a 2D - spatially-arranged array of voxels to the 1D array. + spatially-arranged array of voxels to the 1D array. - :return dict operators: a dictionary containing the derivative - operators: Dij for i, y ∊ (x, y) and Di for i ∊ (x, y). + :return: a dictionary containing the derivative operators: Dij for + i, j ∊ (x, y) and Di for i ∊ (x, y), Dsp and Dsm. This function assumes that all voxels are rectilinear, with their axes aligned to the coordinate axes. Additionally, all voxels are @@ -62,13 +62,18 @@ def generate_derivative_operators(voxel_vertices, grid_index_1d_to_2d_map, D_{xx} \equiv \frac{\partial^2}{\partial x^2}\\ D_{xy} \equiv \frac{\partial^2}{\partial x \partial y} - etc. + etc. It also produces two additional operators, Dsp and Dsm, for + second derivatives on the dy/dx = 1 and dy/dx = -1 diagonals + respectively. Note that the standard 2D laplacian (for isotropic regularisation) - can be trivially calculated as L = Dxx * dx + Dyy * dy, where dx and - dy are the voxel width and height respectively. This expression does - not however produce the 2D laplacian derived from the N-dimensional - case. + can be trivially calculated as follows: + + .. math:: + L = (1 - \alpha) (D_{xx} + D_{yy}) + (\alpha / 2) (D_{sp} + D_{sm}) + + α = 2/3 produces the operator used in Carr et. al. RSI 89, 083506 (2018). + α = 1/3 produces the operator with optimal isotropy. """ # Input argument validation: assume rectilinear voxels voxel_vertices = np.asarray(voxel_vertices) @@ -87,6 +92,8 @@ def generate_derivative_operators(voxel_vertices, grid_index_1d_to_2d_map, Dxx = np.zeros((num_cells, num_cells)) Dxy = np.zeros((num_cells, num_cells)) Dyy = np.zeros((num_cells, num_cells)) + Dsp = np.zeros((num_cells, num_cells)) + Dsm = np.zeros((num_cells, num_cells)) # TODO: for now, we assume all voxels have rectangular cross sections # which are approximately identical. As per Ingesson's notation, we # assume voxels are ordered from top left to bottom right, in column-major @@ -99,6 +106,9 @@ def generate_derivative_operators(voxel_vertices, grid_index_1d_to_2d_map, dx = np.min(abs(dx[dx != 0])).item() dy = np.min(abs(dy[dy != 0])).item() + # Work out how the voxels are ordered: increasing/decreasing in x/y. + xinc, yinc = np.sign(cell_centres[-1] - cell_centres[0]) + # Note that iy increases as y decreases (cells go from top to bottom), # which is the same as Ingesson's notation in equations 37-41 # Use the second version of the second derivative boundary formulae, so @@ -110,8 +120,13 @@ def generate_derivative_operators(voxel_vertices, grid_index_1d_to_2d_map, # get the 2D mesh coordinates of this cell ix, iy = grid_index_1d_to_2d_map[ith_cell] + iright = ix + xinc + ileft = ix - xinc + iabove = iy + yinc + ibelow = iy - yinc + try: - n_left = grid_index_2d_to_1d_map[ix - 1, iy] # left of n0 + n_left = grid_index_2d_to_1d_map[ileft, iy] # left of n0 except KeyError: at_left = True else: @@ -119,7 +134,7 @@ def generate_derivative_operators(voxel_vertices, grid_index_1d_to_2d_map, Dxx[ith_cell, n_left] = 1 try: - n_below_left = grid_index_2d_to_1d_map[ix - 1, iy + 1] # below left of n0 + n_below_left = grid_index_2d_to_1d_map[ileft, ibelow] # below left of n0 except KeyError: # KeyError does not necessarily mean bottom AND left pass @@ -127,7 +142,7 @@ def generate_derivative_operators(voxel_vertices, grid_index_1d_to_2d_map, Dxy[ith_cell, n_below_left] = 1 / 4 try: - n_below = grid_index_2d_to_1d_map[ix, iy + 1] + n_below = grid_index_2d_to_1d_map[ix, ibelow] except KeyError: at_bottom = True else: @@ -135,14 +150,14 @@ def generate_derivative_operators(voxel_vertices, grid_index_1d_to_2d_map, Dyy[ith_cell, n_below] = 1 try: - n_below_right = grid_index_2d_to_1d_map[ix + 1, iy + 1] + n_below_right = grid_index_2d_to_1d_map[iright, ibelow] except KeyError: pass else: Dxy[ith_cell, n_below_right] = -1 / 4 try: - n_right = grid_index_2d_to_1d_map[ix + 1, iy] + n_right = grid_index_2d_to_1d_map[iright, iy] except KeyError: at_right = True else: @@ -150,14 +165,14 @@ def generate_derivative_operators(voxel_vertices, grid_index_1d_to_2d_map, Dxx[ith_cell, n_right] = 1 try: - n_above_right = grid_index_2d_to_1d_map[ix + 1, iy - 1] + n_above_right = grid_index_2d_to_1d_map[iright, iabove] except KeyError: pass else: Dxy[ith_cell, n_above_right] = 1 / 4 try: - n_above = grid_index_2d_to_1d_map[ix, iy - 1] + n_above = grid_index_2d_to_1d_map[ix, iabove] except KeyError: at_top = True else: @@ -165,7 +180,7 @@ def generate_derivative_operators(voxel_vertices, grid_index_1d_to_2d_map, Dyy[ith_cell, n_above] = 1 try: - n_above_left = grid_index_2d_to_1d_map[ix - 1, iy - 1] + n_above_left = grid_index_2d_to_1d_map[ileft, iabove] except KeyError: pass else: @@ -247,14 +262,42 @@ def generate_derivative_operators(voxel_vertices, grid_index_1d_to_2d_map, Dxy[ith_cell, ith_cell] = -1 Dxy[ith_cell, n_above_left] = -1 + if np.isnan(n_above_left) and not np.isnan(n_below_right): + Dsm[ith_cell, ith_cell] = -1 + Dsm[ith_cell, n_below_right] = 1 + elif np.isnan(n_below_right) and not np.isnan(n_above_left): + Dsm[ith_cell, ith_cell] = -1 + Dsm[ith_cell, n_above_left] = 1 + elif np.isnan(n_above_left) and np.isnan(n_below_right): + Dsm[ith_cell, ith_cell] = 0 + else: + Dsm[ith_cell, ith_cell] = -2 + Dsm[ith_cell, n_above_left] = 1 + Dsm[ith_cell, n_below_right] = 1 + + if np.isnan(n_above_right) and not np.isnan(n_below_left): + Dsp[ith_cell, ith_cell] = -1 + Dsp[ith_cell, n_below_left] = 1 + elif np.isnan(n_below_left) and not np.isnan(n_above_right): + Dsp[ith_cell, ith_cell] = -1 + Dsp[ith_cell, n_above_right] = 1 + elif np.isnan(n_below_left) and np.isnan(n_above_right): + Dsp[ith_cell, ith_cell] = 0 + else: + Dsp[ith_cell, ith_cell] = -2 + Dsp[ith_cell, n_above_right] = 1 + Dsp[ith_cell, n_below_left] = 1 + Dx = Dx / dx Dy = Dy / dy Dxx = Dxx / dx**2 Dyy = Dyy / dy**2 Dxy = Dxy / (dx * dy) + Dsp = Dsp / (dx**2 + dy**2) + Dsm = Dsm / (dx**2 + dy**2) # Package all operators up into a dictionary - operators = dict(Dx=Dx, Dy=Dy, Dxx=Dxx, Dyy=Dyy, Dxy=Dxy) + operators = dict(Dx=Dx, Dy=Dy, Dxx=Dxx, Dyy=Dyy, Dxy=Dxy, Dsp=Dsp, Dsm=Dsm) return operators @@ -263,22 +306,21 @@ def calculate_admt(voxel_radii, derivative_operators, psi_at_voxels, dx, dy, ani Calculate the ADMT regularisation operator. :param ndarray voxel_radii: a 1D array of the radius at the centre - of each voxel in the grid + of each voxel in the grid :param tuple derivative_operators: a named tuple with the derivative - operators for the grid, as returned by :func:generate_derivative_operators + operators for the grid, as returned by :func:generate_derivative_operators :param ndarray psi_at_voxels: the magnetic flux at the centre of - each voxel in the grid + each voxel in the grid :param float dx: the width of each voxel. :param float dy: the height of each voxel :param float anisotropy: the ratio of the smoothing in the parallel - and perpendicular directions. - - :return ndarray admt: the ADMT regularisation operator. + and perpendicular directions. + :return: the ADMT regularisation operator. The degree of anisotropy dictates the relative suppression of gradients in the directions parallel and perpendicular to the - magnetic field. For example, `anisotropy=10` implies parallel - gradients in solution are 10 times smaller than perpendicular + magnetic field. For example, ``anisotropy=10`` implies parallel + gradients in the solution are 10 times smaller than perpendicular gradients. This function assumes that all voxels are rectilinear, with their diff --git a/demos/observers/bolometry/admt_tomographic_inversion.py b/demos/observers/bolometry/admt_tomographic_inversion.py new file mode 100644 index 00000000..e27fc518 --- /dev/null +++ b/demos/observers/bolometry/admt_tomographic_inversion.py @@ -0,0 +1,285 @@ +""" +This example demonstrates performing a tomographic reconstruction of a +radiation profile using Cherab's anisotropic diffusion (ADMT) regularisation +utilities. We use the machine geometry, sample bolometers and equilibrium +from Generomak. +""" +import matplotlib.pyplot as plt +import numpy as np + +from raysect.core.math.function.float import Exp2D, Arg2D, Atan4Q2D +from raysect.core.math import translate +from raysect.optical import World +from raysect.optical.material import AbsorbingSurface, VolumeTransform +from raysect.primitive import Cylinder, Subtract + +from cherab.generomak.machine import load_first_wall +from cherab.generomak.equilibrium import load_equilibrium +from cherab.generomak.diagnostics import load_bolometers +from cherab.core.math import sample2d, sample2d_grid, sample2d_points, AxisymmetricMapper +from cherab.tools.emitters import RadiationFunction +from cherab.tools.raytransfer import RayTransferCylinder, RayTransferPipeline0D +from cherab.tools.inversions import admt_utils as admt +from cherab.tools.inversions import invert_regularised_nnls + + +plt.ion() + +################################################################################ +# Define the emissivity profile. +################################################################################ +# The emissivity profile consists of a blob, a ring and part of a ring on the LFS. +# The blob and the ring are Gaussian flux functions. +# The ring is Gaussian in flux and poloidal angle. +# All have equal maximum emissivities, but not necessarily equal total power. +# We use Raysect's function framework to specify an analytic form for the +# emissivity profile, as this is very quick to sample and ray trace. +eq = load_equilibrium() +psin = eq.psi_normalised +axis = eq.magnetic_axis +blob_centre_psin = 0 +blob_width_psin = 0.1 +blob = Exp2D(-0.5 * (psin - blob_centre_psin)**2 / (blob_width_psin**2)) +ring_centre_psin = 0.5 +ring_width_psin = 0.05 +ring = Exp2D(-0.5 * (psin - ring_centre_psin)**2 / (ring_width_psin**2)) +theta = Atan4Q2D(Arg2D('y') - axis.y, Arg2D('x') - axis.x) +lfs_centre_psin = 0.85 +lfs_width_psin = 0.1 +lfs_centre_theta = 0 +lfs_width_theta = 0.5 +lfs = Exp2D(-0.5 * (((psin - lfs_centre_psin) / lfs_width_psin)**2 + + ((theta - lfs_centre_theta) / lfs_width_theta)**2)) +emissivity = blob + ring + lfs +# Assume no emission from these contributors outside the separatrix. +emissivity = emissivity * eq.inside_lcfs + +# Visualise the emissivity profile with the equilibrium overlayed. +plt.figure() +rsamp, zsamp, psisamp = sample2d(psin, (*eq.r_range, 500), (*eq.z_range, 1000)) +plt.contour(rsamp, zsamp, psisamp.T, linewidths=0.5, alpha=0.3, + levels=np.linspace(0, 1, 10), colors=['k']*9 + ['red']) +rsamp, zsamp, emsamp = sample2d(emissivity, (*eq.r_range, 500), (*eq.z_range, 1000)) +im = plt.imshow(emsamp.T, extent=(rsamp[0], rsamp[-1], zsamp[0], zsamp[-1]), cmap='Purples') +plt.xlabel("R[m]") +plt.ylabel("Z[m]") +plt.colorbar(im, label="Model emissivity [W/m3]") +plt.xlim([rsamp[0], rsamp[-1]]) +plt.ylim([zsamp[0], zsamp[-1]]) +plt.gca().set_aspect('equal') +plt.pause(0.5) + + +################################################################################ +# Load the machine wall and diagnostic. +################################################################################ +print("Loading the geometry...") +world = World() +load_first_wall(world, material=AbsorbingSurface()) +bolos = load_bolometers(world) +# Only consider the purely-poloidal cameras for now... +bolos = bolos[:3] + +######################################################################## +# Produce a voxel grid +######################################################################## +print("Producing the voxel grid...") +# Define the centres of each voxel, as an (nx, ny, 2) array. +nx = 40 +ny = 85 +cell_r, cell_dx = np.linspace(0.7, 2.5, nx, retstep=True) +cell_z, cell_dz = np.linspace(-1.8, 1.6, ny, retstep=True) +cell_r_grid, cell_z_grid = np.broadcast_arrays(cell_r[:, None], cell_z[None, :]) +cell_centres = np.stack((cell_r_grid, cell_z_grid), axis=-1) # (nx, ny, 2) array + +# Define the positions of the vertices of the voxels. +cell_vertices_r = np.linspace(cell_r[0] - 0.5 * cell_dx, cell_r[-1] + 0.5 * cell_dx, nx + 1) +cell_vertices_z = np.linspace(cell_z[0] - 0.5 * cell_dz, cell_z[-1] + 0.5 * cell_dz, ny + 1) + +# Build a mask, only including cells within the wall. +mask_2d = sample2d_grid(eq.inside_limiter, cell_r, cell_z) +mask_3d = mask_2d[:, np.newaxis, :] +ncells = mask_3d.sum() + +# We'll use the Ray Transfer frameworks as these voxels are rectangular +# and it's much faster than the Voxel framework for simple cases like this. +ray_transfer_grid = RayTransferCylinder( + radius_outer=cell_vertices_r[-1], + radius_inner=cell_vertices_r[0], + height=cell_vertices_z[-1] - cell_vertices_z[0], + n_radius=nx, n_height=ny, mask=mask_3d, n_polar=1, + transform=translate(0, 0, cell_vertices_z[0]), +) + +######################################################################## +# Calculate the geometry matrix for the grid +######################################################################## +print("Calculating the geometry matrix...") +# The ray transfer object must be in the same world as the bolometers +ray_transfer_grid.parent = world + +sensitivity_matrix = [] +for camera in bolos: + for foil in camera: + # Temporarily override foil pipelines for the sensitivity calculation. + orig_pipelines = foil.pipelines + foil.pipelines = [RayTransferPipeline0D(kind=foil.units)] + # All objects in world have wavelength-independent material properties, + # so it doesn't matter which wavelength range we use (as long as + # max_wavelength - min_wavelength = 1) + foil.min_wavelength = 1 + foil.max_wavelength = 2 + foil.spectral_bins = ray_transfer_grid.bins + foil.observe() + sensitivity_matrix.append(foil.pipelines[0].matrix) + # Restore original pipelines for subsequent observe calls. + foil.pipelines = orig_pipelines +sensitivity_matrix = np.asarray(sensitivity_matrix) + + +################################################################################ +# Generate the regularisation operators. +################################################################################ +print("Generating regularisation operators...") +# generate_derivative_operators requires two mappings, one from a flat list of +# voxels to the original 2D grid, and one for the 2D grid coordinates to the +# flat list of voxels. We could build these by hand, but the RayTransferCylinder +# object helpfully provides the data already so we just need to convert from +# arrays to dictionaries. +grid_index_1d_to_2d_map = {} +for k, idx2d in enumerate(ray_transfer_grid.invert_voxel_map()): + # We want the x and z elements, as the Ray Transfer grid is 3D and this + # inversion is going to be in 2D. + grid_index_1d_to_2d_map[k] = (idx2d[0].item(), idx2d[2].item()) +grid_index_2d_to_1d_map = {} +nx, _, ny = ray_transfer_grid.voxel_map.shape +for i in range(nx): + for j in range(ny): + voxel_index = ray_transfer_grid.voxel_map[i, 0, j] + if voxel_index != -1: + grid_index_2d_to_1d_map[(i, j)] = voxel_index +# We now need an (Nx4x2) array of voxel vertices, which can be easily calculated. +voxel_centres = np.array([cell_centres[grid_index_1d_to_2d_map[i]] + for i in range(ray_transfer_grid.bins)]) +vertex_displacements = np.array([[-cell_dx/2, -cell_dz/2], + [-cell_dx/2, cell_dz/2], + [cell_dx/2, cell_dz/2], + [cell_dx/2, -cell_dz/2]]) +# Combine the (N,2) and (4,2) arrays to get an (N,4,2) array. +voxel_vertices = voxel_centres[:, None, :] + vertex_displacements[None, :, :] +derivative_operators = admt.generate_derivative_operators( + voxel_vertices, grid_index_1d_to_2d_map, grid_index_2d_to_1d_map +) + +# As described in the docstring for generate_derivative_operators, we can +# calculate a 2D laplacian operator for "isotropic" smoothing easily: +alpha = 1/3 # Optimal isotropy +aligned = derivative_operators['Dxx'] * cell_dx**2 + derivative_operators['Dyy'] * cell_dz**2 +skewed = (derivative_operators['Dsp'] + derivative_operators['Dsm']) * (cell_dx**2 + cell_dz**2) +laplacian = (1 - alpha) * aligned + (alpha / 2) * skewed +# We could also use alpha = 2/3, which would produce an operator akin to the one +# used in Carr et. al. RSI 89, 083506 (2018). + +# We can also derive an anistoropic regularisation operator, which calculates the +# amount of un-smoothness parallel and perpendicular to the magnetic field lines. +# For this we need the radii of the voxels and the magnetic flux at each voxel, +# along with a few other inputs. +voxel_radii = voxel_centres[:, 0] +psi_at_voxels = sample2d_points(eq.psi_normalised, voxel_centres) +# We also need to decide on the degree of anisotropy we expect, i.e. how much more +# smooth the radiation is along the field lines vs perpendicular to them. +# The optimal value will depend on the problem at hand. +anisotropy = 50 +admt_operator = admt.calculate_admt( + voxel_radii, derivative_operators, psi_at_voxels, cell_dx, cell_dz, anisotropy +) + +################################################################################ +# Forward model the measurements. +################################################################################ +print("Modelling the measurement values...") +# Create an emitting object whose emission is defined by the analytic form we +# produced earlier. As the emission depends on the equilibrium, this object +# should have an extent no larger than the equilibrium reconstruction extent. +# We actually make the emitter slightly smaller than the equilibrium region to +# avoid numerical precision issues creating attempts to calculate the emissivity +# outside of the equlibrium domain. +CYLINDER_RADIUS = eq.r_range[-1] - 1e-6 +CYLINDER_HEIGHT = eq.z_range[-1] - eq.z_range[0] - 2e-6 +CYLINDER_SHIFT = eq.z_range[0] + 1e-6 +emitter = Cylinder(radius=CYLINDER_RADIUS, height=CYLINDER_HEIGHT, + transform=translate(0, 0, CYLINDER_SHIFT)) +# Cut out middle of cylinder as well: equilibrium not defined here. +emitter = Subtract(emitter, Cylinder(radius=eq.r_range[0] + 1e-6, height=10, + transform=translate(0, 0, -5))) +emission_function_3d = AxisymmetricMapper(emissivity) +emitting_material = VolumeTransform(RadiationFunction(emission_function_3d), + transform=emitter.transform.inverse()) +emitter.material = emitting_material +emitter.parent = world + +# Calculate the line-integral bolometer measurements by observing the emitter +# with all bolometers. The measurements should have the same channel order as +# the sensitivity matrix. +all_measurements = [] +for camera in bolos: + all_measurements.extend(camera.observe()) + + +################################################################################ +# Perform the inversions. +################################################################################ +print("Performing inversions...") +# We'll use NNLS with regularisation. The hyperparameters have been chosen by +# hand but techniques such as the discrepancy principle or L curve optimisation +# could also be used to determine them. That is out of the scope of this demo. +isotropic_alpha = 1e-10 +isotropic_inversion, _ = invert_regularised_nnls( + sensitivity_matrix, all_measurements, alpha=isotropic_alpha, + tikhonov_matrix=laplacian, +) + +admt_alpha = 1e-10 +admt_inversion, _ = invert_regularised_nnls( + sensitivity_matrix, all_measurements, alpha=admt_alpha, + tikhonov_matrix=admt_operator, +) + + +################################################################################ +# Plot the inversion results. +################################################################################ +emiss2d = np.zeros((nx, ny)) + +# Isotropic +for index1d, indices2d in grid_index_1d_to_2d_map.items(): + emiss2d[indices2d] = isotropic_inversion[index1d] +plt.figure() +im = plt.imshow(emiss2d.T, extent=(cell_r[0], cell_r[-1], cell_z[0], cell_z[-1]), cmap='Purples') +plt.contour(rsamp, zsamp, psisamp.T, linewidths=0.5, alpha=0.3, + levels=np.linspace(0, 1, 10), colors=['k']*9 + ['red']) +plt.xlabel("R[m]") +plt.ylabel("Z[m]") +plt.colorbar(im, label="Inverted\nEmissivity [W/m3]") +plt.xlim([rsamp[0], rsamp[-1]]) +plt.ylim([zsamp[0], zsamp[-1]]) +plt.gca().set_aspect('equal') +plt.title("Isotropic regularisation") + +# Anisotropic. +for index1d, indices2d in grid_index_1d_to_2d_map.items(): + emiss2d[indices2d] = admt_inversion[index1d] +plt.figure() +im = plt.imshow(emiss2d.T, extent=(cell_r[0], cell_r[-1], cell_z[0], cell_z[-1]), cmap='Purples') +plt.contour(rsamp, zsamp, psisamp.T, linewidths=0.5, alpha=0.3, + levels=np.linspace(0, 1, 10), colors=['k']*9 + ['red']) +plt.xlabel("R[m]") +plt.ylabel("Z[m]") +plt.colorbar(im, label="Inverted\nEmissivity [W/m3]") +plt.xlim([rsamp[0], rsamp[-1]]) +plt.ylim([zsamp[0], zsamp[-1]]) +plt.gca().set_aspect('equal') +plt.title("Anisotropic regularisation") + +plt.ioff() +plt.show() diff --git a/docs/source/tools/tomography.rst b/docs/source/tools/tomography.rst index 3f3984f6..523f41a1 100644 --- a/docs/source/tools/tomography.rst +++ b/docs/source/tools/tomography.rst @@ -119,3 +119,36 @@ Use spectral pipelines from Raysect if you need these features. .. autoclass:: cherab.tools.raytransfer.pipelines.RayTransferPipeline1D .. autoclass:: cherab.tools.raytransfer.pipelines.RayTransferPipeline2D + + +Regularisation +-------------- + +Some of the inversion methods take a regularisation operator, which provides +additional constraints to help achieved unique solutions to ill-posed +tomography problems. Many regularisation schemes impose constraints on the smoothness +of the resulting solution, with this smoothness quantified by the second derivative +of the solution. Two such regularisation schemes are common in fusion applications: + +#. Isotropic smoothing, where the solution has the same smoothness in all directions. +#. Anisotropic smoothing, so-called "anisotropic diffusion model tomography" (ADMT), + where the solution is smoother parallel to the magnetic field and less smooth + perpendicular to the magnetic field. + +Cherab provides some utility functions to assist in calculating appropriate +operators using these (and other) derivative-based regularisation schemes. These can be used +on inversion grids defined using both the Voxel and Ray Transfer frameworks, and passed +directly to the inversion methods in Cherab which take regularisation operators, such as +cherab.tools.inversions.invert_constrained_sart and cherab.tools.inversions.invert_regularised_nnls. + +The routines to calculate derivative operators for inversion grids, and further to calculate +the ADMT operator for a given set of derivative operators and magnetic field, are taken from +work published by L. C. Ingesson in `JET-R(99)08`_. + + +.. autofunction:: cherab.tools.inversions.admt_utils.generate_derivative_operators + +.. autofunction:: cherab.tools.inversions.admt_utils.calculate_admt + + +.. _JET-R(99)08: http://www.euro-fusionscipub.org/wp-content/uploads/2014/11/JETR99008.pdf From e7cfa3dad7176222495fab3de5d0a55b3dde4da6 Mon Sep 17 00:00:00 2001 From: "Lovell, Jack J" Date: Tue, 19 May 2026 15:56:35 +0200 Subject: [PATCH 2/5] Tidy up the for loop creating generomak bolometer sensors --- cherab/generomak/diagnostics/bolometers.py | 135 ++++++++++----------- 1 file changed, 65 insertions(+), 70 deletions(-) diff --git a/cherab/generomak/diagnostics/bolometers.py b/cherab/generomak/diagnostics/bolometers.py index e4693712..e64d9dec 100644 --- a/cherab/generomak/diagnostics/bolometers.py +++ b/cherab/generomak/diagnostics/bolometers.py @@ -20,6 +20,7 @@ BOX_WIDTH = 0.2 BOX_HEIGHT = 0.07 BOX_DEPTH = 0.2 +THICKNESS = 1e-3 SLIT_WIDTH = 0.004 SLIT_HEIGHT = 0.005 FOIL_WIDTH = 0.0013 @@ -39,13 +40,13 @@ def _make_bolometer_camera(slit_sensor_separation, sensor_angles): """ camera_box = Box(lower=Point3D(-BOX_WIDTH / 2, -BOX_HEIGHT / 2, -BOX_DEPTH), upper=Point3D(BOX_WIDTH / 2, BOX_HEIGHT / 2, 0)) - # Hollow out the box - inside_box = Box(lower=camera_box.lower + Vector3D(1e-5, 1e-5, 1e-5), - upper=camera_box.upper - Vector3D(1e-5, 1e-5, 1e-5)) + # Hollow out the box: it has 1 mm thick walls. + inside_box = Box(lower=camera_box.lower + Vector3D(THICKNESS, THICKNESS, THICKNESS), + upper=camera_box.upper - Vector3D(THICKNESS, THICKNESS, THICKNESS)) camera_box = Subtract(camera_box, inside_box) - # The slit is a hole in the box - aperture = Box(lower=Point3D(-SLIT_WIDTH / 2, -SLIT_HEIGHT / 2, -1e-4), - upper=Point3D(SLIT_WIDTH / 2, SLIT_HEIGHT / 2, 1e-4)) + # The slit is a hole in the box. Make it thicker than the wall. + aperture = Box(lower=Point3D(-SLIT_WIDTH / 2, -SLIT_HEIGHT / 2, -1.1 * THICKNESS), + upper=Point3D(SLIT_WIDTH / 2, SLIT_HEIGHT / 2, 0.1 * THICKNESS)) camera_box = Subtract(camera_box, aperture) camera_box.material = AbsorbingSurface() bolometer_camera = BolometerCamera(camera_geometry=camera_box) @@ -93,75 +94,69 @@ def load_bolometers(parent=None): :return: a list of BolometerCamera instances, one for each of the cameras described above. """ - poloidal_camera_rotations = [ - 30, -30, # Horizontal poloidal, - -90, # Vertical poloidal, - 0, # Tangential, - 0, # Combined poloidal/tangential, - ] - toroidal_camera_rotations = [ - 0, 0, # Horizontal poloidal - 0, # Vertical poloidal - -40, # Tangential - 40, # Combined poloidal/tangential - ] - radial_camera_rotations = [ - 0, 0, # Horizontal poloidal - 0, # Vertical poloidal - 90, # Tangential - 90, # Combined poloidal/tangential - ] - camera_origins = [ - Point3D(2.5, 0.05, 0), Point3D(2.5, -0.05, 0), # Horizontal poloidal - Point3D(1.3, 0, 1.42), # Vertical poloidal - Point3D(2.5, 0, 0), # Midplane tangential horizontal - Point3D(2.5, 0, -0.2), # Combined poloidal/tangential - ] - slit_sensor_separations = [ - 0.08, 0.08, # Horizontal poloidal - 0.05, # Vertical poloidal - 0.1, # Tangential - 0.1, # Combined poloidal/tangential - ] - all_sensor_angles = [ - [-22.5, -7.5, 7.5, 22.5], [-22.5, -7.5, 7.5, 22.5], # Horizontal poloidal - [-36, -12, 12, 36], # Vertical poloidal - [-18, -6, 6, 18], # Tangential - [-18, -6, 6, 18], # Combined poloidal/tangential - ] - toroidal_angles = [ - 10, 10, # Horizontal poloidal, need to avoid LFS limiters. - 0, # Vertical poloidal, happy to hit LFS limiters. - -15, # Tangential, avoid LFS limiters. - 15, # Combined poloidal/tangential, avoid LFS limiters. - ] - names = [ - 'HozPol1', 'HozPol2', - 'VertPol', - 'TanMid1', - 'TanPol1', - ] + camera_properties = { + 'HozPol1': {}, # Horizontal poloidal + 'HozPol2': {}, # Horizontal poloidal, + 'VertPol': {}, # Vertical poloidal + 'TanMid1': {}, # Tangential + 'TanPol1': {} # Combined poloidal/tangential + } + # poloidal rotations + camera_properties['HozPol1']['rotation_poloidal'] = 30 + camera_properties['HozPol2']['rotation_poloidal'] = -30 + camera_properties['VertPol']['rotation_poloidal'] = -90 + camera_properties['TanMid1']['rotation_poloidal'] = 0 + camera_properties['TanPol1']['rotation_poloidal'] = 0 + # toroidal rotation + camera_properties['HozPol1']['rotation_toroidal'] = 0 + camera_properties['HozPol2']['rotation_toroidal'] = 0 + camera_properties['VertPol']['rotation_toroidal'] = 0 + camera_properties['TanMid1']['rotation_toroidal'] = -40 + camera_properties['TanPol1']['rotation_toroidal'] = 40 + # radial rotation + camera_properties['HozPol1']['rotation_radial'] = 0 + camera_properties['HozPol2']['rotation_radial'] = 0 + camera_properties['VertPol']['rotation_radial'] = 0 + camera_properties['TanMid1']['rotation_radial'] = 90 + camera_properties['TanPol1']['rotation_radial'] = 90 + # origins + camera_properties['HozPol1']['origin'] = Point3D(2.5, 0.05, 0) + camera_properties['HozPol2']['origin'] = Point3D(2.5, -0.05, 0) + camera_properties['VertPol']['origin'] = Point3D(1.3, 0, 1.42) + camera_properties['TanMid1']['origin'] = Point3D(2.5, 0, 0) + camera_properties['TanPol1']['origin'] = Point3D(2.5, 0, -0.2) + # slit sensor separations + camera_properties['HozPol1']['slit_sensor_separation'] = 0.08 + camera_properties['HozPol2']['slit_sensor_separation'] = 0.08 + camera_properties['VertPol']['slit_sensor_separation'] = 0.05 + camera_properties['TanMid1']['slit_sensor_separation'] = 0.1 + camera_properties['TanPol1']['slit_sensor_separation'] = 0.1 + # sensor angles + camera_properties['HozPol1']['sensor_angles'] = [-22.5, -7.5, 7.5, 22.5] + camera_properties['HozPol2']['sensor_angles'] = [-22.5, -7.5, 7.5, 22.5] + camera_properties['VertPol']['sensor_angles'] = [-36, -12, 12, 36] + camera_properties['TanMid1']['sensor_angles'] = [-18, -6, 6, 18] + camera_properties['TanPol1']['sensor_angles'] = [-18, -6, 6, 18] + # toroidal angles + camera_properties['HozPol1']['toroidal_angle'] = 10 # need to avoid LFS limiters + camera_properties['HozPol2']['toroidal_angle'] = 10 # need to avoid LFS limiters + camera_properties['VertPol']['toroidal_angle'] = 0 # happy to hit LFS limiters + camera_properties['TanMid1']['toroidal_angle'] = -15 # avoid LFS limiters + camera_properties['TanPol1']['toroidal_angle'] = 15 # avoid LFS limiters + cameras = [] - # FIXME: this for loop definition is really ugly! - for ( - poloidal_rotation, toroidal_rotation, radial_rotation, camera_origin, - slit_sensor_separation, sensor_angles, toroidal_angle, name - ) in zip( - poloidal_camera_rotations, toroidal_camera_rotations, radial_camera_rotations, - camera_origins, - slit_sensor_separations, all_sensor_angles, toroidal_angles, names - ): - camera = _make_bolometer_camera(slit_sensor_separation, sensor_angles) + for name, prop in camera_properties.items(): + camera = _make_bolometer_camera(prop['slit_sensor_separation'], prop['sensor_angles']) # FIXME: work out how to combine tangential and poloidal rotations. camera.transform = ( - rotate_z(toroidal_angle) - * translate(camera_origin.x, camera_origin.y, camera_origin.z) - * rotate_z(toroidal_rotation) - * rotate_x(radial_rotation) - * rotate_y(poloidal_rotation + 90) + rotate_z(prop['toroidal_angle']) + * translate(prop['origin'].x, prop['origin'].y, prop['origin'].z) + * rotate_z(prop['rotation_toroidal']) + * rotate_x(prop['rotation_radial']) + * rotate_y(prop['rotation_poloidal'] + 90) * rotate_basis(-ZAXIS, YAXIS) ) camera.parent = parent - camera.name = "{} at angle {}".format(name, poloidal_rotation) + camera.name = "{} at angle {}".format(name, prop['rotation_poloidal']) cameras.append(camera) return cameras From e6d233d09be497c8372310de6f2221844926a4cc Mon Sep 17 00:00:00 2001 From: "Lovell, Jack J" Date: Wed, 27 May 2026 12:04:58 +0100 Subject: [PATCH 3/5] Remove unused import --- cherab/generomak/diagnostics/bolometers.py | 1 - 1 file changed, 1 deletion(-) diff --git a/cherab/generomak/diagnostics/bolometers.py b/cherab/generomak/diagnostics/bolometers.py index e64d9dec..a2412348 100644 --- a/cherab/generomak/diagnostics/bolometers.py +++ b/cherab/generomak/diagnostics/bolometers.py @@ -3,7 +3,6 @@ """ from raysect.core import (Node, Point3D, Vector3D, rotate_basis, rotate_x, rotate_y, rotate_z, translate) -from raysect.core.math import rotate from raysect.optical.material import AbsorbingSurface from raysect.primitive import Box, Subtract From c7b504f60bf64bb0de768ee11100ce195b52390d Mon Sep 17 00:00:00 2001 From: Jack Lovell Date: Wed, 17 Jun 2026 10:26:43 +0100 Subject: [PATCH 4/5] Improve Generomak bolometer geometry * Fix the series of transforms to position and orient the bolometer cameras. * Update geometry parameters such that the lines of sight all follow a consistant convention and specify the channel ordering in the docstring. This required adding the ability to rotate the sensors about the viewing axis to reverse the channel order of the horizontal-poloidal camera. * Correct the docstring of load_bolometers to specify the correct number of 4-channel sensors in each camera. * Expand docstrings and comments to better describe the intricacies of setting up the bolometer geometry. --- cherab/generomak/diagnostics/bolometers.py | 123 +++++++++++++++------ 1 file changed, 91 insertions(+), 32 deletions(-) diff --git a/cherab/generomak/diagnostics/bolometers.py b/cherab/generomak/diagnostics/bolometers.py index a2412348..e522fb36 100644 --- a/cherab/generomak/diagnostics/bolometers.py +++ b/cherab/generomak/diagnostics/bolometers.py @@ -16,7 +16,7 @@ ORIGIN = Point3D(0, 0, 0) # Bolometer geometry, independent of camera. BOX_WIDTH = 0.05 -BOX_WIDTH = 0.2 +BOX_WIDTH = 0.1 BOX_HEIGHT = 0.07 BOX_DEPTH = 0.2 THICKNESS = 1e-3 @@ -28,14 +28,23 @@ FOIL_SEPARATION = 0.00508 # 0.2 inch between foils -def _make_bolometer_camera(slit_sensor_separation, sensor_angles): +def _make_bolometer_camera(slit_sensor_separation, sensor_angles, sensor_rotations): """ Build a single bolometer camera. The camera consists of a box with a rectangular slit and 4 sensors, each of which has 4 foils. + In its local coordinate system, the camera's slit is located at the - origin and the sensors below the z=0 plane, looking up towards the slit. + origin with its width along the X axis and its height along the y + axis, and the sensors are below the z=0 plane looking up towards the + slit. + + The sensors are rotated by sensor_angles about the y axis to form a + fan, and by sensor_rotations about the axis defined by the line + between the slit and the sensor. A rotation of 180 degrees flips + the sensor upside down and therefore reverses the spatial ordering + of lines of sight relative to a rotation of 0 degrees. """ camera_box = Box(lower=Point3D(-BOX_WIDTH / 2, -BOX_HEIGHT / 2, -BOX_DEPTH), upper=Point3D(BOX_WIDTH / 2, BOX_HEIGHT / 2, 0)) @@ -55,10 +64,14 @@ def _make_bolometer_camera(slit_sensor_separation, sensor_angles): slit = BolometerSlit(slit_id="Example slit", centre_point=ORIGIN, basis_x=XAXIS, dx=SLIT_WIDTH, basis_y=YAXIS, dy=SLIT_HEIGHT, parent=bolometer_camera) - for j, angle in enumerate(sensor_angles): + for j, (angle, rotation) in enumerate(zip(sensor_angles, sensor_rotations)): # 4 bolometer foils, spaced at equal intervals along the local X axis - sensor = Node(name="Bolometer sensor", parent=bolometer_camera, - transform=rotate_y(angle) * translate(0, 0, -slit_sensor_separation)) + sensor = Node(name="Bolometer sensor", parent=bolometer_camera) + sensor.transform = ( + rotate_y(angle) + * rotate_z(rotation) + * translate(0, 0, -slit_sensor_separation) + ) for i, shift in enumerate([-1.5, -0.5, 0.5, 1.5]): # Note that the foils will be parented to the camera rather than the # sensor, so we need to define their transform relative to the camera. @@ -77,8 +90,8 @@ def load_bolometers(parent=None): """ Load the Generomak bolometers. - The Generomak bolometer diagnostic consists of multiple 12-channel - cameras. Each camera has 3 4-channel sensors inside. + The Generomak bolometer diagnostic consists of multiple 16-channel + cameras. Each camera has 4 4-channel sensors inside. * 2 cameras are located at the midplane with purely-poloidal, horizontal views. @@ -89,10 +102,40 @@ def load_bolometers(parent=None): curved lines of sight in the poloidal plane. It looks at the lower divertor. + Channel ordering is as follows: + * Poloidal channels are ordered anti-clockwise by line-of-sight: + channel 1 of HozPol1 views the top of the machine and channel 16 + HozPol2 views the bottom of the machine. Similarly, channel 1 of + VertPol views the high field side and channel 16 views the low + field side. + * Tangential channels are ordered by increasing tangency radius: + channel 1 of TanMid1 has its tangency radius on the high field + side and channel 16 has its tangency radius on the low field side. + * The combined tangential/poloidal channels follow both conventions: + channel 1 views the high field side and channel 16 views the low + field side. + :param parent: the scenegraph node the bolometers will belong to. :return: a list of BolometerCamera instances, one for each of the cameras described above. """ + # The coordinate system conventions are as follows. All angles are in + # degrees and increase clockwise when viewing along the relevant axes: + # y axis for poloidal rotation, z axis for toroidal rotation and x axis + # for radial rotation. + # - rotation_poloidal: viewing angle of the slit in the poloidal plane, + # with 0 being horizontally inwards. + # - rotation_toroidal: viewing angle of the slit in the toroidal plane, + # with 0 being purely radial. + # - rotation_radial: rotation about the radial axis, 0 being vertically upwards. + # - origin: position of the slit relative to the (x, z) poloidal plane i.e. y=0. + # - slit_sensor_separation: distance between slit and each 4-channel sensor. + # - sensor_angles: angle between slit normal and sensor normal. + # - sensor_rotations: rotation angle about the slit-sensor vector, enables + # reversing the order of lines of sight spatially within + # each sensor. + # - toroidal_angle: the angle of the poloidal plane in which the origin is + # definied, with 0 being the (x, z) plane. camera_properties = { 'HozPol1': {}, # Horizontal poloidal 'HozPol2': {}, # Horizontal poloidal, @@ -105,7 +148,7 @@ def load_bolometers(parent=None): camera_properties['HozPol2']['rotation_poloidal'] = -30 camera_properties['VertPol']['rotation_poloidal'] = -90 camera_properties['TanMid1']['rotation_poloidal'] = 0 - camera_properties['TanPol1']['rotation_poloidal'] = 0 + camera_properties['TanPol1']['rotation_poloidal'] = -25 # toroidal rotation camera_properties['HozPol1']['rotation_toroidal'] = 0 camera_properties['HozPol2']['rotation_toroidal'] = 0 @@ -113,30 +156,36 @@ def load_bolometers(parent=None): camera_properties['TanMid1']['rotation_toroidal'] = -40 camera_properties['TanPol1']['rotation_toroidal'] = 40 # radial rotation - camera_properties['HozPol1']['rotation_radial'] = 0 - camera_properties['HozPol2']['rotation_radial'] = 0 - camera_properties['VertPol']['rotation_radial'] = 0 - camera_properties['TanMid1']['rotation_radial'] = 90 - camera_properties['TanPol1']['rotation_radial'] = 90 - # origins - camera_properties['HozPol1']['origin'] = Point3D(2.5, 0.05, 0) - camera_properties['HozPol2']['origin'] = Point3D(2.5, -0.05, 0) + camera_properties['HozPol1']['rotation_radial'] = -90 + camera_properties['HozPol2']['rotation_radial'] = -90 + camera_properties['VertPol']['rotation_radial'] = -90 + camera_properties['TanMid1']['rotation_radial'] = 0 + camera_properties['TanPol1']['rotation_radial'] = 0 + # origins relative to the poloidal (x, z) plane + camera_properties['HozPol1']['origin'] = Point3D(2.45, 0.05, 0) + camera_properties['HozPol2']['origin'] = Point3D(2.45, -0.05, 0) camera_properties['VertPol']['origin'] = Point3D(1.3, 0, 1.42) camera_properties['TanMid1']['origin'] = Point3D(2.5, 0, 0) - camera_properties['TanPol1']['origin'] = Point3D(2.5, 0, -0.2) - # slit sensor separations + camera_properties['TanPol1']['origin'] = Point3D(2.2, 0, -0.8) + # slit-sensor separations camera_properties['HozPol1']['slit_sensor_separation'] = 0.08 camera_properties['HozPol2']['slit_sensor_separation'] = 0.08 camera_properties['VertPol']['slit_sensor_separation'] = 0.05 camera_properties['TanMid1']['slit_sensor_separation'] = 0.1 - camera_properties['TanPol1']['slit_sensor_separation'] = 0.1 - # sensor angles - camera_properties['HozPol1']['sensor_angles'] = [-22.5, -7.5, 7.5, 22.5] - camera_properties['HozPol2']['sensor_angles'] = [-22.5, -7.5, 7.5, 22.5] - camera_properties['VertPol']['sensor_angles'] = [-36, -12, 12, 36] - camera_properties['TanMid1']['sensor_angles'] = [-18, -6, 6, 18] - camera_properties['TanPol1']['sensor_angles'] = [-18, -6, 6, 18] - # toroidal angles + camera_properties['TanPol1']['slit_sensor_separation'] = 0.15 + # sensor angles relative to the slit + camera_properties['HozPol1']['sensor_angles'] = [22.5, 7.5, -7.5, -22.5] + camera_properties['HozPol2']['sensor_angles'] = [22.5, 7.5, -7.5, -22.5] + camera_properties['VertPol']['sensor_angles'] = [36, 12, -12, -36] + camera_properties['TanMid1']['sensor_angles'] = [18, 6, -6, -18] + camera_properties['TanPol1']['sensor_angles'] = [-12, -4, 4, 12] + # sensor rotation relative to the slit + camera_properties['HozPol1']['sensor_rotations'] = [0, 0, 0, 0] + camera_properties['HozPol2']['sensor_rotations'] = [0, 0, 0, 0] + camera_properties['VertPol']['sensor_rotations'] = [0, 0, 0, 0] + camera_properties['TanMid1']['sensor_rotations'] = [0, 0, 0, 0] + camera_properties['TanPol1']['sensor_rotations'] = [180, 180, 180, 180] + # toroidal angles about which to rotate the poloidal plane camera_properties['HozPol1']['toroidal_angle'] = 10 # need to avoid LFS limiters camera_properties['HozPol2']['toroidal_angle'] = 10 # need to avoid LFS limiters camera_properties['VertPol']['toroidal_angle'] = 0 # happy to hit LFS limiters @@ -145,17 +194,27 @@ def load_bolometers(parent=None): cameras = [] for name, prop in camera_properties.items(): - camera = _make_bolometer_camera(prop['slit_sensor_separation'], prop['sensor_angles']) - # FIXME: work out how to combine tangential and poloidal rotations. + camera = _make_bolometer_camera( + prop['slit_sensor_separation'], + prop['sensor_angles'], + prop['sensor_rotations'], + ) + # The transform is applied as follows: + # 1. Point the camera along the inward radial direction in the (x, z) plane. + # 2. Make the poloidal, radial and toroidal rotations while the camera is at + # the origin. + # 3. Move the camera to its position relative to the (x, z) plane. + # 4. Rotate the (x, z) plane to the correct toroidal angle. + # Transforms are applied right-to-left (or bottom-to-top with one per line): camera.transform = ( rotate_z(prop['toroidal_angle']) * translate(prop['origin'].x, prop['origin'].y, prop['origin'].z) * rotate_z(prop['rotation_toroidal']) + * rotate_y(prop['rotation_poloidal']) * rotate_x(prop['rotation_radial']) - * rotate_y(prop['rotation_poloidal'] + 90) - * rotate_basis(-ZAXIS, YAXIS) + * rotate_basis(-XAXIS, ZAXIS) ) camera.parent = parent - camera.name = "{} at angle {}".format(name, prop['rotation_poloidal']) + camera.name = name cameras.append(camera) return cameras From 2d838a6f13a1b31b88c7b8d42c86683fa8ba8ff4 Mon Sep 17 00:00:00 2001 From: Jack Lovell Date: Wed, 17 Jun 2026 10:36:00 +0100 Subject: [PATCH 5/5] Update transform comment with more logical order --- cherab/generomak/diagnostics/bolometers.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cherab/generomak/diagnostics/bolometers.py b/cherab/generomak/diagnostics/bolometers.py index e522fb36..05e2115e 100644 --- a/cherab/generomak/diagnostics/bolometers.py +++ b/cherab/generomak/diagnostics/bolometers.py @@ -201,7 +201,7 @@ def load_bolometers(parent=None): ) # The transform is applied as follows: # 1. Point the camera along the inward radial direction in the (x, z) plane. - # 2. Make the poloidal, radial and toroidal rotations while the camera is at + # 2. Make the radial, poloidal and toroidal rotations while the camera is at # the origin. # 3. Move the camera to its position relative to the (x, z) plane. # 4. Rotate the (x, z) plane to the correct toroidal angle.