Repository: wtbarnes/synthesizAR

Issue:
Rename `Loop` to `Strand`
This is more consistent with the nomenclature that a "strand" is the fundamental, often unresolved building block of the solar atmosphere, constituting the largest cross-sectional area over which the plasma is isothermal while "loops" are observed field-aligned intensity enhancements which can be made up of many unresolved strands.

Discussion:

Patch:
diff --git a/.github/workflows/run-tests.yml b/.github/workflows/ci.yml
similarity index 100%
rename from .github/workflows/run-tests.yml
rename to .github/workflows/ci.yml
diff --git a/README.md b/README.md
index 00e236f..d212537 100644
--- a/README.md
+++ b/README.md
@@ -1,22 +1,22 @@
-<div align="center">
-  <img src="docs/_static/synthesizar_logo.png" width=55%><br>
-</div>
-
 # synthesizAR
+
 [![Powered by SunPy Badge]( http://img.shields.io/badge/powered%20by-SunPy-orange.svg?style=flat)](http://www.sunpy.org)
-![synthesizAR CI status](https://github.com/wtbarnes/synthesizAR/workflows/Run%20tests/badge.svg)
+[![synthesizAR CI status](https://github.com/wtbarnes/synthesizAR/actions/workflows/ci.yml/badge.svg)](https://github.com/wtbarnes/synthesizAR/actions/workflows/ci.yml)
 [![Documentation Status](http://readthedocs.org/projects/synthesizar/badge/?version=latest)](http://synthesizar.readthedocs.io/en/latest/?badge=latest)
-[![codecov](https://codecov.io/gh/wtbarnes/synthesizAR/branch/master/graph/badge.svg)](https://codecov.io/gh/wtbarnes/synthesizAR)
+[![codecov](https://codecov.io/gh/wtbarnes/synthesizAR/graph/badge.svg?token=QUTQZGWQKB)](https://codecov.io/gh/wtbarnes/synthesizAR)
 
-synthesizAR is a Python package for forward modeling emission from solar active regions using hydrodynamic simulations of coronal loops
+synthesizAR is a Python package for forward modeling optically-thin emission from field-aligned hydrodynamic simulations of astrophysical plasmas.
 
 ## Install
+
 To clone the repository and install the package and all of its needed dependencies,
+
 ```shell
-$ git clone https://github.com/wtbarnes/synthesizAR.git
-$ cd synthesizAR
-$ pip install -e .[all]
+git clone https://github.com/wtbarnes/synthesizAR.git
+cd synthesizAR
+pip install -e .[all]
 ```
 
 ## Help
+
 See the [docs](http://synthesizar.readthedocs.io). To report bugs or request features, create an issue or submit a pull request.
diff --git a/examples/arcade-rtv.py b/examples/arcade-rtv.py
index 36200dd..22cc075 100644
--- a/examples/arcade-rtv.py
+++ b/examples/arcade-rtv.py
@@ -6,15 +6,17 @@ This example shows how to model AIA emission from an arcade
 of semi-circular loops who's thermal structure is modeled
 using the RTV scaling laws.
 """
-import astropy.units as u
 import astropy.time
+import astropy.units as u
+
 from astropy.coordinates import SkyCoord
 from sunpy.coordinates import get_earth
 
 import synthesizAR
-from synthesizAR.models import semi_circular_arcade
-from synthesizAR.interfaces import RTVInterface
+
 from synthesizAR.instruments import InstrumentSDOAIA
+from synthesizAR.interfaces import RTVInterface
+from synthesizAR.models import semi_circular_arcade
 
 #########################################################################
 # First, set up the coordinates for loops in the arcade.
@@ -24,7 +26,7 @@ arcade_coords = semi_circular_arcade(80*u.Mm, 5*u.deg, 50, pos, inclination=10*u
 
 #########################################################################
 # Next, assemble the arcade.
-strands = [synthesizAR.Loop(f'strand{i}', c) for i, c in enumerate(arcade_coords)]
+strands = [synthesizAR.Strand(f'strand{i}', c) for i, c in enumerate(arcade_coords)]
 arcade = synthesizAR.Skeleton(strands)
 
 #########################################################################
diff --git a/examples/loop-bundle-rtv.py b/examples/loop-bundle-rtv.py
index 23422c4..cfe48e2 100644
--- a/examples/loop-bundle-rtv.py
+++ b/examples/loop-bundle-rtv.py
@@ -6,17 +6,20 @@ This example shows how to model AIA emission from a
 bundle of semi-circular strands for two different
 viewpoints.
 """
-import matplotlib.pyplot as plt
 import astropy.time
 import astropy.units as u
-from astropy.visualization import quantity_support
+import matplotlib.pyplot as plt
+
 from astropy.coordinates import SkyCoord
+from astropy.visualization import quantity_support
 from sunpy.map import pixelate_coord_path, sample_at_coords
 
 import synthesizAR
+
 from synthesizAR.instruments import InstrumentSDOAIA
 from synthesizAR.interfaces import RTVInterface
 from synthesizAR.models import semi_circular_bundle
+
 # sphinx_gallery_thumbnail_number = -1
 
 ###########################################################################
@@ -32,7 +35,7 @@ bundle_coords = semi_circular_bundle(50 * u.Mm, 1*u.Mm, 500, observer=pos)
 ###########################################################################
 # As in other examples, we then use the coordinates of our strands to
 # construct the `~synthesizAR.Skeleton` object.
-strands = [synthesizAR.Loop(f'strand{i}', c) for i, c in enumerate(bundle_coords)]
+strands = [synthesizAR.Strand(f'strand{i}', c) for i, c in enumerate(bundle_coords)]
 bundle = synthesizAR.Skeleton(strands)
 bundle.peek(observer=pos)
 
diff --git a/examples/multi-instrument.py b/examples/multi-instrument.py
index 4f1ea92..e9d2f8b 100644
--- a/examples/multi-instrument.py
+++ b/examples/multi-instrument.py
@@ -9,16 +9,24 @@ It also demonstrates how to define a custom
 instrument class.
 """
 import astropy.units as u
-from astropy.coordinates import SkyCoord
-import numpy as np
 import matplotlib.pyplot as plt
-from sunpy.coordinates import get_earth, get_horizons_coord, HeliographicStonyhurst, Helioprojective
+import numpy as np
+
+from astropy.coordinates import SkyCoord
 from astropy.visualization import quantity_support
+from sunpy.coordinates import (
+    get_earth,
+    get_horizons_coord,
+    HeliographicStonyhurst,
+    Helioprojective,
+)
 
 import synthesizAR
-from synthesizAR.models import semi_circular_arcade
+
+from synthesizAR.instruments import InstrumentHinodeXRT, InstrumentSDOAIA
 from synthesizAR.interfaces import MartensInterface
-from synthesizAR.instruments import InstrumentSDOAIA, InstrumentHinodeXRT
+from synthesizAR.models import semi_circular_arcade
+
 # sphinx_gallery_thumbnail_number = -1
 
 ###############################################################################
@@ -27,7 +35,7 @@ from synthesizAR.instruments import InstrumentSDOAIA, InstrumentHinodeXRT
 obstime = '2021-10-28T15:00:00'
 loc = SkyCoord(HeliographicStonyhurst(lon=0*u.deg,lat=-30*u.deg, radius=1*u.R_sun, obstime=obstime))
 arcade = semi_circular_arcade(150*u.Mm, 10*u.deg, 50, loc, gamma=90*u.deg, n_points=5000)
-skeleton = synthesizAR.Skeleton([synthesizAR.Loop(f'{i}', c) for i,c in enumerate(arcade)])
+skeleton = synthesizAR.Skeleton([synthesizAR.Strand(f'{i}', c) for i,c in enumerate(arcade)])
 
 ###############################################################################
 # We'll select a few different observer locations for SDO and STERO-A and use
@@ -64,10 +72,10 @@ skeleton.load_loop_simulations(martens)
 with quantity_support():
     plt.figure(figsize=(11, 5))
     ax1 = plt.subplot(121)
-    for l in skeleton.loops:
+    for l in skeleton.strands:
         plt.plot(l.field_aligned_coordinate_center.to('Mm'), l.electron_temperature[0].to('MK'), color='k')
     plt.subplot(122)
-    for l in skeleton.loops:
+    for l in skeleton.strands:
         plt.plot(l.field_aligned_coordinate_center.to('Mm'), l.density[0], color='k')
     plt.yscale('log')
 
diff --git a/licenses/README.rst b/licenses/README.rst
deleted file mode 100644
index 67b82f6..0000000
--- a/licenses/README.rst
+++ /dev/null
@@ -1,9 +0,0 @@
-Licenses
-========
-
-This directory holds license and credit information for the package,
-works the package is derived from, and/or datasets.
-
-Ensure that you pick a package licence which is in this folder and it matches
-the one mentioned in the top level README.rst file. If you are using the
-pre-rendered version of this template check for the word 'Other' in the README.
diff --git a/licenses/TEMPLATE_LICENSE.rst b/licenses/TEMPLATE_LICENSE.rst
deleted file mode 100644
index f29177b..0000000
--- a/licenses/TEMPLATE_LICENSE.rst
+++ /dev/null
@@ -1,31 +0,0 @@
-This project is based upon the Astropy package template
-(https://github.com/astropy/package-template/) which is licenced under the terms
-of the following licence.
-
----
-
-Copyright (c) 2018, Astropy Developers
-All rights reserved.
-
-Redistribution and use in source and binary forms, with or without modification,
-are permitted provided that the following conditions are met:
-
-* Redistributions of source code must retain the above copyright notice, this
-  list of conditions and the following disclaimer.
-* Redistributions in binary form must reproduce the above copyright notice, this
-  list of conditions and the following disclaimer in the documentation and/or
-  other materials provided with the distribution.
-* Neither the name of the Astropy Team nor the names of its contributors may be
-  used to endorse or promote products derived from this software without
-  specific prior written permission.
-
-THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
-ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
-WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
-DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR
-ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
-(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
-LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON
-ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
-(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
-SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
diff --git a/pyproject.toml b/pyproject.toml
index bd3eb28..af8841e 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -11,7 +11,7 @@ authors = [
   {name="Will Barnes", email="will.t.barnes@gmail.com"}
 ]
 license = {file="LICENSE.rst"}
-description = "A Python package for forward-modeling optically-thin emission from field-aligned coronal loop models"
+description = "A Python package for forward modeling optically-thin emission from field-aligned hydrodynamic simulations of astrophysical plasmas."
 readme = {file="README.md", content-type = "text/markdown"}
 requires-python = ">=3.10"
 dependencies = [
diff --git a/synthesizAR/__init__.py b/synthesizAR/__init__.py
index 60b697a..bc65648 100644
--- a/synthesizAR/__init__.py
+++ b/synthesizAR/__init__.py
@@ -3,14 +3,15 @@ synthesizAR is a package for forward modeling emission from solar active regions
 hydrodynamic simulations of coronal loops.
 """
 try:
-    from .version import __version__
+    from synthesizAR.version import __version__
 except ImportError:
     __version__ = "unknown"
 
 
-from .loop import *
-from .skeleton import *
-
 # Set up logger
 from synthesizAR.util.logger import _init_log
+
+from .strand import *  # NOQA
+from .skeleton import *  # NOQA
+
 log = _init_log()
diff --git a/synthesizAR/instruments/base.py b/synthesizAR/instruments/base.py
index 65d6522..63ad4ee 100644
--- a/synthesizAR/instruments/base.py
+++ b/synthesizAR/instruments/base.py
@@ -1,21 +1,21 @@
 """
 Base class for instrument objects.
 """
+import astropy.units as u
 import copy
-import tempfile
+import numpy as np
 import pathlib
-from dataclasses import dataclass
+import tempfile
+import zarr
 
-import numpy as np
+from astropy.coordinates import SkyCoord
+from dataclasses import dataclass
 from scipy.interpolate import interp1d
 from scipy.ndimage import gaussian_filter
-import astropy.units as u
-from astropy.coordinates import SkyCoord
-from sunpy.coordinates.frames import Helioprojective, HeliographicStonyhurst
+from sunpy.coordinates.frames import HeliographicStonyhurst, Helioprojective
 from sunpy.map import make_fitswcs_header, Map
-import zarr
 
-from synthesizAR.util import is_visible, find_minimum_fov
+from synthesizAR.util import find_minimum_fov, is_visible
 from synthesizAR.util.decorators import return_quantity_as_tuple
 
 __all__ = ['ChannelBase', 'InstrumentBase']
@@ -27,7 +27,7 @@ class ChannelBase:
     channel: u.Quantity = None
 
 
-class InstrumentBase(object):
+class InstrumentBase:
     """
     Base class for instruments. This object is not meant to be instantiated directly. Instead,
     specific instruments should subclass this base object and implement a
@@ -170,17 +170,17 @@ class InstrumentBase(object):
             if client:
                 # Parallel
                 kernel_futures = client.map(self.calculate_intensity_kernel,
-                                            skeleton.loops,
+                                            skeleton.strands,
                                             channel=channel,
                                             **kwargs)
                 kernel_interp_futures = client.map(self.interpolate_to_instrument_time,
                                                    kernel_futures,
-                                                   skeleton.loops,
+                                                   skeleton.strands,
                                                    observing_time=(self.observing_time.value, self.observing_time.unit.to_string()))
             else:
                 # Serial
                 kernels_interp = []
-                for l in skeleton.loops:
+                for l in skeleton.strands:
                     k = self.calculate_intensity_kernel(l, channel=channel, **kwargs)
                     k = self.interpolate_to_instrument_time(
                         k, l, observing_time=(self.observing_time.value, self.observing_time.unit.to_string()),
@@ -189,12 +189,12 @@ class InstrumentBase(object):
 
             if kwargs.get('save_kernels_to_disk', False):
                 with tempfile.TemporaryDirectory() as tmpdir:
-                    self._make_stacked_kernel_array(tmpdir, skeleton.loops, channel)
-                    indices = self._find_loop_array_bounds(skeleton.loops)
+                    self._make_stacked_kernel_array(tmpdir, skeleton.strands, channel)
+                    indices = self._find_loop_array_bounds(skeleton.strands)
                     if client:
                         files = client.map(self.write_kernel_to_file,
                                            kernel_interp_futures,
-                                           skeleton.loops,
+                                           skeleton.strands,
                                            indices,
                                            channel=channel,
                                            name=self.name,
@@ -202,9 +202,9 @@ class InstrumentBase(object):
                         # NOTE: block here to avoid pileup of tasks that can overwhelm the scheduler
                         distributed.wait(files)
                     else:
-                        for k, l, i in zip(kernels_interp, skeleton.loops, indices):
+                        for k, l, i in zip(kernels_interp, skeleton.strands, indices):
                             self.write_kernel_to_file(k, l, i, channel, self.name, tmpdir)
-                    self._rechunk_stacked_kernels(tmpdir, skeleton.loops[0].model_results_filename, channel)
+                    self._rechunk_stacked_kernels(tmpdir, skeleton.strands[0].model_results_filename, channel)
                     kernels = self.observing_time.shape[0]*[None]  # placeholder so we know to read from a file
             else:
                 # NOTE: this can really blow up your memory if you are not careful
@@ -336,7 +336,7 @@ class InstrumentBase(object):
         # Compute weights
         if kernels is None:
             i_time = np.where(time == self.observing_time)[0][0]
-            root = skeleton.loops[0].zarr_root
+            root = skeleton.strands[0].zarr_root
             ds = root[f'{self.name}/{channel.name}_stacked_kernels']
             kernels = u.Quantity(ds[i_time, :], ds.attrs['unit'])
         # If a volumetric quantity, integrate over the cell and normalize by pixel area.
diff --git a/synthesizAR/instruments/hinode.py b/synthesizAR/instruments/hinode.py
index 92f2690..6014f3b 100644
--- a/synthesizAR/instruments/hinode.py
+++ b/synthesizAR/instruments/hinode.py
@@ -2,18 +2,18 @@
 Class for Hinode/EIS instrument. Holds information about spectral, temporal, and spatial resolution
 and other instrument-specific information.
 """
-from dataclasses import dataclass
+import astropy.constants as const
+import astropy.units as u
 import json
+import numpy as np
 import pkg_resources
+import xrtpy
 
-import numpy as np
+from dataclasses import dataclass
 from scipy.ndimage import gaussian_filter
-import astropy.units as u
-import astropy.constants as const
-import xrtpy
 
+from synthesizAR.instruments import ChannelBase, InstrumentBase
 from synthesizAR.util import SpatialPair
-from synthesizAR.instruments import InstrumentBase, ChannelBase
 from synthesizAR.util.decorators import return_quantity_as_tuple
 
 __all__ = ['InstrumentHinodeEIS', 'InstrumentHinodeXRT']
@@ -132,7 +132,7 @@ class InstrumentHinodeEIS(InstrumentBase):
         """
         hinode_fn = pkg_resources.resource_filename('synthesizAR',
                                                     'instruments/data/hinode_eis.json')
-        with open(hinode_fn, 'r') as f:
+        with open(hinode_fn) as f:
             eis_info = json.load(f)
 
         self.channels = []
@@ -170,7 +170,7 @@ class InstrumentHinodeEIS(InstrumentBase):
         """
         Build HDF5 files to store detector counts
         """
-        additional_fields = ['{}'.format(line.value) for line in field.loops[0].resolved_wavelengths]
+        additional_fields = [f'{line.value}' for line in field.strands[0].resolved_wavelengths]
         super().build_detector_file(file_template, chunks, additional_fields=additional_fields)
 
     def flatten(self, loop, interp_s, hf, start_index):
@@ -179,7 +179,7 @@ class InstrumentHinodeEIS(InstrumentBase):
         """
         for wavelength in loop.resolved_wavelengths:
             emiss, ion_name = loop.get_emission(wavelength, return_ion_name=True)
-            dset = hf['{}'.format(str(wavelength.value))]
+            dset = hf[f'{str(wavelength.value)}']
             dset.attrs['ion_name'] = ion_name
             self.interpolate_and_store(emiss, loop, interp_s, dset, start_index)
 
@@ -201,6 +201,7 @@ class InstrumentHinodeEIS(InstrumentBase):
         Calculate response of Hinode/EIS detector for given loop object.
         """
         import plasmapy
+
         # trim the instrument response to the appropriate wavelengths
         trimmed_indices = []
         for w in channel['model_wavelengths']:
@@ -215,7 +216,7 @@ class InstrumentHinodeEIS(InstrumentBase):
         counts = np.zeros(temperature.shape+response_x.shape)
         for wavelength in channel['model_wavelengths']:
             # thermal width + instrument width
-            ion_name = hf['{}'.format(str(wavelength.value))].attrs['ion_name']
+            ion_name = hf[f'{str(wavelength.value)}'].attrs['ion_name']
             ion_mass = plasmapy.atomic.ion_mass(ion_name.split(' ')[0].capitalize()).cgs
             thermal_velocity = 2.*const.k_B.cgs*temperature/ion_mass
             thermal_velocity = np.expand_dims(thermal_velocity, axis=2)*thermal_velocity.unit
@@ -225,14 +226,13 @@ class InstrumentHinodeEIS(InstrumentBase):
             doppler_shift = wavelength*los_velocity/const.c.cgs
             doppler_shift = np.expand_dims(doppler_shift, axis=2)*doppler_shift.unit
             # combine emissivity with instrument response function
-            dset = hf['{}'.format(str(wavelength.value))]
+            dset = hf[f'{str(wavelength.value)}']
             hist, edges = np.histogramdd(self.total_coordinates.value,
                                          bins=[self.bins.x, self.bins.y, self.bins.z],
                                          range=[self.bin_range.x, self.bin_range.y, self.bin_range.z],
                                          weights=np.array(dset[i_time, :]))
             emiss = np.dot(hist, np.diff(edges[2])).T
-            emiss = (np.expand_dims(emiss, axis=2)
-                     * u.Unit(get_keys(dset.attrs, ('unit', 'units')))*self.total_coordinates.unit)
+            emiss = np.expand_dims(emiss, axis=2) * u.Unit(dset.attrs.get('unit')) *self.total_coordinates.unit
             intensity = emiss*response_y/np.sqrt(2.*np.pi*line_width)
             intensity *= np.exp(-((response_x - wavelength - doppler_shift)**2)/(2.*line_width))
             if not hasattr(counts, 'unit'):
diff --git a/synthesizAR/interfaces/ebtel/ebtel.py b/synthesizAR/interfaces/ebtel/ebtel.py
index 68fe303..d975490 100644
--- a/synthesizAR/interfaces/ebtel/ebtel.py
+++ b/synthesizAR/interfaces/ebtel/ebtel.py
@@ -1,25 +1,26 @@
 """
 Interface between loop object and ebtel++ simulation
 """
-import os
+import astropy.units as u
 import copy
-import warnings
+import h5py
+import numpy as np
 import toolz
+import warnings
 
-import numpy as np
-import h5py
-import astropy.units as u
 from fiasco import Element
+
 try:
     import distributed
 except ImportError:
     warnings.warn('Dask library required for NEI calculation')
 
 from synthesizAR.atomic import non_equilibrium_ionization
+
 from .util import run_ebtel
 
 
-class EbtelInterface(object):
+class EbtelInterface:
     """
     Interface to the Enthalpy-Based Thermal Evolution of Loops (EBTEL) model
 
@@ -50,7 +51,7 @@ class EbtelInterface(object):
 
         Parameters
         ----------
-        loop : `synthesizAR.Loop` object
+        loop : `synthesizAR.Strand` object
         """
         # Configure run
         output_dict = copy.deepcopy(self.base_config)
@@ -122,8 +123,8 @@ class EbtelInterface(object):
             partial_nei = toolz.curry(EbtelInterface.compute_nei)(el)
             partial_write = toolz.curry(EbtelInterface.write_to_hdf5)(
                 element_name=el_name, savefile=emission_model.ionization_fraction_savefile)
-            y = client.map(partial_nei, skeleton.loops, pure=False)
-            write_y = client.map(partial_write, y, skeleton.loops, pure=False)
+            y = client.map(partial_nei, skeleton.strands, pure=False)
+            write_y = client.map(partial_write, y, skeleton.strands, pure=False)
             distributed.client.wait(write_y)
             futures[el_name] = write_y
 
diff --git a/synthesizAR/models/heating.py b/synthesizAR/models/heating.py
index d75930f..e76f32e 100644
--- a/synthesizAR/models/heating.py
+++ b/synthesizAR/models/heating.py
@@ -20,7 +20,7 @@ def b_over_l_scaling(loop, H_0=0.0738*u.Unit('erg cm-3 s-1'), alpha=0.3, beta=0.
 
     Parameters
     ----------
-    loop: `synthesizAR.Loop`
+    loop: `synthesizAR.Strand`
     H_0: `~astropy.units.Quantity`, optional
         Nominal heating rate.
     alpha: `float`, optional
diff --git a/synthesizAR/skeleton.py b/synthesizAR/skeleton.py
index bc5301b..756df6c 100644
--- a/synthesizAR/skeleton.py
+++ b/synthesizAR/skeleton.py
@@ -10,7 +10,7 @@ from astropy.coordinates import SkyCoord
 from functools import cached_property
 from scipy.interpolate import interp1d, splev, splprep
 
-from synthesizAR import Loop
+from synthesizAR import Strand
 from synthesizAR.visualize import plot_fieldlines
 
 __all__ = ['Skeleton']
@@ -22,19 +22,20 @@ class Skeleton:
 
     Parameters
     ----------
-    loops : `list`
-        List of `Loop` objects
+    strands : `list` of `~synthesizAR.Strand` objects
+        List of objects containing the information about each strand in the magnetic
+        skeleton.
 
     Examples
     --------
     >>> import synthesizAR
     >>> import astropy.units as u
-    >>> loop = synthesizAR.Loop('loop', SkyCoord(x=[1,4]*u.Mm, y=[2,5]*u.Mm, z=[3,6]*u.Mm,frame='heliographic_stonyhurst', representation_type='cartesian'), [1e2,1e3] * u.G)
-    >>> field = synthesizAR.Skeleton([loop,])
+    >>> strand = synthesizAR.Strand('strand', SkyCoord(x=[1,4]*u.Mm, y=[2,5]*u.Mm, z=[3,6]*u.Mm,frame='heliographic_stonyhurst', representation_type='cartesian'), [1e2,1e3] * u.G)
+    >>> field = synthesizAR.Skeleton([strand,])
     """
 
-    def __init__(self, loops):
-        self.loops = loops
+    def __init__(self, strands):
+        self.strands = strands
 
     @classmethod
     def from_coordinates(cls, coordinates, field_strengths=None, **kwargs):
@@ -43,29 +44,29 @@ class Skeleton:
 
         Parameters
         ----------
-        coordinates : `list`
-            List of `~astropy.coordinates.SkyCoord` loop coordinates
-        field_strengths : `list`
-            List of `~astropy.units.Quantity` scalar magnetic field strength along the loop
+        coordinates : `list` of `~astropy.coordinates.SkyCoord` objects
+            Coordinate of each strand
+        field_strengths : `list` of `~astropy.units.Quantity`
+            Scalar magnetic field strength along the strand
         """
-        loops = []
+        strands = []
         if field_strengths is None:
             field_strengths = len(coordinates) * [None]
         for i, (coord, fs) in enumerate(zip(coordinates, field_strengths)):
-            loops.append(Loop(f'loop{i:06d}', coord, field_strength=fs, **kwargs))
-        return cls(loops)
+            strands.append(Strand(f'strand{i:06d}', coord, field_strength=fs, **kwargs))
+        return cls(strands)
 
     def __repr__(self):
         return f'''synthesizAR Skeleton Object
 ------------------------
-Number of loops: {len(self.loops)}'''
+Number of strands: {len(self.strands)}'''
 
     def to_asdf(self, filename):
         """
         Serialize this instance of `Skeleton` to an ASDF file
         """
         tree = {}
-        for l in self.loops:
+        for l in self.strands:
             tree[l.name] = {
                 'field_strength': l.field_strength,
                 'coordinate': l.coordinate,
@@ -86,146 +87,146 @@ Number of loops: {len(self.loops)}'''
         >>> restored_field = synthesizAR.Skeleton.from_asdf('/path/to/skeleton.asdf') # doctest: +SKIP
         """
         exclude_keys = ['asdf_library', 'history']
-        loops = []
+        strands = []
         with asdf.open(filename, mode='r', memmap=False) as af:
             for k in af.keys():
                 if k in exclude_keys:
                     continue
                 model_results_filename = af.tree[k].get('model_results_filename', None)
                 cross_sectional_area = af.tree[k].get('cross_sectional_area', None)
-                loops.append(Loop(
+                strands.append(Strand(
                     k,
                     SkyCoord(af.tree[k]['coordinate']),
                     af.tree[k]['field_strength'],
                     cross_sectional_area=cross_sectional_area,
                     model_results_filename=model_results_filename,
                 ))
-        return cls(loops)
+        return cls(strands)
 
     @u.quantity_input
-    def refine_loops(self, delta_s: u.cm, **kwargs):
+    def refine_strands(self, delta_s: u.cm, **kwargs):
         """
-        Interpolate loop coordinates and field strengths to a specified spatial resolution
+        Interpolate strand coordinates and field strengths to a specified spatial resolution
         and return a new `Skeleton` object.
 
         This can be important in order to ensure that an adequate number of points are used
         to represent each fieldline when binning intensities onto the instrument grid.
         """
-        new_loops = []
-        for l in self.loops:
-            _l = self.refine_loop(l, delta_s, **kwargs)
-            new_loops.append(_l)
+        new_strands = []
+        for l in self.strands:
+            _l = self.refine_strand(l, delta_s, **kwargs)
+            new_strands.append(_l)
 
-        return Skeleton(new_loops)
+        return Skeleton(new_strands)
 
     @staticmethod
     @u.quantity_input
-    def refine_loop(loop, delta_s: u.cm, **kwargs):
+    def refine_strand(strand, delta_s: u.cm, **kwargs):
         evkwargs = kwargs.get('evkwargs', {})
         prepkwargs = kwargs.get('prepkwargs', {})
         try:
-            tck, _ = splprep(loop.coordinate.cartesian.xyz.value,
-                             u=loop.field_aligned_coordinate_norm, **prepkwargs)
-            new_s = np.arange(0, loop.length.to(u.Mm).value, delta_s.to(u.Mm).value) * u.Mm
-            x, y, z = splev((new_s/loop.length).decompose(), tck, **evkwargs)
+            tck, _ = splprep(strand.coordinate.cartesian.xyz.value,
+                             u=strand.field_aligned_coordinate_norm, **prepkwargs)
+            new_s = np.arange(0, strand.length.to(u.Mm).value, delta_s.to(u.Mm).value) * u.Mm
+            x, y, z = splev((new_s/strand.length).decompose(), tck, **evkwargs)
         except (ValueError, TypeError) as e:
-            raise Exception(f'Failed to refine {loop.name}') from e
-        unit = loop.coordinate.cartesian.xyz.unit
+            raise Exception(f'Failed to refine {strand.name}') from e
+        unit = strand.coordinate.cartesian.xyz.unit
         new_coord = SkyCoord(x=x*unit,
                              y=y*unit,
                              z=z*unit,
-                             frame=loop.coordinate.frame,
-                             representation_type=loop.coordinate.representation_type)
-        f_B = interp1d(loop.field_aligned_coordinate.to(u.Mm), loop.field_strength)
-        new_field_strength = f_B(new_s.to(u.Mm)) * loop.field_strength.unit
-        f_A = interp1d(loop.field_aligned_coordinate.to(u.Mm), loop.cross_sectional_area)
-        new_area = f_A(new_s.to(u.Mm)) * loop.cross_sectional_area.unit
-        return Loop(loop.name,
-                    new_coord,
-                    field_strength=new_field_strength,
-                    cross_sectional_area=new_area,
-                    model_results_filename=loop.model_results_filename)
+                             frame=strand.coordinate.frame,
+                             representation_type=strand.coordinate.representation_type)
+        f_B = interp1d(strand.field_aligned_coordinate.to(u.Mm), strand.field_strength)
+        new_field_strength = f_B(new_s.to(u.Mm)) * strand.field_strength.unit
+        f_A = interp1d(strand.field_aligned_coordinate.to(u.Mm), strand.cross_sectional_area)
+        new_area = f_A(new_s.to(u.Mm)) * strand.cross_sectional_area.unit
+        return Strand(strand.name,
+                      new_coord,
+                      field_strength=new_field_strength,
+                      cross_sectional_area=new_area,
+                      model_results_filename=strand.model_results_filename)
 
     @property
     def all_coordinates(self):
         """
-        Coordinates for all loops in the skeleton.
+        Coordinates for all stands in the skeleton.
 
         .. note:: This should be treated as a collection of points and NOT a
                   continuous structure.
         """
-        return SkyCoord([l.coordinate for l in self.loops],
-                        frame=self.loops[0].coordinate.frame,
-                        representation_type=self.loops[0].coordinate.representation_type)
+        return SkyCoord([l.coordinate for l in self.strands],
+                        frame=self.strands[0].coordinate.frame,
+                        representation_type=self.strands[0].coordinate.representation_type)
 
     @property
     def all_coordinates_centers(self):
         """
-        Coordinates for all grid cell centers of all loops in the skeleton
+        Coordinates for all grid cell centers of all strands in the skeleton
 
         .. note:: This should be treated as a collection of points and NOT a
                   continuous structure.
         """
-        return SkyCoord([l.coordinate_center for l in self.loops],
-                        frame=self.loops[0].coordinate_center.frame,
-                        representation_type=self.loops[0].coordinate_center.representation_type)
+        return SkyCoord([l.coordinate_center for l in self.strands],
+                        frame=self.strands[0].coordinate_center.frame,
+                        representation_type=self.strands[0].coordinate_center.representation_type)
 
     @cached_property
     def all_widths(self) -> u.cm:
         """
-        Widths for all loops concatenated together
+        Widths for all strands concatenated together
         """
-        return np.concatenate([l.field_aligned_coordinate_width for l in self.loops])
+        return np.concatenate([l.field_aligned_coordinate_width for l in self.strands])
 
     @cached_property
     def all_cross_sectional_areas(self) -> u.cm**2:
         """
-        Cross-sectional areas for all loops concatenated together.
+        Cross-sectional areas for all strands concatenated together.
 
-        .. note:: These are the cross-sectional areas evaluated at the center of the loop.
+        .. note:: These are the cross-sectional areas evaluated at the center of the strand.
         """
-        return np.concatenate([l.cross_sectional_area_center for l in self.loops])
+        return np.concatenate([l.cross_sectional_area_center for l in self.strands])
 
     def peek(self, **kwargs):
         """
-        Plot loop coordinates on the solar disk.
+        Plot strand coordinates on the solar disk.
 
         See Also
         --------
         synthesizAR.visualize.plot_fieldlines
         """
-        plot_fieldlines(*[_.coordinate for _ in self.loops], **kwargs)
+        plot_fieldlines(*[_.coordinate for _ in self.strands], **kwargs)
 
     def configure_loop_simulations(self, interface, **kwargs):
         """
-        Configure hydrodynamic simulations for each loop object
+        Configure hydrodynamic simulations for each strand object
         """
-        for loop in self.loops:
-            interface.configure_input(loop, **kwargs)
+        for strand in self.strands:
+            interface.configure_input(strand, **kwargs)
 
     @staticmethod
-    def _load_loop_simulation(loop, root=None, interface=None):
+    def _load_loop_simulation(strand, root=None, interface=None):
         # Load in parameters from interface
         (time, electron_temperature, ion_temperature,
-            density, velocity) = interface.load_results(loop)
+            density, velocity) = interface.load_results(strand)
         # If no Zarr file is passed, set the quantites as attributes on the loops
         if root is None:
-            loop._time = time
-            loop._electron_temperature = electron_temperature
-            loop._ion_temperature = ion_temperature
-            loop._density = density
-            loop._velocity = velocity
-            loop._simulation_type = interface.name
+            strand._time = time
+            strand._electron_temperature = electron_temperature
+            strand._ion_temperature = ion_temperature
+            strand._density = density
+            strand._velocity = velocity
+            strand._simulation_type = interface.name
         else:
             # Write to file
-            grp = root.create_group(loop.name)
+            grp = root.create_group(strand.name)
             grp.attrs['simulation_type'] = interface.name
             # time
             dset_time = grp.create_dataset('time', data=time.value)
             dset_time.attrs['unit'] = time.unit.to_string()
             # NOTE: Set the chunk size such that accessing all entries for a given timestep
             # is the most efficient pattern.
-            chunks = (None,) + loop.field_aligned_coordinate_center.shape
+            chunks = (None,) + strand.field_aligned_coordinate_center.shape
             # electron temperature
             dset_electron_temperature = grp.create_dataset(
                 'electron_temperature', data=electron_temperature.value, chunks=chunks)
@@ -244,24 +245,31 @@ Number of loops: {len(self.loops)}'''
 
     def load_loop_simulations(self, interface, filename=None, **kwargs):
         """
-        Load in loop parameters from hydrodynamic results.
+        Load results from hydrodynamic results.
+
+        Parameters
+        ----------
+        interface : model interface object
+            Interface to the hydrodynamic simulation from which to load the results
+        filename : `str` or path-like
+            Path to `zarr` store to write hydrodynamic results to
         """
         if filename is None:
             root = None
         else:
             root = zarr.open(store=filename, mode='w', **kwargs)
-        for loop in self.loops:
-            loop.model_results_filename = filename
+        for strand in self.strands:
+            strand.model_results_filename = filename
         try:
             import distributed
             client = distributed.get_client()
         except (ImportError, ValueError):
-            for l in self.loops:
+            for l in self.strands:
                 self._load_loop_simulation(l, root=root, interface=interface)
         else:
             status = client.map(
                 self._load_loop_simulation,
-                self.loops,
+                self.strands,
                 root=root,
                 interface=interface,
             )
@@ -287,11 +295,11 @@ Number of loops: {len(self.loops)}'''
 
         from synthesizAR.atomic import equilibrium_ionization
 
-        root = zarr.open(store=self.loops[0].model_results_filename, mode='a', **kwargs)
+        root = zarr.open(store=self.strands[0].model_results_filename, mode='a', **kwargs)
         # Check if we can load from the model
         FROM_MODEL = False
         if interface is not None and hasattr(interface, 'load_ionization_fraction'):
-            frac = interface.load_ionization_fraction(self.loops[0], emission_model[0])
+            frac = interface.load_ionization_fraction(self.strands[0], emission_model[0])
             # Some models may optionally output the ionization fractions such that
             # they will have this method, but it may not return anything
             if frac is not None:
@@ -301,17 +309,17 @@ Number of loops: {len(self.loops)}'''
         for el_name in element_names:
             el = Element(el_name, emission_model.temperature)
             ions = [i for i in emission_model if i.element_name == el.element_name]
-            for loop in self.loops:
-                chunks = (None,) + loop.field_aligned_coordinate_center.shape
+            for strand in self.strands:
+                chunks = (None,) + strand.field_aligned_coordinate_center.shape
                 if not FROM_MODEL:
-                    frac_el = equilibrium_ionization(el, loop.electron_temperature)
-                if 'ionization_fraction' in root[loop.name]:
-                    grp = root[f'{loop.name}/ionization_fraction']
+                    frac_el = equilibrium_ionization(el, strand.electron_temperature)
+                if 'ionization_fraction' in root[strand.name]:
+                    grp = root[f'{strand.name}/ionization_fraction']
                 else:
-                    grp = root[loop.name].create_group('ionization_fraction')
+                    grp = root[strand.name].create_group('ionization_fraction')
                 for ion in ions:
                     if FROM_MODEL:
-                        frac = interface.load_ionization_fraction(loop, ion)
+                        frac = interface.load_ionization_fraction(strand, ion)
                         desc = f'{ion.ion_name} ionization fraction computed by {interface.name}'
                     else:
                         frac = frac_el[:, :, ion.charge_state]
diff --git a/synthesizAR/loop.py b/synthesizAR/strand.py
similarity index 95%
rename from synthesizAR/loop.py
rename to synthesizAR/strand.py
index 23b8b2f..f6511ae 100644
--- a/synthesizAR/loop.py
+++ b/synthesizAR/strand.py
@@ -1,5 +1,5 @@
 """
-Loop object for holding field-aligned coordinates and quantities
+Strand object for holding field-aligned coordinates and quantities
 """
 import astropy.units as u
 import numpy as np
@@ -10,10 +10,10 @@ from astropy.coordinates import SkyCoord
 from scipy.interpolate import interp1d, splev, splprep
 from sunpy.coordinates import HeliographicStonyhurst
 
-__all__ = ['Loop']
+__all__ = ['Strand']
 
 
-class Loop:
+class Strand:
     r"""
     Container for geometric and thermodynamic properties of a coronal loop
 
@@ -21,7 +21,7 @@ class Loop:
     ----------
     name : `str`
     coordinate : `astropy.coordinates.SkyCoord`
-        Loop coordinates; should be able to transform to HEEQ
+        Coordinates in the field-aligned direction; should be able to transform to HEEQ
     field_strength : `astropy.units.Quantity`
         Scalar magnetic field strength along the loop. If not specified, defaults
         to NaN with same shape as ``coordinate``.
@@ -39,11 +39,11 @@ class Loop:
     >>> import synthesizAR
     >>> coordinate = SkyCoord(x=[1,4]*u.Mm, y=[2,5]*u.Mm, z=[3,6]*u.Mm, frame='heliographic_stonyhurst', representation_type='cartesian')
     >>> field_strength = u.Quantity([100,200], 'gauss')
-    >>> loop = synthesizAR.Loop('coronal_loop', coordinate, field_strength)
-    >>> loop
-    synthesizAR Loop
-    ----------------
-    Name : coronal_loop
+    >>> strand = synthesizAR.Strand('coronal_strand', coordinate, field_strength)
+    >>> strand
+    synthesizAR Strand
+    ------------------
+    Name : coronal_strand
     Loop full-length, L : 5.196 Mm
     Footpoints : (1 Mm,2 Mm,3 Mm),(4 Mm,5 Mm,6 Mm)
     Maximum field strength : 200.00 G
@@ -76,8 +76,8 @@ class Loop:
     def __repr__(self):
         f0 = f'{self.coordinate.x[0]:.3g},{self.coordinate.y[0]:.3g},{self.coordinate.z[0]:.3g}'
         f1 = f'{self.coordinate.x[-1]:.3g},{self.coordinate.y[-1]:.3g},{self.coordinate.z[-1]:.3g}'
-        return f'''synthesizAR Loop
-----------------
+        return f'''synthesizAR Strand
+------------------
 Name : {self.name}
 Loop full-length, L : {self.length.to(u.Mm):.3f}
 Footpoints : ({f0}),({f1})
@@ -298,7 +298,7 @@ def add_property(name, unit, doc):
         return self._get_quantity(name)
     property_template.__doc__ = doc
     property_template.__name__ = name
-    setattr(Loop, property_template.__name__, property(property_template))
+    setattr(Strand, property_template.__name__, property(property_template))
 
 
 properties = [
diff --git a/synthesizAR/util/util.py b/synthesizAR/util/util.py
index 2f89e4c..4e454eb 100644
--- a/synthesizAR/util/util.py
+++ b/synthesizAR/util/util.py
@@ -2,16 +2,16 @@
 Some basic tools/utilities needed for active region construction. These functions are generally
 peripheral to the actual physics.
 """
-from collections import namedtuple
-import warnings
-
-import numpy as np
-from scipy.interpolate import RegularGridInterpolator
-import astropy.units as u
-from astropy.coordinates import SkyCoord
 import astropy.constants as const
+import astropy.units as u
+import numpy as np
 import sunpy.coordinates
 import sunpy.sun.constants as sun_const
+import warnings
+
+from astropy.coordinates import SkyCoord
+from collections import namedtuple
+from scipy.interpolate import RegularGridInterpolator
 
 import synthesizAR
 
@@ -208,7 +208,7 @@ def from_pfsspy(fieldlines,
                 name_template=None,
                 cross_sectional_area=None):
     """
-    Convert a `pfsspy.fieldline.FieldLines` structure into a list of `~synthesizAR.Loop` objects.
+    Convert a `pfsspy.fieldline.FieldLines` structure into a list of `~synthesizAR.Strand` objects.
 
     Parameters
     ----------
@@ -220,18 +220,18 @@ def from_pfsspy(fieldlines,
         The desired obstime of the coordinates. Use this if the coordinates need to be at a
         different obstime than that of the Carrington map they were traced from.
     length_min : `astropy.units.Quantity`, optional
-        Minimum allowed loop length. All loops with length below this are excluded.
+        Minimum allowed loop length. All strands with length below this are excluded.
     length_max : `astropy.units.Quantity`, optional
-        Maximum allowed loop length. All loops with length above this are excluded.
+        Maximum allowed loop length. All strands with length above this are excluded.
     name_template: `str`, optional
-        Name template to use when building loops. Defaults to 'loop_{:06d}'
+        Name template to use when building strands. Defaults to 'strand_{:06d}'
     cross_sectional_area: `astropy.units.Quantity`, optional
         The cross-sectional area to assign to each loop.
     """
     from synthesizAR import log
     if name_template is None:
-        name_template = 'loop_{:06d}'
-    loops = []
+        name_template = 'strand_{:06d}'
+    strands = []
     for i, f in enumerate(fieldlines):
         if f.coords.shape[0] <= n_min:
             log.debug(f'Dropping {f} as it has less than {n_min} points.')
@@ -258,21 +258,21 @@ def from_pfsspy(fieldlines,
             coord = f.coords
         # Construct the loop here to easily filter on loop length and interpolate
         # NaNs from the field strength.
-        loop = synthesizAR.Loop(name_template.format(i),
+        strand = synthesizAR.Strand(name_template.format(i),
                                 coord,
                                 field_strength=b,
                                 cross_sectional_area=cross_sectional_area)
-        if loop.length < length_min or loop.length > length_max:
-            log.debug(f'Dropping {loop} as it has length {loop.length} outside of the allowed range.')
+        if strand.length < length_min or strand.length > length_max:
+            log.debug(f'Dropping {strand} as it has length {strand.length} outside of the allowed range.')
             continue
-        if np.any(np.isnan(loop.field_strength)):
+        if np.any(np.isnan(strand.field_strength)):
             # There are often NaN values that show up in the interpolated field strengths.
             # Interpolate over these.
-            b = loop.field_strength
-            s = loop.field_aligned_coordinate
+            b = strand.field_strength
+            s = strand.field_aligned_coordinate
             nans = np.isnan(b)
             b[nans] = np.interp(s[nans].value, s[~nans].value, b[~nans].value) * b.unit
-            loop.field_strength = b
-        loops.append(loop)
+            strand.field_strength = b
+        strands.append(strand)
 
-    return loops
+    return strands


Test patch:
diff --git a/synthesizAR/conftest.py b/synthesizAR/conftest.py
index 1973197..9d0853c 100644
--- a/synthesizAR/conftest.py
+++ b/synthesizAR/conftest.py
@@ -1,21 +1,22 @@
 """
 Configure test skeletons here so that they can be used everywhere
 """
+import astropy.units as u
 import pytest
 
-import astropy.units as u
 from sunpy.coordinates import get_earth
 
 import synthesizAR
-from synthesizAR.models import semi_circular_loop, semi_circular_arcade
+
 from synthesizAR.interfaces import MartensInterface
+from synthesizAR.models import semi_circular_arcade, semi_circular_loop
 
 
 @pytest.fixture
 def bare_skeleton():
     observer = get_earth(time='2020-01-01T00:00:00')
     arcade = semi_circular_arcade(100*u.Mm, 20*u.deg, 10, observer)
-    loops = [synthesizAR.Loop(f'{i}', c) for i, c in enumerate(arcade)]
+    loops = [synthesizAR.Strand(f'{i}', c) for i, c in enumerate(arcade)]
     return synthesizAR.Skeleton(loops)
 
 
@@ -29,4 +30,4 @@ def skeleton_with_model(bare_skeleton):
 @pytest.fixture
 def semi_circle_strand(bare_skeleton):
     coords = semi_circular_loop(length=100*u.Mm)
-    return synthesizAR.Loop('test', coords)
+    return synthesizAR.Strand('test', coords)
diff --git a/synthesizAR/tests/test_field.py b/synthesizAR/tests/test_field.py
index 44127ec..6a7ac15 100644
--- a/synthesizAR/tests/test_field.py
+++ b/synthesizAR/tests/test_field.py
@@ -1,28 +1,28 @@
 """
 Tests for Skeleton object
 """
+import astropy.units as u
 import pathlib
-
 import pytest
-import astropy.units as u
+
 from astropy.coordinates import SkyCoord
 
 import synthesizAR
 
 
-def test_skeleton_has_loops(bare_skeleton):
-    assert hasattr(bare_skeleton, 'loops')
-    assert type(bare_skeleton.loops) is list
-    for l in bare_skeleton.loops:
-        assert isinstance(l, synthesizAR.Loop)
+def test_skeleton_has_strands(bare_skeleton):
+    assert hasattr(bare_skeleton, 'strands')
+    assert isinstance(bare_skeleton.strands, list)
+    for l in bare_skeleton.strands:
+        assert isinstance(l, synthesizAR.Strand)
 
 
 def test_create_skeleton_from_coords(bare_skeleton):
-    coords = [l.coordinate for l in bare_skeleton.loops]
-    field_strengths = [l.field_strength for l in bare_skeleton.loops]
+    coords = [l.coordinate for l in bare_skeleton.strands]
+    field_strengths = [l.field_strength for l in bare_skeleton.strands]
     skeleton = synthesizAR.Skeleton.from_coordinates(coords, field_strengths=field_strengths)
-    assert hasattr(skeleton, 'loops')
-    assert type(skeleton.loops) is list
+    assert hasattr(skeleton, 'strands')
+    assert isinstance(skeleton.strands, list)
 
 
 def test_roundtrip(bare_skeleton, tmpdir):
@@ -30,17 +30,17 @@ def test_roundtrip(bare_skeleton, tmpdir):
     filename = pathlib.Path(dirname) / 'test-save.asdf'
     bare_skeleton.to_asdf(filename)
     skeleton_2 = synthesizAR.Skeleton.from_asdf(filename)
-    assert len(bare_skeleton.loops) == len(skeleton_2.loops)
-    for i in range(len(bare_skeleton.loops)):
-        l1 = bare_skeleton.loops[i].coordinate.cartesian.xyz
-        l2 = skeleton_2.loops[i].coordinate.cartesian.xyz
+    assert len(bare_skeleton.strands) == len(skeleton_2.strands)
+    for i in range(len(bare_skeleton.strands)):
+        l1 = bare_skeleton.strands[i].coordinate.cartesian.xyz
+        l2 = skeleton_2.strands[i].coordinate.cartesian.xyz
         assert u.allclose(l2, l1, rtol=1e-9)
 
 
 def test_refine_loops(bare_skeleton):
-    bare_skeleton_refined = bare_skeleton.refine_loops(1*u.Mm)
+    bare_skeleton_refined = bare_skeleton.refine_strands(1*u.Mm)
     assert isinstance(bare_skeleton_refined, synthesizAR.Skeleton)
-    assert len(bare_skeleton_refined.loops) == len(bare_skeleton.loops)
+    assert len(bare_skeleton_refined.strands) == len(bare_skeleton.strands)
 
 
 @pytest.mark.parametrize(
@@ -54,7 +54,7 @@ def test_coordinate_properties(bare_skeleton, name):
 
 
 def test_loops_have_model_type(skeleton_with_model):
-    for l in skeleton_with_model.loops:
+    for l in skeleton_with_model.strands:
         assert hasattr(l, 'simulation_type')
 
 
@@ -69,5 +69,5 @@ def test_loops_have_model_type(skeleton_with_model):
 )
 def test_loops_have_model_quantities(skeleton_with_model, name):
     "These quantities exist only after an interface is defined"
-    for l in skeleton_with_model.loops:
+    for l in skeleton_with_model.strands:
         assert isinstance(getattr(l, name), u.Quantity)
diff --git a/synthesizAR/tests/test_loop.py b/synthesizAR/tests/test_loop.py
index 21c4ce9..a90d2e3 100644
--- a/synthesizAR/tests/test_loop.py
+++ b/synthesizAR/tests/test_loop.py
@@ -1,13 +1,15 @@
 """
-Tests for Loop object
+Tests for Strand object
 """
-import pytest
-import numpy as np
 import astropy.units as u
+import numpy as np
+import pytest
+
 from astropy.coordinates import SkyCoord
 from sunpy.coordinates import get_earth
 
 import synthesizAR
+
 from synthesizAR.models import semi_circular_loop
 
 
@@ -58,7 +60,7 @@ def simple_strand(field_aligned_coordinate):
     observer = get_earth('2020-01-01')
     coord = semi_circular_loop(s=field_aligned_coordinate,
                                observer=observer)
-    return synthesizAR.Loop('test', coord)
+    return synthesizAR.Strand('test', coord)
 
 
 def test_loop_field_aligned_coordinate(simple_strand,