FFT Tutorial with discretisedfield

This Jupyter notebook is a tutorial guide for using the Fast Fourier Transform (FFT) functionality in discretisedfield. FFT is a highly efficient algorithm used to decompose a dataset into its base frequencies. In discretisedfield, our primary focus is on spatial transformations.

In this tutorial we explain the different methods to perform a FFT on a field by looking at a few examples. We will start with a simple scalar field in a region with two spatial dimensions and move on to give a more complex example of a vector field in three spatial dimensions.

Setup field

To get started, let’s import the discretisedfield package.

import discretisedfield as df

As an example we study a rectangular geometry with a sinusoidal wave propagating along the y direction with a \(9\, \mathrm{m}\) wavelength. (For this first example, we use a region and wavelength on the metre scale to avoid large negative/positive multipliers in real/Fourier space)

import numpy as np

mesh = df.Mesh(p1=(0, 0), p2=(20, 45), cell=(1, 1))


def value(p):
    _, y = p  # Ignoring the x-coordinate
    k = 1 / 9  # 1/m
    return np.sin(2 * np.pi * k * y)


field = df.Field(mesh, nvdim=1, value=value)
field

Field

• Mesh
  • Region
    • pmin = [0, 0]
    • pmax = [20, 45]
    • dims = ['x', 'y']
    • units = ['m', 'm']
  • n = [20, 45]
• nvdim = 1

field.mpl()

FFT Methods

The discretisedfield package provides four different methods for performing Fast Fourier Transform (FFT) operations on a discretisedfield.Field object:

1. fftn: This is the “full” FFT method which computes the n-dimensional FFT. It returns the complex-valued FFT result that has the same shape as the input field. It is used when the full frequency spectrum is needed.

2. ifftn: This is the inverse of the “full” FFT method, which takes a complex-valued FFT result and computes the inverse FFT to return a field in the spatial domain. It is used to go back from the spatial frequency domain to the spatial domain.

3. rfftn: This stands for “real FFT”. It is similar to fftn, but it only computes the non-negative frequency terms along the last dimension. The result is a complex-valued array that has the same shape as the input array but half the size of the last dimension compared to the input. It is used when the input field only has real values, and hence its frequency spectrum is symmetric. Thus, it saves computation time and storage.

4. irfftn: This is the inverse of the “real FFT” method. It takes a half-sized complex-valued FFT result (as output by rfftn) and computes the inverse FFT to return a real-valued field in the spatial domain. It is used when the original field is known to be real-valued.

All of these methods are n dimensional and work with any number of value (vector) and mesh (geometric) dimensions. These four methods provide the necessary tools to perform FFT operations and allow the transition between the spatial domain and frequency domain, covering both real and complex-valued fields.

We use scipy.fft for the calculations. Please refer to the SciPy documentation for further details.

FFT

Now we can perform out first FFT.

fft_field = field.fftn()
fft_field

Field

• Mesh
  • Region
    • pmin = [-0.525, -0.5]
    • pmax = [0.47500000000000003, 0.5]
    • dims = ['k_x', 'k_y']
    • units = ['(m)$^{-1}$', '(m)$^{-1}$']
  • n = [20, 45]
• nvdim = 1

With the application of the FFT methods in discretisedfield, notable transformations occur in both the mesh and the field:

1. Mesh Transformation: The mesh is transposed from the regular space to the spatial frequency domain. As a consequence, the units transition from meters (m) to inverse meters \(\left(\mathrm{m}^{-1}\right)\). Furthermore, to signify this transformation, k_ is prepended to each of the dims of the mesh. However, the Fourier transformed mesh maintains the same number of points (n) as the original mesh.

   NOTE: The values of the axis in Fourier space are frequency NOT angular frequency so DO NOT include a factor of \(2\pi\), i.e. \(|{\bf k}| = \frac{1}{\lambda} \neq \frac{2 \pi}{\lambda}\), where \(\bf k\) is the wave vector and \(\lambda\) is the wavelength.

2. Field Transformation: The field, although it continues to possess three value dimensions, undergoes a transformation, highlighted by the prefix ft_ added to each dimension. This indicates that these dimensions have been processed through a Fourier Transform.

These changes not only assist in distinguishing the transformed quantities but also serve a practical purpose by providing the values in appropriate units and dimensions.

To visualise this Fourier transform we can plot the Power Spectral Density given by:

\[PSD(k) = | F(k) |^2\]

where:

• \(PSD(k)\) represents the power spectral density at spatial frequency \(k\).

