diff --git a/CHANGELOG.md b/CHANGELOG.md
index 155eab9..cf37de5 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -33,6 +33,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 - docstrings and documentation files for refactored functions [#160](https://github.com/RECETOX/recetox-aplcms/pull/160)
 - refactored parameter names to keep them more harmonized [#167](https://github.com/RECETOX/recetox-aplcms/pull/167)
 - moved some utility functions to a more suitable locations [#164](https://github.com/RECETOX/recetox-aplcms/pull/164)
+- refactored and documented `prof.to.features.R` [#170](https://github.com/RECETOX/recetox-aplcms/pull/170)
 
 ### Removed
 - `extract_features` and `feature.align` [#154](https://github.com/RECETOX/recetox-aplcms/pull/154)
diff --git a/NAMESPACE b/NAMESPACE
index c46aa30..6387392 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -75,7 +75,7 @@ export(semi.sup)
 export(span)
 export(two.step.hybrid)
 export(unsupervised)
-export(validate_inputs)
+export(validate_model_method_input)
 import(MASS)
 import(arrow)
 import(doParallel)
diff --git a/R/prof.to.features.R b/R/prof.to.features.R
index 383d8a9..a573dd8 100644
--- a/R/prof.to.features.R
+++ b/R/prof.to.features.R
@@ -1,10 +1,10 @@
-
-#' Validate that provided inputs match expected, exit execution otherwise.
+#' @description
+#' Validate that provided model and method inputs match expected, exit execution otherwise.
 #' @param shape_model The mathematical model for the shape of a peak. There are two choices - "bi-Gaussian" and "Gaussian".
 #'  When the peaks are asymmetric, the bi-Gaussian is better. The default is "bi-Gaussian".
 #' @param peak_estim_method The estimation method for the bi-Gaussian peak model. Two possible values: moment and EM.
 #' @export
-validate_inputs <- function(shape_model, peak_estim_method) {
+validate_model_method_input <- function(shape_model, peak_estim_method) {
   if (!shape_model %in% c("Gaussian", "bi-Gaussian")) {
     stop("shape_model argument must be 'Gaussian' or 'bi-Gaussian'")
   }
@@ -13,6 +13,8 @@ validate_inputs <- function(shape_model, peak_estim_method) {
   }
 }
 
+
+#' @description 
 #' Initialize minimum and maximum bandwidth values if none given. Ensure that minimum bandwidth is lower that maximum, else set minimum to 1/4 of maximum value.
 #' @param min_bandwidth The minimum bandwidth to use in the kernel smoother.
 #' @param max_bandwidth The maximum bandwidth to use in the kernel smoother.
@@ -42,7 +44,8 @@ preprocess_bandwidth <- function(min_bandwidth, max_bandwidth, profile) {
   return(list("min_bandwidth" = min_bandwidth, "max_bandwidth" = max_bandwidth))
 }
 
-#' Convert input matrix to a dataframe with column names (see source code for the names).
+#' @description
+#' Convert matrix to dataframe and rename columns.
 #' @param profile Profile table with shape number-of-features*4. The table contains following columns:
 #' \itemize{
 #'   \item float - mass-to-charge ratio of feature
@@ -61,13 +64,12 @@ preprocess_bandwidth <- function(min_bandwidth, max_bandwidth, profile) {
 preprocess_profile <- function(profile) {
   keys <- c("mz", "rt", "intensity", "group_number")
   colnames(profile) <- keys
-
   return(data.frame(profile))
 }
 
+#' @description
 #' Compute parameters of chromatographic peak shape if peaks are considered to be gaussian
 #' @param rt_profile A matrix with two columns: "base.curve" (rt) and "intensity".
-#' @param moment_power The parameter for data transformation when fitting the bi-Gaussian or Gaussian mixture model in an EIC.
 #' @param bw Bandwidth vector to use in the kernel smoother.
 #' @param component_eliminate When a component accounts for a proportion of intensities less than this value, the component will be ignored.
 #' @param BIC_factor The factor that is multiplied on the number of parameters to modify the BIC criterion. If larger than 1,
@@ -81,8 +83,8 @@ preprocess_profile <- function(profile) {
 #'   \item scale - float - estimated total signal strength (total area of the estimated normal curve)
 #'}
 #' @export
-compute_gaussian_peak_shape <- function(rt_profile, moment_power, bw, component_eliminate, BIC_factor, aver_diff) {
-  rt_peak_shape <- normix.bic(rt_profile[, "base.curve"], rt_profile[, 2], moment_power = moment_power, bw = bw, eliminate = component_eliminate, BIC_factor = BIC_factor, aver_diff = aver_diff)$param
+compute_gaussian_peak_shape <- function(rt_profile, bw, component_eliminate, BIC_factor, aver_diff) {
+  rt_peak_shape <- normix.bic(rt_profile[, "base.curve"], rt_profile[, 2], bw = bw, eliminate = component_eliminate, BIC_factor = BIC_factor, aver_diff = aver_diff)$param
   if (nrow(rt_peak_shape) == 1) {
     rt_peak_shape <- c(rt_peak_shape[1, 1:2], rt_peak_shape[1, 2], rt_peak_shape[1, 3])
   } else {
@@ -91,7 +93,14 @@ compute_gaussian_peak_shape <- function(rt_profile, moment_power, bw, component_
   return(rt_peak_shape)
 }
 
-#' This function solves the value of a using the x, t, a from the previous step, and sigma.1, and sigma.2 (original authors' comment).
+#' @description
+#' This function solves peak summit using the x, t, peak summit from the previous step, 
+#' and sigma.1, and sigma.2 (original authors' comment).
+#' @param x A vector of numerical values (intensities).
+#' @param t A vector of numerical values (rt).
+#' @param a A vector of peak summits.
+#' @param sigma.1 left standard deviation of the gaussian curve
+#' @param sigma.2 right standard deviation of the gaussian curve
 #' @export
 solve.a <- function(x, t, a, sigma.1, sigma.2) {
   # This function is a part of bigauss.esti.EM and is not covered by any of test-cases
@@ -99,7 +108,12 @@ solve.a <- function(x, t, a, sigma.1, sigma.2) {
   return(sum(t * w) / sum(w))
 }
 
-#' This function prepares the parameters required for latter compuation. u, v, and sum of x (original authors' comment).
+#' @description
+#' This function prepares the parameters required for latter computation.
+#' u, v, and sum of x (original authors' comment).
+#' @param x A vector of numerical values (intensities).
+#' @param t A vector of numerical values (rt).
+#' @param a A vector of peak summits.
 #' @export
 prep.uv <- function(x, t, a) {
   # This function is a part of bigauss.esti.EM and is not covered by any of test-cases
@@ -113,7 +127,16 @@ prep.uv <- function(x, t, a) {
   ))
 }
 
-#' This function takes the value intensity level x, retention time t and assumed breaking point a, calculates the square estimated of sigma.1 and sigma.2 (original authors' comment).
+#' @description
+#' Calculates the square estimated of left and right standard deviations.
+#' @param x A vector of numerical values (intensities).
+#' @param t A vector of numerical values (rt).
+#' @param a A vector of peak summits.
+#' @return A vector of:
+#' \itemize{
+#'   \item standard deviation at the left side of the gaussian curve
+#'   \item standard deviation at the right side of the gaussian curve
+#' }
 #' @export
 solve.sigma <- function(x, t, a) {
   # This function is a part of bigauss.esti.EM and is not covered by any of test-cases
@@ -127,10 +150,22 @@ solve.sigma <- function(x, t, a) {
 }
 
 #' @description
-#' Function takes into x and t, and then computes the value of sigma.1, sigma.2 and a using iterative method. the returned values include estimated sigmas,
-#' a and a boolean variable on whether the termination criteria is satified upon the end of the program (original authors' comment).
+#' Computes bi-gaussian parameters using an iterative method. 
+#' @param t A vector of numerical values (rt).
+#' @param x A vector of numerical values (intensities).
+#' @param max.iter Maximum number of iterations.
+#' @param epsilon Threshold for continuing the iteration
+#' @param sigma_ratio_lim A vector of two. It enforces the belief of the range of the ratio between the left-standard deviation 
+#' and the right-standard deviation of the bi-Gaussian function.
+#' @return A vector with length 4. The items are as follows going from first to last:
+#' \itemize{
+#'   \item mean of gaussian curve
+#'   \item standard deviation at the left side of the gaussian curve
+#'   \item standard deviation at the right side of the gaussian curve
+#'   \item estimated total signal strength (total area of the estimated normal curve)
+#' }
 #' @export
-bigauss.esti.EM <- function(t, x, max.iter = 50, epsilon = 0.005, moment_power = 1, do.plot = FALSE, truth = NA, sigma_ratio_lim = c(0.3, 1)) {
+bigauss.esti.EM <- function(t, x, max.iter = 50, epsilon = 0.005, do.plot = FALSE, sigma_ratio_lim = c(0.3, 1)) {
   # This function is not covered by any test case
   sel <- which(x > 1e-10)
   if (length(sel) == 0) {
@@ -146,7 +181,7 @@ bigauss.esti.EM <- function(t, x, max.iter = 50, epsilon = 0.005, moment_power =
   ## a smaller than epsilon will terminate the iteration.
   ## epsilon <- min(diff(sort(t)))/2
 
-  ## using the median value of t as the initial value of a.
+  ## using the median value of t as the initial value of a (peak summit)
   a.old <- t[which(x == max(x))[1]]
   a.new <- a.old
   change <- 10 * epsilon
@@ -174,6 +209,7 @@ bigauss.esti.EM <- function(t, x, max.iter = 50, epsilon = 0.005, moment_power =
   return(c(a.new, sigma$sigma.1, sigma$sigma.2, scale))
 }
 
+#' @description
 #' Computes vector of cumulative sums on reversed input. Returns cumulative sum vector going from the sum of all elements to one.
 #' @param x float - vector of numerical values
 #' @return Returns a vector
@@ -183,7 +219,11 @@ rev_cum_sum <- function(x) {
   return(rev(cumsum(x)))
 }
 
-#' TODO: Document
+#' @description
+#' Computes initial bound of set of values.
+#' @param x Cumulative intensity values.
+#' @param left_sigma_ratio_lim. Left-standard deviation of the bi-Gaussian function.
+#' @return Returns end bound.
 #' @export
 compute_start_bound <- function(x, left_sigma_ratio_lim) {
   start_bound <- 1
@@ -196,7 +236,11 @@ compute_start_bound <- function(x, left_sigma_ratio_lim) {
   return (start_bound)
 }
 
-#' TODO: Document
+#' @description
+#' Computes final bound of set of values.
+#' @param x Cumulative intensity values.
+#' @param rigth_sigma_ratio_lim. Right-standard deviation of the bi-Gaussian function.
+#' @return Returns end bound.
 #' @export
 compute_end_bound <- function(x, right_sigma_ratio_lim) {
   len_x <- length(x)
@@ -209,9 +253,11 @@ compute_end_bound <- function(x, right_sigma_ratio_lim) {
   return (end_bound)
 }
 
+#' @description
+#' Computes initial and final bounds of set of values.
 #' @param x Cumulative intensity values.
 #' @param sigma_ratio_lim A vector of two. It enforces the belief of the range of the ratio between the left-standard deviation.
-#'  and the right-standard deviation of the bi-Gaussian function used to fit the data.
+#' and the right-standard deviation of the bi-Gaussian function used to fit the data.
 #' @return Returns a list with bounds with following items:
 #' \itemize{
 #'   \item start - start bound
@@ -224,7 +270,9 @@ compute_bounds <- function(x, sigma_ratio_lim) {
   return(list(start = start, end = end))
 }
 
-#' Compute difference between neighbouring elements of a vector and optionally apply a mask such that the maximum difference is no higher than 4-fold minimum difference.
+#' @description
+#' Compute difference between neighbouring elements of a vector and optionally apply a 
+#' mask such that the maximum difference is no higher than 4-fold minimum difference.
 #' @param x - float - a vector of numerical values.
 #' @param apply_mask - boolean - whether to apply threshold mask to the output vector.
 #' @return Returns vector of numeric differences between neighbouring values.
@@ -248,7 +296,8 @@ compute_dx <- function(x, apply_mask=TRUE) {
   return (dx)
 }
 
-#' Find base.curve RTs that lay within RT range of the whole feature table and append the intensities to these RTs.
+#' @description
+#' Find base.curve RTs that lay within RT range of the whole feature table and append intensities to these RTs.
 #' @param profile Profile table with shape number-of-features*4. The table contains following columns:
 #' \itemize{
 #'   \item mz - float - mass-to-charge ratio of feature
@@ -261,17 +310,17 @@ compute_dx <- function(x, apply_mask=TRUE) {
 #' @export
 compute_chromatographic_profile <- function(profile, base.curve) {
   rt_range <- range(profile[, "rt"])
-  rt_profile <- base.curve[between(base.curve[, "base.curve"], min(rt_range), max(rt_range)), ]
+  rt_profile <- base.curve[dplyr::between(base.curve[, "base.curve"], min(rt_range), max(rt_range)), ]
   rt_profile[rt_profile[, "base.curve"] %in% profile[, "rt"], 2] <- profile[, "intensity"]
   colnames(rt_profile)[2] <- "intensity"
-
   return (rt_profile)
 }
 
-#' Estimate total signal strength (total area of the estimated normal curve).
+#' @description
+#' Estimates total signal strength (total area of the estimated normal curve).
 #' @param y - float - a vector of intensities.
 #' @param d - float - a vector of \emph{y} values in a gaussian curve.
-#' @param scale - float - a vector of scaled intensity values.
+#' @return scale - float - a vector of scaled intensity values.
 #' @export
 compute_scale <- function(y, d) {
   dy_ratio <- d^2 * log(y / d)
@@ -282,6 +331,7 @@ compute_scale <- function(y, d) {
   return (scale)
 }
 
+#' @description
 #' Estimate the parameters of Bi-Gaussian curve.
 #' @param x Vector of RTs that lay in the same RT cluster.
 #' @param y Intensities that belong to x.
@@ -362,7 +412,7 @@ bigauss.esti <- function(x, y, moment_power = 1, do.plot = FALSE, sigma_ratio_li
 
     d1 <- dnorm(x[sel1], sd = s1, mean = m)
     d2 <- dnorm(x[sel2], sd = s2, mean = m)
-    d <- c(d1 * s1, d2 * s2) # notice this "density" doesnt integrate to 1. Rather it integrates to (s1+s2)/2
+    d <- c(d1 * s1, d2 * s2) # notice this "density" does not integrate to 1. Rather it integrates to (s1+s2)/2
     y <- y.0
 
     scale <- compute_scale(y, d)
@@ -382,10 +432,17 @@ bigauss.esti <- function(x, y, moment_power = 1, do.plot = FALSE, sigma_ratio_li
   return(to.return)
 }
 
+#' @description
+#' Calculates the three initial bi-gaussian parameters (sd1, sd2, and scaling factor)
 #' @param rt_profile A matrix with two columns: "base.curve" (rt) and "intensity".
 #' @param vlys A vector of sorted RT-valley values at which the kernel estimate was computed.
 #' @param dx Difference between neighbouring RT values with step 2.
 #' @param pks A vector of sorted RT-peak values at which the kernel estimate was computed.
+#' @return A list. The items are as follows going from first to last:
+#' \itemize{
+#'   \item standard deviation at the left side of the gaussian curve
+#'   \item standard deviation at the right side of the gaussian curve
+#'   \item estimated total signal strength (total area of the estimated normal curve)
 #' @export
 compute_initiation_params <- function(rt_profile, vlys, dx, pks) {
   m <- s1 <- s2 <- delta <- pks
@@ -397,13 +454,16 @@ compute_initiation_params <- function(rt_profile, vlys, dx, pks) {
     sel.2 <- which(rt_profile[, "base.curve"] >= m[i] & rt_profile[, "base.curve"] < min(vlys[vlys > m[i]]))
     s2[i] <- sqrt(sum((rt_profile[sel.2, "base.curve"] - m[i])^2 * rt_profile[sel.2, "intensity"] * dx[sel.2]) / sum(rt_profile[sel.2, "intensity"] * dx[sel.2]))
 
-    delta[i] <- (sum(rt_profile[sel.1, "intensity"] * dx[sel.1]) + sum(rt_profile[sel.2, "intensity"] * dx[sel.2])) / ((sum(dnorm(rt_profile[sel.1, "base.curve"], mean = m[i], sd = s1[i])) * s1[i] / 2) + (sum(dnorm(rt_profile[sel.2, "base.curve"], mean = m[i], sd = s2[i])) * s2[i] / 2))
+    delta[i] <- (sum(rt_profile[sel.1, "intensity"] * dx[sel.1]) + sum(rt_profile[sel.2, "intensity"] * dx[sel.2])) / 
+    ((sum(dnorm(rt_profile[sel.1, "base.curve"], mean = m[i], sd = s1[i])) * s1[i] / 2) + (sum(dnorm(rt_profile[sel.2, "base.curve"], mean = m[i], sd = s2[i])) * s2[i] / 2))
   }
   return (list(s1 = s1,
     s2 = s2,
     delta = delta))
 }
 
+#' @description
+#' Computes the expectation step of the EM method.
 #' @param m A vector of sorted RT-peak values at which the kernel estimate was computed.
 #' @param rt_profile A matrix with two columns: "base.curve" (rt) and "intensity".
 #' @param delta Parameter computed by the initiation step.
@@ -429,6 +489,9 @@ compute_e_step <- function(m, rt_profile, delta, s1, s2) {
   return(list(fit = fit, sum.fit = sum.fit))
 }
 
+#' @description
+#' Estimates the optimal bi-gaussian parameters using the EM method. It accepts two internal computation of parameters for "moment"
+#' and "EM" model input options.
 #' @param rt_profile Dataframe that stores RTs and intensities of features.
 #' @param moment_power The parameter for data transformation when fitting the bi-Gaussian or Gaussian mixture model in an EIC.
 #' @param sigma_ratio_lim A vector of two. It enforces the belief of the range of the ratio between the left-standard deviation
@@ -455,9 +518,7 @@ bigauss.mix <- function(rt_profile, moment_power = 1, do.plot = FALSE, sigma_rat
   last.num.pks <- Inf
 
   rt_profile_unfiltered <- rt_profile
-  rt_profile <- data.frame(rt_profile) |>
-    filter(intensity > 1e-5) |>
-    arrange(base.curve)
+  rt_profile <- data.frame(rt_profile) |> dplyr::filter(intensity > 1e-5) |> dplyr::arrange(base.curve)
 
   for (bw.n in length(all.bw):1)
   {
@@ -527,7 +588,7 @@ bigauss.mix <- function(rt_profile, moment_power = 1, do.plot = FALSE, sigma_rat
           if (peak_estim_method == "moment") {
             this.fit <- bigauss.esti(rt_profile[, "base.curve"], this.y, moment_power = moment_power, do.plot = FALSE, sigma_ratio_lim = sigma_ratio_lim)
           } else {
-            this.fit <- bigauss.esti.EM(rt_profile[, "base.curve"], this.y, moment_power = moment_power, do.plot = FALSE, sigma_ratio_lim = sigma_ratio_lim)
+            this.fit <- bigauss.esti.EM(rt_profile[, "base.curve"], this.y, do.plot = FALSE, sigma_ratio_lim = sigma_ratio_lim)
           }
           m[i] <- this.fit[1]
           s1[i] <- this.fit[2]
@@ -583,7 +644,8 @@ bigauss.mix <- function(rt_profile, moment_power = 1, do.plot = FALSE, sigma_rat
   return(rec)
 }
 
-#' Reevaluate parameters of rtomatographic gaussian curves.
+#' @description
+#' Reevaluate parameters of chromatographic gaussian curves.
 #' @param that.curve Dataframe that stores RTs and intensities of features.
 #' @param pks A vector of sorted RT-peak values at which the kernel estimate was computed.
 #' @param vlys A vector of sorted RT-valley values at which the kernel estimate was computed.
@@ -592,6 +654,12 @@ bigauss.mix <- function(rt_profile, moment_power = 1, do.plot = FALSE, sigma_rat
 #' @param max.iter Maximum number of iterations when reevaluating gaussian curves.
 #' @param aver_diff Average retention time difference across RTs of all features.
 #' @importFrom dplyr between
+#' @return A list. The items are as follows going from first to last:
+#' \itemize{
+#'   \item miu - float - mean value of the gaussian curve
+#'   \item sigma - float - standard deviation of the gaussian curve
+#'   \item scale - float - estimated total signal strength (total area of the estimated normal curve)
+#'}
 #' @export
 normix <- function(that.curve, pks, vlys, ignore = 0.1, max.iter = 50, aver_diff) {
   x <- that.curve[, 1]
@@ -723,16 +791,17 @@ normix <- function(that.curve, pks, vlys, ignore = 0.1, max.iter = 50, aver_diff
   return(rec)
 }
 
+#' @description
+#' Estimates parameters of a gaussian curve.
 #' @param x Vector of RTs that lay in the same RT cluster.
 #' @param y Intensities that belong to x.
-#' @param moment_power The parameter for data transformation when fitting the bi-Gaussian or Gaussian mixture model in an EIC.
 #' @param bw Bandwidth vector to use in the kernel smoother.
 #' @param eliminate When a component accounts for a proportion of intensities less than this value, the component will be ignored.
 #' @param max.iter Maximum number of iterations when executing the E step.
 #' @param BIC_factor The factor that is multiplied on the number of parameters to modify the BIC criterion. If larger than 1,
 #' @param aver_diff Average retention time difference across RTs of all features.
 #' @export
-normix.bic <- function(x, y, moment_power = 2, do.plot = FALSE, bw = c(15, 30, 60), eliminate = .05, max.iter = 50, BIC_factor = 2, aver_diff) {
+normix.bic <- function(x, y, do.plot = FALSE, bw = c(15, 30, 60), eliminate = .05, max.iter = 50, BIC_factor = 2, aver_diff) {
   all.bw <- bw[order(bw)]
   sel <- y > 1e-5
   x <- x[sel]
@@ -799,12 +868,10 @@ normix.bic <- function(x, y, moment_power = 2, do.plot = FALSE, bw = c(15, 30, 6
 }
 
 #' Generate feature table from noise-removed LC/MS profile.
-#'
 #' @description
 #' Each LC/MS profile is first processed by the function proc.cdf() to remove noise and reduce data size. A matrix containing m/z
 #' value, retention time, intensity, and group number is output from proc.cdf(). This matrix is then fed to the function
 #' prof.to.features() to generate a feature list. Every detected feature is summarized into a single row in the output matrix from this function.
-#'
 #' @param profile The matrix output from proc.cdf(). It contains columns of m/z value, retention time, intensity and group number.
 #' @param bandwidth A value between zero and one. Multiplying this value to the length of the signal along the time axis helps
 #'  determine the bandwidth in the kernel smoother used for peak identification.
@@ -837,7 +904,7 @@ prof.to.features <- function(profile,
                              component_eliminate,
                              BIC_factor,
                              do.plot) {
-  validate_inputs(shape_model, peak_estim_method)
+  validate_model_method_input(shape_model, peak_estim_method)
 
   profile <- preprocess_profile(profile)
 
@@ -846,32 +913,40 @@ prof.to.features <- function(profile,
   max_bandwidth <- bws[["max_bandwidth"]]
 
   # base.curve <- compute_base_curve(profile[, "rt"])
-  base.curve <- sort(unique(profile[, "rt"]))
+  base.curve <- sort(unique(profile$rt))
   base.curve <- cbind(base.curve, base.curve * 0)
-  all_rts <- compute_delta_rt(base.curve[, 1])
-  aver_diff <- mean(diff(base.curve))
+  
+  all_diff_mean_rts <- compute_delta_rt(base.curve[, 1]) # computes diff of mean values from consecutive values 
+  aver_diff <- mean(diff(base.curve))  
 
   keys <- c("mz", "rt", "sd1", "sd2", "area")
   peak_parameters <- matrix(0, nrow = 0, ncol = length(keys), dimnames = list(NULL, keys))
-
+  
   feature_groups <- split(profile, profile$group_number)
+
+  # loop over each group
   for (i in seq_along(feature_groups))
   {
-    feature_group <- feature_groups[[i]]
-    feature_group <- feature_group[order(feature_group[, "rt"]), ]
+    # init variables
+    feature_group <- feature_groups[[i]] |> dplyr::arrange_at("rt") 
 
     num_features <- nrow(feature_group)
+    # The estimation procedure for a single peak
+    # Defines the dataframe containing median_mz, median_rt, sd1, sd2, and area
     if (dplyr::between(num_features, 2, 10)) {
-      eic_area <- interpol.area(feature_group[, "rt"], feature_group[, "intensity"], base.curve[, "base.curve"], all_rts)
+      # linear interpolation of  missing intensities and calculate the area for a single EIC
+      eic_area <- interpol.area(feature_group[, "rt"], feature_group[, "intensity"], base.curve[, "base.curve"], all_diff_mean_rts)
       rt_peak_shape <- c(median(feature_group[, "mz"]), median(feature_group[, "rt"]), sd(feature_group[, "rt"]), sd(feature_group[, "rt"]), eic_area)
       peak_parameters <- rbind(peak_parameters, rt_peak_shape)
     }
     if (num_features < 2) {
-      time_weights <- all_rts[which(base.curve[, "base.curve"] %in% feature_group[2])]
+      time_weights <- all_diff_mean_rts[which(base.curve[, "base.curve"] %in% feature_group[2])]
       rt_peak_shape <- c(feature_group[1], feature_group[2], NA, NA, feature_group[3] * time_weights)
       peak_parameters <- rbind(peak_parameters, rt_peak_shape)
     }
+    # application of selected model and method  
     if (num_features > 10) {
+      # find bandwidth for these particular range
       rt_range <- range(feature_group[, "rt"])
       bw <- min(max(bandwidth * (max(rt_range) - min(rt_range)), min_bandwidth), max_bandwidth)
       bw <- seq(bw, 2 * bw, length.out = 3)
@@ -881,7 +956,7 @@ prof.to.features <- function(profile,
 
       rt_profile <- compute_chromatographic_profile(feature_group, base.curve)
       if (shape_model == "Gaussian") {
-        rt_peak_shape <- compute_gaussian_peak_shape(rt_profile, moment_power, bw, component_eliminate, BIC_factor, aver_diff)
+        rt_peak_shape <- compute_gaussian_peak_shape(rt_profile, bw, component_eliminate, BIC_factor, aver_diff)
       } else {
         rt_peak_shape <- bigauss.mix(rt_profile, sigma_ratio_lim = sigma_ratio_lim, bw = bw, moment_power = moment_power, peak_estim_method = peak_estim_method, eliminate = component_eliminate, BIC_factor = BIC_factor)$param[, c(1, 2, 3, 5)]
       }
diff --git a/R/utils.R b/R/utils.R
index 023aace..c936d9f 100644
--- a/R/utils.R
+++ b/R/utils.R
@@ -29,7 +29,7 @@ register_functions_to_cluster <- function(cluster) {
         'bigauss.esti',
         'rev_cum_sum',
         'compute_bounds',
-        'validate_inputs',
+        'validate_model_method_input',
         'preprocess_bandwidth',
         'preprocess_profile',
         'compute_gaussian_peak_shape',