diff --git a/.gitignore b/.gitignore
index fd6a896..7c5f19b 100644
--- a/.gitignore
+++ b/.gitignore
@@ -265,3 +265,4 @@ package.json
 # Repo specific
 examples/using_sunkit_magex_pfss/aia_map.fits
 docs/sg_execution_times.rst
+aia_map.fits
diff --git a/CHANGELOG.rst b/CHANGELOG.rst
new file mode 100644
index 0000000..5c6435e
--- /dev/null
+++ b/CHANGELOG.rst
@@ -0,0 +1,22 @@
+1.0.0 (2024-05-31)
+==================
+
+This is the first release and aims to keep the API the same from `pfsspy` to
+`sunkit_magex`.  The main difference is that you will need to replace the
+import like so:
+
+  .. code-block:: python
+
+    # Before
+    import pfsspy
+    from pfsspy import tracing
+
+    # After
+    from sunkit_magex import pfss as pfsspy
+    from sunkit_magex.pfss import tracing
+
+The main changes from the previous release of `pfsspy` are as follows:
+
+* The ``GongSynopticMap`` class has moved into `sunpy`, note that the new
+  class in ``sunpy`` is slightly different in that it does not modify the
+  header.
diff --git a/changelog/63.breaking.1.rst b/changelog/63.breaking.1.rst
new file mode 100644
index 0000000..d5013c1
--- /dev/null
+++ b/changelog/63.breaking.1.rst
@@ -0,0 +1 @@
+Made streamtracer a hard dependency of sunkit-magex.
diff --git a/changelog/63.breaking.2.rst b/changelog/63.breaking.2.rst
new file mode 100644
index 0000000..f229281
--- /dev/null
+++ b/changelog/63.breaking.2.rst
@@ -0,0 +1,6 @@
+Deprecated the following pfss utils:
+
+- pfss.utils.fix_hmi_meta
+- pfss.utils.load_adapt
+
+These functions are no longer needed and will be removed.
diff --git a/changelog/63.breaking.rst b/changelog/63.breaking.rst
new file mode 100644
index 0000000..28b0e40
--- /dev/null
+++ b/changelog/63.breaking.rst
@@ -0,0 +1 @@
+Increased the minimum required version of ``sunpy``  to v6.0.0.
diff --git a/changelog/README.rst b/changelog/README.rst
new file mode 100644
index 0000000..7d388e3
--- /dev/null
+++ b/changelog/README.rst
@@ -0,0 +1,34 @@
+=========
+Changelog
+=========
+
+.. note::
+
+    This README was adapted from the pytest changelog readme under the terms of the MIT licence.
+
+This directory contains "news fragments" which are short files that contain a small **ReST**-formatted text that will be added to the next ``CHANGELOG``.
+
+The ``CHANGELOG`` will be read by users, so this description should be aimed at SunPy users instead of describing internal changes which are only relevant to the developers.
+
+Make sure to use full sentences with correct case and punctuation, for example::
+
+    Add support for Helioprojective coordinates in `sunpy.coordinates.frames`.
+
+Please try to use Sphinx intersphinx using backticks.
+
+Each file should be named like ``<PULL REQUEST>.<TYPE>[.<COUNTER>].rst``, where ``<PULL REQUEST>`` is a pull request number, ``COUNTER`` is an optional number if a PR needs multiple entries with the same type and ``<TYPE>`` is one of:
+
+* ``breaking``: A change which requires users to change code and is not backwards compatible. (Not to be used for removal of deprecated features.)
+* ``feature``: New user facing features and any new behavior.
+* ``bugfix``: Fixes a reported bug.
+* ``doc``: Documentation addition or improvement, like rewording an entire session or adding missing docs.
+* ``deprecation``: Feature deprecation
+* ``removal``: Feature removal.
+* ``trivial``: A change which has no user facing effect or is tiny change.
+
+So for example: ``123.feature.rst``, ``456.bugfix.rst``.
+
+If you are unsure what pull request type to use, don't hesitate to ask in your PR.
+
+Note that the ``towncrier`` tool will automatically reflow your text, so it will work best if you stick to a single paragraph, but multiple sentences and links are OK and encouraged.
+You can install ``towncrier`` and then run ``towncrier --draft`` if you want to get a preview of how your change will look in the final release notes.
diff --git a/docs/changelog.rst b/docs/changelog.rst
index c22f33b..c9c3609 100644
--- a/docs/changelog.rst
+++ b/docs/changelog.rst
@@ -4,25 +4,7 @@
 Full Changelog
 **************
 
-1.0.0 (2024-05-31)
-==================
-
-This is the first release and aims to keep the API the same from `pfsspy` to
-`sunkit_magex`.  The main difference is that you will need to replace the
-import like so:
-
-  .. code-block:: python
-
-    # Before
-    import pfsspy
-    from pfsspy import tracing
-
-    # After
-    from sunkit_magex import pfss as pfsspy
-    from sunkit_magex.pfss import tracing
-
-The main changes from the previous release of `pfsspy` are as follows:
-
-* The ``GongSynopticMap`` class has moved into `sunpy`, note that the new
-  class in ``sunpy`` is slightly different in that it does not modify the
-  header.
+.. changelog::
+   :towncrier: ../
+   :towncrier-skip-if-empty:
+   :changelog_file: ../CHANGELOG.rst
diff --git a/docs/conf.py b/docs/conf.py
index a82bbe1..9cae120 100644
--- a/docs/conf.py
+++ b/docs/conf.py
@@ -35,6 +35,7 @@
     "sphinx_automodapi.automodapi",
     "sphinx_automodapi.smart_resolver",
     'sphinx_gallery.gen_gallery',
+    'sphinx_changelog',
 ]
 
 # Add any paths that contain templates here, relative to this directory.
diff --git a/docs/installing.rst b/docs/installing.rst
index 9958e7f..970f6e0 100644
--- a/docs/installing.rst
+++ b/docs/installing.rst
@@ -17,8 +17,8 @@ Alternatively ``sunkit-magex`` can be installed from PyPI using:
     pip install sunkit-magex
 
 This will install sunkit_magex and all of its dependencies.
-In addition to the core dependencies, there are optional dependencies (numba, streamtracer) that can improve code performance.
-These can be installed with
+In addition to the core dependencies, numba is an optional dependency that can improve code performance.
+This can be installed with:
 
 .. code:: shell
 
diff --git a/docs/performance.rst b/docs/performance.rst
index b7d282a..3c90dbe 100644
--- a/docs/performance.rst
+++ b/docs/performance.rst
@@ -11,9 +11,8 @@ To enable this simply `install numba`_  and use `sunkit_magex.pfss` as normal.
 Streamline tracing
 ==================
 
-`sunkit_magex.pfss` has two streamline tracers: a pure python `sunkit_magex.pfss.tracing.PythonTracer` and a complied tracer `sunkit_magex.pfss.tracing.FortranTracer`.
-The compiled version is significantly faster, using the `streamtracer`_ package.
+`sunkit_magex.pfss` uses a complied tracer `sunkit_magex.pfss.tracing.PerformanceTracer` using the `streamtracer`_ package.
 
 .. _numba: https://numba.pydata.org
 .. _install numba: http://numba.pydata.org/numba-doc/latest/user/installing.html