• \(F(k)\) represents the Fourier transform of the original signal at spatial frequency \(k\).

• \(|\) \(|\) denotes the absolute value.

In essence, the PSD is the magnitude squared of the Fourier transform of the original signal, providing a measure of the power (or intensity) of the signal across the different spatial frequencies.

psd = np.power(np.abs(fft_field), 2)

psd.mpl()

In this plot, we observe prominent peaks at \(k_y = \pm \frac{1}{9}\,\mathrm{m}^{-1}= \pm 0.11\,\mathrm{m}^{-1}\) aligning with our expectations for a sinusoidal modulation with a wavelength of \(7\,\mathrm{m}\) and confirms the presence of the anticipated spatial frequency in our data.

iFFT

Now we have a Fourier transformed field we can transform it back into real space by using an ifft.

ifft_field = fft_field.ifftn()
ifft_field

Field

• Mesh
  • Region
    • pmin = [-10.0, -22.5]
    • pmax = [10.0, 22.5]
    • dims = ['x', 'y']
    • units = ['m', 'm']
  • n = [20, 45]
• nvdim = 1

Upon transforming our field back into the regular space, we observe several changes. Firstly, our units return to the meter (m) notation. Secondly, both the dimensions and the value dimensions revert to their original identifiers, removing the k_ and ft_ prefixes, respectively. This reflects the transition from the frequency to the spatial domain. As expected, the total number of data points within the field remains constant.

However, one key aspect of Fourier transforms to be aware of is the loss of the original translation during the transformation process. Consequently, the Fourier transformed field or signal is inherently centred around the origin.

This characteristic is due to the Fourier transform operating in the frequency domain, which identifies the frequencies present in the signal but disregards the location or shift of those frequencies in the original signal. Hence, after the Fourier transform, the information about the original position of the frequencies in the spatial domain is lost and the resulting transformed signal is always centred around the origin.

However, we can translate the mesh in order to move the field back to the desired position inplace. In the translate function we provide a vector by which we want to translate our region. Due to the application of the Fourier transforms, the mesh of the new field is centred around the origin. In order to align the new mesh with the original mesh we can provide the center of the original mesh as the vector.

ifft_field.mesh.translate(field.mesh.region.centre, inplace=True)

Mesh

• Region
  • pmin = [0.0, 0.0]
  • pmax = [20.0, 45.0]
  • dims = ['x', 'y']
  • units = ['m', 'm']
• n = [20, 45]

As a check we can see if our original field is the same as the field which has been Fourier transformed and then inverse Fourier transformed. By using the allclose method in the Field class we can see that the two field are almost the same.

field.allclose(ifft_field)

True

While our field has been transformed back into its original space, using the == operator reveals that the initial and final states are not identical.

field == ifft_field

False

This discrepancy stems primarily from two factors:

1. Floating-Point Precision: Arithmetic operations on floating-point numbers can lead to small imprecisions.

2. Complex Values: Initially, our field comprised real values. However, the processing via FFT and inverse FFT results in a field with complex values.

These factors contribute to the subtle differences when comparing the original and transformed fields directly and shows how the allclose method can offer a more suitable comparison.

We can visualise these differences by plotting the real and imaginary parts.

ifft_field.real.mpl()

ifft_field.imag.mpl()

As we can see, the real part of the inverse Fourier transformed field is visually identical to the original field, while the imaginary part is practically zero.

rFFT

As the original field we are using is real we can use a rfft.

rfft_field = field.rfftn()
rfft_field

Field

• Mesh
  • Region
    • pmin = [-0.525, -0.011111111111111112]
    • pmax = [0.47500000000000003, 0.5]
    • dims = ['k_x', 'k_y']
    • units = ['(m)$^{-1}$', '(m)$^{-1}$']
  • n = [20, 23]
• nvdim = 1

Lets also plot the PSD

psd = np.power(np.abs(rfft_field), 2)

psd.mpl()

We can see that the mesh and the plot look different for real Fourier transforms compared to full Fourier transforms. For full Fourier transforms, the mesh is converted from the standard spatial domain to the spatial frequency domain but for real Fourier transforms there are some key differences. As before, each dimension of the mesh (dims) is prepended with k_ to signify this change. However, the Fourier transformed mesh no longer has the the same number of points (n) as the original mesh in all directions. In the final dimension, the number of points is reduced by half (If n is even, the length of the transformed axis is \((n/2)+1\). If n is odd, the length is \((n+1)/2\)).

