From ec88a9c385a787f826ac8e0ccfaf4cac93740e11 Mon Sep 17 00:00:00 2001
From: Connor Moreno <camoreno@Denali.ep.wisc.edu>
Date: Tue, 27 Aug 2024 13:13:10 -0500
Subject: [PATCH 1/7] Change magnet modeling workflow to use CadQuery

---
 Examples/config.yaml          |   5 +-
 Examples/parastell_example.py |   5 +-
 parastell/magnet_coils.py     | 589 ++++++++++++++--------------------
 parastell/parastell.py        |  73 ++---
 tests/test_magnet_coils.py    | 101 +++---
 tests/test_parastell.py       |  15 +-
 6 files changed, 328 insertions(+), 460 deletions(-)

diff --git a/Examples/config.yaml b/Examples/config.yaml
index 263e7707..f7ccef8b 100644
--- a/Examples/config.yaml
+++ b/Examples/config.yaml
@@ -68,10 +68,11 @@ invessel_build:
 
 magnet_coils:
   coils_file: coils.example
-  cross_section: [circle, 20.0]
+  width: 40.0
+  thickness: 50.0
   toroidal_extent: 90.0
   sample_mod: 6
-  step_filename: magnet_coils
+  step_filename: magnet_set
   export_mesh: True
   mesh_filename: magnet_mesh
 
diff --git a/Examples/parastell_example.py b/Examples/parastell_example.py
index b5d2bcc5..65d678a9 100644
--- a/Examples/parastell_example.py
+++ b/Examples/parastell_example.py
@@ -54,11 +54,12 @@
 
 # Define build parameters for magnet coils
 coils_file = "coils.example"
