Skip to content

models.grain.geometries

Concrete grain segment geometries, both analytical and FMM-based.

2D (constant cross-section):

  • BatesSegment — Standard cylindrical grain with axial core (the most common hobby/amateur geometry).
  • TubularSegment — Thin-walled tube, burns from inner and/or outer surfaces.
  • RodAndTubeSegment — Central rod surrounded by an outer tubular section.
  • MultiPortGrainSegment — Multiple circular ports arranged radially in the cross-section.
  • StarGrainSegment — Star-shaped port (FMM-based). Configurable point count, length, and width.
  • DGrainSegment — D-shaped slot port (FMM-based).

3D (varying cross-section):

  • ConicalGrainSegment — Linearly tapered core from upper to lower diameter.
  • WagonWheelSegment — Complex multi-spoke geometry (FMM-based).

All segments implement the GrainSegment interface and can be mixed within a single Grain.

machwave.models.grain.geometries

BatesSegment

Bases: GrainSegment2D

Source code in machwave/models/grain/geometries/bates.py 
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
class BatesSegment(GrainSegment2D):
    INHIBITED_SURFACES = InhibitedSurfaces(
        outer_surface=True,
        inner_surface=False,
        upper_end=False,
        lower_end=False,
    )

    def __init__(
        self,
        outer_diameter: float,
        core_diameter: float,
        length: float,
        density_ratio: float = 1.0,
    ) -> None:
        self.core_diameter = core_diameter

        super().__init__(
            length=length,
            outer_diameter=outer_diameter,
            inhibited_surfaces=self.INHIBITED_SURFACES,
            density_ratio=density_ratio,
        )

    def validate(self) -> None:
        super().validate()

        if not self.outer_diameter > self.core_diameter:
            raise GrainGeometryError(
                f"Outer diameter ({self.outer_diameter}) must be greater than "
                f"core diameter ({self.core_diameter})"
            )
        if not self.core_diameter > 0:
            raise GrainGeometryError(
                f"Core diameter must be positive, got {self.core_diameter}"
            )

    def get_core_diameter(self, web_distance: float) -> float:
        return self.core_diameter + 2 * web_distance

    def get_port_area(self, web_distance: float) -> float:
        return get_circle_area(diameter=self.get_core_diameter(web_distance))

    def get_core_area(self, web_distance: float) -> float:
        length = self.get_length(web_distance=web_distance)
        core_diameter = self.core_diameter + 2 * web_distance
        return get_cylinder_surface_area(length, core_diameter)

    def get_face_area(self, web_distance: float) -> float:
        core_diameter = self.get_core_diameter(web_distance)
        return np.pi * (((self.outer_diameter**2) - (core_diameter) ** 2) / 4)

    def get_web_thickness(self) -> float:
        """
        More details on the web thickness of BATES grains can be found in:
        https://www.nakka-rocketry.net/design1.html
        """
        return 0.5 * (self.outer_diameter - self.core_diameter)

    def get_optimal_length(self) -> float:
        """
        Returns the optimal length for BATES segment.
        More details on the calculation:
        https://www.nakka-rocketry.net/th_grain.html

        :return: Optimal length for neutral burn of BATES segment
        :rtype: float
        """
        return 1e3 * 0.5 * (3 * self.outer_diameter + self.core_diameter)

    def get_center_of_gravity(
        self, web_distance: float = 0.0
    ) -> np.typing.NDArray[np.float64]:
        """
        Calculate the center of gravity of the BATES grain segment.

        BATES is a symmetrical 2D geometry that burns radially. Due to its
        cylindrical symmetry, the center of gravity remains constant at the
        geometric center regardless of web distance burned.

        Args:
            web_distance: Web distance traveled (unused for BATES due to symmetry),
                         in meters. Included for API consistency.

        Returns:
            Center of gravity in 3D space [x, y, z], in meters, measured from
            the aft end (port, closest to nozzle). Always returns [length/2, 0, 0]
            for symmetric BATES grains.
        """
        # For symmetric BATES grains, CoG doesn't change with burn
        # The grain burns radially inward, maintaining axial symmetry
        # CoG is at geometric center, which is length/2 from the aft end (port)
        return np.array([self.length / 2, 0.0, 0.0], dtype=np.float64)

    def get_moment_of_inertia(
        self, ideal_density: float, web_distance: float = 0.0
    ) -> np.typing.NDArray[np.float64]:
        """
        Calculate the moment of inertia tensor of the BATES grain segment (hollow
        cylinder) at its center of gravity.

        Args:
            ideal_density: Propellant ideal density [kg/m^3].
            web_distance: Web distance traveled [m].

        Returns:
            A 3x3 inertia tensor [kg-m^2].
        """
        # Geometry at given web distance
        r_inner = (self.core_diameter + 2 * web_distance) / 2
        r_outer = self.outer_diameter / 2
        current_length = self.get_length(web_distance)

        volume = self.get_volume(web_distance)
        mass = volume * ideal_density * self.density_ratio

        r_sum_sq = r_inner**2 + r_outer**2

        # Ixx: moment about axial axis
        Ixx = mass * r_sum_sq / 2

        # Iyy, Izz: moments about radial axes
        Iyy = mass * (r_sum_sq / 4 + current_length**2 / 12)
        Izz = Iyy

        return np.array(
            [[Ixx, 0.0, 0.0], [0.0, Iyy, 0.0], [0.0, 0.0, Izz]], dtype=np.float64
        )