This halving is a result of the properties of real-valued functions, which have a Hermitian-symmetric Fourier transform. As a result, it is sufficient to store only half of the frequency components (those for non-negative frequencies).

irFFT

As we have an rFFT field we can convert this back to the spatial domain by using an irFFT.

irfft_field = rfft_field.irfftn()
irfft_field

Field

• Mesh
  • Region
    • pmin = [-10.0, -22.500000000000004]
    • pmax = [10.0, 22.500000000000004]
    • dims = ['x', 'y']
    • units = ['m', 'm']
  • n = [20, 44]
• nvdim = 1

As anticipated, the inverse transform function works effectively, yet the total number of points in the final geometric dimension of the mesh is one less than the original. This discrepancy arises due to the characteristics of grids with an odd number of cells.

To reconcile this, we can explicitly specify the desired shape of the output while calling the irfftn function as this information is not stored in the Fourier transformed field. This way, we ensure that the final transformed field aligns with the original field in terms of its dimensionality.

irfft_field = rfft_field.irfftn(shape=field.mesh.n)
irfft_field.mesh.translate(field.mesh.region.centre, inplace=True)
irfft_field

Field

• Mesh
  • Region
    • pmin = [0.0, -3.552713678800501e-15]
    • pmax = [20.0, 45.0]
    • dims = ['x', 'y']
    • units = ['m', 'm']
  • n = [20, 45]
• nvdim = 1

Now we can use the allclose function to see that this is a good comparison to our original field.

field.allclose(irfft_field)

True

Angular frequency

As we noted earlier, the values of the axis in Fourier space are frequency NOT angular frequency so DO NOT include a factor of \(2\pi\), i.e. \(|{\bf k}| = \frac{1}{\lambda} \neq \frac{2 \pi}{\lambda}\), where \(\bf k\) is the wave vector and \(\lambda\) is the wavelength. However, we can convert the mesh into angular frequency.

First we will perform a new Fourier transform and plot the PSD.

fft_field_freq = field.fftn()
fft_field_freq

Field

• Mesh
  • Region
    • pmin = [-0.525, -0.5]
    • pmax = [0.47500000000000003, 0.5]
    • dims = ['k_x', 'k_y']
    • units = ['(m)$^{-1}$', '(m)$^{-1}$']
  • n = [20, 45]
• nvdim = 1

psd = np.power(np.abs(fft_field_freq), 2)

psd.mpl()

The easiest method to convert to angular frequency is by modifying the existing field directly inplace. To achieve this, we can scale the mesh by a factor of \(2\pi\).

fft_field_freq.mesh.scale(2 * np.pi, inplace=True)

Mesh

• Region
  • pmin = [-3.166592653589793, -3.141592653589793]
  • pmax = [3.116592653589793, 3.141592653589793]
  • dims = ['k_x', 'k_y']
  • units = ['(m)$^{-1}$', '(m)$^{-1}$']
• n = [20, 45]

Plotting the PSD of the field again, we can see that the axes labels have expanded by a factor of \(2\pi\).

psd = np.power(np.abs(fft_field_freq), 2)

psd.mpl()

Visualising and analysing FFTs

To demonstrate the power of FFTs and provide some visualisation techniques for their outputs, we have prepared the following example.

We have initialised the field with a complicated three dimensional vector field on a mesh with three spatial dimensions. This example is intricate in order to showcase some of the functionality so don’t worry if you don’t understand all of the details.

We have created a vector field made up of three distinct waves (and additional noise): 1. A helix propagating in the \(y\) direction with a wavelength of \(7\, \mathrm{nm}\). The vector components point in the \(xz\) plane. 2. A density wave propagating in the \(xyz\) direction with a wavelength of \(5\, \mathrm{nm}\). The vector components point in the \(z\) direction. 3. A density wave propagating in the \(xy\) direction with a wavelength of \(10\, \mathrm{nm}\). The vector components point in the \(z\) direction.

mesh = df.Mesh(p1=(0, 0, 0), p2=(30e-9, 30e-9, 30e-9), cell=(0.5e-9, 0.5e-9, 0.5e-9))


def value(p):
    x, y, z = p
    # First wave
    k1 = 1 / 7e-9
    v = np.array([np.sin(2 * np.pi * k1 * y), 0, np.cos(2 * np.pi * k1 * y)])
    # Second wave
    k2 = 1 / 5e-9
    v += 2 * np.array([0, 0, np.cos(2 * np.pi * k2 * (x + y + z) / np.sqrt(3))])
    # Third wave
    k3 = 1 / 10e-9
    v += np.array([0, 0, np.cos(2 * np.pi * k3 * (x + y) / np.sqrt(2))])
    # Random noise
    v += np.random.normal(0, 0.1, 3)
    return v


