diff --git a/NAMESPACE b/NAMESPACE
index a14daa17..27ab76b2 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -2,6 +2,12 @@
 
 S3method(D1,fFIT)
 S3method(D2,fFIT)
+S3method(PhenoDeriv,default)
+S3method(PhenoDeriv,fFIT)
+S3method(PhenoGu,default)
+S3method(PhenoGu,fFIT)
+S3method(PhenoTrs,default)
+S3method(PhenoTrs,fFIT)
 S3method(curvature,fFIT)
 S3method(get_GOF,fFITs)
 S3method(get_GOF,list)
@@ -9,12 +15,14 @@ S3method(get_fitting,fFITs)
 S3method(get_fitting,list)
 S3method(get_param,fFITs)
 S3method(get_param,list)
+S3method(get_pheno,fFITs)
+S3method(get_pheno,list)
+S3method(get_pheno,rfit)
 S3method(plot,fFITs)
 S3method(print,fFIT)
 export("%<-%")
 export("%<>%")
 export("%>%")
-export(.Greenup)
 export(D1)
 export(D2)
 export(FitDL.AG)
@@ -34,6 +42,7 @@ export(PhenoKl)
 export(PhenoTrs)
 export(R2_sign)
 export(add_HeadTail)
+export(brks2rfit)
 export(check_input)
 export(check_season_dt)
 export(check_season_list)
@@ -68,7 +77,6 @@ export(get_fitting)
 export(get_options)
 export(get_param)
 export(get_pheno)
-export(get_pheno.fFITs)
 export(init_AG)
 export(init_AG2)
 export(init_Beck)
@@ -182,6 +190,7 @@ importFrom(lubridate,yday)
 importFrom(lubridate,year)
 importFrom(lubridate,ymd)
 importFrom(methods,is)
+importFrom(purrr,as_mapper)
 importFrom(purrr,map)
 importFrom(purrr,map_df)
 importFrom(purrr,map_dfc)
diff --git a/NEWS.md b/NEWS.md
index e7f54c8c..ae0eb9c1 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,12 +1,4 @@
-# phenofit 0.3.4 (2021-11-14)
-
-- Fix the bug of `findpeaks`, which lead to sharp changed growing season failed to detect.
-- Fix the bug of `PhenoKlos`, where `minpeakheight` not work in previous version.
-- Fix the bug of `check_season_dt`, where `peak` might be able to greater than `end`.
-- Remove the parameter `check_season_dt` in `removeClosedExtreme`, which might eliminate
-  good extreme values.
-  
-# phenofit 0.3.3 (2021-10-26)
+# phenofit 0.3.3 (2021-11-14)
 
 - Fix typo error in curvefits' document.
 
@@ -30,6 +22,19 @@
   used the strict mathematical solution to find the extreme values in the curve
   of curvature's change rate.
 
+** MAJOR updates to improve multi-GS phenology extraction **
+
+- Fix the bug of `findpeaks`, which lead to sharp changed growing season failed to detect.
+
+- Fix the bug of `PhenoKlos`, where `minpeakheight` not work in previous version.
+
+- Fix the bug of `check_season_dt`, where `peak` might be able to greater than `end`.
+
+- Remove the parameter `check_season_dt` in `removeClosedExtreme`, which might eliminate
+  good extreme values.
+  
+- add `get_pheno.rfit` to extract vegetation phenology from rough fitting directly.
+
 # phenofit 0.3.2 (2021-10-15)
 
 - Parameters of `season_mov` and `curvefits`  are wrapped into options. Scripts of phenofit v2.0 will not work anymore.
