diff --git a/README.md b/README.md
index c07ca8a..89a7421 100644
--- a/README.md
+++ b/README.md
@@ -15,7 +15,8 @@ This study aimed at characterizing the rest-activity patterns in the general pop
 
 ## Findings
 
-* Four distinct rest-activity profiles were identified to describe aspects of overall amplitude, early rising time, prolonged daytime activity and biphasic pattern. [insert the profiles figure here. need access to upload figure]
+* Four distinct rest-activity profiles were identified to describe aspects of overall amplitude, early rising time, prolonged daytime activity and biphasic pattern. 
+![fPCs](image/fPCs.png)
 * Rest-activity profiles are highly associated with age, race, education levels and household income levels.
 * Rest-activity profiles differ on weekdays and weekends by demographics and socioeconomics.
 * Lower overall activity is associated with higher odds of self-reported poor or fair health.
@@ -25,7 +26,7 @@ This study aimed at characterizing the rest-activity patterns in the general pop
 In a nationally representative sample of US adults, we identified four distinct profiles for the 24-h rest-activity cycle. We found considerable variation in these profiles across different subgroups by age, gender, race/ethnicity, SES and work status. We also observed associations between rest-activity profiles and self-rated health status.
 
 ## Data
-A sample data of 50 randomly selected subjects from NHANES 2011-2014 with accelerometer data can be found [here](NHANES_sample_data.rds).
+A sample data of 50 randomly selected subjects from NHANES 2011-2014 with accelerometer data can be found [here](data/NHANES_sample_data.rds).
 
 This dataset includes the following variables:
 - SEQN
@@ -36,13 +37,13 @@ This dataset includes the following variables:
 
 ## Codes
 
-The main code to run fpca analysis on the sample data can be found [here](fPCA_analysis.R)
+The main code to run fpca analysis on the sample data can be found [here](code/fPCA_analysis.R)
 
-R function to [plot a functional object](plot_fd_new1.R)
+R function to [plot a functional object](code/plot_fd_new1.R)
 - code adapted from plot.fd{fda}
 - changes are made to present labels of xasis as clock time
 
-R function to generate Figure 1 in the published paper of [functional principal components](plot_fpca_new1.R)
+R function to generate Figure 1 in the published paper of [functional principal components](code/plot_fpca_new1.R)
 - code adapted from plot.pca.fd{fda}
 - add argument 'flip' to decide which component(s) to have sign(s) flipped ('+' to '-' or vice versa)
 - use different colors to represent component's sign
