diff --git a/CHANGELOG b/CHANGELOG
index 51fc324a0..e8bee4fde 100644
--- a/CHANGELOG
+++ b/CHANGELOG
@@ -14,6 +14,7 @@ What's new
 - Support for user-defined isotopes using the ``Isotope.register()`` method. See the simulation gallery for use cases and examples.
 - Shortcuts for frequency contributions, such as ``Shielding``, ``Isotropic``, and ``cross``.
   Sets of contributions can also be excluded by placing an exclamation mark in front of the string, for example ``"!Shielding"`` excludes shielding interactions.
+- New functions for fitting Czjzek and Extended Czjzek tensor distribution models to experimental spectra. See (TODO: ADD) in the examples gallery for more information.
 - Adds ``DelayEvent`` to the events library.
 - Support for python 3.11
 
@@ -21,6 +22,10 @@ What's new
 
 - New instance method for the `Simulator` class -- `.optimize()` -- which pre-computes transition pathways before least-squares minimization. This improves the efficiency of least-squares minimizations.
 
+
+**Czjzek and Extended Czjzek**
+- The Czjzek model now uses an analytical expression for calculating the probability distribution greatly improving quality and calculation speed.
+
 **Bug Fixes**
 
 - Fixed bug where ``MixingEnum`` class had no attribute ``json`` (Issue [#260](https://github.com/deepanshs/mrsimulator/issues/260))
@@ -31,6 +36,7 @@ What's new
 - Fix bug when calculating ppm scale for large reference offsets.
 
 
+
 v0.7.0
 ------
 
diff --git a/docs/user_guide/spin_system_distributions/czjzek.rst b/docs/user_guide/spin_system_distributions/czjzek.rst
index d2250cf3d..53aaf272d 100644
--- a/docs/user_guide/spin_system_distributions/czjzek.rst
+++ b/docs/user_guide/spin_system_distributions/czjzek.rst
@@ -3,13 +3,21 @@
 Czjzek distribution
 -------------------
 
-The Czjzek distribution models random variations of a second-rank traceless
-symmetric tensors about zero, i.e., a tensor with zeta of zero. See :ref:`czjzek_model`
-for a mathematical description of the model as well as references to examples using the Czjzek
-distribution at the bottom of this page.
+The Czjzek distribution models random variations of second-rank traceless
+symmetric tensors about zero, i.e., a tensor with zeta of zero. An analytical expression
+for the Czjzek distribution exists (cite) which follows
 
-Czjzek distribution of symmetric shielding tensors
-''''''''''''''''''''''''''''''''''''''''''''''''''
+.. math::
+    f(\zeta, \eta, \sigma) = \eta \left(1-\frac{\eta^2}{9}\right)\frac{\zeta^4}{32\sigma^5 \sqrt{2 \pi}} \times \exp\left(-\frac{\zeta^2}{8\sigma^2}\left(1+\frac{\eta^2}{3}\right)\right),
+
+where :math:`\zeta` and :math:`\eta` are the Haberlen components of the tensor and :math:`\sigma` is the Czjzek width parameter. See :ref:`czjzek_model` for a further mathematical description of the model.
+
+The remainder of this page quickly describes how to generate Czjzek distributions and generate
+:py:class:`~mrsimulator.spin_system.SpinSystem` objects from these distributions. Also, look at the
+gallery examples using the Czjzek distribution listed at the bottom of this page.
+
+Creating and sampling a Czjzek distribution
+'''''''''''''''''''''''''''''''''''''''''''
 
 To generate a Czjzek distribution, use the :py:class:`~mrsimulator.models.CzjzekDistribution`
 class as follows.
@@ -21,7 +29,7 @@ class as follows.
 
     cz_model = CzjzekDistribution(sigma=0.8)
 
-The **CzjzekDistribution** class accepts a single argument, ``sigma``, which is the standard
+The **CzjzekDistribution** class accepts the argument, ``sigma``, which is the standard
 deviation of the second-rank traceless symmetric tensor parameters. In the above example,
 we create ``cz_model`` as an instance of the CzjzekDistribution class with
 :math:`\sigma=0.8`.
@@ -36,16 +44,13 @@ function. Let's first draw points from this distribution, using the
 
     zeta_dist, eta_dist = cz_model.rvs(size=50000)
 
-In the above example, we draw *50000* random points of the distribution. The output
-``zeta_dist`` and ``eta_dist`` hold the tensor parameter coordinates of the points, defined
-in the Haeberlen convention.
-The scatter plot of these coordinates is shown below.
+In the above example, we draw *50000* random points of the distribution. The output ``zeta_dist`` and ``eta_dist`` hold the tensor parameter coordinates of the points, defined in the Haeberlen convention. It is further assumed that the points in ``zeta_dist`` are in units of ``ppm`` while ``eta_dist`` has values since :math:`\eta` is dimensionless. The scatter plot of these coordinates is shown below.
 
 .. skip: next
 
 .. plot::
     :context: close-figs
-    :caption: Czjzek Distribution of shielding tensors.
+    :caption: Random sampling of Czjzek distribution of shielding tensors.
 
     import matplotlib.pyplot as plt
 
@@ -57,52 +62,211 @@ The scatter plot of these coordinates is shown below.
     plt.tight_layout()
     plt.show()
 
+Creating and sampling a Czjzek distribution in polar coordinates
+''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
+
+The :py:class:`~mrsimulator.models.czjzek.CzjzekDistribution` class also supports sampling tensors in polar coordinates. The logic behind transforming from a :math:`\zeta`-:math:`\eta` Cartesian grid is further described in mrinversion (cite), and the following definitions are used
+
+.. math::
+
+    \begin{split}r_\zeta = \left| \zeta \right| ~~~~\text{and}~~~~
+    \theta = \left\{ \begin{array}{l r}
+                \frac{\pi}{4} \eta      &: \zeta \le 0, \\
+                \frac{\pi}{2} \left(1 - \frac{\eta}{2} \right) &: \zeta > 0.
+            \end{array}
+            \right.\end{split}
+
+Because Cartesian grids are more manageable in computation, the above polar piece-wise grid is re-express as the x-y Cartesian grid following,
+
+.. math::
+
+    x = r_\zeta \cos\theta ~~~~\text{and}~~~~ y = r_\zeta \sin\theta.
+
+Below, we create another instance of the :py:class:`~mrsimulator.models.czjzek.CzjzekDistribution`
+class with the same value of :math:`sigma=0.8`, but we now also include the argument ``polar=True``
+which means the :py:meth:`~mrsimulator.models.CzjzekDistribution.rvs` will sample x and y points.
+
+.. skip: next
+
+.. plot::
+    :context: close-figs
+    :caption: Random sampling of Czjzek distribution of shielding tensors in polar coordinates.
+
+    cz_model_polar = CzjzekDistribution(sigma=0.8, polar=True)
+
+    # Sample (x, y) points
+    x_dist, y_dist = cz_model_polar.rvs(size=50000)
+
+    # Plot the distribution
+    plt.figure(figsize=(4, 4))
+    plt.scatter(x_dist, y_dist, s=4, alpha=0.02)
+    plt.xlabel("$x$ / ppm")
+    plt.ylabel("$y$ / ppm")
+    plt.xlim(0, 8)
+    plt.ylim(0, 8)
+    plt.tight_layout()
+    plt.show()
+
+
 ----
 
-Czjzek distribution of symmetric quadrupolar tensors
-''''''''''''''''''''''''''''''''''''''''''''''''''''
+Generating probability distribution functions from a Czjzek model
+'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
+
+The :py:meth:`~mrsimulator.models.CzjzekDistribution.pdf` instance method will generate a
+probability distribution function on the supplied grid using the analytical function defined above.
+The provided grid -- passed to the ``pos`` keyword argument -- needs to be defined in either
+Cartesian or polar coordinates depending on if the
+:py:attr:`~mrsimulator.models.CzjzekDistribution.polar` attribute is ``True`` or ``False``.
 
-The Czjzek distribution of symmetric quadrupolar tensors follows a similar setup as the
-Czjzek distribution of symmetric shielding tensors, except we assign the outputs to Cq
-and :math:`\eta_q`. In the following example, we generate the probability distribution
-function using the :py:meth:`~mrsimulator.models.CzjzekDistribution.pdf` method.
+Below, we generate and plot a probability distribution on a :math:`\zeta`-:math:`\eta` Cartesian
+grid where ``zeta_range`` and ``eta_range`` define the desired coordinates in each dimension of the
+grid system.
 
 .. plot::
     :context: close-figs
 
     import numpy as np
 
-    Cq_range = np.arange(100) * 0.3 - 15  # pre-defined Cq range in MHz.
-    eta_range = np.arange(21) / 20  # pre-defined eta range.
-    Cq, eta, amp = cz_model.pdf(pos=[Cq_range, eta_range])
+    cz_model = CzjzekDistribution(sigma=1.2, polar=False)  # sample in (zeta, eta)
 
-To generate a probability distribution, we need to define a grid system over which the
-distribution probabilities will be evaluated. We do so by defining the range of coordinates
-along the two dimensions. In the above example, ``Cq_range`` and ``eta_range`` are the
-range of :math:`\text{Cq}` and :math:`\eta_q` coordinates, which is then given as the
-argument to the :py:meth:`~mrsimulator.models.CzjzekDistribution.pdf` method. The output
-``Cq``, ``eta``, and ``amp`` hold the two coordinates and amplitude, respectively.
+    zeta_range = np.linspace(-12, 12, num=200)  # pre-defined zeta range in units of ppm
+    eta_range = np.linspace(0, 1, num=50)       # pre-defined eta range.
+    zeta_grid, eta_grid, amp = cz_model.pdf(pos=[zeta_range, eta_range])
 
-The plot of the Czjzek probability distribution is shown below.
+Here, ``zeta_grid`` and ``eta_grid`` are numpy arrays defining a set of pair-wise points on the
+grid system, and ``amp`` is another numpy array holding the probability density at each point
+on the grid. Below, the distribution is plotted
 
 .. skip: next
 
 .. plot::
     :context: close-figs
-    :caption: Czjzek Distribution of EFG tensors.
+    :caption: Czjzek Distribution of shielding tensors.
 
-    import matplotlib.pyplot as plt
 
-    plt.contourf(Cq, eta, amp, levels=10)
-    plt.xlabel("$C_q$ / MHz")
+    plt.contourf(zeta_grid, eta_grid, amp, levels=10)
+    plt.xlabel("$\zeta$ / ppm")
     plt.ylabel("$\eta$")
     plt.tight_layout()
     plt.show()
 
-.. note::
-    The ``pdf`` method of the instance generates the probability distribution function
-    by first drawing random points from the distribution and then binning it
-    onto a pre-defined grid.
+---
+
+The probability distribution function can also be generated in polar coordinates. The workflow
+is the same, except we now define an (x, y) grid system using the variables ``x_range``
+and ``y_range``. The code to generate and plot the polar Czjzek distribution is shown below.
+
+.. skip: next
+
+.. plot::
+    :context: close-figs
+    :caption: Czjzek Distribution of shielding tensors in polar coordinates.
+
+    cz_model_polar = CzjzekDistribution(sigma=1.2, polar=True)  # sample in (x, y)
+
+    x_range = np.linspace(0, 10, num=150)
+    y_range = np.linspace(0, 10, num=150)
+    x_grid, y_grid, amp = cz_model_polar.pdf(pos=[x_range, y_range])
+
+    plt.figure(figsize=(4, 4))
+    plt.contourf(x_grid, y_grid, amp, levels=10)
+    plt.xlabel("$x$ / ppm")
+    plt.ylabel("$y$ / ppm")
+    plt.tight_layout()
+    plt.show()
+
+
+Distributions of shielding and quadrupolar tensors and a note on units
+''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
+
+The :py:class:`~mrsimulator.models.CzjzekDistribution` class can be used to generate
+distributions for both symmetric chemical shielding tensors and electric field gradient
+tensors. It is important to note the Czjzek model is only aware of the Haberlen components
+of the tensors and not the units of the tensor. In the above code cells, we generated
+distributions for symmetric shielding tensors and assumed all units for :math:`\zeta` were
+in ppm.
+
+Quadrupolar tensors, defined using values of :math:`C_q` in MHz and unitless :math:`\eta`,
+can also be drawn from the Czjzek distribution in the same manner; however, the dimensions
+are assumed to be in units of MHz. The following code draws a distribution of quadrupolar
+tensor parameters.
+
+.. skip: next
+
+.. plot::
+    :context: close-figs
+
+    Cq_range = np.linspace(-12, 12, num=200)  # pre-defined Cq range in units of MHz
+    eta_range = np.linspace(0, 1, num=50)     # pre-defined eta range.
+    Cq_grid, eta_grid, amp = cz_model.pdf(pos=[Cq_range, eta_range])
+
+
+the units for ``Cq_range`` and ``Cq_grid`` are assumed in MHz. Similarly, x and y are assumed to
+be in units of MHz when sampling quadrupolar tensors in polar coordinates.
+
+.. skip: next
+
+.. plot::
+    :context: close-figs
+
+    x_range = np.linspace(0, 10, num=150)  # pre-defined x grid in units of MHz
+    y_range = np.linspace(0, 10, num=150)  # pre-defined y grid in units of MHz
+    x_grid, y_grid, amp = cz_model_polar.pdf(pos=[x_range, y_range])
+
+Generating a list of SpinSystem objects from a Czjzek model
+'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
+
+The utility function :py:meth:`~mrsimulator.utils.collection.single_site_system_generator`, further
+described in :ref:`single_site_system_generator_documentation`, can be used in conjunction with
+the :py:class:`~mrsimulator.models.CzjzekDistribution` class to generate a set of spin systems whose
+tensor parameters follow the Czjzek distribution.
+
+.. plot::
+    :context: close-figs
+
+    from mrsimulator.utils.collection import single_site_system_generator
+
+    # Distribution of quadrupolar tensors
+    cz_model = CzjzekDistribution(sigma=0.7)
+    Cq_range = np.linspace(0, 10, num=100)
+    eta_range = np.linspace(0, 1, num=50)
+
+    # Create (Cq, eta) grid points and amplitude
+    Cq_grid, eta_grid = np.meshgrid(Cq_range, eta_range)
+    _, _, amp = cz_model.pdf(pos=[Cq_range, eta_range])
+
+    sys = single_site_system_generator(
+        isotope="27Al",
+        quadrupolar={"Cq": Cq_grid * 1e6, "eta": eta_grid},  # Cq argument in units of Hz
+        abundance=amp,
+    )
+
+A spin system will be generated for each point on the :math:`\zeta`-:math:`\eta` grid, and the
+abundance of each spin system matches the amplitude of the Czjzek distribution. When working in
+polar coordinates, the set of :math:`\left(x, y\right)` coordinates needs to be transformed into
+a set of :math:`\left(\zeta, \eta\right)` coordinates before being passed to the
+:py:meth:`~mrsimulator.utils.collection.single_site_system_generator` function. The utility
+function :py:meth:`~mrsimulator.utils.x_y_to_zeta_eta` performs this transformation, as shown
+below.
+
+.. plot::
+    :context: close-figs
+
+    from mrsimulator.models.utils import x_y_to_zeta_eta
+
+    # Sample distribution of shielding tensors in polar coords
+    cz_model_polar = CzjzekDistribution(sigma=0.7, polar=True)
+    x_range = np.linspace(0, 10, num=150)
+    y_range = np.linspace(0, 10, num=150)
+
+    # Create (x, y) grid points and amplitude
+    x_grid, y_grid, amp = cz_model_polar.pdf(pos=[x_range, y_range])
+
+    # To transformation (x, y) -> (zeta, eta)
+    zeta_grid, eta_grid = x_y_to_zeta_eta(x_grid, y_grid)
+
+---
 
 .. minigallery:: mrsimulator.models.CzjzekDistribution
     :add-heading: Mini-gallery using the Czjzek distributions
diff --git a/docs/user_guide/spin_system_distributions/extended_czjzek.rst b/docs/user_guide/spin_system_distributions/extended_czjzek.rst
index 6b5b81c13..5af596864 100644
--- a/docs/user_guide/spin_system_distributions/extended_czjzek.rst
+++ b/docs/user_guide/spin_system_distributions/extended_czjzek.rst
@@ -3,9 +3,7 @@
 Extended Czjzek distribution
 ----------------------------
 
-The Extended Czjzek distribution models random variations of a second-rank traceless
-symmetric tensors about a non-zero tensor.  See :ref:`ext_czjzek_model` and
-references within for a brief description of the model.
+The Extended Czjzek distribution models random variations of second-rank traceless symmetric tensors about a non-zero tensor. Unlike the Czjzek distribution, the Extended Czjzek model has no known analytical function for the probability distribution. Therefore, mrsimulator relies on random sampling to approximate the probability distribution function. See :ref:`ext_czjzek_model` and references within for a further description of the model.
 
 Extended Czjzek distribution of symmetric shielding tensors
 '''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
@@ -76,20 +74,24 @@ function using the :py:meth:`~mrsimulator.models.ExtCzjzekDistribution.pdf` meth
 
     import numpy as np
 
-    Cq_range = np.arange(100) * 0.04 + 2  # pre-defined Cq range in MHz.
-    eta_range = np.arange(21) / 20  # pre-defined eta range.
+    Cq_range = np.linspace(2, 6, num=100)  # pre-defined Cq range in MHz.
+    eta_range = np.linspace(0, 1, num=20)  # pre-defined eta range.
 
     quad_tensor = {"Cq": 3.5, "eta": 0.23}  # Cq assumed in MHz
     model_quad = ExtCzjzekDistribution(quad_tensor, eps=0.2)
-    Cq, eta, amp = model_quad.pdf(pos=[Cq_range, eta_range])
+    Cq_grid, eta_grid, amp = model_quad.pdf(pos=[Cq_range, eta_range], size=400000)
 
-As with the case with Czjzek distribution, to generate a probability distribution of the
+As with the case of the Czjzek distribution, to generate a probability distribution of the
 extended Czjzek distribution, we need to define a grid system over which the distribution
 probabilities will be evaluated. We do so by defining the range of coordinates along the
 two dimensions. In the above example, ``Cq_range`` and ``eta_range`` are the
 range of :math:`\text{Cq}` and :math:`\eta_q` coordinates, which is then given as the
-argument to the :py:meth:`~mrsimulator.models.ExtCzjzekDistribution.pdf` method. The output
-``Cq``, ``eta``, and ``amp`` hold the two coordinates and amplitude, respectively.
+argument to the :py:meth:`~mrsimulator.models.ExtCzjzekDistribution.pdf` method. The pdf
+method also accepts the keyword argument ``size`` which defines the number of random samples
+used to approximate the probability distribution. A larger number will create better
+approximations, although this increased quality comes at the expense of computation time.
+The output ``Cq_grid``, ``eta_grid``, and ``amp`` hold the two coordinates and
+amplitude, respectively.
 
 The plot of the extended Czjzek probability distribution is shown below.
 
@@ -101,12 +103,55 @@ The plot of the extended Czjzek probability distribution is shown below.
 
     import matplotlib.pyplot as plt
 
-    plt.contourf(Cq, eta, amp, levels=10)
+    plt.contourf(Cq_grid, eta_grid, amp, levels=10)
     plt.xlabel("$C_q$ / MHz")
     plt.ylabel("$\eta$")
     plt.tight_layout()
     plt.show()
 
+Extended Czjzek distribution in polar coordinates
+'''''''''''''''''''''''''''''''''''''''''''''''''
+
+As with the Czjzek distribution, we can sample an Extended Czjzek distribution on a polar
+(x, y) grid. Below, we construct two equivalent
+:py:class:`~mrsimulator.models.ExtCzjzekDistribution` objects, except one is defined in polar
+coordinates.
+
+.. skip: next
+
+.. plot::
+    :context: close-figs
+    :caption: Two equivalent Extended Czjzek distributions in Cartesian :math:`\left(\zeta, \eta\right)` coordinates (left) and in polar :math:`\left(x, y\right)` coordinates (right).
+
+    quad_tensor = {"Cq": 4.2, "eta": 0.15}  # Cq assumed in MHz
+    ext_cz_model = ExtCzjzekDistribution(quad_tensor, eps=0.4)
+    ext_cz_model_polar = ExtCzjzekDistribution(quad_tensor, eps=0.4, polar=True)
+
+    # Distribution in cartesian (zeta, eta) coordinates
+    Cq_range = np.linspace(2, 8, num=50)
+    eta_range = np.linspace(0, 1, num=20)
+    Cq_grid, eta_grid, amp = ext_cz_model.pdf(pos=[Cq_range, eta_range], size=2000000)
+
+    # Distribution in polar coordinates
+    x_range = np.linspace(0, 6, num=36)
+    y_range = np.linspace(0, 6, num=36)
+    x_grid, y_grid, amp_polar = ext_cz_model_polar.pdf(pos=[x_range, y_range], size=2000000)
+
+    # Plot the distributions
+    fig, ax = plt.subplots(1, 2, figsize=(9, 4), gridspec_kw={"width_ratios": (5, 4)})
+    ax[0].contourf(Cq_grid, eta_grid, amp, levels=10)
+    ax[0].set_xlabel("$C_q$ / MHz")
+    ax[0].set_ylabel("$\eta$")
+    ax[0].set_title("Cartesian coordinates")
+    ax[1].contourf(x_grid, y_grid, amp_polar, levels=10)
+    ax[1].set_xlabel("x / MHz")
+    ax[1].set_ylabel("y / MHz")
+    ax[1].set_title("Polar coordinates")
+
+    plt.tight_layout()
+    plt.show()
+
+
 .. note::
     The ``pdf`` method of the instance generates the probability distribution function
     by first drawing random points from the distribution and then binning it
diff --git a/examples_source/1D_simulation(macro_amorphous)/plot_5_czjzek_distribution.py b/examples_source/1D_simulation(macro_amorphous)/plot_5_czjzek_distribution.py
index 16588cfc3..6d17e9914 100644
--- a/examples_source/1D_simulation(macro_amorphous)/plot_5_czjzek_distribution.py
+++ b/examples_source/1D_simulation(macro_amorphous)/plot_5_czjzek_distribution.py
@@ -33,7 +33,8 @@
 # The range of zeta and eta coordinates over which the distribution is sampled.
 z_range = np.arange(100) - 50  # in ppm
 e_range = np.arange(21) / 20
-z_dist, e_dist, amp = CzjzekDistribution(sigma=3.1415).pdf(pos=[z_range, e_range])
+z_dist, e_dist = np.meshgrid(z_range, e_range)
+_, _, amp = CzjzekDistribution(sigma=3.1415).pdf(pos=[z_range, e_range])
 
 # %%
 # Here ``z_range`` and ``e_range`` are the coordinates along the :math:`\zeta` and
@@ -94,7 +95,8 @@
 # The range of Cq and eta coordinates over which the distribution is sampled.
 cq_range = np.arange(100) * 0.6 - 30  # in MHz
 e_range = np.arange(21) / 20
-cq_dist, e_dist, amp = CzjzekDistribution(sigma=2.3).pdf(pos=[cq_range, e_range])
+cq_dist, e_dist = np.meshgrid(cq_range, e_range)
+_, _, amp = CzjzekDistribution(sigma=2.3).pdf(pos=[cq_range, e_range])
 
 # The following is the contour plot of the Czjzek distribution.
 plt.figure(figsize=(4.25, 3.0))
diff --git a/examples_source/1D_simulation(macro_amorphous)/plot_6_extended_czjzek.py b/examples_source/1D_simulation(macro_amorphous)/plot_6_extended_czjzek.py
index e31693e33..e05c6d457 100644
--- a/examples_source/1D_simulation(macro_amorphous)/plot_6_extended_czjzek.py
+++ b/examples_source/1D_simulation(macro_amorphous)/plot_6_extended_czjzek.py
@@ -36,7 +36,8 @@
 e_lim = np.arange(21) / 20
 
 dominant = {"zeta": 60, "eta": 0.3}
-z_dist, e_dist, amp = ExtCzjzekDistribution(dominant, eps=0.14).pdf(pos=[z_lim, e_lim])
+z_dist, e_dist = np.meshgrid(z_lim, e_lim)
+_, _, amp = ExtCzjzekDistribution(dominant, eps=0.14).pdf(pos=[z_lim, e_lim])
 
 # %%
 # The following is the plot of the extended Czjzek distribution.
@@ -93,9 +94,8 @@
 e_lim = np.arange(21) / 20
 
 dominant = {"Cq": 6.1, "eta": 0.1}
-cq_dist, e_dist, amp = ExtCzjzekDistribution(dominant, eps=0.25).pdf(
-    pos=[cq_lim, e_lim]
-)
+cq_dist, e_dist = np.meshgrid(cq_lim, e_lim)
+_, _, amp = ExtCzjzekDistribution(dominant, eps=0.25).pdf(pos=[cq_lim, e_lim])
 
 # %%
 # The following is the plot of the extended Czjzek distribution.
diff --git a/fitting_source/1D_fitting/plot_6_139La_Ext_Czjzek.py b/fitting_source/1D_fitting/plot_6_139La_Ext_Czjzek.py
index 0ba9d73f4..8f40c8bfb 100644
--- a/fitting_source/1D_fitting/plot_6_139La_Ext_Czjzek.py
+++ b/fitting_source/1D_fitting/plot_6_139La_Ext_Czjzek.py
@@ -22,7 +22,7 @@
 import matplotlib.pyplot as plt
 
 
-# sphinx_gallery_thumbnail_number = 3
+# sphinx_gallery_thumbnail_number = 4
 
 # %%
 # Import the dataset
@@ -33,7 +33,7 @@
 experiment.x[0].to("ppm", "nmr_frequency_ratio")
 experiment /= experiment.max()
 
-plt.figure(figsize=(6, 4))
+plt.figure(figsize=(4, 3))
 ax = plt.subplot(projection="csdm")
 ax.plot(experiment, "k", alpha=0.5)
 plt.grid()
@@ -91,10 +91,11 @@ def make_kernel(pos, method, isotope):
     Cq, eta = np.meshgrid(pos[0], pos[1], indexing="xy")
 
     spin_systems = [
-        SpinSystem(sites=[Site(isotope=isotope, quadrupolar=dict(Cq=cq * 1e6, eta=e))])
-        for cq, e in zip(Cq.ravel(), eta.ravel())
+        SpinSystem(sites=[Site(isotope=isotope, quadrupolar=dict(Cq=cq_, eta=e_))])
+        for cq_, e_ in zip(Cq.ravel(), eta.ravel())
     ]
     sim = Simulator(spin_systems=spin_systems, methods=[method])
+    sim.config.number_of_sidebands = 4
     sim.config.decompose_spectrum = "spin_system"
     sim.run(pack_as_csdm=False)  # Will return spectrum as numpy array, not CSDM object
 
@@ -103,7 +104,7 @@ def make_kernel(pos, method, isotope):
 
 
 # Create ranges to construct cq and eta grid points
-cq_range = np.linspace(0, 100, num=100) * 0.8 + 25
+cq_range = (np.linspace(0, 100, num=100) * 0.8 + 25) * 1e6  # in Hz
 eta_range = np.arange(21) / 20
 pos = [cq_range, eta_range]
 kernel = make_kernel(pos, method, "139La")
@@ -148,13 +149,13 @@ def make_spectrum_from_parameters(params, kernel, processor, pos, distribution):
     iso_shift = values["dist_iso_shift"]
 
     # Setup model object and get the amplitude
-    distribution.symmetric_tensor = {"Cq": Cq, "eta": eta}
+    distribution.symmetric_tensor.Cq = Cq
+    distribution.symmetric_tensor.eta = eta
     distribution.eps = eps
-    # model_quad = ExtCzjzekDistribution({"Cq": Cq, "eta": eta}, eps)
     _, _, amp = distribution.pdf(pos=pos)
 
     # Create spectra by dotting the amplitude distribution with the kernel
-    dist = kernel.dot(amp.ravel())
+    dist = np.dot(kernel, amp.ravel())
 
     # Pack numpy array as csdm object and apply signal processing
     guess_dataset = cp.CSDM(
@@ -197,10 +198,10 @@ def residuals(exp_spectra, simulated_spectra):
 #     If you adapt this example to your own dataset, make sure the initial guess is
 #     decently good, otherwise LMFIT is likely to fall into a local minima.
 params = lmfit.Parameters()
-params.add("dist_Cq", value=49.5)
+params.add("dist_Cq", value=49.5e6)  # Hz
 params.add("dist_eta", value=0.55, min=0, max=1)
-params.add("dist_eps", value=0.24, min=0)
-params.add("dist_iso_shift", value=350)
+params.add("dist_eps", value=0.1, min=0)
+params.add("dist_iso_shift", value=350)  # ppm
 params.add("sp_scale_factor", value=3.8e3, min=0)
 
 # Plot the initial guess spectrum along with the experimental data
@@ -208,18 +209,32 @@ def residuals(exp_spectra, simulated_spectra):
     symmetric_tensor={"Cq": params["dist_Cq"], "eta": params["dist_eta"]},
     eps=params["dist_eps"],
 )
-initial_guess = make_spectrum_from_parameters(
+guess_distribution = distribution.pdf(pos, size=400_000, pack_as_csdm=True)
+
+initial_guess_spectrum = make_spectrum_from_parameters(
     params, kernel, processor, pos, distribution
 )
-residual_spectrum = residuals(experiment, initial_guess)
+residual_spectrum = residuals(experiment, initial_guess_spectrum)
 
-plt.figure(figsize=(6, 4))
+plt.figure(figsize=(4, 3))
 ax = plt.subplot(projection="csdm")
-ax.plot(experiment.real, "k", alpha=0.5)
-ax.plot(initial_guess.real, "r", alpha=0.75, label="Guess")
-ax.plot(residual_spectrum.real, "b", alpha=0.75, label="Residuals")
+ax.plot(experiment.real, "k", alpha=0.5, label="Experiment")
+ax.plot(initial_guess_spectrum.real, "r", alpha=0.3, label="Guess")
+ax.plot(residual_spectrum.real, "b", alpha=0.3, label="Residuals")
 plt.legend()
 plt.grid()
+plt.title("Initial Guess")
+plt.tight_layout()
+plt.show()
+
+
+# %%
+# Guess distribution
+plt.figure(figsize=(4, 3))
+ax = plt.subplot(projection="csdm")
+ax.imshow(guess_distribution, interpolation="none", cmap="gist_ncar_r", aspect="auto")
+ax.set_xlabel("Cq / Hz")
+ax.set_ylabel(r"$\eta$")
 plt.tight_layout()
 plt.show()
 
@@ -242,17 +257,28 @@ def minimization_function(
 
 
 # %%
-# Create the Minimization object with all the functions and variables previously
-# created.
+# Sinice the probabilty distribution is generated from a sparsely sampled
+# from a 5D second rank tensor parameter space, we increase the `diff_step` size from
+# machine precession to avoid approaching local minima from noise.
+scipy_minimization_kwargs = dict(
+    diff_step=1e-4,  # Increase step size from machine precession
+    gtol=1e-10,  # Decrease global convergence requirement (default 1e-8)
+    xtol=1e-10,  # Decrease variable convergence requirement (default 1e-8)
+    verbose=2,  # Print minimization info during each step
+    loss="linear",
+)
+
 minner = Minimizer(
     minimization_function,
     params,
     fcn_args=(experiment, processor, kernel, pos, distribution),
+    **scipy_minimization_kwargs,
 )
-result = minner.minimize(method="powell")
+result = minner.minimize(method="least_squares")
 best_fit_params = result.params  # Grab the Parameters object from the best fit
 result
 
+
 # %%
 # Plot the best-fit solution
 # --------------------------
@@ -263,17 +289,29 @@ def minimization_function(
 )
 residual_spectrum = residuals(experiment, final_fit)
 
-plt.figure(figsize=(6, 4))
+plt.figure(figsize=(4, 3))
 ax = plt.subplot(projection="csdm")
-ax.plot(experiment, "k", alpha=0.5)
-ax.plot(final_fit, "r", alpha=0.75, label="Fit")
-ax.plot(residual_spectrum, "b", alpha=0.75, label="Residuals")
+ax.plot(experiment, "k", alpha=0.5, label="Experiment")
+ax.plot(final_fit, "r", alpha=0.3, label="Fit")
+ax.plot(residual_spectrum, "b", alpha=0.3, label="Residuals")
 plt.legend()
 ax.set_xlim(-11000, 9000)
 plt.grid()
+plt.title("Best Fit")
 plt.tight_layout()
 plt.show()
 
+
+# %%
+# Best fit probability distribution
+prob = distribution.pdf(pos, pack_as_csdm=True)
+plt.figure(figsize=(4, 3))
+ax = plt.subplot(projection="csdm")
+ax.imshow(prob, interpolation="none", cmap="gist_ncar_r", aspect="auto")
+ax.set_xlabel("Cq / Hz")
+ax.set_ylabel(r"$\eta$")
+plt.tight_layout()
+plt.show()
 # %%
 #
 # .. [#f1] A. J. Fernández-Carrión, M. Allix, P. Florian, M. R. Suchomel, and A. I.
diff --git a/fitting_source/1D_fitting/plot_7_27Al_Extended_Czjzek.py b/fitting_source/1D_fitting/plot_7_27Al_Extended_Czjzek.py
new file mode 100644
index 000000000..5aa22de4f
--- /dev/null
+++ b/fitting_source/1D_fitting/plot_7_27Al_Extended_Czjzek.py
@@ -0,0 +1,214 @@
+#!/usr/bin/env python
+"""
+Fitting a Czjzek Model
+^^^^^^^^^^^^^^^^^^^^^^
+"""
+# %%
+# In this example, we illustrate an application of mrsimulator in fitting a lineshape
+# from a Czjzek distribution of quadrupolar tensors to an experimental
+# :math:`^{27}\text{Al}` MAS spectrum of a phospho-aluminosilicate glass. Setting up
+# the least-squares minimization for a distribution of spin systems is slightly
+# different than that of a crystalline solid.
+#
+# There are 4 steps involved in the processs:
+#  - Importing the experimental dataset,
+#  - Generating a pre-computed line shape kernel of subspectra,
+#  - Creating parameters for the Czjzek distribution model from an initial guess,
+#  - Minimizing and visualizing.
+#
+import numpy as np
+import csdmpy as cp
+import matplotlib.pyplot as plt
+import mrsimulator.signal_processor as sp
+from mrsimulator.simulator.config import ConfigSimulator
+import mrsimulator.utils.spectral_fitting as sf
+from lmfit import Minimizer
+
+from mrsimulator.method.lib import BlochDecayCTSpectrum
+from mrsimulator.utils import get_spectral_dimensions
+
+from mrsimulator.models.czjzek import CzjzekDistribution
+from mrsimulator.models.utils import LineShapeKernel
+
+# sphinx_gallery_thumbnail_number = 3
+
+# %%
+# Import the experimental dataset
+# -------------------------------
+#
+# Below we import and visualize the experimental dataset.
+host = "http://ssnmr.org/sites/default/files/mrsimulator/"
+filename = "20K_20Al_10P_50Si_HahnEcho_27Al.csdf"
+exp_data = cp.load(host + filename).real
+
+exp_data.x[0].to("ppm", "nmr_frequency_ratio")
+exp_data /= exp_data.max()
+
+plt.figure(figsize=(4, 3))
+ax = plt.subplot(projection="csdm")
+ax.plot(exp_data)
+plt.tight_layout()
+plt.show()
+
+# %%
+# Generating a line shape kernel
+# ------------------------------
+#
+# The Czjzek distribution is a statistical model that represents a five-dimensional
+# multivariate normal distribution of second-rank tensor components. These random
+# second-rank tensors are converted into corresponding anisotropy and asymmetry
+# parameters using the Haeberlen convention. The resulting parameters are then binned
+# on a two-dimensional grid, forming the Czjzek probability distribution. The Czjzek
+# spectra is which are the probability weighted sum of lineshapes at each grid point.
+#
+# However, simulating the spectra of spin systems at each minimization step for the
+# least-squares fitting is computationally expensive. A more efficient way is to
+# pre-define a grid system for the tensor parameters, simulate a library of sub-spectra
+# for each grid point, and only update the probability distribution during each
+# minimization step.
+#
+# To simulate the spectra of the given experiment, we create a Method object. Using
+# this method, we generate a kernel of lineshape sub-spectra defined on a polar grid by
+# using the LineShapeKernel class. The argument "method" is the mrsimulator method
+# object used for the simulation, "pos" is a tuple of coordinates defining the
+# two-dimensional grid, "polar" defines a polar parameter coordinate space, and
+# "config" is the Simulator config object used in lineshape simulation. Here, we use
+# the "config" argument to set the number of sidebands as 8.
+#
+# The kernel is generated with the "generate_lineshape()" function, where we pass the
+# "tensor_type='quadrupolar'" argument to specify a quadrupolar parameter grid. A
+# symmetric shielding lineshape kernel can also be generated by specifying
+
+# Create a Method object to simulate the spectrum
+spectral_dimensions = get_spectral_dimensions(exp_data)
+method = BlochDecayCTSpectrum(
+    channels=["27Al"],
+    rotor_frequency=14.2e3,
+    spectral_dimensions=spectral_dimensions,
+    experiment=exp_data,
+)
+
+# Define a polar grid for the lineshape kernel
+x = np.linspace(0, 2e7, num=36)
+y = np.linspace(0, 2e7, num=36)
+pos = (x, y)
+
+# Generate the kernel
+sim_config = ConfigSimulator(number_of_sidebands=8)
+quad_kernel = LineShapeKernel(method=method, pos=pos, polar=True, config=sim_config)
+quad_kernel.generate_lineshape(tensor_type="quadrupolar")
+
+print("Desired Kernel shape: ", (x.size * y.size, spectral_dimensions[0]["count"]))
+print("Actual Kernel shape:  ", quad_kernel.kernel.shape)
+
+# %%
+# Create a Parameters object
+# --------------------------
+# Next, we create an instance of the `CzjzekDistribution` class with initial guess
+# values along with a `SignalProcessor` object.
+
+# Create initial guess CzjzekDistribution
+cz_model = CzjzekDistribution(
+    mean_isotropic_chemical_shift=60.0, sigma=1.4e6, polar=True
+)
+all_models = [cz_model]
+
+processor = sp.SignalProcessor(
+    operations=[
+        sp.IFFT(),
+        sp.apodization.Gaussian(FWHM="600 Hz"),
+        sp.FFT(),
+        sp.Scale(factor=30),
+    ]
+)
+
+# %%
+# Make the Parameters object and simulate a guess spectrum.
+# Note that the variable `sf_kwargs` holds some additional keyword arguments that
+# many of the spectral fitting function takes in. This dictionary needs to be updated to
+# reflect any changes made in the minimization.
+params = sf.make_LMFIT_params(spin_system_models=all_models, processors=[processor])
+
+# Additional keyword arguments passed to best-fit and residual functions.
+sf_kwargs = dict(
+    kernel=quad_kernel,
+    spin_system_models=all_models,
+    processor=processor,
+)
+
+# Make a guess and residuals spectrum from the initial guess
+guess = sf.bestfit_dist(params=params, **sf_kwargs)
+residuals = sf.residuals_dist(params=params, **sf_kwargs)
+
+plt.figure(figsize=(4, 3))
+ax = plt.subplot(projection="csdm")
+ax.plot(exp_data, "k", alpha=0.5, label="Experiment")
+ax.plot(guess, "r", alpha=0.3, label="Guess")
+ax.plot(residuals, "b", alpha=0.3, label="Residuals")
+plt.legend()
+plt.grid()
+plt.title("Initial Guess")
+plt.tight_layout()
+plt.show()
+
+# Print the Parameters object
+params
+
+# %%
+# Create and run a minimization
+# -----------------------------
+# Finally, a `Minimizer` object is created and a minimization run using least-squares.
+# The same arguments defined in the `addtl_sf_kwargs` variable are also passed to the
+# minimizer. Sinice the probabilty distribution is generated from a sparsely sampled
+# from a 5D second rank tensor parameter space, we increase the `diff_step` size from
+# machine precession to avoid approaching local minima from noise.
+scipy_minimization_kwargs = dict(
+    diff_step=1e-4,  # Increase step size from machine precesion.
+    gtol=1e-10,  # Decrease global convergence requirement (default 1e-8)
+    xtol=1e-10,  # Decrease variable convergence requirement (default 1e-8)
+    verbose=2,  # Print minimization info during each step
+    loss="linear",
+)
+
+minner = Minimizer(
+    sf.LMFIT_min_function_dist,
+    params,
+    fcn_kws=sf_kwargs,
+    **scipy_minimization_kwargs,
+)
+result = minner.minimize(method="least_squares")
+result
+
+# %%
+# Plot the best-fit spectrum
+# --------------------------
+bestfit = sf.bestfit_dist(params=result.params, **sf_kwargs)
+residuals = sf.residuals_dist(params=result.params, **sf_kwargs)
+
+plt.figure(figsize=(4, 3))
+ax = plt.subplot(projection="csdm")
+ax.plot(exp_data, "k", alpha=0.5, label="Experiment")
+ax.plot(bestfit, "r", alpha=0.3, label="Fit")
+ax.plot(residuals, "b", alpha=0.3, label="Residuals")
+plt.legend()
+plt.grid()
+plt.title("Best Fit")
+plt.tight_layout()
+plt.show()
+
+# %%
+# Plot the best-fit distribution
+# ------------------------------
+#
+for i, model in enumerate(all_models):
+    model.update_lmfit_params(result.params, i)
+
+amp = cz_model.pdf(pos=pos, pack_as_csdm=True)
+
+plt.figure(figsize=(4, 3))
+ax = plt.subplot(projection="csdm")
+ax.imshow(amp, cmap="gist_ncar_r", interpolation="none", aspect="auto")
+ax.set_xlabel("x / Hz")
+ax.set_ylabel("y / Hz")
+plt.tight_layout()
+plt.show()
diff --git a/src/mrsimulator/models/analytical_distributions.py b/src/mrsimulator/models/analytical_distributions.py
index f36e42e7d..b1084064e 100644
--- a/src/mrsimulator/models/analytical_distributions.py
+++ b/src/mrsimulator/models/analytical_distributions.py
@@ -1,5 +1,5 @@
 import numpy as np
-from mrsimulator.models.utils import x_y_from_zeta_eta
+from mrsimulator.models.utils import zeta_eta_to_x_y
 
 
 __author__ = "Deepansh J. Srivastava"
@@ -52,7 +52,7 @@ def czjzek_zeta_eta(sigma: float, pos: list):
     pdf_model = pdf_model.ravel()
     pdf_model[eta_idx] /= 2.0
     pdf_model.shape = zeta.shape
-    return zeta, eta, pdf_model
+    return pos[0], pos[1], pdf_model
 
 
 def czjzek_polar(sigma: float, pos: list):
@@ -73,7 +73,7 @@ def czjzek_polar(sigma: float, pos: list):
     pdf_model = czjzek_distribution(sigma, zeta, eta)
     eta = eta.ravel()
     zeta = zeta.ravel()
-    x, y = x_y_from_zeta_eta(zeta, eta)
+    x, y = zeta_eta_to_x_y(zeta, eta)
 
     eta_idx = np.where(eta == 1)
     pdf_model = pdf_model.ravel()
@@ -89,4 +89,4 @@ def czjzek_polar(sigma: float, pos: list):
     )
     hist_x_y += hist_x_y.T
     hist_x_y /= hist_x_y.sum()
-    return x, y, hist_x_y
+    return pos[0], pos[1], hist_x_y
diff --git a/src/mrsimulator/models/czjzek.py b/src/mrsimulator/models/czjzek.py
index 80883694b..dfca75b14 100644
--- a/src/mrsimulator/models/czjzek.py
+++ b/src/mrsimulator/models/czjzek.py
@@ -4,14 +4,14 @@
 __author__ = "Deepansh J. Srivastava"
 __email__ = "dsrivastava@hyperfine.io"
 """
+import csdmpy as cp
 import mrsimulator.models.analytical_distributions as analytical_dist
 import numpy as np
 from mrsimulator.spin_system.tensors import SymmetricTensor
 
 from .utils import get_Haeberlen_components
 from .utils import get_principal_components
-from .utils import x_y_from_zeta_eta
-
+from .utils import zeta_eta_to_x_y
 
 __author__ = "Deepansh J. Srivastava"
 __email__ = "srivastava.89@osu.edu"
@@ -60,15 +60,15 @@ def _czjzek_random_distribution_tensors(sigma, n):
     sqrt_3_U5 = np.random.normal(0.0, sqrt_3_sigma, n)
 
     # Create N random tensors
-    tensors = np.zeros((n, 3, 3))  # n x 3 x 3 tensors
+    tensors = np.empty((n, 3, 3))  # n x 3 x 3 tensors
 
     tensors[:, 0, 0] = sqrt_3_U5 - U1  # xx
-    tensors[:, 0, 1] = sqrt_3_U4  # xy
-    tensors[:, 0, 2] = sqrt_3_U2  # xz
+    # tensors[:, 0, 1] = sqrt_3_U4  # xy
+    # tensors[:, 0, 2] = sqrt_3_U2  # xz
 
     tensors[:, 1, 0] = sqrt_3_U4  # yx
     tensors[:, 1, 1] = -sqrt_3_U5 - U1  # yy
-    tensors[:, 1, 2] = sqrt_3_U3  # yz
+    # tensors[:, 1, 2] = sqrt_3_U3  # yz
 
     tensors[:, 2, 0] = sqrt_3_U2  # zx
     tensors[:, 2, 1] = sqrt_3_U3  # zy
@@ -78,7 +78,31 @@ def _czjzek_random_distribution_tensors(sigma, n):
 
 
 class AbstractDistribution:
-    def pdf(self, pos, size: int = 400000, analytical: bool = True):
+    """Abstract distribution"""
+
+    model_name = "base"
+
+    def __init__(
+        self,
+        mean_isotropic_chemical_shift=0.0,
+        abundance=100.0,
+        polar=False,
+        cache_tensors=False,
+    ):
+        """Basic class attributes for distributions"""
+        self._cache_tensors = cache_tensors
+        self._tensors = None
+        self.mean_isotropic_chemical_shift = mean_isotropic_chemical_shift
+        self.abundance = abundance
+        self.polar = polar
+
+    def pdf(
+        self,
+        pos,
+        size: int = 400000,
+        analytical: bool = True,
+        pack_as_csdm: bool = False,
+    ):
         """Generates a probability distribution function by binning the random
         variates of length size onto the given grid system.
 
@@ -86,21 +110,30 @@ def pdf(self, pos, size: int = 400000, analytical: bool = True):
             pos: A list of coordinates along the two dimensions given as NumPy arrays.
             size: The number of random variates drawn in generating the pdf. The default
                 is 400000.
+            pack_as_csdm: If true, returns as csdm object.
 
         Returns:
-            A list of x and y coordinates and the corresponding amplitudes.
+            A list of x and y coordinates and the corresponding amplitudes if not packed
+            as csdm object, else a csdm object.
 
         Example:
             >>> import numpy as np
             >>> cq = np.arange(50) - 25
             >>> eta = np.arange(21)/20
-            >>> Cq_dist, eta_dist, amp = cz_model.pdf(pos=[cq, eta])
+            >>> amp = cz_model.pdf(pos=[cq, eta])  # returns amp as a CSDM object.
         """
         if analytical and self.model_name in ANALYTICAL_AVAILABLE:
             analytical_model = ANALYTICAL_AVAILABLE[self.model_name]
-            return analytical_model(self.sigma, pos, self.polar)
+            pos_a, pos_b, data_ = analytical_model(self.sigma, pos, self.polar)
         else:
-            return self.pdf_numerical(pos, size)
+            pos_a, pos_b, data_ = self.pdf_numerical(pos, size)
+
+        if pack_as_csdm:
+            if len(pos_a.shape) == 1:
+                return self.pack_csdm_object([pos_a, pos_b], data_)
+            else:
+                return self.pack_csdm_object([pos_a[0, :], pos_b[:, 0]], data_)
+        return pos_a, pos_b, data_
 
     def pdf_numerical(self, pos, size: int = 400000):
         """Generate distribution numerically"""
@@ -115,9 +148,30 @@ def pdf_numerical(self, pos, size: int = 400000):
         hist, _, _ = np.histogram2d(zeta, eta, bins=[x_size, y_size], range=[x, y])
 
         hist /= hist.sum()
-
-        x_, y_ = np.meshgrid(pos[0], pos[1])
-        return x_, y_, hist.T
+        return pos[0], pos[1], hist.T
+
+    def pack_csdm_object(self, pos, data):
+        """Pack data and coordinates as csdm objects"""
+        dims = [cp.as_dimension(item) for item in pos]
+        dvs = cp.as_dependent_variable(data)
+        return cp.CSDM(dimensions=dims, dependent_variables=[dvs])
+
+    def add_lmfit_params(self, params, i):
+        """Add lmfit params for base class"""
+        prefix = self.model_name
+        params.add(
+            f"{prefix}_{i}_mean_isotropic_chemical_shift",
+            value=self.mean_isotropic_chemical_shift,
+        )
+        params.add(f"{prefix}_{i}_abundance", value=self.abundance, min=0, max=100)
+
+    def update_lmfit_params(self, params, i):
+        """Update lmfit params for base class"""
+        prefix = self.model_name
+        self.mean_isotropic_chemical_shift = params[
+            f"{prefix}_{i}_mean_isotropic_chemical_shift"
+        ].value
+        self.abundance = params[f"{prefix}_{i}_abundance"].value
 
 
 class CzjzekDistribution(AbstractDistribution):
@@ -156,9 +210,21 @@ class CzjzekDistribution(AbstractDistribution):
     """
     model_name = "czjzek"
 
-    def __init__(self, sigma: float, polar=False):
+    def __init__(
+        self,
+        sigma: float,
+        mean_isotropic_chemical_shift: float = 0.0,
+        abundance: float = 100.0,
+        polar=False,
+        cache=True,
+    ):
+        super().__init__(
+            cache_tensors=cache,
+            polar=polar,
+            mean_isotropic_chemical_shift=mean_isotropic_chemical_shift,
+            abundance=abundance,
+        )
         self.sigma = sigma
-        self.polar = polar
 
     def rvs(self, size: int):
         """Draw random variates of length `size` from the distribution.
@@ -174,10 +240,29 @@ def rvs(self, size: int):
         Example:
             >>> Cq_dist, eta_dist = cz_model.rvs(size=1000000)
         """
-        tensors = _czjzek_random_distribution_tensors(self.sigma, size)
+        if self._cache_tensors:
+            if self._tensors is None:
+                self._tensors = _czjzek_random_distribution_tensors(self.sigma, size)
+            tensors = self._tensors
+        else:
+            tensors = _czjzek_random_distribution_tensors(self.sigma, size)
+
         if not self.polar:
             return get_Haeberlen_components(tensors)
-        return x_y_from_zeta_eta(*get_Haeberlen_components(tensors))
+        return zeta_eta_to_x_y(*get_Haeberlen_components(tensors))
+
+    def add_lmfit_params(self, params, i):
+        """Create lmfit params for index i"""
+        prefix = self.model_name
+        params.add(f"{prefix}_{i}_sigma", value=self.sigma, min=0)
+        super().add_lmfit_params(params, i)
+        return params
+
+    def update_lmfit_params(self, params, i):
+        """Update lmfit params for index i"""
+        prefix = self.model_name
+        self.sigma = params[f"{prefix}_{i}_sigma"].value
+        super().update_lmfit_params(params, i)
 
 
 class ExtCzjzekDistribution(AbstractDistribution):
@@ -220,12 +305,28 @@ class ExtCzjzekDistribution(AbstractDistribution):
     >>> S0 = {"Cq": 1e6, "eta": 0.3}
     >>> ext_cz_model = ExtCzjzekDistribution(S0, eps=0.35)
     """
-    model_name = "extended czjzek"
-
-    def __init__(self, symmetric_tensor: SymmetricTensor, eps: float, polar=False):
-        self.symmetric_tensor = symmetric_tensor
+    model_name = "ext_czjzek"
+
+    def __init__(
+        self,
+        symmetric_tensor: SymmetricTensor,
+        eps: float,
+        mean_isotropic_chemical_shift: float = 0.0,
+        abundance: float = 100.0,
+        polar=False,
+        cache=True,
+    ):
+        super().__init__(
+            cache_tensors=cache,
+            polar=polar,
+            mean_isotropic_chemical_shift=mean_isotropic_chemical_shift,
+            abundance=abundance,
+        )
+        if isinstance(symmetric_tensor, dict):
+            self.symmetric_tensor = SymmetricTensor(**symmetric_tensor)
+        else:
+            self.symmetric_tensor = symmetric_tensor
         self.eps = eps
-        self.polar = polar
 
     def rvs(self, size: int):
         """Draw random variates of length `size` from the distribution.
@@ -243,13 +344,15 @@ def rvs(self, size: int):
         """
 
         # czjzek_random_distribution model
-        tensors = _czjzek_random_distribution_tensors(1, size)
+        if self._cache_tensors:
+            if self._tensors is None:
+                self._tensors = _czjzek_random_distribution_tensors(1, size)
+            tensors = self._tensors
+        else:
+            tensors = _czjzek_random_distribution_tensors(1, size)
 
         symmetric_tensor = self.symmetric_tensor
 
-        if isinstance(symmetric_tensor, dict):
-            symmetric_tensor = SymmetricTensor(**symmetric_tensor)
-
         zeta = symmetric_tensor.zeta or symmetric_tensor.Cq
         eta = symmetric_tensor.eta
 
@@ -261,12 +364,44 @@ def rvs(self, size: int):
         # 2-norm of the tensor
         norm_T0 = np.linalg.norm(T0)
 
-        # the perturbation factor
-        rho = self.eps * norm_T0 / np.sqrt(30)
+        # the perturbation factor # np.sqrt(30) = 5.477225575
+        rho = self.eps * norm_T0 / 5.4772255751
 
         # total tensor
         total_tensors = np.diag(T0) + rho * tensors
 
         if not self.polar:
             return get_Haeberlen_components(total_tensors)
-        return x_y_from_zeta_eta(*get_Haeberlen_components(total_tensors))
+        return zeta_eta_to_x_y(*get_Haeberlen_components(total_tensors))
+
+    def add_lmfit_params(self, params, i):
+        """Create lmfit params for index i"""
+        prefix = self.model_name
+        if self.symmetric_tensor.zeta is not None:
+            zeta = self.symmetric_tensor.zeta
+        else:
+            zeta = self.symmetric_tensor.Cq
+        params.add(f"{prefix}_{i}_symmetric_tensor_zeta", value=zeta)
+        params.add(
+            f"{prefix}_{i}_symmetric_tensor_eta",
+            value=self.symmetric_tensor.eta,
+            min=0,
+            max=1,
+        )
+        params.add(f"{prefix}_{i}_eps", value=self.eps, min=1e-6)
+        super().add_lmfit_params(params, i)
+        return params
+
+    def update_lmfit_params(self, params, i):
+        """Create lmfit params for index i"""
+        prefix = self.model_name
+
+        zeta = params[f"{prefix}_{i}_symmetric_tensor_zeta"].value
+        if self.symmetric_tensor.zeta is not None:
+            self.symmetric_tensor.zeta = zeta
+        else:
+            self.symmetric_tensor.Cq = zeta
+
+        self.symmetric_tensor.eta = params[f"{prefix}_{i}_symmetric_tensor_eta"].value
+        self.eps = params[f"{prefix}_{i}_eps"].value
+        super().update_lmfit_params(params, i)
diff --git a/src/mrsimulator/models/tests/test_czjzek_and_extended_czjzek.py b/src/mrsimulator/models/tests/test_czjzek_and_extended_czjzek.py
index 18aa9637a..31e15e9b5 100644
--- a/src/mrsimulator/models/tests/test_czjzek_and_extended_czjzek.py
+++ b/src/mrsimulator/models/tests/test_czjzek_and_extended_czjzek.py
@@ -3,8 +3,8 @@
 import numpy as np
 from mrsimulator.models import CzjzekDistribution
 from mrsimulator.models import ExtCzjzekDistribution
-from mrsimulator.models.utils import x_y_from_zeta_eta
 from mrsimulator.models.utils import x_y_to_zeta_eta
+from mrsimulator.models.utils import zeta_eta_to_x_y
 
 __author__ = "Deepansh J. Srivastava"
 __email__ = "srivastava.89@osu.edu"
@@ -29,7 +29,7 @@ def test_extended_czjzek_eta_distribution_1():
 def test_extended_czjzek_polar():
     S0 = {"zeta": 1, "eta": 0.1}
     x, y = ExtCzjzekDistribution(S0, eps=0.05, polar=True).rvs(size=COUNT)
-    x1, y1 = x_y_from_zeta_eta(*x_y_to_zeta_eta(x, y))
+    x1, y1 = zeta_eta_to_x_y(*x_y_to_zeta_eta(x, y))
     np.testing.assert_almost_equal(x, x1)
     np.testing.assert_almost_equal(y, y1)
 
@@ -78,10 +78,10 @@ def test_czjzek_distribution():
     z_range = (ran_z[1:] + ran_z[:-1]) / 2
 
     # czjzek distribution from analytical formula
-    _, _, res = CzjzekDistribution(sigma).pdf(pos=[z_range, e_range])
+    res = CzjzekDistribution(sigma).pdf(pos=[z_range, e_range], pack_as_csdm=True)
 
-    eta_pro = res.sum(axis=1)
-    zeta_pro = res.sum(axis=0)
+    eta_pro = res.sum(axis=0).y[0].components[0]
+    zeta_pro = res.sum(axis=1).y[0].components[0]
 
     # eta test
     message = "failed to compare eta projection for Czjzek distribution"
@@ -137,6 +137,6 @@ def test_czjzek_pdf_sub_eta():
 
 def test_czjzek_polar():
     x, y = CzjzekDistribution(sigma=0.5, polar=True).rvs(size=COUNT)
-    x1, y1 = x_y_from_zeta_eta(*x_y_to_zeta_eta(x, y))
+    x1, y1 = zeta_eta_to_x_y(*x_y_to_zeta_eta(x, y))
     np.testing.assert_almost_equal(x, x1)
     np.testing.assert_almost_equal(y, y1)
diff --git a/src/mrsimulator/models/tests/test_distribution_params.py b/src/mrsimulator/models/tests/test_distribution_params.py
new file mode 100644
index 000000000..0358d4c24
--- /dev/null
+++ b/src/mrsimulator/models/tests/test_distribution_params.py
@@ -0,0 +1,57 @@
+from lmfit import Parameters
+from mrsimulator.models import CzjzekDistribution
+from mrsimulator.models import ExtCzjzekDistribution
+
+
+def test_dist_1():
+    dist = CzjzekDistribution(sigma=1.5, mean_isotropic_chemical_shift=23.4)
+    prefix = "czjzek"
+
+    params = Parameters()
+    params = dist.add_lmfit_params(params=params, i=0)
+    assert len(params.keys()) == 3
+
+    assert params[f"{prefix}_0_sigma"] == 1.5
+    assert params[f"{prefix}_0_mean_isotropic_chemical_shift"] == 23.4
+    assert params[f"{prefix}_0_abundance"] == 100.0
+
+    params[f"{prefix}_0_sigma"].value = 10
+    params[f"{prefix}_0_mean_isotropic_chemical_shift"].value = -63.4
+    params[f"{prefix}_0_abundance"].value = 0.15
+
+    dist.update_lmfit_params(params=params, i=0)
+    assert dist.sigma == 10.0
+    assert dist.mean_isotropic_chemical_shift == -63.4
+    assert dist.abundance == 0.15
+
+
+def test_dist_2():
+    dist = ExtCzjzekDistribution(
+        symmetric_tensor={"zeta": 4.1, "eta": 0.2},
+        eps=0.3,
+        mean_isotropic_chemical_shift=0.0,
+    )
+    prefix = "ext_czjzek"
+
+    params = Parameters()
+    params = dist.add_lmfit_params(params=params, i=0)
+    assert len(params.keys()) == 5
+
+    assert params[f"{prefix}_0_symmetric_tensor_zeta"] == 4.1
+    assert params[f"{prefix}_0_symmetric_tensor_eta"] == 0.2
+    assert params[f"{prefix}_0_eps"] == 0.3
+    assert params[f"{prefix}_0_mean_isotropic_chemical_shift"] == 0.0
+    assert params[f"{prefix}_0_abundance"] == 100.0
+
+    params[f"{prefix}_0_symmetric_tensor_zeta"].value = -30.2
+    params[f"{prefix}_0_symmetric_tensor_eta"].value = 0.65
+    params[f"{prefix}_0_eps"].value = 0.6
+    params[f"{prefix}_0_mean_isotropic_chemical_shift"].value = -3.4
+    params[f"{prefix}_0_abundance"].value = 0.9
+
+    dist.update_lmfit_params(params=params, i=0)
+    assert dist.symmetric_tensor.zeta == -30.2
+    assert dist.symmetric_tensor.eta == 0.65
+    assert dist.eps == 0.6
+    assert dist.mean_isotropic_chemical_shift == -3.4
+    assert dist.abundance == 0.9
diff --git a/src/mrsimulator/models/tests/test_model_utils.py b/src/mrsimulator/models/tests/test_model_utils.py
index 570a372e4..abb72b389 100644
--- a/src/mrsimulator/models/tests/test_model_utils.py
+++ b/src/mrsimulator/models/tests/test_model_utils.py
@@ -1,8 +1,8 @@
 import numpy as np
 from mrsimulator.models.utils import get_Haeberlen_components
 from mrsimulator.models.utils import get_principal_components
-from mrsimulator.models.utils import x_y_from_zeta_eta
 from mrsimulator.models.utils import x_y_to_zeta_eta
+from mrsimulator.models.utils import zeta_eta_to_x_y
 
 
 __author__ = "Deepansh J. Srivastava"
@@ -49,10 +49,10 @@ def test_get_Haeberlen_components():
     assert components == (30, 0.5)
 
 
-def test_x_y_from_zeta_eta():
+def test_zeta_eta_to_x_y():
     zeta = np.random.rand(20) * 40
     eta = np.random.rand(20)
-    x, y = x_y_from_zeta_eta(zeta, eta)
+    x, y = zeta_eta_to_x_y(zeta, eta)
 
     theta = np.pi * eta / 4.0
     x_ = zeta * np.sin(theta)
diff --git a/src/mrsimulator/models/utils.py b/src/mrsimulator/models/utils.py
index 9bf6a4792..35821cddc 100644
--- a/src/mrsimulator/models/utils.py
+++ b/src/mrsimulator/models/utils.py
@@ -1,4 +1,13 @@
+from dataclasses import dataclass
+from dataclasses import field
+
 import numpy as np
+from mrsimulator import Method
+from mrsimulator import Simulator
+from mrsimulator import Site
+from mrsimulator import SpinSystem
+from mrsimulator.simulator import ConfigSimulator
+from mrsimulator.spin_system.isotope import Isotope
 
 __author__ = "Deepansh J. Srivastava"
 __email__ = "srivastava.89@osu.edu"
@@ -41,9 +50,10 @@ def get_Haeberlen_components(tensors):
     return zeta, eta
 
 
-def x_y_from_zeta_eta(zeta, eta):
-    """Convert the zeta, eta coordinates from the Haeberlen convention to the
-    x-y notation."""
+def zeta_eta_to_x_y(zeta, eta):
+    """Convert a set of (zeta, eta) coordinates from the Haeberlen convention to a set
+    of (x, y) coordinates.
+    """
     xa = np.empty(zeta.size)
     ya = np.empty(zeta.size)
 
@@ -64,7 +74,9 @@ def x_y_from_zeta_eta(zeta, eta):
 
 
 def x_y_to_zeta_eta(x, y):
-    """Same as def x_y_to_zeta_eta, but for ndarrays."""
+    """Convert a set of (x, y) coordinates defined by two numpy arrays into equivalent
+    (zeta, eta) coordinates.
+    """
     x = np.abs(x)
     y = np.abs(y)
     zeta = np.sqrt(x**2 + y**2)  # + offset
@@ -77,3 +89,79 @@ def x_y_to_zeta_eta(x, y):
     eta[index] = (4.0 / np.pi) * np.arctan(x[index] / y[index])
 
     return zeta, eta
+
+
+def _simulate_spectra_over_zeta_and_eta(ZZ, ee, method, tensor_type, config):
+    """Helper function to generate the kernel"""
+    isotope = Isotope.parse(
+        method.channels[0]
+    ).symbol  # Grab isotope from Method object
+
+    spin_systems = [
+        (
+            SpinSystem(sites=[Site(isotope=isotope, quadrupolar=dict(Cq=Z, eta=e))])
+            if tensor_type == "quadrupolar"
+            else SpinSystem(
+                sites=[Site(isotope=isotope, shielding_symmetric=dict(zeta=Z, eta=e))]
+            )
+        )
+        for Z, e in zip(ZZ.ravel(), ee.ravel())
+    ]
+    sim = Simulator(spin_systems=spin_systems, methods=[method])
+    sim.config = config
+    sim.config.decompose_spectrum = "spin_system"
+    sim.run(pack_as_csdm=False)
+
+    amp = sim.methods[0].simulation.real
+    return amp
+
+
+@dataclass
+class LineShapeKernel:
+    """lineshape kernel object
+
+    Arguments:
+            (tuple) pos: A tuple of numpy arrays defining the 2D grid space.
+            (bool) polar: If true, the grid is defined in polar coordinates.
+            (mrsimulator.Method) method: The :py:class:`~mrsimulator.method.Method`
+                used to simulate the spectra.
+            (ConfigSimulator) config: Simulator config to be used in simulation.
+    """
+
+    pos: list
+    method: Method
+    kernel: np.ndarray = None
+    polar: bool = False
+    config: ConfigSimulator = field(default_factory=ConfigSimulator())
+
+    def generate_lineshape(self, tensor_type: str = "shielding") -> np.ndarray:
+        """Pre-compute a lineshape kernel to use for the least-squares fitting of an
+        experimental spectrum. The isotope for the spin system is the isotope
+        at index zero of the method's channel.
+
+        Note: The lineshape kernel is currently limited to simulating the spectrum of
+        an uncoupled spin system for a single Method over a range of nuclear shielding
+        or quadrupolar tensors. Functionality may expand in the future.
+
+        Arguments:
+
+            (str) tensor_type: A string enumeration literal describing the type of
+                tensor to use in the simulation. The allowed values are `shielding` and
+                `quadrupolar`.
+
+        Returns:
+            (np.ndarray) kernel: The simulated lineshape kernel
+        """
+        if tensor_type not in {"shielding", "quadrupolar"}:
+            raise ValueError(f"Unrecognized value of {tensor_type} for `tensor_type.")
+
+        # If in polar coordinates, then (ZZ, ee) right now is really (xx, yy)
+        ZZ, ee = np.meshgrid(self.pos[0], self.pos[1], indexing="xy")
+
+        # Convert polar coords to cartesian coords
+        if self.polar:
+            ZZ, ee = x_y_to_zeta_eta(ZZ.ravel(), ee.ravel())
+
+        self.kernel = _simulate_spectra_over_zeta_and_eta(
+            ZZ, ee, self.method, tensor_type, self.config
+        )
diff --git a/src/mrsimulator/simulator/__init__.py b/src/mrsimulator/simulator/__init__.py
index c8d53fe22..b650f1452 100644
--- a/src/mrsimulator/simulator/__init__.py
+++ b/src/mrsimulator/simulator/__init__.py
@@ -1,6 +1,7 @@
 """Base Simulator class."""
 import json
 from copy import deepcopy
+from typing import Any
 from typing import List
 
 import csdmpy as cp
@@ -130,6 +131,7 @@ class Simulator(Parseable):
     """
 
     spin_systems: List[SpinSystem] = []
+    spin_system_models: List[Any] = []
     methods: List[Method] = []
     config: ConfigSimulator = ConfigSimulator()
 
diff --git a/src/mrsimulator/utils/spectral_fitting.py b/src/mrsimulator/utils/spectral_fitting.py
index 1aeb1e8d8..f32568140 100644
--- a/src/mrsimulator/utils/spectral_fitting.py
+++ b/src/mrsimulator/utils/spectral_fitting.py
@@ -1,7 +1,9 @@
+import csdmpy as cp
 import mrsimulator.signal_processor as sp
 import numpy as np
 from lmfit import Parameters
 from mrsimulator import Simulator
+from mrsimulator.models.utils import LineShapeKernel
 
 __author__ = ["Maxwell C Venetos", "Deepansh Srivastava"]
 __email__ = ["maxvenetos@gmail.com", "srivastava.89@osu.edu"]
@@ -142,7 +144,7 @@ def _traverse_dictionaries(instance, parent="spin_systems"):
     return []
 
 
-def _post_sim_LMFIT_params(params, process, index):
+def _make_params_single_processor(params, process, index):
     """Creates a LMFIT Parameters object for SignalProcessor operations involved in
     spectrum fitting.
 
@@ -179,7 +181,10 @@ def make_signal_processor_params(processors: list):
         raise ValueError("Expecting a list of `SignalProcessor` objects.")
 
     params = Parameters()
-    _ = [_post_sim_LMFIT_params(params, obj, i) for i, obj in enumerate(processors)]
+    _ = [
+        _make_params_single_processor(params, proc, i)
+        for i, proc in enumerate(processors)
+    ]
     return params
 
 
@@ -215,7 +220,9 @@ def _get_simulator_object_value(sim, string):
     return obj
 
 
-def make_simulator_params(sim: Simulator, include={}):
+def make_simulator_params(
+    sim: Simulator = Simulator(), spin_system_models: list = [], include={}
+):
     """Parse the Simulator object for a list of LMFIT parameters.
 
     Args:
@@ -228,13 +235,25 @@ def make_simulator_params(sim: Simulator, include={}):
     temp_list = _traverse_dictionaries(_list_of_dictionaries(sim.spin_systems))
 
     # get total abundance scaling factor
-    length = len(sim.spin_systems)
-    abundance_scale = 100 / sum(sys.abundance for sys in sim.spin_systems)
+    sys_length = len(sim.spin_systems)
+    sys_model_length = len(spin_system_models)
+    total_abundance = 0.0
+    total_abundance += sum(sys.abundance for sys in sim.spin_systems)
+    total_abundance += sum(sys.abundance for sys in spin_system_models)
+    abundance_scale = 100.0 / total_abundance
 
     # expression for the last abundance.
-    last_abundance = f"{length - 1}_abundance"
-    expression = "-".join([f"{START}{i}_abundance" for i in range(length - 1)])
-    expression = "100" if expression == "" else f"100-{expression}"
+    if sys_length > 0:
+        last_abundance = f"{sys_length - 1}_abundance"
+        expression = "-".join([f"{START}{i}_abundance" for i in range(sys_length - 1)])
+        expression = "100" if expression == "" else f"100-{expression}"
+
+    skip_last = sys_length > 0
+    if sys_model_length > 0:
+        param_dist = make_distribution_params(
+            spin_system_models, norm=abundance_scale, skip_last=skip_last
+        )
+        _ = params.update(param_dist)
 
     for items in temp_list:
         value = _get_simulator_object_value(sim, items)
@@ -293,7 +312,12 @@ def make_LMFIT_parameters(sim: Simulator, processors: list = None, include={}):
     return make_LMFIT_params(sim, processors)
 
 
-def make_LMFIT_params(sim: Simulator, processors: list = None, include={}):
+def make_LMFIT_params(
+    sim: Simulator = Simulator(),
+    processors: list = None,
+    spin_system_models: list = [],
+    include={},
+):
     r"""Parse the Simulator and PostSimulator objects for a list of LMFIT parameters.
 
     Args:
@@ -321,7 +345,11 @@ def make_LMFIT_params(sim: Simulator, processors: list = None, include={}):
         LMFIT Parameters object.
     """
     params = Parameters()
-    params.update(make_simulator_params(sim, include))
+    params.update(
+        make_simulator_params(
+            sim=sim, spin_system_models=spin_system_models, include=include
+        )
+    )
 
     proc = make_signal_processor_params(processors) if processors is not None else None
     _ = params.update(proc) if proc is not None else None
@@ -520,3 +548,166 @@ def residuals(sim: Simulator, processors: list = None):
         res.y[0].components[0] *= -1
 
     return residual_
+
+
+def _apply_iso_shift(csdm_obj, iso_shift_ppm, larmor_freq_Hz):
+    """Apply isotropic chemical shift to a CSDM object using the FFT shift theorem."""
+    csdm_obj = csdm_obj.fft()
+    time_coords = csdm_obj.x[0].coordinates.to("s").value
+    iso_shift_Hz = larmor_freq_Hz * iso_shift_ppm
+    csdm_obj.y[0].components[0] *= np.exp(-np.pi * 2j * time_coords * iso_shift_Hz)
+    csdm_obj = csdm_obj.fft()
+
+    return csdm_obj
+
+
+def make_distribution_params(
+    spin_system_models: list, norm: float, skip_last: bool = False
+) -> Parameters:
+    """Generate LMfit Parameters object for spin system distribution model.
+    The distribution has the following parameters:
+    """
+    n_dist = len(spin_system_models)
+
+    expr_terms = []
+    params = Parameters()
+    for i, model in enumerate(spin_system_models):
+        # normalize the abundance
+        model.abundance *= norm
+        params = model.add_lmfit_params(params, i)
+
+        # Set last abundance parameter as a non-free variable
+        model_prefix = model.model_name
+        if i < n_dist - 1 or skip_last:
+            expr_terms.append(f"{model_prefix}_{i}_abundance")
+        else:
+            expr = "-".join(expr_terms)
+            expr = "100" if expr == "" else f"100 - {expr}"
+            params[f"{model_prefix}_{i}_abundance"].vary = False
+            params[f"{model_prefix}_{i}_abundance"].expr = expr
+    return params
+
+
+def _generate_distribution_spectrum(
+    params: Parameters,
+    kernel: LineShapeKernel,
+    spin_system_models: list,
+    processor: sp.SignalProcessor = None,
+) -> cp.CSDM:
+    """Helper function for generating a spectrum from a set of LMfit Parameters and
+    and arguments for defining the grid, kernel, etc. The other functions used in least-
+    squares minimization use this function to reduce code overlap.
+
+    Arguments:
+        (Parameters) params: The LMfit parameters object holding parameters used during
+            the least-squares minimization.
+        (cp.CSDM) exp_spectrum: The experimental spectrum to fit to.
+        (tuple) pos: A tuple of two np.ndarray objects defining the grid on which to
+            sample the distribution.
+        (np.ndarray) The pre-computed lineshape kernel. The kernel needs to be defined
+            on the same grid defined by pos.
+        (float) larmor_freq_Hz: This value is used in conjunction with the FFT shift
+            theorem to apply an isotropic chemical shift to each distribution.
+        (sp.SignalProcessor) processor: A
+            :py:class:~`mrsimulator.signal_processor.Processor` object used to apply
+            post-simulation signal processing to the resulting spectrum.
+            TODO: expand processor to apply to individual distributions (as a list)
+
+    Returns:
+        Guess spectrum as a cp.CSDM object
+
+    """
+    method = kernel.method
+    larmor_freq_Hz = method.channels[0].larmor_freq(B0=method.magnetic_flux_density)
+    exp_spectrum = method.experiment
+
+    guess_spectrum = exp_spectrum.copy()
+    guess_spectrum.y[0].components[:] = 0  # Initialize guess spectrum with zeros
+    dv = cp.as_dependent_variable(np.empty(guess_spectrum.y[0].components.size))
+
+    for i, model in enumerate(spin_system_models):
+        model.update_lmfit_params(params, i)
+        _, _, amp = model.pdf(kernel.pos)
+
+        # Dot amplitude with kernel, then package as CSDM object
+        spec_tmp = cp.CSDM(dimensions=exp_spectrum.x, dependent_variables=[dv])
+        spec_tmp.y[0].components[0] = np.dot(amp.ravel(), kernel.kernel)
+
+        # Apply isotropic chemical shift to distribution using FFT shift theorem
+        spec_tmp = _apply_iso_shift(
+            csdm_obj=spec_tmp,
+            iso_shift_ppm=model.mean_isotropic_chemical_shift,
+            larmor_freq_Hz=larmor_freq_Hz,
+        ).real
+        spec_tmp *= model.abundance
+        guess_spectrum.y[0].components[0] += spec_tmp.y[0].components[0]
+
+    if processor is not None:
+        _update_processors_from_LMFIT_params(params, [processor])
+        guess_spectrum = processor.apply_operations(dataset=guess_spectrum).real
+
+    return guess_spectrum
+
+
+def LMFIT_min_function_dist(
+    params: Parameters,
+    kernel: LineShapeKernel,
+    spin_system_models: list,
+    processor: sp.SignalProcessor = None,
+) -> np.ndarray:
+    """The minimization routine for fitting a set of Czjzek models to an experimental
+    spectrum.
+
+    Arguments:
+        (Parameters) params: The LMfit parameters object holding parameters used during
+            the least-squares minimization.
+        (cp.CSDM) exp_spectrum: The experimental spectrum to fit to.
+        (tuple) pos: A tuple of two np.ndarray objects defining the grid on which to
+            sample the distribution.
+        (bool) polar: True if the sample grid is in polar coordinates. False if the grid
+            is defined using the Haberlen components.
+        (np.ndarray) The pre-computed lineshape kernel. The kernel needs to be defined
+            on the same grid defined by pos.
+        (int) n_dist: The number of Czjzek distributions to fit for.
+        (float) larmor_freq_Hz: This value is used in conjunction with the FFT shift
+            theorem to apply an isotropic chemical shift to each distribution.
+        (sp.SignalProcessor) processor: A
+            :py:class:~`mrsimulator.signal_processor.Processor` object used to apply
+            post-simulation signal processing to the resulting spectrum.
+            TODO: expand processor to apply to individual distributions (as a list)
+
+        Returns:
+            A residual array as a numpy array.
+    """
+    guess_spectrum = _generate_distribution_spectrum(
+        params, kernel, spin_system_models, processor
+    )
+    exp_spectrum = kernel.method.experiment
+    return (exp_spectrum - guess_spectrum).y[0].components[0]
+
+
+def bestfit_dist(
+    params: Parameters,
+    kernel: LineShapeKernel,
+    spin_system_models: list,
+    processor: sp.SignalProcessor = None,
+) -> cp.CSDM:
+    """Returns the best-fit spectrum as a CSDM object"""
+    return _generate_distribution_spectrum(
+        params, kernel, spin_system_models, processor
+    )
+
+
+def residuals_dist(
+    params: Parameters,
+    kernel: LineShapeKernel,
+    spin_system_models: list,
+    processor: sp.SignalProcessor = None,
+) -> cp.CSDM:
+    """Returns the residuals spectrum as a CSDM object"""
+    bestfit_ = _generate_distribution_spectrum(
+        params, kernel, spin_system_models, processor
+    )
+
+    exp_spectrum = kernel.method.experiment
+    return exp_spectrum - bestfit_
diff --git a/src/mrsimulator/utils/tests/test_spectral_fitting.py b/src/mrsimulator/utils/tests/test_spectral_fitting.py
index 9ae49f5c3..031e95d97 100644
--- a/src/mrsimulator/utils/tests/test_spectral_fitting.py
+++ b/src/mrsimulator/utils/tests/test_spectral_fitting.py
@@ -10,6 +10,8 @@
 from mrsimulator.method.lib import BlochDecaySpectrum
 from mrsimulator.method.lib import SSB2D
 from mrsimulator.method.lib import ThreeQ_VAS
+from mrsimulator.models import CzjzekDistribution
+from mrsimulator.models import ExtCzjzekDistribution
 
 
 def test_str_encode():
@@ -357,3 +359,27 @@ def test_array():
     # params = sf.make_LMFIT_params(sim, processor)
     # a = sf.LMFIT_min_function(params, sim, processor)
     # np.testing.assert_almost_equal(-a.sum(), dataset.sum().real, decimal=8)
+
+
+def test_distribution():
+    dist1 = CzjzekDistribution(sigma=0.5)
+    dist2 = CzjzekDistribution(sigma=1.4, mean_isotropic_chemical_shift=4.7)
+    dist3 = ExtCzjzekDistribution(symmetric_tensor={"Cq": 1e6, "eta": 0.1}, eps=0.5)
+
+    models = [dist1, dist2, dist3]
+
+    params = sf.make_LMFIT_params(spin_system_models=models)
+
+    assert params["czjzek_0_sigma"] == 0.5
+    assert params["czjzek_0_mean_isotropic_chemical_shift"] == 0.0
+    assert np.allclose(params["czjzek_0_abundance"], 100.0 / 3.0)
+
+    assert params["czjzek_1_sigma"] == 1.4
+    assert params["czjzek_1_mean_isotropic_chemical_shift"] == 4.7
+    assert np.allclose(params["czjzek_1_abundance"], 100.0 / 3.0)
+
+    assert params["ext_czjzek_2_symmetric_tensor_zeta"] == 1e6
+    assert params["ext_czjzek_2_symmetric_tensor_eta"] == 0.1
+    assert params["ext_czjzek_2_eps"] == 0.5
+    assert params["ext_czjzek_2_mean_isotropic_chemical_shift"] == 0.0
+    assert np.allclose(params["ext_czjzek_2_abundance"], 100.0 / 3.0)