field = df.Field(mesh, nvdim=3, value=value)

As expected due to the complexity, visualising this vector field can be challenging. While mpl can be used to depict the magnetic structure, it struggles to capture the entire picture due to variations in the vector field across all dimensions.

field.sel("z").mpl()

In this circumstance we can use hv to examine the magnetic texture in more detail.

field.hv.scalar(kdims=["x", "y"])

Now we will perform an FFT on the field and calculate the PSD

fft_field = field.fftn()
fft_field

Field

• Mesh
  • Region
    • pmin = [-1016666666.6666667, -1016666666.6666667, -1016666666.6666667]
    • pmax = [983333333.3333334, 983333333.3333334, 983333333.3333334]
    • dims = ['k_x', 'k_y', 'k_z']
    • units = ['(m)$^{-1}$', '(m)$^{-1}$', '(m)$^{-1}$']
  • n = [60, 60, 60]
• nvdim = 3
• vdims:
  • ft_x
  • ft_y
  • ft_z

psd = np.power(np.abs(fft_field), 2)
psd

Field

• Mesh
  • Region
    • pmin = [-1016666666.6666667, -1016666666.6666667, -1016666666.6666667]
    • pmax = [983333333.3333334, 983333333.3333334, 983333333.3333334]
    • dims = ['k_x', 'k_y', 'k_z']
    • units = ['(m)$^{-1}$', '(m)$^{-1}$', '(m)$^{-1}$']
  • n = [60, 60, 60]
• nvdim = 3
• vdims:
  • ft_x
  • ft_y
  • ft_z

We can use mpl to plot two dimensional cross-sections of the PSD in a similar way to the previous example. Given that two of our propagation vectors reside in the \(xy\) plane, we will select a cut at \(k_z=0\) to visualise the \(k_x k_y\) plane. For analysis, it is often more straightforward to plot and examine one component at a time.

The plot of the ft_x component only shows peaks due to the helix (first wave) as it is the only wave with a vector component in the \(x\) direction.

psd.ft_x.sel("k_z").mpl()

There are not waves with vector components in the \(y\) direction, hence all we see is noise.

psd.ft_y.sel("k_z").mpl()

All three waves have vector components in the \(z\) direction. However, only the helix (first wave) and one of the density waves (third wave) have their propagation vectors lying in the \(xy\) plane (i.e. \(k_z=0\)). Hence we see four peaks corresponding to these two waves.

psd.ft_z.sel("k_z").mpl()

But what about the other density wave (second wave) propagating in the \(xyz\) direction?

We cannot access a \(xyz\) plane in discretisedfield without rotations as the mesh is based on a discretised cartesian grid. However, we can plot the \(k_x k_y\) plane at a non-zero \(k_z\) value to see one of the peaks.

psd.ft_z.sel(k_z=1.15e8).mpl()

We can now clearly see the peak related to the density wave propagating in the \(xyz\) \((111)\) direction.

An easier way to visualise this data is again using the hv plotting functionality. This allows us to select the component of the field we want to look at while using a slider to move through the other dimensions.

Take a while to play with the plot. Make sure to look at all of the components at the center plane \((k_z=0\,\mathrm{m}^{-1})\) and also where we expect to see out other propagation vector \(\left(k_z=\pm \frac{1}{5\times 10 ^{-9}\sqrt{3}}\,\mathrm{m}^{-1} =\pm 1.15\times 10 ^{8}\,\mathrm{m}^{-1}\right)\). Further information of plotting using hv can be found in the Documentation and API reference.

psd.hv.scalar(kdims=["k_x", "k_y"], clim=(0, psd.array.max()))

Rather than plotting a plane we can use the line functionality in discretisedfield to examine the values along a one dimensional line in any arbitrary direction. Here we will create a diagonal line from the bottom left corner of our mesh \(\left(p_\mathrm{min}\right)\) to the top right corner of our mesh \(\left(p_\mathrm{max}\right)\). For further information on the line method please see the documentation and API reference.

line = psd.line(p1=psd.mesh.region.pmin, p2=psd.mesh.region.pmax, n=60)
line.data

	r	k_x	k_y	k_z	vft_x	vft_y	vft_z
