diff --git a/autogalaxy/aggregator/agg_util.py b/autogalaxy/aggregator/agg_util.py
index da3fa69e3..035f3947d 100644
--- a/autogalaxy/aggregator/agg_util.py
+++ b/autogalaxy/aggregator/agg_util.py
@@ -11,7 +11,7 @@
 from autogalaxy.analysis.adapt_images.adapt_images import AdaptImages
 
 
-def mask_header_from(fit):
+def mask_header_from(fit, name="dataset"):
     """
     Returns the mask, header and pixel scales of the `PyAutoFit` `Fit` object.
 
@@ -30,7 +30,7 @@ def mask_header_from(fit):
     The mask, header and pixel scales of the `PyAutoFit` `Fit` object.
     """
 
-    header = aa.Header(header_sci_obj=fit.value(name="dataset")[0].header)
+    header = aa.Header(header_sci_obj=fit.value(name=name)[0].header)
     pixel_scales = (
         header.header_sci_obj[Mask2DKeys.PIXSCAY.value],
         header.header_sci_obj[Mask2DKeys.PIXSCAY.value],
@@ -40,7 +40,7 @@ def mask_header_from(fit):
         header.header_sci_obj[Mask2DKeys.ORIGINX.value],
     )
     mask = aa.Mask2D(
-        mask=ndarray_via_hdu_from(fit.value(name="dataset")[0]),
+        mask=ndarray_via_hdu_from(fit.value(name=name)[0]),
         pixel_scales=pixel_scales,
         origin=origin,
     )
