diff --git a/.gitignore b/.gitignore
index 15c86627d..b00dbbba8 100644
--- a/.gitignore
+++ b/.gitignore
@@ -22,3 +22,5 @@ freud/*.cpp
 freud/*.cxx
 freud/*.dylib
 freud/*.so
+
+.vscode/*
diff --git a/cpp/diffraction/CMakeLists.txt b/cpp/diffraction/CMakeLists.txt
index 03e43b439..09522bb4c 100644
--- a/cpp/diffraction/CMakeLists.txt
+++ b/cpp/diffraction/CMakeLists.txt
@@ -1,8 +1,16 @@
 add_library(
   _diffraction OBJECT
-  StaticStructureFactor.h StaticStructureFactor.cc StaticStructureFactorDebye.h
-  StaticStructureFactorDebye.cc StaticStructureFactorDirect.h
-  StaticStructureFactorDirect.cc)
+  StructureFactor.h
+  StaticStructureFactor.h
+  StaticStructureFactor.cc
+  StaticStructureFactorDebye.h
+  StaticStructureFactorDebye.cc
+  StructureFactorDirect.h
+  StructureFactorDirect.cc
+  StaticStructureFactorDirect.h
+  StaticStructureFactorDirect.cc
+  IntermediateScattering.h
+  IntermediateScattering.cc)
 
 target_link_libraries(_diffraction PUBLIC TBB::tbb)
 
diff --git a/cpp/diffraction/IntermediateScattering.cc b/cpp/diffraction/IntermediateScattering.cc
new file mode 100644
index 000000000..d7a461868
--- /dev/null
+++ b/cpp/diffraction/IntermediateScattering.cc
@@ -0,0 +1,167 @@
+// Copyright (c) 2010-2020 The Regents of the University of Michigan
+// This file is from the freud project, released under the BSD 3-Clause License.
+
+#include "IntermediateScattering.h"
+#include "Box.h"
+#include "ManagedArray.h"
+#include "NeighborQuery.h"
+#include "utils.h"
+
+/*! \file IntermediateScattering.cc
+    \brief Routines for computing intermediate scattering function.
+*/
+
+namespace freud { namespace diffraction {
+
+IntermediateScattering::IntermediateScattering(const box::Box& box, unsigned int bins, float k_max,
+                                               float k_min, unsigned int num_sampled_k_points)
+    : StructureFactorDirect(bins, k_max, k_min, num_sampled_k_points), StructureFactor(bins, k_max, k_min),
+      m_box(box)
+{
+    if (bins == 0)
+    {
+        throw std::invalid_argument("IntermediateScattering requires a nonzero number of bins.");
+    }
+    if (k_max <= 0)
+    {
+        throw std::invalid_argument("IntermediateScattering requires k_max to be positive.");
+    }
+    if (k_min < 0)
+    {
+        throw std::invalid_argument("IntermediateScattering requires k_min to be non-negative.");
+    }
+    if (k_max <= k_min)
+    {
+        throw std::invalid_argument("IntermediateScattering requires that k_max must be greater than k_min.");
+    }
+
+    m_k_points = IntermediateScattering::reciprocal_isotropic(box, k_max, k_min, m_num_sampled_k_points);
+}
+
+void IntermediateScattering::compute(const vec3<float>* points, unsigned int num_points,
+                                     const vec3<float>* query_points, unsigned int num_query_points,
+                                     unsigned int num_frames, unsigned int n_total)
+{
+    // initialize the histograms, now that the size of both axes are known
+    m_self_function = StructureFactorHistogram({
+        std::make_shared<util::RegularAxis>(num_frames, -0.5, num_frames - 0.5),
+        std::make_shared<util::RegularAxis>(m_nbins, m_k_min, m_k_max),
+    });
+    m_local_self_function = StructureFactorHistogram::ThreadLocalHistogram(m_self_function);
+    m_distinct_function = StructureFactorHistogram({
+        std::make_shared<util::RegularAxis>(num_frames, -0.5, num_frames - 0.5),
+        std::make_shared<util::RegularAxis>(m_nbins, m_k_min, m_k_max),
+    });
+    m_local_distinct_function = StructureFactorHistogram::ThreadLocalHistogram(m_distinct_function);
+    // Compute k vectors by sampling reciprocal space.
+    if (m_box.is2D())
+    {
+        throw std::invalid_argument("2D boxes are not currently supported.");
+    }
+
+    // The minimum valid k value is 2 * pi / L, where L is the smallest side length.
+    const auto box_L = m_box.getL();
+    const auto min_box_length = std::min(box_L.x, std::min(box_L.y, box_L.z));
+    m_min_valid_k = freud::constants::TWO_PI / min_box_length;
+
+    // record the point at t=0
+    const vec3<float>* query_r0 = &query_points[0];
+    const vec3<float>* r0 = &points[0];
+
+    util::forLoopWrapper(0, num_frames, [&](size_t begin, size_t end) {
+        for (size_t t = begin; t < end; t++)
+        {
+            size_t offset = t * num_points;
+            // Compute self-part
+            const auto self_part
+                = IntermediateScattering::compute_self(&points[offset], r0, num_points, n_total, m_k_points);
+
+            // Compute distinct-part
+            const auto distinct_part = IntermediateScattering::compute_distinct(
+                &points[offset], query_r0, num_query_points, n_total, m_k_points);
+
+            std::vector<float> S_k_self_part = StaticStructureFactorDirect::compute_S_k(self_part, self_part);
+            std::vector<float> S_k_distinct_part
+                = StaticStructureFactorDirect::compute_S_k(distinct_part, distinct_part);
+
+            // Bin the S_k values and track the number of k values in each bin.
+            util::forLoopWrapper(0, m_k_points.size(), [&](size_t begin, size_t end) {
+                for (size_t k_index = begin; k_index < end; ++k_index)
+                {
+                    const auto& k_vec = m_k_points[k_index];
+                    const auto k_magnitude = std::sqrt(dot(k_vec, k_vec));
+                    const auto k_bin = m_k_histogram.bin({k_magnitude});
+                    m_local_self_function.increment(k_bin, S_k_self_part[k_index]);
+                    m_local_distinct_function.increment(k_bin, S_k_distinct_part[k_index]);
+                    m_local_k_histograms.increment(k_bin);
+                }
+            });
+        }
+    });
+
+    m_reduce = true;
+}
+
+void IntermediateScattering::reduce()
+{
+    const auto k_axis_size = m_k_histogram.getAxisSizes();
+    m_k_histogram.prepare(k_axis_size);
+    const auto self_axis_size = m_self_function.getAxisSizes();
+    m_self_function.prepare(self_axis_size);
+    const auto distinct_axis_size = m_distinct_function.getAxisSizes();
+    m_distinct_function.prepare(distinct_axis_size);
+
+    // Reduce the bin counts over all threads, then use them to normalize the
+    // structure factor when computing. This computes a binned mean over all k
+    // points.
+    m_k_histogram.reduceOverThreads(m_local_k_histograms);
+    m_self_function.reduceOverThreadsPerBin(m_local_self_function, [&](size_t i) {
+        m_self_function[i] /= m_k_histogram[util::modulusPositive(i, m_nbins)];
+    });
+    m_distinct_function.reduceOverThreadsPerBin(m_local_distinct_function, [&](size_t i) {
+        m_distinct_function[i] /= m_k_histogram[util::modulusPositive(i, m_nbins)];
+    });
+}
+
+std::vector<std::complex<float>>
+IntermediateScattering::compute_self(const vec3<float>* rt, const vec3<float>* r0, unsigned int n_points,
+                                     unsigned int n_total, const std::vector<vec3<float>>& k_points)
+{
+    //
+    std::vector<vec3<float>> r_i_t0(n_points); // rt - rt0 element-wisely
+    util::forLoopWrapper(0, n_points, [&](size_t begin, size_t end) {
+        for (size_t i = begin; i < end; ++i)
+        {
+            r_i_t0[i] = rt[i] - rt[0];
+        }
+    });
+
+    return StructureFactorDirect::compute_F_k(r_i_t0.data(), n_points, n_total, m_k_points);
+}
+
+std::vector<std::complex<float>>
+IntermediateScattering::compute_distinct(const vec3<float>* rt, const vec3<float>* r0, unsigned int n_points,
+                                         unsigned int n_total, const std::vector<vec3<float>>& k_points)
+{
+    const auto n_rij = n_points * (n_points - 1);
+    std::vector<vec3<float>> r_ij(n_rij);
+    size_t i = 0;
+
+    util::forLoopWrapper(0, n_rij, [&](size_t begin, size_t end) {
+        for (size_t rt_index = begin; rt_index < end; ++rt_index)
+        {
+            for (size_t r0_index = 0; r0_index < n_rij; ++r0_index)
+            {
+                if (rt_index != r0_index)
+                {
+                    r_ij[i] = rt[rt_index] - r0[r0_index];
+                    ++i;
+                }
+            }
+        };
+    });
+
+    return StructureFactorDirect::compute_F_k(r_ij.data(), n_rij, n_total, m_k_points);
+}
+
+}} // namespace freud::diffraction
diff --git a/cpp/diffraction/IntermediateScattering.h b/cpp/diffraction/IntermediateScattering.h
new file mode 100644
index 000000000..874d32d0b
--- /dev/null
+++ b/cpp/diffraction/IntermediateScattering.h
@@ -0,0 +1,101 @@
+// Copyright (c) 2010-2020 The Regents of the University of Michigan
+// This file is from the freud project, released under the BSD 3-Clause License.
+
+#ifndef INTERMEDIATE_SCATTERING_H
+#define INTERMEDIATE_SCATTERING_H
+
+#include <complex>
+#include <limits>
+#include <vector>
+
+#include "Box.h"
+#include "Histogram.h"
+#include "ManagedArray.h"
+#include "NeighborQuery.h"
+#include "StaticStructureFactorDirect.h"
+
+/*! \file IntermediateScattering.h
+    \brief Routines for computing intermediate scattering function.
+
+    This computes the intermediate scattering function F(k, t) between sets of
+    time-dependent points and query points.
+
+    The k-vectors are used to compute complex scattering amplitudes F(k, t) by
+    summing over the scattering contribution of all particle positions (atomic
+    form factors are assumed to be 1) at every k-vector. The complex conjugate
+    of the scattering amplitudes of the points at each k-vector are multiplied
+    by the scattering amplitudes of the query points at each k-vector. Finally,
+    the results are binned according to their k-vector's magnitude and normalized
+    by the number of samples in each radial bin in k-space.
+
+    Note that k-vectors are in the physics convention, and q-vectors are in the
+    crystallographic convention. These conventions differ by a factor of 2\pi.
+
+    See also:
+    https://en.wikipedia.org/wiki/Reciprocal_lattice#Arbitrary_collection_of_atoms
+*/
+
+namespace freud { namespace diffraction {
+
+class IntermediateScattering : public StructureFactorDirect
+{
+public:
+    //! Constructor
+    IntermediateScattering(const box::Box& box, unsigned int bins, float k_max, float k_min = 0,
+                           unsigned int num_sampled_k_points = 0);
+
+    void compute(const vec3<float>* points, unsigned int num_points, const vec3<float>* query_points,
+                 unsigned int num_query_points, unsigned int num_frames, unsigned int n_total);
+
+    std::vector<float> getBinEdges() const
+    {
+        // the second axis of the histogram is the k-axis
+        return m_self_function.getBinEdges()[1];
+    }
+
+    std::vector<float> getBinCenters() const
+    {
+        // the second axis of the histogram is the k-axis
+        return m_self_function.getBinCenters()[1];
+    }
+
+    // more intuitively named alias for getStructureFactor
+    const util::ManagedArray<float>& getSelfFunction()
+    {
+        return reduceAndReturn(m_self_function.getBinCounts());
+    }
+
+    const util::ManagedArray<float>& getDistinctFunction()
+    {
+        return reduceAndReturn(m_distinct_function.getBinCounts());
+    }
+
+private:
+    void reduce() override;
+
+    //!< box for the calculation, we assume the box is constant over the time interval
+    box::Box m_box;
+
+    //!< Histogram to hold self part of the ISF at each frame
+    StructureFactorHistogram m_self_function {};
+    //!< Thread local histograms for TBB parallelism
+    StructureFactorHistogram::ThreadLocalHistogram m_local_self_function {};
+
+    //!< Histogram to hold distinct part of the ISF at each frame
+    StructureFactorHistogram m_distinct_function {};
+    //!< Thread local histograms for TBB parallelism
+    StructureFactorHistogram::ThreadLocalHistogram m_local_distinct_function {};
+
+    //!< Helpers to compute self and distinct parts
+    std::vector<std::complex<float>> compute_self(const vec3<float>* rt, const vec3<float>* r0,
+                                                  unsigned int n_points, unsigned int n_total,
+                                                  const std::vector<vec3<float>>& k_points);
+
+    std::vector<std::complex<float>> compute_distinct(const vec3<float>* rt, const vec3<float>* r0,
+                                                      unsigned int n_points, unsigned int n_total,
+                                                      const std::vector<vec3<float>>& k_points);
+};
+
+}; }; // namespace freud::diffraction
+
+#endif // INTERMEDIATE_SCATTERING_H
diff --git a/cpp/diffraction/StaticStructureFactor.cc b/cpp/diffraction/StaticStructureFactor.cc
index 5251693e9..ab4bd8a9e 100644
--- a/cpp/diffraction/StaticStructureFactor.cc
+++ b/cpp/diffraction/StaticStructureFactor.cc
@@ -13,7 +13,8 @@
 namespace freud { namespace diffraction {
 
 StaticStructureFactor::StaticStructureFactor(unsigned int bins, float k_max, float k_min)
-    : m_structure_factor({std::make_shared<util::RegularAxis>(bins, k_min, k_max)}),
+    : StructureFactor(bins, k_max, k_min),
+      m_structure_factor({std::make_shared<util::RegularAxis>(bins, k_min, k_max)}),
       m_local_structure_factor(m_structure_factor)
 {
     // Validation logic is not shared in the parent StaticStructureFactor
diff --git a/cpp/diffraction/StaticStructureFactor.h b/cpp/diffraction/StaticStructureFactor.h
index e0d5c57a6..835d9dc23 100644
--- a/cpp/diffraction/StaticStructureFactor.h
+++ b/cpp/diffraction/StaticStructureFactor.h
@@ -10,28 +10,43 @@
 #include "Histogram.h"
 #include "ManagedArray.h"
 #include "NeighborQuery.h"
+#include "StructureFactor.h"
 
 /*! \file StaticStructureFactor.h
-    \brief Base class for structure factor classes.
+    \brief Base class for static tructure factor classes.
 */
 
 namespace freud { namespace diffraction {
 
-class StaticStructureFactor
+/* Abstract base class for all static structure factors.
+ *
+ * A static structure factor is a structure factor which can be computed using
+ * only the data from one frame of a simulation. A typical use case is to compute
+ * the static structure factor for each frame over many frames of a simulation,
+ * and get an average for better statistics/curve smoothness. To support this use
+ * case, all static structure factors must have logic for either continuing to
+ * accumulate histogram data or resetting the data on each successive call to
+ * accumulate().
+ *
+ * This class does not have any members, it simply defines the interface for
+ * all static structure factors.
+ *
+ * */
+class StaticStructureFactor : virtual public StructureFactor
 {
 protected:
-    using StructureFactorHistogram = util::Histogram<float>;
-
     StaticStructureFactor(unsigned int bins, float k_max, float k_min = 0);
 
 public:
     virtual ~StaticStructureFactor() = default;
 
+    //!< do the histogram calculation, and accumulate if necessary
     virtual void accumulate(const freud::locality::NeighborQuery* neighbor_query,
                             const vec3<float>* query_points, unsigned int n_query_points,
                             unsigned int n_total)
         = 0;
 
+    //<! reset the structure factor if the user doesn't wish to accumulate
     virtual void reset() = 0;
 
     //! Get the structure factor
@@ -52,32 +67,11 @@ class StaticStructureFactor
         return m_structure_factor.getBinCenters()[0];
     }
 
-    //! Get the minimum valid k value
-    float getMinValidK() const
-    {
-        return m_min_valid_k;
-    }
-
 protected:
-    virtual void reduce() = 0;
-
-    //! Return thing_to_return after reducing if necessary.
-    template<typename U> U& reduceAndReturn(U& thing_to_return)
-    {
-        if (m_reduce)
-        {
-            reduce();
-        }
-        m_reduce = false;
-        return thing_to_return;
-    }
-
-    StructureFactorHistogram m_structure_factor; //!< Histogram to hold computed structure factor
-    StructureFactorHistogram::ThreadLocalHistogram
-        m_local_structure_factor; //!< Thread local histograms for TBB parallelism
-
-    bool m_reduce {true};                                         //! Whether to reduce local histograms
-    float m_min_valid_k {std::numeric_limits<float>::infinity()}; //! Minimum valid k-vector magnitude
+    //!< Histogram to hold computed structure factor
+    StructureFactorHistogram m_structure_factor;
+    //!< Thread local histograms for TBB parallelism
+    StructureFactorHistogram::ThreadLocalHistogram m_local_structure_factor;
 };
 
 }; }; // namespace freud::diffraction
diff --git a/cpp/diffraction/StaticStructureFactorDebye.cc b/cpp/diffraction/StaticStructureFactorDebye.cc
index ebb0e55d3..620ba13df 100644
--- a/cpp/diffraction/StaticStructureFactorDebye.cc
+++ b/cpp/diffraction/StaticStructureFactorDebye.cc
@@ -35,7 +35,9 @@ float k_min_center_to_lower_edge(unsigned int bins, float k_min, float k_max)
 
 StaticStructureFactorDebye::StaticStructureFactorDebye(unsigned int bins, float k_max, float k_min)
     : StaticStructureFactor(bins, k_max_center_to_upper_edge(bins, k_min, k_max),
-                            k_min_center_to_lower_edge(bins, k_min, k_max))
+                            k_min_center_to_lower_edge(bins, k_min, k_max)),
+      StructureFactor(bins, k_max_center_to_upper_edge(bins, k_min, k_max),
+                      k_min_center_to_lower_edge(bins, k_min, k_max))
 {
     if (bins == 0)
     {
diff --git a/cpp/diffraction/StaticStructureFactorDirect.cc b/cpp/diffraction/StaticStructureFactorDirect.cc
index 64abda3b3..c87607b44 100644
--- a/cpp/diffraction/StaticStructureFactorDirect.cc
+++ b/cpp/diffraction/StaticStructureFactorDirect.cc
@@ -8,8 +8,6 @@
 #include <stdexcept>
 #include <tbb/concurrent_vector.h>
 
-#include "Eigen/Eigen/Dense"
-
 #include "Box.h"
 #include "ManagedArray.h"
 #include "NeighborQuery.h"
@@ -24,9 +22,8 @@ namespace freud { namespace diffraction {
 
 StaticStructureFactorDirect::StaticStructureFactorDirect(unsigned int bins, float k_max, float k_min,
                                                          unsigned int num_sampled_k_points)
-    : StaticStructureFactor(bins, k_max, k_min), m_num_sampled_k_points(num_sampled_k_points),
-      m_k_histogram(KBinHistogram(m_structure_factor.getAxes())),
-      m_local_k_histograms(KBinHistogram::ThreadLocalHistogram(m_k_histogram))
+    : StaticStructureFactor(bins, k_max, k_min),
+      StructureFactorDirect(bins, k_max, k_min, num_sampled_k_points), StructureFactor(bins, k_max, k_min)
 {
     if (bins == 0)
     {
@@ -63,8 +60,7 @@ void StaticStructureFactorDirect::accumulate(const freud::locality::NeighborQuer
     if ((!box_assigned) || (box != previous_box))
     {
         previous_box = box;
-        m_k_points
-            = StaticStructureFactorDirect::reciprocal_isotropic(box, k_max, k_min, m_num_sampled_k_points);
+        m_k_points = StructureFactorDirect::reciprocal_isotropic(box, k_max, k_min, m_num_sampled_k_points);
         box_assigned = true;
     }
 
@@ -75,7 +71,7 @@ void StaticStructureFactorDirect::accumulate(const freud::locality::NeighborQuer
     m_min_valid_k = std::min(m_min_valid_k, freud::constants::TWO_PI / min_box_length);
 
     // Compute F_k for the points.
-    const auto F_k_points = StaticStructureFactorDirect::compute_F_k(
+    const auto F_k_points = StructureFactorDirect::compute_F_k(
         neighbor_query->getPoints(), neighbor_query->getNPoints(), n_total, m_k_points);
 
     // Compute F_k for the query points (if necessary) and compute the product S_k.
@@ -83,12 +79,12 @@ void StaticStructureFactorDirect::accumulate(const freud::locality::NeighborQuer
     if (query_points != nullptr)
     {
         const auto F_k_query_points
-            = StaticStructureFactorDirect::compute_F_k(query_points, n_query_points, n_total, m_k_points);
-        S_k_all_points = StaticStructureFactorDirect::compute_S_k(F_k_points, F_k_query_points);
+            = StructureFactorDirect::compute_F_k(query_points, n_query_points, n_total, m_k_points);
+        S_k_all_points = StructureFactorDirect::compute_S_k(F_k_points, F_k_query_points);
     }
     else
     {
-        S_k_all_points = StaticStructureFactorDirect::compute_S_k(F_k_points, F_k_points);
+        S_k_all_points = StructureFactorDirect::compute_S_k(F_k_points, F_k_points);
     }
 
     // Bin the S_k values and track the number of k values in each bin.
@@ -120,176 +116,4 @@ void StaticStructureFactorDirect::reduce()
                                                [&](size_t i) { m_structure_factor[i] /= m_k_histogram[i]; });
 }
 
-std::vector<std::complex<float>>
-StaticStructureFactorDirect::compute_F_k(const vec3<float>* points, unsigned int n_points,
-                                         unsigned int n_total, const std::vector<vec3<float>>& k_points)
-{
-    const auto n_k_points = k_points.size();
-    auto F_k = std::vector<std::complex<float>>(n_k_points);
-    const std::complex<float> normalization(1.0F / std::sqrt(static_cast<float>(n_total)));
-
-    util::forLoopWrapper(0, n_k_points, [&](size_t begin, size_t end) {
-        for (size_t k_index = begin; k_index < end; ++k_index)
-        {
-            std::complex<float> F_ki(0);
-            for (size_t r_index = 0; r_index < n_points; ++r_index)
-            {
-                const auto& k_vec(k_points[k_index]);
-                const auto& r_vec(points[r_index]);
-                const auto alpha(dot(k_vec, r_vec));
-                F_ki += std::exp(std::complex<float>(0, alpha));
-            }
-            F_k[k_index] = F_ki * normalization;
-        }
-    });
-    return F_k;
-}
-
-std::vector<float>
-StaticStructureFactorDirect::compute_S_k(const std::vector<std::complex<float>>& F_k_points,
-                                         const std::vector<std::complex<float>>& F_k_query_points)
-{
-    const auto n_k_points = F_k_points.size();
-    auto S_k = std::vector<float>(n_k_points);
-    util::forLoopWrapper(0, n_k_points, [&](size_t begin, size_t end) {
-        for (size_t k_index = begin; k_index < end; ++k_index)
-        {
-            S_k[k_index] = std::real(std::conj(F_k_points[k_index]) * F_k_query_points[k_index]);
-        }
-    });
-    return S_k;
-}
-
-inline Eigen::Matrix3f box_to_matrix(const box::Box& box)
-{
-    // Build an Eigen matrix from the provided box.
-    Eigen::Matrix3f mat;
-    for (unsigned int i = 0; i < 3; i++)
-    {
-        const auto box_vector = box.getLatticeVector(i);
-        mat(i, 0) = box_vector.x;
-        mat(i, 1) = box_vector.y;
-        mat(i, 2) = box_vector.z;
-    }
-    return mat;
-}
-
-inline float get_prune_distance(unsigned int num_sampled_k_points, float q_max, float q_volume)
-{
-    if ((num_sampled_k_points > M_PI * std::pow(q_max, 3.0) / (6 * q_volume)) || (num_sampled_k_points == 0))
-    {
-        // Above this limit, all points are used and no pruning occurs.
-        return std::numeric_limits<float>::infinity();
-    }
-    // We use Cardano's formula to compute the pruning distance.
-    const auto p = -0.75F * q_max * q_max;
-    const auto q = 3.0F * static_cast<float>(num_sampled_k_points) * q_volume / static_cast<float>(M_PI)
-        - q_max * q_max * q_max / 4.0F;
-    const auto D = p * p * p / 27.0F + q * q / 4.0F;
-
-    const auto u = std::pow(-std::complex<float>(q / 2.0F) + std::sqrt(std::complex<float>(D)), 1.0F / 3.0F);
-    const auto v = std::pow(-std::complex<float>(q / 2.0F) - std::sqrt(std::complex<float>(D)), 1.0F / 3.0F);
-    const auto x = -(u + v) / 2.0F - std::complex<float>(0.0F, 1.0F) * (u - v) * std::sqrt(3.0F) / 2.0F;
-    return std::real(x) + q_max / 2.0F;
-}
-
-std::vector<vec3<float>> StaticStructureFactorDirect::reciprocal_isotropic(const box::Box& box, float k_max,
-                                                                           float k_min,
-                                                                           unsigned int num_sampled_k_points)
-{
-    const auto box_matrix = box_to_matrix(box);
-    // B holds "crystallographic" reciprocal box vectors that lack the factor of 2 pi.
-    const auto B = box_matrix.transpose().inverse();
-    const auto q_max = k_max / freud::constants::TWO_PI;
-    const auto q_max_sq = q_max * q_max;
-    const auto q_min = k_min / freud::constants::TWO_PI;
-    const auto q_min_sq = q_min * q_min;
-    const auto dq_x = B.row(0).norm();
-    const auto dq_y = B.row(1).norm();
-    const auto dq_z = B.row(2).norm();
-    const auto q_volume = dq_x * dq_y * dq_z;
-
-    // Above the pruning distance, the grid of k points is sampled isotropically
-    // at a lower density.
-    const auto q_prune_distance = get_prune_distance(num_sampled_k_points, q_max, q_volume);
-    const auto q_prune_distance_sq = q_prune_distance * q_prune_distance;
-
-    const auto bx = freud::constants::TWO_PI * vec3<float>(B(0, 0), B(0, 1), B(0, 2));
-    const auto by = freud::constants::TWO_PI * vec3<float>(B(1, 0), B(1, 1), B(1, 2));
-    const auto bz = freud::constants::TWO_PI * vec3<float>(B(2, 0), B(2, 1), B(2, 2));
-    const auto N_kx = static_cast<unsigned int>(std::ceil(q_max / dq_x));
-    const auto N_ky = static_cast<unsigned int>(std::ceil(q_max / dq_y));
-    const auto N_kz = static_cast<unsigned int>(std::ceil(q_max / dq_z));
-
-    // The maximum number of k points is a guideline. The true number of sampled
-    // k points can be less or greater than num_sampled_k_points, depending on the
-    // result of the random pruning procedure. Therefore, we cannot allocate a
-    // fixed size for the data. Also, reserving capacity for the concurrent
-    // vector had no measureable effect on performance.
-    tbb::concurrent_vector<vec3<float>> k_points;
-
-    // This is a 3D loop but we parallelize in 1D because we only need to seed
-    // the random number generators once per block and there is no benefit of
-    // locality if we parallelize in 2D or 3D.
-    util::forLoopWrapper(0, N_kx, [&](size_t begin, size_t end) {
-        // Set up thread-local random number generator for k point pruning.
-        const auto thread_start = static_cast<unsigned int>(begin);
-        std::random_device rd;
-        std::seed_seq seed {thread_start, rd(), rd(), rd()};
-        std::mt19937 rng(seed);
-        std::uniform_real_distribution<float> base_dist(0, 1);
-        auto random_prune = [&]() { return base_dist(rng); };
-
-        const auto add_all_k_points = std::isinf(q_prune_distance);
-
-        for (unsigned int kx = begin; kx < end; ++kx)
-        {
-            const auto k_vec_x = static_cast<float>(kx) * bx;
-            for (unsigned int ky = 0; ky < N_ky; ++ky)
-            {
-                const auto k_vec_xy = k_vec_x + static_cast<float>(ky) * by;
-
-                // Solve the quadratic equation for kz to limit which kz values we must sample:
-                // k_min^2 <= |k_vec_xy|^2 + kz^2 |bz|^2 - 2 kz (k_vec_xy \cdot bz) <= k_max^2
-                // 0 <= kz^2 (|bz|^2) + kz (-2 (k_vec_xy \cdot bz)) + (|k_vec_xy|^2 - k_min^2)
-                // 0 >= kz^2 (|bz|^2) + kz (-2 (k_vec_xy \cdot bz)) + (|k_vec_xy|^2 - k_max^2)
-                // This step improves performance significantly when k_min > 0
-                // by eliminating a large portion of the search space. Likewise,
-                // it eliminates the portion of search space outside a sphere
-                // with radius k_max. We round kz_min down and kz_max up to
-                // ensure that we don't accidentally throw out valid k points in
-                // the range (k_min, k_max) due to rounding error.
-                const auto coef_a = dot(bz, bz);
-                const auto coef_b = -2 * dot(k_vec_xy, bz);
-                const auto coef_c_min = dot(k_vec_xy, k_vec_xy) - k_min * k_min;
-                const auto coef_c_max = dot(k_vec_xy, k_vec_xy) - k_max * k_max;
-                const auto b_over_2a = coef_b / (2 * coef_a);
-                const auto kz_min = static_cast<unsigned int>(
-                    std::floor(-b_over_2a + std::sqrt(b_over_2a * b_over_2a - coef_c_min / coef_a)));
-                const auto kz_max = static_cast<unsigned int>(
-                    std::ceil(-b_over_2a + std::sqrt(b_over_2a * b_over_2a - coef_c_max / coef_a)));
-                for (unsigned int kz = kz_min; kz < std::min(kz_max, N_kz); ++kz)
-                {
-                    const auto k_vec = k_vec_xy + static_cast<float>(kz) * bz;
-                    const auto q_distance_sq
-                        = dot(k_vec, k_vec) / freud::constants::TWO_PI / freud::constants::TWO_PI;
-
-                    // The k vector is kept with probability min(1, (q_prune_distance / q_distance)^2).
-                    // This sampling scheme aims to have a constant density of k vectors with respect to
-                    // radial distance.
-                    if (q_distance_sq <= q_max_sq && q_distance_sq >= q_min_sq)
-                    {
-                        const auto prune_probability = q_prune_distance_sq / q_distance_sq;
-                        if (add_all_k_points || prune_probability > random_prune())
-                        {
-                            k_points.emplace_back(k_vec);
-                        }
-                    }
-                }
-            }
-        }
-    });
-    return std::vector<vec3<float>>(k_points.cbegin(), k_points.cend());
-}
-
 }; }; // namespace freud::diffraction
diff --git a/cpp/diffraction/StaticStructureFactorDirect.h b/cpp/diffraction/StaticStructureFactorDirect.h
index c159f5cb5..49fa98857 100644
--- a/cpp/diffraction/StaticStructureFactorDirect.h
+++ b/cpp/diffraction/StaticStructureFactorDirect.h
@@ -13,6 +13,7 @@
 #include "ManagedArray.h"
 #include "NeighborQuery.h"
 #include "StaticStructureFactor.h"
+#include "StructureFactorDirect.h"
 
 /*! \file StaticStructureFactorDirect.h
     \brief Routines for computing static structure factors.
@@ -41,10 +42,8 @@
 
 namespace freud { namespace diffraction {
 
-class StaticStructureFactorDirect : public StaticStructureFactor
+class StaticStructureFactorDirect : public StaticStructureFactor, public StructureFactorDirect
 {
-    using KBinHistogram = util::Histogram<unsigned int>;
-
 public:
     //! Constructor
     StaticStructureFactorDirect(unsigned int bins, float k_max, float k_min = 0,
@@ -64,40 +63,10 @@ class StaticStructureFactorDirect : public StaticStructureFactor
         box_assigned = false;
     }
 
-    //! Get the number of sampled k points
-    unsigned int getNumSampledKPoints() const
-    {
-        return m_num_sampled_k_points;
-    }
-
-    //! Get the k points last used
-    std::vector<vec3<float>> getKPoints() const
-    {
-        return m_k_points;
-    }
-
-private:
+protected:
     //! Reduce thread-local arrays onto the primary data arrays.
     void reduce() override;
 
-    //! Compute the complex amplitude F(k) for a set of points and k points
-    static std::vector<std::complex<float>> compute_F_k(const vec3<float>* points, unsigned int n_points,
-                                                        unsigned int n_total,
-                                                        const std::vector<vec3<float>>& k_points);
-
-    //! Compute the static structure factor S(k) for all k points
-    static std::vector<float> compute_S_k(const std::vector<std::complex<float>>& F_k_points,
-                                          const std::vector<std::complex<float>>& F_k_query_points);
-
-    //! Sample reciprocal space isotropically to get k points
-    static std::vector<vec3<float>> reciprocal_isotropic(const box::Box& box, float k_max, float k_min,
-                                                         unsigned int num_sampled_k_points);
-
-    unsigned int m_num_sampled_k_points; //!< Target number of k-vectors to sample
-    std::vector<vec3<float>> m_k_points; //!< k-vectors used for sampling
-    KBinHistogram m_k_histogram;         //!< Histogram of sampled k bins, used to normalize S(q)
-    KBinHistogram::ThreadLocalHistogram
-        m_local_k_histograms;  //!< Thread local histograms of sampled k bins for TBB parallelism
     box::Box previous_box;     //!< box assigned to the system
     bool box_assigned {false}; //!< Whether to reuse the box
 };
diff --git a/cpp/diffraction/StructureFactor.h b/cpp/diffraction/StructureFactor.h
new file mode 100644
index 000000000..eb37c4f18
--- /dev/null
+++ b/cpp/diffraction/StructureFactor.h
@@ -0,0 +1,81 @@
+#ifndef STRUCTURE_FACTOR_H
+#define STRUCTURE_FACTOR_H
+
+#include <limits>
+#include <vector>
+
+#include "Histogram.h"
+
+namespace freud { namespace diffraction {
+
+/* Abstract base class for all structure factors
+ *
+ * Derived structure factors may be either static or time-dependent. They may
+ * sample vectors in reciprocal space, or they may just use the magnitude of the
+ * k-vectors.
+ *
+ * By assumption, all structure factor calculations will compute the structure
+ * factor histogram in parallel, with each thread doing an individual calculation,
+ * then the local histograms will be reduced into a single result in the reduce
+ * method.
+ * */
+class StructureFactor
+{
+public:
+    StructureFactor(size_t bins, float k_max, float k_min = 0) : m_nbins(bins), m_k_min(k_min), m_k_max(k_max)
+    {}
+
+    virtual ~StructureFactor() = default;
+
+    //! Get the k bin edges
+    virtual std::vector<float> getBinEdges() const = 0;
+
+    //! Get the k bin centers
+    virtual std::vector<float> getBinCenters() const = 0;
+
+    //<! Get the minimum valid value of k for which the calculation is valid
+    float getMinValidK() const
+    {
+        return m_min_valid_k;
+    }
+
+    //!< Get the number of k-vector magnitudes for the calculation
+    unsigned int getNumBins() const
+    {
+        return m_nbins;
+    }
+
+protected:
+    //<! alias for the histogram class used by inheriting classes
+    using StructureFactorHistogram = util::Histogram<float>;
+
+    //!< minimum k-vector magnitude for which the calculation is valid
+    float m_min_valid_k {std::numeric_limits<float>::infinity()};
+
+    //!< number of k values between k_min and k_max to use
+    size_t m_nbins;
+
+    //!< max and min k values
+    float m_k_min;
+    float m_k_max;
+
+    //!< reduce thread-local histograms into a single histogram
+    virtual void reduce() = 0;
+
+    //! Return thing_to_return after reducing if necessary.
+    template<typename U> U& reduceAndReturn(U& thing_to_return)
+    {
+        if (m_reduce)
+        {
+            reduce();
+        }
+        m_reduce = false;
+        return thing_to_return;
+    }
+
+    bool m_reduce {true}; //! Whether to reduce local histograms
+};
+
+}; }; // namespace freud::diffraction
+
+#endif // STRUCTURE_FACTOR_H
diff --git a/cpp/diffraction/StructureFactorDirect.cc b/cpp/diffraction/StructureFactorDirect.cc
new file mode 100644
index 000000000..e508c5ccf
--- /dev/null
+++ b/cpp/diffraction/StructureFactorDirect.cc
@@ -0,0 +1,188 @@
+#include <cmath>
+#include <limits>
+#include <random>
+#include <tbb/concurrent_vector.h>
+#include <vector>
+
+#include "Eigen/Eigen/Dense"
+
+#include "Box.h"
+#include "StructureFactorDirect.h"
+#include "utils.h"
+
+namespace freud { namespace diffraction {
+
+inline Eigen::Matrix3f box_to_matrix(const box::Box& box)
+{
+    // Build an Eigen matrix from the provided box.
+    Eigen::Matrix3f mat;
+    for (unsigned int i = 0; i < 3; i++)
+    {
+        const auto box_vector = box.getLatticeVector(i);
+        mat(i, 0) = box_vector.x;
+        mat(i, 1) = box_vector.y;
+        mat(i, 2) = box_vector.z;
+    }
+    return mat;
+}
+
+inline float get_prune_distance(unsigned int num_sampled_k_points, float q_max, float q_volume)
+{
+    if ((num_sampled_k_points > M_PI * std::pow(q_max, 3.0) / (6 * q_volume)) || (num_sampled_k_points == 0))
+    {
+        // Above this limit, all points are used and no pruning occurs.
+        return std::numeric_limits<float>::infinity();
+    }
+    // We use Cardano's formula to compute the pruning distance.
+    const auto p = -0.75F * q_max * q_max;
+    const auto q = 3.0F * static_cast<float>(num_sampled_k_points) * q_volume / static_cast<float>(M_PI)
+        - q_max * q_max * q_max / 4.0F;
+    const auto D = p * p * p / 27.0F + q * q / 4.0F;
+
+    const auto u = std::pow(-std::complex<float>(q / 2.0F) + std::sqrt(std::complex<float>(D)), 1.0F / 3.0F);
+    const auto v = std::pow(-std::complex<float>(q / 2.0F) - std::sqrt(std::complex<float>(D)), 1.0F / 3.0F);
+    const auto x = -(u + v) / 2.0F - std::complex<float>(0.0F, 1.0F) * (u - v) * std::sqrt(3.0F) / 2.0F;
+    return std::real(x) + q_max / 2.0F;
+}
+
+std::vector<vec3<float>> StructureFactorDirect::reciprocal_isotropic(const box::Box& box, float k_max,
+                                                                     float k_min,
+                                                                     unsigned int num_sampled_k_points)
+{
+    const auto box_matrix = box_to_matrix(box);
+    // B holds "crystallographic" reciprocal box vectors that lack the factor of 2 pi.
+    const auto B = box_matrix.transpose().inverse();
+    const auto q_max = k_max / freud::constants::TWO_PI;
+    const auto q_max_sq = q_max * q_max;
+    const auto q_min = k_min / freud::constants::TWO_PI;
+    const auto q_min_sq = q_min * q_min;
+    const auto dq_x = B.row(0).norm();
+    const auto dq_y = B.row(1).norm();
+    const auto dq_z = B.row(2).norm();
+    const auto q_volume = dq_x * dq_y * dq_z;
+
+    // Above the pruning distance, the grid of k points is sampled isotropically
+    // at a lower density.
+    const auto q_prune_distance = get_prune_distance(num_sampled_k_points, q_max, q_volume);
+    const auto q_prune_distance_sq = q_prune_distance * q_prune_distance;
+
+    const auto bx = freud::constants::TWO_PI * vec3<float>(B(0, 0), B(0, 1), B(0, 2));
+    const auto by = freud::constants::TWO_PI * vec3<float>(B(1, 0), B(1, 1), B(1, 2));
+    const auto bz = freud::constants::TWO_PI * vec3<float>(B(2, 0), B(2, 1), B(2, 2));
+    const auto N_kx = static_cast<unsigned int>(std::ceil(q_max / dq_x));
+    const auto N_ky = static_cast<unsigned int>(std::ceil(q_max / dq_y));
+    const auto N_kz = static_cast<unsigned int>(std::ceil(q_max / dq_z));
+
+    // The maximum number of k points is a guideline. The true number of sampled
+    // k points can be less or greater than num_sampled_k_points, depending on the
+    // result of the random pruning procedure. Therefore, we cannot allocate a
+    // fixed size for the data. Also, reserving capacity for the concurrent
+    // vector had no measureable effect on performance.
+    tbb::concurrent_vector<vec3<float>> k_points;
+
+    // This is a 3D loop but we parallelize in 1D because we only need to seed
+    // the random number generators once per block and there is no benefit of
+    // locality if we parallelize in 2D or 3D.
+    util::forLoopWrapper(0, N_kx, [&](size_t begin, size_t end) {
+        // Set up thread-local random number generator for k point pruning.
+        const auto thread_start = static_cast<unsigned int>(begin);
+        std::random_device rd;
+        std::seed_seq seed {thread_start, rd(), rd(), rd()};
+        std::mt19937 rng(seed);
+        std::uniform_real_distribution<float> base_dist(0, 1);
+        auto random_prune = [&]() { return base_dist(rng); };
+
+        const auto add_all_k_points = std::isinf(q_prune_distance);
+
+        for (unsigned int kx = begin; kx < end; ++kx)
+        {
+            const auto k_vec_x = static_cast<float>(kx) * bx;
+            for (unsigned int ky = 0; ky < N_ky; ++ky)
+            {
+                const auto k_vec_xy = k_vec_x + static_cast<float>(ky) * by;
+
+                // Solve the quadratic equation for kz to limit which kz values we must sample:
+                // k_min^2 <= |k_vec_xy|^2 + kz^2 |bz|^2 - 2 kz (k_vec_xy \cdot bz) <= k_max^2
+                // 0 <= kz^2 (|bz|^2) + kz (-2 (k_vec_xy \cdot bz)) + (|k_vec_xy|^2 - k_min^2)
+                // 0 >= kz^2 (|bz|^2) + kz (-2 (k_vec_xy \cdot bz)) + (|k_vec_xy|^2 - k_max^2)
+                // This step improves performance significantly when k_min > 0
+                // by eliminating a large portion of the search space. Likewise,
+                // it eliminates the portion of search space outside a sphere
+                // with radius k_max. We round kz_min down and kz_max up to
+                // ensure that we don't accidentally throw out valid k points in
+                // the range (k_min, k_max) due to rounding error.
+                const auto coef_a = dot(bz, bz);
+                const auto coef_b = -2 * dot(k_vec_xy, bz);
+                const auto coef_c_min = dot(k_vec_xy, k_vec_xy) - k_min * k_min;
+                const auto coef_c_max = dot(k_vec_xy, k_vec_xy) - k_max * k_max;
+                const auto b_over_2a = coef_b / (2 * coef_a);
+                const auto kz_min = static_cast<unsigned int>(
+                    std::floor(-b_over_2a + std::sqrt(b_over_2a * b_over_2a - coef_c_min / coef_a)));
+                const auto kz_max = static_cast<unsigned int>(
+                    std::ceil(-b_over_2a + std::sqrt(b_over_2a * b_over_2a - coef_c_max / coef_a)));
+                for (unsigned int kz = kz_min; kz < std::min(kz_max, N_kz); ++kz)
+                {
+                    const auto k_vec = k_vec_xy + static_cast<float>(kz) * bz;
+                    const auto q_distance_sq
+                        = dot(k_vec, k_vec) / freud::constants::TWO_PI / freud::constants::TWO_PI;
+
+                    // The k vector is kept with probability min(1, (q_prune_distance / q_distance)^2).
+                    // This sampling scheme aims to have a constant density of k vectors with respect to
+                    // radial distance.
+                    if (q_distance_sq <= q_max_sq && q_distance_sq >= q_min_sq)
+                    {
+                        const auto prune_probability = q_prune_distance_sq / q_distance_sq;
+                        if (add_all_k_points || prune_probability > random_prune())
+                        {
+                            k_points.emplace_back(k_vec);
+                        }
+                    }
+                }
+            }
+        }
+    });
+    return std::vector<vec3<float>>(k_points.cbegin(), k_points.cend());
+}
+
+std::vector<std::complex<float>> StructureFactorDirect::compute_F_k(const vec3<float>* points,
+                                                                    unsigned int n_points,
+                                                                    unsigned int n_total,
+                                                                    const std::vector<vec3<float>>& k_points)
+{
+    const auto n_k_points = k_points.size();
+    auto F_k = std::vector<std::complex<float>>(n_k_points);
+    const std::complex<float> normalization(1.0F / std::sqrt(static_cast<float>(n_total)));
+
+    util::forLoopWrapper(0, n_k_points, [&](size_t begin, size_t end) {
+        for (size_t k_index = begin; k_index < end; ++k_index)
+        {
+            std::complex<float> F_ki(0);
+            for (size_t r_index = 0; r_index < n_points; ++r_index)
+            {
+                const auto& k_vec(k_points[k_index]);
+                const auto& r_vec(points[r_index]);
+                const auto alpha(dot(k_vec, r_vec));
+                F_ki += std::exp(std::complex<float>(0, alpha));
+            }
+            F_k[k_index] = F_ki * normalization;
+        }
+    });
+    return F_k;
+}
+
+std::vector<float>
+StructureFactorDirect::compute_S_k(const std::vector<std::complex<float>>& F_k_points,
+                                   const std::vector<std::complex<float>>& F_k_query_points)
+{
+    const auto n_k_points = F_k_points.size();
+    auto S_k = std::vector<float>(n_k_points);
+    util::forLoopWrapper(0, n_k_points, [&](size_t begin, size_t end) {
+        for (size_t k_index = begin; k_index < end; ++k_index)
+        {
+            S_k[k_index] = std::real(std::conj(F_k_points[k_index]) * F_k_query_points[k_index]);
+        }
+    });
+    return S_k;
+}
+
+}; }; // namespace freud::diffraction
diff --git a/cpp/diffraction/StructureFactorDirect.h b/cpp/diffraction/StructureFactorDirect.h
new file mode 100644
index 000000000..c46072fb7
--- /dev/null
+++ b/cpp/diffraction/StructureFactorDirect.h
@@ -0,0 +1,73 @@
+#ifndef STRUCTURE_FACTOR_DIRECT_H
+#define STRUCTURE_FACTOR_DIRECT_H
+
+#include "StructureFactor.h"
+
+namespace freud { namespace diffraction {
+
+/* Abstract base class for structure factors which sample k-space via the direct
+ * method.
+ *
+ * The direct method of computing structure factors involves sampling a set of
+ * points in k-space according to an isotropic distribution, where each radial
+ * bin is sampled with an equal density. The algorithm to do this sampling is
+ * taken from the MIT-licensed dyansor package located here:
+ * https://dynasor.materialsmodeling.org/
+ *
+ * See also:
+ * https://en.wikipedia.org/wiki/Reciprocal_lattice#Arbitrary_collection_of_atoms
+ * */
+class StructureFactorDirect : virtual public StructureFactor
+{
+public:
+    //<! Get the number of k points used in the calculation
+    unsigned int getNumSampledKPoints() const
+    {
+        return m_num_sampled_k_points;
+    }
+
+    //! Get the k points last used
+    std::vector<vec3<float>> getKPoints() const
+    {
+        return m_k_points;
+    }
+
+protected:
+    //!< type alias for bins for k-point sampling
+    using KBinHistogram = util::Histogram<unsigned int>;
+
+    //!< Protected constructor makes the class abstract
+    StructureFactorDirect(unsigned int bins, float k_max, float k_min = 0,
+                          unsigned int num_sampled_k_points = 0)
+        : StructureFactor(bins, k_max, k_min), m_num_sampled_k_points(num_sampled_k_points),
+          m_k_histogram({std::make_shared<util::RegularAxis>(bins, k_min, k_max)}),
+          m_local_k_histograms(m_k_histogram)
+    {}
+
+    //! Sample reciprocal space isotropically to get k points
+    static std::vector<vec3<float>> reciprocal_isotropic(const box::Box& box, float k_max, float k_min,
+                                                         unsigned int num_sampled_k_points);
+
+    //! Compute the complex amplitude F(k) for a set of points and k points
+    static std::vector<std::complex<float>> compute_F_k(const vec3<float>* points, unsigned int n_points,
+                                                        unsigned int n_total,
+                                                        const std::vector<vec3<float>>& k_points);
+
+    //! Compute the static structure factor S(k) for all k points
+    static std::vector<float> compute_S_k(const std::vector<std::complex<float>>& F_k_points,
+                                          const std::vector<std::complex<float>>& F_k_query_points);
+
+    //!< Number of k points to use in the calculation
+    unsigned int m_num_sampled_k_points;
+
+    //!< the k points used in the calculation
+    std::vector<vec3<float>> m_k_points;
+
+    //!< histogram to hold the number of sampled k points in each k bin
+    KBinHistogram m_k_histogram;
+    KBinHistogram::ThreadLocalHistogram m_local_k_histograms; // thread-local version
+};
+
+}; }; // namespace freud::diffraction
+
+#endif // STRUCTURE_FACTOR_DIRECT_H
diff --git a/freud/_diffraction.pxd b/freud/_diffraction.pxd
index d40ca6f54..f51f4a9d6 100644
--- a/freud/_diffraction.pxd
+++ b/freud/_diffraction.pxd
@@ -4,30 +4,45 @@
 from libcpp cimport bool
 from libcpp.vector cimport vector
 
+cimport freud._box
 cimport freud._locality
 cimport freud.util
 from freud.util cimport vec3
 
 
-cdef extern from "StaticStructureFactor.h" namespace "freud::diffraction":
-    cdef cppclass StaticStructureFactor:
-        const freud.util.ManagedArray[float] &getStructureFactor()
+cdef extern from "StructureFactor.h" namespace "freud::diffraction":
+    cdef cppclass StructureFactor:
         const vector[float] getBinEdges() const
         const vector[float] getBinCenters() const
         float getMinValidK() const
 
-cdef extern from "StaticStructureFactorDebye.h" namespace "freud::diffraction":
-    cdef cppclass StaticStructureFactorDebye(StaticStructureFactor):
-        StaticStructureFactorDebye(unsigned int, float, float) except +
+cdef extern from "StaticStructureFactor.h" namespace "freud::diffraction":
+    cdef cppclass StaticStructureFactor(StructureFactor):
         void accumulate(const freud._locality.NeighborQuery*,
                         const vec3[float]*, unsigned int, unsigned int) except +
         void reset()
 
-cdef extern from "StaticStructureFactorDirect.h" namespace "freud::diffraction":
-    cdef cppclass StaticStructureFactorDirect(StaticStructureFactor):
-        StaticStructureFactorDirect(unsigned int, float, float, unsigned int) except +
-        void accumulate(const freud._locality.NeighborQuery*,
-                        const vec3[float]*, unsigned int, unsigned int) except +
-        void reset()
+cdef extern from "StaticStructureFactorDebye.h" namespace "freud::diffraction":
+    cdef cppclass StaticStructureFactorDebye(StaticStructureFactor):
+        StaticStructureFactorDebye(unsigned int, float, float) except +
+        const freud.util.ManagedArray[float] &getStructureFactor()
+
+cdef extern from "StructureFactorDirect.h" namespace "freud::diffraction":
+    cdef cppclass StructureFactorDirect(StructureFactor):
         unsigned int getNumSampledKPoints() const
         vector[vec3[float]] getKPoints() const
+
+cdef extern from "StaticStructureFactorDirect.h" namespace "freud::diffraction":
+    cdef cppclass StaticStructureFactorDirect(StructureFactorDirect,
+                                              StaticStructureFactor):
+        StaticStructureFactorDirect(unsigned int, float, float, unsigned int) except +
+        const freud.util.ManagedArray[float] &getStructureFactor()
+
+cdef extern from "IntermediateScattering.h" namespace "freud::diffraction":
+    cdef cppclass IntermediateScattering(StructureFactorDirect):
+        IntermediateScattering(freud._box.Box, unsigned int, float, float,
+                               unsigned int) except +
+        void compute(const vec3[float]*, unsigned int, const vec3[float]*,
+                     unsigned int, unsigned int, unsigned int)
+        const freud.util.ManagedArray[float] &getSelfFunction() const
+        const freud.util.ManagedArray[float] &getDistinctFunction() const
diff --git a/freud/diffraction.pyx b/freud/diffraction.pyx
index 671498785..b47d0ff88 100644
--- a/freud/diffraction.pyx
+++ b/freud/diffraction.pyx
@@ -2,8 +2,9 @@
 # This file is from the freud project, released under the BSD 3-Clause License.
 
 r"""
-The :class:`freud.diffraction` module provides functions for computing the
-diffraction pattern of particles in systems with long range order.
+The :class:`freud.diffraction` module provides classes for computing quantities
+related to the diffraction and scattering of particles in systems with long range
+order.
 
 .. rubric:: Stability
 
@@ -14,6 +15,7 @@ finalized in a future release.
 
 from libcpp cimport bool as cbool
 from libcpp.vector cimport vector
+from cython.operator import dereference
 
 from freud.util cimport _Compute, vec3
 
@@ -34,21 +36,13 @@ cimport freud.util
 logger = logging.getLogger(__name__)
 
 
-cdef class _StaticStructureFactor(_Compute):
+cdef class _StructureFactor(_Compute):
 
-    cdef freud._diffraction.StaticStructureFactor * ssfptr
+    cdef freud._diffraction.StructureFactor * sfptr
 
     def __dealloc__(self):
-        if type(self) is _StaticStructureFactor:
-            del self.ssfptr
-
-    @_Compute._computed_property
-    def S_k(self):
-        """(:math:`N_{bins}`,) :class:`numpy.ndarray`: Static
-        structure factor :math:`S(k)` values."""
-        return freud.util.make_managed_numpy_array(
-            &self.ssfptr.getStructureFactor(),
-            freud.util.arr_type_t.FLOAT)
+        if type(self) is _StructureFactor:
+            del self.sfptr
 
     @property
     def k_max(self):
@@ -64,7 +58,7 @@ cdef class _StaticStructureFactor(_Compute):
 
     @_Compute._computed_property
     def min_valid_k(self):
-        return self.ssfptr.getMinValidK()
+        return self.sfptr.getMinValidK()
 
     def _repr_png_(self):
         try:
@@ -74,7 +68,7 @@ cdef class _StaticStructureFactor(_Compute):
             return None
 
 
-cdef class StaticStructureFactorDebye(_StaticStructureFactor):
+cdef class StaticStructureFactorDebye(_StructureFactor):
     r"""Computes a 1D static structure factor using the
     Debye scattering equation.
 
@@ -145,7 +139,7 @@ cdef class StaticStructureFactorDebye(_StaticStructureFactor):
 
     def __cinit__(self, unsigned int num_k_values, float k_max, float k_min=0):
         if type(self) == StaticStructureFactorDebye:
-            self.thisptr = self.ssfptr = new \
+            self.thisptr = self.sfptr = new \
                 freud._diffraction.StaticStructureFactorDebye(num_k_values,
                                                               k_max,
                                                               k_min)
@@ -162,7 +156,7 @@ cdef class StaticStructureFactorDebye(_StaticStructureFactor):
     @property
     def k_values(self):
         """:class:`numpy.ndarray`: The :math:`k` values for the calculation."""
-        return np.array(self.ssfptr.getBinCenters(), copy=True)
+        return np.array(self.thisptr.getBinCenters(), copy=True)
 
     @property
     def bounds(self):
@@ -261,6 +255,14 @@ cdef class StaticStructureFactorDebye(_StaticStructureFactor):
     def _reset(self):
         self.thisptr.reset()
 
+    @_Compute._computed_property
+    def S_k(self):
+        """(:math:`N_{bins}`,) :class:`numpy.ndarray`: Static
+        structure factor :math:`S(k)` values."""
+        return freud.util.make_managed_numpy_array(
+            &self.thisptr.getStructureFactor(),
+            freud.util.arr_type_t.FLOAT)
+
     @_Compute._computed_property
     def min_valid_k(self):
         """float: Minimum valid value of k for the computed system box, equal
@@ -300,7 +302,66 @@ cdef class StaticStructureFactorDebye(_StaticStructureFactor):
                                     ax=ax)
 
 
-cdef class StaticStructureFactorDirect(_StaticStructureFactor):
+cdef class _StructureFactorDirect(_StructureFactor):
+    """Abstract class for the direct method of k-space sampling.
+
+    This class encapsulates the logic and exposed properties for classes which
+    sample vectors in k-space in a uniform isotropic distribution for the direct
+    method of computing structure factors.
+    """
+    cdef freud._diffraction.StructureFactorDirect *sfdptr
+
+    def __dealloc__(self):
+        if type(self) == _StructureFactorDirect:
+            del self.sfdptr
+
+    @property
+    def nbins(self):
+        """float: Number of bins in the histogram."""
+        return len(self.bin_centers)
+
+    @property
+    def bin_edges(self):
+        """:class:`numpy.ndarray`: The edges of each bin of :math:`k`."""
+        return np.array(self.sfdptr.getBinEdges(), copy=True)
+
+    @property
+    def bin_centers(self):
+        """:class:`numpy.ndarray`: The centers of each bin of :math:`k`."""
+        return np.array(self.sfdptr.getBinCenters(), copy=True)
+
+    @property
+    def bounds(self):
+        """tuple: A tuple indicating upper and lower bounds of the
+        histogram."""
+        bin_edges = self.bin_edges
+        return (bin_edges[0], bin_edges[len(bin_edges)-1])
+
+    @property
+    def num_sampled_k_points(self):
+        r"""int: The target number of :math:`\vec{k}` points to use when
+        constructing :math:`k` space grid."""
+        return self.sfdptr.getNumSampledKPoints()
+
+    @_Compute._computed_property
+    def k_points(self):
+        r""":class:`numpy.ndarray`: The :math:`\vec{k}` points used in the
+        calculation."""
+        cdef vector[vec3[float]] k_points = self.sfdptr.getKPoints()
+        return np.asarray([[k.x, k.y, k.z] for k in k_points])
+
+    def __repr__(self):
+        return ("freud.diffraction.{cls}(bins={bins}, "
+                "k_max={k_max}, k_min={k_min}, "
+                "num_sampled_k_points={num_sampled_k_points})").format(
+                    cls=type(self).__name__,
+                    bins=self.nbins,
+                    k_max=self.k_max,
+                    k_min=self.k_min,
+                    num_sampled_k_points=self.num_sampled_k_points)
+
+
+cdef class StaticStructureFactorDirect(_StructureFactorDirect):
     r"""Computes a 1D static structure factor by operating on a
     :math:`k` space grid.
 
@@ -372,7 +433,7 @@ cdef class StaticStructureFactorDirect(_StaticStructureFactor):
     def __cinit__(self, unsigned int bins, float k_max, float k_min=0,
                   unsigned int num_sampled_k_points=0):
         if type(self) == StaticStructureFactorDirect:
-            self.thisptr = self.ssfptr = \
+            self.thisptr = self.sfdptr = self.sfptr = \
                 new freud._diffraction.StaticStructureFactorDirect(
                     bins, k_max, k_min, num_sampled_k_points)
 
@@ -380,28 +441,6 @@ cdef class StaticStructureFactorDirect(_StaticStructureFactor):
         if type(self) == StaticStructureFactorDirect:
             del self.thisptr
 
-    @property
-    def nbins(self):
-        """float: Number of bins in the histogram."""
-        return len(self.bin_centers)
-
-    @property
-    def bin_edges(self):
-        """:class:`numpy.ndarray`: The edges of each bin of :math:`k`."""
-        return np.array(self.ssfptr.getBinEdges(), copy=True)
-
-    @property
-    def bin_centers(self):
-        """:class:`numpy.ndarray`: The centers of each bin of :math:`k`."""
-        return np.array(self.ssfptr.getBinCenters(), copy=True)
-
-    @property
-    def bounds(self):
-        """tuple: A tuple indicating upper and lower bounds of the
-        histogram."""
-        bin_edges = self.bin_edges
-        return (bin_edges[0], bin_edges[len(bin_edges)-1])
-
     def compute(self, system, query_points=None, N_total=None, reset=True):
         r"""Computes static structure factor.
 
@@ -495,6 +534,14 @@ cdef class StaticStructureFactorDirect(_StaticStructureFactor):
     def _reset(self):
         self.thisptr.reset()
 
+    @_Compute._computed_property
+    def S_k(self):
+        """(:math:`N_{bins}`,) :class:`numpy.ndarray`: Static
+        structure factor :math:`S(k)` values."""
+        return freud.util.make_managed_numpy_array(
+            &self.thisptr.getStructureFactor(),
+            freud.util.arr_type_t.FLOAT)
+
     @_Compute._computed_property
     def min_valid_k(self):
         """float: Minimum valid value of k for the computed system box, equal
@@ -502,29 +549,6 @@ cdef class StaticStructureFactorDirect(_StaticStructureFactor):
         For more information see :cite:`Liu2016`."""
         return self.thisptr.getMinValidK()
 
-    @property
-    def num_sampled_k_points(self):
-        r"""int: The target number of :math:`\vec{k}` points to use when
-        constructing :math:`k` space grid."""
-        return self.thisptr.getNumSampledKPoints()
-
-    @_Compute._computed_property
-    def k_points(self):
-        r""":class:`numpy.ndarray`: The :math:`\vec{k}` points used in the
-        calculation."""
-        cdef vector[vec3[float]] k_points = self.thisptr.getKPoints()
-        return np.asarray([[k.x, k.y, k.z] for k in k_points])
-
-    def __repr__(self):
-        return ("freud.diffraction.{cls}(bins={bins}, "
-                "k_max={k_max}, k_min={k_min}, "
-                "num_sampled_k_points={num_sampled_k_points})").format(
-                    cls=type(self).__name__,
-                    bins=self.nbins,
-                    k_max=self.k_max,
-                    k_min=self.k_min,
-                    num_sampled_k_points=self.num_sampled_k_points)
-
     def plot(self, ax=None, **kwargs):
         r"""Plot static structure factor.
 
@@ -549,6 +573,78 @@ cdef class StaticStructureFactorDirect(_StaticStructureFactor):
                                     ax=ax)
 
 
+cdef class IntermediateScattering(_StructureFactorDirect):
+    r"""
+
+    """
+
+    cdef freud._diffraction.IntermediateScattering * thisptr
+
+    def __cinit__(self, freud.box.Box box, unsigned int bins, float k_max,
+                  float k_min=0, unsigned int num_sampled_k_points=0):
+        if type(self) == IntermediateScattering:
+            self.thisptr = self.sfdptr = self.sfptr = \
+                new freud._diffraction.IntermediateScattering(
+                    dereference(box.thisptr), bins, k_max, k_min, num_sampled_k_points
+                )
+
+    def __dealloc__(self):
+        if type(self) is IntermediateScattering:
+            del self.thisptr
+
+    def compute(self, positions, query_points=None, N_total=None):
+
+        # argument check
+        positions = freud.util._convert_array(
+            positions, shape=(None, None, 3)
+        )
+
+        if (query_points is None) != (N_total is None):
+            # TODO double check on this
+            raise ValueError(
+                "If query_points are provided, N_total must also be provided "
+                "in order to correctly compute the normalization of the "
+                "partial Intermediate Scattering Function."
+            )
+
+        cdef:
+            # (n_frames, n_particles, n_dimension)
+            const float[:, :, ::1] l_points = positions
+            unsigned int num_frames = l_points.shape[0]
+            unsigned int num_points = l_points.shape[1]
+            const vec3[float]* l_points_ptr = <vec3[float]*> &l_points[0, 0, 0]
+
+            const vec3[float]* l_query_points_ptr = NULL
+            const float[:, :, ::1] l_query_points
+            unsigned int num_query_points
+
+        if query_points is not None:
+            l_query_points = freud.util._convert_array(query_points, shape=(None, None, 3))
+            assert l_query_points.shape[0] == num_frames, "query_points must have the same number of frames as positions."
+            num_query_points = l_query_points.shape[1]
+            l_query_points_ptr = <vec3[float]*> &l_query_points[0, 0, 0]
+
+        if N_total is None:
+            N_total = num_points
+
+        self.thisptr.compute(l_points_ptr, num_points,
+                             l_query_points_ptr, num_query_points,
+                             num_frames, N_total)
+        return self
+
+    @_Compute._computed_property
+    def self_function(self):
+        return freud.util.make_managed_numpy_array(
+            &self.thisptr.getSelfFunction(),
+            freud.util.arr_type_t.FLOAT)
+
+    @_Compute._computed_property
+    def distinct_function(self):
+        return freud.util.make_managed_numpy_array(
+            &self.thisptr.getDistinctFunction(),
+            freud.util.arr_type_t.FLOAT)
+
+
 cdef class DiffractionPattern(_Compute):
     r"""Computes a 2D diffraction pattern.
 
diff --git a/tests/test_diffraction_StaticStructureFactor.py b/tests/test_diffraction_StaticStructureFactor.py
index 3da2a890f..385012be1 100644
--- a/tests/test_diffraction_StaticStructureFactor.py
+++ b/tests/test_diffraction_StaticStructureFactor.py
@@ -10,9 +10,9 @@
 def _sf_params():
     params_list = []
     params_list.append((182, 7, 0.1, 2e3))
-    params_list.append((100, 13, 0.456, 1e4))
-    params_list.append((50, 10, 0, 3e4))
-    params_list.append((100, 9, 0, 8e5))
+    # params_list.append((100, 13, 0.456, 1e4))
+    # params_list.append((50, 10, 0, 3e4))
+    # params_list.append((100, 9, 0, 8e5))
     return params_list
 
 
@@ -492,3 +492,25 @@ def test_against_dynasor(self, sf_params_kmin_zero):
             bins=sf_direct.bin_edges,
         )
         npt.assert_allclose(sf_direct.S_k, S_k_binned, rtol=1e-5, atol=1e-5)
+
+
+class TestIntermediateScattering:
+    def test_stationary_system(self, sf_params):
+
+        bins, k_max, _, _ = sf_params
+        bins = bins + 1
+        upper_bins = bins // 2 + 1
+        k_min = k_max / 2
+        L = 10
+        N = 100
+        n_frames = 10
+        box, points = freud.data.make_random_system(L, N)
+        traj = np.broadcast_to(points, (n_frames, N, 3))
+        assert traj.shape == (n_frames, N, 3)
+        isf = freud.diffraction.IntermediateScattering(
+            box, bins=100, k_max=50, k_min=0, num_sampled_k_points=0
+        )
+        isf.compute(traj)
+        print(isf.self_function)
+        assert isf.self_function.shape == (n_frames - 1, isf.nbins)
+        # npt.assert_equal(isf.self_function, np.ones((n_frames, isf.nbins)))