diff --git a/autogalaxy/imaging/fit_imaging.py b/autogalaxy/imaging/fit_imaging.py
index f06fab55..667a4eaa 100644
--- a/autogalaxy/imaging/fit_imaging.py
+++ b/autogalaxy/imaging/fit_imaging.py
@@ -1,3 +1,19 @@
+"""
+Fit an imaging (CCD) dataset with a model consisting of one or more galaxies.
+
+`FitImaging` is the central object in a **PyAutoGalaxy** imaging analysis. It:
+
+1. Computes the summed image of all non-linear (standard) galaxy light profiles.
+2. Convolves that image with the PSF to create the ``blurred_image``.
+3. Subtracts the blurred image from the data to create the ``profile_subtracted_image``.
+4. If linear light profiles or a pixelization are present, fits the residual image via a linear
+   inversion (see `autogalaxy.galaxy.to_inversion`).
+5. Combines the blurred image and inversion reconstruction into the ``model_data``.
+6. Computes residuals, chi-squared, log-likelihood (or log-evidence when an inversion is used).
+
+When fitting a model via `AnalysisImaging`, the `figure_of_merit` property of `FitImaging` is
+evaluated inside `log_likelihood_function` and returned to the non-linear search.
+"""
 import numpy as np
 from typing import Dict, List, Optional
 
@@ -116,6 +132,13 @@ def profile_subtracted_image(self) -> aa.Array2D:
 
     @property
     def galaxies_to_inversion(self) -> GalaxiesToInversion:
+        """
+        Returns a `GalaxiesToInversion` object that converts the galaxies containing linear light profiles or
+        pixelizations into the inversion objects required for the linear algebra solve.
+
+        The dataset passed to this object is the ``profile_subtracted_image`` — the imaging data with all
+        standard (non-linear) light profile images subtracted — so the inversion only fits the residual signal.
+        """
         dataset = aa.DatasetInterface(
             data=self.profile_subtracted_image,
             noise_map=self.noise_map,
diff --git a/autogalaxy/imaging/model/analysis.py b/autogalaxy/imaging/model/analysis.py
index dcee26d5..1a2f0e46 100644
--- a/autogalaxy/imaging/model/analysis.py
+++ b/autogalaxy/imaging/model/analysis.py
@@ -1,3 +1,18 @@
+"""
+`AnalysisImaging` — the **PyAutoFit** `Analysis` class for fitting galaxy models to CCD imaging data.
+
+This module provides `AnalysisImaging`, which implements `log_likelihood_function` by:
+
+1. Extracting galaxies from the model instance.
+2. Constructing a `FitImaging` object using those galaxies and the stored `Imaging` dataset.
+3. Returning the `figure_of_merit` of the fit (log-likelihood or log-evidence).
+
+It also handles:
+
+- Adapt images: per-galaxy model images from a previous search that drive adaptive pixelizations.
+- Visualization: automatic figure generation during and after the model-fit.
+- Results: wrapping fit outputs into a `ResultImaging` object for downstream use.
+"""
 import numpy as np
 from typing import Optional
 
diff --git a/autogalaxy/imaging/simulator.py b/autogalaxy/imaging/simulator.py
index 68f19cf2..df762288 100644
--- a/autogalaxy/imaging/simulator.py
+++ b/autogalaxy/imaging/simulator.py
@@ -1,3 +1,14 @@
+"""
+Extends the **PyAutoArray** `SimulatorImaging` class with galaxy-aware simulation methods.
+
+`SimulatorImaging` (from `autoarray`) handles the low-level simulation pipeline: applying a PSF, adding
+Poisson noise, and adding a background sky. This module adds a `via_galaxies_from` method that takes a
+list of `Galaxy` objects and a 2D grid, evaluates the galaxy images, and passes them to the parent
+simulation pipeline.
+
+It also handles `LightProfileSNR` objects, automatically scaling each profile's `intensity` so that its
+peak signal-to-noise ratio matches the requested value given the simulation exposure time and sky level.
+"""
 import numpy as np
 from typing import List