From a07ada3bca7e1894b6cbc9aa2574f4b42154e599 Mon Sep 17 00:00:00 2001
From: Chi Cheng <chicheng85@gmail.com>
Date: Mon, 24 Feb 2025 17:30:00 +0000
Subject: [PATCH] updated correlation documentation, added notation section and
 spectrum code comment

---
 Doc/index.rst                           |  1 +
 Doc/pages/R_notation.rst                | 19 +++++++
 Doc/pages/correlation.rst               | 68 ++++++++++++-------------
 MDANSE/Src/MDANSE/Mathematics/Signal.py |  7 +--
 4 files changed, 56 insertions(+), 39 deletions(-)
 create mode 100644 Doc/pages/R_notation.rst

diff --git a/Doc/index.rst b/Doc/index.rst
index 0e5122dd3..60526b417 100644
--- a/Doc/index.rst
+++ b/Doc/index.rst
@@ -93,5 +93,6 @@ J. Chem. Inf. Model. 2017, 57, 1, 1–5 <https://doi.org/10.1021/acs.jcim.6b0057
    pages/parameters
    pages/R_traj
    pages/R_units
+   pages/R_notation
    pages/R_further
    pages/references
diff --git a/Doc/pages/R_notation.rst b/Doc/pages/R_notation.rst
new file mode 100644
index 000000000..5e701159b
--- /dev/null
+++ b/Doc/pages/R_notation.rst
@@ -0,0 +1,19 @@
+
+.. _notation:
+
+Frequently used symbols and notations
+=====================================
+
+======================== ===============
+**Symbol**               **Description**
+:math:`t`                time
+:math:`t_{\mathrm{cor}}` the total time of the correlation function
+:math:`t_{\mathrm{tot}}` the total time of the MD simulation
+:math:`\Delta t`         the MD time step
+:math:`T`                temperature
+:math:`n_{\mathrm{c}}`   total number of time steps calculated for the correlation function
+:math:`n_{\mathrm{t}}`   total number of time steps in the MD simulation
+:math:`N`                total number of atoms in the unit cell, simulation box or volume :math:`V`
+:math:`\omega`           frequency
+:math:`\Delta \omega`    the frequency step
+======================== ===============
diff --git a/Doc/pages/correlation.rst b/Doc/pages/correlation.rst
index 7b8f16e3d..38b8d6e3b 100644
--- a/Doc/pages/correlation.rst
+++ b/Doc/pages/correlation.rst
@@ -6,65 +6,61 @@ Correlation Functions
 
 Most of the quantities which can be extracted from MD
 simulations are time correlation functions. In MDANSE we use a correlation
-window method to ensure that the time averaging for each time step
+window to ensure that the time averaging for each time step
 is done in a consistent way. Consider two time series
 
 .. math::
    :label: eqn-fca1
 
-   a(k\Delta t)
-   \qquad b(k\Delta t)
-   \qquad k = 0, \ldots, N_t-1,
+   A(k \Delta t) \quad \text{and} \quad B(k \Delta t) \qquad k = 0, \ldots, n_{\mathrm{t}}-1,
 
-of length :math:`T = (N_t -1)\Delta t` which are
-to be correlated. The following the shorthands
-:math:`a(k)` and :math:`b(k)` will be used. Now we choose a specific
-number of correlation frames :math:`N_c` to use which will define
-the length of our correlation function. The correlation function of
-:math:`a(k)` and :math:`b(k)` will be
+of length :math:`t_{\mathrm{tot}} = (n_{\mathrm{t}} -1) \Delta t` which are
+to be correlated. In MDANSE,
+correlation function are calculated by first choosing a specific
+number of correlation time steps :math:`n_{\mathrm{c}}` which will define
+the length of our correlation function :math:`t_{\mathrm{cor}} = (n_{\mathrm{c}} -1) \Delta t`. The correlation function of
+:math:`A(k \Delta t)` and :math:`B(k \Delta t)` will be
 
 .. math::
    :label: eqn-fca2
 
-   c_{ab}(m) = \frac{1}{N_t - N_c + 1} \sum\limits_{k=0}^{N_t - N_c + 1} a^{*}(k)b(k + m) \qquad m = 0, \ldots, N_{c} - 1.
+   C_{AB}(l \Delta t) = \frac{1}{n_{\mathrm{t}} - n_{\mathrm{c}} + 1} \sum\limits_{k=0}^{n_{\mathrm{t}} - n_{\mathrm{c}} + 1} A^{*}(k\Delta t)B([k + l]\Delta t) \qquad l = 0, \ldots, n_{c} - 1.
 
