diff --git a/Code/Source/solver/CMakeLists.txt b/Code/Source/solver/CMakeLists.txt
index c546c2822..30294bf16 100644
--- a/Code/Source/solver/CMakeLists.txt
+++ b/Code/Source/solver/CMakeLists.txt
@@ -156,7 +156,8 @@ set(CSRCS
   PetscLinearAlgebra.h PetscLinearAlgebra.cpp
   TrilinosLinearAlgebra.h TrilinosLinearAlgebra.cpp
   Tensor4.h Tensor4.cpp
-  Vector.h Vector.cpp 
+  Vector.h Vector.cpp
+  fourier_interpolation.cpp
 
   lapack_defs.h
 
diff --git a/Code/Source/solver/CmMod.cpp b/Code/Source/solver/CmMod.cpp
index f3485a2bb..2fbfb1c86 100644
--- a/Code/Source/solver/CmMod.cpp
+++ b/Code/Source/solver/CmMod.cpp
@@ -138,6 +138,11 @@ void cmType::bcast(const CmMod& cm_mod, Vector<int>& data) const
   MPI_Bcast(data.data(), data.size(), cm_mod::mpint, cm_mod.master, com());
 }
 
+/// @brief bcast unsigned int
+void cmType::bcast(const CmMod &cm_mod, unsigned int *data) const {
+  MPI_Bcast(data, 1, cm_mod::mpuint, cm_mod.master, com());
+}
+
 /// @brief gather int array
 void cmType::gather(const CmMod& cm_mod, const int* send_data, int send_count, int* recv_data, int recv_count, int root) const
 {
diff --git a/Code/Source/solver/CmMod.h b/Code/Source/solver/CmMod.h
index 36ac60e18..5d2587b67 100644
--- a/Code/Source/solver/CmMod.h
+++ b/Code/Source/solver/CmMod.h
@@ -24,6 +24,7 @@ using MpiCommWorldType = MPI_Comm;
 const decltype(MPI_CXX_BOOL) mplog = MPI_CXX_BOOL;
 //const decltype(MPI_LOGICAL) mplog  = MPI_LOGICAL;
 const decltype(MPI_INTEGER) mpint = MPI_INTEGER;
+const decltype(MPI_UNSIGNED) mpuint = MPI_UNSIGNED;
 const decltype(MPI_DOUBLE_PRECISION) mpreal = MPI_DOUBLE_PRECISION;
 const decltype(MPI_CHARACTER) mpchar = MPI_CHARACTER;
 };
@@ -83,6 +84,8 @@ class cmType {
     void bcast(const CmMod& cm_mod, int* data) const;
     void bcast(const CmMod& cm_mod, Vector<int>& data) const;
 
+    void bcast(const CmMod &cm_mod, unsigned int *data) const;
+
     void bcast(const CmMod& cm_mod, Array<int>& data, const std::string& name="") const;
 
     // Gather operations
diff --git a/Code/Source/solver/ComMod.h b/Code/Source/solver/ComMod.h
index 86db290f4..1d00fd230 100644
--- a/Code/Source/solver/ComMod.h
+++ b/Code/Source/solver/ComMod.h
@@ -8,19 +8,20 @@
 // are defined here.
 
 #ifndef COMMOD_H 
-#define COMMOD_H 
+#define COMMOD_H
 
 #include "Array.h"
 #include "Array3.h"
-#include "SolutionStates.h"
 #include "CepMod.h"
 #include "ChnlMod.h"
 #include "CmMod.h"
+#include "CoupledBoundaryCondition.h"
 #include "Parameters.h"
 #include "RobinBoundaryCondition.h"
-#include "CoupledBoundaryCondition.h"
+#include "SolutionStates.h"
 #include "Timer.h"
 #include "Vector.h"
+#include "fourier_interpolation.h"
 
 #include "DebugMsg.h"
 
@@ -43,42 +44,6 @@
 
 class LinearAlgebra;
 
-/// @brief Fourier coefficients that are used to specify unsteady BCs
-//
-class fcType
-{
-  public:
-
-    bool defined() { return n != 0; };
-
-    // If this is a ramp function
-    bool lrmp = false;
-
-    // Number of Fourier coefficient
-    int n = 0;
-
-    // No. of dimensions (scalar or vector)
-    int d = 0;
-   
-    // Initial value
-    Vector<double> qi;
-
-    // Time derivative of linear part
-    Vector<double> qs;
-
-    // Period
-    double T = 0.0;
-
-    // Initial time
-    double ti = 0.0;
-
-    // Imaginary part of coefficint
-    Array<double> i;
-
-    // Real part of coefficint
-    Array<double> r;
-};
-
 /// @brief Moving boundary data structure (used for general BC)
 //
 class MBType
@@ -194,7 +159,7 @@ class bcType
     //
     // This is declare ALLOCATABLE in MOD.f. 
     //
-    fcType gt;
+    FourierInterpolation gt;
 
     // Neu: RCR
     rcrType RCR;
@@ -306,7 +271,7 @@ class bfType
     Array<double> bx;
 
     // Time dependant (unsteady imposed value)
-    fcType bt;
+    FourierInterpolation bt;
 
     // General (unsteady and spatially dependent combination)
     MBType bm;
@@ -334,7 +299,7 @@ class fibStrsType
     double eta_n = 0.0;
 
     // Unsteady time-dependent values
-    fcType gt;
+    FourierInterpolation gt;
 };
 
 /// @brief Structural domain type
diff --git a/Code/Source/solver/Core/Exception.h b/Code/Source/solver/Core/Exception.h
index 80e6968ec..c60e92d9d 100644
--- a/Code/Source/solver/Core/Exception.h
+++ b/Code/Source/solver/Core/Exception.h
@@ -346,6 +346,23 @@ class DependencyException : public CoreException {
     }
 };
 
