diff --git a/.gitmodules b/.gitmodules
index ed0482f1..96217218 100644
--- a/.gitmodules
+++ b/.gitmodules
@@ -1,3 +1,3 @@
 [submodule "scripts"]
 	path = scripts
-	url = https://github.com/kongdd/phenofit-scripts
+	url = https://github.com/eco-hydro/phenofit-scripts
diff --git a/NAMESPACE b/NAMESPACE
index 0b098b65..a57f58c1 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -74,6 +74,7 @@ export(lambda_cv_jl)
 export(lambda_vcurve)
 export(lambda_vcurve_jl)
 export(logistic)
+export(make_date)
 export(melt_list)
 export(movmean)
 export(opt_nlm)
@@ -164,6 +165,7 @@ importFrom(lubridate,ddays)
 importFrom(lubridate,dyears)
 importFrom(lubridate,is.Date)
 importFrom(lubridate,leap_year)
+importFrom(lubridate,make_date)
 importFrom(lubridate,month)
 importFrom(lubridate,yday)
 importFrom(lubridate,year)
diff --git a/R/PhenoExtractMeth.R b/R/PhenoExtractMeth.R
index caa6183d..03d13a13 100644
--- a/R/PhenoExtractMeth.R
+++ b/R/PhenoExtractMeth.R
@@ -1,3 +1,7 @@
+# colors    <- c("blue", "green3", "orange", "red")
+colors <- c("green3", "darkgreen", "darkorange3", "red")
+linewidth <- 1.2
+
 #' @name PhenoExtractMeth
 #' @title Phenology Extraction methods
 #'
@@ -16,7 +20,7 @@
 #' # 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)