@@ -79,7 +79,7 @@ def adapt_images_from(
     adapt_images_list = []
 
     for fit in fit_list:
-        mask, header = mask_header_from(fit=fit)
+        mask, header = mask_header_from(fit=fit, name="adapt_images")
 
         galaxy_name_image_dict = {}
 
diff --git a/autogalaxy/aggregator/subplot.py b/autogalaxy/aggregator/subplot.py
index 30aa03763..54fdd3546 100644
--- a/autogalaxy/aggregator/subplot.py
+++ b/autogalaxy/aggregator/subplot.py
@@ -6,7 +6,7 @@ class FITSFit(Enum):
     The HDUs that can be extracted from the fit.fits file.
     """
 
-    model_image = "MODEL_IMAGE"
+    model_data = "MODEL_DATA"
     residual_map = "RESIDUAL_MAP"
     normalized_residual_map = "NORMALIZED_RESIDUAL_MAP"
     chi_squared_map = "CHI_SQUARED_MAP"
@@ -38,7 +38,7 @@ class SubplotFit(Enum):
 
     data = (0, 0)
     signal_to_noise_map = (1, 0)
-    model_image = (2, 0)
+    model_data = (2, 0)
     normalized_residual_map = (0, 1)
     normalized_residual_map_one_sigma = (1, 1)
     chi_squared_map = (2, 1)
diff --git a/autogalaxy/analysis/analysis/analysis.py b/autogalaxy/analysis/analysis/analysis.py
index ceeecbbb6..7b87a910e 100644
--- a/autogalaxy/analysis/analysis/analysis.py
+++ b/autogalaxy/analysis/analysis/analysis.py
@@ -192,33 +192,33 @@ def profile_log_likelihood_function(
 
         try:
             info_dict["image_pixels"] = self.dataset.grids.lp.shape_slim
-            info_dict[
-                "sub_total_light_profiles"
-            ] = self.dataset.grids.lp.over_sampler.sub_total
+            info_dict["sub_total_light_profiles"] = (
+                self.dataset.grids.lp.over_sampler.sub_total
+            )
         except AttributeError:
             pass
 
         if fit.model_obj.has(cls=aa.Pixelization):
             info_dict["use_w_tilde"] = fit.inversion.settings.use_w_tilde
             try:
-                info_dict[
-                    "sub_total_pixelization"
-                ] = self.dataset.grids.pixelization.over_sampler.sub_total
+                info_dict["sub_total_pixelization"] = (
+                    self.dataset.grids.pixelization.over_sampler.sub_total
+                )
             except AttributeError:
                 pass
-            info_dict[
-                "use_positive_only_solver"
-            ] = fit.inversion.settings.use_positive_only_solver
-            info_dict[
-                "force_edge_pixels_to_zeros"
-            ] = fit.inversion.settings.force_edge_pixels_to_zeros
+            info_dict["use_positive_only_solver"] = (
+                fit.inversion.settings.use_positive_only_solver
+            )
+            info_dict["force_edge_pixels_to_zeros"] = (
+                fit.inversion.settings.force_edge_pixels_to_zeros
+            )
             info_dict["use_w_tilde_numpy"] = fit.inversion.settings.use_w_tilde_numpy
             info_dict["source_pixels"] = len(fit.inversion.reconstruction)
 
             if hasattr(fit.inversion, "w_tilde"):
-                info_dict[
-                    "w_tilde_curvature_preload_size"
-                ] = fit.inversion.w_tilde.curvature_preload.shape[0]
+                info_dict["w_tilde_curvature_preload_size"] = (
+                    fit.inversion.w_tilde.curvature_preload.shape[0]
+                )
 
         self.output_profiling_info(
             paths=paths, run_time_dict=run_time_dict, info_dict=info_dict
diff --git a/autogalaxy/analysis/chaining_util.py b/autogalaxy/analysis/chaining_util.py
index 930b6a5f6..95306e45f 100644
--- a/autogalaxy/analysis/chaining_util.py
+++ b/autogalaxy/analysis/chaining_util.py
@@ -215,39 +215,29 @@ def extra_galaxies_from(
 
         for extra_galaxy_index in range(len(result.instance.extra_galaxies)):
             if hasattr(result.instance.extra_galaxies[extra_galaxy_index], "mass"):
-                extra_galaxies[
-                    extra_galaxy_index
-                ].mass.centre = result.instance.extra_galaxies[
-                    extra_galaxy_index
-                ].mass.centre
-                extra_galaxies[
-                    extra_galaxy_index
-                ].mass.einstein_radius = result.model.extra_galaxies[
-                    extra_galaxy_index
-                ].mass.einstein_radius
+                extra_galaxies[extra_galaxy_index].mass.centre = (
+                    result.instance.extra_galaxies[extra_galaxy_index].mass.centre
+                )
+                extra_galaxies[extra_galaxy_index].mass.einstein_radius = (
+                    result.model.extra_galaxies[extra_galaxy_index].mass.einstein_radius
+                )
                 if free_centre:
-                    extra_galaxies[
-                        extra_galaxy_index
-                    ].mass.centre = result.model.extra_galaxies[
-                        extra_galaxy_index
-                    ].mass.centre
+                    extra_galaxies[extra_galaxy_index].mass.centre = (
+                        result.model.extra_galaxies[extra_galaxy_index].mass.centre
+                    )
 
     elif light_as_model:
         extra_galaxies = result.instance.extra_galaxies.as_model((LightProfile,))
 
         for extra_galaxy_index in range(len(result.instance.extra_galaxies)):
             if extra_galaxies[extra_galaxy_index].bulge is not None:
-                extra_galaxies[
-                    extra_galaxy_index
-                ].bulge.centre = result.instance.extra_galaxies[
-                    extra_galaxy_index
-                ].bulge.centre
+                extra_galaxies[extra_galaxy_index].bulge.centre = (
+                    result.instance.extra_galaxies[extra_galaxy_index].bulge.centre
+                )
                 if free_centre:
-                    extra_galaxies[
-                        extra_galaxy_index
-                    ].bulge.centre = result.model.extra_galaxies[
-                        extra_galaxy_index
-                    ].bulge.centre
+                    extra_galaxies[extra_galaxy_index].bulge.centre = (
+                        result.model.extra_galaxies[extra_galaxy_index].bulge.centre
+                    )
 
     else:
         extra_galaxies = result.instance.extra_galaxies.as_model(())
diff --git a/autogalaxy/analysis/plotter_interface.py b/autogalaxy/analysis/plotter_interface.py
index f7285bd2a..ddec7258f 100644
--- a/autogalaxy/analysis/plotter_interface.py
+++ b/autogalaxy/analysis/plotter_interface.py
@@ -1,3 +1,4 @@
+import csv
 import os
 from pathlib import Path
 from typing import List, Union
@@ -11,7 +12,6 @@
 from autogalaxy.analysis.adapt_images.adapt_images import AdaptImages
 from autogalaxy.galaxy.galaxy import Galaxy
 from autogalaxy.galaxy.galaxies import Galaxies
-from autogalaxy.galaxy.plot.galaxy_plotters import GalaxyPlotter
 from autogalaxy.galaxy.plot.galaxies_plotters import GalaxiesPlotter
 from autogalaxy.galaxy.plot.adapt_plotters import AdaptPlotter
 
@@ -180,9 +180,13 @@ def should_plot(name):
             pass
 
         if should_plot("fits_galaxy_images"):
+
+            image_list = [
+                galaxy.image_2d_from(grid=grid).native_for_fits for galaxy in galaxies
+            ]
+
             hdu_list = hdu_list_for_output_from(
-                values_list=[grid.mask.astype("float")]
-                + [galaxy.image_2d_from(grid=grid) for galaxy in galaxies],
+                values_list=[image_list[0].mask.astype("float")] + image_list,
                 ext_name_list=["mask"] + [f"galaxy_{i}" for i in range(len(galaxies))],
                 header_dict=grid.mask.header_dict,
             )
@@ -221,8 +225,32 @@ def should_plot(name):
         if should_plot("subplot_inversion"):
             mapper_list = inversion.cls_list_from(cls=aa.AbstractMapper)
 
-            for mapper_index in range(len(mapper_list)):
-                inversion_plotter.subplot_of_mapper(mapper_index=mapper_index)
+            for i in range(len(mapper_list)):
+                suffix = "" if len(mapper_list) == 1 else f"_{i}"
+
+                inversion_plotter.subplot_of_mapper(
+                    mapper_index=i, auto_filename=f"subplot_inversion{suffix}"
+                )
+
+        if should_plot("csv_reconstruction"):
+            mapper_list = inversion.cls_list_from(cls=aa.AbstractMapper)
+
+            for i, mapper in enumerate(mapper_list):
+                y = mapper.mapper_grids.source_plane_mesh_grid[:, 0]
+                x = mapper.mapper_grids.source_plane_mesh_grid[:, 1]
+                reconstruction = inversion.reconstruction_dict[mapper]
+                noise_map = inversion.reconstruction_noise_map_dict[mapper]
+
+                with open(
+                    self.image_path / f"inversion_reconstruction_{i}.csv",
+                    mode="w",
+                    newline="",
+                ) as file:
+                    writer = csv.writer(file)
+                    writer.writerow(["y", "x", "reconstruction", "noise_map"])  # header
+
+                    for i in range(len(x)):
+                        writer.writerow([y[i], x[i], reconstruction[i], noise_map[i]])
 
     def adapt_images(
         self,
@@ -260,16 +288,16 @@ def should_plot(name):
             )
 
         if should_plot("fits_adapt_images"):
-            values_list = [
-                adapt_images.galaxy_name_image_dict[name].native
+            image_list = [
+                adapt_images.galaxy_name_image_dict[name].native_for_fits
                 for name in adapt_images.galaxy_name_image_dict.keys()
             ]
 
             hdu_list = hdu_list_for_output_from(
                 values_list=[
-                    adapt_images.mask.astype("float"),
+                    image_list[0].mask.astype("float"),
                 ]
-                + values_list,
+                + image_list,
                 ext_name_list=["mask"]
                 + list(adapt_images.galaxy_name_image_dict.keys()),
                 header_dict=adapt_images.mask.header_dict,
diff --git a/autogalaxy/config/visualize/plots.yaml b/autogalaxy/config/visualize/plots.yaml
index c66b85f45..5a72a6e52 100644
--- a/autogalaxy/config/visualize/plots.yaml
+++ b/autogalaxy/config/visualize/plots.yaml
@@ -11,6 +11,7 @@
 # These can be disabled to save on hard-disk space but will lead to certain database functionality being disabled.
 
 subplot_format: [png]                      # Output format of all subplots, can be png, pdf or both (e.g. [png, pdf])
+fits_are_zoomed: true                      # If true, output .fits files are zoomed in on the center of the unmasked region image, saving hard-disk space.
 
 dataset:                                   # Settings for plots of all datasets (e.g. ImagingPlotter, InterferometerPlotter).
   subplot_dataset: true                    # Plot subplot containing all dataset quantities (e.g. the data, noise-map, etc.)?
@@ -36,7 +37,7 @@ galaxies:                                  # Settings for plots of galaxies (e.g
 inversion:                                 # Settings for plots of inversions (e.g. InversionPlotter).
   subplot_inversion: true                  # Plot subplot of all quantities in each inversion (e.g. reconstrucuted image, reconstruction)?
   subplot_mappings: true                   # Plot subplot of the image-to-source pixels mappings of each pixelization?
-  fits_reconstruction: false               # output reconstruction.fits containing the reconstructed pixelization and noise map on the adaptive mesh?
+  csv_reconstruction: false                # output reconstruction.csv containing the source-plane mesh y, x, reconstruction and noise map values.
 
 adapt:                                     # Settings for plots of adapt images used by adaptive pixelizations.
   subplot_adapt_images: true               # Plot subplot showing each adapt image used for adaptive pixelization?
diff --git a/autogalaxy/ellipse/model/plotter_interface.py b/autogalaxy/ellipse/model/plotter_interface.py
index 1c6767eac..646e2c3b3 100644
--- a/autogalaxy/ellipse/model/plotter_interface.py
+++ b/autogalaxy/ellipse/model/plotter_interface.py
@@ -44,12 +44,13 @@ def should_plot(name):
         if should_plot("subplot_dataset"):
             dataset_plotter.subplot_dataset()
 
+        image_list = [
+            dataset.data.native_for_fits,
+            dataset.noise_map.native_for_fits,
+        ]
+
         hdu_list = hdu_list_for_output_from(
-            values_list=[
-                dataset.mask.astype("float"),
-                dataset.data.native,
-                dataset.noise_map.native,
-            ],
+            values_list=[image_list[0].mask.astype("float")] + image_list,
             ext_name_list=[
                 "mask",
                 "data",
diff --git a/autogalaxy/imaging/fit_imaging.py b/autogalaxy/imaging/fit_imaging.py
index 4fba3a0de..08819dfe1 100644
--- a/autogalaxy/imaging/fit_imaging.py
+++ b/autogalaxy/imaging/fit_imaging.py
@@ -220,9 +220,9 @@ def subtracted_images_of_galaxies_dict(self) -> Dict[Galaxy, aa.Array2D]:
 
             subtracted_image = self.data - sum(other_galaxies_model_images)
 
-            subtracted_images_of_galaxies_dict[
-                self.galaxies[galaxy_index]
-            ] = subtracted_image
+            subtracted_images_of_galaxies_dict[self.galaxies[galaxy_index]] = (
+                subtracted_image
+            )
 
         return subtracted_images_of_galaxies_dict
 
diff --git a/autogalaxy/imaging/model/plotter_interface.py b/autogalaxy/imaging/model/plotter_interface.py
index 299c8f478..97b91c42e 100644
--- a/autogalaxy/imaging/model/plotter_interface.py
+++ b/autogalaxy/imaging/model/plotter_interface.py
@@ -35,14 +35,19 @@ def fits_to_fits(
     """
 
     if should_plot("fits_fit"):
+
+        image_list = [
+            fit.model_data.native_for_fits,
+            fit.residual_map.native_for_fits,
+            fit.normalized_residual_map.native_for_fits,
+            fit.chi_squared_map.native_for_fits,
+        ]
+
         hdu_list = hdu_list_for_output_from(
             values_list=[
-                fit.mask.astype("float"),
-                fit.model_data,
-                fit.residual_map,
-                fit.normalized_residual_map,
-                fit.chi_squared_map,
-            ],
+                image_list[0].mask.astype("float"),
+            ]
+            + image_list,
             ext_name_list=[
                 "mask",
                 "model_data",
@@ -58,9 +63,12 @@ def fits_to_fits(
     if should_plot("fits_model_galaxy_images"):
         number_plots = len(fit.galaxy_model_image_dict.keys()) + 1
 
+        image_list = [
+            image.native_for_fits for image in fit.galaxy_model_image_dict.values()
+        ]
+
         hdu_list = hdu_list_for_output_from(
-            values_list=[fit.mask.astype("float")]
-            + [image.native for image in fit.galaxy_model_image_dict.values()],
+            values_list=[image_list[0].mask.astype("float")] + image_list,
             ext_name_list=[
                 "mask",
             ]
@@ -104,15 +112,16 @@ def should_plot(name):
             dataset_plotter.subplot_dataset()
 
         if should_plot("fits_dataset"):
+            image_list = [
+                dataset.data.native_for_fits,
+                dataset.noise_map.native_for_fits,
+                dataset.psf.native_for_fits,
+                dataset.grids.lp.over_sample_size.native_for_fits,
+                dataset.grids.pixelization.over_sample_size.native_for_fits,
+            ]
+
             hdu_list = hdu_list_for_output_from(
-                values_list=[
-                    dataset.mask.astype("float"),
-                    dataset.data.native,
-                    dataset.noise_map.native,
-                    dataset.psf.native,
-                    dataset.grids.lp.over_sample_size.native,
-                    dataset.grids.pixelization.over_sample_size.native,
-                ],
+                values_list=[image_list[0].mask.astype("float")] + image_list,
                 ext_name_list=[
                     "mask",
                     "data",
diff --git a/autogalaxy/interferometer/model/plotter_interface.py b/autogalaxy/interferometer/model/plotter_interface.py
index bbf478650..a37096fbc 100644
--- a/autogalaxy/interferometer/model/plotter_interface.py
+++ b/autogalaxy/interferometer/model/plotter_interface.py
@@ -36,9 +36,13 @@ def fits_to_fits(
     """
 
     if should_plot("fits_model_galaxy_images"):
+
+        image_list = [
+            image.native_for_fits for image in fit.galaxy_model_image_dict.values()
+        ]
+
         hdu_list = hdu_list_for_output_from(
-            values_list=[fit.dataset.real_space_mask.astype("float")]
-            + [image.native for image in fit.galaxy_model_image_dict.values()],
+            values_list=[image_list[0].mask.astype("float")] + image_list,
             ext_name_list=["mask"]
             + [f"galaxy_{i}" for i in range(len(fit.galaxy_model_image_dict.values()))],
             header_dict=fit.dataset.real_space_mask.header_dict,
@@ -47,16 +51,18 @@ def fits_to_fits(
         hdu_list.writeto(image_path / "model_galaxy_images.fits", overwrite=True)
 
     if should_plot("fits_dirty_images"):
+
+        image_list = [
+            fit.dirty_image.native_for_fits,
+            fit.dirty_noise_map.native_for_fits,
+            fit.dirty_model_image.native_for_fits,
+            fit.dirty_residual_map.native_for_fits,
+            fit.dirty_normalized_residual_map.native_for_fits,
+            fit.dirty_chi_squared_map.native_for_fits,
+        ]
+
         hdu_list = hdu_list_for_output_from(
-            values_list=[
-                fit.dataset.real_space_mask.astype("float"),
-                fit.dirty_image.native,
-                fit.dirty_noise_map.native,
-                fit.dirty_model_image.native,
-                fit.dirty_residual_map.native,
-                fit.dirty_normalized_residual_map.native,
-                fit.dirty_chi_squared_map.native,
-            ],
+            values_list=[image_list[0].mask.astype("float")] + image_list,
             ext_name_list=["mask"]
             + [
                 "dirty_image",
@@ -109,6 +115,7 @@ def should_plot(name):
             dataset_plotter.subplot_dataset()
 
         if should_plot("fits_dataset"):
+
             hdu_list = hdu_list_for_output_from(
                 values_list=[
                     dataset.real_space_mask.astype("float"),
diff --git a/autogalaxy/operate/deflections.py b/autogalaxy/operate/deflections.py
index 28f01958c..002737e9c 100644
--- a/autogalaxy/operate/deflections.py
+++ b/autogalaxy/operate/deflections.py
@@ -56,7 +56,9 @@ def wrapper(
 
         pixel_scale_ratio = grid.pixel_scale / pixel_scale
 
-        zoom_shape_native = grid.mask.zoom_shape_native
+        zoom = aa.Zoom2D(mask=grid.mask)
+
+        zoom_shape_native = zoom.shape_native
         shape_native = (
             int(pixel_scale_ratio * zoom_shape_native[0]),
             int(pixel_scale_ratio * zoom_shape_native[1]),
@@ -78,7 +80,7 @@ def wrapper(
         grid = aa.Grid2D.uniform(
             shape_native=shape_native,
             pixel_scales=(pixel_scale, pixel_scale),
-            origin=grid.mask.zoom_offset_scaled,
+            origin=zoom.offset_scaled,
         )
 
         grid.is_evaluation_grid = True
@@ -662,9 +664,7 @@ def einstein_mass_angular_list_from(
         einstein_radius_list = self.einstein_radius_list_from(
             grid=grid, pixel_scale=pixel_scale
         )
-        return [
-            np.pi * einstein_radius**2 for einstein_radius in einstein_radius_list
-        ]
+        return [np.pi * einstein_radius**2 for einstein_radius in einstein_radius_list]
 
     @evaluation_grid
     def einstein_mass_angular_from(
diff --git a/autogalaxy/profiles/basis.py b/autogalaxy/profiles/basis.py
index 63747b22c..f71777a00 100644
--- a/autogalaxy/profiles/basis.py
+++ b/autogalaxy/profiles/basis.py
@@ -125,9 +125,11 @@ def image_2d_list_from(
         The image of the light profiles in the basis summed together.
         """
         return [
-            light_profile.image_2d_from(grid=grid, operated_only=operated_only)
-            if not isinstance(light_profile, lp_linear.LightProfileLinear)
-            else np.zeros((grid.shape[0],))
+            (
+                light_profile.image_2d_from(grid=grid, operated_only=operated_only)
+                if not isinstance(light_profile, lp_linear.LightProfileLinear)
+                else np.zeros((grid.shape[0],))
+            )
             for light_profile in self.light_profile_list
         ]
 
diff --git a/autogalaxy/profiles/mass/dark/gnfw_virial_mass_conc.py b/autogalaxy/profiles/mass/dark/gnfw_virial_mass_conc.py
index 9247c9c9d..838b67453 100644
--- a/autogalaxy/profiles/mass/dark/gnfw_virial_mass_conc.py
+++ b/autogalaxy/profiles/mass/dark/gnfw_virial_mass_conc.py
@@ -43,9 +43,7 @@ def kappa_s_and_scale_radius(
 
     ##############################
     def integrand(r):
-        return (r**2 / r**inner_slope) * (1 + r / scale_radius_kpc) ** (
-            inner_slope - 3
-        )
+        return (r**2 / r**inner_slope) * (1 + r / scale_radius_kpc) ** (inner_slope - 3)
 
     de_c = (
         (overdens / 3.0)
diff --git a/autogalaxy/profiles/mass/dark/gnfw_virial_mass_gnfw_conc.py b/autogalaxy/profiles/mass/dark/gnfw_virial_mass_gnfw_conc.py
index ad466df1a..5be2e7840 100644
--- a/autogalaxy/profiles/mass/dark/gnfw_virial_mass_gnfw_conc.py
+++ b/autogalaxy/profiles/mass/dark/gnfw_virial_mass_gnfw_conc.py
@@ -49,9 +49,7 @@ def kappa_s_and_scale_radius(
 
     ##############################
     def integrand(r):
-        return (r**2 / r**inner_slope) * (1 + r / scale_radius_kpc) ** (
-            inner_slope - 3
-        )
+        return (r**2 / r**inner_slope) * (1 + r / scale_radius_kpc) ** (inner_slope - 3)
 
     de_c = (
         (overdens / 3.0)
diff --git a/autogalaxy/profiles/mass/dark/nfw_hk24_util.py b/autogalaxy/profiles/mass/dark/nfw_hk24_util.py
index 0c46b4cfa..71bdce62c 100644
--- a/autogalaxy/profiles/mass/dark/nfw_hk24_util.py
+++ b/autogalaxy/profiles/mass/dark/nfw_hk24_util.py
@@ -3,6 +3,7 @@
 
 @author: felixvecchi
 """
+
 import numpy as np
 from astropy.cosmology import Planck15
 
@@ -206,16 +207,11 @@ def g1_g2_from(x1, x2, e, k_s):
 
     # Prefactors for g1
     pre_f0_g1 = (
-        ((x1 - e) ** 2 + x2**2)
-        * ((x1 + e) ** 2 + x2**2)
-        * (x1**2 - x2**2 - e**2)
+        ((x1 - e) ** 2 + x2**2) * ((x1 + e) ** 2 + x2**2) * (x1**2 - x2**2 - e**2)
     )
 
     pre_f1_g1 = (
-        2
-        * e**2
-        * (x1**2 - 1)
-        * ((x1**2 - x2**2 - e**2) ** 2 - 4 * x1**2 * x2**2)
+        2 * e**2 * (x1**2 - 1) * ((x1**2 - x2**2 - e**2) ** 2 - 4 * x1**2 * x2**2)
         - (x1**2 - x2**2 - e**2) ** 3 * (3 + e**2) / 2
         + 6 * x1**2 * x2**2 * (e**2 - 1) * (x1**2 - x2**2 - e**2)
     )
@@ -225,12 +221,7 @@ def g1_g2_from(x1, x2, e, k_s):
         - 8 * e**2 * x1**2 * x2**2
     )
 
-    pre_f3_g1 = (
-        2
-        * x1
-        * x2
-        * ((x1**2 + x2**2 + e**2) ** 2 - 4 * e**2 * (x2**2 + e**2))
-    )
+    pre_f3_g1 = 2 * x1 * x2 * ((x1**2 + x2**2 + e**2) ** 2 - 4 * e**2 * (x2**2 + e**2))
 
     # First component shear
     g1 = full_pre_factor * (
@@ -245,21 +236,12 @@ def g1_g2_from(x1, x2, e, k_s):
         * x2
         * (
             (x1**2 + x2**2 + e**2)
-            * (
-                (5 * e**2 - 3) * x1**2
-                - 3 * (1 + e**2) * x2**2
-                + (5 - 3 * e**2) * e**2
-            )
+            * ((5 * e**2 - 3) * x1**2 - 3 * (1 + e**2) * x2**2 + (5 - 3 * e**2) * e**2)
             - 4 * e**2 * x1**2 * (1 + e**2)
         )
     )
 
-    pre_f2_g2 = (
-        -2
-        * x1
-        * x2
-        * ((x1**2 + x2**2 + e**2) ** 2 - 4 * e**2 * (x2**2 + e**2))
-    )
+    pre_f2_g2 = -2 * x1 * x2 * ((x1**2 + x2**2 + e**2) ** 2 - 4 * e**2 * (x2**2 + e**2))
 
     pre_f3_g2 = -(
         (x1**2 - x2**2 - e**2) * ((x1**2 + x2**2) ** 2 - e**4)
diff --git a/autogalaxy/quantity/model/plotter_interface.py b/autogalaxy/quantity/model/plotter_interface.py
index 1d925b40c..d15741191 100644
--- a/autogalaxy/quantity/model/plotter_interface.py
+++ b/autogalaxy/quantity/model/plotter_interface.py
@@ -27,12 +27,17 @@ def dataset_quantity(self, dataset: DatasetQuantity):
         dataset
             The imaging dataset which is visualized.
         """
+
+        image_list = [
+            dataset.data.native_for_fits,
+            dataset.noise_map.native_for_fits,
+        ]
+
         hdu_list = hdu_list_for_output_from(
             values_list=[
-                dataset.mask.astype("float"),
-                dataset.data.native,
-                dataset.noise_map.native,
-            ],
+                image_list[0].mask.astype("float"),
+            ]
+            + image_list,
             ext_name_list=[
                 "mask",
                 "data",
diff --git a/autogalaxy/util/shear_field.py b/autogalaxy/util/shear_field.py
index b0f5f02d2..179e8f09a 100644
--- a/autogalaxy/util/shear_field.py
+++ b/autogalaxy/util/shear_field.py
@@ -94,7 +94,6 @@ class ShearYX2D(aa.VectorYX2D, AbstractShearField):
 
 
 class ShearYX2DIrregular(aa.VectorYX2DIrregular, AbstractShearField):
-
     """
     An irregular shear field, which is collection of (y,x) vectors which are located on an irregular grid
     of (y,x) coordinates.
diff --git a/test_autogalaxy/analysis/test_plotter_interface.py b/test_autogalaxy/analysis/test_plotter_interface.py
index ef7e96878..9577f7345 100644
--- a/test_autogalaxy/analysis/test_plotter_interface.py
+++ b/test_autogalaxy/analysis/test_plotter_interface.py
@@ -1,4 +1,6 @@
+import csv
 import shutil
+import numpy as np
 from os import path
 import pytest
 
@@ -35,10 +37,10 @@ def test__galaxies(
     )
 
     image = ag.ndarray_via_fits_from(
-        file_path=path.join(plot_path, "galaxy_images.fits"), hdu=0
+        file_path=path.join(plot_path, "galaxy_images.fits"), hdu=1
     )
 
-    assert image.shape == (7, 7)
+    assert image.shape == (5, 5)
 
 
 def test__inversion(
@@ -59,6 +61,22 @@ def test__inversion(
 
     assert path.join(plot_path, "subplot_inversion_0.png") in plot_patch.paths
 
+    with open(path.join(plot_path, "inversion_reconstruction_0.csv"), mode="r") as file:
+        reader = csv.reader(file)
+        header_list = next(reader)  # ['y', 'x', 'reconstruction', 'noise_map']
+
+        reconstruction_dict = {header: [] for header in header_list}
+
+        for row in reader:
+            for key, value in zip(header_list, row):
+                reconstruction_dict[key].append(float(value))
+
+        # Convert lists to NumPy arrays
+        for key in reconstruction_dict:
+            reconstruction_dict[key] = np.array(reconstruction_dict[key])
+
+    assert reconstruction_dict["x"][0] == pytest.approx(-0.8333333333333334, rel=1.0e-2)
+
 
 def test__adapt_images(
     masked_imaging_7x7,
@@ -83,7 +101,7 @@ def test__adapt_images(
     assert path.join(plot_path, "subplot_adapt_images.png") in plot_patch.paths
 
     image = ag.ndarray_via_fits_from(
-        file_path=path.join(plot_path, "adapt_images.fits"), hdu=0
+        file_path=path.join(plot_path, "adapt_images.fits"), hdu=1
     )
 
-    assert image.shape == (7, 7)
+    assert image.shape == (5, 5)
diff --git a/test_autogalaxy/config/visualize.yaml b/test_autogalaxy/config/visualize.yaml
index 53a1f0700..1c0b5d16c 100644
--- a/test_autogalaxy/config/visualize.yaml
+++ b/test_autogalaxy/config/visualize.yaml
@@ -378,6 +378,7 @@ plots:
     subplot_adapt_images: true
   inversion:
     subplot_inversion: true
+    csv_reconstruction: true          # output reconstruction_mesh.fits containing the reconstructed pixelization and noise map on the source-plane mesh?
   galaxies:
     subplot_galaxies: true
     subplot_galaxy_images: true
diff --git a/test_autogalaxy/imaging/model/test_plotter_interface_imaging.py b/test_autogalaxy/imaging/model/test_plotter_interface_imaging.py
index 0f86cd413..d24358750 100644
--- a/test_autogalaxy/imaging/model/test_plotter_interface_imaging.py
+++ b/test_autogalaxy/imaging/model/test_plotter_interface_imaging.py
@@ -24,7 +24,7 @@ def test__imaging(imaging_7x7, include_2d_all, plot_path, plot_patch):
     assert path.join(plot_path, "subplot_dataset.png") in plot_patch.paths
 
     image = ag.ndarray_via_fits_from(
-        file_path=path.join(plot_path, "dataset.fits"), hdu=0
+        file_path=path.join(plot_path, "dataset.fits"), hdu=1
     )
 
     assert image.shape == (7, 7)
@@ -59,15 +59,15 @@ def test__fit_imaging(
 
     assert path.join(plot_path, "subplot_fit.png") in plot_patch.paths
 
-    image = ag.ndarray_via_fits_from(file_path=path.join(plot_path, "fit.fits"), hdu=0)
+    image = ag.ndarray_via_fits_from(file_path=path.join(plot_path, "fit.fits"), hdu=1)
 
-    assert image.shape == (7, 7)
+    assert image.shape == (5, 5)
 
     image = ag.ndarray_via_fits_from(
-        file_path=path.join(plot_path, "model_galaxy_images.fits"), hdu=0
+        file_path=path.join(plot_path, "model_galaxy_images.fits"), hdu=1
     )
 
-    assert image.shape == (7, 7)
+    assert image.shape == (5, 5)
 
 
 def test__fit_imaging_combined(
diff --git a/test_autogalaxy/interferometer/model/test_plotter_interface_interferometer.py b/test_autogalaxy/interferometer/model/test_plotter_interface_interferometer.py
index 6bf737532..a9f6caca8 100644
--- a/test_autogalaxy/interferometer/model/test_plotter_interface_interferometer.py
+++ b/test_autogalaxy/interferometer/model/test_plotter_interface_interferometer.py
@@ -23,10 +23,10 @@ def test__interferometer(interferometer_7, include_2d_all, plot_path, plot_patch
     assert path.join(plot_path, "subplot_dataset.png") in plot_patch.paths
 
     image = ag.ndarray_via_fits_from(
-        file_path=path.join(plot_path, "dataset.fits"), hdu=0
+        file_path=path.join(plot_path, "dataset.fits"), hdu=1
     )
 
-    assert image.shape == (7, 7)
+    assert image.shape == (7, 2)
 
 
 def test__fit_interferometer(
@@ -45,13 +45,13 @@ def test__fit_interferometer(
     assert path.join(plot_path, "subplot_fit.png") in plot_patch.paths
 
     image = ag.ndarray_via_fits_from(
-        file_path=path.join(plot_path, "model_galaxy_images.fits"), hdu=0
+        file_path=path.join(plot_path, "model_galaxy_images.fits"), hdu=1
     )
 
-    assert image.shape == (7, 7)
+    assert image.shape == (5, 5)
 
     image = ag.ndarray_via_fits_from(
-        file_path=path.join(plot_path, "fit_dirty_images.fits"), hdu=0
+        file_path=path.join(plot_path, "fit_dirty_images.fits"), hdu=1
     )
 
-    assert image.shape == (7, 7)
+    assert image.shape == (5, 5)
diff --git a/test_autogalaxy/quantity/model/test_plotter_interface_quantity.py b/test_autogalaxy/quantity/model/test_plotter_interface_quantity.py
index 96e3c5484..eeb52afc5 100644
--- a/test_autogalaxy/quantity/model/test_plotter_interface_quantity.py
+++ b/test_autogalaxy/quantity/model/test_plotter_interface_quantity.py
@@ -28,7 +28,7 @@ def test__dataset(
     PlotterInterface.dataset_quantity(dataset=dataset_quantity_7x7_array_2d)
 
     image = ag.ndarray_via_fits_from(
-        file_path=path.join(plot_path, "dataset.fits"), hdu=0
+        file_path=path.join(plot_path, "dataset.fits"), hdu=1
     )
 
     assert image.shape == (7, 7)