-.. _streamtracer: https://github.com/sunpy/streamtracer
+.. _streamtracer: https://docs.sunpy.org/projects/streamtracer/en/stable/
diff --git a/examples/finding_data/README.rst b/examples/finding_data/README.rst
deleted file mode 100644
index 0cfdfce..0000000
--- a/examples/finding_data/README.rst
+++ /dev/null
@@ -1,3 +0,0 @@
-Finding data
-------------
-Examples showing how to find, download, and load magnetograms.
diff --git a/examples/finding_data/adapt_map_sequence.py b/examples/finding_data/adapt_map_sequence.py
deleted file mode 100644
index 0de7ed4..0000000
--- a/examples/finding_data/adapt_map_sequence.py
+++ /dev/null
@@ -1,65 +0,0 @@
-"""
-Parsing ADAPT Ensemble .fits files
-==================================
-
-Parse an ADAPT FITS file into a `sunpy.map.MapSequence`.
-"""
-import matplotlib.pyplot as plt
-from matplotlib import gridspec
-
-from astropy.io import fits
-
-import sunpy.map
-
-from sunkit_magex import pfss
-
-###############################################################################
-# Load an example ADAPT fits file.
-
-adapt_fname = pfss.sample_data.get_adapt_map()
-
-###############################################################################
-# ADAPT synoptic magnetograms contain 12 realizations of synoptic magnetograms
-# output as a result of varying model assumptions. See
-# [here](https://www.swpc.noaa.gov/sites/default/files/images/u33/SWW_2012_Talk_04_27_2012_Arge.pdf))
-#
-# Because the fits data is 3D, it cannot be passed directly to `sunpy.map.Map`,
-# because this will take the first slice only and the other realizations are
-# lost. We want to end up with a `sunpy.map.MapSequence` containing all these
-# realiations as individual maps. These maps can then be individually accessed
-# and PFSS solutions generated from them.
-#
-# We first read in the fits file:
-
-adapt_fits = fits.open(adapt_fname)
-
-###############################################################################
-# ``adapt_fits`` is a list of ``HDPair`` objects. The first of these contains
-# the 12 realizations data and a header with sufficient information to build
-# the `~sunpy.map.MapSequence`. We unpack this ``HDPair`` into a list of
-# ``(data,header)`` tuples where ``data`` are the different adapt realizations.
-
-data_header_pairs = [(map_slice, adapt_fits[0].header) for map_slice in adapt_fits[0].data]
-
-###############################################################################
-# Next, pass this list of tuples as the argument to `sunpy.map.Map` to create
-# the map sequence:
-
-adapt_maps = sunpy.map.Map(data_header_pairs, sequence=True)
-
-###############################################################################
-# ``adapt_map_sequence`` is now a list of our individual adapt realizations.
-# Note the ``.peek()` and ``.plot()`` methods of `~sunpy.map.MapSequence`
-# returns instances of
-# ``sunpy.visualization.MapSequenceAnimator`` and
-# ``matplotlib.animation.FuncAnimation1``. Here, we generate a static
-# plot accessing the individual maps in turn:
-
-fig = plt.figure(figsize=(7, 8))
-gs = gridspec.GridSpec(4, 3, figure=fig)
-for i, a_map in enumerate(adapt_maps):
-    ax = fig.add_subplot(gs[i], projection=a_map)
-    a_map.plot(axes=ax, cmap='bwr', vmin=-2, vmax=2, title=f"Realization {1+i:02d}")
-plt.tight_layout(pad=5, h_pad=2)
-
-plt.show()
diff --git a/examples/internals/plot_computation_grid.py b/examples/internals/computation_grid.py
similarity index 98%
rename from examples/internals/plot_computation_grid.py
rename to examples/internals/computation_grid.py
index 1d7020e..cf0730c 100644
--- a/examples/internals/plot_computation_grid.py
+++ b/examples/internals/computation_grid.py
@@ -1,5 +1,5 @@
 r"""
-Magnetic field grid
+Magnetic Field Grid
 ===================
 
 A plot of the grid corners which the magnetic field values are taken from
diff --git a/examples/internals/tracer_performance.py b/examples/internals/tracer_performance.py
index 057b3b0..7809f87 100644
--- a/examples/internals/tracer_performance.py
+++ b/examples/internals/tracer_performance.py
@@ -1,8 +1,8 @@
 """
-Tracer performance
+Tracer Performance
 ==================
 
-Comparing the performance of Python and Rust tracers.
+Comparing the performance of the Python and Compiled tracers.
 """
 import timeit
 
@@ -19,47 +19,45 @@
 ###############################################################################
 # Create a dipole map.
 
+def dipole_Br(r, theta):
+    return 2 * np.sin(theta) / r**3
+
+
 ntheta = 180
 nphi = 360
 nr = 50
 rss = 2.5
-
 phi = np.linspace(0, 2 * np.pi, nphi)
 theta = np.linspace(-np.pi / 2, np.pi / 2, ntheta)
 theta, phi = np.meshgrid(theta, phi)
-
-def dipole_Br(r, theta):
-    return 2 * np.sin(theta) / r**3
-
-
 br = dipole_Br(1, theta)
 br = sunpy.map.Map(br.T, pfss.utils.carr_cea_wcs_header('2010-01-01', br.shape))
 pfss_input = pfss.Input(br, nr, rss)
 pfss_output = pfss.pfss(pfss_input)
-print('Computed PFSS solution')
 
 ###############################################################################
 # Trace some field lines.
 
 seed0 = np.atleast_2d(np.array([1, 1, 0]))
-tracers = [pfss.tracing.PythonTracer(),
-           pfss.tracing.FortranTracer()]
+tracers = [
+    pfss.tracing.PythonTracer(),
+    pfss.tracing.PerformanceTracer(),
+]
 nseeds = 2**np.arange(14)
-times = [[], []]
+times = [[] for _ in tracers]
 
 for nseed in nseeds:
-    print(nseed)
     seeds = np.repeat(seed0, nseed, axis=0)
     r, lat, lon = pfss.coords.cart2sph(seeds[:, 0], seeds[:, 1], seeds[:, 2])
     r = r * astropy.constants.R_sun
     lat = (lat - np.pi / 2) * u.rad
     lon = lon * u.rad
     seeds = astropy.coordinates.SkyCoord(lon, lat, r, frame=pfss_output.coordinate_frame)
-
     for i, tracer in enumerate(tracers):
+        # Skip the Python tracer for large numbers of seeds.
+        # It is too slow.
         if nseed > 64 and i == 0:
             continue
-
         t = timeit.timeit(lambda: tracer.trace(seeds, pfss_output), number=1)
         times[i].append(t)
 
@@ -67,24 +65,23 @@ def dipole_Br(r, theta):
 # Plot the results.
 
 fig, ax = plt.subplots()
-ax.scatter(nseeds[1:len(times[0])], times[0][1:], label='python')
-ax.scatter(nseeds[1:], times[1][1:], label='compiled')
 
-pydt = (times[0][4] - times[0][3]) / (nseeds[4] - nseeds[3])
-ax.plot([1, 1e5], [pydt, 1e5 * pydt])
-
-fort0 = times[1][1]
-fordt = (times[1][-1] - times[1][-2]) / (nseeds[-1] - nseeds[-2])
-ax.plot(np.logspace(0, 5, 100), fort0 + fordt * np.logspace(0, 5, 100))
+for i, tracer in enumerate(tracers):
+    if i == 0:
+        ax.scatter(nseeds[1:len(times[i])], times[i][1:], label=tracer.__class__.__name__)
+        pydt = (times[0][4] - times[0][3]) / (nseeds[4] - nseeds[3])
+        ax.plot([1, 1e5], [pydt, 1e5 * pydt])
+    else:
+        ax.scatter(nseeds[1:], times[i][1:], label=tracer.__class__.__name__)
+        t0 = times[i][1]
+        dt = (times[i][-1] - times[i][-2]) / (nseeds[-1] - nseeds[-2])
+        ax.plot(np.logspace(0, 5, 100), t0 + dt * np.logspace(0, 5, 100))
 
 ax.set_xscale('log')
 ax.set_yscale('log')
-
 ax.set_xlabel('Number of seeds')
 ax.set_ylabel('Seconds')
-
 ax.axvline(180 * 360, color='k', linestyle='--', label='180x360 seed points')
-
 ax.legend()
 
 plt.show()
diff --git a/examples/testing/helpers.py b/examples/testing/_helpers.py
similarity index 99%
rename from examples/testing/helpers.py
rename to examples/testing/_helpers.py
index 0e5c726..d8b048f 100644
--- a/examples/testing/helpers.py
+++ b/examples/testing/_helpers.py
@@ -16,7 +16,6 @@
 import sunkit_magex.pfss.utils
 
 result_dir = (pathlib.Path(__file__) / '..' / 'results').resolve()
-
 pi = np.pi * u.rad
 
 
@@ -83,8 +82,6 @@ def open_flux_numeric(l, m, zss, nrho):
     return np.sum(np.abs(br)) * (4 * np.pi) / nphi / ns
 
 
-######################
-# Field line helpers #
 @u.quantity_input
 def fr(r: u.m, rss: u.m, l):
     rho = r / rss
diff --git a/examples/testing/plot_error_map.py b/examples/testing/field_line_error_map.py
similarity index 90%
rename from examples/testing/plot_error_map.py
rename to examples/testing/field_line_error_map.py
index 06d8531..c68ead7 100644
--- a/examples/testing/plot_error_map.py
+++ b/examples/testing/field_line_error_map.py
@@ -7,13 +7,13 @@
 """
 import matplotlib.pyplot as plt
 import numpy as np
+from _helpers import pffspy_output, phi_fline_coords, theta_fline_coords
 
 import astropy.constants as const
 import astropy.units as u
 from astropy.coordinates import SkyCoord
 from astropy.visualization import quantity_support
 
-from examples.testing.helpers import pffspy_output, phi_fline_coords, theta_fline_coords
 from sunkit_magex import pfss
 
 quantity_support()
@@ -48,7 +48,7 @@
 dthetas = []
 print(f'Tracing {step_size}...')
 # Trace
-tracer = pfss.tracing.FortranTracer(step_size=step_size)
+tracer = pfss.tracing.PerformanceTracer(step_size=step_size)
 flines = tracer.trace(seeds, pfss_out)
 # Set a mask of open field lines
 mask = flines.connectivities.astype(bool).reshape(theta.shape)