get_center_of_gravity(web_distance=0.0)

Calculate the center of gravity of the BATES grain segment.

BATES is a symmetrical 2D geometry that burns radially. Due to its cylindrical symmetry, the center of gravity remains constant at the geometric center regardless of web distance burned.

Parameters:

Name Type Description Default
web_distance float

Web distance traveled (unused for BATES due to symmetry), in meters. Included for API consistency.

 0.0

Returns:

Type Description
NDArray[float64]

Center of gravity in 3D space [x, y, z], in meters, measured from
NDArray[float64]

the aft end (port, closest to nozzle). Always returns [length/2, 0, 0]
NDArray[float64]

for symmetric BATES grains.
Source code in machwave/models/grain/geometries/bates.py 
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
def get_center_of_gravity(
    self, web_distance: float = 0.0
) -> np.typing.NDArray[np.float64]:
    """
    Calculate the center of gravity of the BATES grain segment.

    BATES is a symmetrical 2D geometry that burns radially. Due to its
    cylindrical symmetry, the center of gravity remains constant at the
    geometric center regardless of web distance burned.

    Args:
        web_distance: Web distance traveled (unused for BATES due to symmetry),
                     in meters. Included for API consistency.

    Returns:
        Center of gravity in 3D space [x, y, z], in meters, measured from
        the aft end (port, closest to nozzle). Always returns [length/2, 0, 0]
        for symmetric BATES grains.
    """
    # For symmetric BATES grains, CoG doesn't change with burn
    # The grain burns radially inward, maintaining axial symmetry
    # CoG is at geometric center, which is length/2 from the aft end (port)
    return np.array([self.length / 2, 0.0, 0.0], dtype=np.float64)

get_moment_of_inertia(ideal_density, web_distance=0.0)

Calculate the moment of inertia tensor of the BATES grain segment (hollow cylinder) at its center of gravity.

Parameters:

Name Type Description Default
ideal_density float

Propellant ideal density [kg/m^3].

 required
web_distance float

Web distance traveled [m].

 0.0

Returns:

Type Description
NDArray[float64]

A 3x3 inertia tensor [kg-m^2].
Source code in machwave/models/grain/geometries/bates.py 
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
def get_moment_of_inertia(
    self, ideal_density: float, web_distance: float = 0.0
) -> np.typing.NDArray[np.float64]:
    """
    Calculate the moment of inertia tensor of the BATES grain segment (hollow
    cylinder) at its center of gravity.

    Args:
        ideal_density: Propellant ideal density [kg/m^3].
        web_distance: Web distance traveled [m].

    Returns:
        A 3x3 inertia tensor [kg-m^2].
    """
    # Geometry at given web distance
    r_inner = (self.core_diameter + 2 * web_distance) / 2
    r_outer = self.outer_diameter / 2
    current_length = self.get_length(web_distance)

    volume = self.get_volume(web_distance)
    mass = volume * ideal_density * self.density_ratio

    r_sum_sq = r_inner**2 + r_outer**2

    # Ixx: moment about axial axis
    Ixx = mass * r_sum_sq / 2

    # Iyy, Izz: moments about radial axes
    Iyy = mass * (r_sum_sq / 4 + current_length**2 / 12)
    Izz = Iyy

    return np.array(
        [[Ixx, 0.0, 0.0], [0.0, Iyy, 0.0], [0.0, 0.0, Izz]], dtype=np.float64
    )