-cross_section = ["circle", 20]
+width = 40.0
+thickness = 50.0
 toroidal_extent = 90.0
 # Construct magnets
 stellarator.construct_magnets(
-    coils_file, cross_section, toroidal_extent, sample_mod=6
+    coils_file, width, thickness, toroidal_extent, sample_mod=6
 )
 # Export magnet files
 stellarator.export_magnets(
diff --git a/parastell/magnet_coils.py b/parastell/magnet_coils.py
index 0bd5ff95..4ac72947 100644
--- a/parastell/magnet_coils.py
+++ b/parastell/magnet_coils.py
@@ -1,7 +1,9 @@
 import argparse
+from pathlib import Path
 
 import numpy as np
 
+import cadquery as cq
 import cubit
 
 from . import log
@@ -11,17 +13,46 @@
 export_allowed_kwargs = ["step_filename", "export_mesh", "mesh_filename"]
 
 
+def compute_tangent(prev_line, next_line):
+    """Computes tangent at "current" filament point using central difference
+    approximation and previous and next lines in coil filament input data text
+    file.
+
+    Arguments:
+        prev_line (str): line in input data file representing coordinates of
+            filament point previous to the current point.
+        next_line (str): line in input data file representing coordinates of
+            filament point next to the current point.
+
+    Returns:
+        tangent (array of float): tangent vector at filament point.
+    """
+    prev_columns = prev_line.strip().split()
+    prev_x = float(prev_columns[0])
+    prev_y = float(prev_columns[1])
+    prev_z = float(prev_columns[2])
+    prev_pt = np.array([prev_x, prev_y, prev_z])
+
+    next_columns = next_line.strip().split()
+    next_x = float(next_columns[0])
+    next_y = float(next_columns[1])
+    next_z = float(next_columns[2])
+    next_pt = np.array([next_x, next_y, next_z])
+
+    tangent = next_pt - prev_pt
+    tangent = normalize(tangent)
+
+    return tangent
+
+
 class MagnetSet(object):
     """An object representing a set of modular stellarator magnet coils.
 
     Arguments:
         coils_file (str): path to coil filament data file.
-        cross_section (list): coil cross-section definiton. The cross-section
-            shape must be either a circle or rectangle. For a circular
-            cross-section, the list format is
-            ['circle' (str), radius [cm](float)]
-            For a rectangular cross-section, the list format is
-            ['rectangle' (str), width [cm](float), thickness [cm](float)]
+        width (float): width of coil cross-section in toroidal direction [cm].
+        thickness (float): thickness of coil cross-section in radial direction
+            [cm].
         toroidal_extent (float): toroidal extent to model [deg].
         logger (object): logger object (optional, defaults to None). If no
             logger is supplied, a default logger will be instantiated.
@@ -38,12 +69,13 @@ class MagnetSet(object):
     """
 
     def __init__(
-        self, coils_file, cross_section, toroidal_extent, logger=None, **kwargs
+        self, coils_file, width, thickness, toroidal_extent, logger=None, **kwargs
     ):
 
         self.logger = logger
         self.coils_file = coils_file
-        self.cross_section = cross_section
+        self.width = width
+        self.thickness = thickness
         self.toroidal_extent = toroidal_extent
 
         self.start_line = 3
@@ -59,16 +91,32 @@ def __init__(
         ):
             self.__setattr__(name, kwargs[name])
 
-        cubit_io.init_cubit()
+        # Define maximum length of coil cross-section
+        self.max_cs_len = max(self._width, self._thickness)
 
     @property
-    def cross_section(self):
-        return self._cross_section
+    def width(self):
+        return self._width
+
+    @width.setter
+    def width(self, value):
+        self._width = value
+        if self._width < 0.0:
+            e = ValueError("Coil cross-section width cannot be negative.")
+            self._logger.error(e.args[0])
+            raise e
 
-    @cross_section.setter
-    def cross_section(self, shape):
-        self._cross_section = shape
-        self._extract_cross_section()
+    @property
+    def thickness(self):
+        return self._thickness
+
+    @thickness.setter
+    def thickness(self, value):
+        self._thickness = value
+        if self._thickness < 0.0:
+            e = ValueError("Coil cross-section thickness cannot be negative.")
+            self._logger.error(e.args[0])
+            raise e
 
     @property
     def toroidal_extent(self):
@@ -90,247 +138,164 @@ def logger(self):
     def logger(self, logger_object):
         self._logger = log.check_init(logger_object)
 
-    def _extract_cross_section(self):
-        """Extract coil cross-section parameters.
-        (Internal function not intended to be called externally)
-
-        Arguments:
-            cross_section (list or tuple of str, float, float): coil
-                cross-section definition. Note that the cross-section shape
-                must be either a circle or rectangle.
-                For a circular cross-section, the list format is
-                ['circle' (str), radius (float, cm)]
-                For a rectangular cross-section, the list format is
-                ['rectangle' (str), width (float, cm), thickness (float, cm)]
-            logger (object): logger object.
-
-        Returns:
-            shape (str): cross-section shape.
-            shape_str (str): string to pass to Cubit for cross-section
-                generation. For a circular cross-section, the string format is
-                '{shape} radius {radius}'
-                For a rectangular cross-section, the string format is
-                '{shape} width {thickness} height {width}'
-            mag_len (float): characteristic length of magnets.
-        """
-        # Extract coil cross-section shape
-        shape = self._cross_section[0]
-
-        # Conditionally extract parameters for circular cross-section
-        if shape == "circle":
-            # Check that list format is correct
-            if len(self._cross_section) == 1:
-                e = ValueError(
-                    "Format of list defining circular cross-section must be\n"
-                    '["circle" (str), radius (float, cm)]'
-                )
-                self._logger.error(e.args[0])
-                raise e
-            elif len(self._cross_section) > 2:
-                w = Warning(
-                    "More than one length dimension has been defined for "
-                    "cross_section. Interpreting the first as the circle's"
-                    'radius; did you mean to use "rectangle"?'
-                )
-                self._logger.warning(w.args[0])
-
-            # Extract parameters
-            mag_len = self._cross_section[1]
-            # Define string to pass to Cubit for cross-section generation
-            shape_str = f"{shape} radius {mag_len}"
-        # Conditinally extract parameters for rectangular cross-section
-        elif shape == "rectangle":
-            # Check that list format is correct
-            if len(self._cross_section) != 3:
-                e = ValueError(
-                    "Format of list defining rectangular cross-section must \n"
-                    'be ["rectangle" (str), width (float, cm), thickness '
-                    "(float, cm)]"
-                )
-                self._logger.error(e.args[0])
-                raise e
-            # Extract parameters
-            width = self._cross_section[1]
-            thickness = self._cross_section[2]
-            # Detemine largest parameter
-            mag_len = max(width, thickness)
-            # Define string to pass to Cubit for cross-section generation
-            shape_str = f"{shape} width {thickness} height {width}"
-        # Otherwise, if input string is neither 'circle' nor 'rectangle',
-        #  raise an exception
-        else:
-            e = ValueError(
-                "Magnet cross-section must be either a circle or rectangle. "
-                "The first entry of the list defining the cross-section must be"
-                " the shape, with the following entries defining the shape"
-                "parameters.\n"
-                "\n"
-                "For a circular cross-section, the list format is\n"
-                '["circle" (str), radius (float, cm)]\n'
-                "\n"
-                "For a rectangular cross-section, the list format is\n"
-                '["rectangle" (str), width (float, cm),'
-                "thickness (float, cm)]"
-            )
-            self._logger.error(e.args[0])
-            raise e
-
-        self.shape = shape
-        self.shape_str = shape_str
-        self.mag_len = mag_len
-
-    def _extract_filaments(self):
-        """Extracts filament data from magnet coil data file.
+    def _extract_filament_data(self):
+        """Extracts filament coordinate data from input data file, sampling
+        filament point-loci according to the input sampling modifier and
+        computing tangents at each sampled point.
         (Internal function not intended to be called externally)
         """
         with open(self.coils_file, "r") as file:
             data = file.readlines()[self.start_line :]
 
         coords = []
-        filaments = []
+        filament_coords = []
+        tangents = []
+        filament_tangents = []
 
         # Ensure that sampling always starts on the first line of each filament
         sample_counter = 0
 
-        for line in data:
+        for i, line in enumerate(data):
             columns = line.strip().split()
 
             if columns[0] == "end":
                 break
 
-            x = float(columns[0]) * self.scale
-            y = float(columns[1]) * self.scale
-            z = float(columns[2]) * self.scale
-
             # Coil current
             s = float(columns[3])
 
             # s = 0 signals end of filament
             if s != 0:
                 if sample_counter % self.sample_mod == 0:
+                    x = float(columns[0]) * self.scale
+                    y = float(columns[1]) * self.scale
+                    z = float(columns[2]) * self.scale
                     coords.append([x, y, z])
+
+                    # Compute tangent
+                    if sample_counter == 0:
+                        # To compute tangent at initial point, store
+                        # corresponding "next" point
+                        next_line_init = data[i + 1]
+                    else:
+                        prev_line = data[i - 1]
+                        next_line = data[i + 1]
+                        tangent = compute_tangent(prev_line, next_line)
+                        tangents.append(tangent)
+
                 sample_counter += 1
             else:
-                coords.append([x, y, z])
-                filaments.append(np.array(coords))
-                sample_counter = 0
+                coords.append(coords[0])
+                filament_coords.append(np.array(coords))
                 coords.clear()
 
-        self.filaments = filaments
+                # To compute tangent at initial point, store point "previous"
+                # to final point (initial and final points are the same since
+                # filaments are closed loops)
+                prev_line_init = data[i - 1]
+                tangent = compute_tangent(prev_line_init, next_line_init)
+                tangents.insert(0, tangent)
+                tangents.append(tangent)
+                filament_tangents.append(np.array(tangents))
+                tangents.clear()
+
+                sample_counter = 0
 
-    def _set_average_radial_distance(self):
-        """Computes average radial distance of filament points.
+        self.filament_coords = filament_coords
+        self.filament_tangents = filament_tangents
+
+    def _compute_radial_distance_data(self):
+        """Computes average and maximum radial distance of filament points.
         (Internal function not intended to be called externally)
         """
         radial_distance = []
-        for f in self.filaments:
+
+        for f in self.filament_coords:
             radial_distance.extend(list(np.linalg.norm(f[:, :2], axis=1)))
+
         self.average_radial_distance = np.average(radial_distance)
+        self.max_radial_distance = np.max(radial_distance)
 
-    def _set_filtered_filaments(self):
+    def _filter_filaments(self):
         """Cleans filament data such that only filaments within the toroidal
         extent of the model are included and filaments are sorted by toroidal
         angle.
         (Internal function not intended to be called externally)
-
-        Arguments:
-            filaments (np array of list of list of float): list of filament
-                coordinates. Each filament is a list of coordinates.
-            r_avg (float): average radial distance of magnets (cm).
-            mag_len (float): characteristic length of magnets.
-
-        Returns:
-            filtered_filaments (list of list of list of float): sorted list
-            of filament coordinates.
         """
         # Initialize data for filaments within toroidal extent of model
-        reduced_fils = []
-        # Initialize list of filament centers of mass for those within toroidal
-        # extent of model
+        filtered_coords = []
+        filtered_tangents = []
+        # Initialize list of filament centers of mass
         com_list = []
 
-        # Define tolerance of toroidal extent to account for width of coils
+        # Define tolerance of toroidal extent to account for dimensionality of
+        # coil cross-section
         # Multiply by factor of 2 to be conservative
-        tol = 2 * np.arctan2(self.mag_len, self.average_radial_distance)
+        tol = 2 * np.arctan2(self.max_cs_len, self.average_radial_distance)
 
         # Compute lower and upper bounds of toroidal extent within tolerance
-        min_rad = 2 * np.pi - tol
-        max_rad = self._toroidal_extent + tol
+        lower_bound = 2 * np.pi - tol
+        upper_bound = self._toroidal_extent + tol
 
-        for fil in self.filaments:
+        for coords, tangents in zip(self.filament_coords, self.filament_tangents):
             # Compute filament center of mass
-            com = np.average(fil, axis=0)
+            com = np.average(coords, axis=0)
             # Compute toroidal angle of each point in filament
-            phi_pts = np.arctan2(fil[:, 1], fil[:, 0])
+            toroidal_angles = np.arctan2(coords[:, 1], coords[:, 0])
             # Ensure angles are positive
-            phi_pts = (phi_pts + 2 * np.pi) % (2 * np.pi)
+            toroidal_angles = (toroidal_angles + 2 * np.pi) % (2 * np.pi)
             # Compute bounds of toroidal extent of filament
-            min_phi = np.min(phi_pts)
-            max_phi = np.max(phi_pts)
+            min_tor_ang = np.min(toroidal_angles)
+            max_tor_ang = np.max(toroidal_angles)
 
             # Determine if filament toroidal extent overlaps with that of model
-            if (min_phi >= min_rad or min_phi <= max_rad) or (
-                max_phi >= min_rad or max_phi <= max_rad
+            if (min_tor_ang >= lower_bound or min_tor_ang <= upper_bound) or (
+                max_tor_ang >= lower_bound or max_tor_ang <= upper_bound
             ):
-                reduced_fils.append(fil)
+                filtered_coords.append(coords)
+                filtered_tangents.append(tangents)
                 com_list.append(com)
 
-        reduced_fils = np.array(reduced_fils)
+        filtered_coords = np.array(filtered_coords)
+        filtered_tangents = np.array(filtered_tangents)
         com_list = np.array(com_list)
 
         # Compute toroidal angles of filament centers of mass
-        phi_arr = np.arctan2(com_list[:, 1], com_list[:, 0])
-        phi_arr = (phi_arr + 2 * np.pi) % (2 * np.pi)
+        com_toroidal_angles = np.arctan2(com_list[:, 1], com_list[:, 0])
+        com_toroidal_angles = (com_toroidal_angles + 2 * np.pi) % (2 * np.pi)
 
-        # Sort filaments by toroidal angle
-        self.filtered_filaments = [
-            x for _, x in sorted(zip(phi_arr, reduced_fils))
-        ]
+        # Sort filaments by toroidal angle and overwrite respective arrays
+        self.filament_coords = np.array(
+            [x for _, x in sorted(zip(com_toroidal_angles, filtered_coords))]
+        )
+        self.filament_tangents = np.array(
+            [x for _, x in sorted(zip(com_toroidal_angles, filtered_tangents))]
+        )
+        self.filament_com = np.array(
+            [x for _, x in sorted(zip(com_toroidal_angles, com_list))]
+        )
 
-    def _cut_magnets(self, volume_ids):
-        """Cleanly cuts the magnets at the planes defining the toriodal extent.
+    def _cut_magnets(self):
+        """Cuts the magnets at the planes defining the toriodal extent.
         (Internal function not intended to be called externally)
-
-        Arguments:
-            volume_ids (list): volume ids corresponding to each magnet volume
-
-        Returns:
-            volume_ids (range): new volume ids corresponding to magnet volumes
-                following cutting operation
         """
-        # Define sweeping surface width
-        # Multiply by factor of 2 to be conservative
-        rec_width = 2 * self.average_radial_distance
-
-        cubit.cmd(f"create surface rectangle width {rec_width} yplane")
-        surf_id = cubit.get_last_id("surface")
-
-        # Shift surface to positive x axis
-        cubit.cmd(f"move Surface {surf_id} x {rec_width/2}")
-
-        # Revolve surface to create wedge spanning toroidal extent
-        cubit.cmd(
-            (
-                f"sweep surface {surf_id} zaxis angle "
-                f"{np.rad2deg(self._toroidal_extent)}"
-            )
+        toroidal_region = cq.Workplane("XZ")
+        toroidal_region = toroidal_region.transformed(
+            offset=(1.25 * self.max_radial_distance / 2, 0)
         )
-        sweep_id = cubit.get_last_id("volume")
-
-        # Remove magnets and magnet portions not within toroidal extent
-        cubit.cmd(
-            "intersect volume "
-            + " ".join(str(i) for i in volume_ids)
-            + f" {sweep_id}"
+        toroidal_region = toroidal_region.rect(
+            1.25 * self.max_radial_distance, 1.25 * self.max_radial_distance
         )
+        toroidal_region = toroidal_region.revolve(
+            np.rad2deg(self._toroidal_extent),
+            (-1.25 * self.max_radial_distance / 2, 0),
+            (-1.25 * self.max_radial_distance / 2, 1),
+        )
+        toroidal_region = toroidal_region.val()
 
-        # Renumber volume ids from 1 to N
-        cubit.cmd("compress all")
-
-        # Extract new volume ids
-        volume_ids = cubit.get_entities("volume")
-
-        return volume_ids
+        for coil in self.magnet_coils:
+            cut_coil = coil.solid.intersect(toroidal_region)
+            coil.solid = cut_coil
 
     def build_magnet_coils(self):
         """Builds each filament in self.filtered_filaments in cubit, then cuts
@@ -338,44 +303,50 @@ def build_magnet_coils(self):
         """
         self._logger.info("Constructing magnet coils...")
 
-        self._extract_filaments()
-        self._set_average_radial_distance()
-        self._set_filtered_filaments()
+        self._extract_filament_data()
+        self._compute_radial_distance_data()
+        self._filter_filaments()
 
         self.magnet_coils = [
-            MagnetCoil(filament, self.shape, self.shape_str)
-            for filament in self.filtered_filaments
+            MagnetCoil(coords, tangents, center_of_mass, self._width, self._thickness)
+            for coords, tangents, center_of_mass in zip(
+                self.filament_coords, self.filament_tangents, self.filament_com
+            )
         ]
 
-        volume_ids = []
-
-        for magnet_coil in self.magnet_coils:
-            volume_id = magnet_coil.create_magnet()
-            volume_ids.append(volume_id)
+        [magnet_coil.create_magnet() for magnet_coil in self.magnet_coils]
 
-        volume_ids = self._cut_magnets(volume_ids)
+        self._cut_magnets()
 
-        self.volume_ids = volume_ids
-
-    def export_step(self, step_filename="magnets", export_dir=""):
-        """Export CAD solids as a STEP file via Coreform Cubit.
+    def export_step(self, step_filename="magnet_set", export_dir=""):
+        """Export CAD solids as a STEP file via CadQuery.
 
         Arguments:
             step_filename (str): name of STEP output file, excluding '.step'
-                extension (optional, defaults to 'magnets').
+                extension (optional, defaults to 'magnet_set').
             export_dir (str): directory to which to export the STEP output file
                 (optional, defaults to empty string).
         """
         self._logger.info("Exporting STEP file for magnet coils...")
 
-        cubit_io.export_step_cubit(
-            filename=step_filename, export_dir=export_dir
+        self.export_dir = export_dir
+        self.step_filename = step_filename
+
+        export_path = Path(self.export_dir) / Path(self.step_filename).with_suffix(
+            ".step"
         )
 
+        coil_set = cq.Compound.makeCompound([coil.solid for coil in self.magnet_coils])
+        cq.exporters.export(coil_set, str(export_path))
+
     def mesh_magnets(self):
         """Creates tetrahedral mesh of magnet volumes via Coreform Cubit."""
         self._logger.info("Generating tetrahedral mesh of magnet coils...")
 
+        last_vol_id = cubit_io.import_step_cubit(self.step_filename, self.export_dir)
+
+        self.volume_ids = range(1, last_vol_id + 1)
+
         for vol in self.volume_ids:
             cubit.cmd(f"volume {vol} scheme tetmesh")
             cubit.cmd(f"mesh volume {vol}")
@@ -392,173 +363,95 @@ def export_mesh(self, mesh_filename="magnet_mesh", export_dir=""):
         """
         self._logger.info("Exporting mesh H5M file for magnet coils...")
 
-        cubit_io.export_mesh_cubit(
-            filename=mesh_filename, export_dir=export_dir
-        )
+        cubit_io.export_mesh_cubit(filename=mesh_filename, export_dir=export_dir)
 
 
 class MagnetCoil(object):
     """An object representing a single modular stellarator magnet coil.
 
     Arguments:
-        filament (np.ndarray(double)): set of Cartesian coordinates defining
+        coords (2-D array of float): set of Cartesian coordinates defining
             magnet filament location.
-        shape (str): shape of coil cross-section.
-        shape_str (str): string defining cross-section shape for Coreform Cubit.
+        tangents (2-D array of float): set of tangent vectors at each filament
+            location.
+        center_of_mass (1-D array of float): Cartesian coordinates of filament
+            center of mass.
+        width (float): width of coil cross-section in toroidal direction [cm].
+        thickness (float): thickness of coil cross-section in radial direction
+            [cm].
     """
 
-    def __init__(self, filament, shape, shape_str):
-
-        self.filament = filament
-        self.shape = shape
-        self.shape_str = shape_str
+    def __init__(self, coords, tangents, center_of_mass, width, thickness):
 
-    def _orient_rectangle(
-        self, path_origin, surf_id, t_vec, norm, rot_axis, rot_ang_norm
-    ):
-        """Orients rectangular cross-section in the normal plane such that its
-        thickness direction faces the origin.
-        (Internal function not intended to be called externally)
-
-        Arguments:
-            path_origin (int): index of initial point in filament path.
-            surf_id (int): index of cross-section surface.
-            t_vec (list of float): cross-section thickness vector.
-            norm (list of float): cross-section normal vector.
-            rot_axis (list of float): axis about which to rotate the
-                cross-section.
-            rot_ang_norm (float): angle by which cross-section was rotated to
-                align its normal with the initial point tangent (deg).
-        """
-        # Determine orientation of thickness vector after cross-section was
-        # oriented along filament origin tangent
-
-        # Compute part of thickness vector parallel to rotation axis
-        t_vec_par = normalize(np.inner(t_vec, rot_axis) * rot_axis)
-        # Compute part of thickness vector orthogonal to rotation axis
-        t_vec_perp = normalize(t_vec - t_vec_par)
-
-        # Compute vector othogonal to both rotation axis and orthogonal
-        # part of thickness vector
-        orth = normalize(np.cross(rot_axis, t_vec_perp))
-
-        # Determine part of rotated vector parallel to original
-        rot_par = np.cos(rot_ang_norm)
-        # Determine part of rotated vector orthogonal to original
-        rot_perp = np.sin(rot_ang_norm)
-
-        # Compute orthogonal part of thickness vector after rotation
-        t_vec_perp_rot = rot_par * t_vec_perp + rot_perp * orth
-        # Compute thickness vector after rotation
-        t_vec_rot = normalize(t_vec_perp_rot + t_vec_par)
-
-        # Orient cross-section in its plane such that it faces the global origin
-
-        # Extract initial path point
-        pos = cubit.vertex(path_origin).coordinates()
-
-        # Project position vector onto cross-section
-        pos_proj = normalize(pos - np.inner(pos, norm) * norm)
-
-        # Compute angle by which to rotate cross-section such that it faces the
-        # origin
-        rot_ang_orig = np.arccos(np.inner(pos_proj, t_vec_rot))
-
-        # Re-orient rotated cross-section such that thickness vector faces
-        # origin
-        cubit.cmd(
-            f"rotate Surface {surf_id} angle {np.rad2deg(rot_ang_orig)} about "
-            "origin 0 0 0 direction " + " ".join(str(i) for i in norm)
-        )
+        self.coords = coords
+        self.tangents = tangents
+        self.center_of_mass = center_of_mass
+        self.width = width
+        self.thickness = thickness
 
     def create_magnet(self):
-        """Creates magnet coil volumes in cubit.
+        """Creates a single magnet coil CAD solid in CadQuery.
 
         Returns:
-            volume_id (int): magnet volume ids in cubit
+            coil (object): cq.Solid object representing a single magnet coil.
         """
-        # Cross-section inititally populated with thickness vector
-        # oriented along x axis
-        t_vec = np.array([1, 0, 0])
+        tangent_vectors = [cq.Vector(tuple(tangent)) for tangent in self.tangents]
 
-        # Create cross-section for sweep
-        cubit.cmd(f"create surface " + self.shape_str + " zplane")
+        # Define coil filament path normals such that they face the filament
+        # center of mass
+        normal_dirs = np.array([i - self.center_of_mass for i in self.coords])
+        normal_dirs = normal_dirs / np.linalg.norm(normal_dirs, axis=1)[:, np.newaxis]
 
-        # Store cross-section index
-        cs_id = cubit.get_last_id("surface")
-        # Cross-section initially populated with normal oriented along z
-        # axis
-        cs_axis = np.array([0, 0, 1])
+        # Project normal directions onto desired coil cross-section (CS) plane
+        # at each filament position to define true filament path normals
+        parallel_parts = []
+        for dir, tangent in zip(normal_dirs, self.tangents):
+            parallel_parts.append(np.dot(dir, tangent))
+        parallel_parts = np.array(parallel_parts)
 
-        # Initialize path list
-        path = []
+        normals = normal_dirs - parallel_parts[:, np.newaxis] * self.tangents
+        normals = normals / np.linalg.norm(normals, axis=1)[:, np.newaxis]
 
-        # Create vertices in filament path
-        for x, y, z in self.filament:
-            cubit.cmd(f"create vertex {x} {y} {z}")
-            path += [cubit.get_last_id("vertex")]
+        # Compute binormals projected onto CS plane at each position
+        binormals = np.cross(self.tangents, normals)
 
-        # Ensure final vertex in path is the same as the first
-        path += [path[0]]
+        # Compute coordinates of edges of rectangular coils
+        coil_edge_coords = []
 
-        cubit.cmd(
-            f"create curve spline location vertex "
-            + " ".join(str(i) for i in path)
-        )
-        curve_id = cubit.get_last_id("curve")
-
-        # Define new surface normal vector as that pointing between path
-        # points adjacent to initial point
-
-        # Extract path points adjacent to initial point
-        next_pt = np.array(cubit.vertex(path[1]).coordinates())
-        last_pt = np.array(cubit.vertex(path[-2]).coordinates())
-        # Compute direction in which to align surface normal
-        tang = normalize(np.subtract(next_pt, last_pt))
-
-        # Define axis and angle of rotation to orient cross-section along
-        # defined normal
-
-        # Define axis of rotation as orthogonal to both z axis and surface
-        # normal
-        rot_axis = normalize(np.cross(cs_axis, tang))
-        # Compute angle by which to rotate cross-section to orient along
-        # defined surface normal
-        rot_ang_norm = np.arccos(np.inner(cs_axis, tang))
-
-        # Copy cross-section for sweep
-        cubit.cmd(f"surface {cs_id} copy")
-        surf_id = cubit.get_last_id("surface")
-
-        # Orient cross-section along defined normal
-        cubit.cmd(
-            f"rotate Surface {surf_id} angle {np.rad2deg(rot_ang_norm)} about "
-            "origin 0 0 0 direction " + " ".join(str(i) for i in rot_axis)
+        coil_edge1_coords = (
+            self.coords - (self.width / 2) * binormals - (self.thickness / 2) * normals
         )
+        coil_edge_coords.append([cq.Vector(tuple(pos)) for pos in coil_edge1_coords])
 
-        # Conditionally orients rectangular cross-section
-        if self.shape == "rectangle":
-            self._orient_rectangle(
-                path[0], surf_id, t_vec, tang, rot_axis, rot_ang_norm
-            )
+        coil_edge2_coords = (
+            self.coords - (self.width / 2) * binormals + (self.thickness / 2) * normals
+        )
+        coil_edge_coords.append([cq.Vector(tuple(pos)) for pos in coil_edge2_coords])
 
-        # Move cross-section to initial path point
-        cubit.cmd(f"move Surface {surf_id} location vertex {path[0]}")
+        coil_edge3_coords = (
+            self.coords + (self.width / 2) * binormals + (self.thickness / 2) * normals
+        )
+        coil_edge_coords.append([cq.Vector(tuple(pos)) for pos in coil_edge3_coords])
 
-        # Sweep cross-section to create magnet coil
-        cubit.cmd(
-            f"sweep surface {surf_id} along curve {curve_id} " f"individual"
+        coil_edge4_coords = (
+            self.coords + (self.width / 2) * binormals - (self.thickness / 2) * normals
         )
-        volume_id = cubit.get_last_id("volume")
+        coil_edge_coords.append([cq.Vector(tuple(pos)) for pos in coil_edge4_coords])
+        # Append first edge once again
+        coil_edge_coords.append([cq.Vector(tuple(pos)) for pos in coil_edge1_coords])
 
-        # Delete extraneous curves and vertices
-        cubit.cmd(f"delete curve {curve_id}")
-        cubit.cmd("delete vertex all")
+        coil_edges = [
+            cq.Edge.makeSpline(coord_vectors, tangents=tangent_vectors).close()
+            for coord_vectors in coil_edge_coords
+        ]
 
-        # Delete original cross-section
-        cubit.cmd(f"delete surface {cs_id}")
+        face_list = [
+            cq.Face.makeRuledSurface(edge1, edge2)
+            for edge1, edge2 in zip(coil_edges[:-1], coil_edges[1:])
+        ]
 
-        return volume_id
+        shell = cq.Shell.makeShell(face_list)
+        self.solid = cq.Solid.makeSolid(shell)
 
 
 def parse_args():
diff --git a/parastell/parastell.py b/parastell/parastell.py
index cb5d2bea..c3650b5f 100644
--- a/parastell/parastell.py
+++ b/parastell/parastell.py
@@ -193,18 +193,16 @@ def export_invessel_build(
             )
 
     def construct_magnets(
-        self, coils_file, cross_section, toroidal_extent, **kwargs
+        self, coils_file, width, thickness, toroidal_extent, **kwargs
     ):
         """Constructs MagnetSet class object.
 
         Arguments:
             coils_file (str): path to coil filament data file.
-            cross_section (list): coil cross-section definiton. The
-                cross-section shape must be either a circle or rectangle. For a
-                circular cross-section, the list format is
-                ['circle' (str), radius [cm](float)]
-                For a rectangular cross-section, the list format is
-                ['rectangle' (str), width [cm](float), thickness [cm](float)]
+            width (float): width of coil cross-section in toroidal direction
+                [cm].
+            thickness (float): thickness of coil cross-section in radial
+                direction [cm].
             toroidal_extent (float): toroidal extent to model [deg].
 
         Optional attributes:
@@ -219,7 +217,8 @@ def construct_magnets(
         """
         self.magnet_set = mc.MagnetSet(
             coils_file,
-            cross_section,
+            width,
+            thickness,
             toroidal_extent,
             logger=self._logger,
             **kwargs,
@@ -229,7 +228,7 @@ def construct_magnets(
 
     def export_magnets(
         self,
-        step_filename="magnets",
+        step_filename="magnet_set",
         export_mesh=False,
         mesh_filename="magnet_mesh",
         export_dir="",
@@ -238,7 +237,8 @@ def export_magnets(
 
         Arguments:
             step_filename (str): name of STEP export output file, excluding
-                '.step' extension (optional, optional, defaults to 'magnets').
+                '.step' extension (optional, optional, defaults to
+                'magnet_set').
             export_mesh (bool): flag to indicate tetrahedral mesh generation
                 for magnet volumes (optional, defaults to False).
             mesh_filename (str): name of tetrahedral mesh H5M file, excluding
@@ -246,9 +246,7 @@ def export_magnets(
             export_dir (str): directory to which to export output files
                 (optional, defaults to empty string).
         """
-        self.magnet_set.export_step(
-            step_filename=step_filename, export_dir=export_dir
-        )
+        self.magnet_set.export_step(step_filename=step_filename, export_dir=export_dir)
 
         if export_mesh:
             self.magnet_set.mesh_magnets()
@@ -311,29 +309,31 @@ def _import_ivb_step(self):
             name,
             data,
         ) in self.invessel_build.radial_build.radial_build.items():
-            vol_id = cubit_io.import_step_cubit(
-                name, self.invessel_build.export_dir
-            )
+            vol_id = cubit_io.import_step_cubit(name, self.invessel_build.export_dir)
             data["vol_id"] = vol_id
 
+    def _import_magnets_step(self):
+        """Import STEP file for magnet set into Coreform Cubit.
+        (Internal function not intended to be called externally)
+        """
+        last_vol_id = cubit_io.import_step_cubit(
+            self.magnet_set.step_filename, self.magnet_set.export_dir
+        )
+
+        self.magnet_set.volume_ids = range(1, last_vol_id + 1)
+
     def _tag_materials_legacy(self):
         """Applies material tags to corresponding CAD volumes for legacy DAGMC
         neutronics model export.
         (Internal function not intended to be called externally)
         """
         if self.magnet_set:
-            vol_id_str = " ".join(
-                str(i) for i in list(self.magnet_set.volume_ids)
-            )
-            cubit.cmd(
-                f'group "mat:{self.magnet_set.mat_tag}" add volume {vol_id_str}'
-            )
+            vol_id_str = " ".join(str(i) for i in list(self.magnet_set.volume_ids))
+            cubit.cmd(f'group "mat:{self.magnet_set.mat_tag}" add volume {vol_id_str}')
 
         if self.invessel_build:
             for data in self.invessel_build.radial_build.radial_build.values():
-                cubit.cmd(
-                    f'group "mat:{data["mat_tag"]}" add volume {data["vol_id"]}'
-                )
+                cubit.cmd(f'group "mat:{data["mat_tag"]}" add volume {data["vol_id"]}')
 
     def _tag_materials_native(self):
         """Applies material tags to corresponding CAD volumes for native DAGMC
@@ -367,9 +367,15 @@ def build_cubit_model(self, skip_imprint=False, legacy_faceting=True):
         """
         self.legacy_faceting = legacy_faceting
 
-        self._logger.info(
-            "Building DAGMC neutronics model via Coreform Cubit..."
-        )
+        self._logger.info("Building DAGMC neutronics model via Coreform Cubit...")
+
+        if cubit_io.initialized:
+            cubit.cmd("new")
+        else:
+            cubit_io.init_cubit()
+
+        if self.magnet_set:
+            self._import_magnets_step()
 
         if self.invessel_build:
             self._import_ivb_step()
@@ -466,8 +472,7 @@ def parse_args():
         "--logger",
         action="store_true",
         help=(
-            "flag to indicate the instantiation of a logger object (default: "
-            "False)"
+            "flag to indicate the instantiation of a logger object (default: " "False)"
         ),
     )
 
@@ -485,9 +490,7 @@ def parse_args():
         "-m",
         "--magnets",
         action="store_true",
-        help=(
-            "flag to indicate the creation of magnet geometry (default: False)"
-        ),
+        help=("flag to indicate the creation of magnet geometry (default: False)"),
     )
 
     parser.add_argument(
@@ -513,9 +516,7 @@ def parse_args():
     return parser.parse_args()
 
 
-def check_inputs(
-    invessel_build, magnet_coils, source_mesh, dagmc_export, logger
-):
+def check_inputs(invessel_build, magnet_coils, source_mesh, dagmc_export, logger):
     """Checks inputs for consistency across ParaStell classes.
 
     Arguments:
diff --git a/tests/test_magnet_coils.py b/tests/test_magnet_coils.py
index 99d999a7..57887aa8 100644
--- a/tests/test_magnet_coils.py
+++ b/tests/test_magnet_coils.py
@@ -1,107 +1,80 @@
 from pathlib import Path
 
 import pytest
+import numpy as np
 
 import parastell.magnet_coils as magnet_coils
 
 
 def remove_files():
 
-    if Path("magnets.step").exists():
-        Path.unlink("magnets.step")
+    if Path("magnet_set.step").exists():
+        Path.unlink("magnet_set.step")
     if Path("magnet_mesh.exo").exists():
         Path.unlink("magnet_mesh.exo")
     if Path("magnet_mesh.h5m").exists():
         Path.unlink("magnet_mesh.h5m")
     if Path("stellarator.log").exists():
         Path.unlink("stellarator.log")
-    if Path("step_export.log").exists():
-        Path.unlink("step_export.log")
+    if Path("step_import.log").exists():
+        Path.unlink("step_import.log")
 
 
 @pytest.fixture
-def rect_coil_set():
+def coil_set():
 
     coils_file = Path("files_for_tests") / "coils.example"
-    rect_cross_section = ["rectangle", 20, 60]
+    width = 40.0
+    thickness = 50.0
     toroidal_extent = 90.0
-    sample_mod = 6
+    sample_mod = 10
 
-    rect_coil_obj = magnet_coils.MagnetSet(
-        coils_file, rect_cross_section, toroidal_extent, sample_mod=sample_mod
+    coil_set_obj = magnet_coils.MagnetSet(
+        coils_file, width, thickness, toroidal_extent, sample_mod=sample_mod
     )
 
-    return rect_coil_obj
+    return coil_set_obj
 
 
-@pytest.fixture
-def circ_coil_set():
-
-    coils_file = Path("files_for_tests") / "coils.example"
-    circ_cross_section = ["circle", 25]
-    toroidal_extent = 90.0
-    sample_mod = 6
-
-    circ_coil_obj = magnet_coils.MagnetSet(
-        coils_file, circ_cross_section, toroidal_extent, sample_mod=sample_mod
-    )
-
-    return circ_coil_obj
-
-
-def test_rectangular_magnets(rect_coil_set):
-
-    shape_exp = "rectangle"
-    shape_str_exp = "rectangle width 60 height 20"
-    mag_len_exp = 60
-
-    remove_files()
-
-    assert rect_coil_set.shape == shape_exp
-    assert rect_coil_set.shape_str == shape_str_exp
-    assert rect_coil_set.mag_len == mag_len_exp
-
-    remove_files()
-
-
-def test_circular_magnets(circ_coil_set):
+def test_magnet_construction(coil_set):
 
-    len_filaments_exp = 2
-    average_radial_distance_exp = 1068.3010006669892
-    len_filtered_filaments_exp = 1
-    shape_exp = "circle"
-    shape_str_exp = "circle radius 25"
-    mag_len_exp = 25
-    len_test_coil_filament_exp = 23
+    width_exp = 40.0
+    thickness_exp = 50.0
+    toroidal_extent_exp = np.deg2rad(90.0)
+    max_cs_len_exp = 50.0
+    average_radial_distance_exp = 1075.9559376755847
+    max_radial_distance_exp = 1641.2890540431476
+    len_filaments_exp = 1
+    len_coords_exp = 14
 
     remove_files()
 
-    circ_coil_set.build_magnet_coils()
+    coil_set.build_magnet_coils()
 
-    assert len(circ_coil_set.filaments) == len_filaments_exp
-    assert circ_coil_set.average_radial_distance == average_radial_distance_exp
-    assert len(circ_coil_set.filtered_filaments) == len_filtered_filaments_exp
-    assert circ_coil_set.shape == shape_exp
-    assert circ_coil_set.shape_str == shape_str_exp
-    assert circ_coil_set.mag_len == mag_len_exp
+    assert len(coil_set.filament_coords) == len_filaments_exp
+    assert coil_set.width == width_exp
+    assert coil_set.thickness == thickness_exp
+    assert coil_set.toroidal_extent == toroidal_extent_exp
+    assert coil_set.max_cs_len == max_cs_len_exp
+    assert coil_set.average_radial_distance == average_radial_distance_exp
+    assert coil_set.max_radial_distance == max_radial_distance_exp
 
-    test_coil = circ_coil_set.magnet_coils[0]
-    test_coil_filament = test_coil.filament
-    assert len(test_coil_filament) == len_test_coil_filament_exp
+    test_coil = coil_set.magnet_coils[0]
+    assert len(test_coil.coords) == len_coords_exp
 
     remove_files()
 
 
-def test_magnet_exports(circ_coil_set):
+def test_magnet_exports(coil_set):
 
     remove_files()
 
-    circ_coil_set.build_magnet_coils()
-    circ_coil_set.export_step()
-    assert Path("magnets.step").exists()
+    coil_set.build_magnet_coils()
+    coil_set.export_step()
+    assert Path("magnet_set.step").exists()
 
-    circ_coil_set.mesh_magnets()
-    circ_coil_set.export_mesh()
+    coil_set.mesh_magnets()
+    coil_set.export_mesh()
     assert Path("magnet_mesh.h5m").exists()
 
     remove_files()
diff --git a/tests/test_parastell.py b/tests/test_parastell.py
index 90bc7343..5def3c9a 100644
--- a/tests/test_parastell.py
+++ b/tests/test_parastell.py
@@ -12,8 +12,8 @@ def remove_files():
         Path.unlink("chamber.step")
     if Path("component.step").exists():
         Path.unlink("component.step")
-    if Path("magnets.step").exists():
-        Path.unlink("magnets.step")
+    if Path("magnet_set.step").exists():
+        Path.unlink("magnet_set.step")
     if Path("magnet_mesh.exo").exists():
         Path.unlink("magnet_mesh.exo")
     if Path("magnet_mesh.h5m").exists():
@@ -54,9 +54,7 @@ def test_parastell(stellarator):
     component_name_exp = "component"
     radial_build_dict = {
         component_name_exp: {
-            "thickness_matrix": np.ones(
-                (len(toroidal_angles), len(poloidal_angles))
-            )
+            "thickness_matrix": np.ones((len(toroidal_angles), len(poloidal_angles)))
             * 10
         }
     }
@@ -81,15 +79,16 @@ def test_parastell(stellarator):
     # Magnet Coils
 
     coils_file = Path("files_for_tests") / "coils.example"
-    cross_section = ["circle", 25]
+    width = 40.0
+    thickness = 50.0
     toroidal_extent = 90.0
     sample_mod = 6
 
     stellarator.construct_magnets(
-        coils_file, cross_section, toroidal_extent, sample_mod=sample_mod
+        coils_file, width, thickness, toroidal_extent, sample_mod=sample_mod
     )
 
-    step_filename_exp = "magnets"
+    step_filename_exp = "magnet_set"
     export_mesh = True
     mesh_filename_exp = "magnet_mesh"
 

From 38abfb82aa337e12f48e9818ee85265245ee42d2 Mon Sep 17 00:00:00 2001
From: Connor Moreno <camoreno@Denali.ep.wisc.edu>
Date: Tue, 27 Aug 2024 13:55:02 -0500
Subject: [PATCH 2/7] Fix black formatting

---
 parastell/magnet_coils.py | 78 +++++++++++++++++++++++++++++----------
 parastell/parastell.py    | 35 +++++++++++++-----
 tests/test_parastell.py   |  4 +-
 3 files changed, 87 insertions(+), 30 deletions(-)

diff --git a/parastell/magnet_coils.py b/parastell/magnet_coils.py
index 4ac72947..cb2bd3fb 100644
--- a/parastell/magnet_coils.py
+++ b/parastell/magnet_coils.py
@@ -69,7 +69,13 @@ class MagnetSet(object):
     """
 
     def __init__(
-        self, coils_file, width, thickness, toroidal_extent, logger=None, **kwargs
+        self,
+        coils_file,
+        width,
+        thickness,
+        toroidal_extent,
+        logger=None,
+        **kwargs,
     ):
 
         self.logger = logger
@@ -237,7 +243,9 @@ def _filter_filaments(self):
         lower_bound = 2 * np.pi - tol
         upper_bound = self._toroidal_extent + tol
 
-        for coords, tangents in zip(self.filament_coords, self.filament_tangents):
+        for coords, tangents in zip(
+            self.filament_coords, self.filament_tangents
+        ):
             # Compute filament center of mass
             com = np.average(coords, axis=0)
             # Compute toroidal angle of each point in filament
@@ -308,7 +316,9 @@ def build_magnet_coils(self):
         self._filter_filaments()
 
         self.magnet_coils = [
-            MagnetCoil(coords, tangents, center_of_mass, self._width, self._thickness)
+            MagnetCoil(
+                coords, tangents, center_of_mass, self._width, self._thickness
+            )
             for coords, tangents, center_of_mass in zip(
                 self.filament_coords, self.filament_tangents, self.filament_com
             )
@@ -332,18 +342,22 @@ def export_step(self, step_filename="magnet_set", export_dir=""):
         self.export_dir = export_dir
         self.step_filename = step_filename
 
-        export_path = Path(self.export_dir) / Path(self.step_filename).with_suffix(
-            ".step"
-        )
+        export_path = Path(self.export_dir) / Path(
+            self.step_filename
+        ).with_suffix(".step")
 
-        coil_set = cq.Compound.makeCompound([coil.solid for coil in self.magnet_coils])
+        coil_set = cq.Compound.makeCompound(
+            [coil.solid for coil in self.magnet_coils]
+        )
         cq.exporters.export(coil_set, str(export_path))
 
     def mesh_magnets(self):
         """Creates tetrahedral mesh of magnet volumes via Coreform Cubit."""
         self._logger.info("Generating tetrahedral mesh of magnet coils...")
 
-        last_vol_id = cubit_io.import_step_cubit(self.step_filename, self.export_dir)
+        last_vol_id = cubit_io.import_step_cubit(
+            self.step_filename, self.export_dir
+        )
 
         self.volume_ids = range(1, last_vol_id + 1)
 
@@ -363,7 +377,9 @@ def export_mesh(self, mesh_filename="magnet_mesh", export_dir=""):
         """
         self._logger.info("Exporting mesh H5M file for magnet coils...")
 
-        cubit_io.export_mesh_cubit(filename=mesh_filename, export_dir=export_dir)
+        cubit_io.export_mesh_cubit(
+            filename=mesh_filename, export_dir=export_dir
+        )
 
 
 class MagnetCoil(object):
@@ -395,12 +411,16 @@ def create_magnet(self):
         Returns:
             coil (object): cq.Solid object representing a single magnet coil.
         """
-        tangent_vectors = [cq.Vector(tuple(tangent)) for tangent in self.tangents]
+        tangent_vectors = [
+            cq.Vector(tuple(tangent)) for tangent in self.tangents
+        ]
 
         # Define coil filament path normals such that they face the filament
         # center of mass
         normal_dirs = np.array([i - self.center_of_mass for i in self.coords])
-        normal_dirs = normal_dirs / np.linalg.norm(normal_dirs, axis=1)[:, np.newaxis]
+        normal_dirs = (
+            normal_dirs / np.linalg.norm(normal_dirs, axis=1)[:, np.newaxis]
+        )
 
         # Project normal directions onto desired coil cross-section (CS) plane
         # at each filament position to define true filament path normals
@@ -419,26 +439,44 @@ def create_magnet(self):
         coil_edge_coords = []
 
         coil_edge1_coords = (
-            self.coords - (self.width / 2) * binormals - (self.thickness / 2) * normals
+            self.coords
+            - (self.width / 2) * binormals
+            - (self.thickness / 2) * normals
+        )
+        coil_edge_coords.append(
+            [cq.Vector(tuple(pos)) for pos in coil_edge1_coords]
         )
-        coil_edge_coords.append([cq.Vector(tuple(pos)) for pos in coil_edge1_coords])
 
         coil_edge2_coords = (
-            self.coords - (self.width / 2) * binormals + (self.thickness / 2) * normals
+            self.coords
+            - (self.width / 2) * binormals
+            + (self.thickness / 2) * normals
+        )
+        coil_edge_coords.append(
+            [cq.Vector(tuple(pos)) for pos in coil_edge2_coords]
         )
-        coil_edge_coords.append([cq.Vector(tuple(pos)) for pos in coil_edge2_coords])
 
         coil_edge3_coords = (
-            self.coords + (self.width / 2) * binormals + (self.thickness / 2) * normals
+            self.coords
+            + (self.width / 2) * binormals
+            + (self.thickness / 2) * normals
+        )
+        coil_edge_coords.append(
+            [cq.Vector(tuple(pos)) for pos in coil_edge3_coords]
         )
-        coil_edge_coords.append([cq.Vector(tuple(pos)) for pos in coil_edge3_coords])
 
         coil_edge4_coords = (
-            self.coords + (self.width / 2) * binormals - (self.thickness / 2) * normals
+            self.coords
+            + (self.width / 2) * binormals
+            - (self.thickness / 2) * normals
+        )
+        coil_edge_coords.append(
+            [cq.Vector(tuple(pos)) for pos in coil_edge4_coords]
         )
-        coil_edge_coords.append([cq.Vector(tuple(pos)) for pos in coil_edge4_coords])
         # Append first edge once again
-        coil_edge_coords.append([cq.Vector(tuple(pos)) for pos in coil_edge1_coords])
+        coil_edge_coords.append(
+            [cq.Vector(tuple(pos)) for pos in coil_edge1_coords]
+        )
 
         coil_edges = [
             cq.Edge.makeSpline(coord_vectors, tangents=tangent_vectors).close()
diff --git a/parastell/parastell.py b/parastell/parastell.py
index c3650b5f..fcfbd83f 100644
--- a/parastell/parastell.py
+++ b/parastell/parastell.py
@@ -246,7 +246,9 @@ def export_magnets(
             export_dir (str): directory to which to export output files
                 (optional, defaults to empty string).
         """
-        self.magnet_set.export_step(step_filename=step_filename, export_dir=export_dir)
+        self.magnet_set.export_step(
+            step_filename=step_filename, export_dir=export_dir
+        )
 
         if export_mesh:
             self.magnet_set.mesh_magnets()
@@ -309,7 +311,9 @@ def _import_ivb_step(self):
             name,
             data,
         ) in self.invessel_build.radial_build.radial_build.items():
-            vol_id = cubit_io.import_step_cubit(name, self.invessel_build.export_dir)
+            vol_id = cubit_io.import_step_cubit(
+                name, self.invessel_build.export_dir
+            )
             data["vol_id"] = vol_id
 
     def _import_magnets_step(self):
@@ -328,12 +332,18 @@ def _tag_materials_legacy(self):
         (Internal function not intended to be called externally)
         """
         if self.magnet_set:
-            vol_id_str = " ".join(str(i) for i in list(self.magnet_set.volume_ids))
-            cubit.cmd(f'group "mat:{self.magnet_set.mat_tag}" add volume {vol_id_str}')
+            vol_id_str = " ".join(
+                str(i) for i in list(self.magnet_set.volume_ids)
+            )
+            cubit.cmd(
+                f'group "mat:{self.magnet_set.mat_tag}" add volume {vol_id_str}'
+            )
 
         if self.invessel_build:
             for data in self.invessel_build.radial_build.radial_build.values():
-                cubit.cmd(f'group "mat:{data["mat_tag"]}" add volume {data["vol_id"]}')
+                cubit.cmd(
+                    f'group "mat:{data["mat_tag"]}" add volume {data["vol_id"]}'
+                )
 
     def _tag_materials_native(self):
         """Applies material tags to corresponding CAD volumes for native DAGMC
@@ -367,7 +377,9 @@ def build_cubit_model(self, skip_imprint=False, legacy_faceting=True):
         """
         self.legacy_faceting = legacy_faceting
 
-        self._logger.info("Building DAGMC neutronics model via Coreform Cubit...")
+        self._logger.info(
+            "Building DAGMC neutronics model via Coreform Cubit..."
+        )
 
         if cubit_io.initialized:
             cubit.cmd("new")
@@ -472,7 +484,8 @@ def parse_args():
         "--logger",
         action="store_true",
         help=(
-            "flag to indicate the instantiation of a logger object (default: " "False)"
+            "flag to indicate the instantiation of a logger object (default: "
+            "False)"
         ),
     )
 
@@ -490,7 +503,9 @@ def parse_args():
         "-m",
         "--magnets",
         action="store_true",
-        help=("flag to indicate the creation of magnet geometry (default: False)"),
+        help=(
+            "flag to indicate the creation of magnet geometry (default: False)"
+        ),
     )
 
     parser.add_argument(
@@ -516,7 +531,9 @@ def parse_args():
     return parser.parse_args()
 
 
-def check_inputs(invessel_build, magnet_coils, source_mesh, dagmc_export, logger):
+def check_inputs(
+    invessel_build, magnet_coils, source_mesh, dagmc_export, logger
+):
     """Checks inputs for consistency across ParaStell classes.
 
     Arguments:
diff --git a/tests/test_parastell.py b/tests/test_parastell.py
index 5def3c9a..77acbe2a 100644
--- a/tests/test_parastell.py
+++ b/tests/test_parastell.py
@@ -54,7 +54,9 @@ def test_parastell(stellarator):
     component_name_exp = "component"
     radial_build_dict = {
         component_name_exp: {
-            "thickness_matrix": np.ones((len(toroidal_angles), len(poloidal_angles)))
+            "thickness_matrix": np.ones(
+                (len(toroidal_angles), len(poloidal_angles))
+            )
             * 10
         }
     }

From aaa0d51962036434316921ae09c75cd4d36cc70e Mon Sep 17 00:00:00 2001
From: Connor Moreno <camoreno@Denali.ep.wisc.edu>
Date: Wed, 4 Sep 2024 13:00:35 -0500
Subject: [PATCH 3/7] Modify data structure design

---
 parastell/magnet_coils.py  | 167 +++++++++++++------------------------
 tests/test_magnet_coils.py |   4 +-
 2 files changed, 60 insertions(+), 111 deletions(-)

diff --git a/parastell/magnet_coils.py b/parastell/magnet_coils.py
index cb2bd3fb..779e8b0b 100644
--- a/parastell/magnet_coils.py
+++ b/parastell/magnet_coils.py
@@ -28,16 +28,10 @@ def compute_tangent(prev_line, next_line):
         tangent (array of float): tangent vector at filament point.
     """
     prev_columns = prev_line.strip().split()
-    prev_x = float(prev_columns[0])
-    prev_y = float(prev_columns[1])
-    prev_z = float(prev_columns[2])
-    prev_pt = np.array([prev_x, prev_y, prev_z])
+    prev_pt = np.array([float(coord) for coord in prev_columns[0:3]])
 
     next_columns = next_line.strip().split()
-    next_x = float(next_columns[0])
-    next_y = float(next_columns[1])
-    next_z = float(next_columns[2])
-    next_pt = np.array([next_x, next_y, next_z])
+    next_pt = np.array([float(coord) for coord in next_columns[0:3]])
 
     tangent = next_pt - prev_pt
     tangent = normalize(tangent)
@@ -145,9 +139,7 @@ def logger(self, logger_object):
         self._logger = log.check_init(logger_object)
 
     def _extract_filament_data(self):
-        """Extracts filament coordinate data from input data file, sampling
-        filament point-loci according to the input sampling modifier and
-        computing tangents at each sampled point.
+        """Extracts filament coordinate numerical data from input data file.
         (Internal function not intended to be called externally)
         """
         with open(self.coils_file, "r") as file:
@@ -155,13 +147,8 @@ def _extract_filament_data(self):
 
         coords = []
         filament_coords = []
-        tangents = []
-        filament_tangents = []
 
-        # Ensure that sampling always starts on the first line of each filament
-        sample_counter = 0
-
-        for i, line in enumerate(data):
+        for line in data:
             columns = line.strip().split()
 
             if columns[0] == "end":
@@ -172,43 +159,17 @@ def _extract_filament_data(self):
 
             # s = 0 signals end of filament
             if s != 0:
-                if sample_counter % self.sample_mod == 0:
-                    x = float(columns[0]) * self.scale
-                    y = float(columns[1]) * self.scale
-                    z = float(columns[2]) * self.scale
-                    coords.append([x, y, z])
-
-                    # Compute tangent
-                    if sample_counter == 0:
-                        # To compute tangent at initial point, store
-                        # corresponding "next" point
-                        next_line_init = data[i + 1]
-                    else:
-                        prev_line = data[i - 1]
-                        next_line = data[i + 1]
-                        tangent = compute_tangent(prev_line, next_line)
-                        tangents.append(tangent)
-
-                sample_counter += 1
+                x = float(columns[0]) * self.scale
+                y = float(columns[1]) * self.scale
+                z = float(columns[2]) * self.scale
+                coords.append([x, y, z])
+            
             else:
                 coords.append(coords[0])
                 filament_coords.append(np.array(coords))
                 coords.clear()
 
-                # To compute tangent at initial point, store point "previous"
-                # to final point (initial and final points are the same since
-                # filaments are closed loops)
-                prev_line_init = data[i - 1]
-                tangent = compute_tangent(prev_line_init, next_line_init)
-                tangents.insert(0, tangent)
-                tangents.append(tangent)
-                filament_tangents.append(np.array(tangents))
-                tangents.clear()
-
-                sample_counter = 0
-
         self.filament_coords = filament_coords
-        self.filament_tangents = filament_tangents
 
     def _compute_radial_distance_data(self):
         """Computes average and maximum radial distance of filament points.
@@ -228,9 +189,9 @@ def _filter_filaments(self):
         angle.
         (Internal function not intended to be called externally)
         """
-        # Initialize data for filaments within toroidal extent of model
+        # Initialize coordinate data for filaments within toroidal extent of
+        # model
         filtered_coords = []
-        filtered_tangents = []
         # Initialize list of filament centers of mass
         com_list = []
 
@@ -243,9 +204,7 @@ def _filter_filaments(self):
         lower_bound = 2 * np.pi - tol
         upper_bound = self._toroidal_extent + tol
 
-        for coords, tangents in zip(
-            self.filament_coords, self.filament_tangents
-        ):
+        for coords in self.filament_coords:
             # Compute filament center of mass
             com = np.average(coords, axis=0)
             # Compute toroidal angle of each point in filament
@@ -261,24 +220,20 @@ def _filter_filaments(self):
                 max_tor_ang >= lower_bound or max_tor_ang <= upper_bound
             ):
                 filtered_coords.append(coords)
-                filtered_tangents.append(tangents)
                 com_list.append(com)
 
         filtered_coords = np.array(filtered_coords)
-        filtered_tangents = np.array(filtered_tangents)
         com_list = np.array(com_list)
 
         # Compute toroidal angles of filament centers of mass
         com_toroidal_angles = np.arctan2(com_list[:, 1], com_list[:, 0])
         com_toroidal_angles = (com_toroidal_angles + 2 * np.pi) % (2 * np.pi)
 
-        # Sort filaments by toroidal angle and overwrite respective arrays
+        # Sort filaments by toroidal angle and overwrite coordinate data
         self.filament_coords = np.array(
             [x for _, x in sorted(zip(com_toroidal_angles, filtered_coords))]
         )
-        self.filament_tangents = np.array(
-            [x for _, x in sorted(zip(com_toroidal_angles, filtered_tangents))]
-        )
+        # Store filtered centers of mass since they have already been computed
         self.filament_com = np.array(
             [x for _, x in sorted(zip(com_toroidal_angles, com_list))]
         )
@@ -287,17 +242,17 @@ def _cut_magnets(self):
         """Cuts the magnets at the planes defining the toriodal extent.
         (Internal function not intended to be called externally)
         """
+        side_length = 1.25 * self.max_radial_distance
+        
         toroidal_region = cq.Workplane("XZ")
-        toroidal_region = toroidal_region.transformed(
-            offset=(1.25 * self.max_radial_distance / 2, 0)
-        )
-        toroidal_region = toroidal_region.rect(
-            1.25 * self.max_radial_distance, 1.25 * self.max_radial_distance
+        toroidal_region = toroidal_region.transformed(offset=(
+            side_length / 2, 0)
         )
+        toroidal_region = toroidal_region.rect(side_length, side_length)
         toroidal_region = toroidal_region.revolve(
             np.rad2deg(self._toroidal_extent),
-            (-1.25 * self.max_radial_distance / 2, 0),
-            (-1.25 * self.max_radial_distance / 2, 1),
+            (-side_length / 2, 0),
+            (-side_length / 2, 1),
         )
         toroidal_region = toroidal_region.val()
 
@@ -317,10 +272,11 @@ def build_magnet_coils(self):
 
         self.magnet_coils = [
             MagnetCoil(
-                coords, tangents, center_of_mass, self._width, self._thickness
+                coords, center_of_mass, self._width, self._thickness,
+                self.sample_mod
             )
-            for coords, tangents, center_of_mass in zip(
-                self.filament_coords, self.filament_tangents, self.filament_com
+            for coords, center_of_mass in zip(
+                self.filament_coords, self.filament_com
             )
         ]
 
@@ -397,14 +353,32 @@ class MagnetCoil(object):
             [cm].
     """
 
-    def __init__(self, coords, tangents, center_of_mass, width, thickness):
+    def __init__(self, coords, center_of_mass, width, thickness, sample_mod):
 
+        self.sample_mod = sample_mod
         self.coords = coords
-        self.tangents = tangents
         self.center_of_mass = center_of_mass
         self.width = width
         self.thickness = thickness
 
+    @property
+    def coords(self):
+        return self._coords
+
+    @coords.setter
+    def coords(self, data):
+        # Compute tangents
+        tangents = np.subtract(
+            np.append(data[1:], [data[1]], axis=0),
+            np.append([data[-2]], data[0:-1], axis=0)
+        )
+        
+        # Sample filament coordinates and tangents by modifier
+        self._coords = data[0:-1:self.sample_mod]
+        self._coords = np.append(self._coords, [self._coords[0]], axis=0)
+        self.tangents = tangents[0:-1:self.sample_mod]
+        self.tangents = np.append(self.tangents, [self.tangents[0]], axis=0)
+
     def create_magnet(self):
         """Creates a single magnet coil CAD solid in CadQuery.
 
@@ -417,7 +391,7 @@ def create_magnet(self):
 
         # Define coil filament path normals such that they face the filament
         # center of mass
-        normal_dirs = np.array([i - self.center_of_mass for i in self.coords])
+        normal_dirs = self._coords - self.center_of_mass
         normal_dirs = (
             normal_dirs / np.linalg.norm(normal_dirs, axis=1)[:, np.newaxis]
         )
@@ -436,47 +410,22 @@ def create_magnet(self):
         binormals = np.cross(self.tangents, normals)
 
         # Compute coordinates of edges of rectangular coils
-        coil_edge_coords = []
-
-        coil_edge1_coords = (
-            self.coords
-            - (self.width / 2) * binormals
-            - (self.thickness / 2) * normals
-        )
-        coil_edge_coords.append(
-            [cq.Vector(tuple(pos)) for pos in coil_edge1_coords]
-        )
+        edge_offsets = np.array([[-1, -1], [-1, 1], [1, 1], [1, -1]])
 
-        coil_edge2_coords = (
-            self.coords
-            - (self.width / 2) * binormals
-            + (self.thickness / 2) * normals
-        )
-        coil_edge_coords.append(
-            [cq.Vector(tuple(pos)) for pos in coil_edge2_coords]
-        )
+        coil_edge_coords = []
+        for edge_offset in edge_offsets:
+            coil_edge = (
+                self._coords
+                + edge_offset[0] * binormals * (self.width / 2)
+                + edge_offset[1] * normals * (self.thickness / 2)
+            )
 
-        coil_edge3_coords = (
-            self.coords
-            + (self.width / 2) * binormals
-            + (self.thickness / 2) * normals
-        )
-        coil_edge_coords.append(
-            [cq.Vector(tuple(pos)) for pos in coil_edge3_coords]
-        )
+            coil_edge_coords.append(
+                [cq.Vector(tuple(pos)) for pos in coil_edge]
+            )            
 
-        coil_edge4_coords = (
-            self.coords
-            + (self.width / 2) * binormals
-            - (self.thickness / 2) * normals
-        )
-        coil_edge_coords.append(
-            [cq.Vector(tuple(pos)) for pos in coil_edge4_coords]
-        )
         # Append first edge once again
-        coil_edge_coords.append(
-            [cq.Vector(tuple(pos)) for pos in coil_edge1_coords]
-        )
+        coil_edge_coords.append(coil_edge_coords[0])
 
         coil_edges = [
             cq.Edge.makeSpline(coord_vectors, tangents=tangent_vectors).close()
diff --git a/tests/test_magnet_coils.py b/tests/test_magnet_coils.py
index 57887aa8..312eb337 100644
--- a/tests/test_magnet_coils.py
+++ b/tests/test_magnet_coils.py
@@ -42,8 +42,8 @@ def test_magnet_construction(coil_set):
     thickness_exp = 50.0
     toroidal_extent_exp = np.deg2rad(90.0)
     max_cs_len_exp = 50.0
-    average_radial_distance_exp = 1075.9559376755847
-    max_radial_distance_exp = 1641.2890540431476
+    average_radial_distance_exp = 1028.5044183872055
+    max_radial_distance_exp = 1646.3258131460148
     len_filaments_exp = 1
     len_coords_exp = 14
 

From 89892d39f2ff0e920bca2efa2e07a84f2614c49d Mon Sep 17 00:00:00 2001
From: Connor Moreno <camoreno@Denali.ep.wisc.edu>
Date: Wed, 4 Sep 2024 13:14:12 -0500
Subject: [PATCH 4/7] Ensure tangents are stored as unit vectors

---
 parastell/magnet_coils.py | 27 +--------------------------
 1 file changed, 1 insertion(+), 26 deletions(-)

diff --git a/parastell/magnet_coils.py b/parastell/magnet_coils.py
index 779e8b0b..206db6bb 100644
--- a/parastell/magnet_coils.py
+++ b/parastell/magnet_coils.py
@@ -13,32 +13,6 @@
 export_allowed_kwargs = ["step_filename", "export_mesh", "mesh_filename"]
 
 
-def compute_tangent(prev_line, next_line):
-    """Computes tangent at "current" filament point using central difference
-    approximation and previous and next lines in coil filament input data text
-    file.
-
-    Arguments:
-        prev_line (str): line in input data file representing coordinates of
-            filament point previous to the current point.
-        next_line (str): line in input data file representing coordinates of
-            filament point next to the current point.
-
-    Returns:
-        tangent (array of float): tangent vector at filament point.
-    """
-    prev_columns = prev_line.strip().split()
-    prev_pt = np.array([float(coord) for coord in prev_columns[0:3]])
-
-    next_columns = next_line.strip().split()
-    next_pt = np.array([float(coord) for coord in next_columns[0:3]])
-
-    tangent = next_pt - prev_pt
-    tangent = normalize(tangent)
-
-    return tangent
-
-
 class MagnetSet(object):
     """An object representing a set of modular stellarator magnet coils.
 
@@ -372,6 +346,7 @@ def coords(self, data):
             np.append(data[1:], [data[1]], axis=0),
             np.append([data[-2]], data[0:-1], axis=0)
         )
+        tangents = tangents/np.linalg.norm(tangents, axis=1)[:, np.newaxis]
         
         # Sample filament coordinates and tangents by modifier
         self._coords = data[0:-1:self.sample_mod]

From ed3db0fda9edefad06a8ce78210affbed06d0daa Mon Sep 17 00:00:00 2001
From: Connor Moreno <camoreno@Denali.ep.wisc.edu>
Date: Wed, 4 Sep 2024 13:15:38 -0500
Subject: [PATCH 5/7] Use black formatting

---
 parastell/magnet_coils.py | 27 +++++++++++++++------------
 1 file changed, 15 insertions(+), 12 deletions(-)

diff --git a/parastell/magnet_coils.py b/parastell/magnet_coils.py
index 206db6bb..9c390050 100644
--- a/parastell/magnet_coils.py
+++ b/parastell/magnet_coils.py
@@ -137,7 +137,7 @@ def _extract_filament_data(self):
                 y = float(columns[1]) * self.scale
                 z = float(columns[2]) * self.scale
                 coords.append([x, y, z])
-            
+
             else:
                 coords.append(coords[0])
                 filament_coords.append(np.array(coords))
@@ -217,10 +217,10 @@ def _cut_magnets(self):
         (Internal function not intended to be called externally)
         """
         side_length = 1.25 * self.max_radial_distance
-        
+
         toroidal_region = cq.Workplane("XZ")
-        toroidal_region = toroidal_region.transformed(offset=(
-            side_length / 2, 0)
+        toroidal_region = toroidal_region.transformed(
+            offset=(side_length / 2, 0)
         )
         toroidal_region = toroidal_region.rect(side_length, side_length)
         toroidal_region = toroidal_region.revolve(
@@ -246,8 +246,11 @@ def build_magnet_coils(self):
 
         self.magnet_coils = [
             MagnetCoil(
-                coords, center_of_mass, self._width, self._thickness,
-                self.sample_mod
+                coords,
+                center_of_mass,
+                self._width,
+                self._thickness,
+                self.sample_mod,
             )
             for coords, center_of_mass in zip(
                 self.filament_coords, self.filament_com
@@ -344,14 +347,14 @@ def coords(self, data):
         # Compute tangents
         tangents = np.subtract(
             np.append(data[1:], [data[1]], axis=0),
-            np.append([data[-2]], data[0:-1], axis=0)
+            np.append([data[-2]], data[0:-1], axis=0),
         )
-        tangents = tangents/np.linalg.norm(tangents, axis=1)[:, np.newaxis]
-        
+        tangents = tangents / np.linalg.norm(tangents, axis=1)[:, np.newaxis]
+
         # Sample filament coordinates and tangents by modifier
-        self._coords = data[0:-1:self.sample_mod]
+        self._coords = data[0 : -1 : self.sample_mod]
         self._coords = np.append(self._coords, [self._coords[0]], axis=0)
-        self.tangents = tangents[0:-1:self.sample_mod]
+        self.tangents = tangents[0 : -1 : self.sample_mod]
         self.tangents = np.append(self.tangents, [self.tangents[0]], axis=0)
 
     def create_magnet(self):
@@ -397,7 +400,7 @@ def create_magnet(self):
 
             coil_edge_coords.append(
                 [cq.Vector(tuple(pos)) for pos in coil_edge]
-            )            
+            )
 
         # Append first edge once again
         coil_edge_coords.append(coil_edge_coords[0])

From cd32c3c80a4480bd36516e605434888dc80fa85d Mon Sep 17 00:00:00 2001
From: Connor Moreno <camoreno@Denali.ep.wisc.edu>
Date: Mon, 9 Sep 2024 09:33:05 -0500
Subject: [PATCH 6/7] Streamline code

---
 parastell/magnet_coils.py | 41 +++++++++++++++++++++++++--------------
 1 file changed, 26 insertions(+), 15 deletions(-)

diff --git a/parastell/magnet_coils.py b/parastell/magnet_coils.py
index 9c390050..2e7a13fd 100644
--- a/parastell/magnet_coils.py
+++ b/parastell/magnet_coils.py
@@ -133,10 +133,9 @@ def _extract_filament_data(self):
 
             # s = 0 signals end of filament
             if s != 0:
-                x = float(columns[0]) * self.scale
-                y = float(columns[1]) * self.scale
-                z = float(columns[2]) * self.scale
-                coords.append([x, y, z])
+                coords.append(
+                    [float(ord) * self.scale for ord in columns[0:3]]
+                )
 
             else:
                 coords.append(coords[0])
@@ -149,13 +148,19 @@ def _compute_radial_distance_data(self):
         """Computes average and maximum radial distance of filament points.
         (Internal function not intended to be called externally)
         """
-        radial_distance = []
+        radii_count = 0
+        self.average_radial_distance = 0
+        self.max_radial_distance = -1
 
         for f in self.filament_coords:
-            radial_distance.extend(list(np.linalg.norm(f[:, :2], axis=1)))
+            radii = np.linalg.norm(f[:, :2], axis=1)
+            radii_count += len(radii)
+            self.average_radial_distance += np.sum(radii)
+            self.max_radial_distance = max(
+                self.max_radial_distance, np.max(radii)
+            )
 
-        self.average_radial_distance = np.average(radial_distance)
-        self.max_radial_distance = np.max(radial_distance)
+        self.average_radial_distance /= radii_count
 
     def _filter_filaments(self):
         """Cleans filament data such that only filaments within the toroidal
@@ -369,19 +374,25 @@ def create_magnet(self):
 
         # Define coil filament path normals such that they face the filament
         # center of mass
-        normal_dirs = self._coords - self.center_of_mass
-        normal_dirs = (
-            normal_dirs / np.linalg.norm(normal_dirs, axis=1)[:, np.newaxis]
+        # Compute "outward" direction as difference between filament positions
+        # and filament center of mass
+        outward_dirs = self._coords - self.center_of_mass
+        outward_dirs = (
+            outward_dirs / np.linalg.norm(outward_dirs, axis=1)[:, np.newaxis]
         )
 
-        # Project normal directions onto desired coil cross-section (CS) plane
-        # at each filament position to define true filament path normals
+        # Project outward directions onto desired coil cross-section (CS) plane
+        # at each filament position to define filament path normals
         parallel_parts = []
-        for dir, tangent in zip(normal_dirs, self.tangents):
+        for dir, tangent in zip(outward_dirs, self.tangents):
             parallel_parts.append(np.dot(dir, tangent))
         parallel_parts = np.array(parallel_parts)
 
-        normals = normal_dirs - parallel_parts[:, np.newaxis] * self.tangents
+        parallel_parts = np.diagonal(
+            np.matmul(outward_dirs, self.tangents.transpose())
+        )
+
+        normals = outward_dirs - parallel_parts[:, np.newaxis] * self.tangents
         normals = normals / np.linalg.norm(normals, axis=1)[:, np.newaxis]
 
         # Compute binormals projected onto CS plane at each position

From 3fbecafc5a339d544ec4c57e9f59161124b60c3f Mon Sep 17 00:00:00 2001
From: Connor Moreno <camoreno@Denali.ep.wisc.edu>
Date: Mon, 9 Sep 2024 09:34:53 -0500
Subject: [PATCH 7/7] Streamline code

---
 parastell/magnet_coils.py | 5 -----
 1 file changed, 5 deletions(-)

diff --git a/parastell/magnet_coils.py b/parastell/magnet_coils.py
index 2e7a13fd..4201b7e6 100644
--- a/parastell/magnet_coils.py
+++ b/parastell/magnet_coils.py
@@ -383,11 +383,6 @@ def create_magnet(self):
 
         # Project outward directions onto desired coil cross-section (CS) plane
         # at each filament position to define filament path normals
-        parallel_parts = []
-        for dir, tangent in zip(outward_dirs, self.tangents):
-            parallel_parts.append(np.dot(dir, tangent))
-        parallel_parts = np.array(parallel_parts)
-
         parallel_parts = np.diagonal(
             np.matmul(outward_dirs, self.tangents.transpose())
         )