@@ -73,7 +73,6 @@
 
 fig, axs = plt.subplots(nrows=2, sharex=True, sharey=True)
 
-
 def plot_map(field, ax, label, title):
     kwargs = dict(cmap='RdBu', vmin=-0.5, vmax=0.5, shading='nearest', edgecolors='face')
     im = ax.pcolormesh(phi.to_value(u.deg), np.sin(theta).value,
@@ -85,10 +84,10 @@ def plot_map(field, ax, label, title):
 
 
 plot_map(dtheta.to_value(u.deg), axs[0],
-         r'$\theta_{sunkit_magex.pfss} - \theta_{analytic}$ (deg)',
+         r'$\theta_{pfss} - \theta_{analytic}$ (deg)',
          'Error in latitude')
 plot_map(dphi.to_value(u.deg), axs[1],
-         r'$\phi_{sunkit_magex.pfss} - \phi_{analytic}$ (deg)',
+         r'$\phi_{pfss} - \phi_{analytic}$ (deg)',
          'Error in longitude')
 
 ax = axs[1]
diff --git a/examples/testing/plot_open_flux_harmonics.py b/examples/testing/open_flux_harmonics_comparison.py
similarity index 74%
rename from examples/testing/plot_open_flux_harmonics.py
rename to examples/testing/open_flux_harmonics_comparison.py
index 9af13ec..7c739e6 100644
--- a/examples/testing/plot_open_flux_harmonics.py
+++ b/examples/testing/open_flux_harmonics_comparison.py
@@ -14,8 +14,7 @@
 import matplotlib.cm as cm
 import matplotlib.colors as mcolor
 import matplotlib.pyplot as plt
-
-from examples.testing.helpers import LMAxes, result_dir
+from _helpers import LMAxes, result_dir
 
 with open(result_dir / "open_flux_harmonics.json") as f:
     results = json.load(f, parse_int=int)
@@ -27,15 +26,11 @@
 for lstr in results['analytic']:
     for mstr in results['analytic'][lstr]:
         l, m = int(lstr), int(mstr)
-
         ax = axs[l, m]
-        ax.set_facecolor(cmap(norm(results['numeric'][lstr][mstr] /
-                                   results['analytic'][lstr][mstr])))
+        ax.set_facecolor(cmap(norm(results['numeric'][lstr][mstr] / results['analytic'][lstr][mstr])))
         ax.set_aspect('equal')
 
-cbar = axs.fig.colorbar(cm.ScalarMappable(norm=norm, cmap=cmap),
-                        ax=axs.all_axs)
-cbar.ax.set_ylabel(r'$\frac{\Phi_{sunkit_magex.pfss}}{\Phi_{analytic}}$', rotation=0,
-                   size=18, labelpad=27, va='center')
+cbar = axs.fig.colorbar(cm.ScalarMappable(norm=norm, cmap=cmap), ax=axs.all_axs)
+cbar.ax.set_ylabel(r'$\frac{\Phi_{pfss}}{\Phi_{analytic}}$', rotation=0, size=18, labelpad=27, va='center')
 
 plt.show()
diff --git a/examples/testing/plot_open_flux_nr.py b/examples/testing/open_flux_nr_grid_points.py
similarity index 92%
rename from examples/testing/plot_open_flux_nr.py
rename to examples/testing/open_flux_nr_grid_points.py
index b5002c3..0e2c8d6 100644
--- a/examples/testing/plot_open_flux_nr.py
+++ b/examples/testing/open_flux_nr_grid_points.py
@@ -10,8 +10,7 @@
 import matplotlib.pyplot as plt
 import matplotlib.ticker as mticker
 import pandas as pd
-
-from examples.testing.helpers import LMAxes, result_dir
+from _helpers import LMAxes, result_dir
 
 df = pd.read_csv(result_dir / 'open_flux_results.csv', index_col=0)
 axs = LMAxes(nl=5)
@@ -41,7 +40,7 @@
 
         ax.yaxis.set_ticks([1, 1.05, 1.1])
         ax.yaxis.tick_right()
-        ax.set_ylabel(r'$\frac{\Phi_{sunkit_magex.pfss}}{\Phi_{analytic}}$',
+        ax.set_ylabel(r'$\frac{\Phi_{pfss}}{\Phi_{analytic}}$',
                       rotation=0, labelpad=30, fontsize=16, loc='center')
         ax.yaxis.set_label_position('right')
         ax.yaxis.set_major_formatter(mticker.ScalarFormatter())
diff --git a/examples/testing/skip_open_flux_harmonics.py b/examples/testing/skip_open_flux_harmonics.py
index 5be1edc..028d64f 100644
--- a/examples/testing/skip_open_flux_harmonics.py
+++ b/examples/testing/skip_open_flux_harmonics.py
@@ -12,7 +12,7 @@
 import json
 from collections import defaultdict
 
-from examples.testing.helpers import open_flux_analytic, open_flux_numeric, result_dir
+from _helpers import open_flux_analytic, open_flux_numeric, result_dir
 
 ###############################################################################
 # Set the source surface height, and the (l, m) values to investigate:
@@ -20,8 +20,7 @@
 zss = 2
 nrho = 40
 
-results = {'numeric': defaultdict(dict),
-           'analytic': defaultdict(dict)}
+results = {'numeric': defaultdict(dict), 'analytic': defaultdict(dict)}
 
 for l in range(1, 6):
     for m in range(-l, l + 1):
diff --git a/examples/testing/skip_open_flux_nr.py b/examples/testing/skip_open_flux_nr.py
index 4cd4ed8..f217293 100644
--- a/examples/testing/skip_open_flux_nr.py
+++ b/examples/testing/skip_open_flux_nr.py
@@ -8,14 +8,11 @@
 fluxes in PFSS solutions of spherical harmonics, as a function of the number of
 radial grid cells in the sunkit_magex.pfss grid.
 """
-
-
 import numpy as np
 import pandas as pd
+from _helpers import open_flux_analytic, open_flux_numeric, result_dir
 from tqdm import tqdm
 
-from examples.testing.helpers import open_flux_analytic, open_flux_numeric, result_dir
-
 ###############################################################################
 # Set the source surface height and range of radial grid points.
 
diff --git a/examples/testing/tracer_step_size.py b/examples/testing/skip_tracer_step_size.py
similarity index 93%
rename from examples/testing/tracer_step_size.py
rename to examples/testing/skip_tracer_step_size.py
index b641b3c..d05383d 100644
--- a/examples/testing/tracer_step_size.py
+++ b/examples/testing/skip_tracer_step_size.py
@@ -4,12 +4,12 @@
 """
 import numpy as np
 import pandas as pd
+from _helpers import pffspy_output, phi_fline_coords, result_dir, theta_fline_coords
 
 import astropy.constants as const
 import astropy.units as u
 from astropy.coordinates import SkyCoord
 
-from examples.testing.helpers import pffspy_output, phi_fline_coords, result_dir, theta_fline_coords
 from sunkit_magex.pfss import tracing
 
 l = 3
@@ -19,14 +19,15 @@
 nr = 40
 rss = 2
 
-
 ###############################################################################
 # Calculate PFSS solution
+
 pfss_out = pffspy_output(nphi, ns, nr, rss, l, m)
 
 ###############################################################################
 # Trace an array of field lines from the source surface down to the solar
 # surface
+
 n = 90
 # Create 1D theta, phi arrays
 phi = np.linspace(0, 360, n * 2)
@@ -45,11 +46,10 @@
 for step_size in step_sizes:
     print(f'Tracing {step_size}...')
     # Trace
-    tracer = tracing.FortranTracer(step_size=step_size)
+    tracer = tracing.PerformanceTracer(step_size=step_size)
     flines = tracer.trace(seeds, pfss_out)
     # Set a mask of open field lines
     mask = flines.connectivities.astype(bool).reshape(theta.shape)
-
     # Get solar surface latitude
     phi_solar = np.ones_like(phi) * np.nan
     phi_solar[mask] = flines.open_field_lines.solar_feet.lon
@@ -57,7 +57,6 @@
     theta_solar[mask] = flines.open_field_lines.solar_feet.lat
     r_out = np.ones_like(theta.value) * const.R_sun * np.nan
     r_out[mask] = flines.open_field_lines.solar_feet.radius
-
     # Calculate analytical solution
     theta_analytic = theta_fline_coords(r_out, rss, l, m, theta)
     dtheta = (theta_solar - theta_analytic).to_value(u.deg)
@@ -66,14 +65,14 @@
     # Wrap phi values
     dphi[dphi > 180] -= 360
     dphi[dphi < -180] += 360
-
     dthetas.append(dtheta.ravel())
     dphis.append(dphi.ravel())
 
 
 ###############################################################################
 # Save results. This saves the maximum error in both phi and theta as a
-# function of thet tracer step size.
+# function of the tracer step size.
+
 dthetas = pd.DataFrame(data=np.array(dthetas), index=step_sizes)
 dthetas = dthetas.mask(np.abs(dthetas) > 30).max(axis=1)
 
diff --git a/examples/testing/plot_harmonic_comparison.py b/examples/testing/spherical_harmonic_comparison.py
similarity index 94%
rename from examples/testing/plot_harmonic_comparison.py
rename to examples/testing/spherical_harmonic_comparison.py
index 19b3864..7074208 100644
--- a/examples/testing/plot_harmonic_comparison.py
+++ b/examples/testing/spherical_harmonic_comparison.py
@@ -6,8 +6,7 @@
 """
 import matplotlib.pyplot as plt
 import matplotlib.ticker as mticker
-
-from examples.testing.helpers import LMAxes, brss_analytic, brss_pfss
+from _helpers import LMAxes, brss_analytic, brss_pfss
 
 ###############################################################################
 # Compare the the `sunkit_magex.pfss` solution to the analytic solutions.
diff --git a/examples/testing/plot_tracer_step_size.py b/examples/testing/tracer_step_size_comparision.py
similarity index 89%
rename from examples/testing/plot_tracer_step_size.py
rename to examples/testing/tracer_step_size_comparision.py
index ca702f7..174f3ba 100644
--- a/examples/testing/plot_tracer_step_size.py
+++ b/examples/testing/tracer_step_size_comparision.py
@@ -5,8 +5,7 @@
 import matplotlib.pyplot as plt
 import matplotlib.ticker as mticker
 import pandas as pd
-
-from examples.testing.helpers import LMAxes, result_dir
+from _helpers import LMAxes, result_dir
 
 nl = 3
 
@@ -16,10 +15,8 @@
     for m in range(-l, l+1):
         ax = axs[l, m]
         try:
-            dphis = pd.read_csv(result_dir / f'flines/dphis_{l}{m}.csv',
-                                header=None, index_col=0)
-            dthetas = pd.read_csv(result_dir / f'flines/dthetas_{l}{m}.csv',
-                                  header=None, index_col=0)
+            dphis = pd.read_csv(result_dir / f'flines/dphis_{l}{m}.csv', header=None, index_col=0)
+            dthetas = pd.read_csv(result_dir / f'flines/dthetas_{l}{m}.csv', header=None, index_col=0)
             print(l, m)
         except FileNotFoundError:
             print(f'❌ Files not found for l={l}, m={m}')
diff --git a/examples/using_sunkit_magex_pfss/plot_dipole.py b/examples/using_sunkit_magex_pfss/dipole_source_solution.py
similarity index 94%
rename from examples/using_sunkit_magex_pfss/plot_dipole.py
rename to examples/using_sunkit_magex_pfss/dipole_source_solution.py
index 5714963..3aa78f0 100644
--- a/examples/using_sunkit_magex_pfss/plot_dipole.py
+++ b/examples/using_sunkit_magex_pfss/dipole_source_solution.py
@@ -5,6 +5,8 @@
 A simple example showing how to use PFSS to compute the solution to a dipole
 source field.
 """
+# sphinx_gallery_thumbnail_number = 3
+
 import matplotlib.patches as mpatch
 import matplotlib.pyplot as plt
 import numpy as np
@@ -26,7 +28,6 @@
 
 nphi = 360
 ns = 180
-
 phi = np.linspace(0, 2 * np.pi, nphi)
 s = np.linspace(-1, 1, ns)
 s, phi = np.meshgrid(s, phi)
@@ -59,15 +60,18 @@ def dipole_Br(r, s):
 ###############################################################################
 # Using the Input object, plot the input field.
 
-m = pfss_in.map
+in_map = pfss_in.map
+
 fig = plt.figure()
-ax = plt.subplot(projection=m)
-m.plot()
+ax = fig.add_subplot(projection=in_map)
+
+in_map.plot(axes=ax)
 plt.colorbar()
 ax.set_title('Input dipole field')
 
 ###############################################################################
 # Now calculate the PFSS solution.
+
 pfss_out = pfss.pfss(pfss_in)
 
 ###############################################################################
@@ -76,15 +80,15 @@ def dipole_Br(r, s):
 
 ss_br = pfss_out.source_surface_br
 
-# Create the figure and axes
 fig = plt.figure()
-ax = plt.subplot(projection=ss_br)
+ax = fig.add_subplot(projection=ss_br)
 
 # Plot the source surface map
-ss_br.plot()
+ss_br.plot(axes=ax)
+
 # Plot the polarity inversion line
 ax.plot_coord(pfss_out.source_surface_pils[0])
-# Plot formatting
+
 plt.colorbar()
 ax.set_title('Source surface magnetic field')
 
@@ -102,7 +106,7 @@ def dipole_Br(r, s):
 lat = np.linspace(-np.pi / 2, np.pi / 2, 33) * u.rad
 seeds = SkyCoord(lon, lat, r, frame=pfss_out.coordinate_frame)
 
-tracer = pfss.tracing.FortranTracer()
+tracer = pfss.tracing.PerformanceTracer()
 field_lines = tracer.trace(seeds, pfss_out)
 
 for field_line in field_lines:
diff --git a/examples/using_sunkit_magex_pfss/plot_gong.py b/examples/using_sunkit_magex_pfss/gong_pfss.py
similarity index 88%
rename from examples/using_sunkit_magex_pfss/plot_gong.py
rename to examples/using_sunkit_magex_pfss/gong_pfss.py
index 3a0ef1f..c8929e6 100644
--- a/examples/using_sunkit_magex_pfss/plot_gong.py
+++ b/examples/using_sunkit_magex_pfss/gong_pfss.py
@@ -37,22 +37,15 @@
 
 pfss_in = pfss.Input(gong_map, nrho, rss)
 
-
 ###############################################################################
 # Using the Input object, plot the input field.
 
-def set_axes_lims(ax):
-    ax.set_xlim(0, 360)
-    ax.set_ylim(0, 180)
-
-
-m = pfss_in.map
+input_map = pfss_in.map
 fig = plt.figure()
-ax = plt.subplot(projection=m)
-m.plot()
+ax = fig.add_subplot(projection=input_map)
+input_map.plot(axes=ax)
 plt.colorbar()
 ax.set_title('Input field')
-set_axes_lims(ax)
 
 ###############################################################################
 # Now calculate the PFSS solution.
@@ -64,18 +57,15 @@ def set_axes_lims(ax):
 # polarity inversion line.
 
 ss_br = pfss_out.source_surface_br
-# Create the figure and axes
+
 fig = plt.figure()
-ax = plt.subplot(projection=ss_br)
+ax = fig.add_subplot(projection=ss_br)
 
-# Plot the source surface map
-ss_br.plot()
+ss_br.plot(axes=ax)
 # Plot the polarity inversion line
 ax.plot_coord(pfss_out.source_surface_pils[0])
-# Plot formatting
 plt.colorbar()
 ax.set_title('Source surface magnetic field')
-set_axes_lims(ax)
 
 ###############################################################################
 # It is also easy to plot the magnetic field at an arbitrary height within
@@ -89,16 +79,12 @@ def set_axes_lims(ax):
 # Get the radial coordinate
 r = np.exp(pfss_out.grid.rc[ridx])
 
-# Create the figure and axes
 fig = plt.figure()
-ax = plt.subplot(projection=br)
+ax = fig.add_subplot(projection=br)
 
-# Plot the source surface map
 br.plot(cmap='RdBu')
-# Plot formatting
 plt.colorbar()
 ax.set_title('$B_{r}$ ' + f'at r={r:.2f}' + '$r_{\\odot}$')
-set_axes_lims(ax)
 
 ###############################################################################
 # Finally, using the 3D magnetic field solution we can trace some field lines.
@@ -108,7 +94,7 @@ def set_axes_lims(ax):
 fig = plt.figure()
 ax = fig.add_subplot(111, projection='3d')
 
-tracer = pfss.tracing.FortranTracer()
+tracer = pfss.tracing.PerformanceTracer()
 r = 1.2 * const.R_sun
 lat = np.linspace(-np.pi / 2, np.pi / 2, 8, endpoint=False)
 lon = np.linspace(0, 2 * np.pi, 8, endpoint=False)
@@ -116,9 +102,7 @@ def set_axes_lims(ax):
 lat, lon = lat.ravel() * u.rad, lon.ravel() * u.rad
 
 seeds = SkyCoord(lon, lat, r, frame=pfss_out.coordinate_frame)
-
 field_lines = tracer.trace(seeds, pfss_out)
-
 for field_line in field_lines:
     color = {0: 'black', -1: 'tab:blue', 1: 'tab:red'}.get(field_line.polarity)
     coords = field_line.coords
@@ -127,7 +111,6 @@ def set_axes_lims(ax):
             coords.y / const.R_sun,
             coords.z / const.R_sun,
             color=color, linewidth=1)
-
 ax.set_title('PFSS solution')
 
 plt.show()
diff --git a/examples/using_sunkit_magex_pfss/plot_hmi.py b/examples/using_sunkit_magex_pfss/hmi_pfss.py
similarity index 93%
rename from examples/using_sunkit_magex_pfss/plot_hmi.py
rename to examples/using_sunkit_magex_pfss/hmi_pfss.py
index da3990f..3467cd8 100644
--- a/examples/using_sunkit_magex_pfss/plot_hmi.py
+++ b/examples/using_sunkit_magex_pfss/hmi_pfss.py
@@ -46,12 +46,12 @@
 # object.
 
 hmi_map = sunpy.map.Map(files[0])
-print('Data shape: ', hmi_map.data.shape)
 
 ###############################################################################
 # Since this map is far to big to calculate a PFSS solution quickly, lets
 # resample it down to a smaller size.
 
+print('Old shape: ', hmi_map.data.shape)
 hmi_map = hmi_map.resample([360, 180] * u.pix)
 print('New shape: ', hmi_map.data.shape)
 
@@ -68,15 +68,13 @@
 # polarity inversion line.
 
 ss_br = pfss_out.source_surface_br
-# Create the figure and axes
+
 fig = plt.figure()
-ax = plt.subplot(projection=ss_br)
+ax = fig.add_subplot(projection=ss_br)
 
-# Plot the source surface map
-ss_br.plot()
+ss_br.plot(axes=ax)
 # Plot the polarity inversion line
 ax.plot_coord(pfss_out.source_surface_pils[0])
-# Plot formatting
 plt.colorbar()
 ax.set_title('Source surface magnetic field')
 
diff --git a/examples/using_sunkit_magex_pfss/plot_field_line_magnetic_field.py b/examples/using_sunkit_magex_pfss/magnetic_field_along_field_line.py
similarity index 98%
rename from examples/using_sunkit_magex_pfss/plot_field_line_magnetic_field.py
rename to examples/using_sunkit_magex_pfss/magnetic_field_along_field_line.py
index 7c18b6a..5aa4acf 100644
--- a/examples/using_sunkit_magex_pfss/plot_field_line_magnetic_field.py
+++ b/examples/using_sunkit_magex_pfss/magnetic_field_along_field_line.py
@@ -40,11 +40,10 @@
 # Now take a seed point, and trace a magnetic field line through the PFSS
 # solution from this point
 
-tracer = pfss.tracing.FortranTracer()
+tracer = pfss.tracing.PerformanceTracer()
 r = 1.2 * const.R_sun
 lat = 70 * u.deg
 lon = 0 * u.deg
-
 seeds = SkyCoord(lon, lat, r, frame=pfss_out.coordinate_frame)
 field_lines = tracer.trace(seeds, pfss_out)
 
@@ -59,6 +58,7 @@
 field_line = field_lines[0]
 B = field_line.b_along_fline
 r = field_line.coords.radius
+
 fig, ax = plt.subplots()
 
 ax.plot(r.to(const.R_sun), B[:, 0], label=r'$B_{r}$')
diff --git a/examples/using_sunkit_magex_pfss/plot_open_closed_map.py b/examples/using_sunkit_magex_pfss/open_closed_field_map.py
similarity index 89%
rename from examples/using_sunkit_magex_pfss/plot_open_closed_map.py
rename to examples/using_sunkit_magex_pfss/open_closed_field_map.py
index a5c41bd..af2f20c 100644
--- a/examples/using_sunkit_magex_pfss/plot_open_closed_map.py
+++ b/examples/using_sunkit_magex_pfss/open_closed_field_map.py
@@ -41,22 +41,19 @@
 #
 # First, set up the tracing seeds.
 
-r = const.R_sun
 # Number of steps in cos(latitude)
 nsteps = 45
 lon_1d = np.linspace(0, 2 * np.pi, nsteps * 2 + 1)
 lat_1d = np.arcsin(np.linspace(-1, 1, nsteps + 1))
 lon, lat = np.meshgrid(lon_1d, lat_1d, indexing='ij')
 lon, lat = lon*u.rad, lat*u.rad
-seeds = SkyCoord(lon.ravel(), lat.ravel(), r, frame=pfss_out.coordinate_frame)
+seeds = SkyCoord(lon.ravel(), lat.ravel(), const.R_sun, frame=pfss_out.coordinate_frame)
 
 ###############################################################################
 # Trace the field lines.
 
-print('Tracing field lines...')
-tracer = pfss.tracing.FortranTracer(max_steps=2000)
+tracer = pfss.tracing.PerformanceTracer()
 field_lines = tracer.trace(seeds, pfss_out)
-print('Finished tracing field lines')
 
 ###############################################################################
 # Plot the result. The to plot is the input magnetogram, and the bottom plot
@@ -64,9 +61,9 @@
 # field regions and 0 for closed field regions.
 
 fig = plt.figure()
-m = pfss_in.map
-ax = fig.add_subplot(2, 1, 1, projection=m)
-m.plot()
+input_map = pfss_in.map
+ax = fig.add_subplot(2, 1, 1, projection=input_map)
+input_map.plot(axes=ax)
 ax.set_title('Input GONG magnetogram')
 
 ax = fig.add_subplot(2, 1, 2)
@@ -75,8 +72,9 @@
 pols = field_lines.polarities.reshape(2 * nsteps + 1, nsteps + 1).T
 ax.contourf(np.rad2deg(lon_1d), np.sin(lat_1d), pols, norm=norm, cmap=cmap)
 ax.set_ylabel('sin(latitude)')
-
 ax.set_title('Open (blue/red) and closed (black) field')
 ax.set_aspect(0.5 * 360 / 2)
 
+fig.tight_layout()
+
 plt.show()
diff --git a/examples/using_sunkit_magex_pfss/plot_aia_overplotting.py b/examples/using_sunkit_magex_pfss/overplotting_onto_aia.py
similarity index 76%
rename from examples/using_sunkit_magex_pfss/plot_aia_overplotting.py
rename to examples/using_sunkit_magex_pfss/overplotting_onto_aia.py
index 164542b..8a5329d 100644
--- a/examples/using_sunkit_magex_pfss/plot_aia_overplotting.py
+++ b/examples/using_sunkit_magex_pfss/overplotting_onto_aia.py
@@ -5,7 +5,7 @@
 This example shows how to take a PFSS solution, trace some field lines, and
 overplot the traced field lines on an AIA 193 map.
 """
-# sphinx_gallery_thumbnail_number = 5
+# sphinx_gallery_thumbnail_number = 4
 
 import os
 
@@ -36,8 +36,7 @@
         'http://jsoc2.stanford.edu/data/aia/synoptic/2020/09/01/H1300/AIA20200901_1300_0193.fits',
         'aia_map.fits')
 
-aia = sunpy.map.Map('aia_map.fits')
-dtime = aia.date
+aia_map = sunpy.map.Map('aia_map.fits')
 
 ###############################################################################
 # The PFSS solution is calculated on a regular 3D grid in (phi, s, rho), where
@@ -53,26 +52,12 @@
 
 pfss_in = sunkit_magex.pfss.Input(gong_map, nrho, rss)
 
-###############################################################################
-# Using the `~sunkit_magex.pfss.Input` object, plot the input photospheric magnetic field.
-
-m = pfss_in.map
-fig = plt.figure()
-ax = plt.subplot(projection=m)
-m.plot()
-plt.colorbar()
-ax.set_title('Input field')
-
 ###############################################################################
 # We can also plot the AIA map to give an idea of the global picture. There
 # is a nice active region in the top right of the AIA plot, that can also
 # be seen in the top left of the photospheric field plot above.
-
-ax = plt.subplot(1, 1, 1, projection=aia)
-aia.plot(axes=ax)
-
-###############################################################################
-# Now we construct a 5 x 5 grid of footpoints to trace some magnetic field
+#
+# We construct a 5 x 5 grid of footpoints to trace some magnetic field
 # lines from. These coordinates are defined in the helioprojective frame of the
 # AIA image.
 
@@ -80,28 +65,28 @@
 hp_lat = np.linspace(250, 500, 5) * u.arcsec
 # Make a 2D grid from these 1D points
 lon, lat = np.meshgrid(hp_lon, hp_lat)
-seeds = SkyCoord(lon.ravel(), lat.ravel(), frame=aia.coordinate_frame)
+seeds = SkyCoord(lon.ravel(), lat.ravel(), frame=aia_map.coordinate_frame)
+
 fig = plt.figure()
-ax = plt.subplot(projection=aia)
-aia.plot(axes=ax)
+ax = fig.add_subplot(projection=aia_map)
+
+aia_map.plot(axes=ax)
 ax.plot_coord(seeds, color='white', marker='o', linewidth=0)
 
 ###############################################################################
-# Plot the magnetogram and the seed footpoints The footpoints are centered
+# Using the `~sunkit_magex.pfss.Input` object, we can plot the input photospheric
+# magnetic field and the seed footpoints The footpoints are centered
 # around the active region mentioned above.
 
-m = pfss_in.map
-fig = plt.figure()
-ax = plt.subplot(projection=m)
-m.plot()
-plt.colorbar()
+input_map = pfss_in.map
 
+fig = plt.figure()
+ax = fig.add_subplot(projection=input_map)
+input_map.plot(axes=ax)
 ax.plot_coord(seeds, color='black', marker='o', linewidth=0, markersize=2)
 
-# Set the axes limits. These limits have to be in pixel values
-# ax.set_xlim(0, 180)
-# ax.set_ylim(45, 135)
-ax.set_title('Field line footpoints')
+plt.colorbar()
+ax.set_title('Input field and seed footpoints')
 ax.set_ylim(bottom=0)
 
 #######################################################################
@@ -112,22 +97,20 @@
 ###############################################################################
 # Trace field lines from the footpoints defined above.
 
-tracer = tracing.FortranTracer()
+tracer = tracing.PythonTracer()
 flines = tracer.trace(seeds, pfss_out)
 
 ###############################################################################
 # Plot the input GONG magnetic field map, along with the traced magnetic field
 # lines.
 
-m = pfss_in.map
 fig = plt.figure()
-ax = plt.subplot(projection=m)
-m.plot()
-plt.colorbar()
+ax = fig.add_subplot(projection=input_map)
 
+input_map.plot(axes=ax)
+plt.colorbar()
 for fline in flines:
     ax.plot_coord(fline.coords, color='black', linewidth=1)
-
 ax.set_title('Photospheric field and traced field lines')
 
 ###############################################################################
@@ -136,8 +119,8 @@
 # and then plotted on top of the map.
 
 fig = plt.figure()
-ax = plt.subplot(1, 1, 1, projection=aia)
-aia.plot(axes=ax)
+ax = fig.add_subplot(1, 1, 1, projection=aia_map)
+aia_map.plot(axes=ax)
 for fline in flines:
     ax.plot_coord(fline.coords, alpha=0.8, linewidth=1, color='white')
 
diff --git a/examples/utils/README.rst b/examples/utils/README.rst
index a5bd0b3..18d41d0 100644
--- a/examples/utils/README.rst
+++ b/examples/utils/README.rst
@@ -1,3 +1,4 @@
 Utilities
 ---------
+
 Useful code that doesn't involve doing a PFSS extrapolation.
diff --git a/examples/utils/plot_car_reproject.py b/examples/utils/reproject_car_to_cea.py
similarity index 76%
rename from examples/utils/plot_car_reproject.py
rename to examples/utils/reproject_car_to_cea.py
index d7139e9..f77b892 100644
--- a/examples/utils/plot_car_reproject.py
+++ b/examples/utils/reproject_car_to_cea.py
@@ -12,13 +12,14 @@
 """
 import matplotlib.pyplot as plt
 
+import sunpy.map
+
 from sunkit_magex import pfss
 
 ###############################################################################
 # Load a sample ADAPT map, which has a CAR projection.
 
-adapt_maps = pfss.utils.load_adapt(pfss.sample_data.get_adapt_map())
-adapt_map_car = adapt_maps[0]
+adapt_map_car = sunpy.map.Map(pfss.sample_data.get_adapt_map(), hdus=0)
 
 ###############################################################################
 # Re-project into a CEA projection.
@@ -28,9 +29,13 @@
 ###############################################################################
 # Plot the original map and the reprojected map.
 
-plt.figure()
-adapt_map_car.plot()
-plt.figure()
-adapt_map_cea.plot()
+fig = plt.figure(figsize=(12, 6))
+ax = fig.add_subplot(1, 2, 1, projection=adapt_map_car)
+ax1 = fig.add_subplot(1, 2, 2, projection=adapt_map_cea)
+
+adapt_map_car.plot(axes=ax)
+adapt_map_cea.plot(axes=ax1)
+
+fig.tight_layout()
 
 plt.show()
diff --git a/pyproject.toml b/pyproject.toml
index aa4af83..3f1084a 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -13,8 +13,8 @@ requires-python = ">=3.10"
 readme = { file = "README.rst", content-type = "text/x-rst" }
 license = { file = "licenses/LICENSE.rst" }
 authors = [
-  { name = "David Stansby" },
   { name = "The SunPy Community", email = "sunpy@googlegroups.com" },
+  { name = "David Stansby" },
 ]
 dependencies = [
   "astropy>=5.3.0",
@@ -35,17 +35,17 @@ tests = [
   "pytest-arraydiff",
   "parfive",
   "reproject",
-  "sympy",
+  "sunkit-magex[all]",
 ]
 docs = [
   "sphinx",
   "sphinx-automodapi",
-  "pillow",
   "reproject",
   "sphinx-gallery",
-  "sunpy[net,map]",
-  "sympy",
+  "sunpy[net]>=6.0",
   "sunpy-sphinx-theme",
+  "sunkit-magex[all]",
+  "sphinx-changelog",
 ]
 performance = [
   "numba",
@@ -60,7 +60,7 @@ dev = ["sunkit-magex[all,tests,docs]"]
 Homepage = "https://sunpy.org"
 Download = "https://pypi.org/project/sunkit-magex/"
 "Source Code" = "https://github.com/sunpy/sunkit-magex/"
-Documentation = "https://docs.sunpy.org/"
+Documentation = "https://docs.sunpy.org/projects/sunkit-magex/en/stable/"
 Changelog = "https://docs.sunpy.org/en/stable/whatsnew/changelog.html"
 "Issue Tracker" = "https://github.com/sunpy/sunkit-magex/issues"
 
@@ -73,3 +73,47 @@ exclude = ["sunkit_magex._dev*"]
 
 [tool.setuptools_scm]
 write_to = "sunkit_magex/_version.py"
+
+# TODO: This should be in towncrier.toml but Giles currently only works looks in
+# pyproject.toml we should move this back when it's fixed.
+[tool.towncrier]
+  package = "sunkit_magex"
+  filename = "CHANGELOG.rst"
+  directory = "changelog/"
+  issue_format = "`#{issue} <https://github.com/sunkit_magex/sunkit_magex/pull/{issue}>`__"
+  title_format = "{version} ({project_date})"
+
+  [[tool.towncrier.type]]
+  directory = "breaking"
+  name = "Breaking Changes"
+  showcontent = true
+
+  [[tool.towncrier.type]]
+  directory = "deprecation"
+  name = "Deprecations"
+  showcontent = true
+
+  [[tool.towncrier.type]]
+  directory = "removal"
+  name = "Removals"
+  showcontent = true
+
+  [[tool.towncrier.type]]
+  directory = "feature"
+  name = "New Features"
+  showcontent = true
+
+  [[tool.towncrier.type]]
+  directory = "bugfix"
+  name = "Bug Fixes"
+  showcontent = true
+
+  [[tool.towncrier.type]]
+  directory = "doc"
+  name = "Documentation"
+  showcontent = true
+
+  [[tool.towncrier.type]]
+  directory = "trivial"
+  name = "Internal Changes"
+  showcontent = true
diff --git a/sunkit_magex/pfss/__init__.py b/sunkit_magex/pfss/__init__.py
index 56f5c3c..8ae25ab 100644
--- a/sunkit_magex/pfss/__init__.py
+++ b/sunkit_magex/pfss/__init__.py
@@ -1,5 +1,3 @@
-# Import this to register map sources
-import sunkit_magex.pfss.map as _
 from sunkit_magex.pfss import coords, fieldline, sample_data, tracing, utils
 from sunkit_magex.pfss.input import Input
 from sunkit_magex.pfss.output import Output
diff --git a/sunkit_magex/pfss/map.py b/sunkit_magex/pfss/map.py
deleted file mode 100644
index ce3ddeb..0000000
--- a/sunkit_magex/pfss/map.py
+++ /dev/null
@@ -1,24 +0,0 @@
-"""
-Custom `sunpy.map.GenericMap` sub-classes for different magnetogram sources.
-"""
-import sunpy.map
-
-__all__ = ['ADAPTMap']
-
-
-class ADAPTMap(sunpy.map.GenericMap):
-    def __init__(self, data, header, **kwargs):
-        if 'date-obs' not in header:
-            header['date-obs'] = header['maptime']
-        # Fix CTYPE
-        if header['ctype1'] == 'Long':
-            header['ctype1'] = 'CRLN-CAR'
-        if header['ctype2'] == 'Lat':
-            header['ctype2'] = 'CRLT-CAR'
-
-        super().__init__(data, header, **kwargs)
-
-    @classmethod
-    def is_datasource_for(cls, data, header, **kwargs):
-        """Determines if header corresponds to an ADAPT map."""
-        return header.get('model') == 'ADAPT'
diff --git a/sunkit_magex/pfss/tests/conftest.py b/sunkit_magex/pfss/tests/conftest.py
index 3ff7915..d935f8a 100644
--- a/sunkit_magex/pfss/tests/conftest.py
+++ b/sunkit_magex/pfss/tests/conftest.py
@@ -6,7 +6,6 @@
 from sunpy.map import Map
 
 import sunkit_magex.pfss
-from sunkit_magex.pfss import utils
 from sunkit_magex.tests.helpers import get_dummy_map_from_header, get_fitsfile_from_header
 
 
@@ -65,10 +64,12 @@ def adapt_test_file(tmp_path):
         package="sunkit_magex.pfss.tests.data"
     )
 
-
 @pytest.fixture
 def adapt_map(adapt_test_file):
-    return utils.load_adapt(adapt_test_file)
+    import warnings
+    with warnings.catch_warnings():
+        warnings.simplefilter("ignore")
+        return Map(adapt_test_file, hdus=0)
 
 
 @pytest.fixture
diff --git a/sunkit_magex/pfss/tests/test_map.py b/sunkit_magex/pfss/tests/test_map.py
deleted file mode 100644
index 0b994d4..0000000
--- a/sunkit_magex/pfss/tests/test_map.py
+++ /dev/null
@@ -1,12 +0,0 @@
-import astropy.io
-
-import sunpy.map
-
-import sunkit_magex.pfss.map
-
-
-def test_adapt_map(adapt_test_file):
-    with astropy.io.fits.open(adapt_test_file) as adapt_fits:
-        for map_slice in adapt_fits[0].data:
-            m = sunpy.map.Map((map_slice, adapt_fits[0].header))
-            assert isinstance(m, sunkit_magex.pfss.map.ADAPTMap)
diff --git a/sunkit_magex/pfss/tests/test_pfss.py b/sunkit_magex/pfss/tests/test_pfss.py
index c06252e..2f3a948 100644
--- a/sunkit_magex/pfss/tests/test_pfss.py
+++ b/sunkit_magex/pfss/tests/test_pfss.py
@@ -55,7 +55,7 @@ def test_bunit(gong_map):
 
 
 def test_expansion_factor(dipole_result):
-    inp, out = dipole_result
+    _, out = dipole_result
     out_frame = out.coordinate_frame
 
     tracer = tracing.PythonTracer()
@@ -80,7 +80,7 @@ def test_expansion_factor(dipole_result):
 
 
 def test_field_line_polarity(dipole_result):
-    input, out = dipole_result
+    _, out = dipole_result
     out_frame = out.coordinate_frame
 
     tracer = tracing.PythonTracer()
@@ -99,7 +99,7 @@ def test_field_line_polarity(dipole_result):
 
 
 def test_footpoints(dipole_result):
-    input, out = dipole_result
+    _, out = dipole_result
     out_frame = out.coordinate_frame
 
     tracer = tracing.PythonTracer(atol=1e-8, rtol=1e-8)
@@ -128,12 +128,12 @@ def check_open_fline(fline):
 
 def test_shape(zero_map):
     # Test output map shapes
-    input, out = zero_map
-    nr = input.grid.nr
-    nphi = input.grid.nphi
-    ns = input.grid.ns
+    _input, out = zero_map
+    nr = _input.grid.nr
+    nphi = _input.grid.nphi
+    ns = _input.grid.ns
 
-    out = sunkit_magex.pfss.pfss(input)
+    out = sunkit_magex.pfss.pfss(_input)
     alr, als, alp = out._al
     for comp in (alr, als, alp):
         assert np.all(comp == 0)
diff --git a/sunkit_magex/pfss/tests/test_tracers.py b/sunkit_magex/pfss/tests/test_tracers.py
index d617aea..ac8b8de 100644
--- a/sunkit_magex/pfss/tests/test_tracers.py
+++ b/sunkit_magex/pfss/tests/test_tracers.py
@@ -10,7 +10,7 @@
 
 
 @pytest.fixture(params=[tracing.PythonTracer(),
-                        tracing.FortranTracer()],
+                        tracing.PerformanceTracer()],
                 ids=['python', 'compiled'])
 def flines(dipole_result, request):
     tracer = request.param
@@ -43,13 +43,13 @@ def test_fline_step_size(dipole_result):
     seed = coord.SkyCoord(2*u.deg, -45*u.deg, 1.01*const.R_sun,
                           frame=out_frame)
 
-    tracer = tracing.FortranTracer(step_size=0.5)
+    tracer = tracing.PerformanceTracer(step_size=0.5)
     flines = tracer.trace(seed, out)
     assert out.grid.nr == 10
     # With a step size of 0.5, this should be ~20
     assert len(flines[0]) == 21
 
-    tracer = tracing.FortranTracer(step_size=0.2)
+    tracer = tracing.PerformanceTracer(step_size=0.2)
     flines = tracer.trace(seed, out)
     assert out.grid.nr == 10
     # With a step size of 0.2, this should be ~50
@@ -57,7 +57,7 @@ def test_fline_step_size(dipole_result):
 
 
 def test_rot_warning(dipole_result):
-    tracer = tracing.FortranTracer(max_steps=2)
+    tracer = tracing.PerformanceTracer(max_steps=2)
     input, out = dipole_result
     out_frame = out.coordinate_frame
     seed = coord.SkyCoord(0*u.deg, -45*u.deg, 1.01*const.R_sun,
diff --git a/sunkit_magex/pfss/tests/test_utils.py b/sunkit_magex/pfss/tests/test_utils.py
index e66b574..c7e1302 100644
--- a/sunkit_magex/pfss/tests/test_utils.py
+++ b/sunkit_magex/pfss/tests/test_utils.py
@@ -7,17 +7,13 @@
 
 import sunpy.map
 
-import sunkit_magex.pfss
-from sunkit_magex.pfss import utils
-
-# Ignore missing metadata warnings
-pytestmark = [pytest.mark.filterwarnings('ignore:Missing metadata for observer')]
-
-def test_load_adapt(adapt_test_file):
-    adaptMapSequence = utils.load_adapt(adapt_test_file)
-    assert isinstance(adaptMapSequence, sunpy.map.MapSequence)
-    for map_ in adaptMapSequence:
-        assert map_.meta['model'] == "ADAPT"
+from sunkit_magex.pfss.utils import (
+    car_to_cea,
+    carr_cea_wcs_header,
+    is_cea_map,
+    is_full_sun_synoptic_map,
+    roll_map,
+)
 
 
 def test_header_generation():
@@ -25,7 +21,7 @@ def test_header_generation():
     nphi = 90
     dtime = '2001-01-01 00:00:00'
     shape = [nphi, ntheta]
-    header = sunkit_magex.pfss.utils.carr_cea_wcs_header(dtime, shape)
+    header = carr_cea_wcs_header(dtime, shape)
     assert header['LONPOLE'] == 0
     assert header['CTYPE1'] == 'CRLN-CEA'
     assert header['CTYPE2'] == 'CRLT-CEA'
@@ -58,20 +54,19 @@ def test_header_generation():
 
 @pytest.mark.parametrize('error', [True, False])
 def test_validation(dipole_map, error):
-    assert utils.is_cea_map(dipole_map, error)
-    assert utils.is_full_sun_synoptic_map(dipole_map, error)
+    assert is_cea_map(dipole_map, error)
+    assert is_full_sun_synoptic_map(dipole_map, error)
 
 
 def test_validation_not_full_map(dipole_map):
     dipole_map.meta['cdelt1'] = 0.001
-    assert not utils.is_full_sun_synoptic_map(dipole_map)
+    assert not is_full_sun_synoptic_map(dipole_map)
     with pytest.raises(ValueError, match='Number of points in phi direction times'):
-        utils.is_full_sun_synoptic_map(dipole_map, error=True)
+        is_full_sun_synoptic_map(dipole_map, error=True)
 
 
-def test_car_reproject(adapt_test_file):
-    adapt_map = utils.load_adapt(adapt_test_file)[0]
-    adapt_reproj = utils.car_to_cea(adapt_map)
+def test_car_reproject(adapt_map):
+    adapt_reproj = car_to_cea(adapt_map)
 
     assert np.all(np.isfinite(adapt_map.data))
     assert np.all(np.isfinite(adapt_reproj.data))
@@ -81,13 +76,13 @@ def test_car_reproject(adapt_test_file):
         assert adapt_reproj.meta[f'CTYPE{i}'][5:8] == 'CEA'
 
     with pytest.raises(ValueError, match='method must be one of'):
-        utils.car_to_cea(adapt_map, method='gibberish')
+        car_to_cea(adapt_map, method='gibberish')
 
 
 def test_roll_map(adapt_map, gong_map):
     lh_edge_test = 0.0 * u.deg
     gong_map = sunpy.map.Map(gong_map)
-    rolled_map = utils.roll_map(gong_map,
+    rolled_map = roll_map(gong_map,
                                 lh_edge_lon=lh_edge_test)
 
     # Test ref pixel rolled correctly
@@ -99,38 +94,38 @@ def test_roll_map(adapt_map, gong_map):
     assert np.all(np.isfinite(rolled_map.data))
 
     # Test output map is full sun synoptic
-    assert utils.is_full_sun_synoptic_map(rolled_map, error=True)
+    assert is_full_sun_synoptic_map(rolled_map, error=True)
 
     # Test reproject method error handling
     with pytest.raises(ValueError, match='method must be one of'):
-        utils.roll_map(gong_map, method='gibberish')
+        roll_map(gong_map, method='gibberish')
 
     # Test left hand edge input type validation
     # 1. No Units
     with pytest.raises(TypeError,
                        match="has no 'unit' attribute"):
-        utils.roll_map(adapt_map, lh_edge_lon=0)
+        roll_map(adapt_map, lh_edge_lon=0)
     # 2. Incompatible units
     with pytest.raises(u.UnitsError,
                        match="must be in units convertible to 'deg'"):
-        utils.roll_map(adapt_map, lh_edge_lon=0 * u.m)
+        roll_map(adapt_map, lh_edge_lon=0 * u.m)
 
     # Test left hand edge input range validation
     with pytest.raises(ValueError,
                        match='lh_edge_lon must be in'):
-        utils.roll_map(adapt_map, lh_edge_lon=361 * u.deg)
+        roll_map(adapt_map, lh_edge_lon=361 * u.deg)
 
 
 def test_cea_header():
     # Assert default reference pixel is at 0 deg lon
-    cea_default = utils.carr_cea_wcs_header(
+    cea_default = carr_cea_wcs_header(
         datetime.datetime(2020, 1, 1),
         [360, 180]
     )
     assert cea_default['crval1'] == 0.0
 
     # Assert custom reference pixel is expected lon
-    cea_shift = utils.carr_cea_wcs_header(
+    cea_shift = carr_cea_wcs_header(
         datetime.datetime(2020, 1, 1),
         [360, 180],
         map_center_longitude=10.0*u.deg
@@ -140,14 +135,14 @@ def test_cea_header():
     # Test reference pixel shift error handling
     # 1: No units
     with pytest.raises(u.UnitTypeError):
-        cea_default = utils.carr_cea_wcs_header(
+        cea_default = carr_cea_wcs_header(
             datetime.datetime(2020, 1, 1),
             [360, 180],
             map_center_longitude=0.0
         )
     # 2: Wrong Units
     with pytest.raises(u.UnitTypeError):
-        cea_default = utils.carr_cea_wcs_header(
+        cea_default = carr_cea_wcs_header(
             datetime.datetime(2020, 1, 1),
             [360, 180],
             map_center_longitude=0.0*u.m
diff --git a/sunkit_magex/pfss/tracing.py b/sunkit_magex/pfss/tracing.py
index 046483c..8ab80fc 100644
--- a/sunkit_magex/pfss/tracing.py
+++ b/sunkit_magex/pfss/tracing.py
@@ -2,14 +2,18 @@
 import warnings
 
 import numpy as np
+from streamtracer import StreamTracer, VectorGrid
 
 import astropy.constants as const
 import astropy.coordinates as astrocoords
 import astropy.units as u
 
+from sunpy.util.decorators import deprecated
+
 import sunkit_magex.pfss
 import sunkit_magex.pfss.fieldline as fieldline
 
+__all__ = ['Tracer', 'PerformanceTracer', 'FortranTracer', 'PythonTracer']
 
 class Tracer(abc.ABC):
     """
@@ -37,8 +41,9 @@ def validate_seeds(seeds):
         Check that *seeds* has the right shape and is the correct type.
         """
         if not isinstance(seeds, astrocoords.SkyCoord):
-            raise ValueError('seeds must be SkyCoord object '
-                             f'(got {type(seeds)} instead)')
+            raise ValueError(
+                f'seeds must be SkyCoord object (got {type(seeds)} instead)'
+            )
 
     @staticmethod
     def cartesian_to_coordinate():
@@ -73,7 +78,7 @@ def coords_to_xyz(seeds, output):
         return x, y, z
 
 
-class FortranTracer(Tracer):
+class PerformanceTracer(Tracer):
     r"""
     Tracer using compiled code via streamtracer.
 
@@ -99,12 +104,6 @@ class FortranTracer(Tracer):
     directly on the poles will not go anywhere.
     """
     def __init__(self, max_steps='auto', step_size=1):
-        try:
-            from streamtracer import StreamTracer
-        except ModuleNotFoundError as e:
-            raise RuntimeError(
-                'Using this tracer requires the streamtracer module, '
-                'but streamtracer could not be loaded') from e
         self.max_steps = max_steps
         self.step_size = step_size
         self.max_steps = max_steps
@@ -116,8 +115,6 @@ def vector_grid(output):
         """
         Create a `streamtracer.VectorGrid` object from an `~sunkit_magex.pfss.Output`.
         """
-        from streamtracer import VectorGrid
-
         # The indexing order on the last index is (phi, s, r)
         vectors = output.bg.copy()
 
@@ -172,9 +169,7 @@ def trace(self, seeds, output):
 
         # Filter out of bounds points out
         rho_ss = np.log(output.grid.rss)
-        xs = [x[(x[:, 2] <= rho_ss) & (x[:, 2] >= 0) &
-                (np.abs(x[:, 1]) < 1), :]
-              for x in xs]
+        xs = [x[(x[:, 2] <= rho_ss) & (x[:, 2] >= 0) & (np.abs(x[:, 1]) < 1), :] for x in xs]
 
         rots = self.tracer.ROT
         if np.any(rots == 1):
@@ -191,7 +186,6 @@ def trace(self, seeds, output):
 class PythonTracer(Tracer):
     """
     Tracer using native python code.
-
     Uses `scipy.integrate.solve_ivp`, with an LSODA method.
     """
     def __init__(self, atol=1e-4, rtol=1e-4):
@@ -219,3 +213,9 @@ def trace(self, seeds, output):
 
             flines.append(fline)
         return fieldline.FieldLines(flines)
+
+
+class FortranTracer(PerformanceTracer):
+    @deprecated('1.0', message='This class has been renamed to PerformanceTracer as this now does not use Fortran.')
+    def __init__(*args, **kwargs):
+        super().__init__(*args, **kwargs)
diff --git a/sunkit_magex/pfss/utils.py b/sunkit_magex/pfss/utils.py
index 2e850f2..4f528bc 100644
--- a/sunkit_magex/pfss/utils.py
+++ b/sunkit_magex/pfss/utils.py
@@ -12,6 +12,7 @@
 import sunpy.coordinates
 import sunpy.map
 import sunpy.time
+from sunpy.util.decorators import deprecated
 
 __all__ = ['fix_hmi_meta', 'load_adapt', 'carr_cea_wcs_header', 'is_cea_map', 'is_car_map', 'is_full_sun_synoptic_map', 'car_to_cea', 'roll_map']
 
@@ -36,6 +37,7 @@ def _earth_obs_coord_meta(obstime):
     return _observer_coord_meta(sunpy.coordinates.get_earth(obstime))
 
 
+@deprecated('1.0', message="This is now fixed within sunpy when you load a HMI file.")
 def fix_hmi_meta(hmi_map):
     """
     Fix non-compliant FITS metadata in HMI maps.
@@ -79,6 +81,7 @@ def fix_hmi_meta(hmi_map):
         hmi_map.meta.update(_earth_obs_coord_meta(hmi_map.meta['date-obs']))
 
 
+@deprecated('1.0', message="This is not required anymore. An example of how to load an ADAPT file is available in the documentation.")
 def load_adapt(adapt_path):
     """
     Parse adapt .fts file as a `sunpy.map.MapSequence`
@@ -150,7 +153,7 @@ def carr_cea_wcs_header(dtime, shape, *, map_center_longitude=0*u.deg):
 
 
 def _get_projection(m, i):
-    return m.meta[f'ctype{i}'][5:8]
+    return m.coordinate_system.axis1.split("-")[-1]
 
 
 def _check_projection(m, proj_code, error=False):