@@ -57,7 +61,7 @@ PhenoTrs <- function(fFIT, approach = c("White", "Trs"), trs = 0.5, #, min.mean
 
     # get peak of season position
     half.season <- median(which.max(values)) %>% round() # + 20, half season + 20 was unreasonable
-    pop <- t[half.season]
+    pos <- t[half.season]
 
     if (all(is.na(values))) return(metrics)
     if (half.season < 5 || half.season > (n - 5)) return(metrics)
@@ -107,8 +111,8 @@ PhenoTrs <- function(fFIT, approach = c("White", "Trs"), trs = 0.5, #, min.mean
     # 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 < pop], na.rm = TRUE))
-    eos <- round(median(t[!greenup & bool & t > pop], 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))
     # los <- eos - sos#los < 0 indicate that error
     # los[los < 0] <- n + (eos[los < 0] - sos[los < 0])
 
@@ -125,7 +129,7 @@ PhenoTrs <- function(fFIT, approach = c("White", "Trs"), trs = 0.5, #, min.mean
     #     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, pop = pop, mgs = mgs,
+    # 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
 
@@ -134,15 +138,14 @@ PhenoTrs <- function(fFIT, approach = c("White", "Trs"), trs = 0.5, #, min.mean
         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)
-
+        # 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, POP, MGS, rsp, rau, PEAK, MSP and MAU. }
+    ### The function returns a vector with SOS, EOS, LOS, POS, MGS, rsp, rau, PEAK, MSP and MAU. }
 }
 
 # identify greenup or dormancy(brown) period
@@ -164,16 +167,16 @@ PhenoDeriv <- function(fFIT,
     analytical = TRUE, smoothed.spline = FALSE,
     IsPlot = TRUE, show.lgd = TRUE, ...)
 {
-    PhenoNames <- c("SOS", "POP", "EOS")
-    metrics <- setNames(rep(NA, 3), c("sos", "pop", "eos")) # template
+    PhenoNames <- c("SOS", "POS", "EOS")
+    metrics <- setNames(rep(NA, 3), c("sos", "pos", "eos")) # template
 
     t      <- fFIT$tout
     values <- last(fFIT$zs)
     n      <- length(t)
 
     # get peak of season position
-    half.season <- median(which.max(values)) # deal with multiple pop values
-    pop <- t[half.season]
+    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)
@@ -197,7 +200,7 @@ PhenoDeriv <- function(fFIT,
     if (is_empty(sos)) sos <- NA
     if (is_empty(eos)) eos <- NA
 
-    metrics <- c(sos = sos, pop = pop, eos = eos)#, los = los
+    metrics <- c(sos = sos, pos = pos, eos = eos)#, los = los
 
     if (IsPlot){
         main  <- ifelse(all(par("mar") == 0), "", "DER")
@@ -205,21 +208,21 @@ PhenoDeriv <- function(fFIT,
         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 = pop, col ="darkgreen", lty = 1, 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, -1, 1)*0.1*A
+            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("darkgreen", colors[3]), adj = c(0, 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)
@@ -245,8 +248,8 @@ PhenoGu <- function(fFIT,
     n      <- length(t)
 
     # get peak of season position
-    half.season <- median(which.max(values)) # deal with multiple pop values
-    pop <- t[half.season]
+    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)
@@ -329,6 +332,7 @@ PhenoGu <- function(fFIT,
         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))
     }
@@ -352,7 +356,7 @@ PhenoKl <- function(fFIT,
 
     # get peak of season position
     half.season <- median(which.max(values)) # + 20, half season + 20 was unreasonable
-    pop <- t[half.season]
+    pos <- t[half.season]
 
     if (all(is.na(values))) return(metrics)
     if (half.season < 5 || half.season > (n - 5)) return(metrics)
@@ -429,12 +433,15 @@ PhenoKl <- function(fFIT,
             legend('topright', c("K'"), lty = c(3), col = c("black"), bty='n') ##pch =c(20, 1),
         }
 
-        pop     <- t[half.season]
+        pos     <- t[half.season]
         abline(v = metrics, col = colors, lwd = linewidth)
-        # abline(v = pop, col ="darkgreen", lty = 1, lwd = linewidth)
-        # abline(v = pop + 20, col ="darkgreen", lty = 2, lwd = linewidth)
-        I <- c(1, 3); text(xlons[I], ylons[I], PhenoNames[I], col = colors[I], adj = c(0, 0))
-        I <- c(2, 4); text(xlons[I], ylons[I], PhenoNames[I], col = colors[I], adj = c(1, 0))
+        # 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",
diff --git a/R/S3_fFITs.R b/R/S3_fFITs.R
index 5e79f2eb..25f59036 100644
--- a/R/S3_fFITs.R
+++ b/R/S3_fFITs.R
@@ -72,7 +72,7 @@ plot.fFITs <- function(x, method, ...){
 
         pred   <- last(fFIT$zs)
 
-        pop    <- t[which.max(pred)]
+        pos    <- t[which.max(pred)]
         derivs <- curvature(fFIT)
 
         # plot for every method
@@ -84,7 +84,7 @@ plot.fFITs <- function(x, method, ...){
             type= "b", pch = 20, cex = 1.3, col = "grey",
             main = "curve fitting VI", xlab = "Index", ylab = "VI")
         lines(t, pred); grid()
-        abline(v = pop, col ="green")
+        abline(v = pos, col ="green")
 
         maxd_der1 <- t[which.max(derivs$der1)]
         mind_der1 <- t[which.min(derivs$der1)]
@@ -92,14 +92,14 @@ plot.fFITs <- function(x, method, ...){
         plot(t, derivs$der1, main = "D1"); grid()
         abline(v = maxd_der1, col ="blue")
         abline(v = mind_der1, col ="red")
-        abline(v = pop, col ="green")
+        abline(v = pos, col ="green")
 
         plot(t, derivs$der2, main = "D2"); grid()
         plot(t, derivs$k, main = "k")    ; grid()
         abline(v = maxd_der1, col ="blue")
         abline(v = mind_der1, col ="red")
-        abline(v = pop, col ="green")
-        abline(v = pop + 20, col ="green", lty = 2)
+        abline(v = pos, col ="green")
+        abline(v = pos + 20, col ="green", lty = 2)
         # plot(diff(der1_diff), main = "diff2")
 
         # k <- derivs$k
diff --git a/R/curvefit.R b/R/curvefit.R
index 5ac5f1f9..6c8f7be5 100644
--- a/R/curvefit.R
+++ b/R/curvefit.R
@@ -1,5 +1,5 @@
 phenonames <- c('TRS2.SOS', 'TRS2.EOS', 'TRS5.SOS', 'TRS5.EOS', 'TRS6.SOS', 'TRS6.EOS',
-    'DER.SOS', 'DER.POP', 'DER.EOS',
+    'DER.SOS', 'DER.POS', 'DER.EOS',
     'UD', 'SD', 'DD','RD',
     'GreenUp', 'Maturity', 'Senescence', 'Dormancy')
 
diff --git a/R/findpeaks.R b/R/findpeaks.R
index 4c33811b..de0485cf 100644
--- a/R/findpeaks.R
+++ b/R/findpeaks.R
@@ -176,7 +176,7 @@ findpeaks_season_jl <- function(
     nyear = 1)
 {
     A = max(ypred) - min(ypred)
-    ans = JuliaCall::julia_call("findpeaks_season", ypred,
+    ans = JuliaCall::julia_call("phenofit.findpeaks_season", ypred,
         r_max = r_max, r_min = r_min,
         # r_max = y_max/A, r_min = y_min/A,
         minpeakdistance = as.integer(minpeakdistance), minpeakheight = minpeakheight,
@@ -184,4 +184,3 @@ findpeaks_season_jl <- function(
     ans$threshold <- data.table(y_max = r_max*A, y_min = r_min*A, r_max, r_min)
     ans
 }
-
diff --git a/R/get_pheno.R b/R/get_pheno.R
index c944ba4d..dcca2ed9 100644
--- a/R/get_pheno.R
+++ b/R/get_pheno.R
@@ -1,6 +1,3 @@
-colors    <- c("blue", "green3", "orange", "red")
-linewidth <- 1.2
-
 # ' PhenoPlot
 # '
 # ' @inheritParams check_input
@@ -11,8 +8,8 @@ linewidth <- 1.2
 PhenoPlot <- function(t, y, main = "", ...){
     plot(t, y, main = main, ...,
              type= "l", cex = 2, col = "black", lwd = linewidth) #pch = 20,
-    # grid(nx = NA)
-    grid(ny = 4, nx = NA)
+    grid(nx = NA)
+    # grid(ny = 4, nx = NA)
 }
 
 #' get_pheno
@@ -52,9 +49,13 @@ get_pheno <- function(fits, method,
     res <- lapply(set_names(seq_along(methods), methods), function(k){
         method <- methods[k]
         if (IsPlot){
-            oma <- if (show_title) c(1, 2, 4, 1) else c(1, 2, 2, 1)
-            op <- par(mfrow = c(length(fits), 5), oma = oma,
-                mar = rep(0, 4), yaxt = "n", xaxt = "n")
+            op <- par(mfrow = c(length(fits), 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)
+                par(oma = oma)
+            }
         }
 
         # fFITs
@@ -121,7 +122,7 @@ 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)
+        PhenoPlot(fFITs$tout, ypred, ylim = ylim, yaxt = "s")
         lines(ti, yi, lwd = 1, col = "grey60")
 
         QC_flag <- fFITs$data$QC_flag
@@ -153,9 +154,10 @@ 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 = 0.2)
+        mtext(title_left, side = 2, line = 1.8)
     }
-    if (showName_pheno && IsPlot) mtext("Fitting", line = 0.2)
+    if (showName_pheno && IsPlot) mtext("Fine fitting", line = 0.2)
 
     p_TRS <- lapply(TRS, function(trs) {
         PhenoTrs(fFIT, approach = "White", trs = trs, IsPlot = FALSE)
@@ -171,9 +173,14 @@ get_pheno.fFITs <- function(fFITs, method,
         IsPlot, ylim = ylim)
     param_common2 <- c(param_common, list(show.lgd = show.lgd))
 
-    der   <- do.call(PhenoDeriv, param_common2);  if (showName_pheno && IsPlot) mtext("DER", line = 0.2)
-    gu    <- do.call(PhenoGu, param_common)[1:4]; if (showName_pheno && IsPlot) mtext("GU", line = 0.2)
-    zhang <- do.call(PhenoKl, param_common2);     if (showName_pheno && IsPlot) mtext("ZHANG", line = 0.2)
+    der <- do.call(PhenoDeriv, param_common2)
+    if (showName_pheno && IsPlot) mtext("DER", line = 0.2)
+
+    zhang <- do.call(PhenoKl, param_common2)
+    if (showName_pheno && IsPlot) mtext("Inflexion ", line = 0.2)
+
+    gu <- do.call(PhenoGu, param_common)[1:4]
+    if (showName_pheno && IsPlot) mtext("Gu", line = 0.2)
 
     c(p_TRS, list(der, gu, zhang)) %>% set_names(methods)
 }
diff --git a/R/plot_curvefits.R b/R/plot_curvefits.R
index 5e35ec99..4fc643fe 100644
--- a/R/plot_curvefits.R
+++ b/R/plot_curvefits.R
@@ -1,3 +1,5 @@
+# ' @param theme ggplot theme to be applied
+
 #' plot_curvefits
 #'
 #' @param d_fit data.frame of curve fittings returned by [get_fitting()].
@@ -10,11 +12,12 @@
 #' @param yticks ticks of y axis
 #' @param font.size Font size of axis.text
 #' @param show.legend Boolean
-#' @param theme ggplot theme to be applied
 #' @param shape the shape of input VI observation? `line` or `point`
 #' @param cex point size for VI observation.
 #' @param angle `text.x` angle
-#'
+#' @param layer_extra (not used) extra ggplot layers
+#' @param ... ignored
+#' 
 #' @example inst/examples/ex-curvefits.R
 #' 
 #' @export
@@ -29,7 +32,9 @@ plot_curvefits <- function(
     theme = NULL,
     cex = 2,
     shape = "point", angle = 30,
-    show.legend = TRUE)
+    show.legend = TRUE, 
+    layer_extra = NULL, 
+    ...)
 {
     methods <- d_fit$meth %>% unique() %>% rm_empty() # in case of NA
     nmethod <- length(methods) # how many curve fitting methods?
@@ -69,8 +74,7 @@ plot_curvefits <- function(
         p <- p + geom_point(
             data = d_obs,
             aes_string("t", "y", shape = "QC_flag", color = "QC_flag", fill = "QC_flag"),
-            size = cex, alpha = 0.7
-        )
+            size = cex, alpha = 0.7)
     } else {
         p <- if (shape == "point") {
             p + geom_point(data = d_obs, aes_string("t", "y"), size = cex, alpha = 0.6, color = "grey60")
@@ -92,13 +96,13 @@ plot_curvefits <- function(
         scale_shape_manual(values = qc_shapes, drop = F) +
         coord_cartesian(xlim = xlim)
 
-    if (!is.null(theme)) p <- p + theme
     if (!is.null(yticks)) p <- p + scale_y_continuous(breaks = yticks)
     if  (is.null(title)) p <- p + theme(plot.title = element_blank())
-
+    if (!is.null(layer_extra)) p <- p + layer_extra
+    
     if (show.legend) {
-        iters_name_fine = c("Rough fitting", "Fine fitting")
-        lines_colors = c("black", "red")
+        iters_name_fine = c("Rough fitting", "iter1", "iter2")
+        lines_colors = c("black", "blue", "red")
         lgd <- make_legend_nmax(iters_name_fine, lines_colors, d_obs$QC_flag)
         p <- p + theme(legend.position = "none")
         p <- arrangeGrob(p, lgd,
diff --git a/R/plot_input.R b/R/plot_input.R
index 00df13cd..18370cda 100644
--- a/R/plot_input.R
+++ b/R/plot_input.R
@@ -13,7 +13,7 @@
 #' IsPlot = TRUE
 #' nptperyear = 23
 #' ypeak_min  = 0.05
-#' 
+#'
 #' INPUT    <- check_input(d$t, d$y, d$w, d$QC_flag, nptperyear,
 #'                         maxgap = nptperyear/4, alpha = 0.02, wmin = 0.2)
 #' plot_input(INPUT)
@@ -65,7 +65,11 @@ plot_input <- function(INPUT, wmin = 0.2, show.y0 = TRUE, ylab = "VI", ...){
         plot(t[I], y[I], type = "l", xaxt = "n", ylab = ylab, xlab = "Time", main = main, ...) # , ann = FALSE
         axis.Date(1, at=seq(min(t), max(t), by="month"), format="%Y-%m")
     } else {
-        plot(t[I], y[I], type = "l", ylab = ylab, xlab = "Time", main = main, ...)
+        plot(t[I], y[I], type = "l", ylab = ylab, xlab = "Time", main = main, xlim = xlim, ...)
+        year_min = min(years)
+        year_max = max(years)
+        axis.Date(1, at=seq(make_date(year_min), make_date(year_max), by="year"),
+                  format="%Y")
     }
 
     show.goodPoints = INPUT$nptperyear < 90
@@ -87,7 +91,7 @@ plot_input <- function(INPUT, wmin = 0.2, show.y0 = TRUE, ylab = "VI", ...){
         t_grids  <- seq.Date(date_beg, date_end, by = "year")
         abline(v = t_grids, col = "grey60", lty = 3)
     }
-    grid(nx = NA, NULL)
+    # grid(nx = NA, NULL)
     ylu <- INPUT$ylu
     if (!is.null(ylu)) abline(h=ylu, col="red", lty=2) # show ylims
 }
diff --git a/R/plot_season.R b/R/plot_season.R
index 9ae10d89..d5d3e2d5 100644
--- a/R/plot_season.R
+++ b/R/plot_season.R
@@ -14,18 +14,22 @@
 #' Only if `IsPlot=TRUE`, [plot_input()] will be used to plot.
 #' Known that y and w in `INPUT` have been changed, we suggest using the
 #' original data.table.
-
+#' 
+#' @param show.shade Boolean, period inside growing cycle colored as shade?
+#' @param margin `ylim = c(ymin, ymax + margin * A); A = ymax - ymin`.
+#' 
 #' @importFrom grid viewport pushViewport grid.draw
 #' @export
 plot_season <- function(
     INPUT, brks, plotdat, 
     # ylu,
     IsPlot.OnlyBad = FALSE, show.legend = TRUE,
-    ylab = "VI", title = NULL)
+    ylab = "VI", title = NULL, 
+    show.shade = TRUE,
+    margin = 0.35)
 {
-    if (missing(plotdat)) {
-        plotdat <- INPUT
-    }
+    if (missing(plotdat)) plotdat <- INPUT
+    
     if (is.data.frame(brks$dt[[1]])) {
         brks$dt %<>% do.call(rbind, .)
     }
@@ -64,10 +68,10 @@ plot_season <- function(
     ymin = min(INPUT$y0, na.rm = TRUE)
     ymax = max(INPUT$y0, na.rm = TRUE)
     # ylim = c(ymin, (ymax - ymin)*0.12 + ymax)
-    ylim = c(ymin, ylu[2] + 0.35 * A)
+    ylim = c(ymin, ylu[2] + margin * A)
 
     plot_input(plotdat, ylab = ylab, ylim = ylim)
-    plot_season_boundary(dt)
+    if (show.shade) plot_season_boundary(dt)
     
     NITER  <- ncol(zs)
     lines_colors <- iter_colors(NITER)
diff --git a/R/process_phenofit.R b/R/process_phenofit.R
index 5e48ceed..d1fc68f7 100644
--- a/R/process_phenofit.R
+++ b/R/process_phenofit.R
@@ -5,7 +5,8 @@
 #' @inheritParams curvefits
 #' @inheritParams season
 #' @param ... other parameters to [curvefits()]
-#' 
+#'
+#' @keywords internal
 #' @export
 process_phenofit <- function(
     d_obs,
@@ -14,7 +15,7 @@ process_phenofit <- function(
     .v_curve = FALSE,
     options_season = list(
         # rFUN = "smooth_wWHIT",
-        wFUN = "wTSM", 
+        wFUN = "wTSM",
         # wmin = 0.1,
         # iters = 2,
         # .lambda_vcurve = TRUE, lambda = NULL,
@@ -33,14 +34,17 @@ process_phenofit <- function(
         use.y0 = FALSE
     ),
     brks = NULL,
-    # TRS = c(0.1, 0.2, 0.5, 0.6, 0.8, 0.9),
+    TRS = c(0.1, 0.2, 0.5, 0.6, 0.8, 0.9),
     # ymin = 0.1, used for check_input
     # wsnow = 0.8,
     # use.y0 = FALSE,
     # overwrite = FALSE,
     run.curvefit = TRUE,
-    ...) 
+    ...)
 {
+    options_season %<>% modifyList(list(...))
+    options_fitting %<>% modifyList(list(...))
+
     ## 2.1 load site data
     # d_obs <- listk(t, y, w, QC_flag) %>% as.data.table()
     if (!("QC_flag" %in% colnames(d_obs))) {
@@ -60,7 +64,6 @@ process_phenofit <- function(
         r <- v_curve(INPUT, lg_lambdas, d = 2, IsPlot = FALSE)
         lambda <- r$lambda
     }
-    print(lambda)
     # wFUN <- "wBisquare", "wTSM", threshold_max = 0.1, IGBP = CSH
 
     brks2 <- season_mov(INPUT, options_season, ...)
diff --git a/R/process_season.R b/R/process_season.R
index 05a4149e..479d5632 100644
--- a/R/process_season.R
+++ b/R/process_season.R
@@ -7,8 +7,8 @@
 #' parameter, lambda.
 #'
 #' @note site-year may be not continuous.
-#'
-#' @rdname season_mov
+#' 
+#' @keywords internal
 #' @export
 process_season <- function(
     d_obs,
diff --git a/R/tools.R b/R/tools.R
index 3cdcb12e..6a8d25e2 100644
--- a/R/tools.R
+++ b/R/tools.R
@@ -117,3 +117,7 @@ magrittr::`%>%`
 
 #' @export
 magrittr::`%<>%`
+
+#' @importFrom lubridate make_date
+#' @export
+lubridate::make_date
diff --git a/data-raw/dat_CA-NS6.R b/data-raw/dat_CA-NS6.R
index efff5efe..c416f973 100644
--- a/data-raw/dat_CA-NS6.R
+++ b/data-raw/dat_CA-NS6.R
@@ -4,8 +4,11 @@ data("MOD13A1")
 df <- tidy_MOD13(MOD13A1$dt)
 st <- MOD13A1$st
 
-date_start <- as.Date("2013-01-01")
+date_start <- as.Date("2010-01-01")
 date_end <- as.Date("2016-12-31")
 
 sitename <- "CA-NS6" # df$site[1]
 d <- df[site == sitename & (date >= date_start & date <= date_end), ]
+
+CA_NS6 = d
+use_data(CA_NS6, overwrite = TRUE)
diff --git a/data/CA_NS6.rda b/data/CA_NS6.rda
index ac47854d..573c595f 100644
Binary files a/data/CA_NS6.rda and b/data/CA_NS6.rda differ
diff --git a/man/CA_NS6.Rd b/man/CA_NS6.Rd
index f3efe9b0..9fade5e8 100644
--- a/man/CA_NS6.Rd
+++ b/man/CA_NS6.Rd
@@ -5,7 +5,7 @@
 \alias{CA_NS6}
 \title{MOD13A1 EVI observations at flux site CA-NS6}
 \format{
-An object of class \code{data.table} (inherits from \code{data.frame}) with 92 rows and 6 columns.
+An object of class \code{data.table} (inherits from \code{data.frame}) with 161 rows and 6 columns.
 }
 \usage{
 data('CA_NS6')
diff --git a/man/plot_curvefits.Rd b/man/plot_curvefits.Rd
index e4861e43..7435560d 100644
--- a/man/plot_curvefits.Rd
+++ b/man/plot_curvefits.Rd
@@ -17,7 +17,9 @@ plot_curvefits(
   cex = 2,
   shape = "point",
   angle = 30,
-  show.legend = TRUE
+  show.legend = TRUE,
+  layer_extra = NULL,
+  ...
 )
 }
 \arguments{
@@ -37,8 +39,6 @@ of \code{t}, \code{y} and \code{QC_flag}. If not specified, it will be determine
 
 \item{font.size}{Font size of axis.text}
 
-\item{theme}{ggplot theme to be applied}
-
 \item{cex}{point size for VI observation.}
 
 \item{shape}{the shape of input VI observation? \code{line} or \code{point}}
@@ -46,6 +46,10 @@ of \code{t}, \code{y} and \code{QC_flag}. If not specified, it will be determine
 \item{angle}{\code{text.x} angle}
 
 \item{show.legend}{Boolean}
+
+\item{layer_extra}{(not used) extra ggplot layers}
+
+\item{...}{ignored}
 }
 \description{
 plot_curvefits
diff --git a/man/plot_season.Rd b/man/plot_season.Rd
index 4bcabe27..5cf94cc3 100644
--- a/man/plot_season.Rd
+++ b/man/plot_season.Rd
@@ -11,7 +11,9 @@ plot_season(
   IsPlot.OnlyBad = FALSE,
   show.legend = TRUE,
   ylab = "VI",
-  title = NULL
+  title = NULL,
+  show.shade = TRUE,
+  margin = 0.35
 )
 }
 \arguments{
@@ -32,6 +34,10 @@ original data.table.}
 \item{ylab}{y axis title}
 
 \item{title}{The main title (on top)}
+
+\item{show.shade}{Boolean, period inside growing cycle colored as shade?}
+
+\item{margin}{\verb{ylim = c(ymin, ymax + margin * A); A = ymax - ymin}.}
 }
 \description{
 Plot growing season divding result.
diff --git a/man/process_phenofit.Rd b/man/process_phenofit.Rd
index 6c8d165e..e2097d9a 100644
--- a/man/process_phenofit.Rd
+++ b/man/process_phenofit.Rd
@@ -14,6 +14,7 @@ process_phenofit(
   options_fitting = list(methods = c("AG", "Zhang", "Beck", "Elmore", "Gu"), wFUN =
     "wTSM", maxExtendMonth = 12, minExtendMonth = 0.5, use.y0 = FALSE),
   brks = NULL,
+  TRS = c(0.1, 0.2, 0.5, 0.6, 0.8, 0.9),
   run.curvefit = TRUE,
   ...
 )
@@ -35,3 +36,4 @@ dividing information.}
 \description{
 Extract Vegetation Phenology at site scale
 }
+\keyword{internal}
diff --git a/man/process_season.Rd b/man/process_season.Rd
new file mode 100644
index 00000000..42004ef9
--- /dev/null
+++ b/man/process_season.Rd
@@ -0,0 +1,38 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/process_season.R
+\name{process_season}
+\alias{process_season}
+\title{divide_seasons}
+\usage{
+process_season(
+  d_obs,
+  options = list(wFUN = "wTSM", maxExtendMonth = 12, MaxPeaksPerYear = 3,
+    MaxTroughsPerYear = 4),
+  nptperyear = 36,
+  south = FALSE,
+  .v_curve = FALSE,
+  ...
+)
+}
+\arguments{
+\item{nptperyear}{Integer, number of images per year.}
+
+\item{south}{Boolean. In south hemisphere, growing year is 1 July to the
+following year 31 June; In north hemisphere, growing year is 1 Jan to 31 Dec.}
+
+\item{.v_curve}{If true, it will use V-curve theory to optimize Whittaker
+parameter, lambda.}
+
+\item{...}{Others will be ignored.}
+
+\item{d}{data.frame, with the columns of \code{t}, \code{y} and \code{w}.}
+
+\item{options_season}{options of \code{\link[=season_mov]{season_mov()}}}
+}
+\description{
+divide_seasons
+}
+\note{
+site-year may be not continuous.
+}
+\keyword{internal}
diff --git a/man/reexports.Rd b/man/reexports.Rd
index 4ab626ef..beb7efc2 100644
--- a/man/reexports.Rd
+++ b/man/reexports.Rd
@@ -6,6 +6,7 @@
 \alias{\%<-\%}
 \alias{\%>\%}
 \alias{\%<>\%}
+\alias{make_date}
 \title{Objects exported from other packages}
 \keyword{internal}
 \description{
@@ -13,6 +14,8 @@ These objects are imported from other packages. Follow the links
 below to see their documentation.
 
 \describe{
+  \item{lubridate}{\code{\link[lubridate:make_datetime]{make_date}}}
+
   \item{magrittr}{\code{\link[magrittr:compound]{\%<>\%}}, \code{\link[magrittr:pipe]{\%>\%}}}
 
   \item{zeallot}{\code{\link[zeallot:operator]{\%<-\%}}}
diff --git a/man/season_mov.Rd b/man/season_mov.Rd
index c8f255e8..d64abdf8 100644
--- a/man/season_mov.Rd
+++ b/man/season_mov.Rd
@@ -1,53 +1,25 @@
 % Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/process_season.R, R/season_mov.R
-\name{process_season}
-\alias{process_season}
+% Please edit documentation in R/season_mov.R
+\name{season_mov}
 \alias{season_mov}
-\title{divide_seasons}
+\title{Moving growing season division}
 \usage{
-process_season(
-  d_obs,
-  options = list(wFUN = "wTSM", maxExtendMonth = 12, MaxPeaksPerYear = 3,
-    MaxTroughsPerYear = 4),
-  nptperyear = 36,
-  south = FALSE,
-  .v_curve = FALSE,
-  ...
-)
-
 season_mov(INPUT, options = list(r_min = 0), ..., years.run = NULL)
 }
 \arguments{
-\item{options}{see details}
-
-\item{nptperyear}{Integer, number of images per year.}
-
-\item{south}{Boolean. In south hemisphere, growing year is 1 July to the
-following year 31 June; In north hemisphere, growing year is 1 Jan to 31 Dec.}
+\item{INPUT}{A list object with the elements of \code{t}, \code{y}, \code{w},
+\code{Tn} (optional) and \code{ylu}, returned by \code{\link[=check_input]{check_input()}}.}
 
-\item{.v_curve}{If true, it will use V-curve theory to optimize Whittaker
-parameter, lambda.}
+\item{options}{see details}
 
 \item{...}{others to \code{\link[=season]{season()}}}
 
-\item{INPUT}{A list object with the elements of \code{t}, \code{y}, \code{w},
-\code{Tn} (optional) and \code{ylu}, returned by \code{\link[=check_input]{check_input()}}.}
-
 \item{years.run}{Numeric vector. Which years to run? If not specified, it is
 all years.}
-
-\item{d}{data.frame, with the columns of \code{t}, \code{y} and \code{w}.}
-
-\item{options_season}{options of \code{\link[=season_mov]{season_mov()}}}
 }
 \description{
-divide_seasons
-
 Moving growing season division
 }
-\note{
-site-year may be not continuous.
-}
 \section{options}{
 
 \itemize{
diff --git a/scripts b/scripts
index 34fa8643..19a2ede0 160000
--- a/scripts
+++ b/scripts
@@ -1 +1 @@
-Subproject commit 34fa864380ea089fc36d9de25158999c32f78767
+Subproject commit 19a2ede0faaa3034bef10bfb3e1a72e407544cb5