0	0.000000e+00	-1.016667e+09	-1.016667e+09	-1.016667e+09	1.533123e+02	967.133826	4.288057e+03
1	5.871359e+07	-9.827684e+08	-9.827684e+08	-9.827684e+08	1.847779e+03	1158.211500	7.352008e+03
2	1.174272e+08	-9.488701e+08	-9.488701e+08	-9.488701e+08	7.947446e+02	4409.870244	3.265669e+02
3	1.761408e+08	-9.149718e+08	-9.149718e+08	-9.149718e+08	1.304691e+03	562.386254	2.033242e+02
4	2.348543e+08	-8.810734e+08	-8.810734e+08	-8.810734e+08	1.980885e+02	682.398721	9.901236e+02
5	2.935679e+08	-8.471751e+08	-8.471751e+08	-8.471751e+08	3.188206e+02	2177.488013	2.383882e+02
6	3.522815e+08	-8.132768e+08	-8.132768e+08	-8.132768e+08	9.382068e+03	1556.538399	5.629643e+02
7	4.109951e+08	-7.793785e+08	-7.793785e+08	-7.793785e+08	1.283730e+03	1158.909728	7.189914e+02
8	4.697087e+08	-7.454802e+08	-7.454802e+08	-7.454802e+08	3.769435e+01	1683.641055	2.109827e+03
9	5.284223e+08	-7.115819e+08	-7.115819e+08	-7.115819e+08	3.637889e+03	3599.612375	1.371286e+03
10	5.871359e+08	-6.776836e+08	-6.776836e+08	-6.776836e+08	4.212858e+03	81.808281	3.167806e+03
11	6.458495e+08	-6.437853e+08	-6.437853e+08	-6.437853e+08	7.499259e+02	1453.503092	1.171385e+03
12	7.045630e+08	-6.098870e+08	-6.098870e+08	-6.098870e+08	3.549571e+02	3877.647581	6.529764e+03
13	7.632766e+08	-5.759887e+08	-5.759887e+08	-5.759887e+08	2.796901e+02	6730.064855	1.402915e+03
14	8.219902e+08	-5.420904e+08	-5.420904e+08	-5.420904e+08	1.278498e+03	3391.607301	3.417761e+03
15	8.807038e+08	-5.081921e+08	-5.081921e+08	-5.081921e+08	3.345894e+03	529.881944	4.418821e+03
16	9.394174e+08	-4.742938e+08	-4.742938e+08	-4.742938e+08	2.493517e+03	428.868239	3.075442e+03
17	9.981310e+08	-4.403955e+08	-4.403955e+08	-4.403955e+08	1.094941e+03	838.860841	4.024590e+03
18	1.056845e+09	-4.064972e+08	-4.064972e+08	-4.064972e+08	8.956192e+02	1430.065794	6.318671e+03
19	1.115558e+09	-3.725989e+08	-3.725989e+08	-3.725989e+08	1.263918e+03	1234.678006	8.904181e+03
20	1.174272e+09	-3.387006e+08	-3.387006e+08	-3.387006e+08	1.192933e+03	2959.557627	5.674751e+03
21	1.232985e+09	-3.048023e+08	-3.048023e+08	-3.048023e+08	6.265783e+03	4182.300782	3.237929e+03
22	1.291699e+09	-2.709040e+08	-2.709040e+08	-2.709040e+08	1.921405e+03	2750.472104	8.690839e+03
23	1.350412e+09	-2.370056e+08	-2.370056e+08	-2.370056e+08	1.059326e+03	96.690581	2.470224e+04
24	1.409126e+09	-2.031073e+08	-2.031073e+08	-2.031073e+08	1.051780e+03	2841.657226	2.013952e+05
25	1.467840e+09	-1.692090e+08	-1.692090e+08	-1.692090e+08	1.454145e+03	2922.751065	3.601955e+06
26	1.526553e+09	-1.353107e+08	-1.353107e+08	-1.353107e+08	1.188864e+03	2110.892584	1.974571e+09
27	1.585267e+09	-1.014124e+08	-1.014124e+08	-1.014124e+08	6.252956e+03	2065.062863	4.671110e+09
28	1.643980e+09	-6.751412e+07	-6.751412e+07	-6.751412e+07	1.059485e+03	7540.272416	4.641762e+06
29	1.702694e+09	-3.361582e+07	-3.361582e+07	-3.361582e+07	4.920847e+03	363.959313	1.525691e+05
30	1.761408e+09	2.824859e+05	2.824859e+05	2.824859e+05	9.768406e+07	0.709852	6.876681e+07
31	1.820121e+09	3.418079e+07	3.418079e+07	3.418079e+07	4.920847e+03	363.959313	1.525691e+05
32	1.878835e+09	6.807910e+07	6.807910e+07	6.807910e+07	1.059485e+03	7540.272416	4.641762e+06
33	1.937548e+09	1.019774e+08	1.019774e+08	1.019774e+08	6.252956e+03	2065.062863	4.671110e+09
34	1.996262e+09	1.358757e+08	1.358757e+08	1.358757e+08	1.188864e+03	2110.892584	1.974571e+09
35	2.054976e+09	1.697740e+08	1.697740e+08	1.697740e+08	1.454145e+03	2922.751065	3.601955e+06
36	2.113689e+09	2.036723e+08	2.036723e+08	2.036723e+08	1.051780e+03	2841.657226	2.013952e+05
37	2.172403e+09	2.375706e+08	2.375706e+08	2.375706e+08	1.059326e+03	96.690581	2.470224e+04
38	2.231116e+09	2.714689e+08	2.714689e+08	2.714689e+08	1.921405e+03	2750.472104	8.690839e+03
39	2.289830e+09	3.053672e+08	3.053672e+08	3.053672e+08	6.265783e+03	4182.300782	3.237929e+03
40	2.348543e+09	3.392655e+08	3.392655e+08	3.392655e+08	1.192933e+03	2959.557627	5.674751e+03
41	2.407257e+09	3.731638e+08	3.731638e+08	3.731638e+08	1.263918e+03	1234.678006	8.904181e+03
42	2.465971e+09	4.070621e+08	4.070621e+08	4.070621e+08	8.956192e+02	1430.065794	6.318671e+03
43	2.524684e+09	4.409605e+08	4.409605e+08	4.409605e+08	1.094941e+03	838.860841	4.024590e+03
44	2.583398e+09	4.748588e+08	4.748588e+08	4.748588e+08	2.493517e+03	428.868239	3.075442e+03
45	2.642111e+09	5.087571e+08	5.087571e+08	5.087571e+08	3.345894e+03	529.881944	4.418821e+03
46	2.700825e+09	5.426554e+08	5.426554e+08	5.426554e+08	1.278498e+03	3391.607301	3.417761e+03
47	2.759539e+09	5.765537e+08	5.765537e+08	5.765537e+08	2.796901e+02	6730.064855	1.402915e+03
48	2.818252e+09	6.104520e+08	6.104520e+08	6.104520e+08	3.549571e+02	3877.647581	6.529764e+03
49	2.876966e+09	6.443503e+08	6.443503e+08	6.443503e+08	7.499259e+02	1453.503092	1.171385e+03
50	2.935679e+09	6.782486e+08	6.782486e+08	6.782486e+08	4.212858e+03	81.808281	3.167806e+03
51	2.994393e+09	7.121469e+08	7.121469e+08	7.121469e+08	3.637889e+03	3599.612375	1.371286e+03
52	3.053107e+09	7.460452e+08	7.460452e+08	7.460452e+08	3.769435e+01	1683.641055	2.109827e+03
53	3.111820e+09	7.799435e+08	7.799435e+08	7.799435e+08	1.283730e+03	1158.909728	7.189914e+02
54	3.170534e+09	8.138418e+08	8.138418e+08	8.138418e+08	9.382068e+03	1556.538399	5.629643e+02
55	3.229247e+09	8.477401e+08	8.477401e+08	8.477401e+08	3.188206e+02	2177.488013	2.383882e+02
56	3.287961e+09	8.816384e+08	8.816384e+08	8.816384e+08	1.980885e+02	682.398721	9.901236e+02
57	3.346674e+09	9.155367e+08	9.155367e+08	9.155367e+08	1.304691e+03	562.386254	2.033242e+02
58	3.405388e+09	9.494350e+08	9.494350e+08	9.494350e+08	7.947446e+02	4409.870244	3.265669e+02
59	3.464102e+09	9.833333e+08	9.833333e+08	9.833333e+08	1.847779e+03	1158.211500	7.352008e+03