-In case that :math:`a(k)` and
-:math:`b(k)` are identical, the corresponding correlation function
-:math:`c_{aa}(m)` is called an *autocorrelation* function. Notice that
-the prefactor is the same for all time steps :math:`m`, in previous
-versions of MDANSE this was not the case. This meant that for different
+In case that :math:`A(k \Delta t)` and
+:math:`B(k \Delta t)` are identical, the corresponding correlation function
+:math:`C_{AA}(l \Delta t)` is called an *autocorrelation* function. Notice that
+the prefactor is the same for all :math:`l \Delta t` time steps, this was
+not the case in previous versions of MDANSE. This meant that for different
 time steps a different number of configurations were used to obtain the
 average correlation; leading to spuriously large correlations for some
-time intervals.
+time intervals. However in MDANSE 2 your correlation functions will be
+truncated since :math:`t_{\mathrm{tot}} \geq t_{\mathrm{cor}}`.
 
 In many cases not only is the computation of a correlation function
 required, but also the computation of its Fourier spectrum. In
-MDANSE the spectra can be smoothed by applying a window in the time
-domain [Ref45]_ 
+MDANSE the spectra can be smoothed by applying a instrument resolution
+function
 
 .. math::
    :label: eqn-fca11
 
-   P_{ab}\left(\frac{n}{2N_c}\right) =
-   \Delta t \sum_{m=-(N_c-1)}^{N_c-1}
-   \exp\left[-2\pi i\left(\frac{n m}{2N_c}\right)\right]
-   \,W(m)\,\frac{1}{N_c-|m|}c_{ab}(m).
+   P_{AB}\left(m \Delta \omega \right) = \frac{\Delta t}{2 \pi}\sum_{l=-(n_{\mathrm{c}}-1)}^{n_{\mathrm{c}}-1}
+   \exp\left[- 2 \pi i \frac{l \Delta t }{2n_{\mathrm{c}} - 1} m \Delta \omega \right] \frac{W(l \Delta t)}{W(0)} C_{AB}( \vert l \Delta t \vert )
 
-The time step :math:`\Delta t` in front of the sum yields the proper
-normalization of the spectrum. For example, a Gaussian window
-[Ref46]_ can be used where:
+here :math:`m = -(n_{\mathrm{c}}-1), \ldots, n_{\mathrm{c}}-1` and :math:`\Delta \omega` is the
+frequency step. Notice that the we assume that the correlation is symmetric so
+that :math:`C_{AB}(l \Delta t) = C_{AB}( |l \Delta t| )` which should
+approximately be the case for all the correlation functions calculated
+in MDANSE assuming good (equilibrated, of a sufficient length/size and
+etc) MD trajectories are used. A Gaussian window is one example
+instrument resolution function that can be used can be used where:
 
 .. math::
    :label: eqn-fca12
 
-   W(m) = \exp\left[
-   -\frac{1}{2}\left(\alpha\frac{|m|}{N_c-1}\right)^2
+   W(l \Delta t ) = \frac{1}{2n_{\mathrm{c}} - 1} \frac{\sqrt{2 \pi}}{ \sigma \Delta t} \sum_{m=-(n_{\mathrm{c}}-1)}^{n_{\mathrm{c}}-1} \exp\left[ 2 \pi i \frac{m \Delta \omega}{2n_{\mathrm{c}} - 1} l \Delta t \right] \exp\left[
+   -\frac{1}{2}\left(\frac{ m \Delta \omega }{\sigma}\right)^2
    \right]
-   \qquad m = -(N_c-1), \ldots, N_c-1.
-
-Its widths in the time and frequency domain are :math:`\sigma_t = \alpha/T`
-and :math:`\sigma_\nu = 1/(2\pi\sigma_t)`, respectively.
-:math:`\sigma_\nu` corresponds to the width of the resolution
-function of the Fourier spectrum and
-:math:`T_{c} =(N_c-1)\Delta t` is the length of the correlation
-function.
+
+here :math:`l = -(n_{\mathrm{c}}-1), \ldots, n_{\mathrm{c}}-1` and :math:`\sigma` corresponds to the width of the resolution
+function of the Fourier spectrum.
diff --git a/MDANSE/Src/MDANSE/Mathematics/Signal.py b/MDANSE/Src/MDANSE/Mathematics/Signal.py
index 9d53ef6d0..51851c7eb 100644
--- a/MDANSE/Src/MDANSE/Mathematics/Signal.py
+++ b/MDANSE/Src/MDANSE/Mathematics/Signal.py
@@ -279,9 +279,10 @@ def get_spectrum(signal, window=None, timeStep=1.0, axis=0, fft="fft"):
 
     s = tuple(s)
 
-    # We compute the non-unitary fourier transform with angular
-    # frequencies with a 1/2pi factor applied to the forward transform.
-    # This is done for some historical reason see the git history.
+    # We compute the non-unitary fourier transform with a 1/2pi factor
+    # applied to the forward transform and angular frequencies.
+    # See the derivation of S(q,w) from QM in Principles of Neutron
+    # Scattering from Condensed Matter, Chap. 3.
 
     # For information about the manipulation around fftshift and ifftshift
     # http://www.mathworks.com/matlabcentral/newsreader/view_thread/285244