get_optimal_length()

Returns the optimal length for BATES segment. More details on the calculation: https://www.nakka-rocketry.net/th_grain.html

:return: Optimal length for neutral burn of BATES segment :rtype: float

Source code in machwave/models/grain/geometries/bates.py 
70
71
72
73
74
75
76
77
78
79
def get_optimal_length(self) -> float:
    """
    Returns the optimal length for BATES segment.
    More details on the calculation:
    https://www.nakka-rocketry.net/th_grain.html

    :return: Optimal length for neutral burn of BATES segment
    :rtype: float
    """
    return 1e3 * 0.5 * (3 * self.outer_diameter + self.core_diameter)

get_web_thickness()

More details on the web thickness of BATES grains can be found in: https://www.nakka-rocketry.net/design1.html

Source code in machwave/models/grain/geometries/bates.py 
63
64
65
66
67
68
def get_web_thickness(self) -> float:
    """
    More details on the web thickness of BATES grains can be found in:
    https://www.nakka-rocketry.net/design1.html
    """
    return 0.5 * (self.outer_diameter - self.core_diameter)

MultiPortGrainSegment

Bases: FMMGrainSegment2D

Source code in machwave/models/grain/geometries/multi_port.py 
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
class MultiPortGrainSegment(FMMGrainSegment2D):
    def __init__(
        self,
        length: float,
        outer_diameter: float,
        port_diameter: float,
        port_radial_count: float,
        port_level_count: float,
        inhibited_surfaces: InhibitedSurfaces | None = None,
        density_ratio: float = 1.0,
    ) -> None:
        self.port_diameter = port_diameter
        self.port_radial_count = int(port_radial_count)
        self.port_level_count = int(port_level_count)

        super().__init__(
            length=length,
            outer_diameter=outer_diameter,
            inhibited_surfaces=inhibited_surfaces,
            density_ratio=density_ratio,
        )

    def validate(self) -> None:
        super().validate()

        if not self.port_diameter > 0:
            raise GrainGeometryError(
                f"Port diameter must be positive, got {self.port_diameter}"
            )
        if not self.port_level_count > 0:
            raise GrainGeometryError(
                f"Port level count must be positive, got {self.port_level_count}"
            )
        max_port_size = self.outer_diameter / 2
        total_port_size = self.port_level_count * self.port_diameter
        if not total_port_size < max_port_size:
            raise GrainGeometryError(
                f"Total port size ({total_port_size}) must be less than "
                f"half the outer diameter ({max_port_size})"
            )
        if not self.port_radial_count > 0:
            raise GrainGeometryError(
                f"Port radial count must be positive, got {self.port_radial_count}"
            )

    def get_initial_face_map(self) -> np.typing.NDArray[np.int_]:
        """
        NOTE: Still needs to correctly implement wagon wheel ports.
        """
        map_x, map_y = self.get_maps()
        core_map = self.get_empty_face_map()

        od_norm = self.normalize(self.outer_diameter)
        port_od_norm = self.normalize(self.port_diameter)

        for radius in range(self.port_radial_count):
            angle = np.pi * 2 * radius / self.port_radial_count

            for level in range(self.port_level_count):
                radial_distance = od_norm * level / (self.port_level_count) / 2

                x_offset = radial_distance * np.cos(angle)
                y_offset = radial_distance * np.sin(angle)

                radius = np.sqrt((map_x - x_offset) ** 2 + (map_y - y_offset) ** 2)
                core_map[radius < port_od_norm / 2] = 0

        return core_map

get_initial_face_map()

NOTE: Still needs to correctly implement wagon wheel ports.

