diff --git a/src/struphy/geometry/domains.py b/src/struphy/geometry/domains.py
index 48473dd23..52cce7908 100644
--- a/src/struphy/geometry/domains.py
+++ b/src/struphy/geometry/domains.py
@@ -565,6 +565,56 @@ def __init__(
         super().__init__()
 
 
+class MagnetotailSlab(Domain):
+    r"""Elongated slab with a pinched center.
+
+    This mapping is intended as a simple geometry for studying reconnection in
+    the Earth's night-side magnetotail. The domain is a long slab in the
+    tail-aligned direction with a smooth reduction of the cross-tail width near
+    the center plane.
+
+    Parameters
+    ----------
+    Lx : float
+        Length in the tail-aligned direction (default: 12.0).
+    Ly : float
+        Full width of the slab in the cross-tail direction away from the pinch (default: 4.0).
+    Lz : float
+        Thickness in the normal direction (default: 2.0).
+    pinch : float
+        Relative pinch strength in [0, 1). Larger values give a narrower center (default: 0.6).
+    pinch_width : float
+        Characteristic width of the pinched region in physical x-units (default: 2.0).
+    x_center : float
+        Center position of the pinch in physical x-coordinate (default: 0.0).
+    """
+
+    def __init__(
+        self,
+        Lx: float = 12.0,
+        Ly: float = 4.0,
+        Lz: float = 2.0,
+        pinch: float = 0.6,
+        pinch_width: float = 2.0,
+        x_center: float = 0.0,
+    ):
+        self.kind_map = 13
+
+        self.params = copy.deepcopy(locals())
+        self.params_numpy = self.get_params_numpy()
+
+        assert Lx > 0.0, f"Need positive slab length, got {Lx =}"
+        assert Ly > 0.0, f"Need positive slab width, got {Ly =}"
+        assert Lz > 0.0, f"Need positive slab thickness, got {Lz =}"
+        assert 0.0 <= pinch < 1.0, f"Pinch strength must satisfy 0 <= pinch < 1, got {pinch =}"
+        assert pinch_width > 0.0, f"Need positive pinch width, got {pinch_width =}"
+
+        self.periodic_eta3 = False
+        self.pole = False
+
+        super().__init__()
+
+
 class HollowCylinder(Domain):
     r""" Cylinder with possible hole around the axis.
 
@@ -794,6 +844,219 @@ def inverse_map(self, x, y, z, bounded=True, change_out_order=False):
             return eta1, eta2, eta3
 
 
+class GeospaceFluxDomain(Domain):
+    r"""Global magnetospheric flux-surface inspired mapping.
+
+    This domain maps the logical unit cube onto an Earth-centered, solar-wind
+    distorted geometry with compressed dayside and elongated nightside tail.
+    It provides a single smooth coordinate mapping useful for idealized studies
+    of bow-shock/magnetosheath/magnetotail coupling and tail reconnection setup.
+
+    Coordinates
+    -----------
+    * :math:`\eta_1`: radial/normal-like coordinate from ionosphere to bow shock
+    * :math:`\eta_2`: dayside-to-nightside poloidal angle coordinate
+    * :math:`\eta_3`: clock angle around Earth-Sun axis
+
+    Parameters
+    ----------
+    r_iono : float
+        Inner reference radius (ionosphere-like boundary).
+    r_mp_dayside : float
+        Magnetopause radius on dayside.
+    r_mp_tail : float
+        Magnetopause radius in the nightside tail.
+    r_bs_dayside : float
+        Bow-shock radius on dayside.
+    r_bs_tail : float
+        Bow-shock radius in the nightside tail.
+    mp_eta1 : float
+        Logical location of magnetopause in :math:`\eta_1` (0 < mp_eta1 < 1).
+    sheet_flatten : float
+        Tail current-sheet flattening strength in [0, 1).
+    """
+
+    def __init__(
+        self,
+        r_iono: float = 1.0,
+        r_mp_dayside: float = 8.0,
+        r_mp_tail: float = 30.0,
+        r_bs_dayside: float = 12.0,
+        r_bs_tail: float = 45.0,
+        mp_eta1: float = 0.72,
+        sheet_flatten: float = 0.45,
+    ):
+        self.kind_map = 24
+
+        self.params = copy.deepcopy(locals())
+        self.params_numpy = self.get_params_numpy()
+
+        assert r_iono > 0.0, f"Need positive ionosphere radius, got {r_iono =}"
+        assert r_mp_dayside > r_iono, f"Need r_mp_dayside > r_iono, got {r_mp_dayside =}, {r_iono =}"
+        assert r_mp_tail > r_iono, f"Need r_mp_tail > r_iono, got {r_mp_tail =}, {r_iono =}"
+        assert r_bs_dayside > r_mp_dayside, f"Need r_bs_dayside > r_mp_dayside, got {r_bs_dayside =}, {r_mp_dayside =}"
+        assert r_bs_tail > r_mp_tail, f"Need r_bs_tail > r_mp_tail, got {r_bs_tail =}, {r_mp_tail =}"
+        assert 0.0 < mp_eta1 < 1.0, f"Need 0 < mp_eta1 < 1, got {mp_eta1 =}"
+        assert 0.0 <= sheet_flatten < 1.0, f"Need 0 <= sheet_flatten < 1, got {sheet_flatten =}"
+
+        self.periodic_eta3 = True
+        self.pole = False
+
+        super().__init__()
+
+
+class WarpedAccretionDisk(Domain):
+    r"""Warped accretion disk mapping.
+
+    This mapping describes a cylindrical disk with finite thickness and a smooth,
+    radius-dependent vertical warp of the disk midplane.
+
+    Coordinates
+    -----------
+    * :math:`\eta_1`: radial coordinate from inner to outer disk radius
+    * :math:`\eta_2`: azimuthal coordinate around the central object
+    * :math:`\eta_3`: vertical coordinate across disk thickness
+
+    Parameters
+    ----------
+    r_in : float
+        Inner disk radius.
+    r_out : float
+        Outer disk radius.
+    thickness : float
+        Half-thickness scaling in the vertical direction.
+    warp_amp : float
+        Warp amplitude factor.
+    warp_power : float
+        Power-law exponent for radial growth of the warp.
+    node_angle : float
+        Azimuthal node angle of the warp in radians.
+    tor_period : int
+        Azimuthal periodicity: :math:`\phi = 2\pi\eta_2/\mathrm{tor\_period}`.
+    """
+
+    def __init__(
+        self,
+        r_in: float = 2.0,
+        r_out: float = 12.0,
+        thickness: float = 0.3,
+        warp_amp: float = 0.15,
+        warp_power: float = 1.5,
+        node_angle: float = 0.0,
+        tor_period: int = 1,
+    ):
+        self.kind_map = 25
+
+        self.params = copy.deepcopy(locals())
+        self.params_numpy = self.get_params_numpy()
+
+        assert r_in > 0.0, f"Need positive inner radius, got {r_in =}"
+        assert r_out > r_in, f"Need r_out > r_in, got {r_out =}, {r_in =}"
+        assert thickness > 0.0, f"Need positive thickness, got {thickness =}"
+        assert warp_amp >= 0.0, f"Need non-negative warp amplitude, got {warp_amp =}"
+        assert warp_power >= 0.0, f"Need non-negative warp power, got {warp_power =}"
+        assert tor_period > 0, f"Need positive toroidal periodicity, got {tor_period =}"
+
+        self.periodic_eta3 = True
+        self.pole = False
+
+        super().__init__()
+
+
+class Spheromak(Domain):
+    r"""Compact toroidal plasma proxy in a spherical-like shell.
+
+    This mapping provides a simple geometry for spheromak studies using nested
+    closed flux-surface-like shells with optional vertical elongation.
+
+    Coordinates
+    -----------
+    * :math:`\eta_1`: radial shell coordinate
+    * :math:`\eta_2`: poloidal angle coordinate
+    * :math:`\eta_3`: toroidal/azimuthal angle coordinate
+
+    Parameters
+    ----------
+    r0 : float
+        Inner radius.
+    a : float
+        Plasma minor-size scale (outer radius is r0 + a).
+    kappa : float
+        Vertical elongation factor.
+    tor_period : int
+        Azimuthal periodicity: :math:`\phi=2\pi\eta_3/\mathrm{tor\_period}`.
+    """
+
+    def __init__(
+        self,
+        r0: float = 0.0,
+        a: float = 1.0,
+        kappa: float = 1.0,
+        tor_period: int = 1,
+    ):
+        self.kind_map = 26
+
+        self.params = copy.deepcopy(locals())
+        self.params_numpy = self.get_params_numpy()
+
+        assert r0 >= 0.0, f"Need non-negative inner radius, got {r0 =}"
+        assert a > 0.0, f"Need positive minor size a, got {a =}"
+        assert kappa > 0.0, f"Need positive elongation kappa, got {kappa =}"
+        assert tor_period > 0, f"Need positive toroidal periodicity, got {tor_period =}"
+
+        self.periodic_eta3 = True
+        self.pole = r0 == 0.0
+
+        super().__init__()
+
+
+class HallEffectThrusterChannel(Domain):
+    r"""Coaxial annular Hall thruster channel with axial end packing.
+
+    The mapping models an annular channel where :math:`\eta_3` is the axial
+    coordinate (anode to exit), with increased point packing near both ends
+    controlled by a smooth high-frequency modulation.
+
+    Parameters
+    ----------
+    r_in : float
+        Inner channel radius.
+    r_out : float
+        Outer channel radius.
+    length : float
+        Axial channel length.
+    pack_strength : float
+        End-packing strength in [0, 1). Higher values cluster more points at
+        the anode and exit.
+    tor_period : int
+        Azimuthal periodicity: :math:`\phi=2\pi\eta_2/\mathrm{tor\_period}`.
+    """
+
+    def __init__(
+        self,
+        r_in: float = 1.0,
+        r_out: float = 2.0,
+        length: float = 5.0,
+        pack_strength: float = 0.7,
+        tor_period: int = 1,
+    ):
+        self.kind_map = 27
+
+        self.params = copy.deepcopy(locals())
+        self.params_numpy = self.get_params_numpy()
+
+        assert r_in > 0.0, f"Need positive inner radius, got {r_in =}"
+        assert r_out > r_in, f"Need r_out > r_in, got {r_out =}, {r_in =}"
+        assert length > 0.0, f"Need positive channel length, got {length =}"
+        assert 0.0 <= pack_strength < 1.0, f"Need 0 <= pack_strength < 1, got {pack_strength =}"
+        assert tor_period > 0, f"Need positive toroidal periodicity, got {tor_period =}"
+
+        self.periodic_eta3 = False
+        self.pole = False
+
+        super().__init__()
+
+
 class ShafranovShiftCylinder(Domain):
     r""" Cylinder with quadratic Shafranov shift.
 
@@ -972,3 +1235,121 @@ def __init__(
         self.pole = True
 
         super().__init__()
+
+class MoebiusStrip(Domain):
+    r"""Thickened Moebius strip domain.
+    
+    A toroidal ribbon that performs a 180-degree twist over one revolution.
+    This serves as a test for non-standard periodicity and metric tensor shifts.
+
+    Parameters
+    ----------
+    R : float
+        Major radius of the loop (default: 3.0).
+    width : float
+        Width of the ribbon (default: 1.0).
+    thickness : float
+        Thickness of the ribbon (default: 0.1).
+    """
+
+    def __init__(
+        self,
+        R: float = 3.0,
+        width: float = 1.0,
+        thickness: float = 0.1,
+    ):
+        self.kind_map = 50  # Unique ID for the new domain
+
+        # use params setter
+        self.params = copy.deepcopy(locals())
+        self.params_numpy = self.get_params_numpy()
+
+        # Periodicity logic: 
+        # Technically, eta3 is periodic, but the boundary at eta3=1 
+        # maps to a flipped version of eta3=0.
+        self.periodic_eta3 = True
+        self.pole = False
+
+        super().__init__()
+
+    def map(self, eta1, eta2, eta3):
+        """Analytical mapping for MoebiusStrip"""
+        
+        # Center the coordinates for width and thickness
+        w_hat = self.params['width'] * (eta1 - 0.5)
+        t_hat = self.params['thickness'] * (eta2 - 0.5)
+        
+        # Half-twist angle
+        alpha = xp.pi * eta3
+        # Standard toroidal angle
+        phi = 2 * xp.pi * eta3
+
+        # Effective radius in the XY plane
+        r_eff = self.params['R'] + w_hat * xp.cos(alpha) - t_hat * xp.sin(alpha)
+
+        x = r_eff * xp.cos(phi)
+        y = r_eff * xp.sin(phi)
+        z = w_hat * xp.sin(alpha) + t_hat * xp.cos(alpha)
+
+        return x, y, z
+
+
+class DiagnosticPortHoleTorus(Domain):
+    r"""Torus with a localized diagnostic port deformation.
+
+    This mapping models a toroidal vessel with a smooth, localized outward bulge
+    attached to the torus tube. The bulge approximates the effect of a small
+    radial/poloidal cylindrical diagnostic port while remaining a single smooth
+    patch, making it suitable for studying perturbations near reactor openings.
+
+    Parameters
+    ----------
+    a1 : float
+        Inner minor radius of the torus tube (default: 0.1).
+    a2 : float
+        Outer minor radius of the torus tube (default: 1.0).
+    R0 : float
+        Major radius of the torus (default: 3.0).
+    tor_period : int
+        Toroidal periodicity built into the mapping: :math:`\phi=2\pi\,\eta_3/\mathrm{torperiod}` (default: 1).
+    port_depth : float
+        Maximum outward radial deformation of the localized port (default: 0.25).
+    port_eta2_center : float
+        Center of the port in logical poloidal coordinate :math:`\eta_2` (default: 0.0).
+    port_eta3_center : float
+        Center of the port in logical toroidal coordinate :math:`\eta_3` (default: 0.0).
+    port_eta2_width : float
+        Width parameter of the port in logical poloidal coordinate (default: 0.08).
+    port_eta3_width : float
+        Width parameter of the port in logical toroidal coordinate (default: 0.08).
+    """
+
+    def __init__(
+        self,
+        a1: float = 0.1,
+        a2: float = 1.0,
+        R0: float = 3.0,
+        tor_period: int = 1,
+        port_depth: float = 0.25,
+        port_eta2_center: float = 0.0,
+        port_eta3_center: float = 0.0,
+        port_eta2_width: float = 0.08,
+        port_eta3_width: float = 0.08,
+    ):
+        self.kind_map = 23
+
+        # use params setter
+        self.params = copy.deepcopy(locals())
+        self.params_numpy = self.get_params_numpy()
+
+        assert a2 <= R0, f"The minor radius must be smaller or equal than the major radius! {a2 =}, {R0 =}"
+        assert a2 + port_depth <= R0, (
+            f"The localized port deformation must keep the torus non-self-intersecting! {a2 =}, {port_depth =}, {R0 =}"
+        )
+        assert port_eta2_width > 0.0, f"Need positive port width in eta2, got {port_eta2_width =}"
+        assert port_eta3_width > 0.0, f"Need positive port width in eta3, got {port_eta3_width =}"
+
+        self.periodic_eta3 = True
+        self.pole = a1 == 0.0
+
+        super().__init__()
\ No newline at end of file
diff --git a/src/struphy/geometry/evaluation_kernels.py b/src/struphy/geometry/evaluation_kernels.py
index 4f97b9ce9..a64d411e4 100644
--- a/src/struphy/geometry/evaluation_kernels.py
+++ b/src/struphy/geometry/evaluation_kernels.py
@@ -106,6 +106,19 @@ def f(
             args.params[3],
             f_out,
         )
+    elif args.kind_map == 13:
+        mappings_kernels.magnetotail_slab(
+            eta1,
+            eta2,
+            eta3,
+            args.params[0],
+            args.params[1],
+            args.params[2],
+            args.params[3],
+            args.params[4],
+            args.params[5],
+            f_out,
+        )
     elif args.kind_map == 20:
         mappings_kernels.hollow_cyl(
             eta1,
@@ -141,6 +154,73 @@ def f(
             args.params[5],
             f_out,
         )
+    elif args.kind_map == 23:
+        mappings_kernels.diagnostic_port_hole_torus(
+            eta1,
+            eta2,
+            eta3,
+            args.params[0],
+            args.params[1],
+            args.params[2],
+            args.params[3],
+            args.params[4],
+            args.params[5],
+            args.params[6],
+            args.params[7],
+            args.params[8],
+            f_out,
+        )
+    elif args.kind_map == 24:
+        mappings_kernels.geospace_flux_domain(
+            eta1,
+            eta2,
+            eta3,
+            args.params[0],
+            args.params[1],
+            args.params[2],
+            args.params[3],
+            args.params[4],
+            args.params[5],
+            args.params[6],
+            f_out,
+        )
+    elif args.kind_map == 25:
+        mappings_kernels.warped_accretion_disk(
+            eta1,
+            eta2,
+            eta3,
+            args.params[0],
+            args.params[1],
+            args.params[2],
+            args.params[3],
+            args.params[4],
+            args.params[5],
+            args.params[6],
+            f_out,
+        )
+    elif args.kind_map == 26:
+        mappings_kernels.spheromak(
+            eta1,
+            eta2,
+            eta3,
+            args.params[0],
+            args.params[1],
+            args.params[2],
+            args.params[3],
+            f_out,
+        )
+    elif args.kind_map == 27:
+        mappings_kernels.hall_effect_thruster_channel(
+            eta1,
+            eta2,
+            eta3,
+            args.params[0],
+            args.params[1],
+            args.params[2],
+            args.params[3],
+            args.params[4],
+            f_out,
+        )
     elif args.kind_map == 30:
         mappings_kernels.shafranov_shift(
             eta1,
@@ -177,6 +257,16 @@ def f(
             args.params[6],
             f_out,
         )
+    elif args.kind_map == 50:
+        mappings_kernels.moebius_strip(
+            eta1,
+            eta2,
+            eta3,
+            args.params[0],
+            args.params[1],
+            args.params[2],
+            f_out,
+        )
 
 
 def df(
@@ -265,6 +355,19 @@ def df(
             args.params[3],
             df_out,
         )
+    elif args.kind_map == 13:
+        mappings_kernels.magnetotail_slab_df(
+            eta1,
+            eta2,
+            eta3,
+            args.params[0],
+            args.params[1],
+            args.params[2],
+            args.params[3],
+            args.params[4],
+            args.params[5],
+            df_out,
+        )
     elif args.kind_map == 20:
         mappings_kernels.hollow_cyl_df(
             eta1,
@@ -299,6 +402,73 @@ def df(
             args.params[5],
             df_out,
         )
+    elif args.kind_map == 23:
+        mappings_kernels.diagnostic_port_hole_torus_df(
+            eta1,
+            eta2,
+            eta3,
+            args.params[0],
+            args.params[1],
+            args.params[2],
+            args.params[3],
+            args.params[4],
+            args.params[5],
+            args.params[6],
+            args.params[7],
+            args.params[8],
+            df_out,
+        )
+    elif args.kind_map == 24:
+        mappings_kernels.geospace_flux_domain_df(
+            eta1,
+            eta2,
+            eta3,
+            args.params[0],
+            args.params[1],
+            args.params[2],
+            args.params[3],
+            args.params[4],
+            args.params[5],
+            args.params[6],
+            df_out,
+        )
+    elif args.kind_map == 25:
+        mappings_kernels.warped_accretion_disk_df(
+            eta1,
+            eta2,
+            eta3,
+            args.params[0],
+            args.params[1],
+            args.params[2],
+            args.params[3],
+            args.params[4],
+            args.params[5],
+            args.params[6],
+            df_out,
+        )
+    elif args.kind_map == 26:
+        mappings_kernels.spheromak_df(
+            eta1,
+            eta2,
+            eta3,
+            args.params[0],
+            args.params[1],
+            args.params[2],
+            args.params[3],
+            df_out,
+        )
+    elif args.kind_map == 27:
+        mappings_kernels.hall_effect_thruster_channel_df(
+            eta1,
+            eta2,
+            eta3,
+            args.params[0],
+            args.params[1],
+            args.params[2],
+            args.params[3],
+            args.params[4],
+            df_out,
+        )
     elif args.kind_map == 30:
         mappings_kernels.shafranov_shift_df(
             eta1,
@@ -335,6 +505,16 @@ def df(
             args.params[6],
             df_out,
         )
+    elif args.kind_map == 50:
+        mappings_kernels.moebius_strip_df(
+            eta1,
+            eta2,
+            eta3,
+            args.params[0],
+            args.params[1],
+            args.params[2],
+            df_out,
+        )
 
 
 def det_df(
@@ -435,11 +615,22 @@ def df_inv(
             dfinv_out[1, 2] = 0.0
             dfinv_out[2, 0] = 0.0
             dfinv_out[2, 1] = 0.0
+        elif args.kind_map == 13:
+            dfinv_out[0, 1] = 0.0
+            dfinv_out[0, 2] = 0.0
+            dfinv_out[1, 2] = 0.0
+            dfinv_out[2, 0] = 0.0
+            dfinv_out[2, 1] = 0.0
         elif args.kind_map == 20:
             dfinv_out[0, 2] = 0.0
             dfinv_out[1, 2] = 0.0
             dfinv_out[2, 0] = 0.0
             dfinv_out[2, 1] = 0.0
+        elif args.kind_map == 27:
+            dfinv_out[0, 2] = 0.0
+            dfinv_out[1, 2] = 0.0
+            dfinv_out[2, 0] = 0.0
+            dfinv_out[2, 1] = 0.0
         elif args.kind_map == 21:
             dfinv_out[0, 2] = 0.0
             dfinv_out[1, 2] = 0.0
@@ -534,6 +725,11 @@ def g(
             g_out[1, 2] = 0.0
             g_out[2, 0] = 0.0
             g_out[2, 1] = 0.0
+        elif args.kind_map == 13:
+            g_out[0, 2] = 0.0
+            g_out[1, 2] = 0.0
+            g_out[2, 0] = 0.0
+            g_out[2, 1] = 0.0
         elif args.kind_map == 20:
             g_out[0, 1] = 0.0
             g_out[0, 2] = 0.0
@@ -541,6 +737,13 @@ def g(
             g_out[1, 2] = 0.0
             g_out[2, 0] = 0.0
             g_out[2, 1] = 0.0
+        elif args.kind_map == 27:
+            g_out[0, 1] = 0.0
+            g_out[0, 2] = 0.0
+            g_out[1, 0] = 0.0
+            g_out[1, 2] = 0.0
+            g_out[2, 0] = 0.0
+            g_out[2, 1] = 0.0
         elif args.kind_map == 21:
             g_out[0, 2] = 0.0
             g_out[1, 2] = 0.0
@@ -653,6 +856,11 @@ def g_inv(
             ginv_out[1, 2] = 0.0
             ginv_out[2, 0] = 0.0
             ginv_out[2, 1] = 0.0
+        elif args.kind_map == 13:
+            ginv_out[0, 2] = 0.0
+            ginv_out[1, 2] = 0.0
+            ginv_out[2, 0] = 0.0
+            ginv_out[2, 1] = 0.0
         elif args.kind_map == 20:
             ginv_out[0, 1] = 0.0
             ginv_out[0, 2] = 0.0
@@ -660,6 +868,13 @@ def g_inv(
             ginv_out[1, 2] = 0.0
             ginv_out[2, 0] = 0.0
             ginv_out[2, 1] = 0.0
+        elif args.kind_map == 27:
+            ginv_out[0, 1] = 0.0
+            ginv_out[0, 2] = 0.0
+            ginv_out[1, 0] = 0.0
+            ginv_out[1, 2] = 0.0
+            ginv_out[2, 0] = 0.0
+            ginv_out[2, 1] = 0.0
         elif args.kind_map == 21:
             ginv_out[0, 2] = 0.0
             ginv_out[1, 2] = 0.0
diff --git a/src/struphy/geometry/mappings_kernels.py b/src/struphy/geometry/mappings_kernels.py
index c67771dc0..e7b11f797 100644
--- a/src/struphy/geometry/mappings_kernels.py
+++ b/src/struphy/geometry/mappings_kernels.py
@@ -1,4 +1,4 @@
-from numpy import arcsin, arctan, cos, pi, sin, sqrt, tan, zeros
+from numpy import arcsin, arctan, cos, exp, pi, sin, sqrt, tan, zeros
 from pyccel.decorators import pure, stack_array
 
 import struphy.bsplines.bsplines_kernels as bsplines_kernels
@@ -646,6 +646,62 @@ def colella_df(eta1: float, eta2: float, lx: float, ly: float, alpha: float, lz:
     df_out[2, 2] = lz
 
 
+@pure
+def magnetotail_slab(
+    eta1: float,
+    eta2: float,
+    eta3: float,
+    lx: float,
+    ly: float,
+    lz: float,
+    pinch: float,
+    pinch_width: float,
+    x_center: float,
+    f_out: "float[:]",
+):
+    """Point-wise evaluation of an elongated slab with a smooth central pinch."""
+
+    x = lx * (eta1 - 0.5)
+    x_shift = x - x_center
+    profile = 1.0 - pinch * exp(-(x_shift / pinch_width) ** 2)
+
+    f_out[0] = x
+    f_out[1] = ly * profile * (eta2 - 0.5)
+    f_out[2] = lz * (eta3 - 0.5)
+
+
+@pure
+def magnetotail_slab_df(
+    eta1: float,
+    eta2: float,
+    eta3: float,
+    lx: float,
+    ly: float,
+    lz: float,
+    pinch: float,
+    pinch_width: float,
+    x_center: float,
+    df_out: "float[:,:]",
+):
+    """Jacobian matrix for :meth:`struphy.geometry.mappings_kernels.magnetotail_slab`."""
+
+    x = lx * (eta1 - 0.5)
+    x_shift = x - x_center
+    gauss = exp(-(x_shift / pinch_width) ** 2)
+    profile = 1.0 - pinch * gauss
+    dprofile_dx = 2.0 * pinch * x_shift * gauss / (pinch_width * pinch_width)
+
+    df_out[0, 0] = lx
+    df_out[0, 1] = 0.0
+    df_out[0, 2] = 0.0
+    df_out[1, 0] = ly * (eta2 - 0.5) * dprofile_dx * lx
+    df_out[1, 1] = ly * profile
+    df_out[1, 2] = 0.0
+    df_out[2, 0] = 0.0
+    df_out[2, 1] = 0.0
+    df_out[2, 2] = lz
+
+
 @pure
 def hollow_cyl(eta1: float, eta2: float, eta3: float, a1: float, a2: float, lz: float, poc: float, f_out: "float[:]"):
     r"""Point-wise evaluation of
@@ -927,6 +983,404 @@ def hollow_torus_df(
         df_out[2, 2] = 0.0
 
 
+@pure
+def geospace_flux_domain(
+    eta1: float,
+    eta2: float,
+    eta3: float,
+    r_iono: float,
+    r_mp_dayside: float,
+    r_mp_tail: float,
+    r_bs_dayside: float,
+    r_bs_tail: float,
+    mp_eta1: float,
+    sheet_flatten: float,
+    f_out: "float[:]",
+):
+    """Point-wise evaluation of the GeospaceFluxDomain mapping."""
+
+    theta = pi * eta2
+    phi = 2 * pi * eta3
+
+    tail = 0.5 * (1.0 - cos(theta))
+
+    r_mp = r_mp_dayside + (r_mp_tail - r_mp_dayside) * tail
+    r_bs = r_bs_dayside + (r_bs_tail - r_bs_dayside) * tail
+
+    if eta1 <= mp_eta1:
+        u = eta1 / mp_eta1
+        r = r_iono + (r_mp - r_iono) * u
+    else:
+        u = (eta1 - mp_eta1) / (1.0 - mp_eta1)
+        r = r_mp + (r_bs - r_mp) * u
+
+    scale_z = 1.0 - sheet_flatten * tail
+
+    x = r * cos(theta)
+    rho = r * sin(theta)
+
+    f_out[0] = x
+    f_out[1] = rho * cos(phi)
+    f_out[2] = rho * sin(phi) * scale_z
+
+
+@pure
+def geospace_flux_domain_df(
+    eta1: float,
+    eta2: float,
+    eta3: float,
+    r_iono: float,
+    r_mp_dayside: float,
+    r_mp_tail: float,
+    r_bs_dayside: float,
+    r_bs_tail: float,
+    mp_eta1: float,
+    sheet_flatten: float,
+    df_out: "float[:,:]",
+):
+    """Jacobian matrix for :meth:`struphy.geometry.mappings_kernels.geospace_flux_domain`."""
+
+    theta = pi * eta2
+    phi = 2 * pi * eta3
+
+    sin_theta = sin(theta)
+    cos_theta = cos(theta)
+    sin_phi = sin(phi)
+    cos_phi = cos(phi)
+
+    tail = 0.5 * (1.0 - cos_theta)
+    dtail_deta2 = 0.5 * pi * sin_theta
+
+    r_mp = r_mp_dayside + (r_mp_tail - r_mp_dayside) * tail
+    r_bs = r_bs_dayside + (r_bs_tail - r_bs_dayside) * tail
+    drmp_deta2 = (r_mp_tail - r_mp_dayside) * dtail_deta2
+    drbs_deta2 = (r_bs_tail - r_bs_dayside) * dtail_deta2
+
+    if eta1 <= mp_eta1:
+        u = eta1 / mp_eta1
+        r = r_iono + (r_mp - r_iono) * u
+        dr_deta1 = (r_mp - r_iono) / mp_eta1
+        dr_deta2 = drmp_deta2 * u
+    else:
+        u = (eta1 - mp_eta1) / (1.0 - mp_eta1)
+        r = r_mp + (r_bs - r_mp) * u
+        dr_deta1 = (r_bs - r_mp) / (1.0 - mp_eta1)
+        dr_deta2 = drmp_deta2 * (1.0 - u) + drbs_deta2 * u
+
+    scale_z = 1.0 - sheet_flatten * tail
+    dscale_deta2 = -sheet_flatten * dtail_deta2
+
+    # x = r*cos(theta)
+    df_out[0, 0] = dr_deta1 * cos_theta
+    df_out[0, 1] = dr_deta2 * cos_theta - r * sin_theta * pi
+    df_out[0, 2] = 0.0
+
+    # rho = r*sin(theta)
+    drho_deta1 = dr_deta1 * sin_theta
+    drho_deta2 = dr_deta2 * sin_theta + r * cos_theta * pi
+
+    # y = rho*cos(phi)
+    df_out[1, 0] = drho_deta1 * cos_phi
+    df_out[1, 1] = drho_deta2 * cos_phi
+    df_out[1, 2] = -rho * sin_phi * 2 * pi
+
+    # z = rho*sin(phi)*scale_z
+    df_out[2, 0] = drho_deta1 * sin_phi * scale_z
+    df_out[2, 1] = drho_deta2 * sin_phi * scale_z + rho * sin_phi * dscale_deta2
+    df_out[2, 2] = rho * cos_phi * 2 * pi * scale_z
+
+
+@pure
+def warped_accretion_disk(
+    eta1: float,
+    eta2: float,
+    eta3: float,
+    r_in: float,
+    r_out: float,
+    thickness: float,
+    warp_amp: float,
+    warp_power: float,
+    node_angle: float,
+    tor_period: float,
+    f_out: "float[:]",
+):
+    """Point-wise evaluation of a warped accretion disk."""
+
+    dr = r_out - r_in
+    r = r_in + dr * eta1
+    s = (r - r_in) / dr
+    phi = 2 * pi * eta2 / tor_period
+    psi = phi - node_angle
+
+    z_local = thickness * (eta3 - 0.5)
+    w = warp_amp * r * (s**warp_power)
+
+    f_out[0] = r * cos(phi)
+    f_out[1] = r * sin(phi)
+    f_out[2] = z_local + w * sin(psi)
+
+
+@pure
+def warped_accretion_disk_df(
+    eta1: float,
+    eta2: float,
+    eta3: float,
+    r_in: float,
+    r_out: float,
+    thickness: float,
+    warp_amp: float,
+    warp_power: float,
+    node_angle: float,
+    tor_period: float,
+    df_out: "float[:,:]",
+):
+    """Jacobian matrix for :meth:`struphy.geometry.mappings_kernels.warped_accretion_disk`."""
+
+    dr = r_out - r_in
+    r = r_in + dr * eta1
+    s = (r - r_in) / dr
+    phi = 2 * pi * eta2 / tor_period
+    psi = phi - node_angle
+
+    dphi_deta2 = 2 * pi / tor_period
+
+    w = warp_amp * r * (s**warp_power)
+
+    if s == 0.0 and warp_power < 1.0:
+        dw_dr = 0.0
+    else:
+        dw_dr = warp_amp * (s**warp_power + r * warp_power * (s ** (warp_power - 1.0)) / dr)
+
+    dw_deta1 = dw_dr * dr
+
+    df_out[0, 0] = dr * cos(phi)
+    df_out[0, 1] = -r * sin(phi) * dphi_deta2
+    df_out[0, 2] = 0.0
+
+    df_out[1, 0] = dr * sin(phi)
+    df_out[1, 1] = r * cos(phi) * dphi_deta2
+    df_out[1, 2] = 0.0
+
+    df_out[2, 0] = dw_deta1 * sin(psi)
+    df_out[2, 1] = w * cos(psi) * dphi_deta2
+    df_out[2, 2] = thickness
+
+
+@pure
+def spheromak(
+    eta1: float,
+    eta2: float,
+    eta3: float,
+    r0: float,
+    a: float,
+    kappa: float,
+    tor_period: float,
+    f_out: "float[:]",
+):
+    """Point-wise evaluation of a compact spherical-shell spheromak proxy."""
+
+    r = r0 + a * eta1
+    theta = pi * eta2
+    phi = 2 * pi * eta3 / tor_period
+
+    sin_theta = sin(theta)
+    cos_theta = cos(theta)
+
+    f_out[0] = r * sin_theta * cos(phi)
+    f_out[1] = r * sin_theta * sin(phi)
+    f_out[2] = kappa * r * cos_theta
+
+
+@pure
+def spheromak_df(
+    eta1: float,
+    eta2: float,
+    eta3: float,
+    r0: float,
+    a: float,
+    kappa: float,
+    tor_period: float,
+    df_out: "float[:,:]",
+):
+    """Jacobian matrix for :meth:`struphy.geometry.mappings_kernels.spheromak`."""
+
+    r = r0 + a * eta1
+    theta = pi * eta2
+    phi = 2 * pi * eta3 / tor_period
+
+    sin_theta = sin(theta)
+    cos_theta = cos(theta)
+    sin_phi = sin(phi)
+    cos_phi = cos(phi)
+
+    dtheta_deta2 = pi
+    dphi_deta3 = 2 * pi / tor_period
+
+    df_out[0, 0] = a * sin_theta * cos_phi
+    df_out[0, 1] = r * cos_theta * dtheta_deta2 * cos_phi
+    df_out[0, 2] = -r * sin_theta * sin_phi * dphi_deta3
+
+    df_out[1, 0] = a * sin_theta * sin_phi
+    df_out[1, 1] = r * cos_theta * dtheta_deta2 * sin_phi
+    df_out[1, 2] = r * sin_theta * cos_phi * dphi_deta3
+
+    df_out[2, 0] = kappa * a * cos_theta
+    df_out[2, 1] = -kappa * r * sin_theta * dtheta_deta2
+    df_out[2, 2] = 0.0
+
+
+@pure
+def hall_effect_thruster_channel(
+    eta1: float,
+    eta2: float,
+    eta3: float,
+    r_in: float,
+    r_out: float,
+    length: float,
+    pack_strength: float,
+    tor_period: float,
+    f_out: "float[:]",
+):
+    """Point-wise evaluation of a coaxial Hall thruster channel."""
+
+    dr = r_out - r_in
+    r = r_in + dr * eta1
+    phi = 2 * pi * eta2 / tor_period
+
+    # Axial clustering near eta3=0 and eta3=1 (anode and exit)
+    z = length * (eta3 - pack_strength * sin(2 * pi * eta3) / (2 * pi))
+
+    f_out[0] = r * cos(phi)
+    f_out[1] = r * sin(phi)
+    f_out[2] = z
+
+
+@pure
+def hall_effect_thruster_channel_df(
+    eta1: float,
+    eta2: float,
+    eta3: float,
+    r_in: float,
+    r_out: float,
+    length: float,
+    pack_strength: float,
+    tor_period: float,
+    df_out: "float[:,:]",
+):
+    """Jacobian matrix for :meth:`struphy.geometry.mappings_kernels.hall_effect_thruster_channel`."""
+
+    dr = r_out - r_in
+    r = r_in + dr * eta1
+    phi = 2 * pi * eta2 / tor_period
+
+    dphi_deta2 = 2 * pi / tor_period
+    dz_deta3 = length * (1.0 - pack_strength * cos(2 * pi * eta3))
+
+    df_out[0, 0] = dr * cos(phi)
+    df_out[0, 1] = -r * sin(phi) * dphi_deta2
+    df_out[0, 2] = 0.0
+
+    df_out[1, 0] = dr * sin(phi)
+    df_out[1, 1] = r * cos(phi) * dphi_deta2
+    df_out[1, 2] = 0.0
+
+    df_out[2, 0] = 0.0
+    df_out[2, 1] = 0.0
+    df_out[2, 2] = dz_deta3
+
+
+@pure
+def diagnostic_port_hole_torus(
+    eta1: float,
+    eta2: float,
+    eta3: float,
+    a1: float,
+    a2: float,
+    r0: float,
+    tor_period: float,
+    port_depth: float,
+    port_eta2_center: float,
+    port_eta3_center: float,
+    port_eta2_width: float,
+    port_eta3_width: float,
+    f_out: "float[:]",
+):
+    """Point-wise evaluation of a torus with a localized diagnostic port deformation."""
+
+    da = a2 - a1
+    theta = 2 * pi * eta2
+    phi = 2 * pi * eta3 / tor_period
+
+    sigma2_sq = port_eta2_width * port_eta2_width
+    sigma3_sq = port_eta3_width * port_eta3_width
+
+    bump_theta = exp((cos(2 * pi * (eta2 - port_eta2_center)) - 1.0) / sigma2_sq)
+    bump_phi = exp((cos(2 * pi * (eta3 - port_eta3_center)) - 1.0) / sigma3_sq)
+    bump = bump_theta * bump_phi
+
+    r_minor = a1 + da * eta1 + port_depth * eta1 * eta1 * bump
+    r_major = r0 + r_minor * cos(theta)
+
+    f_out[0] = r_major * cos(phi)
+    f_out[1] = -r_major * sin(phi)
+    f_out[2] = r_minor * sin(theta)
+
+
+@pure
+def diagnostic_port_hole_torus_df(
+    eta1: float,
+    eta2: float,
+    eta3: float,
+    a1: float,
+    a2: float,
+    r0: float,
+    tor_period: float,
+    port_depth: float,
+    port_eta2_center: float,
+    port_eta3_center: float,
+    port_eta2_width: float,
+    port_eta3_width: float,
+    df_out: "float[:,:]",
+):
+    """Jacobian matrix for :meth:`struphy.geometry.mappings_kernels.diagnostic_port_hole_torus`."""
+
+    da = a2 - a1
+    theta = 2 * pi * eta2
+    phi = 2 * pi * eta3 / tor_period
+
+    sigma2_sq = port_eta2_width * port_eta2_width
+    sigma3_sq = port_eta3_width * port_eta3_width
+
+    bump_theta = exp((cos(2 * pi * (eta2 - port_eta2_center)) - 1.0) / sigma2_sq)
+    bump_phi = exp((cos(2 * pi * (eta3 - port_eta3_center)) - 1.0) / sigma3_sq)
+    bump = bump_theta * bump_phi
+
+    dbump_deta2 = bump * (-2 * pi * sin(2 * pi * (eta2 - port_eta2_center)) / sigma2_sq)
+    dbump_deta3 = bump * (-2 * pi * sin(2 * pi * (eta3 - port_eta3_center)) / sigma3_sq)
+
+    r_minor = a1 + da * eta1 + port_depth * eta1 * eta1 * bump
+    dr_deta1 = da + 2 * port_depth * eta1 * bump
+    dr_deta2 = port_depth * eta1 * eta1 * dbump_deta2
+    dr_deta3 = port_depth * eta1 * eta1 * dbump_deta3
+
+    r_major = r0 + r_minor * cos(theta)
+    d_rmajor_deta1 = dr_deta1 * cos(theta)
+    d_rmajor_deta2 = dr_deta2 * cos(theta) - 2 * pi * r_minor * sin(theta)
+    d_rmajor_deta3 = dr_deta3 * cos(theta)
+
+    df_out[0, 0] = d_rmajor_deta1 * cos(phi)
+    df_out[0, 1] = d_rmajor_deta2 * cos(phi)
+    df_out[0, 2] = d_rmajor_deta3 * cos(phi) - 2 * pi / tor_period * r_major * sin(phi)
+
+    df_out[1, 0] = -d_rmajor_deta1 * sin(phi)
+    df_out[1, 1] = -d_rmajor_deta2 * sin(phi)
+    df_out[1, 2] = -d_rmajor_deta3 * sin(phi) - 2 * pi / tor_period * r_major * cos(phi)
+
+    df_out[2, 0] = dr_deta1 * sin(theta)
+    df_out[2, 1] = dr_deta2 * sin(theta) + 2 * pi * r_minor * cos(theta)
+    df_out[2, 2] = dr_deta3 * sin(theta)
+
+
 @pure
 def shafranov_shift(
     eta1: float,
@@ -1160,3 +1614,70 @@ def shafranov_dshaped_df(
     df_out[2, 0] = 0.0
     df_out[2, 1] = 0.0
     df_out[2, 2] = lz
+
+
+@pure
+def moebius_strip(
+    eta1: float,
+    eta2: float,
+    eta3: float,
+    r0: float,
+    width: float,
+    thickness: float,
+    f_out: "float[:]",
+):
+    """Point-wise evaluation of the thickened Moebius strip mapping."""
+
+    w_hat = width * (eta1 - 0.5)
+    t_hat = thickness * (eta2 - 0.5)
+
+    alpha = pi * eta3
+    phi = 2 * pi * eta3
+
+    r_eff = r0 + w_hat * cos(alpha) - t_hat * sin(alpha)
+
+    f_out[0] = r_eff * cos(phi)
+    f_out[1] = r_eff * sin(phi)
+    f_out[2] = w_hat * sin(alpha) + t_hat * cos(alpha)
+
+
+@pure
+def moebius_strip_df(
+    eta1: float,
+    eta2: float,
+    eta3: float,
+    r0: float,
+    width: float,
+    thickness: float,
+    df_out: "float[:,:]",
+):
+    """Jacobian matrix for :meth:`struphy.geometry.mappings_kernels.moebius_strip`."""
+
+    w_hat = width * (eta1 - 0.5)
+    t_hat = thickness * (eta2 - 0.5)
+
+    alpha = pi * eta3
+    phi = 2 * pi * eta3
+
+    sin_alpha = sin(alpha)
+    cos_alpha = cos(alpha)
+    sin_phi = sin(phi)
+    cos_phi = cos(phi)
+
+    r_eff = r0 + w_hat * cos_alpha - t_hat * sin_alpha
+
+    dr_deta1 = width * cos_alpha
+    dr_deta2 = -thickness * sin_alpha
+    dr_deta3 = -(w_hat * sin_alpha + t_hat * cos_alpha) * pi
+
+    df_out[0, 0] = dr_deta1 * cos_phi
+    df_out[0, 1] = dr_deta2 * cos_phi
+    df_out[0, 2] = dr_deta3 * cos_phi - 2 * pi * r_eff * sin_phi
+
+    df_out[1, 0] = dr_deta1 * sin_phi
+    df_out[1, 1] = dr_deta2 * sin_phi
+    df_out[1, 2] = dr_deta3 * sin_phi + 2 * pi * r_eff * cos_phi
+
+    df_out[2, 0] = width * sin_alpha
+    df_out[2, 1] = thickness * cos_alpha
+    df_out[2, 2] = (w_hat * cos_alpha - t_hat * sin_alpha) * pi
diff --git a/src/struphy/geometry/tests/test_domain.py b/src/struphy/geometry/tests/test_domain.py
index 6ee8a5ee3..35acf100c 100644
--- a/src/struphy/geometry/tests/test_domain.py
+++ b/src/struphy/geometry/tests/test_domain.py
@@ -141,10 +141,16 @@ def a_vec(e1, e2, e3):
     "mapping",
     [
         "Cuboid",
+        "MagnetotailSlab",
         "HollowCylinder",
         "Colella",
         "Orthogonal",
         "HollowTorus",
+        "DiagnosticPortHoleTorus",
+        "GeospaceFluxDomain",
+        "WarpedAccretionDisk",
+        "Spheromak",
+        "HallEffectThrusterChannel",
         "PoweredEllipticCylinder",
         "ShafranovShiftCylinder",
         "ShafranovSqrtCylinder",