diff --git a/fpca_NHANES.Rproj b/code/fpca_NHANES.Rproj
similarity index 100%
rename from fpca_NHANES.Rproj
rename to code/fpca_NHANES.Rproj
diff --git a/fpca_analysis.R b/code/fpca_analysis.R
similarity index 100%
rename from fpca_analysis.R
rename to code/fpca_analysis.R
diff --git a/plot_fd_new1.R b/code/plot_fd_new1.R
old mode 100755
new mode 100644
similarity index 96%
rename from plot_fd_new1.R
rename to code/plot_fd_new1.R
index 6f2f6c8..039f56b
--- a/plot_fd_new1.R
+++ b/code/plot_fd_new1.R
@@ -1,205 +1,205 @@
-# 7/11/21 change color
-plot.fd.new <- function (x, y, Lfdobj = 0, href = TRUE, titles = NULL, xlim = NULL, 
-          ylim = NULL, xlab = NULL, ylab = NULL, ask = FALSE, nx = NULL, 
-          axes = NULL, family = NULL, ...) {
-  fdobj <- x
-  if (!(is.fd(fdobj) || is.fdPar(fdobj))) 
-    stop(paste("First argument is neither a functional data or a ", 
-               "functional parameter object."))
-  if (is.fdPar(fdobj)) 
-    fdobj <- fdobj$fd
-  if (is.null(axes)) {
-    if (is.null(fdobj$basis$axes)) {
-      Axes <- TRUE
-      axFun <- FALSE
-    }
-    else {
-      if (!inherits(fdobj$basis$axes, "list")) 
-        stop("fdobj$basis$axes must be a list;  ", 
-             "class(fdobj$basis$axes) = ", class(fdobj$basis$axes))
-      if (!(inherits(fdobj$basis$axes[[1]], "character") || 
-            inherits(fdobj$basis$axes[[1]], "function"))) 
-        stop("fdobj$basis$axes[[1]] must be either a function or the ", 
-             "name of a function;  class(fdobj$basis$axes[[1]]) = ", 
-             class(fdobj$basis$axes[[1]]))
-      Axes <- FALSE
-      axFun <- TRUE
-      axList <- c(fdobj$basis$axes, ...)
-    }
-  }
-  else {
-    if (is.logical(axes)) {
-      Axes <- axes
-      axFun <- FALSE
-    }
-    else {
-      if (!inherits(axes, "list")) 
-        stop("axes must be a logical or a list;  class(axes) = ", 
-             class(axes))
-      if (!(inherits(axes[[1]], "character") || inherits(axes[[1]], 
-                                                         "function"))) 
-        stop("axes[[1]] must be either a function or the ", 
-             "name of a function;  class(axes[[1]]) = ", 
-             class(axes[[1]]))
-      Axes <- FALSE
-      axFun <- TRUE
-      axList <- c(axes, ...)
-    }
-  }
-  Lfdobj <- int2Lfd(Lfdobj)
-  if (!inherits(Lfdobj, "Lfd")) 
-    stop("Second argument is not a linear differential operator.")
-  coef <- fdobj$coefs
-  coefd <- dim(coef)
-  ndim <- length(coefd)
-  nbasis <- coefd[1]
-  if (is.null(nx)) 
-    nx <- max(c(501, 10 * nbasis + 1))
-  nrep <- coefd[2]
-  if (ndim > 2) 
-    nvar <- coefd[3]
-  else nvar <- 1
-  basisobj <- fdobj$basis
-  rangex <- basisobj$rangeval
-  if (missing(y)) {
-    y <- nx
-  }
-  else {
-    if (is.numeric(y)) 
-      y <- as.vector(y)
-  }
-  Y <- y
-  if (length(y) == 1) {
-    if (y >= 1) {
-      y <- seq(rangex[1], rangex[2], len = round(y))
-    }
-    else {
-      stop("'y' a single number less than one.")
-    }
-  }
-  if (min(y) < rangex[1] || max(y) > rangex[2]) 
-    stop("Values in Y are outside the basis range.")
-  if (is.null(xlim)) {
-    xlim <- rangex
-  }
-  else {
-    rangex[1] <- max(rangex[1], xlim[1])
-    rangex[2] <- min(rangex[2], xlim[2])
-    if (length(Y) == 1) 
-      y <- seq(rangex[1], rangex[2], len = round(Y))
-  }
-  fdmat <- eval.fd(y, fdobj, Lfdobj)
-  rangey <- range(fdmat)
-  if (is.null(ylim)) 
-    ylim <- rangey
-  fdnames = fdobj$fdnames
-  fdlabelslist = fdlabels(fdnames, nrep, nvar)
-  xlabel = fdlabelslist$xlabel
-  ylabel = fdlabelslist$ylabel
-  casenames = fdlabelslist$casenames
-  varnames = fdlabelslist$varnames
-  if (is.null(xlab)) 
-    xlab <- xlabel
-  if (is.null(ylab)) 
-    ylab <- ylabel
-  if (ndim < 2) {
-    plot(y, fdmat, type = "l", xlim = xlim, ylim = ylim, 
-         xlab = xlab, ylab = ylab, axes = Axes
-         , xaxt="n", ...)
-    # xaxis labels and locations
-    x.axis <- seq(from=0, to=24, by=1)
-    x.axis.lab <- paste0(x.axis,":00")
-    x.axis.loc <- x.axis/24*288
-    axis(1, at=x.axis.loc, labels=x.axis.lab, las=0) 
-    
-    if (axFun) 
-      do.call(axList[[1]], axList[-1])
-    if (zerofind(fdmat) && href) 
-      abline(h = 0, lty = 2)
-  }
-  if (ndim == 2) {
-    if (!ask) {
-      matplot(y, fdmat, type = "l", xlim = xlim, 
-              ylim = ylim, xlab = xlab, ylab = ylab, axes = Axes
-              , xaxt="n"
-              #, col="blue1"
-              , family=family, ps=11
-              , ...
-              )
-      # xaxis labels and locations
-      # x.axis <- seq(from=0, to=24, by=1)
-      x.axis <- seq(from=0, to=24, by=4)
-      x.axis.lab <- paste0(x.axis,":00")
-      x.axis.loc <- x.axis/24*288
-      axis(1, at=x.axis.loc, labels=x.axis.lab, las=0) 
-      
-      if (axFun) 
-        do.call(axList[[1]], axList[-1])
-      if (zerofind(fdmat) && href) 
-        abline(h = 0, lty = 2)
-    }
-    else {
-      op <- par(ask = FALSE)
-      on.exit(par(op))
-      cat("Multiple plots:  Click in the plot to advance to the next")
-      for (irep in 1:nrep) {
-        plot(y, fdmat[, irep], type = "l", xlim = xlim, 
-             ylim = ylim, xlab = xlab, ylab = ylab, axes = Axes, 
-             ...)
-        if (axFun) 
-          do.call(axList[[1]], axList[-1])
-        if (irep < 2) 
-          par(ask = ask)
-        if (!is.null(casenames)) 
-          title(casenames[irep])
-        else title(paste("Case", irep))
-        if (zerofind(ylim) && href) 
-          abline(h = 0, lty = 2)
-      }
-    }
-  }
-  if (ndim == 3) {
-    if (!ask) {
-      for (ivar in 1:nvar) {
-        matplot(y, fdmat[, , ivar], type = "l", 
-                xlim = xlim, ylim = ylim, xlab = xlab, ylab = ylab, 
-                ask = FALSE, axes = Axes, ...)
-        if (axFun) 
-          do.call(axList[[1]], axList[-1])
-        if (!is.null(varnames)) 
-          title(varnames[ivar])
-        else title(paste("Variable", ivar))
-        if (zerofind(ylim) && href) 
-          abline(h = 0, lty = 2)
-      }
-    }
-    else {
-      op <- par(ask = FALSE)
-      on.exit(par(op))
-      cat("Multiple plots:  Click in the plot to advance to the next")
-      for (irep in 1:nrep) {
-        for (ivar in 1:nvar) {
-          plot(y, fdmat[, irep, ivar], type = "l", 
-               xlim = xlim, ylim = ylim, xlab = xlab, ylab = ylab, 
-               axes = Axes, ...)
-          if (axFun) 
-            do.call(axList[[1]], axList[-1])
-          if (!is.null(casenames)) 
-            titlestr = casenames[irep]
-          else titlestr = paste("Case", irep)
-          if (!is.null(varnames)) {
-            titlestr = paste(titlestr, "  ", varnames[ivar])
-          }
-          else {
-            titlestr = paste(titlestr, "  ", "Variable", 
-                             ivar)
-          }
-          title(titlestr)
-          if (zerofind(ylim) && href) 
-            abline(h = 0, lty = 2)
-        }
-      }
-    }
-  }
-  "done"
-}
+# 7/11/21 change color
+plot.fd.new <- function (x, y, Lfdobj = 0, href = TRUE, titles = NULL, xlim = NULL, 
+          ylim = NULL, xlab = NULL, ylab = NULL, ask = FALSE, nx = NULL, 
+          axes = NULL, family = NULL, ...) {
+  fdobj <- x
+  if (!(is.fd(fdobj) || is.fdPar(fdobj))) 
+    stop(paste("First argument is neither a functional data or a ", 
+               "functional parameter object."))
+  if (is.fdPar(fdobj)) 
+    fdobj <- fdobj$fd
+  if (is.null(axes)) {
+    if (is.null(fdobj$basis$axes)) {
+      Axes <- TRUE
+      axFun <- FALSE
+    }
+    else {
+      if (!inherits(fdobj$basis$axes, "list")) 
+        stop("fdobj$basis$axes must be a list;  ", 
+             "class(fdobj$basis$axes) = ", class(fdobj$basis$axes))
+      if (!(inherits(fdobj$basis$axes[[1]], "character") || 
+            inherits(fdobj$basis$axes[[1]], "function"))) 
+        stop("fdobj$basis$axes[[1]] must be either a function or the ", 
+             "name of a function;  class(fdobj$basis$axes[[1]]) = ", 
+             class(fdobj$basis$axes[[1]]))
+      Axes <- FALSE
+      axFun <- TRUE
+      axList <- c(fdobj$basis$axes, ...)
+    }
+  }
+  else {
+    if (is.logical(axes)) {
+      Axes <- axes
+      axFun <- FALSE
+    }
+    else {
+      if (!inherits(axes, "list")) 
+        stop("axes must be a logical or a list;  class(axes) = ", 
+             class(axes))
+      if (!(inherits(axes[[1]], "character") || inherits(axes[[1]], 
+                                                         "function"))) 
+        stop("axes[[1]] must be either a function or the ", 
+             "name of a function;  class(axes[[1]]) = ", 
+             class(axes[[1]]))
+      Axes <- FALSE
+      axFun <- TRUE
+      axList <- c(axes, ...)
+    }
+  }
+  Lfdobj <- int2Lfd(Lfdobj)
+  if (!inherits(Lfdobj, "Lfd")) 
+    stop("Second argument is not a linear differential operator.")
+  coef <- fdobj$coefs
+  coefd <- dim(coef)
+  ndim <- length(coefd)
+  nbasis <- coefd[1]
+  if (is.null(nx)) 
+    nx <- max(c(501, 10 * nbasis + 1))
+  nrep <- coefd[2]
+  if (ndim > 2) 
+    nvar <- coefd[3]
+  else nvar <- 1
+  basisobj <- fdobj$basis
+  rangex <- basisobj$rangeval
+  if (missing(y)) {
+    y <- nx
+  }
+  else {
+    if (is.numeric(y)) 
+      y <- as.vector(y)
+  }
+  Y <- y
+  if (length(y) == 1) {
+    if (y >= 1) {
+      y <- seq(rangex[1], rangex[2], len = round(y))
+    }
+    else {
+      stop("'y' a single number less than one.")
+    }
+  }
+  if (min(y) < rangex[1] || max(y) > rangex[2]) 
+    stop("Values in Y are outside the basis range.")
+  if (is.null(xlim)) {
+    xlim <- rangex
+  }
+  else {
+    rangex[1] <- max(rangex[1], xlim[1])
+    rangex[2] <- min(rangex[2], xlim[2])
+    if (length(Y) == 1) 
+      y <- seq(rangex[1], rangex[2], len = round(Y))
+  }
+  fdmat <- eval.fd(y, fdobj, Lfdobj)
+  rangey <- range(fdmat)
+  if (is.null(ylim)) 
+    ylim <- rangey
+  fdnames = fdobj$fdnames
+  fdlabelslist = fdlabels(fdnames, nrep, nvar)
+  xlabel = fdlabelslist$xlabel
+  ylabel = fdlabelslist$ylabel
+  casenames = fdlabelslist$casenames
+  varnames = fdlabelslist$varnames
+  if (is.null(xlab)) 
+    xlab <- xlabel
+  if (is.null(ylab)) 
+    ylab <- ylabel
+  if (ndim < 2) {
+    plot(y, fdmat, type = "l", xlim = xlim, ylim = ylim, 
+         xlab = xlab, ylab = ylab, axes = Axes
+         , xaxt="n", ...)
+    # xaxis labels and locations
+    x.axis <- seq(from=0, to=24, by=1)
+    x.axis.lab <- paste0(x.axis,":00")
+    x.axis.loc <- x.axis/24*288
+    axis(1, at=x.axis.loc, labels=x.axis.lab, las=0) 
+    
+    if (axFun) 
+      do.call(axList[[1]], axList[-1])
+    if (zerofind(fdmat) && href) 
+      abline(h = 0, lty = 2)
+  }
+  if (ndim == 2) {
+    if (!ask) {
+      matplot(y, fdmat, type = "l", xlim = xlim, 
+              ylim = ylim, xlab = xlab, ylab = ylab, axes = Axes
+              , xaxt="n"
+              #, col="blue1"
+              , family=family, ps=11
+              , ...
+              )
+      # xaxis labels and locations
+      # x.axis <- seq(from=0, to=24, by=1)
+      x.axis <- seq(from=0, to=24, by=4)
+      x.axis.lab <- paste0(x.axis,":00")
+      x.axis.loc <- x.axis/24*288
+      axis(1, at=x.axis.loc, labels=x.axis.lab, las=0) 
+      
+      if (axFun) 
+        do.call(axList[[1]], axList[-1])
+      if (zerofind(fdmat) && href) 
+        abline(h = 0, lty = 2)
+    }
+    else {
+      op <- par(ask = FALSE)
+      on.exit(par(op))
+      cat("Multiple plots:  Click in the plot to advance to the next")
+      for (irep in 1:nrep) {
+        plot(y, fdmat[, irep], type = "l", xlim = xlim, 
+             ylim = ylim, xlab = xlab, ylab = ylab, axes = Axes, 
+             ...)
+        if (axFun) 
+          do.call(axList[[1]], axList[-1])
+        if (irep < 2) 
+          par(ask = ask)
+        if (!is.null(casenames)) 
+          title(casenames[irep])
+        else title(paste("Case", irep))
+        if (zerofind(ylim) && href) 
+          abline(h = 0, lty = 2)
+      }
+    }
+  }
+  if (ndim == 3) {
+    if (!ask) {
+      for (ivar in 1:nvar) {
+        matplot(y, fdmat[, , ivar], type = "l", 
+                xlim = xlim, ylim = ylim, xlab = xlab, ylab = ylab, 
+                ask = FALSE, axes = Axes, ...)
+        if (axFun) 
+          do.call(axList[[1]], axList[-1])
+        if (!is.null(varnames)) 
+          title(varnames[ivar])
+        else title(paste("Variable", ivar))
+        if (zerofind(ylim) && href) 
+          abline(h = 0, lty = 2)
+      }
+    }
+    else {
+      op <- par(ask = FALSE)
+      on.exit(par(op))
+      cat("Multiple plots:  Click in the plot to advance to the next")
+      for (irep in 1:nrep) {
+        for (ivar in 1:nvar) {
+          plot(y, fdmat[, irep, ivar], type = "l", 
+               xlim = xlim, ylim = ylim, xlab = xlab, ylab = ylab, 
+               axes = Axes, ...)
+          if (axFun) 
+            do.call(axList[[1]], axList[-1])
+          if (!is.null(casenames)) 
+            titlestr = casenames[irep]
+          else titlestr = paste("Case", irep)
+          if (!is.null(varnames)) {
+            titlestr = paste(titlestr, "  ", varnames[ivar])
+          }
+          else {
+            titlestr = paste(titlestr, "  ", "Variable", 
+                             ivar)
+          }
+          title(titlestr)
+          if (zerofind(ylim) && href) 
+            abline(h = 0, lty = 2)
+        }
+      }
+    }
+  }
+  "done"
+}
diff --git a/plot_fpca_new1.R b/code/plot_fpca_new1.R
old mode 100755
new mode 100644
similarity index 96%
rename from plot_fpca_new1.R
rename to code/plot_fpca_new1.R
index 6c1016c..e9e1638
--- a/plot_fpca_new1.R
+++ b/code/plot_fpca_new1.R
@@ -1,90 +1,90 @@
-### based on source code of plot.pca.fd
-### revisions:
-### 1) colored add/subtract pca components
-### 2) change xaxis labels to clock time of 0:00 - 23:00
-### 3) add the argument of 'flip' to determine which components to flip the sign
-### warning: depends on the results and interpretation, may not suitable for other data or analysis. 
-
-plot.fpca.new <- function (x, nx = 128, pointplot = TRUE, harm = 0, expand = 0, 
-          cycle = FALSE, flip=c(1,2,4)
-          , ...) {
-  pcafd <- x
-  
-  # xaxis labels and locations
-  x.axis <- seq(from=0, to=24, by=1)
-  x.axis.lab <- paste0(x.axis,":00")
-  x.axis.loc <- x.axis/24*288
-  
-  if (!(inherits(pcafd, "pca.fd"))) 
-    stop("Argument 'x' is not a pca.fd object.")
-  harmfd <- pcafd[[1]]
-  basisfd <- harmfd$basis
-  rangex <- basisfd$rangeval
-  {
-    if (length(nx) > 1) {
-      x <- nx
-      nx <- length(x)
-    }
-    else x <- seq(rangex[1], rangex[2], length = nx)
-  }
-  fdmat <- eval.fd(x, harmfd)
-  meanmat <- eval.fd(x, pcafd$meanfd)
-  dimfd <- dim(fdmat)
-  nharm <- dimfd[2]
-  plotsPerPg <- sum(par("mfrow"))
-  harm <- as.vector(harm)
-  if (harm[1] == 0) 
-    harm <- (1:nharm)
-  if (length(dimfd) == 2) {
-    for (jharm in 1:length(harm)) {
-      #jharm=1
-      if (jharm == 2) {
-        op <- par(ask = TRUE)
-        on.exit(par(op))
-      }
-      iharm <- harm[jharm]
-      if (expand == 0) {
-        fac <- sqrt(pcafd$values[iharm])
-      }
-      else {
-        fac <- expand
-      }
-      vecharm <- fdmat[, iharm]
-      #if (jharm == 3){ # fPCA3
-      if (jharm %in% flip){
-        pcmat <- cbind(meanmat + fac * (-vecharm)
-                     , meanmat - fac * (-vecharm)
-                     )
-      }
-      #else{ # fPCA1,2,4 flip the sign 
-      else{
-        pcmat <- cbind(meanmat + fac * vecharm
-                       , meanmat - fac * vecharm
-        )
-      }
-      
-      if (pointplot) 
-        plottype <- "p"
-      else plottype <- "l"
-      percentvar <- round(100 * pcafd$varprop[iharm], 1)
-      plot(x, meanmat, type = "l", ylim = c(min(pcmat), 
-                                            max(pcmat)), ylab = paste("Harmonic", iharm), 
-           main = paste("PCA function", iharm, "(Percentage of variability", 
-                        percentvar, ")")
-           , xaxt="n", xlab=NA , ...
-           )
-      axis(1, at=x.axis.loc, labels=x.axis.lab, las=0) 
-      
-      if (pointplot) {
-        points(x, pcmat[, 1], pch = "+", col="darkgoldenrod1")
-        points(x, pcmat[, 2], pch = "-", col="blue1")
-      }
-      else {
-        lines(x, pcmat[, 1], lty = 2, col="darkgoldenrod1")
-        lines(x, pcmat[, 2], lty = 3, col="blue1")
-      }
-    }
-  }
-#  this is not a complete version as original plot.pca.fd function, the rest starting from here is deleted.
-  invisible(NULL)
-}
+### based on source code of plot.pca.fd
+### revisions:
+### 1) colored add/subtract pca components
+### 2) change xaxis labels to clock time of 0:00 - 23:00
+### 3) add the argument of 'flip' to determine which components to flip the sign
+### warning: depends on the results and interpretation, may not suitable for other data or analysis. 
+
+plot.fpca.new <- function (x, nx = 128, pointplot = TRUE, harm = 0, expand = 0, 
+          cycle = FALSE, flip=c(1,2,4)
+          , ...) {
+  pcafd <- x
+  
+  # xaxis labels and locations
+  x.axis <- seq(from=0, to=24, by=1)
+  x.axis.lab <- paste0(x.axis,":00")
+  x.axis.loc <- x.axis/24*288
+  
+  if (!(inherits(pcafd, "pca.fd"))) 
+    stop("Argument 'x' is not a pca.fd object.")
+  harmfd <- pcafd[[1]]
+  basisfd <- harmfd$basis
+  rangex <- basisfd$rangeval
+  {
+    if (length(nx) > 1) {
+      x <- nx
+      nx <- length(x)
+    }
+    else x <- seq(rangex[1], rangex[2], length = nx)
+  }
+  fdmat <- eval.fd(x, harmfd)
+  meanmat <- eval.fd(x, pcafd$meanfd)
+  dimfd <- dim(fdmat)
+  nharm <- dimfd[2]
+  plotsPerPg <- sum(par("mfrow"))
+  harm <- as.vector(harm)
+  if (harm[1] == 0) 
+    harm <- (1:nharm)
+  if (length(dimfd) == 2) {
+    for (jharm in 1:length(harm)) {
+      #jharm=1
+      if (jharm == 2) {
+        op <- par(ask = TRUE)
+        on.exit(par(op))
+      }
+      iharm <- harm[jharm]
+      if (expand == 0) {
+        fac <- sqrt(pcafd$values[iharm])
+      }
+      else {
+        fac <- expand
+      }
+      vecharm <- fdmat[, iharm]
+      #if (jharm == 3){ # fPCA3
+      if (jharm %in% flip){
+        pcmat <- cbind(meanmat + fac * (-vecharm)
+                     , meanmat - fac * (-vecharm)
+                     )
+      }
+      #else{ # fPCA1,2,4 flip the sign 
+      else{
+        pcmat <- cbind(meanmat + fac * vecharm
+                       , meanmat - fac * vecharm
+        )
+      }
+      
+      if (pointplot) 
+        plottype <- "p"
+      else plottype <- "l"
+      percentvar <- round(100 * pcafd$varprop[iharm], 1)
+      plot(x, meanmat, type = "l", ylim = c(min(pcmat), 
+                                            max(pcmat)), ylab = paste("Harmonic", iharm), 
+           main = paste("PCA function", iharm, "(Percentage of variability", 
+                        percentvar, ")")
+           , xaxt="n", xlab=NA , ...
+           )
+      axis(1, at=x.axis.loc, labels=x.axis.lab, las=0) 
+      
+      if (pointplot) {
+        points(x, pcmat[, 1], pch = "+", col="darkgoldenrod1")
+        points(x, pcmat[, 2], pch = "-", col="blue1")
+      }
+      else {
+        lines(x, pcmat[, 1], lty = 2, col="darkgoldenrod1")
+        lines(x, pcmat[, 2], lty = 3, col="blue1")
+      }
+    }
+  }
+#  this is not a complete version as original plot.pca.fd function, the rest starting from here is deleted.
+  invisible(NULL)
+}
diff --git a/NHANES_sample_data.rds b/data/NHANES_sample_data.rds
similarity index 100%
rename from NHANES_sample_data.rds
rename to data/NHANES_sample_data.rds
diff --git a/NHANES_sample_data_day_average.rds b/data/NHANES_sample_data_day_average.rds
similarity index 100%
rename from NHANES_sample_data_day_average.rds
rename to data/NHANES_sample_data_day_average.rds
diff --git a/image/fPCs.png b/image/fPCs.png
new file mode 100644
index 0000000..d9e8463
Binary files /dev/null and b/image/fPCs.png differ