Source code in machwave/models/grain/geometries/multi_port.py 
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
def get_initial_face_map(self) -> np.typing.NDArray[np.int_]:
    """
    NOTE: Still needs to correctly implement wagon wheel ports.
    """
    map_x, map_y = self.get_maps()
    core_map = self.get_empty_face_map()

    od_norm = self.normalize(self.outer_diameter)
    port_od_norm = self.normalize(self.port_diameter)

    for radius in range(self.port_radial_count):
        angle = np.pi * 2 * radius / self.port_radial_count

        for level in range(self.port_level_count):
            radial_distance = od_norm * level / (self.port_level_count) / 2

            x_offset = radial_distance * np.cos(angle)
            y_offset = radial_distance * np.sin(angle)

            radius = np.sqrt((map_x - x_offset) ** 2 + (map_y - y_offset) ** 2)
            core_map[radius < port_od_norm / 2] = 0

    return core_map

RodAndTubeGrainSegment

Bases: FMMGrainSegment2D

Source code in machwave/models/grain/geometries/rod_and_tube.py 
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
class RodAndTubeGrainSegment(FMMGrainSegment2D):
    def __init__(
        self,
        length: float,
        outer_diameter: float,
        rod_outer_diameter: float,
        tube_inner_diameter: float,
        inhibited_surfaces: InhibitedSurfaces | None = None,
        density_ratio: float = 1.0,
    ) -> None:
        self.rod_outer_diameter = rod_outer_diameter
        self.tube_inner_diameter = tube_inner_diameter

        super().__init__(
            length=length,
            outer_diameter=outer_diameter,
            inhibited_surfaces=inhibited_surfaces,
            density_ratio=density_ratio,
        )

    def validate(self) -> None:
        super().validate()

        if not self.rod_outer_diameter > 0:
            raise GrainGeometryError(
                f"Rod outer diameter must be positive, got {self.rod_outer_diameter}"
            )
        if not self.tube_inner_diameter > self.rod_outer_diameter:
            raise GrainGeometryError(
                f"Tube inner diameter ({self.tube_inner_diameter}) must be greater than "
                f"rod outer diameter ({self.rod_outer_diameter})"
            )
        if not self.tube_inner_diameter < self.outer_diameter:
            raise GrainGeometryError(
                f"Tube inner diameter ({self.tube_inner_diameter}) must be less than "
                f"outer diameter ({self.outer_diameter})"
            )

    def get_initial_face_map(self) -> np.typing.NDArray[np.int_]:
        """
        NOTE: Still needs to correctly implement wagon wheel ports.
        """
        map_x, map_y = self.get_maps()
        core_map = self.get_empty_face_map()

        rod_od_norm = self.normalize(self.rod_outer_diameter)
        tube_id_norm = self.normalize(self.tube_inner_diameter)

        radius = np.sqrt(map_x**2 + map_y**2)

        # Create the ring:
        core_map[(radius > rod_od_norm / 2) & (radius < tube_id_norm / 2)] = 0

        return core_map

get_initial_face_map()

NOTE: Still needs to correctly implement wagon wheel ports.

Source code in machwave/models/grain/geometries/rod_and_tube.py 
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
def get_initial_face_map(self) -> np.typing.NDArray[np.int_]:
    """
    NOTE: Still needs to correctly implement wagon wheel ports.
    """
    map_x, map_y = self.get_maps()
    core_map = self.get_empty_face_map()

    rod_od_norm = self.normalize(self.rod_outer_diameter)
    tube_id_norm = self.normalize(self.tube_inner_diameter)

    radius = np.sqrt(map_x**2 + map_y**2)

    # Create the ring:
    core_map[(radius > rod_od_norm / 2) & (radius < tube_id_norm / 2)] = 0

    return core_map

StarGrainSegment

Bases: FMMGrainSegment2D