We can process the data to help it to be more meaningful. We will change the length along the line r to have its zero point at the origin.

line.data["r"] = np.sqrt(3) * line.data["k_x"]
line.data

	r	k_x	k_y	k_z	vft_x	vft_y	vft_z
0	-1.760918e+09	-1.016667e+09	-1.016667e+09	-1.016667e+09	1.533123e+02	967.133826	4.288057e+03
1	-1.702205e+09	-9.827684e+08	-9.827684e+08	-9.827684e+08	1.847779e+03	1158.211500	7.352008e+03
2	-1.643491e+09	-9.488701e+08	-9.488701e+08	-9.488701e+08	7.947446e+02	4409.870244	3.265669e+02
3	-1.584778e+09	-9.149718e+08	-9.149718e+08	-9.149718e+08	1.304691e+03	562.386254	2.033242e+02
4	-1.526064e+09	-8.810734e+08	-8.810734e+08	-8.810734e+08	1.980885e+02	682.398721	9.901236e+02
5	-1.467350e+09	-8.471751e+08	-8.471751e+08	-8.471751e+08	3.188206e+02	2177.488013	2.383882e+02
6	-1.408637e+09	-8.132768e+08	-8.132768e+08	-8.132768e+08	9.382068e+03	1556.538399	5.629643e+02
7	-1.349923e+09	-7.793785e+08	-7.793785e+08	-7.793785e+08	1.283730e+03	1158.909728	7.189914e+02
8	-1.291210e+09	-7.454802e+08	-7.454802e+08	-7.454802e+08	3.769435e+01	1683.641055	2.109827e+03
9	-1.232496e+09	-7.115819e+08	-7.115819e+08	-7.115819e+08	3.637889e+03	3599.612375	1.371286e+03
10	-1.173782e+09	-6.776836e+08	-6.776836e+08	-6.776836e+08	4.212858e+03	81.808281	3.167806e+03
11	-1.115069e+09	-6.437853e+08	-6.437853e+08	-6.437853e+08	7.499259e+02	1453.503092	1.171385e+03
12	-1.056355e+09	-6.098870e+08	-6.098870e+08	-6.098870e+08	3.549571e+02	3877.647581	6.529764e+03
13	-9.976417e+08	-5.759887e+08	-5.759887e+08	-5.759887e+08	2.796901e+02	6730.064855	1.402915e+03
14	-9.389281e+08	-5.420904e+08	-5.420904e+08	-5.420904e+08	1.278498e+03	3391.607301	3.417761e+03
15	-8.802145e+08	-5.081921e+08	-5.081921e+08	-5.081921e+08	3.345894e+03	529.881944	4.418821e+03
16	-8.215009e+08	-4.742938e+08	-4.742938e+08	-4.742938e+08	2.493517e+03	428.868239	3.075442e+03
17	-7.627873e+08	-4.403955e+08	-4.403955e+08	-4.403955e+08	1.094941e+03	838.860841	4.024590e+03
18	-7.040738e+08	-4.064972e+08	-4.064972e+08	-4.064972e+08	8.956192e+02	1430.065794	6.318671e+03
19	-6.453602e+08	-3.725989e+08	-3.725989e+08	-3.725989e+08	1.263918e+03	1234.678006	8.904181e+03
20	-5.866466e+08	-3.387006e+08	-3.387006e+08	-3.387006e+08	1.192933e+03	2959.557627	5.674751e+03
21	-5.279330e+08	-3.048023e+08	-3.048023e+08	-3.048023e+08	6.265783e+03	4182.300782	3.237929e+03
22	-4.692194e+08	-2.709040e+08	-2.709040e+08	-2.709040e+08	1.921405e+03	2750.472104	8.690839e+03
23	-4.105058e+08	-2.370056e+08	-2.370056e+08	-2.370056e+08	1.059326e+03	96.690581	2.470224e+04
24	-3.517922e+08	-2.031073e+08	-2.031073e+08	-2.031073e+08	1.051780e+03	2841.657226	2.013952e+05
25	-2.930787e+08	-1.692090e+08	-1.692090e+08	-1.692090e+08	1.454145e+03	2922.751065	3.601955e+06
26	-2.343651e+08	-1.353107e+08	-1.353107e+08	-1.353107e+08	1.188864e+03	2110.892584	1.974571e+09
27	-1.756515e+08	-1.014124e+08	-1.014124e+08	-1.014124e+08	6.252956e+03	2065.062863	4.671110e+09
28	-1.169379e+08	-6.751412e+07	-6.751412e+07	-6.751412e+07	1.059485e+03	7540.272416	4.641762e+06
29	-5.822431e+07	-3.361582e+07	-3.361582e+07	-3.361582e+07	4.920847e+03	363.959313	1.525691e+05
30	4.892799e+05	2.824859e+05	2.824859e+05	2.824859e+05	9.768406e+07	0.709852	6.876681e+07
31	5.920287e+07	3.418079e+07	3.418079e+07	3.418079e+07	4.920847e+03	363.959313	1.525691e+05
32	1.179165e+08	6.807910e+07	6.807910e+07	6.807910e+07	1.059485e+03	7540.272416	4.641762e+06
33	1.766300e+08	1.019774e+08	1.019774e+08	1.019774e+08	6.252956e+03	2065.062863	4.671110e+09
34	2.353436e+08	1.358757e+08	1.358757e+08	1.358757e+08	1.188864e+03	2110.892584	1.974571e+09
35	2.940572e+08	1.697740e+08	1.697740e+08	1.697740e+08	1.454145e+03	2922.751065	3.601955e+06
36	3.527708e+08	2.036723e+08	2.036723e+08	2.036723e+08	1.051780e+03	2841.657226	2.013952e+05
37	4.114844e+08	2.375706e+08	2.375706e+08	2.375706e+08	1.059326e+03	96.690581	2.470224e+04
38	4.701980e+08	2.714689e+08	2.714689e+08	2.714689e+08	1.921405e+03	2750.472104	8.690839e+03
39	5.289116e+08	3.053672e+08	3.053672e+08	3.053672e+08	6.265783e+03	4182.300782	3.237929e+03
40	5.876251e+08	3.392655e+08	3.392655e+08	3.392655e+08	1.192933e+03	2959.557627	5.674751e+03
41	6.463387e+08	3.731638e+08	3.731638e+08	3.731638e+08	1.263918e+03	1234.678006	8.904181e+03
42	7.050523e+08	4.070621e+08	4.070621e+08	4.070621e+08	8.956192e+02	1430.065794	6.318671e+03
43	7.637659e+08	4.409605e+08	4.409605e+08	4.409605e+08	1.094941e+03	838.860841	4.024590e+03
44	8.224795e+08	4.748588e+08	4.748588e+08	4.748588e+08	2.493517e+03	428.868239	3.075442e+03
45	8.811931e+08	5.087571e+08	5.087571e+08	5.087571e+08	3.345894e+03	529.881944	4.418821e+03
46	9.399067e+08	5.426554e+08	5.426554e+08	5.426554e+08	1.278498e+03	3391.607301	3.417761e+03
47	9.986203e+08	5.765537e+08	5.765537e+08	5.765537e+08	2.796901e+02	6730.064855	1.402915e+03
48	1.057334e+09	6.104520e+08	6.104520e+08	6.104520e+08	3.549571e+02	3877.647581	6.529764e+03
49	1.116047e+09	6.443503e+08	6.443503e+08	6.443503e+08	7.499259e+02	1453.503092	1.171385e+03
50	1.174761e+09	6.782486e+08	6.782486e+08	6.782486e+08	4.212858e+03	81.808281	3.167806e+03
51	1.233475e+09	7.121469e+08	7.121469e+08	7.121469e+08	3.637889e+03	3599.612375	1.371286e+03
52	1.292188e+09	7.460452e+08	7.460452e+08	7.460452e+08	3.769435e+01	1683.641055	2.109827e+03
53	1.350902e+09	7.799435e+08	7.799435e+08	7.799435e+08	1.283730e+03	1158.909728	7.189914e+02
54	1.409615e+09	8.138418e+08	8.138418e+08	8.138418e+08	9.382068e+03	1556.538399	5.629643e+02
55	1.468329e+09	8.477401e+08	8.477401e+08	8.477401e+08	3.188206e+02	2177.488013	2.383882e+02
56	1.527043e+09	8.816384e+08	8.816384e+08	8.816384e+08	1.980885e+02	682.398721	9.901236e+02
57	1.585756e+09	9.155367e+08	9.155367e+08	9.155367e+08	1.304691e+03	562.386254	2.033242e+02
58	1.644470e+09	9.494350e+08	9.494350e+08	9.494350e+08	7.947446e+02	4409.870244	3.265669e+02
59	1.703183e+09	9.833333e+08	9.833333e+08	9.833333e+08	1.847779e+03	1158.211500	7.352008e+03

To plot this line we can use mpl along with adding some axis to the plot.

import matplotlib.pyplot as plt

fig, ax = plt.subplots()
line.mpl(ax=ax)

# Provide titles and labels
ax.set_xlabel("Distance from origin (nm)$^{-1}$")
ax.set_ylabel("PSD (arb.)")

plt.show()

We can clearly see the two peaks originating from the the density wave propagating in the \(xyz\) direction. As expected from the initial wavelength that was set, the peaks occur at \(k = \pm \frac{1}{5}\,(\mathrm{nm})^{-1}= \pm 0.2\,(\mathrm{nm})^{-1}\).