diff --git a/DESCRIPTION b/DESCRIPTION
index 9608873d2fbf56e8de2ec253a9481ced579d8497..0d901fb5a1dd2c0f0650027458f33c12c3c69b91 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -2,8 +2,8 @@ Package: gtseries
 Encoding: UTF-8
 Type: Package
 Title: Time Series Analysis for Paleoclimate and Cyclostratigraphy Purposes
-Version: 1.05
-Date: 2023-12-15
+Version: 1.3
+Date: 2025-02-19
 Author: Michel Crucifix [aut, cre] (<https://orcid.org/0000-0002-3437-4911>)
 Maintainer: Michel Crucifix <michel.crucifix@uclouvain.be>
 Description: Collection of time series analysis tools, focused on spectral decomposition and reconstruction, with an original focus on astronomical theory of paleoclimates and cyclostratigraphy. 
diff --git a/R/develop.R b/R/develop.R
index 311c55a3081e4459db5ff49270b00c282b98dcab..d599ebefac337092371be548138606f9bd9028c7 100644
--- a/R/develop.R
+++ b/R/develop.R
@@ -64,8 +64,8 @@ develop.discreteSpectrum  <- function(M, start = NULL, end = NULL, deltat = NULL
  if (!is.null(maxfreq)) nfreq <- min(nfreq, maxfreq) 
 
  if (timesIsATseries){
-   reconstructed <- lapply(seq(nfreq), function(i) ts( M$Amp[i] * dfunction(M$Freq[i] * times + M$Phase[i]), start=start, deltat=deltat) )} 
- else {
+   reconstructed <- lapply(seq(nfreq), function(i) ts( M$Amp[i] * dfunction(M$Freq[i] * times + M$Phase[i]), start=start, deltat=deltat) )
+ } else {
  reconstructed <- sapply(seq(nfreq), function(i) M$Amp[i] * dfunction(M$Freq[i] * times + M$Phase[i]) )
  }
 
@@ -78,11 +78,12 @@ if ( sum ) {
    else
      reconstructed <- apply(reconstructed, 1 , sum)  + trend * times + shift
  } else if (timesIsATseries) { 
-   reconstructed <- apply(reconstructed, 2, function(x) ts(x,start=start, deltat=deltat)) 
+   reconstructed <- lapply(reconstructed, function(x) ts(x,start=start, deltat=deltat)) 
 } 
  return(reconstructed)
 }
 
+#' Discrete spectrum class
 #' @rdname discreteSpectrum
 #' @export
 as.data.frame.discreteSpectrum <- function(x) {data.frame(Freq=x$Freq, Amp=x$Amp, Phases=x$Phases)}
@@ -125,8 +126,9 @@ lines.discreteSpectrum <- function (M,...){
 #' @rdname discreteSpectrum
 #' @export
 print.discreteSpectrum <- function (M,...){
-  N <- nrows(as.data.frame(M))
-  print.data.frame(cbind(as.data.frame(M[:,min(10,N)]), Period=2*pi/M$Freq))
+  N <- nrow(as.data.frame(M))
+  to_print <- seq(min(10,N))
+  print.data.frame(cbind(as.data.frame(M)[to_print, ], Period=2*pi/M$Freq[to_print]))
   if (N > 10) print(sprintf("... + %d other rows \n", N-10))
 }
 
diff --git a/man/develop.discreteSpectrum.Rd b/man/develop.discreteSpectrum.Rd
index 9a5f4a1fae1ffc5eb26936d7d73b58bed9ff5269..69016dcea573c2dbea6d637962c5e8a4c6808383 100644
--- a/man/develop.discreteSpectrum.Rd
+++ b/man/develop.discreteSpectrum.Rd
@@ -9,9 +9,11 @@
   start = NULL,
   end = NULL,
   deltat = NULL,
-  times,
+  times = NULL,
   dfunction = cos,
-  sum = TRUE
+  maxfreq = NULL,
+  sum = TRUE,
+  trendshift = TRUE
 )
 }
 \arguments{
@@ -23,6 +25,8 @@
 
 \item{sum}{: TRUE if user wants to sum components %in% the reconstruction}
 
+\item{trendshift}{: TRUE if user wants to account for trend and shift encoded in the discreteSpectrum object}
+
 \item{times:}{if supplied, times of the decomposition}
 
 \item{start:}{if supplied, overrides time and will generate a time series with start and deltat, which must then