cherab interface: Calculate heat fluxes to walls

`xbout.wall.AxisymmetricWall` represents an axisymmetric wall

`xbout.wall.read_geqdsk` reads a GEQDSK file, returning an AxisymmetricWall

`xbout.cherab.triangulate.TriangularData.wall_flux`
takes an AxisymmetricWall and calculates heat fluxes to each element.

So to calculate heat fluxes:

```
wall = xbout.wall.read_geqdsk("geqdsk")

ds = ds.bout.with_cherab_grid()

da = -bd['Rd+_ex'].isel(t=1, zeta=0)

da[:2,:] = 0.0
da[-2:,:] = 0.0

data = da.bout.as_cherab_data()

result = data.wall_flux(cat_wall, pixel_samples=1000)
```

xbout/__init__.py

@@ -1,6 +1,6 @@
 from .load import open_boutdataset, collect
 
-from . import geometries
+from . import geometries, wall
 from .geometries import register_geometry, REGISTERED_GEOMETRIES
 
 from .boutdataset import BoutDatasetAccessor

xbout/cherab/triangulate.py

@@ -7,7 +7,6 @@
 """
 
 import numpy as np
-import xarray as xr
 
 
 class TriangularData:
@@ -74,7 +73,7 @@ def to_emitter(
             material=emitting_material,
         )
 
-    def plot_2d(self, ax=None, nr: int = 150, nz: int = 150):
+    def plot_2d(self, ax=None, nr: int = 150, nz: int = 150, colorbar: bool = True):
         """
         Make a 2D plot of the data
 
@@ -95,12 +94,136 @@ def plot_2d(self, ax=None, nr: int = 150, nz: int = 150):
         image = ax.imshow(
             emiss_sampled.T, origin="lower", extent=(r.min(), r.max(), z.min(), z.max())
         )
-        fig.colorbar(image)
+        if colorbar:
+            fig.colorbar(image)
         ax.set_xlabel("r")
         ax.set_ylabel("z")
 
         return ax
 