diff --git a/R/PhenoDeriv.R b/R/PhenoDeriv.R
new file mode 100644
index 00000000..b2f95d76
--- /dev/null
+++ b/R/PhenoDeriv.R
@@ -0,0 +1,108 @@
+#' Phenology extraction in Derivative method (DER)
+#' 
+#' @inheritParams D
+#' @inheritParams PhenoTrs
+#' @param der1 the first order difference
+#' 
+#' @references 
+#' 1. Filippa, G., Cremonese, E., Migliavacca, M., Galvagno, M., Forkel, M.,
+#'    Wingate, L., … Richardson, A. D. (2016). Phenopix: A R package for
+#'    image-based vegetation phenology. Agricultural and Forest Meteorology,
+#'    220, 141–150. \doi{10.1016/j.agrformet.2016.01.006}
+#' 
+#' @seealso [PhenoTrs()], [PhenoGu()], [PhenoKl()]
+#' @export
+PhenoDeriv <- function(x, t, ...) UseMethod("PhenoDeriv", x)
+
+#' @rdname PhenoDeriv
+#' @export
+PhenoDeriv.fFIT <- function(x, t = NULL,
+    analytical = TRUE, smoothed.spline = FALSE, ...) {
+    if (!is.null(x$tout)) t <- x$tout
+    values <- last2(x$zs)
+    der1 <- D1.fFIT(x, t, analytical, smoothed.spline)
+    PhenoDeriv.default(values, t, der1, ...)
+}
+
+#' @rdname PhenoDeriv
+#' @param show.legend whether show figure lelend?
+#' @export
+PhenoDeriv.default <- function(x, t, der1,
+    IsPlot = TRUE, show.legend = TRUE, ...)
+{
+    PhenoNames <- c("SOS", "POS", "EOS")
+    metrics <- setNames(rep(NA, 3), c("sos", "pos", "eos")) # template
+    n      <- length(t)
+
+    # get peak of season position
+    half.season <- median(which.max(x)) # deal with multiple pos x
+    pos <- t[half.season]
+
+    if (all(is.na(x))) return(metrics)
+    if (half.season < 5 || half.season > (n - 5)) return(metrics)
+
+    # get SOS and EOS according to first order derivative
+    # fixed 20180510, ±5 to make sure sos and eos are not in the POP.
+    # I_sos <- median(which.max(der1[1:(half.season - 5)]))
+    # I_eos <- median(which.min(der1[(half.season+5):length(der1)])) + half.season
+
+    # der.sos (eos) is impossible to occur near the pop.
+    nmin  <- 5
+    I_sos <- findpeaks( der1[1:(half.season - 5)],
+        nups=nmin, ndowns=nmin, npeaks=1, sortstr=TRUE)$X$pos
+    I_eos <- findpeaks(-der1[(half.season+5):length(der1)],
+        nups=nmin, ndowns=nmin, npeaks=1, sortstr=TRUE)$X$pos + half.season + 4
+
+    # if half.season > length(der1), error will be occur
+    sos <- t[I_sos]
+    eos <- t[I_eos]
+    if (is_empty(sos)) {
+        I_sos <- findpeaks( der1[1:(half.season - 5)],
+            nups=nmin, ndowns=0, npeaks=1, sortstr=TRUE)$X$pos
+        sos <- t[I_sos]
+        if (is_empty(sos)) sos <- round((t[1] + pos)/2)
+    }
+    if (is_empty(eos)) {
+        I_eos <- findpeaks(-der1[(half.season+5):length(der1)],
+            nups=nmin, ndowns=0, npeaks=1, sortstr=TRUE)$X$pos + half.season + 4
+        eos <- t[I_eos]
+        if (is_empty(eos)) eos <- round((t[n] + pos)/2)
+    }
+    # greenup <- .Greenup(x)
+    # if (is.na(sos)) {
+    #     if (all(!rm_empty(greenup))) sos = pos
+    # }
+    # if (is.na(eos)) {
+    #     if (all(rm_empty(greenup)))  eos = pos
+    # }
+    metrics <- c(sos = sos, pos = pos, eos = eos)#, los = los
+
+    if (IsPlot){
+        main  <- ifelse(all(par("mar") == 0), "", "DER")
+        PhenoPlot(t, x, main = main, ...)
+        if (show.legend) legend('topright', c("f(t)'"), lty = 2, col = "black", bty='n')
+
+        abline(v = c(sos, eos), col = colors[c(1, 4)], lwd = linewidth)
+        abline(v = pos, col ="blue", lty = 1, lwd = linewidth)
+
+        A <- diff(range(x))
+        I_metrics <- match(metrics, t)
+        if (all(is.na(I_metrics))) {
+            ylons <- I_metrics
+        }else{
+            ylons <- x[I_metrics] + c(1, -2, 1)*0.1*A
+        }
+        xlons <- metrics + c(-1, 1, 1)*5
+        xlons[xlons < min(t)] <- min(t)
+        xlons[xlons > max(t)] <- max(t)
+
+        I <- c(1); text(xlons[I], ylons[I], PhenoNames[I], col = colors[I], adj = c(1, 0))
+        I <- 2:3 ; text(xlons[I], ylons[I], PhenoNames[I], col = c("blue", colors[4]), adj = c(0, 0))
+
+        #der1 last plot
+        op <- par(new = TRUE)
+        plot(t, der1, type= "l", lty = 2, lwd = linewidth,
+             col = "black", axes = FALSE)
+    }
+    return(metrics)
+}
diff --git a/R/PhenoExtractMeth.R b/R/PhenoExtractMeth.R
deleted file mode 100644
index 247901dd..00000000
--- a/R/PhenoExtractMeth.R
+++ /dev/null
@@ -1,490 +0,0 @@
-# colors    <- c("blue", "green3", "orange", "red")
-colors <- c("green3", "darkgreen", "darkorange3", "red")
-linewidth <- 1.2
-
-#' @name PhenoExtractMeth
-#' @title Phenology Extraction methods
-#'
-#' @inheritParams D
-#'
-#' @param fFIT `fFIT` object returned by [optim_pheno()].
-#' @param t `date` or `doy` vector, with the same length as `ypred`. This parameter
-#' is for the Julia version `curvefits`.
-#'
-#' @param approach to be used to calculate phenology metrics.
-#' 'White' (White et al. 1997) or 'Trs' for simple threshold.
-#' @param trs threshold to be used for approach "Trs", in (0, 1).
-#' @param IsPlot whether to plot?
-#' @param show.lgd whether show figure lelend?
-#' @param ... other parameters to PhenoPlot
-#'
-#' @examples
-#' library(phenofit)
-#' # simulate vegetation time-series
-#' fFUN = doubleLog.Beck
-#' par  = c( mn  = 0.1 , mx  = 0.7 , sos = 50 , rsp = 0.1 , eos = 250, rau = 0.1)
-#'
-#' t    <- seq(1, 365, 8)
-#' tout <- seq(1, 365, 1)
-#' y <- fFUN(par, t)
-#'
-#' methods <- c("AG", "Beck", "Elmore", "Gu", "Zhang") # "Klos" too slow
-#' fFITs <- curvefit(y, t, tout, methods)
-#' fFIT  <- fFITs$model$AG
-#'
-#' par(mfrow = c(2, 2))
-#' PhenoTrs(fFIT)
-#' PhenoDeriv(fFIT)
-#' PhenoGu(fFIT)
-#' PhenoKl(fFIT)
-NULL
-
-#' @description
-#' \itemize{
-#' \item `PhenoTrs` Threshold method
-#' \item `PhenoDeriv` Derivative method
-#' \item `PhenoGu` Gu method
-#' \item `PhenoKl` Inflection method
-#' }
-#'
-#' @param asymmetric If true, background value in spring season and autumn season
-#' is regarded as different.
-#'
-#' @rdname PhenoExtractMeth
-#' @export
-PhenoTrs <- function(fFIT, t = NULL, approach = c("White", "Trs"), trs = 0.5, #, min.mean = 0.1
-    asymmetric = TRUE,
-    IsPlot = TRUE, ...)
-{
-    metrics <- c(sos = NA, eos = NA)
-    if (!is.null(fFIT$tout)) t <- fFIT$tout
-
-    values <- last2(fFIT$zs)
-    n      <- length(t)
-
-    # get peak of season position
-    half.season <- median(which.max(values)) %>% round() # + 20, half season + 20 was unreasonable
-    pos <- t[half.season]
-
-    if (all(is.na(values))) return(metrics)
-    if (half.season < 5 || half.season > (n - 5)) return(metrics)
-
-    # get statistical values
-    # n    <- t[length(t)]
-    # avg  <- mean(x, na.rm = TRUE)
-
-    # avg2 <- mean(x2[x2 > min.mean], na.rm = TRUE)
-    peak <- max(values, na.rm = TRUE)
-
-    if (asymmetric) {
-        mn_a <- min(values[1:half.season], na.rm = T)
-        mn_b <- min(values[-(1:half.season)], na.rm = T)
-
-        mn <- c(rep(mn_a, half.season),
-            rep(mn_b, n - half.season))
-    } else {
-        mn   <- min(values, na.rm = TRUE)
-    }
-
-    ampl <- peak - mn
-
-    # select (or scale) values and thresholds for different methods
-    approach <- approach[1]
-    if (approach == "White") {
-        # scale annual time series to 0-1
-        ratio   <- (values - mn)/ampl
-        # trs   <- 0.5
-        trs.low <- trs - 0.05
-        trs.up  <- trs + 0.05
-    }
-    if (approach == "Trs") {
-        ratio   <- values
-        a       <- diff(range(ratio, na.rm = TRUE)) * 0.05
-        trs.low <- trs - a
-        trs.up  <- trs + a
-    }
-
-    greenup <- .Greenup(ratio)
-    ## distinguish the first half year and second year
-    # select time where SOS and EOS are located (around trs value)
-    bool <- ratio >= trs.low & ratio <= trs.up
-
-    # get SOS, EOS, LOS
-    # fixed 2017-01-04, according to TP phenology property
-    # sos <- round(median(sose os[greenup & bool], na.rm = TRUE))
-    # eos <- round(median(soseos[!greenup & bool], na.rm = TRUE))
-
-    sos <- round(median(t[ greenup & bool & t < pos], na.rm = TRUE))
-    eos <- round(median(t[!greenup & bool & t > pos], na.rm = TRUE))
-
-    # greenup <- .Greenup(ratio)
-    if (is.na(sos + eos)) {
-        bool2 <- ratio >= (trs-0.1) & ratio <= (trs+0.1)
-        if (is.na(sos)) sos <- round(median(t[ greenup & bool2 & t < pos], na.rm = TRUE))
-        if (is.na(eos)) eos <- round(median(t[!greenup & bool2 & t > pos], na.rm = TRUE))
-
-        # if still is NA
-        if (is.na(sos)) sos <- round((t[1] + pos)/2)
-        if (is.na(eos)) eos <- round((t[n] + pos)/2)
-    }
-    # monotonous = length(unique(greenup[2:n]) %>% rm_empty()) == 1
-    # los <- eos - sos#los < 0 indicate that error
-    # los[los < 0] <- n + (eos[los < 0] - sos[los < 0])
-
-    # get MGS, MSP, MAU
-    # mgs <- mean(x[ratio > trs], na.rm = TRUE)
-    # msp <- mau <- NA
-    # if (!is.na(sos)) {
-    #     id <- (sos - 10):(sos + 10)
-    #     id <- id[(id > 0) & (id < n)]
-    #     msp <- mean(x[which(index(x) %in% id == TRUE)], na.rm = TRUE)
-    # }
-    # if (!is.na(eos)) {
-    #     id <- (eos - 10):(eos + 10)
-    #     id <- id[(id > 0) & (id < n)]
-    #     mau <- mean(x[which(index(x) %in% id == TRUE)], na.rm = TRUE)
-    # }
-    # metrics <- c(sos = sos, eos = eos, los = los, pos = pos, mgs = mgs,
-    #   rsp = NA, rau = NA, peak = peak, msp = msp, mau = mau)
-    metrics <- c(sos = sos, eos = eos)#, los = los
-
-    if (IsPlot){
-        main   <- ifelse(all(par("mar") == 0), "", sprintf("TRS%d", trs*10))
-        PhenoPlot(t, values, main = main, ...)
-
-        lines(t, trs*ampl + mn, lwd = linewidth)
-        # lines(t, trs.low*ampl + mn, lty = 2, lwd = linewidth)
-        # lines(t, trs.up*ampl + mn, lty = 2, lwd = linewidth)
-        abline(v = metrics, col = colors[c(1, 4)], lwd = linewidth)
-        text(metrics[1] - 5, min(trs + 0.15, 1)*ampl[1] + mn[1], "SOS", col = colors[1], adj = c(1, 0))
-        text(metrics[2] + 5, min(trs + 0.15, 1)*last(ampl) + last(mn), "EOS", col = colors[4], adj = c(0, 0))
-    }
-    return(metrics)
-    ### The function returns a vector with SOS, EOS, LOS, POS, MGS, rsp, rau, PEAK, MSP and MAU. }
-}
-
-# identify greenup or dormancy(brown) period
-#' @export
-.Greenup <- function(x, ...) {
-    ratio.deriv <- c(NA, diff(x))
-    greenup     <- rep(NA, length(x))
-    greenup[ratio.deriv > 0] <- TRUE
-    greenup[ratio.deriv < 0] <- FALSE
-    return(greenup)
-}
-
-#' @inheritParams D
-#' @inheritParams PhenoTrs
-#'
-#' @rdname PhenoExtractMeth
-#' @export
-PhenoDeriv <- function(fFIT, t = NULL,
-    analytical = TRUE, smoothed.spline = FALSE,
-    IsPlot = TRUE, show.lgd = TRUE, ...)
-{
-    PhenoNames <- c("SOS", "POS", "EOS")
-    metrics <- setNames(rep(NA, 3), c("sos", "pos", "eos")) # template
-
-    # t      <- fFIT$tout
-    if (!is.null(fFIT$tout)) t <- fFIT$tout
-    values <- last2(fFIT$zs)
-    n      <- length(t)
-
-    # get peak of season position
-    half.season <- median(which.max(values)) # deal with multiple pos values
-    pos <- t[half.season]
-
-    if (all(is.na(values))) return(metrics)
-    if (half.season < 5 || half.season > (n - 5)) return(metrics)
-
-    der1   <- D1.fFIT(fFIT, t, analytical, smoothed.spline)
-    # get SOS and EOS according to first order derivative
-    # fixed 20180510, ±5 to make sure sos and eos are not in the POP.
-    # I_sos <- median(which.max(der1[1:(half.season - 5)]))
-    # I_eos <- median(which.min(der1[(half.season+5):length(der1)])) + half.season
-
-    # der.sos (eos) is impossible to occur near the pop.
-    nmin  <- 5
-    I_sos <- findpeaks( der1[1:(half.season - 5)],
-        nups=nmin, ndowns=nmin, npeaks=1, sortstr=TRUE)$X$pos
-    I_eos <- findpeaks(-der1[(half.season+5):length(der1)],
-        nups=nmin, ndowns=nmin, npeaks=1, sortstr=TRUE)$X$pos + half.season + 4
-
-    # if half.season > length(der1), error will be occur
-    sos <- t[I_sos]
-    eos <- t[I_eos]
-    if (is_empty(sos)) {
-        I_sos <- findpeaks( der1[1:(half.season - 5)],
-            nups=nmin, ndowns=0, npeaks=1, sortstr=TRUE)$X$pos
-        sos <- t[I_sos]
-        if (is_empty(sos)) sos <- round((t[1] + pos)/2)
-    }
-    if (is_empty(eos)) {
-        I_eos <- findpeaks(-der1[(half.season+5):length(der1)],
-            nups=nmin, ndowns=0, npeaks=1, sortstr=TRUE)$X$pos + half.season + 4
-        eos <- t[I_eos]
-        if (is_empty(eos)) eos <- round((t[n] + pos)/2)
-    }
-    # greenup <- .Greenup(values)
-    # if (is.na(sos)) {
-    #     if (all(!rm_empty(greenup))) sos = pos
-    # }
-    # if (is.na(eos)) {
-    #     if (all(rm_empty(greenup)))  eos = pos
-    # }
-    metrics <- c(sos = sos, pos = pos, eos = eos)#, los = los
-
-    if (IsPlot){
-        main  <- ifelse(all(par("mar") == 0), "", "DER")
-        PhenoPlot(t, values, main = main, ...)
-        if (show.lgd) legend('topright', c("f(t)'"), lty = 2, col = "black", bty='n')
-
-        abline(v = c(sos, eos), col = colors[c(1, 4)], lwd = linewidth)
-        abline(v = pos, col ="blue", lty = 1, lwd = linewidth)
-
-        A <- diff(range(values))
-        I_metrics <- match(metrics, t)
-        if (all(is.na(I_metrics))) {
-            ylons <- I_metrics
-        }else{
-            ylons <- values[I_metrics] + c(1, -2, 1)*0.1*A
-        }
-        xlons <- metrics + c(-1, 1, 1)*5
-        xlons[xlons < min(t)] <- min(t)
-        xlons[xlons > max(t)] <- max(t)
-
-        I <- c(1); text(xlons[I], ylons[I], PhenoNames[I], col = colors[I], adj = c(1, 0))
-        I <- 2:3 ; text(xlons[I], ylons[I], PhenoNames[I], col = c("blue", colors[4]), adj = c(0, 0))
-
-        #der1 last plot
-        op <- par(new = TRUE)
-        plot(t, der1, type= "l", lty = 2, lwd = linewidth,
-             col = "black", axes = FALSE)
-    }
-    return(metrics)
-}
-
-
-#' @inheritParams PhenoTrs
-#' @rdname PhenoExtractMeth
-#' @export
-PhenoGu <- function(fFIT, t = NULL,
-    analytical = TRUE, smoothed.spline = FALSE,
-    IsPlot = TRUE, ...)
-{
-    PhenoNames <- c("UD", "SD", "DD", "RD")
-    metrics <- setNames(rep(NA, 4), c("UD", "SD", "DD", "RD"))
-
-    # t      <- fFIT$tout
-    if (!is.null(fFIT$tout)) t <- fFIT$tout
-    values <- last2(fFIT$zs)
-    n      <- length(t)
-
-    # get peak of season position
-    half.season <- median(which.max(values)) # deal with multiple pos values
-    pos <- t[half.season]
-
-    if (all(is.na(values))) return(metrics)
-    if (half.season < 5 || half.season > (n - 5)) return(metrics)
-
-    der1   <- D1.fFIT(fFIT, t, analytical, smoothed.spline)
-    # get SOS and EOS according to first order derivative
-    sos.index <- median(which.max(der1[1:(half.season-5)]))
-    eos.index <- median(which.min(der1[(half.season+5):length(der1)])) + half.season
-    sos <- t[sos.index]
-    eos <- t[eos.index]
-
-    rsp <- der1[sos.index] # rate of spring, also known as peak recovery rate
-    rau <- der1[eos.index] # rate of autumn, peak senescence rate
-    VI.rsp <- values[sos.index] # VI index of corresponding date
-    VI.rau <- values[eos.index]
-
-    # Gu Phenology Extraction method also quite rely on first order derivative
-    rl.b <- VI.rsp - rsp*sos # interception of recovery line
-    sl.b <- VI.rau - rau*eos # interception of recovery line
-
-    baseline <- min(values, na.rm = TRUE)
-    maxline  <- max(values, na.rm = TRUE)
-
-    ## y = kx + b; x = (y - b)/k
-    UD <- (baseline - rl.b)/rsp # upturn day, intersection of rl and x axis
-    SD <- (maxline  - rl.b)/rsp # stabilization day, intersection of maxline and rl
-    DD <- (maxline  - sl.b)/rau # downturn day, intersection of maxline and sl
-    RD <- (baseline - sl.b)/rau # recession day, intersection of sl and x axis
-
-    ## subset data between SD and DD
-    sub.time <- t[t >= SD & t <= DD]
-    sub.gcc  <- values[t >= SD & t <= DD]
-
-    ## compute a linear fit
-    if (length(sub.time) > 3) {
-        X <- cbind(1, sub.time)
-        plateau.lm        <- .lm.fit(X, sub.gcc)
-        plateau.slope     <- plateau.lm$coefficients[2]
-        plateau.intercept <- plateau.lm$coefficients[1]
-        # y1 = rau*t + sl.b
-        # y2 = k2t + b2
-        # so, t = (sl.b - b2)/(k2 -rau)
-        DD <- (sl.b - plateau.intercept) / (plateau.slope - rau)
-    }else{
-        plateau.slope     <- NA
-        plateau.intercept <- NA
-    }
-    ## calculate area under the curve
-    # cut.x <- days[which(days>=UD & days<=RD)]
-    # cut.y <- offset.y[which(days>=UD & days<=RD)]
-    # the.fun <- function(t) {eval(retrieved.formula, envir=as.list(params))}
-
-    ## Avoid NA values by using `clamp` instead
-    # sapply(function(x){ ifelse(x < min(t) || x > max(t), NA, x) }) %>%
-    metrics   <- round(c(UD, SD, DD, RD)) %>%
-        sapply(clamp, lims = c(t[1], t[n])) %>%
-        set_names(PhenoNames)
-    # c("UD", "SD", "DD", "RD", "maxline", "baseline", "rsp", "rau", "plateau.slope")
-    # c(pheno, maxline, baseline, rsp, rau, plateau.slope)
-
-    if (IsPlot){
-        main   <- ifelse(all(par("mar") == 0), "", "Gu")
-        PhenoPlot(t, values, main = main, ...)
-
-        abline(rl.b, rsp, col = "blue", lty=  2, lwd = linewidth,)
-        abline(sl.b, rau, col = "red" , lty=  2, lwd = linewidth,)
-        # any na values input to abline will lead to error
-        if (all(!is.na(c(plateau.slope, plateau.intercept))))
-            abline(a = plateau.intercept, b = plateau.slope,
-                lty = 2, lwd = linewidth, col = "darkgreen")
-
-        abline(h = c(maxline, baseline), lty = 2, lwd = linewidth,)
-        abline(v = metrics[1:4],  col = colors, lwd = linewidth,)
-
-        A <- diff(range(values))
-        I_metrics <- match(metrics, t)
-        if (all(is.na(I_metrics))) {
-            ylons <- I_metrics
-        }else{
-            ylons <- values[I_metrics] + c(1, -3, -3, 1)*0.1*A
-        }
-
-        xlons <- metrics[1:4] + c(-1, -1, 1, 1)*5
-        xlons[xlons < min(t)] <- min(t)
-        xlons[xlons > max(t)] <- max(t)
-
-        I <- c(1, 2); text(xlons[I], ylons[I], PhenoNames[I], col = colors[I], adj = c(1, 0))
-        I <- c(3, 4); text(xlons[I], ylons[I], PhenoNames[I], col = colors[I], adj = c(0, 0))
-    }
-    return(metrics)
-}
-
-#' @inheritParams PhenoTrs
-#' @rdname PhenoExtractMeth
-#' @export
-PhenoKl <- function(fFIT, t = NULL,
-    analytical = TRUE, smoothed.spline = FALSE,
-    IsPlot = TRUE, show.lgd = TRUE, ...)
-{
-    PhenoNames <- c("Greenup", "Maturity", "Senescence", "Dormancy")
-    metrics <- setNames(rep(NA, 4), PhenoNames)
-
-    if (!is.null(fFIT$tout)) t <- fFIT$tout
-    values <- last2(fFIT$zs)
-    n      <- length(t)
-    xlim   <- range(t)
-
-    # get peak of season position
-    half.season <- median(which.max(values)) # + 20, half season + 20 was unreasonable
-    pos <- t[half.season]
-
-    if (all(is.na(values))) return(metrics)
-    if (half.season < 5 || half.season > (n - 5)) return(metrics)
-
-    derivs <- curvature.fFIT(fFIT, t, analytical, smoothed.spline)
-    k      <- derivs$k
-    # define cutoff date for spline functions
-
-    # x <- ifelse(uncert==TRUE, x$uncertainty, x$fit)
-    # if (length(which(is.na(k) == TRUE)) != 0 | length(which(is.infinite(k) == TRUE)) != 0) {
-
-    # If have no NA and infinite values
-    if (!(any(is.na(k)) || any(is.infinite(k)))) {
-        spline.k <- smooth.spline(k, df = 0.1 * length(k))
-        der.k    <- predict(spline.k, d = 1)$y
-        der.k2   <- predict(spline.k, d = 2)$y
-
-        # der.k  <- c(NA, diff(k))
-        # der.k2 <- c(NA, NA, diff(k, differences = 2))
-        ## find maxima of derivative of k ## split season
-        dist_fromPeak <- 0 # days
-        asc.k   <- try(der.k[1:(half.season - dist_fromPeak)]) # k   of first half year
-        asc.k.d <- try(t[1:(half.season - dist_fromPeak)])
-        # doy of first half year
-        des.k   <- try(der.k[(half.season + dist_fromPeak):length(k)])
-        des.k.d <- try(t[(half.season + dist_fromPeak):length(k)])
-
-        A <- range(der.k, na.rm = TRUE) %>% diff()
-        minpeakheight = A*0.025
-
-        # first half maximum local values of k'
-        # minpeakdistance in the unit of day
-        pos <- findpeaks(asc.k, npeaks = 2, ndowns = 0, sortstr = TRUE,
-                         minpeakdistance = 15, minpeakheight = minpeakheight)$X$pos
-        pos <- sort(pos)
-        pos <- c(rep(NA, 2 - length(pos)), pos) #at least two values
-        I_asc <- asc.k.d[pos]
-        if (all(is.na(pos))) {
-            I_asc <- pos
-        } else {
-            I_asc <- asc.k.d[pos]
-        }
-
-        # second half minimum local values of k'
-        pos <- findpeaks(-des.k, npeaks = 2, sortstr = TRUE, 
-            minpeakdistance = 15, minpeakheight = minpeakheight)$X$pos
-        pos <- sort(pos);
-        pos <- c(pos, rep(NA, 2 - length(pos)))
-        if (all(is.na(pos))) {
-            I_dec <- pos
-        } else {
-            I_dec <- des.k.d[pos]
-        }
-        metrics <- c(I_asc, I_dec)
-    }
-
-    if (IsPlot){
-        main   <- ifelse(all(par("mar") == 0), "", "Zhang (Curvature Rate)")
-        A <- diff(range(values))
-
-        I_metrics <- match(metrics, t)
-        if (all(is.na(I_metrics))) {
-            ylons <- I_metrics
-        }else{
-            ylons <- values[I_metrics] + c(-1, -1, -1, 1)*0.1*A
-        }
-
-        xlons  <- metrics + c(1, -1, 1, -1) * 5
-        xlons[xlons < min(t)] <- min(t)
-        xlons[xlons > max(t)] <- max(t)
-        # plotrix::twoord.plot(t, values, t, der.k,
-        #                      main = main, type= c("p", "b"), xlim = xlim,
-        #                      rcol = "black", lcol = "grey60", lpch = 20, rpch = 1,
-        #                      rytickpos = NULL)
-        PhenoPlot(t, values, main = main, ...)
-        if (show.lgd){
-            legend('topright', c("K'"), lty = c(3), col = c("black"), bty='n') ##pch =c(20, 1),
-        }
-
-        pos     <- t[half.season]
-        abline(v = metrics, col = colors, lwd = linewidth)
-        # abline(v = pos, col ="darkgreen", lty = 1, lwd = linewidth)
-        # abline(v = pos + 20, col ="darkgreen", lty = 2, lwd = linewidth)
-        PhenoNames2 <- c("Greenup", "Maturity", "Senescence", "Dormancy")
-        # PhenoNames2 <- c("G", "M", "S", "D")
-
-        I <- c(1, 3); text(xlons[I], ylons[I], PhenoNames2[I], col = colors[I], adj = c(0, 0))
-        I <- c(2, 4); text(xlons[I], ylons[I], PhenoNames2[I], col = colors[I], adj = c(1, 0))
-        # der.k last plot
-        op <- par(new = TRUE)
-        plot(t, der.k, xlim = xlim, type= "l",
-             lty = 3, lwd = linewidth, col = "black", axes = TRUE) # cex = 1, pch = 20,
-    }
-    return(setNames(metrics, PhenoNames))
-}
diff --git a/R/PhenoGu.R b/R/PhenoGu.R
new file mode 100644
index 00000000..cecb2bdb
--- /dev/null
+++ b/R/PhenoGu.R
@@ -0,0 +1,141 @@
+#' Phenology extraction in GU method (GU)
+#'
+#' @inheritParams PhenoDeriv 
+#' @inherit PhenoTrs examples
+#' @param ... other parameters to [PhenoGu.default()] or [PhenoGu.fFIT()]
+#' 
+#' @references
+#' 1. Gu, L., Post, W. M., Baldocchi, D. D., Black, T. A., Suyker, A. E., Verma,
+#'    S. B., … Wofsy, S. C. (2009). Characterizing the Seasonal Dynamics of
+#'    Plant Community Photosynthesis Across a Range of Vegetation Types. In A.
+#'    Noormets (Ed.), Phenology of Ecosystem Processes: Applications in Global
+#'    Change Research (pp. 35–58). New York, NY: Springer New York.
+#'    \doi{10.1007/978-1-4419-0026-5_2}
+#' 2. Filippa, G., Cremonese, E., Migliavacca, M., Galvagno, M., Forkel, M.,
+#'    Wingate, L., … Richardson, A. D. (2016). Phenopix: A R package for
+#'    image-based vegetation phenology. Agricultural and Forest Meteorology,
+#'    220, 141–150. \doi{10.1016/j.agrformet.2016.01.006}
+#' 
+#' @return A numeric vector, with the elements of:
+#' - `UD`: upturn date
+#' - `SD`: stabilisation date
+#' - `DD`: downturn date
+#' - `RD`: recession date 
+#' @export
+PhenoGu <- function(x, t, ...) UseMethod("PhenoGu", x)
+
+#' @rdname PhenoGu
+#' @export
+PhenoGu.fFIT <- function(x, t = NULL,
+    analytical = TRUE, smoothed.spline = FALSE, ...) 
+{
+    if (!is.null(x$tout)) t <- x$tout
+    values <- last2(x$zs)
+    der1 <- D1.fFIT(x, t, analytical, smoothed.spline)
+    PhenoGu.default(values, t, der1, ...)
+}
+
+#' @inheritParams PhenoTrs
+#' @rdname PhenoGu
+#' @export
+PhenoGu.default <- function(x, t, der1,
+    IsPlot = TRUE, ...)
+{
+    PhenoNames <- c("UD", "SD", "DD", "RD")
+    metrics <- setNames(rep(NA, 4), c("UD", "SD", "DD", "RD"))
+
+    n      <- length(t)
+
+    # get peak of season position
+    half.season <- median(which.max(x)) # deal with multiple pos x
+    pos <- t[half.season]
+
+    if (all(is.na(x))) return(metrics)
+    if (half.season < 5 || half.season > (n - 5)) return(metrics)
+
+    # get SOS and EOS according to first order derivative
+    sos.index <- median(which.max(der1[1:(half.season-5)]))
+    eos.index <- median(which.min(der1[(half.season+5):length(der1)])) + half.season
+    sos <- t[sos.index]
+    eos <- t[eos.index]
+
+    rsp <- der1[sos.index] # rate of spring, also known as peak recovery rate
+    rau <- der1[eos.index] # rate of autumn, peak senescence rate
+    VI.rsp <- x[sos.index] # VI index of corresponding date
+    VI.rau <- x[eos.index]
+
+    # Gu Phenology Extraction method also quite rely on first order derivative
+    rl.b <- VI.rsp - rsp*sos # interception of recovery line
+    sl.b <- VI.rau - rau*eos # interception of recovery line
+
+    baseline <- min(x, na.rm = TRUE)
+    maxline  <- max(x, na.rm = TRUE)
+
+    ## y = kx + b; x = (y - b)/k
+    UD <- (baseline - rl.b)/rsp # upturn day, intersection of rl and x axis
+    SD <- (maxline  - rl.b)/rsp # stabilization day, intersection of maxline and rl
+    DD <- (maxline  - sl.b)/rau # downturn day, intersection of maxline and sl
+    RD <- (baseline - sl.b)/rau # recession day, intersection of sl and x axis
+
+    ## subset data between SD and DD
+    sub.time <- t[t >= SD & t <= DD]
+    sub.gcc  <- x[t >= SD & t <= DD]
+
+    ## compute a linear fit
+    if (length(sub.time) > 3) {
+        X <- cbind(1, sub.time)
+        plateau.lm        <- .lm.fit(X, sub.gcc)
+        plateau.slope     <- plateau.lm$coefficients[2]
+        plateau.intercept <- plateau.lm$coefficients[1]
+        # y1 = rau*t + sl.b
+        # y2 = k2t + b2
+        # so, t = (sl.b - b2)/(k2 -rau)
+        DD <- (sl.b - plateau.intercept) / (plateau.slope - rau)
+    }else{
+        plateau.slope     <- NA
+        plateau.intercept <- NA
+    }
+    ## calculate area under the curve
+    # cut.x <- days[which(days>=UD & days<=RD)]
+    # cut.y <- offset.y[which(days>=UD & days<=RD)]
+    # the.fun <- function(t) {eval(retrieved.formula, envir=as.list(params))}
+
+    ## Avoid NA x by using `clamp` instead
+    # sapply(function(x){ ifelse(x < min(t) || x > max(t), NA, x) }) %>%
+    metrics   <- round(c(UD, SD, DD, RD)) %>%
+        sapply(clamp, lims = c(t[1], t[n])) %>%
+        set_names(PhenoNames)
+    # c("UD", "SD", "DD", "RD", "maxline", "baseline", "rsp", "rau", "plateau.slope")
+    # c(pheno, maxline, baseline, rsp, rau, plateau.slope)
+
+    if (IsPlot){
+        main   <- ifelse(all(par("mar") == 0), "", "Gu")
+        PhenoPlot(t, x, main = main, ...)
+
+        abline(rl.b, rsp, col = "blue", lty=  2, lwd = linewidth,)
+        abline(sl.b, rau, col = "red" , lty=  2, lwd = linewidth,)
+        # any na x input to abline will lead to error
+        if (all(!is.na(c(plateau.slope, plateau.intercept))))
+            abline(a = plateau.intercept, b = plateau.slope,
+                lty = 2, lwd = linewidth, col = "darkgreen")
+
+        abline(h = c(maxline, baseline), lty = 2, lwd = linewidth,)
+        abline(v = metrics[1:4],  col = colors, lwd = linewidth,)
+
+        A <- diff(range(x))
+        I_metrics <- match(metrics, t)
+        if (all(is.na(I_metrics))) {
+            ylons <- I_metrics
+        }else{
+            ylons <- x[I_metrics] + c(1, -3, -3, 1)*0.1*A
+        }
+
+        xlons <- metrics[1:4] + c(-1, -1, 1, 1)*5
+        xlons[xlons < min(t)] <- min(t)
+        xlons[xlons > max(t)] <- max(t)
+
+        I <- c(1, 2); text(xlons[I], ylons[I], PhenoNames[I], col = colors[I], adj = c(1, 0))
+        I <- c(3, 4); text(xlons[I], ylons[I], PhenoNames[I], col = colors[I], adj = c(0, 0))
+    }
+    return(metrics)
+}
diff --git a/R/PhenoKl.R b/R/PhenoKl.R
new file mode 100644
index 00000000..a74cb032
--- /dev/null
+++ b/R/PhenoKl.R
@@ -0,0 +1,138 @@
+
+#' Phenology extraction in inflexion method (Zhang)
+#'
+#' @inheritParams PhenoDeriv
+#' @param fFIT object return by [curvefit()]
+#'
+#' @inherit PhenoTrs examples 
+#' @references
+#' 1. Zhang, X., Friedl, M. A., Schaaf, C. B., Strahler, A. H., Hodges, J. C. F.
+#'    F., Gao, F., … Huete, A. (2003). Monitoring vegetation phenology using
+#'    MODIS. Remote Sensing of Environment, 84(3), 471–475.
+#'    \doi{10.1016/S0034-4257(02)00135-9}
+#' @return A numeric vector, with the elements of: 
+#' `Greenup`, `Maturity`, `Senescence`, `Dormancy`.
+#' 
+#' @export
+PhenoKl <- function(fFIT, t = NULL,
+                    analytical = TRUE, smoothed.spline = FALSE,
+                    IsPlot = TRUE, show.legend = TRUE, ...) {
+    
+    PhenoNames <- c("Greenup", "Maturity", "Senescence", "Dormancy")
+    metrics <- setNames(rep(NA, 4), PhenoNames)
+
+    if (!is.null(fFIT$tout)) t <- fFIT$tout
+    x <- last2(fFIT$zs)
+    n <- length(t)
+    xlim <- range(t)
+
+    # get peak of season position
+    half.season <- median(which.max(x)) # + 20, half season + 20 was unreasonable
+    pos <- t[half.season]
+
+    if (all(is.na(x))) {
+        return(metrics)
+    }
+    if (half.season < 5 || half.season > (n - 5)) {
+        return(metrics)
+    }
+
+    derivs <- curvature.fFIT(fFIT, t, analytical, smoothed.spline)
+    k <- derivs$k
+    # define cutoff date for spline functions
+
+    # x <- ifelse(uncert==TRUE, x$uncertainty, x$fit)
+    # if (length(which(is.na(k) == TRUE)) != 0 | length(which(is.infinite(k) == TRUE)) != 0) {
+
+    # If have no NA and infinite x
+    if (!(any(is.na(k)) || any(is.infinite(k)))) {
+        spline.k <- smooth.spline(k, df = 0.1 * length(k))
+        der.k <- predict(spline.k, d = 1)$y
+        der.k2 <- predict(spline.k, d = 2)$y
+
+        # der.k  <- c(NA, diff(k))
+        # der.k2 <- c(NA, NA, diff(k, differences = 2))
+        ## find maxima of derivative of k ## split season
+        dist_fromPeak <- 0 # days
+        asc.k <- try(der.k[1:(half.season - dist_fromPeak)]) # k   of first half year
+        asc.k.d <- try(t[1:(half.season - dist_fromPeak)])
+        # doy of first half year
+        des.k <- try(der.k[(half.season + dist_fromPeak):length(k)])
+        des.k.d <- try(t[(half.season + dist_fromPeak):length(k)])
+
+        A <- range(der.k, na.rm = TRUE) %>% diff()
+        minpeakheight <- A * 0.025
+
+        # first half maximum local x of k'
+        # minpeakdistance in the unit of day
+        pos <- findpeaks(asc.k,
+            npeaks = 2, ndowns = 0, sortstr = TRUE,
+            minpeakdistance = 15, minpeakheight = minpeakheight
+        )$X$pos
+        pos <- sort(pos)
+        pos <- c(rep(NA, 2 - length(pos)), pos) # at least two x
+        I_asc <- asc.k.d[pos]
+        if (all(is.na(pos))) {
+            I_asc <- pos
+        } else {
+            I_asc <- asc.k.d[pos]
+        }
+
+        # second half minimum local x of k'
+        pos <- findpeaks(-des.k,
+            npeaks = 2, sortstr = TRUE,
+            minpeakdistance = 15, minpeakheight = minpeakheight
+        )$X$pos
+        pos <- sort(pos)
+        pos <- c(pos, rep(NA, 2 - length(pos)))
+        if (all(is.na(pos))) {
+            I_dec <- pos
+        } else {
+            I_dec <- des.k.d[pos]
+        }
+        metrics <- c(I_asc, I_dec)
+    }
+
+    if (IsPlot) {
+        main <- ifelse(all(par("mar") == 0), "", "Zhang (Curvature Rate)")
+        A <- diff(range(x))
+
+        I_metrics <- match(metrics, t)
+        if (all(is.na(I_metrics))) {
+            ylons <- I_metrics
+        } else {
+            ylons <- x[I_metrics] + c(-1, -1, -1, 1) * 0.1 * A
+        }
+
+        xlons <- metrics + c(1, -1, 1, -1) * 5
+        xlons[xlons < min(t)] <- min(t)
+        xlons[xlons > max(t)] <- max(t)
+        # plotrix::twoord.plot(t, x, t, der.k,
+        #                      main = main, type= c("p", "b"), xlim = xlim,
+        #                      rcol = "black", lcol = "grey60", lpch = 20, rpch = 1,
+        #                      rytickpos = NULL)
+        PhenoPlot(t, x, main = main, ...)
+        if (show.legend) {
+            legend("topright", c("K'"), lty = c(3), col = c("black"), bty = "n") ## pch =c(20, 1),
+        }
+
+        pos <- t[half.season]
+        abline(v = metrics, col = colors, lwd = linewidth)
+        # abline(v = pos, col ="darkgreen", lty = 1, lwd = linewidth)
+        # abline(v = pos + 20, col ="darkgreen", lty = 2, lwd = linewidth)
+        PhenoNames2 <- c("Greenup", "Maturity", "Senescence", "Dormancy")
+        # PhenoNames2 <- c("G", "M", "S", "D")
+
+        I <- c(1, 3)
+        text(xlons[I], ylons[I], PhenoNames2[I], col = colors[I], adj = c(0, 0))
+        I <- c(2, 4)
+        text(xlons[I], ylons[I], PhenoNames2[I], col = colors[I], adj = c(1, 0))
+        # der.k last plot
+        op <- par(new = TRUE)
+        plot(t, der.k,
+            xlim = xlim, type = "l",
+            lty = 3, lwd = linewidth, col = "black", axes = TRUE
+        ) # cex = 1, pch = 20,
+    }
+    return(setNames(metrics, PhenoNames))
+}
diff --git a/R/PhenoTrs.R b/R/PhenoTrs.R
new file mode 100644
index 00000000..1f9e16ab
--- /dev/null
+++ b/R/PhenoTrs.R
@@ -0,0 +1,158 @@
+# colors    <- c("blue", "green3", "orange", "red")
+colors <- c("green3", "darkgreen", "darkorange3", "red")
+linewidth <- 1.2
+
+
+# identify greenup or dormancy(brown) period
+.Greenup <- function(x, ...) {
+    ratio.deriv <- c(NA, diff(x))
+    greenup <- rep(NA, length(x))
+    greenup[ratio.deriv > 0] <- TRUE
+    greenup[ratio.deriv < 0] <- FALSE
+    return(greenup)
+}
+
+#' Phenology extraction in Threshold method (TRS)
+#' 
+#' @param x numeric vector, or `fFIT` object returned by [curvefit()].
+#' @param t `doy` vector, corresponding doy of vegetation index.
+#' @param approach to be used to calculate phenology metrics.
+#' 'White' (White et al. 1997) or 'Trs' for simple threshold.
+#' @param trs threshold to be used for approach "Trs", in (0, 1).
+#' @param IsPlot whether to plot?
+#' @param ... other parameters to PhenoPlot
+#' @param asymmetric If true, background value in spring season and autumn season
+#' is regarded as different.
+#' 
+#' @example R/examples/ex-PhenoTrs.R
+#' @seealso [PhenoDeriv()], [PhenoGu()], [PhenoKl()]
+#' @export
+PhenoTrs <- function(x, t = NULL, 
+    approach = c("White", "Trs"), trs = 0.5, #, min.mean = 0.1
+    asymmetric = TRUE, IsPlot = TRUE, ...) 
+{
+    UseMethod("PhenoTrs", x)
+}
+
+#' @rdname PhenoTrs
+#' @export
+PhenoTrs.fFIT <- function(x, t = NULL, ...) {
+    if (!is.null(x$tout)) t <- x$tout
+    values <- last2(x$zs)
+    PhenoTrs.default(values, t, ...)
+}
+
+#' @export
+#' @rdname PhenoTrs
+PhenoTrs.default <- function(x, t = NULL,
+                             approach = c("White", "Trs"), trs = 0.5, # , min.mean = 0.1
+                             asymmetric = TRUE,
+                             IsPlot = TRUE, ...) {
+    metrics <- c(sos = NA, eos = NA)
+    n <- length(x)
+
+    # get peak of season position
+    half.season <- median(which.max(x)) %>% round() # + 20, half season + 20 was unreasonable
+    pos <- t[half.season]
+
+    if (all(is.na(x))) {
+        return(metrics)
+    }
+    if (half.season < 5 || half.season > (n - 5)) {
+        return(metrics)
+    }
+
+    # get statistical x
+    # avg  <- mean(x, na.rm = TRUE)
+
+    # avg2 <- mean(x2[x2 > min.mean], na.rm = TRUE)
+    peak <- max(x, na.rm = TRUE)
+
+    if (asymmetric) {
+        mn_a <- min(x[1:half.season], na.rm = T)
+        mn_b <- min(x[-(1:half.season)], na.rm = T)
+
+        mn <- c(
+            rep(mn_a, half.season),
+            rep(mn_b, n - half.season)
+        )
+    } else {
+        mn <- min(x, na.rm = TRUE)
+    }
+
+    ampl <- peak - mn
+
+    # select (or scale) x and thresholds for different methods
+    approach <- approach[1]
+    if (approach == "White") {
+        # scale annual time series to 0-1
+        ratio <- (x - mn) / ampl
+        # trs   <- 0.5
+        trs.low <- trs - 0.05
+        trs.up <- trs + 0.05
+    }
+    if (approach == "Trs") {
+        ratio <- x
+        a <- diff(range(ratio, na.rm = TRUE)) * 0.05
+        trs.low <- trs - a
+        trs.up <- trs + a
+    }
+
+    greenup <- .Greenup(ratio)
+    ## distinguish the first half year and second year
+    # select time where SOS and EOS are located (around trs value)
+    bool <- ratio >= trs.low & ratio <= trs.up
+
+    # get SOS, EOS, LOS
+    # fixed 2017-01-04, according to TP phenology property
+    # sos <- round(median(sose os[greenup & bool], na.rm = TRUE))
+    # eos <- round(median(soseos[!greenup & bool], na.rm = TRUE))
+
+    sos <- round(median(t[greenup & bool & t < pos], na.rm = TRUE))
+    eos <- round(median(t[!greenup & bool & t > pos], na.rm = TRUE))
+
+    # greenup <- .Greenup(ratio)
+    if (is.na(sos + eos)) {
+        bool2 <- ratio >= (trs - 0.1) & ratio <= (trs + 0.1)
+        if (is.na(sos)) sos <- round(median(t[greenup & bool2 & t < pos], na.rm = TRUE))
+        if (is.na(eos)) eos <- round(median(t[!greenup & bool2 & t > pos], na.rm = TRUE))
+
+        # if still is NA
+        if (is.na(sos)) sos <- round((t[1] + pos) / 2)
+        if (is.na(eos)) eos <- round((t[n] + pos) / 2)
+    }
+    # monotonous = length(unique(greenup[2:n]) %>% rm_empty()) == 1
+    # los <- eos - sos#los < 0 indicate that error
+    # los[los < 0] <- n + (eos[los < 0] - sos[los < 0])
+
+    # get MGS, MSP, MAU
+    # mgs <- mean(x[ratio > trs], na.rm = TRUE)
+    # msp <- mau <- NA
+    # if (!is.na(sos)) {
+    #     id <- (sos - 10):(sos + 10)
+    #     id <- id[(id > 0) & (id < n)]
+    #     msp <- mean(x[which(index(x) %in% id == TRUE)], na.rm = TRUE)
+    # }
+    # if (!is.na(eos)) {
+    #     id <- (eos - 10):(eos + 10)
+    #     id <- id[(id > 0) & (id < n)]
+    #     mau <- mean(x[which(index(x) %in% id == TRUE)], na.rm = TRUE)
+    # }
+    # metrics <- c(sos = sos, eos = eos, los = los, pos = pos, mgs = mgs,
+    #   rsp = NA, rau = NA, peak = peak, msp = msp, mau = mau)
+    metrics <- c(sos = sos, eos = eos) # , los = los
+
+    if (IsPlot) {
+        main <- ifelse(all(par("mar") == 0), "", sprintf("TRS%d", trs * 10))
+        PhenoPlot(t, x, main = main, ...)
+
+        lines(t, trs * ampl + mn, lwd = linewidth)
+        # lines(t, trs.low*ampl + mn, lty = 2, lwd = linewidth)
+        # lines(t, trs.up*ampl + mn, lty = 2, lwd = linewidth)
+        abline(v = metrics, col = colors[c(1, 4)], lwd = linewidth)
+        text(metrics[1] - 5, min(trs + 0.15, 1) * ampl[1] + mn[1], "SOS", col = colors[1], adj = c(1, 0))
+        text(metrics[2] + 5, min(trs + 0.15, 1) * last(ampl) + last(mn), "EOS", col = colors[4], adj = c(0, 0))
+    }
+    return(metrics)
+    ### The function returns a vector with SOS, EOS, LOS, POS, MGS, rsp, rau, PEAK, MSP and MAU. }
+}
diff --git a/R/curvefit.R b/R/curvefit.R
index da5c3883..f9299304 100644
--- a/R/curvefit.R
+++ b/R/curvefit.R
@@ -19,7 +19,6 @@ phenonames <- c('TRS2.SOS', 'TRS2.EOS', 'TRS5.SOS', 'TRS5.EOS', 'TRS6.SOS', 'TRS
 #' @note 'Klos' have too many parameters. It will be slow and not stable.
 #'
 #' @return fFITs S3 object, see [fFITs()] for details.
-#' 
 #' @seealso [fFITs()]
 #' 
 #' @examples
@@ -66,6 +65,7 @@ curvefit <- function(y, t = index(y), tout = t,
     structure(list(data = data.table(y, t), tout = tout, model = fFITs), class = 'fFITs')
 }
 
+finefit <- curvefit
 
 #' @rdname curvefit
 #' @export
diff --git a/R/examples/ex-PhenoTrs.R b/R/examples/ex-PhenoTrs.R
new file mode 100644
index 00000000..650a94bd
--- /dev/null
+++ b/R/examples/ex-PhenoTrs.R
@@ -0,0 +1,17 @@
+# simulate vegetation time-series
+fFUN = doubleLog.Beck
+par  = c( mn  = 0.1 , mx  = 0.7 , sos = 50 , rsp = 0.1 , eos = 250, rau = 0.1)
+
+t    <- seq(1, 365, 8)
+tout <- seq(1, 365, 1)
+y <- fFUN(par, t)
+
+methods <- c("AG", "Beck", "Elmore", "Gu", "Zhang") # "Klos" too slow
+fFITs <- curvefit(y, t, tout, methods)
+fFIT  <- fFITs$model$AG
+
+par(mfrow = c(2, 2))
+PhenoTrs(fFIT)
+PhenoDeriv(fFIT)
+PhenoGu(fFIT)
+PhenoKl(fFIT)
diff --git a/R/findpeaks.R b/R/findpeaks.R
index e7aae422..178ac0cb 100644
--- a/R/findpeaks.R
+++ b/R/findpeaks.R
@@ -30,8 +30,10 @@
 #' just the first npeaks are returned in the order of index.
 #' @param sortstr Boolean, Should the peaks be returned sorted in decreasing oreder of
 #' their maximum value?
-#' @param IsPlot Boolean.
-#'
+#' @param include_gregexpr Boolean (default `FALSE`), whether to include the 
+#' matched `gregexpr`?
+#' @param IsPlot Boolean, whether to plot?
+#' 
 #' @note
 #' In versions before v0.3.4, `findpeaks(c(1, 2, 3, 4, 4, 3, 1))` failed to detect
 #' peaks when a flat pattern exit in the middle.
diff --git a/R/get_pheno.R b/R/get_pheno.R
index 6ee57d46..c58544ca 100644
--- a/R/get_pheno.R
+++ b/R/get_pheno.R
@@ -18,54 +18,79 @@ PhenoPlot <- function(t, y, main = "", ...){
 #'
 #' @inheritParams get_GOF
 #' @inheritParams D
-#' @param fits A list of [fFITs()] object, for a single curve fitting
-#' method.
+#' @param x One of:
+#' - `rfit` (rought fitting object), returned by [brks2rfit()].
+#' - `fFITs` (fine fitting object), return by multiple curve fitting methods by [curvefit()] for 
+#'    a growing season.
+#' - list of [fFITs()] object, for multiple growing seasons.
 #' @param method Which fine curve fitting method to be extracted?
 #' @param TRS Threshold for `PhenoTrs`.
 #' @param IsPlot Boolean. Whether to plot figure?
-#' @param show_title Whether to show the name of fine curve fitting method
+#' @param show.title Whether to show the name of fine curve fitting method
 #' in top title?
 #' @param ... ignored.
 #'
-#' @note
-#' Please note that only a single fine curve fitting method allowed here!
-#'
 #' @return List of every year phenology metrics
 #'
 #' @example inst/examples/ex-get_fitting_param_GOF.R
 #' @export
-get_pheno <- function(fits, method,
+get_pheno <- function(x, ...) UseMethod("get_pheno", x)
+
+#' @inheritParams PhenoTrs
+#' @rdname get_pheno
+#' @export
+get_pheno.rfit <- function(x, TRS = c(0.2, 0.5), asymmetric = TRUE, ...) {
+    dt = x$dt
+    tout = x$tout
+    doys = as.numeric(difftime(tout, date.origin, units = "day"))
+    TRS %<>% set_names(., paste0("TRS", .*10))
+    pheno = dt %>% group_by(flag) %>% group_map(function(d, .y) {
+        ind = which(tout >= d$beg & tout <= d$end)
+        t = doys[ind]
+        values = last(x$zs)[ind]
+        der1 = diff(values) %>% c(NA, .)
+
+        p_TRS = map(TRS, ~ PhenoTrs.default(values, t, trs = ., IsPlot = FALSE, asymmetric = asymmetric))
+        p_DER = PhenoDeriv.default(values, t, der1, IsPlot = FALSE)
+        p_GU  = PhenoGu.default(values, t, der1, IsPlot = FALSE)
+        c(p_TRS, list(DER = p_DER, p_GU)) %>% unlist()
+    }) %>% set_names(dt$flag)
+    tidy_pheno(pheno)
+}
+
+#' @rdname get_pheno
+#' @export
+get_pheno.list <- function(x, method,
     TRS = c(0.2, 0.5, 0.6),
     analytical = TRUE, smoothed.spline = FALSE,
-    IsPlot = FALSE, show_title = TRUE, ...)
+    IsPlot = FALSE, show.title = TRUE, ...)
 {
-    if (!is.list(fits)) stop("Unsupported input type!")
-    if (class(fits) == 'fFITs')  fits <- list(fits)
+    if (!is.list(x)) stop("Unsupported input type!")
+    if (class(x) == 'fFITs')  x <- list(x)
 
-    names   <- names(fits) # methods
-    methods <- if (missing(method)) names(fits[[1]]$model) else method
+    names   <- names(x) # methods
+    methods <- if (missing(method)) names(x[[1]]$model) else method
 
     # pheno_list
     res <- lapply(set_names(seq_along(methods), methods), function(k){
         method <- methods[k]
         if (IsPlot){
-            op <- par(mfrow = c(length(fits), 5),
+            op <- par(mfrow = c(length(x), 5),
                 mgp = c(3, 0.6, 0), mar = rep(0, 4), yaxt = "n", xaxt = "n")
             if (isTRUE(all.equal(par("oma"), c(0, 0, 0, 0)))) {
                 margin_l = 5.5
-                oma <- if (show_title) c(1, margin_l, 4, 1) else c(1, margin_l, 2, 1)
+                oma <- if (show.title) c(1, margin_l, 4, 1) else c(1, margin_l, 2, 1)
                 par(oma = oma)
             }
         }
 
-        # fFITs
         pheno_list <- lapply(seq_along(names) %>% set_names(names), function(i){
-            fFITs <- fits[[i]]
+            fFITs <- x[[i]]
             .params = listk(
-                fFITs, method, TRS,
+                x = fFITs, method, TRS,
                 analytical, smoothed.spline, IsPlot,
-                showName_pheno = ifelse(i == 1, TRUE, FALSE),
-                title_left = names[i]
+                show.PhenoName = ifelse(i == 1, TRUE, FALSE),
+                title.left = names[i]
             )
             do.call(get_pheno.fFITs, .params)
         })
@@ -73,7 +98,7 @@ get_pheno <- function(fits, method,
         if (IsPlot){
             par(new = TRUE, mfrow = c(1, 1), oma = rep(0, 4), mar = rep(0, 4))
             plot(0, axes = F, type = "n", xaxt = "n", yaxt = "n") #
-            if (show_title) {
+            if (show.title) {
                 text(1, 1, method, font = 2, cex = 1.3)
             }
         }
@@ -82,19 +107,19 @@ get_pheno <- function(fits, method,
     lapply(res, tidy_pheno) %>% purrr::transpose()
 }
 
-#' @param fFITs `fFITs` object returned by [curvefit()]
-#' @param title_left String of growing season flag.
-#' @param showName_pheno Whether to show phenological methods names in the top panel?
+#' @param fFITs `fFITs` object returned by [curvefits()]
+#' @param title.left String of growing season flag.
+#' @param show.PhenoName Whether to show phenological methods names in the top panel?
 #'
 #' @rdname get_pheno
 #' @export
-get_pheno.fFITs <- function(fFITs, method,
+get_pheno.fFITs <- function(x, method,
     TRS = c(0.2, 0.5),
     analytical = TRUE, smoothed.spline = FALSE,
     IsPlot = FALSE,
-    title_left = "", showName_pheno = TRUE)
+    title.left = "", show.PhenoName = TRUE, ...)
 {
-    meths <- names(fFITs$model)
+    meths <- names(x$model)
     if (missing(method)) {
         method <- meths[1]
         warning(sprintf("method is missing and set to %s!", method))
@@ -103,7 +128,7 @@ get_pheno.fFITs <- function(fFITs, method,
         warning(sprintf("%s not in methods and set to %s!", method, meths[1]))
         method <- meths[1]
     }
-    model <- fFITs$model[[method]]
+    model <- x$model[[method]]
 
     # get_pheno methods
     methods  <- c(paste0("TRS", TRS*10),"DER","GU", "ZHANG")
@@ -112,10 +137,10 @@ get_pheno.fFITs <- function(fFITs, method,
     ypred  <- last2(model$zs)
     all_na <- all(is.na(ypred))
 
-    show.lgd = FALSE
+    show.legend = FALSE
     if (IsPlot && !all_na){
-        ti <- fFITs$data$t
-        yi <- fFITs$data$y # may have NA values.
+        ti <- x$data$t
+        yi <- x$data$y # may have NA values.
         # constrain plot ylims
         ylim0    <- c( pmin(min(yi, na.rm = TRUE), min(ypred)),
                        pmax(max(yi, na.rm = TRUE), max(ypred)))
@@ -123,10 +148,10 @@ get_pheno.fFITs <- function(fFITs, method,
         ylim     <- ylim0 + c(-1, 0.2) * 0.05 *A
         ylim_trs <- (ylim - ylim0) / A # TRS:0-1
 
-        PhenoPlot(fFITs$tout, ypred, ylim = ylim, yaxt = "s")
+        PhenoPlot(x$tout, ypred, ylim = ylim, yaxt = "s")
         lines(ti, yi, lwd = 1, col = "grey60")
 
-        QC_flag <- fFITs$data$QC_flag
+        QC_flag <- x$data$QC_flag
         # Different quality points with different color and shape.
         if (is.null(QC_flag)){
             points(ti, yi)
@@ -141,11 +166,11 @@ get_pheno.fFITs <- function(fFITs, method,
         }
 
         # show legend in first subplot
-        if (showName_pheno){
-            show.lgd <- TRUE
+        if (show.PhenoName){
+            show.legend <- TRUE
             legend('topright', c('y', "f(t)"), lty = c(1, 1), pch =c(1, NA), bty='n')
         }
-        stat <- get_GOF.fFITs(fFITs)
+        stat <- get_GOF.fFITs(x)
         stat <- subset(stat, meth == method) %>% select(-meth) %>% unlist() %>% round(3)
 
         exprs = list(
@@ -155,38 +180,37 @@ get_pheno.fFITs <- function(fFITs, method,
         )
 
         legend('topleft', do.call(expression, exprs), adj = c(0.2, 0.2), bty='n', text.col = "red")
-        # mtext(title_left, side = 2, line = 0.2)
-        mtext(title_left, side = 2, line = 1.8)
+        # mtext(title.left, side = 2, line = 0.2)
+        mtext(title.left, side = 2, line = 1.8)
     }
-    if (showName_pheno && IsPlot) mtext("Fine fitting", line = 0.2)
+    if (show.PhenoName && IsPlot) mtext("Fine fitting", line = 0.2)
 
     p_TRS <- lapply(TRS, function(trs) {
-        PhenoTrs(model, t = fFITs$tout, approach = "White", trs = trs, IsPlot = FALSE)
+        PhenoTrs(model, t = x$tout, approach = "White", trs = trs, IsPlot = FALSE)
     })
 
     if (IsPlot && !all_na) {
         p_TRS_last <- PhenoTrs(model, approach="White", trs=TRS_last, IsPlot = IsPlot, ylim = ylim)
-        if (showName_pheno) mtext(sprintf('TRS%d', TRS_last*10))
+        if (show.PhenoName) mtext(sprintf('TRS%d', TRS_last*10))
     }
 
-    param_common  <- list(model, t = fFITs$tout,
+    param_common  <- list(model, t = x$tout,
         analytical = analytical, smoothed.spline = smoothed.spline,
         IsPlot, ylim = ylim)
-    param_common2 <- c(param_common, list(show.lgd = show.lgd))
+    param_common2 <- c(param_common, list(show.legend = show.legend))
 
     der <- do.call(PhenoDeriv, param_common2)
-    if (showName_pheno && IsPlot) mtext("DER", line = 0.2)
+    if (show.PhenoName && IsPlot) mtext("DER", line = 0.2)
 
     zhang <- do.call(PhenoKl, param_common2)
-    if (showName_pheno && IsPlot) mtext("Inflexion ", line = 0.2)
+    if (show.PhenoName && IsPlot) mtext("Inflexion ", line = 0.2)
 
     gu <- do.call(PhenoGu, param_common)[1:4]
-    if (showName_pheno && IsPlot) mtext("Gu", line = 0.2)
+    if (show.PhenoName && IsPlot) mtext("Gu", line = 0.2)
 
     c(p_TRS, list(der, gu, zhang)) %>% set_names(methods)
 }
 
-
 #' tidy_pheno
 #'
 #' Tidy for every method with multiple years phenology data
@@ -213,17 +237,6 @@ get_pheno.fFITs <- function(fFITs, method,
 #' @export
 tidy_pheno <- function(pheno) {
     doy2date <- function(datenum) as.Date(unlist(datenum), origin = date.origin)
-    # phenonames <- c('TRS2.sos', 'TRS2.eos', 'TRS5.sos', 'TRS5.eos', 'TRS6.sos', 'TRS6.eos',
-    #                 'DER.sos', 'DER.pop', 'DER.eos',
-    #                 'GU.UD', 'GU.SD', 'GU.DD', 'GU.RD',
-    #                 'ZHANG.Greenup', 'ZHANG.Maturity', 'ZHANG.Senescence', 'ZHANG.Dormancy')
-
-    # fix the error of \code{pheno} without name
-    # names <- names(pheno)
-    # if (is.null(names)){
-    #     years <- seq(year(origin), by = 1, length.out = length(pheno))
-    #     names(pheno) <- years
-    # }
     names <- unlist(pheno[[1]]) %>% names()
 
     p_date <- map_df(pheno, function(x) {
@@ -242,6 +255,16 @@ tidy_pheno <- function(pheno) {
     list(doy = p_doy[, ..vars], date = p_date[, ..vars])
 }
 
+# phenonames <- c('TRS2.sos', 'TRS2.eos', 'TRS5.sos', 'TRS5.eos', 'TRS6.sos', 'TRS6.eos',
+#                 'DER.sos', 'DER.pop', 'DER.eos',
+#                 'GU.UD', 'GU.SD', 'GU.DD', 'GU.RD',
+#                 'ZHANG.Greenup', 'ZHANG.Maturity', 'ZHANG.Senescence', 'ZHANG.Dormancy')
+# fix the error of \code{pheno} without name
+# names <- names(pheno)
+# if (is.null(names)){
+#     years <- seq(year(origin), by = 1, length.out = length(pheno))
+#     names(pheno) <- years
+# }
 date2doy <- function(p_date){
     p_date %>% mutate(across(3:ncol(.), ~ as.numeric(.x - origin + 1)))
 }
diff --git a/R/global_options.R b/R/global_options.R
index 260a76e2..17f5c0f6 100644
--- a/R/global_options.R
+++ b/R/global_options.R
@@ -61,8 +61,8 @@ options_fitting <- list(
     south              = FALSE,
     ymin               = 0.1,
     wFUN               = "wTSM",
-    wmin               = 0.1, 
-    
+    wmin               = 0.1,
+
     ws                 = c(0.2, 0.5, 0.8), # initial weights
 
     # methods
@@ -73,7 +73,7 @@ options_fitting <- list(
     debug              = FALSE,
     # parameters
     season             = options_season,
-    fitting            = options_fitting, 
+    fitting            = options_fitting,
     initialized        = FALSE
 ))
 
@@ -82,26 +82,29 @@ options_fitting <- list(
 #' @param ... list of phenofit options
 #' FUN_season: character, `season_mov` or `season`
 #' rFUN: character, rough fitting function. `smooth_wWHIT`, `smooth_wSG` or `smooth_wHANTs`.
-#' @param options If not NULL, `options` will be used and `...` will be ignored.
+#' @param options If not NULL, `options` will be used and `...` patched.
 #' @examples
 #' set_options(verbose = FALSE)
 #' get_options("season") %>% str()
 #' @export
 set_options <- function(..., options = NULL) {
-    if (is.null(options)) options = list(...)
-    # rm unrelated parameters
+    if (is.null(options)) {
+        options = list(...)
+    } else {
+        options %<>% modifyList(list(...))
+    }
+    # find valid options, and rm unrelated parameters
     ind = match(names(options), names(.options)) %>% which.notna()
-    # This step might lead to error, but will improve performance
     if (.options$initialized && length(ind) == 0) return()
-    
+
     .options %<>% modifyList(options[ind])
     # `season` and `fitting` share the same parameter
-    pars_comm = c("wFUN", "wmin", "verbose") 
+    pars_comm = c("wFUN", "wmin", "verbose")
     for (par in pars_comm) {
         if (is.null(.options$fitting[[par]])) .options$fitting[[par]] <- .options[[par]]
         if (is.null(.options$season[[par]])) .options$season[[par]] <- .options[[par]]
     }
-    
+
     .options$fitting$wFUN %<>% check_function()
     .options$season$wFUN %<>% check_function()
     # .options$season$rFUN %<>% check_function()
diff --git a/R/phenofit-package.R b/R/phenofit-package.R
index 6ab33bf4..2ca89513 100644
--- a/R/phenofit-package.R
+++ b/R/phenofit-package.R
@@ -47,7 +47,7 @@ NULL
               "val", "type", "flag", "peak", # season, 
               "i", "qc", "y", "sitename",  # phenofit_TS.avhrr
               season_input, 
-              "years", "nyear", "d_fit", "info_peak"
+              "years", "nyear", "d_fit", "rfit", "info_peak"
             )
         )
     }
diff --git a/R/qcFUN.R b/R/qcFUN.R
index 78d6be7e..0155862f 100644
--- a/R/qcFUN.R
+++ b/R/qcFUN.R
@@ -241,6 +241,7 @@ qc_SPOT <- function (QA, wmin = 0.2, wmid = 0.5, wmax = 1) {
 #' | 10            | Cirrus                                | Good        | \eqn{w_{mid}}  |
 #' | 11            | Snow / Ice                            | Bad         | \eqn{w_{mid}}  |
 #' @references 
+#' https://developers.google.com/earth-engine/datasets/catalog/COPERNICUS_S2_SR
 #' 
 #' @examples
 #' qc_sentinel2(1:11)
diff --git a/R/roughFit.R b/R/roughFit.R
index 985bbf65..0d63288e 100644
--- a/R/roughFit.R
+++ b/R/roughFit.R
@@ -190,3 +190,32 @@ adjustRoughParam <- function(lambda, nf, frame,
     }
     listk(lambda, nf, frame, status = status)
 }
+
+#' get rough fitting
+#'
+#' @param brks returned by function [season_mov()]
+#'
+#' @keywords internal
+#' @return
+#' - `data`:
+#'     + t
+#'     + y
+#'     + QC_flag
+#' - `tout`:
+#' - `zs`: list of iter1, ..., itern
+#' - `ws`: list of iter1, ..., itern
+#' @export
+brks2rfit <- function(brks) {
+    dt = brks$dt
+    fit = brks$fit
+    data = fit[, .(t, y)]
+
+    # doys = difftime(data$t, data$t[1])
+    t = data$t #%>% as.integer()
+    tout = seq(data$t[1], last(data$t), by = "day") #%>% as.integer()
+    zs = ldply(fit %>% select(starts_with("ziter")),
+        ~approx2(t, .x, tout, na.rm = FALSE)$y)
+    structure(listk(data, dt, tout, zs), class = "rfit")
+}
+
+approx2 <- function(...) suppressWarnings(approx(...))
diff --git a/R/season_input.R b/R/season_input.R
index 2b9f6e2a..82d3d2f1 100644
--- a/R/season_input.R
+++ b/R/season_input.R
@@ -151,4 +151,3 @@ season_input <- function(INPUT, options = NULL, verbose = FALSE, ...)
     d_season = find_season.peaks(rfit, info_peak)
     listk(fit = rfit, dt = d_season)
 }
-
diff --git a/R/tools_plyr.R b/R/tools_plyr.R
new file mode 100644
index 00000000..72ee15c5
--- /dev/null
+++ b/R/tools_plyr.R
@@ -0,0 +1,11 @@
+# , .progress = "none", .inform = FALSE, .parallel = FALSE, .paropts = NULL
+#' @importFrom purrr as_mapper
+llply <- function(.data, .f = NULL, ...) {
+    if (is_empty(.data)) return(.data)
+    .f <- as_mapper(.f, ...)
+    lapply(.data, .f, ...)
+}
+
+ldply <- function(.data, .f = NULL, ...) {
+    llply(.data, .f, ...) %>% as.data.table()
+}
diff --git a/inst/examples/ex-curvefits.R b/inst/examples/ex-curvefits.R
index 36fcecaf..d236eca5 100644
--- a/inst/examples/ex-curvefits.R
+++ b/inst/examples/ex-curvefits.R
@@ -18,7 +18,7 @@ brks2 <- season_mov(INPUT,
     ))
 # plot_season(INPUT, brks2, d)
 # Fine fitting
-fit <- curvefits(
+fFITs <- curvefits(
     INPUT, brks2,
     options = list(
         methods = c("AG", "Beck", "Elmore", "Zhang"), #,"klos", "Gu"
@@ -27,10 +27,10 @@ fit <- curvefits(
     )
 )
 
-r_param = get_param(fit)
-r_pheno = get_pheno(fit)
-r_gof = get_GOF(fit)
-d_fit = get_fitting(fit)
+r_param = get_param(fFITs)
+r_pheno = get_pheno(fFITs)
+r_gof = get_GOF(fFITs)
+d_fit = get_fitting(fFITs)
 
 g <- plot_curvefits(d_fit, brks2)
 grid::grid.newpage(); grid::grid.draw(g)
diff --git a/man/PhenoDeriv.Rd b/man/PhenoDeriv.Rd
new file mode 100644
index 00000000..92f3bc75
--- /dev/null
+++ b/man/PhenoDeriv.Rd
@@ -0,0 +1,46 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/PhenoDeriv.R
+\name{PhenoDeriv}
+\alias{PhenoDeriv}
+\alias{PhenoDeriv.fFIT}
+\alias{PhenoDeriv.default}
+\title{Phenology extraction in Derivative method (DER)}
+\usage{
+PhenoDeriv(x, t, ...)
+
+\method{PhenoDeriv}{fFIT}(x, t = NULL, analytical = TRUE, smoothed.spline = FALSE, ...)
+
+\method{PhenoDeriv}{default}(x, t, der1, IsPlot = TRUE, show.legend = TRUE, ...)
+}
+\arguments{
+\item{x}{numeric vector, or \code{fFIT} object returned by \code{\link[=curvefit]{curvefit()}}.}
+
+\item{t}{\code{doy} vector, corresponding doy of vegetation index.}
+
+\item{...}{Other parameters will be ignored.}
+
+\item{analytical}{If true, \code{numDeriv} package \code{grad} and \code{hess}
+will be used; if false, \code{D1} and \code{D2} will be used.}
+
+\item{smoothed.spline}{Whether apply \code{smooth.spline} first?}
+
+\item{der1}{the first order difference}
+
+\item{IsPlot}{whether to plot?}
+
+\item{show.legend}{whether show figure lelend?}
+}
+\description{
+Phenology extraction in Derivative method (DER)
+}
+\references{
+\enumerate{
+\item Filippa, G., Cremonese, E., Migliavacca, M., Galvagno, M., Forkel, M.,
+Wingate, L., … Richardson, A. D. (2016). Phenopix: A R package for
+image-based vegetation phenology. Agricultural and Forest Meteorology,
+220, 141–150. \doi{10.1016/j.agrformet.2016.01.006}
+}
+}
+\seealso{
+\code{\link[=PhenoTrs]{PhenoTrs()}}, \code{\link[=PhenoGu]{PhenoGu()}}, \code{\link[=PhenoKl]{PhenoKl()}}
+}
diff --git a/man/PhenoExtractMeth.Rd b/man/PhenoExtractMeth.Rd
deleted file mode 100644
index 3d4baa3b..00000000
--- a/man/PhenoExtractMeth.Rd
+++ /dev/null
@@ -1,102 +0,0 @@
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/PhenoExtractMeth.R
-\name{PhenoExtractMeth}
-\alias{PhenoExtractMeth}
-\alias{PhenoTrs}
-\alias{PhenoDeriv}
-\alias{PhenoGu}
-\alias{PhenoKl}
-\title{Phenology Extraction methods}
-\usage{
-PhenoTrs(
-  fFIT,
-  t = NULL,
-  approach = c("White", "Trs"),
-  trs = 0.5,
-  asymmetric = TRUE,
-  IsPlot = TRUE,
-  ...
-)
-
-PhenoDeriv(
-  fFIT,
-  t = NULL,
-  analytical = TRUE,
-  smoothed.spline = FALSE,
-  IsPlot = TRUE,
-  show.lgd = TRUE,
-  ...
-)
-
-PhenoGu(
-  fFIT,
-  t = NULL,
-  analytical = TRUE,
-  smoothed.spline = FALSE,
-  IsPlot = TRUE,
-  ...
-)
-
-PhenoKl(
-  fFIT,
-  t = NULL,
-  analytical = TRUE,
-  smoothed.spline = FALSE,
-  IsPlot = TRUE,
-  show.lgd = TRUE,
-  ...
-)
-}
-\arguments{
-\item{fFIT}{\code{fFIT} object returned by \code{\link[=optim_pheno]{optim_pheno()}}.}
-
-\item{t}{\code{date} or \code{doy} vector, with the same length as \code{ypred}. This parameter
-is for the Julia version \code{curvefits}.}
-
-\item{approach}{to be used to calculate phenology metrics.
-'White' (White et al. 1997) or 'Trs' for simple threshold.}
-
-\item{trs}{threshold to be used for approach "Trs", in (0, 1).}
-
-\item{asymmetric}{If true, background value in spring season and autumn season
-is regarded as different.}
-
-\item{IsPlot}{whether to plot?}
-
-\item{...}{other parameters to PhenoPlot}
-
-\item{analytical}{If true, \code{numDeriv} package \code{grad} and \code{hess}
-will be used; if false, \code{D1} and \code{D2} will be used.}
-
-\item{smoothed.spline}{Whether apply \code{smooth.spline} first?}
-
-\item{show.lgd}{whether show figure lelend?}
-}
-\description{
-\itemize{
-\item \code{PhenoTrs} Threshold method
-\item \code{PhenoDeriv} Derivative method
-\item \code{PhenoGu} Gu method
-\item \code{PhenoKl} Inflection method
-}
-}
-\examples{
-library(phenofit)
-# simulate vegetation time-series
-fFUN = doubleLog.Beck
-par  = c( mn  = 0.1 , mx  = 0.7 , sos = 50 , rsp = 0.1 , eos = 250, rau = 0.1)
-
-t    <- seq(1, 365, 8)
-tout <- seq(1, 365, 1)
-y <- fFUN(par, t)
-
-methods <- c("AG", "Beck", "Elmore", "Gu", "Zhang") # "Klos" too slow
-fFITs <- curvefit(y, t, tout, methods)
-fFIT  <- fFITs$model$AG
-
-par(mfrow = c(2, 2))
-PhenoTrs(fFIT)
-PhenoDeriv(fFIT)
-PhenoGu(fFIT)
-PhenoKl(fFIT)
-}
diff --git a/man/PhenoGu.Rd b/man/PhenoGu.Rd
new file mode 100644
index 00000000..185e1074
--- /dev/null
+++ b/man/PhenoGu.Rd
@@ -0,0 +1,75 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/PhenoGu.R
+\name{PhenoGu}
+\alias{PhenoGu}
+\alias{PhenoGu.fFIT}
+\alias{PhenoGu.default}
+\title{Phenology extraction in GU method (GU)}
+\usage{
+PhenoGu(x, t, ...)
+
+\method{PhenoGu}{fFIT}(x, t = NULL, analytical = TRUE, smoothed.spline = FALSE, ...)
+
+\method{PhenoGu}{default}(x, t, der1, IsPlot = TRUE, ...)
+}
+\arguments{
+\item{x}{numeric vector, or \code{fFIT} object returned by \code{\link[=curvefit]{curvefit()}}.}
+
+\item{t}{\code{doy} vector, corresponding doy of vegetation index.}
+
+\item{...}{other parameters to \code{\link[=PhenoGu.default]{PhenoGu.default()}} or \code{\link[=PhenoGu.fFIT]{PhenoGu.fFIT()}}}
+
+\item{analytical}{If true, \code{numDeriv} package \code{grad} and \code{hess}
+will be used; if false, \code{D1} and \code{D2} will be used.}
+
+\item{smoothed.spline}{Whether apply \code{smooth.spline} first?}
+
+\item{der1}{the first order difference}
+
+\item{IsPlot}{whether to plot?}
+}
+\value{
+A numeric vector, with the elements of:
+\itemize{
+\item \code{UD}: upturn date
+\item \code{SD}: stabilisation date
+\item \code{DD}: downturn date
+\item \code{RD}: recession date
+}
+}
+\description{
+Phenology extraction in GU method (GU)
+}
+\examples{
+# simulate vegetation time-series
+fFUN = doubleLog.Beck
+par  = c( mn  = 0.1 , mx  = 0.7 , sos = 50 , rsp = 0.1 , eos = 250, rau = 0.1)
+
+t    <- seq(1, 365, 8)
+tout <- seq(1, 365, 1)
+y <- fFUN(par, t)
+
+methods <- c("AG", "Beck", "Elmore", "Gu", "Zhang") # "Klos" too slow
+fFITs <- curvefit(y, t, tout, methods)
+fFIT  <- fFITs$model$AG
+
+par(mfrow = c(2, 2))
+PhenoTrs(fFIT)
+PhenoDeriv(fFIT)
+PhenoGu(fFIT)
+PhenoKl(fFIT)
+}
+\references{
+\enumerate{
+\item Gu, L., Post, W. M., Baldocchi, D. D., Black, T. A., Suyker, A. E., Verma,
+S. B., … Wofsy, S. C. (2009). Characterizing the Seasonal Dynamics of
+Plant Community Photosynthesis Across a Range of Vegetation Types. In A.
+Noormets (Ed.), Phenology of Ecosystem Processes: Applications in Global
+Change Research (pp. 35–58). New York, NY: Springer New York.
+\doi{10.1007/978-1-4419-0026-5_2}
+\item Filippa, G., Cremonese, E., Migliavacca, M., Galvagno, M., Forkel, M.,
+Wingate, L., … Richardson, A. D. (2016). Phenopix: A R package for
+image-based vegetation phenology. Agricultural and Forest Meteorology,
+220, 141–150. \doi{10.1016/j.agrformet.2016.01.006}
+}
+}
diff --git a/man/PhenoKl.Rd b/man/PhenoKl.Rd
new file mode 100644
index 00000000..e72b4a11
--- /dev/null
+++ b/man/PhenoKl.Rd
@@ -0,0 +1,66 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/PhenoKl.R
+\name{PhenoKl}
+\alias{PhenoKl}
+\title{Phenology extraction in inflexion method (Zhang)}
+\usage{
+PhenoKl(
+  fFIT,
+  t = NULL,
+  analytical = TRUE,
+  smoothed.spline = FALSE,
+  IsPlot = TRUE,
+  show.legend = TRUE,
+  ...
+)
+}
+\arguments{
+\item{fFIT}{object return by \code{\link[=curvefit]{curvefit()}}}
+
+\item{t}{\code{doy} vector, corresponding doy of vegetation index.}
+
+\item{analytical}{If true, \code{numDeriv} package \code{grad} and \code{hess}
+will be used; if false, \code{D1} and \code{D2} will be used.}
+
+\item{smoothed.spline}{Whether apply \code{smooth.spline} first?}
+
+\item{IsPlot}{whether to plot?}
+
+\item{show.legend}{whether show figure lelend?}
+
+\item{...}{Other parameters will be ignored.}
+}
+\value{
+A numeric vector, with the elements of:
+\code{Greenup}, \code{Maturity}, \code{Senescence}, \code{Dormancy}.
+}
+\description{
+Phenology extraction in inflexion method (Zhang)
+}
+\examples{
+# simulate vegetation time-series
+fFUN = doubleLog.Beck
+par  = c( mn  = 0.1 , mx  = 0.7 , sos = 50 , rsp = 0.1 , eos = 250, rau = 0.1)
+
+t    <- seq(1, 365, 8)
+tout <- seq(1, 365, 1)
+y <- fFUN(par, t)
+
+methods <- c("AG", "Beck", "Elmore", "Gu", "Zhang") # "Klos" too slow
+fFITs <- curvefit(y, t, tout, methods)
+fFIT  <- fFITs$model$AG
+
+par(mfrow = c(2, 2))
+PhenoTrs(fFIT)
+PhenoDeriv(fFIT)
+PhenoGu(fFIT)
+PhenoKl(fFIT)
+}
+\references{
+\enumerate{
+\item Zhang, X., Friedl, M. A., Schaaf, C. B., Strahler, A. H., Hodges, J. C. F.
+F., Gao, F., … Huete, A. (2003). Monitoring vegetation phenology using
+MODIS. Remote Sensing of Environment, 84(3), 471–475.
+\doi{10.1016/S0034-4257(02)00135-9}
+}
+}
diff --git a/man/PhenoTrs.Rd b/man/PhenoTrs.Rd
new file mode 100644
index 00000000..018dbcb1
--- /dev/null
+++ b/man/PhenoTrs.Rd
@@ -0,0 +1,72 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/PhenoTrs.R
+\name{PhenoTrs}
+\alias{PhenoTrs}
+\alias{PhenoTrs.fFIT}
+\alias{PhenoTrs.default}
+\title{Phenology extraction in Threshold method (TRS)}
+\usage{
+PhenoTrs(
+  x,
+  t = NULL,
+  approach = c("White", "Trs"),
+  trs = 0.5,
+  asymmetric = TRUE,
+  IsPlot = TRUE,
+  ...
+)
+
+\method{PhenoTrs}{fFIT}(x, t = NULL, ...)
+
+\method{PhenoTrs}{default}(
+  x,
+  t = NULL,
+  approach = c("White", "Trs"),
+  trs = 0.5,
+  asymmetric = TRUE,
+  IsPlot = TRUE,
+  ...
+)
+}
+\arguments{
+\item{x}{numeric vector, or \code{fFIT} object returned by \code{\link[=curvefit]{curvefit()}}.}
+
+\item{t}{\code{doy} vector, corresponding doy of vegetation index.}
+
+\item{approach}{to be used to calculate phenology metrics.
+'White' (White et al. 1997) or 'Trs' for simple threshold.}
+
+\item{trs}{threshold to be used for approach "Trs", in (0, 1).}
+
+\item{asymmetric}{If true, background value in spring season and autumn season
+is regarded as different.}
+
+\item{IsPlot}{whether to plot?}
+
+\item{...}{other parameters to PhenoPlot}
+}
+\description{
+Phenology extraction in Threshold method (TRS)
+}
+\examples{
+# simulate vegetation time-series
+fFUN = doubleLog.Beck
+par  = c( mn  = 0.1 , mx  = 0.7 , sos = 50 , rsp = 0.1 , eos = 250, rau = 0.1)
+
+t    <- seq(1, 365, 8)
+tout <- seq(1, 365, 1)
+y <- fFUN(par, t)
+
+methods <- c("AG", "Beck", "Elmore", "Gu", "Zhang") # "Klos" too slow
+fFITs <- curvefit(y, t, tout, methods)
+fFIT  <- fFITs$model$AG
+
+par(mfrow = c(2, 2))
+PhenoTrs(fFIT)
+PhenoDeriv(fFIT)
+PhenoGu(fFIT)
+PhenoKl(fFIT)
+}
+\seealso{
+\code{\link[=PhenoDeriv]{PhenoDeriv()}}, \code{\link[=PhenoGu]{PhenoGu()}}, \code{\link[=PhenoKl]{PhenoKl()}}
+}
diff --git a/man/brks2rfit.Rd b/man/brks2rfit.Rd
new file mode 100644
index 00000000..e7be7311
--- /dev/null
+++ b/man/brks2rfit.Rd
@@ -0,0 +1,28 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/roughFit.R
+\name{brks2rfit}
+\alias{brks2rfit}
+\title{get rough fitting}
+\usage{
+brks2rfit(brks)
+}
+\arguments{
+\item{brks}{returned by function \code{\link[=season_mov]{season_mov()}}}
+}
+\value{
+\itemize{
+\item \code{data}:
+\itemize{
+\item t
+\item y
+\item QC_flag
+}
+\item \code{tout}:
+\item \code{zs}: list of iter1, ..., itern
+\item \code{ws}: list of iter1, ..., itern
+}
+}
+\description{
+get rough fitting
+}
+\keyword{internal}
diff --git a/man/curvefits.Rd b/man/curvefits.Rd
index deb57791..b06f6891 100644
--- a/man/curvefits.Rd
+++ b/man/curvefits.Rd
@@ -74,7 +74,7 @@ brks2 <- season_mov(INPUT,
     ))
 # plot_season(INPUT, brks2, d)
 # Fine fitting
-fit <- curvefits(
+fFITs <- curvefits(
     INPUT, brks2,
     options = list(
         methods = c("AG", "Beck", "Elmore", "Zhang"), #,"klos", "Gu"
@@ -83,10 +83,10 @@ fit <- curvefits(
     )
 )
 
-r_param = get_param(fit)
-r_pheno = get_pheno(fit)
-r_gof = get_GOF(fit)
-d_fit = get_fitting(fit)
+r_param = get_param(fFITs)
+r_pheno = get_pheno(fFITs)
+r_gof = get_GOF(fFITs)
+d_fit = get_fitting(fFITs)
 
 g <- plot_curvefits(d_fit, brks2)
 grid::grid.newpage(); grid::grid.draw(g)
diff --git a/man/findpeaks.Rd b/man/findpeaks.Rd
index d8283955..1e164df3 100644
--- a/man/findpeaks.Rd
+++ b/man/findpeaks.Rd
@@ -55,7 +55,10 @@ just the first npeaks are returned in the order of index.}
 \item{sortstr}{Boolean, Should the peaks be returned sorted in decreasing oreder of
 their maximum value?}
 
-\item{IsPlot}{Boolean.}
+\item{include_gregexpr}{Boolean (default \code{FALSE}), whether to include the
+matched \code{gregexpr}?}
+
+\item{IsPlot}{Boolean, whether to plot?}
 }
 \description{
 Find peaks (maxima) in a time series. This function is modified from
diff --git a/man/get_pheno.Rd b/man/get_pheno.Rd
index 9b1e81e0..706cfd22 100644
--- a/man/get_pheno.Rd
+++ b/man/get_pheno.Rd
@@ -2,39 +2,56 @@
 % Please edit documentation in R/get_pheno.R
 \name{get_pheno}
 \alias{get_pheno}
+\alias{get_pheno.rfit}
+\alias{get_pheno.list}
 \alias{get_pheno.fFITs}
 \title{get_pheno}
 \usage{
-get_pheno(
-  fits,
+get_pheno(x, ...)
+
+\method{get_pheno}{rfit}(x, TRS = c(0.2, 0.5), asymmetric = TRUE, ...)
+
+\method{get_pheno}{list}(
+  x,
   method,
   TRS = c(0.2, 0.5, 0.6),
   analytical = TRUE,
   smoothed.spline = FALSE,
   IsPlot = FALSE,
-  show_title = TRUE,
+  show.title = TRUE,
   ...
 )
 
-get_pheno.fFITs(
-  fFITs,
+\method{get_pheno}{fFITs}(
+  x,
   method,
   TRS = c(0.2, 0.5),
   analytical = TRUE,
   smoothed.spline = FALSE,
   IsPlot = FALSE,
-  title_left = "",
-  showName_pheno = TRUE
+  title.left = "",
+  show.PhenoName = TRUE,
+  ...
 )
 }
 \arguments{
-\item{fits}{A list of \code{\link[=fFITs]{fFITs()}} object, for a single curve fitting
-method.}
+\item{x}{One of:
+\itemize{
+\item \code{rfit} (rought fitting object), returned by \code{\link[=brks2rfit]{brks2rfit()}}.
+\item \code{fFITs} (fine fitting object), return by multiple curve fitting methods by \code{\link[=curvefit]{curvefit()}} for
+a growing season.
+\item list of \code{\link[=fFITs]{fFITs()}} object, for multiple growing seasons.
+}}
 
-\item{method}{Which fine curve fitting method to be extracted?}
+\item{...}{ignored.}
 
 \item{TRS}{Threshold for \code{PhenoTrs}.}
 
+\item{asymmetric}{If true, background value in spring season and autumn season
+is regarded as different.}
+
+\item{method}{Which fine curve fitting method to be extracted?}
+
 \item{analytical}{If true, \code{numDeriv} package \code{grad} and \code{hess}
 will be used; if false, \code{D1} and \code{D2} will be used.}
 
@@ -42,16 +59,14 @@ will be used; if false, \code{D1} and \code{D2} will be used.}
 
 \item{IsPlot}{Boolean. Whether to plot figure?}
 
-\item{show_title}{Whether to show the name of fine curve fitting method
+\item{show.title}{Whether to show the name of fine curve fitting method
 in top title?}
 
-\item{...}{ignored.}
+\item{title.left}{String of growing season flag.}
 
-\item{fFITs}{\code{fFITs} object returned by \code{\link[=curvefit]{curvefit()}}}
+\item{show.PhenoName}{Whether to show phenological methods names in the top panel?}
 
-\item{title_left}{String of growing season flag.}
-
-\item{showName_pheno}{Whether to show phenological methods names in the top panel?}
+\item{fFITs}{\code{fFITs} object returned by \code{\link[=curvefits]{curvefits()}}}
 }
 \value{
 List of every year phenology metrics
@@ -59,9 +74,6 @@ List of every year phenology metrics
 \description{
 Get yearly vegetation phenological metrics of a curve fitting method
 }
-\note{
-Please note that only a single fine curve fitting method allowed here!
-}
 \examples{
 library(phenofit)
 # simulate vegetation time-series
diff --git a/man/plot_curvefits.Rd b/man/plot_curvefits.Rd
index ae8efdb7..b224dc7e 100644
--- a/man/plot_curvefits.Rd
+++ b/man/plot_curvefits.Rd
@@ -77,7 +77,7 @@ brks2 <- season_mov(INPUT,
     ))
 # plot_season(INPUT, brks2, d)
 # Fine fitting
-fit <- curvefits(
+fFITs <- curvefits(
     INPUT, brks2,
     options = list(
         methods = c("AG", "Beck", "Elmore", "Zhang"), #,"klos", "Gu"
@@ -86,10 +86,10 @@ fit <- curvefits(
     )
 )
 
-r_param = get_param(fit)
-r_pheno = get_pheno(fit)
-r_gof = get_GOF(fit)
-d_fit = get_fitting(fit)
+r_param = get_param(fFITs)
+r_pheno = get_pheno(fFITs)
+r_gof = get_GOF(fFITs)
+d_fit = get_fitting(fFITs)
 
 g <- plot_curvefits(d_fit, brks2)
 grid::grid.newpage(); grid::grid.draw(g)
diff --git a/man/qc_sentinel2.Rd b/man/qc_sentinel2.Rd
index 6d9af99d..fe865bc6 100644
--- a/man/qc_sentinel2.Rd
+++ b/man/qc_sentinel2.Rd
@@ -35,5 +35,5 @@ qc_sentinel2(SCL, wmin = 0.2, wmid = 0.5, wmax = 1)
 qc_sentinel2(1:11)
 }
 \references{
-
+https://developers.google.com/earth-engine/datasets/catalog/COPERNICUS_S2_SR
 }
diff --git a/man/set_options.Rd b/man/set_options.Rd
index f3d31aa9..d6978f9e 100644
--- a/man/set_options.Rd
+++ b/man/set_options.Rd
@@ -14,7 +14,7 @@ get_options(names = NULL)
 FUN_season: character, \code{season_mov} or \code{season}
 rFUN: character, rough fitting function. \code{smooth_wWHIT}, \code{smooth_wSG} or \code{smooth_wHANTs}.}
 
-\item{options}{If not NULL, \code{options} will be used and \code{...} will be ignored.}
+\item{options}{If not NULL, \code{options} will be used and \code{...} patched.}
 
 \item{names}{vector of character, names of options}
 }
diff --git a/scripts b/scripts
index 36c39e14..d139ddf2 160000
--- a/scripts
+++ b/scripts
@@ -1 +1 @@
-Subproject commit 36c39e14f373a5e13cd8ed9fc36a7769455c3d6a
+Subproject commit d139ddf216996f9650129069fceed4a1e93bbfce
diff --git a/tests/testthat/test-PhenoExtract.R b/tests/testthat/test-PhenoExtract.R
index 81ce9d0a..f3e04c9a 100644
--- a/tests/testthat/test-PhenoExtract.R
+++ b/tests/testthat/test-PhenoExtract.R
@@ -5,11 +5,11 @@ wFUN = wTSM # wBisquare #
 
 # The `maxExtendMonth` in season_mov and curvefits is different
 # lambda   <- init_lambda(INPUT$y) # lambda for whittaker
-brks2 <- season_mov(INPUT,
+brks <- season_mov(INPUT,
     options = list(rFUN = "smooth_wWHIT", wFUN = wFUN)
 )
 param <- list(
-    INPUT, brks2,
+    INPUT, brks,
     options = list(
         methods = c("AG", "Beck", "Elmore", "Gu", "Zhang"), # ,"klos",
         wFUN = wFUN, nextend = 2, maxExtendMonth = 3, minExtendMonth = 1,
@@ -17,7 +17,7 @@ param <- list(
     )
 )
 
-fit  <- do.call(curvefits, param)
+fit  <- do.call(curvefits, param)[1:6]
 
 # test plot.fFIT
 expect_silent({
@@ -35,29 +35,33 @@ expect_silent({
         stat  <- get_GOF(fit)
         d_fit <- get_fitting(fit)
 
-        g <- plot_curvefits(d_fit, brks2)
+        g <- plot_curvefits(d_fit, brks)
         grid::grid.newpage(); grid::grid.draw(g)
     })
 })
 
-
-# analytical
-expect_silent({
-    p <- get_pheno(fit,
-        analytical = TRUE, smoothed.spline = FALSE,
-        IsPlot = TRUE)
-})
-
-# numDeriv
-expect_silent({
-    p <- get_pheno(fit,
-        analytical = FALSE, smoothed.spline = FALSE,
-        IsPlot = FALSE)
+test_that("get_pheno.fFIT works", {
+    # analytical
+    expect_silent({
+        p <- get_pheno(fit,
+            analytical = TRUE, smoothed.spline = FALSE, IsPlot = TRUE)
+    })
+    # numDeriv
+    expect_silent({
+        p <- get_pheno(fit,
+            analytical = FALSE, smoothed.spline = FALSE, IsPlot = FALSE)
+    })
+    # spline
+    expect_silent({
+        p <- get_pheno(fit,
+            analytical = FALSE, smoothed.spline = TRUE, IsPlot = FALSE)
+    })
 })
 
-# spline
-expect_silent({
-    p <- get_pheno(fit,
-        analytical = FALSE, smoothed.spline = TRUE,
-        IsPlot = FALSE)
+test_that("get_pheno.rfit works", {
+    rfit = brks2rfit(brks)
+    r = get_pheno(rfit)
+    expect_equal(names(r), c("doy", "date"))
+    expect_equal(nrow(r$doy), nrow(brks$dt))
+    expect_equal(colnames(r$doy)[1:2], c("flag", "origin"))    
 })
diff --git a/vignettes/phenofit-procedures.Rmd b/vignettes/phenofit-procedures.Rmd
index 9a0f4d5a..150ed197 100644
--- a/vignettes/phenofit-procedures.Rmd
+++ b/vignettes/phenofit-procedures.Rmd
@@ -124,7 +124,7 @@ grid::grid.draw(g)
 
 # extract phenology metrics, only the first 3 year showed at here
 # write_fig({
-l_pheno <- get_pheno(fit[1:7], method = "AG", TRS = TRS, IsPlot = TRUE, show_title = FALSE)
+l_pheno <- get_pheno(fit[1:7], method = "AG", TRS = TRS, IsPlot = TRUE, show.title = FALSE)
 # }, "Figure6_phenology_metrics.pdf", 8, 6, show = TRUE)
 ```