+class FileNotFoundException : public CoreException {
+public:
+  FileNotFoundException(const std::string &file_name, const char *file = "",
+                        int line = 0, const char *function = "")
+      : CoreException("Could not open file " + file_name, StatusCode::IOError,
+                      file, line, function) {}
+};
+
+class FileFormatException : public CoreException {
+public:
+  FileFormatException(const std::string &file_name, const std::string &message,
+                      const char *file = "", int line = 0,
+                      const char *function = "")
+      : CoreException("Error parsing file " + file_name + ". " + message,
+                      StatusCode::IOError, file, line, function) {}
+};
+
 inline void ExceptionRuntime::install_terminate_handler()
 {
     std::set_terminate([]() {
diff --git a/Code/Source/solver/bf.cpp b/Code/Source/solver/bf.cpp
index 1574b5df0..2832db0ee 100644
--- a/Code/Source/solver/bf.cpp
+++ b/Code/Source/solver/bf.cpp
@@ -151,8 +151,7 @@ void set_bf_l(ComMod& com_mod, bfType& lBf, mshType& lM, const SolutionStates& s
     f = lBf.b;
 
   } else if (utils::btest(lBf.bType, enum_int(BodyForceType::bfType_ustd))) {
-    Vector<double> rtmp(1);
-    ifft(com_mod, lBf.bt, f, rtmp);
+    f = lBf.bt.value(com_mod.time);
 
   } else if (utils::btest(lBf.bType, enum_int(BodyForceType::bfType_gen))) {
     bfl.resize(idof,nNo); 
diff --git a/Code/Source/solver/distribute.cpp b/Code/Source/solver/distribute.cpp
index 64acce552..742f01d31 100644
--- a/Code/Source/solver/distribute.cpp
+++ b/Code/Source/solver/distribute.cpp
@@ -644,44 +644,10 @@ void dist_bc(ComMod& com_mod, const CmMod& cm_mod, const cmType& cm, bcType& lBc
   }
 
   // Communicating time-dependent BC data
-  //
-  // lBc.gt is declare ALLOCATABLE in MOD.f but we don't
-  // want to use pointers so use the define() method
-  // to check if it has data define.
-  //
-  bool flag = lBc.gt.defined(); 
-  cm.bcast(cm_mod, &flag);
-
-  if (flag) {
-    if (is_slave) {
-      // [NOTE] This is allocated in ComMod.
-      //lBc.gt = new fcType;
-    }
-
-    cm.bcast(cm_mod, &lBc.gt.lrmp);
-    cm.bcast(cm_mod, &lBc.gt.d);
-    cm.bcast(cm_mod, &lBc.gt.n);
-
-    int j = lBc.gt.d;
-    int i = lBc.gt.n;
-
-    if (is_slave) { 
-      lBc.gt.qi.resize(j);
-      lBc.gt.qs.resize(j);
-      lBc.gt.r.resize(j,i);
-      lBc.gt.i.resize(j,i);
-    }
-
-    cm.bcast(cm_mod, &lBc.gt.ti);
-    cm.bcast(cm_mod, &lBc.gt.T);
-    cm.bcast(cm_mod, lBc.gt.qi);
-    cm.bcast(cm_mod, lBc.gt.qs);
-    cm.bcast(cm_mod, lBc.gt.r);
-    cm.bcast(cm_mod, lBc.gt.i);
-  }
+  lBc.gt.distribute(cm_mod, cm);
 
   // Communicating moving BC data
-  flag = lBc.gm.defined();
+  bool flag = lBc.gm.defined();
   cm.bcast(cm_mod, &flag);
   
   if (flag) {
@@ -917,31 +883,7 @@ void dist_bf(ComMod& com_mod, const CmMod& cm_mod, const cmType& cm, bfType& lBf
   }
 
   //  Communicating time-dependent BF data
-  flag = lBf.bt.defined();
-  cm.bcast(cm_mod, &flag);
-  if (flag) {
-    if (is_slave) {
-      //ALLOCATE(lBf.bt)
-    }
-    cm.bcast(cm_mod, &lBf.bt.lrmp);
-    cm.bcast(cm_mod, &lBf.bt.d);
-    cm.bcast(cm_mod, &lBf.bt.n);
-
-    if (is_slave) {
-      int j = lBf.bt.d;
-      int i = lBf.bt.n;
-      lBf.bt.qi.resize(j);
-      lBf.bt.qs.resize(j);
-      lBf.bt.r.resize(j,i);
-      lBf.bt.i.resize(j,i);
-    }
-    cm.bcast(cm_mod, &lBf.bt.ti);
-    cm.bcast(cm_mod, &lBf.bt.T);
-    cm.bcast(cm_mod, lBf.bt.qi);
-    cm.bcast(cm_mod, lBf.bt.qs);
-    cm.bcast(cm_mod, lBf.bt.r);
-    cm.bcast(cm_mod, lBf.bt.i);
-  }
+  lBf.bt.distribute(cm_mod, cm);
 
   // Communicating moving BF data
   //
@@ -1731,26 +1673,7 @@ void dist_mat_consts(const ComMod& com_mod, const CmMod& cm_mod, const cmType& c
     cm.bcast(cm_mod, &lStM.Tf.g);
 
   } else if (utils::btest(lStM.Tf.fType, static_cast<int>(BoundaryConditionType::bType_ustd))) {
-    cm.bcast(cm_mod, &lStM.Tf.gt.lrmp);
-    cm.bcast(cm_mod, &lStM.Tf.gt.d);
-    cm.bcast(cm_mod, &lStM.Tf.gt.n);
-
-    if (cm.slv(cm_mod)) {
-      int j = lStM.Tf.gt.d;
-      int i = lStM.Tf.gt.n;
-      lStM.Tf.gt.qi.resize(j);
-      lStM.Tf.gt.qs.resize(j);
-      lStM.Tf.gt.r.resize(j,i);
-      lStM.Tf.gt.i.resize(j,i);
-   } 
-
-   cm.bcast(cm_mod, &lStM.Tf.gt.ti);
-   cm.bcast(cm_mod, &lStM.Tf.gt.T);
-
-   cm.bcast(cm_mod, lStM.Tf.gt.qi, "lStM.Tf.gt.qi");
-   cm.bcast(cm_mod, lStM.Tf.gt.qs, "lStM.Tf.gt.qs");
-   cm.bcast(cm_mod, lStM.Tf.gt.r, "lStM.Tf.gt.r");
-   cm.bcast(cm_mod, lStM.Tf.gt.i, "lStM.Tf.gt.i");
+    lStM.Tf.gt.distribute(cm_mod, cm);
   }
   
   // Broadcast directional stress distribution parameters
diff --git a/Code/Source/solver/fft.cpp b/Code/Source/solver/fft.cpp
index 5a7c79bf9..5a362ebd0 100644
--- a/Code/Source/solver/fft.cpp
+++ b/Code/Source/solver/fft.cpp
@@ -6,159 +6,6 @@
 #include "fft.h"
 #include <math.h>
 
-/// @brief Replicates Fortran 'SUBROUTINE FFT(fid, np, gt)'.
-///
-/// temporal_values contains all of the values read in from
-/// the temporal values file.
-//
-void fft(const int np, const std::vector<std::vector<double>>& temporal_values, fcType& gt)
-{
-  using namespace consts;
-
-  #define n_debug_fft
-  #ifdef debug_fft
-  DebugMsg dmsg(__func__, 0);
-  dmsg.banner();
-  dmsg << "np: " << np;
-  #endif
-
-  Vector<double> t(np);
-  Array<double> q(gt.d, np);
-
-  for (int i = 0; i < np; i++) {
-    t[i] = temporal_values[i][0];
-    for (int j = 0; j < gt.d; j++) {
-      q(j,i) = temporal_values[i][j+1];
-    }
-  }
-
-
-  gt.ti = t[0];
-  gt.T = t[np-1] - t[0];
-  #ifdef debug_fft
-  dmsg << "pi: " << pi;
-  dmsg << "gt.ti: " << gt.ti;
-  dmsg << "gt.T: " << gt.T;
-  dmsg << "gt.d: " << gt.d;
-  dmsg << "gt.n: " << gt.n;
-  dmsg << "np: " << np;
-  dmsg << "gt.r.nrows(): " << gt.r.nrows();
-  #endif
-
-  for (int j = 0; j < gt.d; j++) {
-    gt.qi[j] = q(j,0);
-    gt.qs[j] = (q(j,np-1) - q(j,0)) / gt.T;
-  }
-
-  for (int i = 0; i < np; i++) {
-    t[i] = t[i] - gt.ti;
-    for (int j = 0; j < gt.d; j++) {
-      q(j,i) = q(j,i) - gt.qi[j] - gt.qs[j]*t[i];
-    }
-  }
-
-  int ns = 0;
-  Vector<double> ns_array(512);
-
-  for (int n = 0; n < gt.n; n++) {
-    double tmp = static_cast<double>(n);
-    gt.r.set_col(n, 0.0);
-    gt.i.set_col(n, 0.0);
-
-    for (int i = 0; i < np-1; i++) {
-      double ko = 2.0 * pi * tmp * t[i] / gt.T;
-      double kn = 2.0 * pi * tmp * t[i+1] / gt.T;
-
-      for (int j = 0; j < gt.d; j++) {
-        double s = (q(j,i+1) - q(j,i)) / (t[i+1] - t[i]);
-        if (n == 0) {
-          gt.r(j,n) = gt.r(j,n) + 0.5*(t[i+1] - t[i]) * (q(j,i+1) + q(j,i));
-        } else {
-          gt.r(j,n) = gt.r(j,n) + s*(cos(kn) - cos(ko));
-          gt.i(j,n) = gt.i(j,n) - s*(sin(kn) - sin(ko));
-        }
-      }
-    }
-
-    if (n == 0) {
-      for (int k = 0; k < gt.d; k++) {
-        gt.r(k,n) = gt.r(k,n) / gt.T;
-      }
-    } else { 
-      for (int k = 0; k < gt.d; k++) {
-        gt.r(k,n) = 0.5 * gt.r(k,n) * gt.T / (pi*pi*tmp*tmp);
-        gt.i(k,n) = 0.5 * gt.i(k,n) * gt.T / (pi*pi*tmp*tmp);
-      }
-    }
-  }
-}
-
-/// @brief This is to calculate flow rate and flow acceleration (IFFT)
-//
-void ifft(const ComMod& com_mod, const fcType& gt, Vector<double>& Y, Vector<double>& dY) 
-{
-  using namespace consts;
-  #define n_debug_ifft
-  #ifdef debug_ifft
-  DebugMsg dmsg(__func__, com_mod.cm.idcm());
-  dmsg.banner();
-  #endif
-
-  double time = com_mod.time;
-  #ifdef debug_ifft
-  dmsg << "time: " << time;
-  dmsg << "gt.lrmp: " << gt.lrmp;
-  dmsg << "gt.ti: " << gt.ti;
-  dmsg << "gt.n: " << gt.n;
-  dmsg << "pi: " << pi;
-  #endif
-
-  if (gt.lrmp) {
-    double t = time - gt.ti;
-    if (t <= 0.0) {
-      t = fmax(t, 0.0);
-    } else {
-      t = fmin(t, gt.T);
-    }
-    #ifdef debug_ifft
-    dmsg << "t: " <<  t;
-    #endif
-
-    for (int i = 0; i < gt.d; i++) { 
-      #ifdef debug_ifft
-      dmsg << " gt.qi(i): " <<  gt.qi(i);
-      dmsg << " gt.qs(i): " <<  gt.qs(i);
-      #endif
-      Y(i) = gt.qi(i) + t*gt.qs(i);
-      dY(i) = gt.qs(i);
-    }
-
-  } else {
-    double t  = fmod(time - gt.ti, gt.T);
-    double tmp = 2.0 * pi / gt.T;
-    #ifdef debug_ifft
-    dmsg << "t: " << t;
-    dmsg << "tmp: " << tmp;
-    dmsg << "gt.d: " << gt.d;
-    #endif
-
-    for (int j = 0; j < gt.d; j++) { 
-      Y(j) = gt.qi(j) + t*gt.qs(j);
-      dY(j) = gt.qs(j);
-    }
-
-    for (int i = 0; i < gt.n; i++) { 
-      double dk = tmp * static_cast<double>(i);
-      double K = t*dk;
-
-      for (int j = 0; j < gt.d; j++) { 
-        Y(j) = Y(j) +  gt.r(j,i)*cos(K) - gt.i(j,i)*sin(K);
-        dY(j) = dY(j) - (gt.r(j,i)*sin(K) + gt.i(j,i)*cos(K))*dk;
-      }
-    }
-  }
-}
-
 /// @brief This routine is for calculating values by the inverse of general BC
 //
 void igbc(const ComMod& com_mod, const MBType& gm, Array<double>& Y, Array<double>& dY)
@@ -183,7 +30,5 @@ void igbc(const ComMod& com_mod, const MBType& gm, Array<double>& Y, Array<doubl
       Y(j,a) = tmp*gm.d(j,a,i+1) + gm.d(j,a,i)*(1.0-tmp);
       dY(j,a) = (gm.d(j,a,i+1) - gm.d(j,a,i)) / delT;
     }
-  } 
-}
-
-
+  }
+}
\ No newline at end of file
diff --git a/Code/Source/solver/fft.h b/Code/Source/solver/fft.h
index 5ffeb0fa2..9fbc1c45a 100644
--- a/Code/Source/solver/fft.h
+++ b/Code/Source/solver/fft.h
@@ -8,10 +8,6 @@
 
 #include <vector>
 
-void fft(const int np, const std::vector<std::vector<double>>& temporal_values, fcType& gt);
-
-void ifft(const ComMod& com_mod, const fcType& gt, Vector<double>& Y, Vector<double>& dY);
-
 void igbc(const ComMod& com_mod, const MBType& gm, Array<double>& Y, Array<double>& dY);
 
 #endif
diff --git a/Code/Source/solver/fourier_interpolation.cpp b/Code/Source/solver/fourier_interpolation.cpp
new file mode 100644
index 000000000..12aaefaff
--- /dev/null
+++ b/Code/Source/solver/fourier_interpolation.cpp
@@ -0,0 +1,622 @@
+// SPDX-FileCopyrightText: Copyright (c) Stanford University, The Regents of the
+// University of California, and others. SPDX-License-Identifier: BSD-3-Clause
+
+#include "fourier_interpolation.h"
+#include "consts.h"
+
+#include <algorithm>
+#include <cmath>
+#include <fstream>
+#include <regex>
+#include <sstream>
+
+namespace {
+/// Clean a line by removing carriage returns and leading/trailing spaces, and
+/// replacing multiple spaces with a single space.
+std::string clean_line(std::string line) {
+  line.erase(std::remove(line.begin(), line.end(), '\r'), line.end());
+  return std::regex_replace(line, std::regex("^ +| +$|( ) +"), "$1");
+}
+} // namespace
+
+FourierInterpolation FourierInterpolation::from_time_series(
+    const unsigned int n_fourier_coefficients,
+    const std::vector<std::vector<double>> &temporal_values,
+    const bool use_ramp) {
+  const unsigned int n_time_points = temporal_values.size();
+
+  if (n_time_points < 2) {
+    svmp::raise<svmp::FE::InvalidArgumentException>(
+        SVMP_HERE, "At least two time points are needed to construct "
+                   "a FourierInterpolation object.");
+  }
+
+  const unsigned int n_components = temporal_values[0].size() - 1;
+
+  // Check that all entries of temporal_values have the same number of elements.
+  // @todo[michelebucelli] This should probably be an Array, in which this is
+  //   enforced a priori, rather than a vector of vectors. This change would
+  //   also save us the extraction for loop below.
+  for (const auto &row : temporal_values) {
+    if (row.size() != n_components + 1) {
+      svmp::raise<svmp::FE::InvalidArgumentException>(
+          SVMP_HERE, "All rows of temporal_values must have the same number of "
+                     "elements.");
+    }
+  }
+
+  Vector<double> times(n_time_points);
+  Array<double> values(n_components, n_time_points);
+
+  for (unsigned int i = 0; i < n_time_points; ++i) {
+    times[i] = temporal_values[i][0];
+    for (unsigned int j = 0; j < n_components; ++j) {
+      values(j, i) = temporal_values[i][j + 1];
+    }
+  }
+
+  FourierInterpolation result;
+
+  result.use_ramp = use_ramp;
+  result.n_components = n_components;
+  result.n_fourier_coefficients = use_ramp ? 1 : n_fourier_coefficients;
+  result.initial_time = times[0];
+  result.period = times[n_time_points - 1] - times[0];
+  result.linear_trend_initial_values.resize(n_components);
+  result.linear_trend_slopes.resize(n_components);
+  result.fourier_coefficients_real.resize(n_components,
+                                          result.n_fourier_coefficients);
+  result.fourier_coefficients_imaginary.resize(n_components,
+                                               result.n_fourier_coefficients);
+
+  // Compute the linear trend part.
+  for (unsigned int j = 0; j < n_components; ++j) {
+    result.linear_trend_initial_values[j] = values(j, 0);
+    result.linear_trend_slopes[j] =
+        (values(j, n_time_points - 1) - values(j, 0)) / result.period;
+  }
+
+  // Subtract the linear trend part from the values.
+  for (unsigned int i = 0; i < n_time_points; ++i) {
+    times[i] -= result.initial_time;
+    for (unsigned int j = 0; j < n_components; ++j) {
+      values(j, i) = values(j, i) - result.linear_trend_initial_values[j] -
+                     result.linear_trend_slopes[j] * times[i];
+    }
+  }
+
+  // Compute the Fourier coefficients.
+  for (unsigned int n = 0; n < result.n_fourier_coefficients; ++n) {
+    const double tmp = static_cast<double>(n);
+    result.fourier_coefficients_real.set_col(n, 0.0);
+    result.fourier_coefficients_imaginary.set_col(n, 0.0);
+
+    for (unsigned int i = 0; i < n_time_points - 1; ++i) {
+      const double ko = 2.0 * consts::pi * tmp * times[i] / result.period;
+      const double kn = 2.0 * consts::pi * tmp * times[i + 1] / result.period;
+
+      for (unsigned int j = 0; j < result.n_components; j++) {
+        const double s =
+            (values(j, i + 1) - values(j, i)) / (times[i + 1] - times[i]);
+
+        if (n == 0) {
+          result.fourier_coefficients_real(j, n) +=
+              0.5 * (times[i + 1] - times[i]) *
+              (values(j, i + 1) + values(j, i));
+        } else {
+          result.fourier_coefficients_real(j, n) +=
+              s * (std::cos(kn) - std::cos(ko));
+          result.fourier_coefficients_imaginary(j, n) -=
+              s * (std::sin(kn) - std::sin(ko));
+        }
+      }
+    }
+
+    if (n == 0) {
+      for (unsigned int k = 0; k < result.n_components; k++) {
+        result.fourier_coefficients_real(k, n) /= result.period;
+      }
+    } else {
+      const double tmp_2 = (consts::pi * consts::pi * tmp * tmp);
+
+      for (unsigned int k = 0; k < result.n_components; k++) {
+        result.fourier_coefficients_real(k, n) =
+            0.5 * result.fourier_coefficients_real(k, n) * result.period /
+            tmp_2;
+        result.fourier_coefficients_imaginary(k, n) =
+            0.5 * result.fourier_coefficients_imaginary(k, n) * result.period /
+            tmp_2;
+      }
+    }
+  }
+
+  return result;
+}
+
+FourierInterpolation
+FourierInterpolation::from_time_series_file(const std::string &file_name,
+                                            const unsigned int n_components,
+                                            const bool use_ramp) {
+  std::ifstream file(file_name);
+
+  if (!file.is_open()) {
+    svmp::raise<svmp::FileNotFoundException>(SVMP_HERE, file_name);
+  }
+
+  // Read the header of the file.
+  int n_time_points, n_fourier_coefficients;
+  file >> n_time_points >> n_fourier_coefficients;
+
+  if (n_time_points < 2) {
+    svmp::raise<svmp::FileFormatException>(
+        SVMP_HERE, file_name,
+        "At least 2 time points are required to set up Fourier interpolation. "
+        "Only " +
+            std::to_string(n_time_points) + " time points were provided.");
+  }
+
+  if (n_fourier_coefficients == 0) {
+    svmp::raise<svmp::FileFormatException>(
+        SVMP_HERE, file_name,
+        "At least 1 Fourier coefficient is required to set up Fourier "
+        "interpolation. 0 Fourier coefficients were provided.");
+  }
+
+  // Read the time-value pairs.
+  std::vector<std::vector<double>> values;
+  double tmp;
+
+  std::string line;
+  int line_number = 1;
+
+  while (std::getline(file, line)) {
+    line = clean_line(line);
+    if (line.empty())
+      continue;
+
+    std::istringstream line_string_stream(line);
+    std::vector<double> line_values;
+
+    while (!line_string_stream.eof()) {
+      line_string_stream >> tmp;
+
+      if (line_string_stream.fail()) {
+        svmp::raise<svmp::FileFormatException>(
+            SVMP_HERE, file_name,
+            "Error reading values for the temporal values file for line " +
+                std::to_string(line_number) + ": '" + line +
+                "'; value number " + std::to_string(line_values.size() + 1) +
+                " is not a double.");
+      }
+
+      line_values.push_back(tmp);
+    }
+
+    if (line_values.size() != 1 + n_components) {
+      svmp::raise<svmp::FileFormatException>(
+          SVMP_HERE, file_name,
+          "Error reading values for the temporal values file for line " +
+              std::to_string(line_number) + ": '" + line + "'; expected " +
+              std::to_string(1 + n_components) + " values, but got " +
+              std::to_string(line_values.size()) + ".");
+    }
+
+    values.push_back(line_values);
+    ++line_number;
+  }
+
+  return FourierInterpolation::from_time_series(n_fourier_coefficients, values,
+                                                use_ramp);
+}
+
+FourierInterpolation FourierInterpolation::from_fourier_coefficients(
+    const Vector<double> &linear_trend_initial_values,
+    const Vector<double> &linear_trend_slopes,
+    const Array<double> &fourier_coefficients_real,
+    const Array<double> &fourier_coefficients_imaginary,
+    const double initial_time, const double period) {
+  // Linear trend initial values and slopes must have the same size (which will
+  // determine the number of components).
+  if (linear_trend_initial_values.size() != linear_trend_slopes.size()) {
+    svmp::raise<svmp::FE::InvalidArgumentException>(
+        SVMP_HERE,
+        "linear_trend_initial_values and linear_trend_slopes must have the "
+        "same size, but their sizes are " +
+            std::to_string(linear_trend_initial_values.size()) + " and " +
+            std::to_string(linear_trend_slopes.size()) + ".");
+  }
+
+  // The number of rows of the Fourier coefficients must match the number of
+  // components (i.e. the size of the linear trend vectors).
+  if (fourier_coefficients_real.nrows() != linear_trend_initial_values.size()) {
+    svmp::raise<svmp::FE::InvalidArgumentException>(
+        SVMP_HERE, "The number of rows of fourier_coefficients_real must match "
+                   "the size of linear_trend_initial_values.");
+  }
+
+  // The number of rows of the Fourier coefficients must match the number of
+  // components (i.e. the size of the linear trend vectors).
+  if (fourier_coefficients_imaginary.nrows() !=
+      linear_trend_initial_values.size()) {
+    svmp::raise<svmp::FE::InvalidArgumentException>(
+        SVMP_HERE, "The number of rows of fourier_coefficients_imaginary must "
+                   "match the size of linear_trend_initial_values.");
+  }
+
+  // The number of columns of the real and imaginary Fourier coefficients must
+  // match.
+  if (fourier_coefficients_real.ncols() !=
+      fourier_coefficients_imaginary.ncols()) {
+    svmp::raise<svmp::FE::InvalidArgumentException>(
+        SVMP_HERE, "The number of columns of fourier_coefficients_real and "
+                   "fourier_coefficients_imaginary must match.");
+  }
+
+  FourierInterpolation result;
+
+  result.use_ramp = false;
+  result.n_fourier_coefficients = fourier_coefficients_real.ncols();
+  result.n_components = linear_trend_initial_values.size();
+  result.linear_trend_initial_values = linear_trend_initial_values;
+  result.linear_trend_slopes = linear_trend_slopes;
+  result.fourier_coefficients_real = fourier_coefficients_real;
+  result.fourier_coefficients_imaginary = fourier_coefficients_imaginary;
+  result.initial_time = initial_time;
+  result.period = period;
+
+  return result;
+}
+
+FourierInterpolation FourierInterpolation::from_fourier_coefficients_file(
+    const std::string &file_name, const unsigned int n_components) {
+  std::ifstream file(file_name);
+
+  if (!file.is_open()) {
+    svmp::raise<svmp::FileNotFoundException>(SVMP_HERE, file_name);
+  }
+
+  double initial_time, period;
+  file >> initial_time >> period;
+
+  if (file.fail()) {
+    svmp::raise<svmp::FileFormatException>(
+        SVMP_HERE, file_name,
+        "Could not read the initial time and period from the first line.");
+  }
+
+  // Read the linear trend part. Each line is expected to contain the initial
+  // value and slope of the linear trend for one component, so it must hold at
+  // least 2 values.
+  Vector<double> linear_trend_initial_values(n_components);
+  Vector<double> linear_trend_slopes(n_components);
+
+  std::string line;
+  double tmp;
+  unsigned int current_component = 0;
+
+  while (current_component < n_components && std::getline(file, line)) {
+    line = clean_line(line);
+    if (line.empty())
+      continue;
+
+    std::istringstream line_string_stream(line);
+    std::vector<double> values;
+
+    while (line_string_stream >> tmp) {
+      values.push_back(tmp);
+    }
+
+    if (values.size() != 2) {
+      svmp::raise<svmp::FileFormatException>(
+          SVMP_HERE, file_name,
+          "Error reading the linear trend for component " +
+              std::to_string(current_component) + ": '" + line +
+              "'; expected exactly 2 values (initial value and slope), but "
+              "got " +
+              std::to_string(values.size()) + ".");
+    }
+
+    linear_trend_initial_values[current_component] = values[0];
+    linear_trend_slopes[current_component] = values[1];
+    ++current_component;
+  }
+
+  if (current_component != n_components) {
+    svmp::raise<svmp::FileFormatException>(
+        SVMP_HERE, file_name,
+        "Expected " + std::to_string(n_components) +
+            " lines of linear trend coefficients (one per component), but only "
+            "found " +
+            std::to_string(current_component) + ".");
+  }
+
+  // Read the Fourier coefficients.
+  unsigned int n_fourier_coefficients;
+  file >> n_fourier_coefficients;
+
+  if (file.fail()) {
+    svmp::raise<svmp::FileFormatException>(
+        SVMP_HERE, file_name,
+        "Could not read the number of Fourier coefficients.");
+  }
+
+  Array<double> fourier_coefficients_real(n_components, n_fourier_coefficients);
+  Array<double> fourier_coefficients_imaginary(n_components,
+                                               n_fourier_coefficients);
+
+  // Each line is expected to hold the real parts of the coefficients for all
+  // components followed by the imaginary parts, so it must hold exactly
+  // 2 * n_components values.
+  const unsigned int n_values_per_line = 2 * n_components;
+
+  unsigned int current_coefficient = 0;
+  while (current_coefficient < n_fourier_coefficients &&
+         std::getline(file, line)) {
+    line = clean_line(line);
+    if (line.empty())
+      continue;
+
+    std::istringstream line_string_stream(line);
+    std::vector<double> values;
+
+    while (line_string_stream >> tmp) {
+      values.push_back(tmp);
+    }
+
+    if (values.size() != n_values_per_line) {
+      svmp::raise<svmp::FileFormatException>(
+          SVMP_HERE, file_name,
+          "Error reading Fourier coefficient " +
+              std::to_string(current_coefficient) + ": '" + line +
+              "'; expected " + std::to_string(n_values_per_line) +
+              " values (real and imaginary parts for " +
+              std::to_string(n_components) + " components), but got " +
+              std::to_string(values.size()) + ".");
+    }
+
+    for (unsigned int j = 0; j < n_components; ++j) {
+      fourier_coefficients_real(j, current_coefficient) = values[j];
+      fourier_coefficients_imaginary(j, current_coefficient) =
+          values[j + n_components];
+    }
+
+    ++current_coefficient;
+  }
+
+  if (current_coefficient != n_fourier_coefficients) {
+    svmp::raise<svmp::FileFormatException>(
+        SVMP_HERE, file_name,
+        "Expected " + std::to_string(n_fourier_coefficients) +
+            " lines of Fourier coefficients, but only found " +
+            std::to_string(current_coefficient) + ".");
+  }
+
+  return FourierInterpolation::from_fourier_coefficients(
+      linear_trend_initial_values, linear_trend_slopes,
+      fourier_coefficients_real, fourier_coefficients_imaginary, initial_time,
+      period);
+}
+
+void FourierInterpolation::distribute(const CmMod &cm_mod, const cmType &cm) {
+  // Only the master knows whether the object has been initialized. Therefore,
+  // we broadcast the initialization flag to all ranks, so that the following if
+  // statement can be run correctly by all.
+  bool initialized = defined();
+  cm.bcast(cm_mod, &initialized);
+
+  if (initialized) {
+    cm.bcast(cm_mod, &use_ramp);
+    cm.bcast(cm_mod, &n_fourier_coefficients);
+    cm.bcast(cm_mod, &n_components);
+
+    // All ranks but the master need to allocate the arrays before receiving
+    // data.
+    if (cm.slv(cm_mod)) {
+      linear_trend_initial_values.resize(n_components);
+      linear_trend_slopes.resize(n_components);
+      fourier_coefficients_real.resize(n_components, n_fourier_coefficients);
+      fourier_coefficients_imaginary.resize(n_components,
+                                            n_fourier_coefficients);
+    }
+
+    cm.bcast(cm_mod, &initial_time);
+    cm.bcast(cm_mod, &period);
+    cm.bcast(cm_mod, linear_trend_initial_values);
+    cm.bcast(cm_mod, linear_trend_slopes);
+    cm.bcast(cm_mod, fourier_coefficients_real);
+    cm.bcast(cm_mod, fourier_coefficients_imaginary);
+  }
+}
+
+Vector<double> FourierInterpolation::value(const double time) const {
+  Vector<double> result(n_components);
+  static Vector<double> dummy;
+
+  evaluate_internal(time, /* evaluate_derivative = */ false, result, dummy);
+
+  return result;
+}
+
+std::pair<Vector<double>, Vector<double>>
+FourierInterpolation::value_and_derivative(const double time) const {
+  Vector<double> value(n_components);
+  Vector<double> derivative(n_components);
+
+  evaluate_internal(time, /* evaluate_derivative = */ true, value, derivative);
+
+  return std::make_pair(value, derivative);
+}
+
+bool FourierInterpolation::defined() const {
+  return n_fourier_coefficients != 0;
+}
+
+unsigned int FourierInterpolation::get_n_components() const {
+  return n_components;
+}
+
+unsigned int FourierInterpolation::get_n_fourier_coefficients() const {
+  return n_fourier_coefficients;
+}
+
+double FourierInterpolation::get_linear_trend_initial_value(
+    const unsigned int component) const {
+  if (!defined()) {
+    svmp::raise<svmp::FE::NotInitializedException>(
+        SVMP_HERE,
+        "Cannot get linear trend initial value of FourierInterpolation "
+        "instance that has not been defined.");
+  }
+
+  if (component >= n_components) {
+    svmp::raise<svmp::FE::InvalidArgumentException>(
+        SVMP_HERE, "Component index " + std::to_string(component) +
+                       " is out of bounds for FourierInterpolation instance "
+                       "with " +
+                       std::to_string(n_components) + " components.");
+  }
+
+  return linear_trend_initial_values[component];
+}
+
+double FourierInterpolation::get_linear_trend_slope(
+    const unsigned int component) const {
+  if (!defined()) {
+    svmp::raise<svmp::FE::NotInitializedException>(
+        SVMP_HERE,
+        "Cannot get linear trend initial value of FourierInterpolation "
+        "instance that has not been defined.");
+  }
+
+  if (component >= n_components) {
+    svmp::raise<svmp::FE::InvalidArgumentException>(
+        SVMP_HERE, "Component index " + std::to_string(component) +
+                       " is out of bounds for FourierInterpolation instance "
+                       "with " +
+                       std::to_string(n_components) + " components.");
+  }
+
+  return linear_trend_slopes[component];
+}
+
+double
+FourierInterpolation::get_coefficient_real(const unsigned int component,
+                                           const unsigned int frequency) const {
+  if (!defined()) {
+    svmp::raise<svmp::FE::NotInitializedException>(
+        SVMP_HERE,
+        "Cannot get linear trend initial value of FourierInterpolation "
+        "instance that has not been defined.");
+  }
+
+  if (component >= n_components) {
+    svmp::raise<svmp::FE::InvalidArgumentException>(
+        SVMP_HERE, "Component index " + std::to_string(component) +
+                       " is out of bounds for FourierInterpolation instance "
+                       "with " +
+                       std::to_string(n_components) + " components.");
+  }
+
+  if (frequency >= n_fourier_coefficients) {
+    svmp::raise<svmp::FE::InvalidArgumentException>(
+        SVMP_HERE, "Frequency index " + std::to_string(frequency) +
+                       " is out of bounds for FourierInterpolation instance "
+                       "with " +
+                       std::to_string(n_fourier_coefficients) +
+                       " Fourier coefficients.");
+  }
+
+  return fourier_coefficients_real(component, frequency);
+}
+
+double FourierInterpolation::get_coefficient_imaginary(
+    const unsigned int component, const unsigned int frequency) const {
+
+  if (!defined()) {
+    svmp::raise<svmp::FE::NotInitializedException>(
+        SVMP_HERE,
+        "Cannot get linear trend initial value of FourierInterpolation "
+        "instance that has not been defined.");
+  }
+
+  if (component >= n_components) {
+    svmp::raise<svmp::FE::InvalidArgumentException>(
+        SVMP_HERE, "Component index " + std::to_string(component) +
+                       " is out of bounds for FourierInterpolation instance "
+                       "with " +
+                       std::to_string(n_components) + " components.");
+  }
+
+  if (frequency >= n_fourier_coefficients) {
+    svmp::raise<svmp::FE::InvalidArgumentException>(
+        SVMP_HERE, "Frequency index " + std::to_string(frequency) +
+                       " is out of bounds for FourierInterpolation instance "
+                       "with " +
+                       std::to_string(n_fourier_coefficients) +
+                       " Fourier coefficients.");
+  }
+
+  return fourier_coefficients_imaginary(component, frequency);
+}
+
+void FourierInterpolation::evaluate_internal(const double time,
+                                             const bool evaluate_derivative,
+                                             Vector<double> &value,
+                                             Vector<double> &derivative) const {
+  if (!defined()) {
+    svmp::raise<svmp::FE::NotInitializedException>(
+        SVMP_HERE, "Cannot evaluate FourierInterpolation instance that has not "
+                   "been defined.");
+  }
+
+  // Shifted and rescaled time.
+  // The input time is shifted by ti. Then, if using the ramp function, it is
+  // clamped to the interval [0, T]. Otherwise, the time is wrapped to the
+  // interval [0, T], to enable periodicity.
+  const double t = use_ramp
+                       ? std::max(std::min(time - initial_time, period), 0.0)
+                       : std::fmod(time - initial_time, period);
+
+  // Linear trend.
+  for (unsigned int i = 0; i < n_components; ++i) {
+    value[i] = linear_trend_initial_values[i] + t * linear_trend_slopes[i];
+
+    // @todo[michelebucelli] This is consistent with the old code, but is
+    //   incorrect when use_ramp = true. In that case, the derivative should be
+    //   zero outside the interpolation interval [initial_time, initial_time +
+    //   period].
+    if (evaluate_derivative)
+      derivative[i] = linear_trend_slopes[i];
+  }
+
+  // Fourier series.
+  if (!use_ramp) {
+    const double tmp = 2.0 * consts::pi / period;
+
+    // Fourier series.
+    for (int i = 0; i < n_fourier_coefficients; ++i) {
+      const double dk = tmp * i;
+      const double K = t * dk;
+
+      for (int j = 0; j < n_components; ++j) {
+        // Using value[j] = value[j] + ... instead of value[j] += ..., because
+        // the latter changes the order of operations enough to break some of
+        // the tests.
+        // @todo[michelebucelli] This seems pretty fragile! It happens in a few
+        //   other places in this file as well, where using an increment
+        //   operator (*= or +=) changes the order of operations resulting in
+        //   changes in the test values. This should be investigated and the
+        //   best (i.e. most accurate) operation order should be chosen.
+        value[j] = value[j] + fourier_coefficients_real(j, i) * std::cos(K) -
+                   fourier_coefficients_imaginary(j, i) * std::sin(K);
+
+        if (evaluate_derivative) {
+          derivative[j] -=
+              (fourier_coefficients_real(j, i) * std::sin(K) +
+               fourier_coefficients_imaginary(j, i) * std::cos(K)) *
+              dk;
+        }
+      }
+    }
+  }
+}
\ No newline at end of file
diff --git a/Code/Source/solver/fourier_interpolation.h b/Code/Source/solver/fourier_interpolation.h
new file mode 100644
index 000000000..cf898cf66
--- /dev/null
+++ b/Code/Source/solver/fourier_interpolation.h
@@ -0,0 +1,397 @@
+// SPDX-FileCopyrightText: Copyright (c) Stanford University, The Regents of the
+// University of California, and others. SPDX-License-Identifier: BSD-3-Clause
+
+#ifndef FOURIER_INTERPOLATION_H
+#define FOURIER_INTERPOLATION_H
+
+#include "Array.h"
+#include "CmMod.h"
+#include "Vector.h"
+
+#include "Core/Exception.h"
+#include "FE/Common/FEException.h"
+
+#include <string>
+#include <utility>
+#include <vector>
+
+/**
+ * @brief Fourier interpolation of time dependent data.
+ *
+ * This class implements interpolation of time dependent data through either a
+ * Fourier series or a linear clamped ramp. Which of the two is used is
+ * determined by the boolean member variable @ref use_ramp.
+ *
+ * In what follows, let @f$\{t_i, \mathbf{v}_i\}_{i=0}^{N-1}@f$ be the
+ * interpolated time series, and @f$T = t_{N-1} - t_0@f$ be the period of the
+ * time series. The time series is assumed to be ordered, i.e. @f$t_i <
+ * t_{i+1}@f$ for all @f$i@f$.
+ *
+ * ## Using this class
+ *
+ * The FourierInterpolation class is not meant to be constructed directly, but
+ * rather through one of the static methods @ref from_time_series, @ref
+ * from_time_series_file, @ref from_fourier_coefficients, or @ref
+ * from_fourier_coefficients_file.
+ *
+ * If data is only read from the master rank in a parallel setting, the method
+ * @ref distribute can be used to broadcast the data to all ranks.
+ *
+ * Once the object has been constructed, the interpolated value at a given time
+ * can be obtained through the method @ref value, and the interpolated value and
+ * its time derivative can be obtained through the method @ref
+ * value_and_derivative.
+ *
+ * ## Fourier series interpolation
+ *
+ * This gives rise to a periodic interpolation, with period @f$T@f$. Let @f$M@f$
+ * be the user-defined number of Fourier modes. The interpolated value at time
+ * @f$t@f$ is given by
+ * @f[ \begin{aligned}
+ * \tilde{\mathbf{v}}(t) &=
+ *   \underbrace{\mathbf{q}_i + \mathbf{q}_s \tau(t)}_{\text{linear trend}} +
+ *   \underbrace{
+ *     \sum_{k=0}^{M-1} \left(
+ *       \mathbf{c}_k^{\text{real}} \cos\left(\frac{2 \pi k \tau(t)}{T}\right) -
+ *       \mathbf{c}_k^{\text{imag}} \sin\left(\frac{2 \pi k \tau(t)}{T}\right)
+ *     \right)
+ *   }_{\text{Fourier series}} \\
+ *   \tau(t) &= (t - t_0)\;\text{mod}\;T
+ * \end{aligned} @f]
+ * The coefficients of the linear trend and Fourier series are determined from
+ * the input time series as follows:
+ * @f[ \begin{aligned}
+ * \mathbf{q}_i &= \mathbf{v}_0 \\
+ * \mathbf{q}_s &= \frac{\mathbf{v}_{N-1} - \mathbf{v}_0}{T} \\
+ * \mathbf{c}_0^{\text{real}} &=
+ *   \frac{1}{2T} \sum_{i=0}^{N-2} (\hat{t}_{i+1} - \hat{t}_i)
+ *     (\hat{\mathbf{v}}_{i+1} + \hat{\mathbf{v}}_i) \\
+ * \mathbf{c}_0^{\text{imag}} &= 0 \\
+ * \mathbf{c}_k^{\text{real}} &= \frac{T}{2 \pi^2 k^2}
+ *   \sum_{i=0}^{N-2} \frac{\hat{\mathbf{v}}_{i+1} - \hat{\mathbf{v}}_i}
+ *   {\hat{t}_{i+1} - \hat{t}_i} \left(
+ *     \cos\left(\frac{2 \pi k \hat{t}_{i+1}}{T}\right) -
+ *     \cos\left(\frac{2 \pi k \hat{t}_i}{T}\right)
+ *   \right) \\
+ * \mathbf{c}_k^{\text{imag}} &= \frac{T}{2 \pi^2 k^2}
+ *   -\sum_{i=0}^{N-2} \frac{\hat{\mathbf{v}}_{i+1} - \hat{\mathbf{v}}_i}
+ *   {\hat{t}_{i+1} - \hat{t}_i} \left(
+ *     \sin\left(\frac{2 \pi k \hat{t}_{i+1}}{T}\right) -
+ *     \sin\left(\frac{2 \pi k \hat{t}_i}{T}\right)
+ *   \right)
+ * \end{aligned} @f]
+ * The quantities @f$\hat{t}_i@f$ and @f$\hat{\mathbf{v}}_i@f$ are the time and
+ * value series after subtracting the linear trend:
+ * @f[ \begin{aligned}
+ * \hat{t}_i &= t_i - t_0 \\
+ * \hat{\mathbf{v}}_i &= \mathbf{v}_i - \mathbf{q}_i - \mathbf{q}_s \hat{t}_i
+ * \end{aligned} @f]
+ *
+ * ## Ramp interpolation
+ *
+ * The interpolated value is equal to @f$\mathbf{v}_0@f$ until time @f$t_0@f$,
+ * then it follows a linear ramp from @f$\mathbf{v}_0@f$ to
+ * @f$\mathbf{v}_{N-1}@f$ at @f$t_{N-1}@f$, and then it remains constant at
+ * @f$\mathbf{v}_{N-1}@f$ for all times greater than @f$t_{N-1}@f$. Notice in
+ * particular that this interpolation is not periodic.
+ *
+ * The interpolated value is given by
+ * @f[
+ * \tilde{\mathbf{v}}(t) = \begin{cases}
+ *   \mathbf{v}_0 & t < t_0 \\
+ *   \mathbf{v}_0 + \frac{\mathbf{v}_{N-1}
+ *     - \mathbf{v}_0}{t_{N-1} - t_0} (t - t_0)
+ *     & t_0 \leq t < t_{N-1} \\
+ *   \mathbf{v}_{N-1} & t \geq t_{N-1}
+ * \end{cases}
+ * @f]
+ *
+ * This is equivalent to only taking the linear trend part of the Fourier series
+ * interpolation, and clamping the time to the interval @f$[t_0, t_{N-1}]@f$.
+ *
+ */
+class FourierInterpolation {
+public:
+  /**
+   * @brief Default constructor.
+   *
+   * This constructor is not meant to be used directly, and is only here to
+   * facilitate storing objects of this type in STL containers. The object
+   * constructed this way is not initialized and will not be usable until it is
+   * assigned to a valid FourierInterpolation instance. Use the static methods
+   * @ref from_time_series, @ref from_time_series_file, @ref
+   * from_fourier_coefficients, or @ref from_fourier_coefficients_file to
+   * construct a valid instance.
+   */
+  FourierInterpolation() = default;
+
+  /**
+   * @brief Construct a FourierInterpolation from a time series.
+   *
+   * @param[in] n_fourier_coefficients The number @f$M@f$ of Fourier modes to
+   *   use in the interpolation.
+   * @param[in] temporal_values The time series to interpolate. This is a vector
+   *   of vectors, where each inner vector contains the time and the values of
+   *   all components at that time. In other words,
+   *   <kbd>temporal_values[i][0]</kbd> is the time @f$t_i@f$, and
+   *   <kbd>temporal_values[i][j]</kbd> is the value the @f$j@f$-th component of
+   *   @f$\mathbf{v}_i@f$.
+   * @param[in] use_ramp Whether to use a ramp function for the interpolation.
+   *   See the general class documentation for the precise meaning of this
+   *   choice.
+   */
+  static FourierInterpolation
+  from_time_series(const unsigned int n_fourier_coefficients,
+                   const std::vector<std::vector<double>> &temporal_values,
+                   const bool use_ramp);
+
+  /**
+   * @brief Read a time series from file and return the corresponding instance
+   * of FourierInterpolation.
+   *
+   * The input file is expected to have the following format:
+   * ```
+   * <number of time points> <number of Fourier coefficients>
+   * <time 0> <value 0, component 0> ... <value 0, component d-1>
+   * <time 1> <value 1, component 0> ... <value 1, component d-1>
+   * ...
+   * <time N-1> <value N-1, component 0> ... <value N-1, component d-1>
+   * ```
+   *
+   * @param[in] file_name The name of the file to read. An exception will be
+   *   thrown if the file cannot be opened.
+   * @param[in] n_components The number of components of the data to be
+   *   interpolated. If any row in the file does not have this number of entries
+   *   (plus one entry for the time), an exception will be thrown.
+   * @param[in] use_ramp Whether to use a ramp function for the interpolation.
+   *   See the general class documentation for the precise meaning of this
+   *   choice. If this parameter is set to true, then the number of Fourier
+   *   coefficients in the file will be ignored, and the time points except the
+   *   first and last will have no effect.
+   */
+  static FourierInterpolation
+  from_time_series_file(const std::string &file_name,
+                        const unsigned int n_components, const bool use_ramp);
+
+  /**
+   * @brief Construct a FourierInterpolation from Fourier coefficients.
+   *
+   * This method bypasses the computation of the Fourier coefficients from a
+   * time series, and construct the instance directly from precomputed
+   * coefficients.
+   *
+   * This construction method always sets @ref use_ramp to false, and therefore
+   * always constructs a periodic interpolation.
+   *
+   * @param[in] linear_trend_initial_values The initial values
+   *   @f$\mathbf{q}_i@f$ of the linear trend part of the interpolation, with
+   *   one entry for each component.
+   * @param[in] linear_trend_slopes The slopes @f$\mathbf{q}_s@f$ of the linear
+   *   trend part of the interpolation, with one entry for each component.
+   * @param[in] fourier_coefficients_real The real part
+   *   @f$\mathbf{c}_k^\text{real}@f$ of the Fourier interpolation. This is a 2D
+   *   array, for which the first index selects the component and the second
+   *   index selects the Fourier mode.
+   * @param[in] fourier_coefficients_imaginary The imaginary part
+   *   @f$\mathbf{c}_k^\text{imag}@f$ of the Fourier interpolation. This is a 2D
+   *   array, for which the first index selects the component and the second
+   *   index selects the Fourier mode.
+   * @param[in] initial_time The initial time @f$t_0@f$ of the interpolation.
+   * @param[in] period The period @f$T@f$ of the interpolation.
+   */
+  static FourierInterpolation
+  from_fourier_coefficients(const Vector<double> &linear_trend_initial_values,
+                            const Vector<double> &linear_trend_slopes,
+                            const Array<double> &fourier_coefficients_real,
+                            const Array<double> &fourier_coefficients_imaginary,
+                            const double initial_time, const double period);
+
+  /**
+   * @brief Read Fourier coefficients from file and return the corresponding
+   * FourierInterpolation instance.
+   *
+   * The input file is expected to have the following format:
+   * ```
+   * <initial time> <period>
+   * <q_i, component 0> <q_s, component 0>
+   * ...
+   * <q_i, component d-1> <q_s, component d-1>
+   * <number of Fourier coefficients>
+   * <c_r[0][0]> <c_r[1][0]> ... <c_r[d-1][0]> <c_i[0][0]> <c_i[1][0]> ... <c_i[d-1][0]>
+   * <c_r[0][1]> <c_r[1][1]> ... <c_r[d-1][1]> <c_i[0][1]> <c_i[1][1]> ... <c_i[d-1][1]>
+   * ...
+   * <c_r[0][M-1]> <c_r[1][M-1]> ... <c_r[d-1][M-1]> <c_i[0][M-1]> <c_i[1][M-1]> ... <c_i[d-1][M-1]>
+   * ```
+   * where <kbd>c_r</kbd> and <kbd>c_i</kbd> are the real and imaginary parts of
+   * the Fourier coefficients.
+   *
+   * @param[in] file_name The name of the file to read. An exception will be
+   *   thrown if the file cannot be opened.
+   * @param[in] n_components The number of components of the data to be
+   *   interpolated. This is used to check the correctness of the file format,
+   *   and an exception will be thrown if the check fails.
+   */
+  static FourierInterpolation
+  from_fourier_coefficients_file(const std::string &file_name,
+                                 const unsigned int n_components);
+
+  /**
+   * @brief Distribute the data to all parallel processes.
+   *
+   * Broadcasts the data contained in this object from the master rank to all
+   * other ranks. This is necessary if only the master rank initializes the
+   * object (e.g. by reading its data from a file), but all ranks need to use
+   * it.
+   *
+   * @param[in] cm_mod The communication module to use for the broadcast.
+   * @param[in] cm The communicator to use for the broadcast.
+   */
+  void distribute(const CmMod &cm_mod, const cmType &cm);
+
+  /**
+   * @brief Return the interpolated value at a given time.
+   *
+   * Refer to the general class documentation for details on how the returned
+   * value is computed.
+   *
+   * @param[in] time The time at which to evaluate the interpolation.
+   */
+  Vector<double> value(const double time) const;
+
+  /**
+   * @brief Return the interpolated value and its time derivative at a given
+   * time.
+   *
+   * Refer to the general class documentation for details on how the returned
+   * value is computed.
+   *
+   * @param[in] time The time at which to evaluate the interpolation.
+   *
+   * @return A pair in which the first element is the interpolated value and the
+   *   second element is its time derivative.
+   */
+  std::pair<Vector<double>, Vector<double>>
+  value_and_derivative(const double time) const;
+
+  /// @name Data member access.
+  /// @{
+
+  /// @brief Return whether this object has been initialized.
+  bool defined() const;
+
+  /// @brief Get the dimension of the data interpolated by this object.
+  unsigned int get_n_components() const;
+
+  /// @brief Get the number of Fourier coefficients used by this object.
+  unsigned int get_n_fourier_coefficients() const;
+
+  /// @brief Get the initial value of the linear trend part for one component.
+  double get_linear_trend_initial_value(const unsigned int component) const;
+
+  /// @brief Get the slope of the linear trend part for one component.
+  double get_linear_trend_slope(const unsigned int component) const;
+
+  /// @brief Get the real part of the Fourier coefficients for one component.
+  double get_coefficient_real(const unsigned int component,
+                              const unsigned int frequency) const;
+
+  /// @brief Get the imaginary part of the Fourier coefficients for one
+  /// component.
+  double get_coefficient_imaginary(const unsigned int component,
+                                   const unsigned int frequency) const;
+
+  /// @}
+
+private:
+  /** @brief Internal evaluation function.
+   *
+   * Uses the inverse Fourier transform to evaluate the value, and optionally
+   * the derivative, of the interpolated data. This function is not meant to be
+   * used directly, but only as a backend to @ref value and @ref
+   * value_and_derivative.
+   *
+   * The vectors value and derivative are assumed to be of size @ref d. This is
+   * not checked by this function.
+   *
+   * @throws svmp::FE::NotInitializedException if this FourierInterpolation
+   * instance has not been initialized (i.e. if @ref defined returns false).
+   *
+   * @param[in] time The time at which to evaluate the interpolation.
+   * @param[in] evaluate_derivative Whether to also evaluate the time
+   *   derivative of the interpolation.
+   * @param[out] value The interpolated value at the given time.
+   * @param[out] derivative The time derivative of the interpolated value at
+   *   the given time. If evaluated_derivative is false, this will not be
+   *   modified or accessed.
+   */
+  void evaluate_internal(const double time, const bool evaluate_derivative,
+                         Vector<double> &value,
+                         Vector<double> &derivative) const;
+
+  /**
+   * @brief Toggle whether this is a ramp function or not.
+   *
+   * See the general class documentation for details on the difference between
+   * ramp and periodic interpolation.
+   */
+  bool use_ramp = false;
+
+  /**
+   * @brief Number of Fourier coefficients.
+   */
+  unsigned int n_fourier_coefficients = 0;
+
+  /**
+   * @brief Number of components of the interpolated data.
+   */
+  unsigned int n_components = 0;
+
+  /**
+   * @brief Initial value for the linear trend.
+   *
+   * This is a vector with n_components entries, where each entry is the initial
+   * value (i.e. the value for time = ti) of the linear trend part of the
+   * interpolated data for that component.
+   */
+  Vector<double> linear_trend_initial_values;
+
+  /**
+   * @brief Time derivative for the linear trend.
+   *
+   * This is a vector with n_components entries, where each entry is the slope
+   * of the linear trend part of the interpolated data for that component.
+   */
+  Vector<double> linear_trend_slopes;
+
+  /**
+   * @brief Period of the interpolated data.
+   *
+   * This is disregarded if use_ramp is true. See the general class
+   * documentation for details on the difference between ramp and periodic
+   * interpolation.
+   */
+  double period = 0.0;
+
+  /**
+   * @brief Initial time.
+   */
+  double initial_time = 0.0;
+
+  /**
+   * @brief Real part of the Fourier series coefficients.
+   *
+   * This is a 2D array with n_components rows and n_fourier_coefficients
+   * columns.
+   */
+  Array<double> fourier_coefficients_real;
+
+  /**
+   * @brief Imaginary part of the Fourier series coefficients.
+   *
+   * This is a 2D array with n_components rows and n_fourier_coefficients
+   * columns.
+   */
+  Array<double> fourier_coefficients_imaginary;
+};
+
+#endif
\ No newline at end of file
diff --git a/Code/Source/solver/mat_models.cpp b/Code/Source/solver/mat_models.cpp
index dd8a608ce..58658ac54 100644
--- a/Code/Source/solver/mat_models.cpp
+++ b/Code/Source/solver/mat_models.cpp
@@ -216,10 +216,8 @@ void compute_fib_stress(const ComMod& com_mod, const CepMod& cep_mod, const fibS
 
   if (utils::btest(Tfl.fType, iBC_std)) {
     g = Tfl.g;
-  } else if (utils::btest(Tfl.fType, iBC_ustd)) { 
-    Vector<double> gv(1), tv(1);
-    ifft(com_mod, Tfl.gt, gv, tv);
-    g = gv[0];
+  } else if (utils::btest(Tfl.fType, iBC_ustd)) {
+    g = Tfl.gt.value(com_mod.time)[0];
   }
 }
 
diff --git a/Code/Source/solver/read_files.cpp b/Code/Source/solver/read_files.cpp
index 6c6723bce..4c70dbce3 100644
--- a/Code/Source/solver/read_files.cpp
+++ b/Code/Source/solver/read_files.cpp
@@ -222,27 +222,26 @@ void read_bc(Simulation* simulation, EquationParameters* eq_params, eqType& lEq,
       lBc.g = bc_params->value.value();
     }
 
-  } else if (ctmp == "Unsteady") { 
-    lBc.bType = utils::ibset(lBc.bType, enum_int(BoundaryConditionType::bType_ustd)); 
-
-    if (utils::btest(lBc.bType, enum_int(BoundaryConditionType::bType_trac))) { 
-      lBc.gt.d = com_mod.nsd;
-    } else { 
-      lBc.gt.d = 1;
-    }
-
-    if (bc_params->temporal_values_file_path.defined()) { 
-      lBc.gt.lrmp = bc_params->ramp_function.value();
-      auto file_name = bc_params->temporal_values_file_path.value();
-      read_temporal_values(file_name, lBc);
-
-    } else { 
-      if (!bc_params->fourier_coefficients_file_path.defined()) { 
-        throw std::runtime_error("[read_bc] Undefined data for boundary condition type '" + bc_type + "'.");
-      }
-      lBc.gt.lrmp = false;
-      auto file_name = bc_params->fourier_coefficients_file_path.value(); 
-      read_fourier_coeff_values_file(file_name, lBc);
+  } else if (ctmp == "Unsteady") {
+    lBc.bType =
+        utils::ibset(lBc.bType, enum_int(BoundaryConditionType::bType_ustd));
+
+    const unsigned int n_dimensions =
+        utils::btest(lBc.bType, enum_int(BoundaryConditionType::bType_trac))
+            ? com_mod.nsd
+            : 1;
+
+    if (bc_params->temporal_values_file_path.defined()) {
+      lBc.gt = FourierInterpolation::from_time_series_file(
+          bc_params->temporal_values_file_path.value(), n_dimensions,
+          bc_params->ramp_function.value());
+    } else if (bc_params->fourier_coefficients_file_path.defined()) {
+      lBc.gt = FourierInterpolation::from_fourier_coefficients_file(
+          bc_params->fourier_coefficients_file_path.value(), n_dimensions);
+    } else {
+      throw std::runtime_error(
+          "[read_bc] Undefined data for boundary condition type '" + bc_type +
+          "'.");
     }
 
   // Coupling to a 0D model:
@@ -291,7 +290,7 @@ void read_bc(Simulation* simulation, EquationParameters* eq_params, eqType& lEq,
       }
     }
 
-  } else if (ctmp == "Resistance") { 
+  } else if (ctmp == "Resistance") {
     lBc.bType = utils::ibset(lBc.bType, enum_int(BoundaryConditionType::bType_res)); 
     if (!utils::btest(lBc.bType, enum_int(BoundaryConditionType::bType_Neu))) { 
       throw std::runtime_error("[read_bc] Resistance is only defined for Neu BC.");
@@ -337,7 +336,7 @@ void read_bc(Simulation* simulation, EquationParameters* eq_params, eqType& lEq,
   //   S p a t i a l  ////
   ////////////////////////
 
-  } else if (ctmp == "Spatial") { 
+  } else if (ctmp == "Spatial") {
     lBc.bType = utils::ibset(lBc.bType, enum_int(BoundaryConditionType::bType_gen)); 
     if (!utils::btest(lBc.bType, enum_int(BoundaryConditionType::bType_Neu))) { 
       throw std::runtime_error("[read_bc] Spatial BC is only defined for Neu BC.");
@@ -361,7 +360,7 @@ void read_bc(Simulation* simulation, EquationParameters* eq_params, eqType& lEq,
   //   G e n e r a l  ////
   ////////////////////////
 
-  } else if (ctmp == "General") { 
+  } else if (ctmp == "General") {
     lBc.bType = utils::ibset(lBc.bType, enum_int(BoundaryConditionType::bType_gen)); 
     int iM = lBc.iM;
     int iFa = lBc.iFa;
@@ -947,16 +946,18 @@ void read_bf(ComMod& com_mod, BodyForceParameters* bf_params, bfType& lBf)
   // Unsteady //
 
   } else if (time_dependence == "unsteady") {
-    lBf.bType = utils::ibset(lBf.bType, enum_int(BodyForceType::bfType_ustd)); 
-    lBf.bt.d = lBf.dof;
+    lBf.bType = utils::ibset(lBf.bType, enum_int(BodyForceType::bfType_ustd));
+
     if (bf_params->temporal_values_file_path.defined()) {
-      lBf.bt.lrmp = bf_params->ramp_function.value();
-      auto file_name = bf_params->temporal_values_file_path.value();
-      read_temporal_values(file_name, lBf);
+      lBf.bt = FourierInterpolation::from_time_series_file(
+          bf_params->temporal_values_file_path.value(), lBf.dof,
+          bf_params->ramp_function.value());
+    } else if (bf_params->fourier_coefficients_file_path.defined()) {
+      lBf.bt = FourierInterpolation::from_fourier_coefficients_file(
+          bf_params->fourier_coefficients_file_path.value(), lBf.dof);
     } else {
-      lBf.bt.lrmp = false;
-      auto fTmp = bf_params->fourier_coefficients_file_path.value();
-      read_fourier_coeff_values_file(fTmp, lBf);
+      throw std::runtime_error("No temporal values or Fourier coefficients "
+                               "file provided for unsteady body force.");
     }
 
   // Spatial //
@@ -1580,62 +1581,6 @@ void read_eq(Simulation* simulation, EquationParameters* eq_params, eqType& lEq)
   }
 }
 
-//---------------------------------
-// read_fiber_temporal_values_file
-//---------------------------------
-// Read fiber reinforcement stress the values from a file.
-//
-// Note: There is no equivalent Fortran function.
-//
-void read_fiber_temporal_values_file(FiberReinforcementStressParameters& fiber_params, dmnType& lDmn)
-{
-  auto file_name = fiber_params.temporal_values_file_path.value();
-  std::ifstream temporal_values_file;
-  temporal_values_file.open(file_name);
-  if (!temporal_values_file.is_open()) {
-    throw std::runtime_error("Failed to open the fiber reinforcement stress the values file '" + file_name + "'.");
-  }
-
-  int i, j;
-  temporal_values_file >> i >> j; 
-  if ((i == 0) || (j == 0) || (i < 2)) {
-    throw std::runtime_error("Error reading the first line of the fiber reinforcement stress the values file '" + file_name + "' has an incorrect format.");
-  }
-
-  lDmn.stM.Tf.gt.d = 1;
-  lDmn.stM.Tf.gt.n = j;
-
-  // Read time/value pairs.
-  //
-  std::vector<std::vector<double>> temporal_values;
-  double time, value;
-  std::string line;
-
-  while (std::getline(temporal_values_file, line)) { 
-    line.erase(std::remove(line.begin(), line.end(), '\r'), line.end());
-    if (line == "") {
-      continue;
-    }
-    std::istringstream line_input(line);
-    std::vector<double> values;
-    while (line_input >> value) {
-      values.push_back(value);
-    }
-    temporal_values.push_back(values);
-  }
-
-  if (lDmn.stM.Tf.gt.lrmp) {
-    lDmn.stM.Tf.gt.n = 1;
-  }
-
-  lDmn.stM.Tf.gt.qi.resize(lDmn.stM.Tf.gt.d);
-  lDmn.stM.Tf.gt.qs.resize(lDmn.stM.Tf.gt.d);
-  lDmn.stM.Tf.gt.r.resize(lDmn.stM.Tf.gt.d,lDmn.stM.Tf.gt.n);
-  lDmn.stM.Tf.gt.i.resize(lDmn.stM.Tf.gt.d,lDmn.stM.Tf.gt.n);
-
-  fft(i, temporal_values, lDmn.stM.Tf.gt);
-}
-
 //------------
 // read_files
 //------------
@@ -1890,151 +1835,6 @@ void read_files(Simulation* simulation, const std::string& file_name)
   #endif
 }
 
-//--------------------------------
-// read_fourier_coeff_values_file 
-//--------------------------------
-// Set boundary condition Fourier coefficients read in from a file. 
-//
-// Note: There is no equivalent Fortran function.
-//
-// [TODO:DaveP] this is not fully implemented. 
-//
-void read_fourier_coeff_values_file(const std::string& file_name, bcType& lBc) 
-{
-  std::ifstream temporal_values_file;
-  temporal_values_file.open(file_name);
-  if (!temporal_values_file.is_open()) {
-    throw std::runtime_error("Failed to open the Fourier coefficients values file '" + file_name + "'.");
-  }
-
-  temporal_values_file >> lBc.gt.ti >> lBc.gt.T; 
-
-  // Read time/value pairs.
-  //
-  lBc.gt.qi.resize(lBc.gt.d); 
-  lBc.gt.qs.resize(lBc.gt.d);
-
-  double time, value;
-  std::string line;
-  int n = 0;
-
-  while (std::getline(temporal_values_file, line)) { 
-    line.erase(std::remove(line.begin(), line.end(), '\r'), line.end());
-    if (line == "") {
-      continue;
-    }
-    std::istringstream line_input(line);
-    std::vector<double> values;
-    while (line_input >> value) {
-      values.push_back(value);
-    }
-    lBc.gt.qi[n] = values[0];
-    lBc.gt.qs[n] = values[1];
-    n += 1;
-
-    if (n == lBc.gt.d) {
-      break;
-    }
-  }
-
-  temporal_values_file >> lBc.gt.n; 
-
-  lBc.gt.r.resize(lBc.gt.d, lBc.gt.n);
-  lBc.gt.i.resize(lBc.gt.d, lBc.gt.n);
-  
-  int j = 0;
-
-  while (std::getline(temporal_values_file, line)) { 
-    line.erase(std::remove(line.begin(), line.end(), '\r'), line.end());
-    if (line == "") {
-      continue;
-    }
-    std::istringstream line_input(line);
-    std::vector<double> values;
-    while (line_input >> value) {
-      values.push_back(value);
-    }
-
-    for (int i = 0; i < lBc.gt.d; i++) { 
-      lBc.gt.r(i,j) = values[i];
-      lBc.gt.i(i,j) = values[i + lBc.gt.d];
-    }
-
-    j += 1;
-    if (j == lBc.gt.n) {
-      break;
-    }
-  }
-}
-
-void read_fourier_coeff_values_file(const std::string& file_name, bfType& lBf) 
-{
-  std::ifstream temporal_values_file;
-  temporal_values_file.open(file_name);
-  if (!temporal_values_file.is_open()) {
-    throw std::runtime_error("Failed to open the Fourier coefficients values file '" + file_name + "'.");
-  }
-
-  temporal_values_file >> lBf.bt.ti >> lBf.bt.T; 
-
-  // Read time/value pairs.
-  //
-  lBf.bt.qi.resize(lBf.bt.d); 
-  lBf.bt.qs.resize(lBf.bt.d);
-
-  double time, value;
-  std::string line;
-  int n = 0;
-
-  while (std::getline(temporal_values_file, line)) { 
-    line.erase(std::remove(line.begin(), line.end(), '\r'), line.end());
-    if (line == "") {
-      continue;
-    }
-    std::istringstream line_input(line);
-    std::vector<double> values;
-    while (line_input >> value) {
-      values.push_back(value);
-    }
-    lBf.bt.qi[n] = values[0];
-    lBf.bt.qs[n] = values[1];
-    n += 1;
-
-    if (n == lBf.bt.d) {
-      break;
-    }
-  }
-
-  temporal_values_file >> lBf.bt.n; 
-
-  lBf.bt.r.resize(lBf.bt.d, lBf.bt.n);
-  lBf.bt.i.resize(lBf.bt.d, lBf.bt.n);
-  
-  int j = 0;
-
-  while (std::getline(temporal_values_file, line)) { 
-    line.erase(std::remove(line.begin(), line.end(), '\r'), line.end());
-    if (line == "") {
-      continue;
-    }
-    std::istringstream line_input(line);
-    std::vector<double> values;
-    while (line_input >> value) {
-      values.push_back(value);
-    }
-
-    for (int i = 0; i < lBf.bt.d; i++) { 
-      lBf.bt.r(i,j) = values[i];
-      lBf.bt.i(i,j) = values[i + lBf.bt.d];
-    }
-
-    j += 1;
-    if (j == lBf.bt.n) {
-      break;
-    }
-  }
-}
-
 //---------
 // read_ls
 //---------
@@ -2230,9 +2030,13 @@ void read_mat_model(Simulation* simulation, EquationParameters* eq_params, Domai
       lDmn.stM.Tf.g = fiber_params.value.value();
 
     } else if (fiber_stress == "unsteady") {
-      lDmn.stM.Tf.fType = utils::ibset(lDmn.stM.Tf.fType, static_cast<int>(BoundaryConditionType::bType_ustd));
-      lDmn.stM.Tf.gt.lrmp = fiber_params.ramp_function.value();
-      read_fiber_temporal_values_file(fiber_params, lDmn);
+      lDmn.stM.Tf.fType =
+          utils::ibset(lDmn.stM.Tf.fType,
+                       static_cast<int>(BoundaryConditionType::bType_ustd));
+
+      lDmn.stM.Tf.gt = FourierInterpolation::from_time_series_file(
+          fiber_params.temporal_values_file_path.value(),
+          /* n_dimensions = */ 1, fiber_params.ramp_function.value());
     }
 
     // Read directional stress distribution parameters
@@ -2664,134 +2468,6 @@ void read_temp_spat_values(const ComMod& com_mod, const mshType& msh, const std:
   } 
 }
 
-//----------------------
-// read_temporal_values
-//----------------------
-// Set boundary condition temporal values read in from a file. 
-//
-// Data modified:
-//
-//  Set lBc.gt Fourier coefficients data (fcType) 
-//     lBc.gt.n 
-//     lBc.gt.qi
-//     lBc.gt.qs
-//     lBc.gt.r
-//     lBc.gt.i
-//
-// Note: There is no equivalent Fortran function.
-//
-void read_temporal_values(const std::string& file_name, bcType& lBc) 
-{
-  std::ifstream temporal_values_file;
-  temporal_values_file.open(file_name);
-  if (!temporal_values_file.is_open()) {
-    throw std::runtime_error("Failed to open the temporal values file '" + file_name + "'.");
-  }
-
-  int i, j;
-  temporal_values_file >> i >> j; 
-  if ((i == 0) || (j == 0) || (i < 2)) {
-    throw std::runtime_error("Error reading the first line of the temporal values file '" + file_name + "'.");
-  }
-  lBc.gt.n = j;
-
-  // Read time/value pairs.
-  //
-  std::vector<std::vector<double>> temporal_values;
-  double time, value;
-  std::string line;
-  int line_number = 1;
-  int num_values_per_line = lBc.gt.d + 1;
-
-  while (std::getline(temporal_values_file, line)) { 
-    line.erase(std::remove(line.begin(), line.end(), '\r'), line.end());
-    if (line == "") {
-      continue;
-    }
-    // Remove leading and trailing spaces.
-    auto cleaned_line = std::regex_replace(line, std::regex("^ +| +$|( ) +"), "$1");
-    std::istringstream line_input(cleaned_line);
-    std::vector<double> values;
-
-    while (!line_input.eof()) {
-      line_input >> value;
-
-      if (line_input.fail()) { 
-        throw std::runtime_error("Error reading values for the temporal values file '" + file_name + "' for line " +
-            std::to_string(line_number) + ": '" + line + "'; value number " + std::to_string(values.size()+1) + " is not a double.");
-      }
-      values.push_back(value);
-    }
-
-    if (values.size() != num_values_per_line) { 
-      throw std::runtime_error("Error reading values for the temporal values file '" + file_name + "' for line " +
-          std::to_string(line_number) + ": '" + line + "'; expected " + std::to_string(num_values_per_line) + " values per line.");
-    }
-
-    temporal_values.push_back(values);
-    line_number += 1;
-  }
-
-  if (lBc.gt.lrmp) {
-    lBc.gt.n = 1;
-  }
-
-  lBc.gt.qi.resize(lBc.gt.d); 
-  lBc.gt.qs.resize(lBc.gt.d);
-  lBc.gt.r.resize(lBc.gt.d,lBc.gt.n);
-  lBc.gt.i.resize(lBc.gt.d,lBc.gt.n);
-
-  fft(i, temporal_values, lBc.gt);
-}
-
-// [NOTE] This is a hack, should really just have a single function for this..
-//
-void read_temporal_values(const std::string& file_name, bfType& lBf) 
-{
-  std::ifstream temporal_values_file;
-  temporal_values_file.open(file_name);
-  if (!temporal_values_file.is_open()) {
-    throw std::runtime_error("Failed to open the temporal values file '" + file_name + "'.");
-  }
-
-  int i, j;
-  temporal_values_file >> i >> j; 
-  if (i < 2) {
-    throw std::runtime_error("The temporal values file '" + file_name + "' has an incorrect format.");
-  }
-  lBf.bt.n = j;
-
-  // Read time/value pairs.
-  //
-  std::vector<std::vector<double>> temporal_values;
-  double time, value;
-  std::string line;
-
-  while (std::getline(temporal_values_file, line)) { 
-    line.erase(std::remove(line.begin(), line.end(), '\r'), line.end());
-    if (line == "") {
-      continue;
-    }
-    std::istringstream line_input(line);
-    std::vector<double> values;
-    while (line_input >> value) {
-      values.push_back(value);
-    }
-    temporal_values.push_back(values);
-  }
-
-  if (lBf.bt.lrmp) {
-    lBf.bt.n = 1;
-  }
-
-  lBf.bt.qi.resize(lBf.bt.d); 
-  lBf.bt.qs.resize(lBf.bt.d);
-  lBf.bt.r.resize(lBf.bt.d, lBf.bt.n);
-  lBf.bt.i.resize(lBf.bt.d, lBf.bt.n);
-
-  fft(i, temporal_values, lBf.bt);
-}
-
 //----------------
 // read_trac_bcff
 //----------------
diff --git a/Code/Source/solver/read_files.h b/Code/Source/solver/read_files.h
index 3cccfc469..197f83fdb 100644
--- a/Code/Source/solver/read_files.h
+++ b/Code/Source/solver/read_files.h
@@ -35,9 +35,6 @@ namespace read_files_ns {
 
   void read_files(Simulation* simulation, const std::string& file_name);
 
-  void read_fourier_coeff_values_file(const std::string& file_name, bcType& lBc);
-  void read_fourier_coeff_values_file(const std::string& file_name, bfType& lBf);
-
   void read_ls(Simulation* simulation, EquationParameters* eq_params, consts::SolverType solver_type, eqType& lEq);
 
   void read_mat_model(Simulation* simulation, EquationParameters* eq_params, DomainParameters* domain_params, dmnType& lDmn);
@@ -48,9 +45,6 @@ namespace read_files_ns {
 
   void read_spatial_values(const ComMod& com_mod, const mshType& msh, const faceType& lFa, const std::string& file_name, bcType& lBc);
 
-  void read_temporal_values(const std::string& file_name, bcType& lBc);
-  void read_temporal_values(const std::string& file_name, bfType& lBf);
-
   void read_temp_spat_values(const ComMod& com_mod, const mshType& msh, const faceType& lFa, 
       const std::string& file_name, bcType& lBc);
   void read_temp_spat_values(const ComMod& com_mod, const mshType& msh, const std::string& file_name, bfType& lBf);
diff --git a/Code/Source/solver/set_bc.cpp b/Code/Source/solver/set_bc.cpp
index 01d361ceb..c9f92735f 100644
--- a/Code/Source/solver/set_bc.cpp
+++ b/Code/Source/solver/set_bc.cpp
@@ -1075,10 +1075,9 @@ void set_bc_dir_l(ComMod& com_mod, const bcType& lBc, const faceType& lFa, Array
     #ifdef debug_set_bc_dir_l
     dmsg << "bType_ustd";
     #endif
-    Vector<double> dirY_v(1), dirA_v(1);
-    ifft(com_mod, lBc.gt, dirY_v, dirA_v);
-    dirY = dirY_v(0);
-    dirA = dirA_v(0);
+    const auto value_and_derivative = lBc.gt.value_and_derivative(com_mod.time);
+    dirY = value_and_derivative.first[0];
+    dirA = value_and_derivative.second[0];
 
   } else { 
     dirA = 0.0;
@@ -1427,7 +1426,7 @@ void set_bc_neu_l(ComMod& com_mod, const CmMod& cm_mod, const bcType& lBc, const
   int nsd = com_mod.nsd;
 
   int nNo = lFa.nNo;
-  Vector<double> h(1), rtmp(1);
+  Vector<double> h(1);
   Vector<double> tmpA(nNo); 
 
   // Geting the contribution of Neu BC
@@ -1466,7 +1465,7 @@ void set_bc_neu_l(ComMod& com_mod, const CmMod& cm_mod, const bcType& lBc, const
        h(0) = lBc.g;
 
      } else if (utils::btest(lBc.bType,iBC_ustd)) {
-       ifft(com_mod, lBc.gt, h, rtmp);
+       h = lBc.gt.value(com_mod.time);
 
      } else {
        throw std::runtime_error("[set_bc_neu_l] Correction in SETBCNEU is needed");
@@ -1827,12 +1826,11 @@ void set_bc_trac_l(ComMod& com_mod, const CmMod& cm_mod, const bcType& lBc, cons
     }
 
   } else if (utils::btest(lBc.bType,iBC_ustd)) {
-     if (lBc.gt.d != nsd) {
-       throw std::runtime_error("[set_bc_trac_l]  Traction dof not initialized properly");
-     }
-     Vector<double> h(nsd);
-     Vector<double> tmp(nsd);
-     ifft(com_mod, lBc.gt, h, tmp);
+    if (lBc.gt.get_n_components() != nsd) {
+      throw std::runtime_error(
+          "[set_bc_trac_l]  Traction dof not initialized properly");
+    }
+     const Vector<double> h = lBc.gt.value(com_mod.time);
      for (int a = 0; a < nNo; a++) {
        int Ac = lFa.gN(a);
        hg.set_col(Ac, h);
diff --git a/tests/unitTests/math_tests/test_fft.cpp b/tests/unitTests/math_tests/test_fft.cpp
deleted file mode 100644
index 044472f29..000000000
--- a/tests/unitTests/math_tests/test_fft.cpp
+++ /dev/null
@@ -1,93 +0,0 @@
-/* Copyright (c) Stanford University, The Regents of the University of California, and others.
- *
- * All Rights Reserved.
- *
- * See Copyright-SimVascular.txt for additional details.
- *
- * Permission is hereby granted, free of charge, to any person obtaining
- * a copy of this software and associated documentation files (the
- * "Software"), to deal in the Software without restriction, including
- * without limitation the rights to use, copy, modify, merge, publish,
- * distribute, sublicense, and/or sell copies of the Software, and to
- * permit persons to whom the Software is furnished to do so, subject
- * to the following conditions:
- *
- * The above copyright notice and this permission notice shall be included
- * in all copies or substantial portions of the Software.
- *
- * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS
- * IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
- * TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
- * PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER
- * OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
- * EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
- * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
- * PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
- * LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
- * NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
- * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
- */
-
-#include "fft.h"
-#include "ComMod.h"
-#include "../test_common.h"
-#include <cmath>
-
-class FFTTest : public ::testing::Test {
-protected:
-    void SetUp() override {}
-    
-    void TearDown() override {}
-
-    // Generate temporal values for function f(t) = sin(t) + cos(t) + 0.1*t
-    void CreateTemporalValues(int N, double x_start, double x_end, 
-                             std::vector<std::vector<double>>& temporal_values) {
-        temporal_values.clear();
-        temporal_values.reserve(N);
-        
-        double step = (x_end - x_start) / (N - 1);
-        for (int i = 0; i < N; ++i) {
-            double t = x_start + i * step;
-            double y = std::sin(t) + std::cos(t) + 0.1 * t;
-            temporal_values.push_back({t, y});
-        }
-    }
-    
-    void InitializeFourierCoefficients(fcType& gt, int d = 1, int n = 16) {
-        gt.d = d;
-        gt.n = n;
-        gt.qi.resize(gt.d);
-        gt.qs.resize(gt.d);
-        gt.r.resize(gt.d, gt.n);
-        gt.i.resize(gt.d, gt.n);
-    }
-};
-
-TEST_F(FFTTest, SinCosLinearCombination) {
-    // Creates a temporal values function: f(t) = sin(t) + cos(t) + 0.1*t 
-    // Test finds the interpolated fourier coefficients of this function using fft.cpp
-    int N = 100;            // 100 timesteps
-    double x_start = 0.0;   // start time
-    double x_end = 10.0;    // end time 
-    std::vector<std::vector<double>> temporal_values;
-    
-    CreateTemporalValues(N, x_start, x_end, temporal_values);
-    
-    fcType gt;
-    // Apply interpolation to 1 dimensional data with n = 16 Fourier modes
-    InitializeFourierCoefficients(gt, 1, 16);
-    
-    // Compute the Fourier coefficients
-    fft(N, temporal_values, gt);
-
-    // Check the slope (first Fourier coefficient)
-    ASSERT_NEAR(gt.qs[0], -0.13830, 1e-2) << "Expected slope ~-0.13830";
-
-    // Check the real and imaginary components of the first three Fourier coefficients
-    ASSERT_NEAR(gt.r(0, 0), 0.32094, 1e-2) << "Expected first real coefficient to be close to 0.32094";
-    ASSERT_NEAR(gt.i(0, 0), 0.0, 1e-2) << "Expected first imaginary coefficient to be close to 0.0";
-    ASSERT_NEAR(gt.r(0, 1), 0.42759, 1e-2) << "Expected second real coefficient to be close to 0.42759";
-    ASSERT_NEAR(gt.i(0, 1), 1.25295, 1e-2) << "Expected second imaginary coefficient to be close to 1.25295";
-    ASSERT_NEAR(gt.r(0, 2), -0.44685, 1e-2) << "Expected third real coefficient to be close to -0.44685";
-    ASSERT_NEAR(gt.i(0, 2), -0.65403, 1e-2) << "Expected third imaginary coefficient to be close to -0.65403";
-}
\ No newline at end of file
diff --git a/tests/unitTests/math_tests/test_fft.h b/tests/unitTests/math_tests/test_fft.h
deleted file mode 100644
index 33c1c2611..000000000
--- a/tests/unitTests/math_tests/test_fft.h
+++ /dev/null
@@ -1,57 +0,0 @@
-/* Copyright (c) Stanford University, The Regents of the University of California, and others.
- *
- * All Rights Reserved.
- *
- * See Copyright-SimVascular.txt for additional details.
- *
- * Permission is hereby granted, free of charge, to any person obtaining
- * a copy of this software and associated documentation files (the
- * "Software"), to deal in the Software without restriction, including
- * without limitation the rights to use, copy, modify, merge, publish,
- * distribute, sublicense, and/or sell copies of the Software, and to
- * permit persons to whom the Software is furnished to do so, subject
- * to the following conditions:
- *
- * The above copyright notice and this permission notice shall be included
- * in all copies or substantial portions of the Software.
- *
- * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS
- * IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
- * TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
- * PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER
- * OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
- * EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
- * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
- * PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
- * LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
- * NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
- * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
- */
-
-#ifndef FFT_TEST_H
-#define FFT_TEST_H
-
-#include "../test_common.h"
-#include <vector>
-#include <cmath>
-
-// Forward declarations for types used in FFT testing
-struct fcType;
-
-// Test fixture class for FFT functionality
-class FFTTest : public ::testing::Test {
-protected:
-    // Set up test environment before each test
-    void SetUp() override;
-    
-    // Clean up test environment after each test  
-    void TearDown() override;
-    
-    // Helper methods for common test operations
-    void CreateTemporalValues(int N, double x_start, double x_end, 
-                             std::vector<std::vector<double>>& temporal_values);
-
-    void InitializeFourierCoefficients(fcType& gt, int d = 1, int n = 16);
-};
-
-#endif // FFT_TEST_H
\ No newline at end of file
diff --git a/tests/unitTests/math_tests/test_fourier_interpolation.cpp b/tests/unitTests/math_tests/test_fourier_interpolation.cpp
new file mode 100644
index 000000000..74eb9303a
--- /dev/null
+++ b/tests/unitTests/math_tests/test_fourier_interpolation.cpp
@@ -0,0 +1,297 @@
+/* Copyright (c) Stanford University, The Regents of the University of
+ * California, and others.
+ *
+ * All Rights Reserved.
+ *
+ * See Copyright-SimVascular.txt for additional details.
+ *
+ * Permission is hereby granted, free of charge, to any person obtaining
+ * a copy of this software and associated documentation files (the
+ * "Software"), to deal in the Software without restriction, including
+ * without limitation the rights to use, copy, modify, merge, publish,
+ * distribute, sublicense, and/or sell copies of the Software, and to
+ * permit persons to whom the Software is furnished to do so, subject
+ * to the following conditions:
+ *
+ * The above copyright notice and this permission notice shall be included
+ * in all copies or substantial portions of the Software.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS
+ * IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
+ * TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
+ * PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER
+ * OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
+ * EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
+ * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
+ * PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
+ * LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
+ * NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+ * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ */
+
+#include "../test_common.h"
+#include "fourier_interpolation.h"
+#include <cmath>
+
+class FourierInterpolationTest : public ::testing::Test {
+protected:
+  void SetUp() override {}
+
+  void TearDown() override {}
+
+  // Generate temporal values for function f(t) = sin(t) + cos(t) + 0.1*t
+  void CreateTemporalValues(int N, double x_start, double x_end,
+                            std::vector<std::vector<double>> &temporal_values) {
+    temporal_values.clear();
+    temporal_values.reserve(N);
+
+    const double step = (x_end - x_start) / (N - 1);
+    for (int i = 0; i < N; ++i) {
+      const double t = x_start + i * step;
+      const double y = interpolated(t, 0);
+      const double z = interpolated(t, 1);
+
+      temporal_values.push_back({t, y, z});
+    }
+  }
+
+  // Interpolated function.
+  static double interpolated(const double t, const unsigned int component) {
+    if (component == 0) {
+      return std::sin(t) + std::cos(t) + 0.1 * t;
+    } else {
+      return 2.0 * std::sin(t) + 3.0 * std::cos(t) + 0.2 * t;
+    }
+  }
+
+  // Derivative of the interpolated function.
+  static double interpolated_derivative(const double t,
+                                        const unsigned int component) {
+    if (component == 0) {
+      return std::cos(t) - std::sin(t) + 0.1;
+    } else {
+      return 2.0 * std::cos(t) - 3.0 * std::sin(t) + 0.2;
+    }
+  }
+};
+
+TEST_F(FourierInterpolationTest, FromTimeSeries) {
+  // Construct a FourierInterpolation instance from a time series and check that
+  // the dimensions of the resulting object are correct.
+
+  int N = 100;                   // 100 timesteps
+  double x_start = 0.0;          // start time
+  double x_end = 2 * consts::pi; // end time
+  std::vector<std::vector<double>> temporal_values;
+
+  CreateTemporalValues(N, x_start, x_end, temporal_values);
+
+  // Compute the Fourier interpolation.
+  const FourierInterpolation gt = FourierInterpolation::from_time_series(
+      /* n_fourier_coefficients = */ 3, temporal_values,
+      /* use_ramp = */ false);
+
+  // Check the number of components and Fourier coefficients.
+  ASSERT_EQ(gt.get_n_components(), 2);
+  ASSERT_EQ(gt.get_n_fourier_coefficients(), 3);
+}
+
+TEST_F(FourierInterpolationTest, FromFourierCoefficients) {
+  // Construct a FourierInterpolation instance from Fourier coefficients and
+  // check that the dimensions of the resulting object are correct, that the
+  // coefficients are set correctly, and that the object evaluates to the
+  // correct function.
+
+  const Vector<double> q_i = {0.0, 0.0};
+  const Vector<double> q_s = {0.1, 0.2};
+  const Array<double> c_real = {{0.0, 1.0}, {0.0, 3.0}};
+  const Array<double> c_imag = {{0.0, -1.0}, {0.0, -2.0}};
+  const double initial_time = 0.0;
+  const double period = 2 * consts::pi;
+
+  const FourierInterpolation gt =
+      FourierInterpolation::from_fourier_coefficients(q_i, q_s, c_real, c_imag,
+                                                      initial_time, period);
+
+  // Check the number of components and Fourier coefficients.
+  ASSERT_EQ(gt.get_n_components(), 2);
+  ASSERT_EQ(gt.get_n_fourier_coefficients(), 2);
+
+  // Check the linear trend coefficients.
+  ASSERT_EQ(gt.get_linear_trend_initial_value(0), 0.0);
+  ASSERT_EQ(gt.get_linear_trend_initial_value(1), 0.0);
+  ASSERT_EQ(gt.get_linear_trend_slope(0), 0.1);
+  ASSERT_EQ(gt.get_linear_trend_slope(1), 0.2);
+
+  // Check the Fourier coefficients.
+  ASSERT_EQ(gt.get_coefficient_real(0, 0), 0.0);
+  ASSERT_EQ(gt.get_coefficient_real(0, 1), 1.0);
+  ASSERT_EQ(gt.get_coefficient_real(1, 0), 0.0);
+  ASSERT_EQ(gt.get_coefficient_real(1, 1), 3.0);
+  ASSERT_EQ(gt.get_coefficient_imaginary(0, 0), 0.0);
+  ASSERT_EQ(gt.get_coefficient_imaginary(0, 1), -1.0);
+  ASSERT_EQ(gt.get_coefficient_imaginary(1, 0), 0.0);
+  ASSERT_EQ(gt.get_coefficient_imaginary(1, 1), -2.0);
+
+  // Check that the interpolation evaluates to the correct function.
+  Vector<double> value, derivative;
+  for (double t = initial_time; t <= initial_time + period; t += 0.5) {
+    std::tie(value, derivative) = gt.value_and_derivative(t);
+
+    const double expected_y = interpolated(t, 0);
+    const double expected_dy = interpolated_derivative(t, 0);
+    const double expected_z = interpolated(t, 1);
+    const double expected_dz = interpolated_derivative(t, 1);
+
+    ASSERT_NEAR(value[0], expected_y, 1e-6);
+    ASSERT_NEAR(value[1], expected_z, 1e-6);
+    ASSERT_NEAR(derivative[0], expected_dy, 1e-6);
+    ASSERT_NEAR(derivative[1], expected_dz, 1e-6);
+  }
+}
+
+TEST_F(FourierInterpolationTest, ValueAndDerivative) {
+  // Evaluate the interpolation and its derivative at a few points and compare
+  // with the original function. Since the interpolated function is a sin-cos
+  // combination plus a linear term, the interpolation is expected to match it
+  // exactly (as long as we use at least two Fourier coefficients).
+
+  int N = 100;                   // 100 timesteps
+  double x_start = 0.0;          // start time
+  double x_end = 2 * consts::pi; // end time
+  std::vector<std::vector<double>> temporal_values;
+
+  CreateTemporalValues(N, x_start, x_end, temporal_values);
+
+  // Compute the Fourier interpolation.
+  const FourierInterpolation gt = FourierInterpolation::from_time_series(
+      /* n_fourier_coefficients = */ 3, temporal_values,
+      /* use_ramp = */ false);
+
+  Vector<double> value, derivative;
+  for (double t = x_start; t <= x_end; t += 0.5) {
+    const double expected_y = interpolated(t, 0);
+    const double expected_dy = interpolated_derivative(t, 0);
+
+    const double expected_z = interpolated(t, 1);
+    const double expected_dz = interpolated_derivative(t, 1);
+
+    std::tie(value, derivative) = gt.value_and_derivative(t);
+
+    ASSERT_NEAR(value[0], expected_y, 1e-2) << "Value mismatch at t = " << t;
+    ASSERT_NEAR(derivative[0], expected_dy, 1e-2)
+        << "Derivative mismatch at t = " << t;
+    ASSERT_NEAR(value[1], expected_z, 1e-2) << "Value mismatch at t = " << t;
+    ASSERT_NEAR(derivative[1], expected_dz, 1e-2)
+        << "Derivative mismatch at t = " << t;
+  }
+
+  // Check periodicity.
+  Vector<double> value_2, derivative_2;
+  std::tie(value, derivative) = gt.value_and_derivative(x_start + 0.1);
+  std::tie(value_2, derivative_2) = gt.value_and_derivative(x_end + 0.1);
+  ASSERT_NEAR(value[0], value_2[0], 1e-2) << "Periodic value mismatch";
+  ASSERT_NEAR(value[1], value_2[1], 1e-2) << "Periodic value mismatch";
+  ASSERT_NEAR(derivative[0], derivative_2[0], 1e-2)
+      << "Periodic derivative mismatch";
+  ASSERT_NEAR(derivative[1], derivative_2[1], 1e-2)
+      << "Periodic derivative mismatch";
+}
+
+TEST_F(FourierInterpolationTest, Coefficients) {
+  // Construct the interpolation and compare the computed Fourier and linear
+  // trend coefficients with the expected values. Notice that we're using a
+  // period here that is different from that of the interpolated function (i.e.
+  // not 2pi). Therefore, the Fourier coefficients will not be the same as those
+  // of the original function.
+
+  int N = 100;          // 100 timesteps
+  double x_start = 0.0; // start time
+  double x_end = 10.0;  // end time
+  std::vector<std::vector<double>> temporal_values;
+
+  CreateTemporalValues(N, x_start, x_end, temporal_values);
+
+  // Compute the Fourier coefficients
+  const FourierInterpolation gt = FourierInterpolation::from_time_series(
+      /* n_fourier_coefficients = */ 16, temporal_values,
+      /* use_ramp = */ false);
+
+  // Check the linear trend slope
+  ASSERT_NEAR(gt.get_linear_trend_slope(0), -0.13830, 1e-2)
+      << "Expected slope ~-0.13830";
+
+  // Check the real and imaginary components of the first three Fourier
+  // coefficients
+  ASSERT_NEAR(gt.get_coefficient_real(0, 0), 0.32094, 1e-2)
+      << "Expected first real coefficient to be close to 0.32094";
+  ASSERT_NEAR(gt.get_coefficient_imaginary(0, 0), 0.0, 1e-2)
+      << "Expected first imaginary coefficient to be close to 0.0";
+  ASSERT_NEAR(gt.get_coefficient_real(0, 1), 0.42759, 1e-2)
+      << "Expected second real coefficient to be close to 0.42759";
+  ASSERT_NEAR(gt.get_coefficient_imaginary(0, 1), 1.25295, 1e-2)
+      << "Expected second imaginary coefficient to be close to 1.25295";
+  ASSERT_NEAR(gt.get_coefficient_real(0, 2), -0.44685, 1e-2)
+      << "Expected third real coefficient to be close to -0.44685";
+  ASSERT_NEAR(gt.get_coefficient_imaginary(0, 2), -0.65403, 1e-2)
+      << "Expected third imaginary coefficient to be close to -0.65403";
+}
+
+TEST_F(FourierInterpolationTest, Ramp) {
+  // Construct a FourierInterpolation instance from a time series and check that
+  // the ramp option behaves correctly.
+
+  int N = 100;          // 100 timesteps
+  double x_start = 0.0; // start time
+  double x_end = 10.0;  // end time
+  std::vector<std::vector<double>> temporal_values;
+
+  CreateTemporalValues(N, x_start, x_end, temporal_values);
+
+  // Compute the Fourier coefficients
+  const FourierInterpolation gt = FourierInterpolation::from_time_series(
+      /* n_fourier_coefficients = */ 16, temporal_values,
+      /* use_ramp = */ true);
+
+  // Check that the number of Fourier coefficients is 1 (since ramp is used)
+  ASSERT_EQ(gt.get_n_fourier_coefficients(), 1)
+      << "Expected number of Fourier coefficients to be 1 when using ramp";
+
+  Vector<double> value, derivative;
+
+  // Check that before the initial time the value is equal to the initial value
+  // of the linear trend.
+  std::tie(value, derivative) = gt.value_and_derivative(x_start - 1.0);
+  ASSERT_NEAR(value[0], temporal_values[0][1], 1e-6)
+      << "Expected value before initial time to be equal to initial value of "
+         "linear trend";
+
+  // Check that after the final time the value is equal to the final value of
+  // the linear trend.
+  std::tie(value, derivative) = gt.value_and_derivative(x_end + 1.0);
+  ASSERT_NEAR(value[0], temporal_values.back()[1], 1e-6)
+      << "Expected value after final time to be equal to final value of linear "
+         "trend";
+
+  // Check that the value in the interpolation interval is a linear
+  // interpolation of the initial and final value.
+  for (double t = x_start; t <= x_end; t += 0.5) {
+    std::tie(value, derivative) = gt.value_and_derivative(t);
+
+    const double expected_y =
+        temporal_values[0][1] +
+        (temporal_values.back()[1] - temporal_values[0][1]) * (t - x_start) /
+            (x_end - x_start);
+    const double expected_z =
+        temporal_values[0][2] +
+        (temporal_values.back()[2] - temporal_values[0][2]) * (t - x_start) /
+            (x_end - x_start);
+
+    ASSERT_NEAR(value[0], expected_y, 1e-6)
+        << "Expected value in interpolation interval to be linear "
+           "interpolation of initial and final value";
+    ASSERT_NEAR(value[1], expected_z, 1e-6)
+        << "Expected value in interpolation interval to be linear "
+           "interpolation of initial and final value";
+  }
+}
\ No newline at end of file