diff --git a/R/scova.R b/R/scova.R
index e3fb5d5..8e89b52 100644
--- a/R/scova.R
+++ b/R/scova.R
@@ -152,7 +152,7 @@ scova <- R6::R6Class(
       dt_out <- merge(
         dt_samples_wide, dt_times, by = "t_id", allow.cartesian = TRUE)
 
-      dt_out[, mu := private$simulate_trajectory(
+      dt_out[, mu := scova_simulate_trajectory(
         t, t0_pop, tp_pop, ts_pop, m1_pop, m2_pop, m3_pop),
                by = c("t", "p", "k", ".draw")]
 
@@ -166,20 +166,6 @@ scova <- R6::R6Class(
 
       data.table::setcolorder(dt_out, c("t", "p", "k"))
     },
-    simulate_trajectory = function(t, t0, tp, ts, m1, m2, m3) {
-      mu <- t0
-      if (t < tp) {
-        mu <- mu + m1 * t
-      } else if (t <= ts) {
-        mu <- mu + m1 * tp + m2 * (t - tp)
-      } else if (t > ts) {
-        mu <- mu + m1 * tp + m2 * (ts - tp) + m3 * (t - ts)
-      }
-      max(mu, 0)
-    },
-    simulate_trajectories = function(dat) {
-      data.table::setDT(simulate_trajectories_cpp(dat))
-    },
     prepare_stan_data = function() {
       stan_id <- titre <- censored <- titre_type_num <- titre_type <- obs_id <- t_since_last_exp <- t_since_min_date <- NULL
       stan_data <- list(
@@ -258,7 +244,7 @@ scova <- R6::R6Class(
 
       dt_out
     },
-    adjust_parameters <- function(dt) {
+    adjust_parameters = function(dt) {
       params_to_adjust <- c(
         "t0_pop", "tp_pop", "ts_pop", "m1_pop", "m2_pop", "m3_pop")
       # Loop through the parameters you want to adjust
@@ -396,11 +382,11 @@ scova <- R6::R6Class(
 
       logger::log_info("Calculating peak and switch titre values")
       dt_peak_switch[, `:=`(
-        mu_0 = private$simulate_trajectory(
+        mu_0 = scova_simulate_trajectory(
           0, t0_pop, tp_pop, ts_pop, m1_pop, m2_pop, m3_pop),
-        mu_p = private$simulate_trajectory(
+        mu_p = scova_simulate_trajectory(
           tp_pop, t0_pop, tp_pop, ts_pop, m1_pop, m2_pop, m3_pop),
-        mu_s = private$simulate_trajectory(
+        mu_s = scova_simulate_trajectory(
           ts_pop, t0_pop, tp_pop, ts_pop, m1_pop, m2_pop, m3_pop)),
                        by = c("p", "k", ".draw")]
 
@@ -465,7 +451,7 @@ scova <- R6::R6Class(
       # Running the C++ code to simulate trajectories for each parameter sample
       # for each individual
       logger::log_info("Simulating individual trajectories")
-      dt_params_ind_traj <- private$simulate_trajectories(dt_params_ind_trim)
+      dt_params_ind_traj <- scova_simulate_trajectories(dt_params_ind_trim)
 
       dt_params_ind_traj <- data.table::setDT(convert_log_scale_inverse_cpp(
           dt_params_ind_traj, vars_to_transform = "mu"))
diff --git a/R/simulate.R b/R/simulate.R
new file mode 100644
index 0000000..65e3987
--- /dev/null
+++ b/R/simulate.R
@@ -0,0 +1,15 @@
+scova_simulate_trajectory <- function(t, t0, tp, ts, m1, m2, m3) {
+  mu <- t0
+  if (t < tp) {
+    mu <- mu + m1 * t
+  } else if (t <= ts) {
+    mu <- mu + m1 * tp + m2 * (t - tp)
+  } else if (t > ts) {
+    mu <- mu + m1 * tp + m2 * (ts - tp) + m3 * (t - ts)
+  }
+  max(mu, 0)
+}
+
+scova_simulate_trajectories <- function(dat) {
+  data.table::setDT(simulate_trajectories_cpp(dat))
+}
diff --git a/tests/testthat/test-simulate-individual-trajectories.R b/tests/testthat/test-simulate-individual-trajectories.R
new file mode 100644
index 0000000..4c1037a
--- /dev/null
+++ b/tests/testthat/test-simulate-individual-trajectories.R
@@ -0,0 +1,26 @@
+test_that("Cpp and R function produce same values", {
+  t_max <- 6L
+  dat_pop <- data.table(t = 0:t_max,
+                        t0 = 0,
+                        tp = 3,
+                        ts = 5,
+                        m1 = 0.5,
+                        m2 = 1,
+                        m3 = 1.5
+  )
+  res_pop <- dat_pop[, mu := scova_simulate_trajectory(t, t0, tp, ts, m1, m2, m3), by = "t"]
+  dat_ind <- data.table(stan_id = 1L,
+                        t_max = t_max,
+                        titre_type = 3L,
+                        draw = 10L,
+                        t0_ind = 0,
+                        tp_ind = 3,
+                        ts_ind = 5,
+                        m1_ind = 0.5,
+                        m2_ind = 1,
+                        m3_ind = 1.5
+  )
+
+  res_ind <- scova_simulate_trajectories(dat_ind)
+  expect_equal(res_pop$mu, res_ind$mu)
+})