Source code in machwave/models/grain/geometries/star.py 
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
class StarGrainSegment(FMMGrainSegment2D):
    def __init__(
        self,
        length: float,
        outer_diameter: float,
        number_of_points: int,
        point_length: float,
        point_width: float,
        inhibited_surfaces: InhibitedSurfaces | None = None,
        density_ratio: float = 1.0,
    ) -> None:
        self.number_of_points = int(number_of_points)
        self.point_length = point_length
        self.point_width = point_width

        super().__init__(
            length=length,
            outer_diameter=outer_diameter,
            inhibited_surfaces=inhibited_surfaces,
            density_ratio=density_ratio,
        )

    def validate(self) -> None:
        super().validate()

        if not self.number_of_points > 0:
            raise GrainGeometryError(
                f"Number of points must be positive, got {self.number_of_points}"
            )
        if not self.number_of_points < 12:
            raise GrainGeometryError(
                f"Number of points must be less than 12, got {self.number_of_points}"
            )
        if not isinstance(self.number_of_points, int):
            raise GrainGeometryError(
                f"Number of points must be an integer, got {type(self.number_of_points).__name__}"
            )
        if not self.point_length > 0:
            raise GrainGeometryError(
                f"Point length must be positive, got {self.point_length}"
            )
        if not self.point_width > 0:
            raise GrainGeometryError(
                f"Point width must be positive, got {self.point_width}"
            )

    def get_initial_face_map(self) -> np.typing.NDArray[np.int_]:
        """
        This method returns the initial face map for a star grain segment.

        References:
        openMotor, https://github.com/reilleya/openMotor
        """
        map_x, map_y = self.get_maps()
        core_map = self.get_empty_face_map()

        point_length_norm = self.normalize(self.point_length)
        point_width_norm = self.normalize(self.point_width)

        radius = (map_x**2 + map_y**2) ** 0.5

        for i in range(0, self.number_of_points):
            theta = 2 * np.pi / self.number_of_points * i
            rect = abs(np.cos(theta) * map_x + np.sin(theta) * map_y)

            width = point_width_norm / 2 * (1 - (radius / point_length_norm))
            vect = rect < width
            near = np.sin(theta) * map_x - np.cos(theta) * map_y > -0.025

            core_map[np.logical_and(vect, near)] = 0

        return core_map

get_initial_face_map()

This method returns the initial face map for a star grain segment.

References: openMotor, https://github.com/reilleya/openMotor

Source code in machwave/models/grain/geometries/star.py 
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
def get_initial_face_map(self) -> np.typing.NDArray[np.int_]:
    """
    This method returns the initial face map for a star grain segment.

    References:
    openMotor, https://github.com/reilleya/openMotor
    """
    map_x, map_y = self.get_maps()
    core_map = self.get_empty_face_map()

    point_length_norm = self.normalize(self.point_length)
    point_width_norm = self.normalize(self.point_width)

    radius = (map_x**2 + map_y**2) ** 0.5

    for i in range(0, self.number_of_points):
        theta = 2 * np.pi / self.number_of_points * i
        rect = abs(np.cos(theta) * map_x + np.sin(theta) * map_y)

        width = point_width_norm / 2 * (1 - (radius / point_length_norm))
        vect = rect < width
        near = np.sin(theta) * map_x - np.cos(theta) * map_y > -0.025

        core_map[np.logical_and(vect, near)] = 0

    return core_map

WagonWheelGrainSegment

Bases: FMMGrainSegment2D