+    def wall_flux(
+        self,
+        wall,
+        pixel_samples: int = 2000,
+        wall_detector_offset: float = 0.001,
+        toroidal_width: float = 0.01,
+        step: float = 0.01,
+    ):
+        """
+        Calculate power onto segments of an axisymmetric wall
+
+        Based on the Cherab manual here:
+        https://www.cherab.info/demonstrations/radiation_loads/symmetric_power_load.html#symmetric-power-load
+        
+        Parameters
+        ----------
+
+        wall : AxisymmetricWall
+            Defines a closed loop in R-Z, that defines an axisymmetric wall.
+
+        pixel_samples : The number of samples per wall segment
+
+        wall_detector_offset
+            Distance of detector pixels from wall [meters].
+            Slightly displaced to avoid numerical overlap (ray trapping)
+
+        toroidal_width
+            toroidal width of the detectors [meters]
+
+        step
+            Ray path step length [meters]
+
+        """
+        from raysect.core import Point3D, Vector3D, translate, rotate_basis
+        from raysect.optical import World
+        from raysect.optical.observer import PowerPipeline0D
+        from raysect.optical.material import AbsorbingSurface
+        from raysect.optical.observer.nonimaging.pixel import Pixel
+        from cherab.tools.primitives import axisymmetric_mesh_from_polygon
+
+        world = World()
+
+        # Create a plasma emitter
+        emitter = self.to_emitter(parent=world, step=step)
+
+        # Create the walls around the plasma
+        wall_polygon = wall.to_polygon()  # [npoints, 2] array
+
+        # create a 3D mesh from the 2D polygon outline using symmetry
+        wall_mesh = axisymmetric_mesh_from_polygon(wall_polygon)
+        wall_mesh.parent = world
+        wall_mesh.material = AbsorbingSurface()
+
+        result = []
+        for rz1, rz2 in wall:
+            p1 = Point3D(rz1[0], 0, rz1[1])
+            p2 = Point3D(rz2[0], 0, rz2[1])
+
+            # evaluate y_vector
+            y_vector_full = p1.vector_to(p2)
+
+            if y_vector_full.length < 1e-5:
+                continue  # Skip; points p1, p2 are the same
+
+            y_vector = y_vector_full.normalise()
+            y_width = y_vector_full.length
+
+            # evaluate normal_vector (inward pointing)
+            normal_vector = y_vector.cross(Vector3D(0, 1, 0)).normalise()
+
+            # evaluate the central point of the detector
+            # Note: Displaced in the normal direction by wall_detector_offset
+            detector_center = (
+                p1 + y_vector_full * 0.5 + wall_detector_offset * normal_vector
+            )
+
+            # Use the power pipeline to record total power arriving at the surface
+            power_data = PowerPipeline0D()
+
+            # Note: Displace the pixel location
+            pixel_transform = translate(
+                detector_center.x, detector_center.y, detector_center.z
+            ) * rotate_basis(normal_vector, y_vector)
+
+            # Use pixel_samples argument to increase amount of sampling and reduce noise
+            pixel = Pixel(
+                [power_data],
+                x_width=toroidal_width,
+                y_width=y_width,
+                name=f"wall-{rz1}-{rz2}",
+                spectral_bins=1,
+                transform=pixel_transform,
+                parent=world,
+                pixel_samples=pixel_samples,
+            )
+
+            pixel_area = toroidal_width * y_width
+
+            # Start collecting samples
+            pixel.observe()
+
+            # Estimate power density
+            power_density = power_data.value.mean / pixel_area
+
+            # For checking energy conservation.
+            # Revolve this tile around the CYLINDRICAL z-axis to get total power collected by these tiles.
+            # Add up all the tile contributions to get total power collected.
+            detector_radius = np.sqrt(detector_center.x ** 2 + detector_center.y ** 2)
+            observed_total_power = power_density * (
+                y_width * 2 * np.pi * detector_radius
+            )
+
+            result.append(
+                {
+                    "rz1": rz1,
+                    "rz2": rz2,
+                    "power_density": power_density,
+                    "total_power": observed_total_power,
+                }
+            )
+
+        return result
+
 
 class Triangulate:
     """

xbout/wall.py

@@ -0,0 +1,116 @@
+"""
+Routines to read and represent wall geometries
+"""
+
+import numpy as np
+
+class AxisymmetricWall:
+    def __init__(self, Rs, Zs):
+        """
+        Defines a 2D (R,Z) axisymmetric wall
+
+        Parameters
+        ----------
+
+        Rs : list or 1D array
+            Major radius coordinates [meters]
+        Zs : list or 1D array
+            Vertical coordinates [meters]
+
+        """
+        if len(Rs) != len(Zs):
+            raise ValueError("Rs and Zs arrays have different lengths")
+
+        # Ensure that the members are numpy arrays
+        self.Rs = np.array(Rs)
+        self.Zs = np.array(Zs)
+
+    def __iter__(self):
+        """
+        Iterate over wall elements
+
+        Each iteration returns a pair of (R,Z) pairs:
+            ((R1, Z1), (R2, Z2))
+
+        These pairs define wall segment.
+        """
+        return iter(zip(zip(self.Rs, self.Zs),
+                        zip(np.roll(self.Rs, -1), np.roll(self.Zs, -1))))
+
+    def to_polygon(self):
+        """
+        Returns a 2D Numpy array [npoints, 2]
+        Index 0 is major radius (R) in meters
+        Index 1 is height (Z) in meters
+        """
+        return np.stack((self.Rs, self.Zs), axis=-1)
+
+    def plot(self, linestyle='k-', ax = None):
+        """
+        Plot the wall on given axis. If no axis
+        is given then a new figure is created.
+
+        Returns
+        -------
+
+        The matplotlib axis containing the plot
+        """
+
+        import matplotlib.pyplot as plt
+
+        if ax is None:
+            fig, ax = plt.subplots()
+
+        ax.plot(self.Rs, self.Zs, linestyle)
+        return ax
+
+def read_geqdsk(filehandle):
+    """
+    Read wall geometry from a GEQDSK file.
+
+    Note: Requires the freeqdsk package
+    """
+
+    if isinstance(filehandle, str):
+        with open(filehandle, "r") as f:
+            return read_geqdsk(f)
+
+    from freeqdsk import geqdsk
+
+    data = geqdsk.read(filehandle)
+    # rlim and zlim should be 1D arrays of wall coordinates
+
+    if not (hasattr(data, "rlim") and hasattr(data, "zlim")):
+        raise ValueError(f"Wall coordinates not found in GEQDSK file")
+
+    return AxisymmetricWall(data["rlim"], data["zlim"])
+
+
+def read_csv(filehandle, delimiter=','):
+    """
+    Parameters
+    ----------
+
+    filehandle: File handle
+        Must contain two columns, for R and Z coordinates [meters]
+    
+    delimier : character
+        A single character that separates fields
+
+    Notes:
+    - Uses the python `csv` module 
+    """
+    import csv
+    reader = csv.reader(filehandle, delimiter=delimiter)
+    Rs = []
+    Zs = []
+    for row in reader:
+        if len(row) == 0:
+            continue # Skip empty rows
+        if len(row) != 2:
+            raise ValueError(f"CSV row should contain two columns: {row}")
+        Rs.append(float(row[0]))
+        Zs.append(float(row[1]))
+
+    return AxisymmetricWall(Rs, Zs)
+