Source code in machwave/models/grain/geometries/wagon_wheel.py 
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
class WagonWheelGrainSegment(FMMGrainSegment2D):
    def __init__(
        self,
        length: float,
        outer_diameter: float,
        core_diameter: float,
        number_of_ports: int,
        port_inner_diameter: float,
        port_outer_diameter: float,
        port_angular_width: float,
        inhibited_surfaces: InhibitedSurfaces | None = None,
        density_ratio: float = 1.0,
    ) -> None:
        self.core_diameter = core_diameter
        self.number_of_ports = int(number_of_ports)
        self.port_inner_diameter = port_inner_diameter
        self.port_outer_diameter = port_outer_diameter
        self.port_angular_width = port_angular_width

        super().__init__(
            length=length,
            outer_diameter=outer_diameter,
            inhibited_surfaces=inhibited_surfaces,
            density_ratio=density_ratio,
        )

    def validate(self) -> None:
        super().validate()

        if not self.number_of_ports > 0:
            raise GrainGeometryError(
                f"Number of ports must be positive, got {self.number_of_ports}"
            )
        if not self.number_of_ports < 12:
            raise GrainGeometryError(
                f"Number of ports must be less than 12, got {self.number_of_ports}"
            )
        if not self.number_of_ports % 2 == 0:
            raise GrainGeometryError(
                f"Number of ports must be even, got {self.number_of_ports}"
            )
        if not isinstance(self.number_of_ports, int):
            raise GrainGeometryError(
                f"Number of ports must be an integer, got {type(self.number_of_ports).__name__}"
            )
        if not self.port_inner_diameter > self.core_diameter:
            raise GrainGeometryError(
                f"Port inner diameter ({self.port_inner_diameter}) must be greater than "
                f"core diameter ({self.core_diameter})"
            )
        if not self.port_outer_diameter > self.port_inner_diameter:
            raise GrainGeometryError(
                f"Port outer diameter ({self.port_outer_diameter}) must be greater than "
                f"port inner diameter ({self.port_inner_diameter})"
            )
        if not self.port_angular_width > 0:
            raise GrainGeometryError(
                f"Port angular width must be positive, got {self.port_angular_width}"
            )
        max_angular_width = 360 / self.number_of_ports
        if not self.port_angular_width < max_angular_width:
            raise GrainGeometryError(
                f"Port angular width ({self.port_angular_width}) must be less than "
                f"{max_angular_width} (360 / {self.number_of_ports})"
            )

    def get_initial_face_map(self) -> np.typing.NDArray[np.int_]:
        """
        NOTE: Still needs to correctly implement wagon wheel ports.
        """
        map_x, map_y = self.get_maps()
        core_map = self.get_empty_face_map()

        core_diameter_norm = self.normalize(self.core_diameter)
        port_inner_diameter_norm = self.normalize(self.port_inner_diameter)
        port_outer_diameter_norm = self.normalize(self.port_outer_diameter)

        radius = np.sqrt(map_x**2 + map_y**2)

        # Create the core:
        core_map[radius < core_diameter_norm / 2] = 0

        # Create the ports:
        for port_index in range(int(self.number_of_ports)):
            displacement_angle = 2 * np.pi / self.number_of_ports * (port_index)

            theta_2 = np.deg2rad(self.port_angular_width / 2) + displacement_angle
            theta_1 = displacement_angle - np.deg2rad(self.port_angular_width / 2)

            map_x_y_arctan = np.arctan(map_y / map_x)

            core_map[
                (radius < port_outer_diameter_norm / 2)
                & (radius > port_inner_diameter_norm / 2)
                & (np.abs(map_x_y_arctan) < theta_2)
                & (np.abs(map_x_y_arctan) > theta_1)
            ] = 0

        return core_map

get_initial_face_map()

NOTE: Still needs to correctly implement wagon wheel ports.

Source code in machwave/models/grain/geometries/wagon_wheel.py 
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
def get_initial_face_map(self) -> np.typing.NDArray[np.int_]:
    """
    NOTE: Still needs to correctly implement wagon wheel ports.
    """
    map_x, map_y = self.get_maps()
    core_map = self.get_empty_face_map()

    core_diameter_norm = self.normalize(self.core_diameter)
    port_inner_diameter_norm = self.normalize(self.port_inner_diameter)
    port_outer_diameter_norm = self.normalize(self.port_outer_diameter)

    radius = np.sqrt(map_x**2 + map_y**2)

    # Create the core:
    core_map[radius < core_diameter_norm / 2] = 0

    # Create the ports:
    for port_index in range(int(self.number_of_ports)):
        displacement_angle = 2 * np.pi / self.number_of_ports * (port_index)

        theta_2 = np.deg2rad(self.port_angular_width / 2) + displacement_angle
        theta_1 = displacement_angle - np.deg2rad(self.port_angular_width / 2)

        map_x_y_arctan = np.arctan(map_y / map_x)

        core_map[
            (radius < port_outer_diameter_norm / 2)
            & (radius > port_inner_diameter_norm / 2)
            & (np.abs(map_x_y_arctan) < theta_2)
            & (np.abs(map_x_y_arctan) > theta_1)
        ] = 0

    return core_map