diff --git a/.Rbuildignore b/.Rbuildignore
index d651831..5d5ca31 100644
--- a/.Rbuildignore
+++ b/.Rbuildignore
@@ -4,3 +4,5 @@
 ^tests$
 ^cran-comments\.md$
 ^CRAN-SUBMISSION$
+^.*\.Rproj$
+^\.Rproj\.user$
diff --git a/.gitignore b/.gitignore
index f304eec..a977dd1 100644
--- a/.gitignore
+++ b/.gitignore
@@ -76,3 +76,4 @@ inst/doc
 
 # Development
 experiment/
+.Rproj.user
diff --git a/R/CheckOptions.R b/R/CheckOptions.R
index c24a6e6..edeac87 100644
--- a/R/CheckOptions.R
+++ b/R/CheckOptions.R
@@ -90,7 +90,7 @@ CheckOptions = function(t,optns,n){
   # if (optns[['methodXi']] == 'IN' && optns[['dataType']] != 'Dense') {
     # stop("integration method can only be applied on dense data now!")
   # }
-  if(!(any(optns[['kernel']] == c('epan','gauss','rect','quar','gausvar')))){ 
+  if(!(any(optns[['kernel']] == c('epan','gauss','rect','quar','gausvar','triangular','triweight','tricube','cosine','logistic','sigmoid','silverman')))){ 
     #method to estimate the PC scores
     stop("FPCA is aborted because the argument: kernel is invalid!\n");   
   }
diff --git a/R/FCCor.R b/R/FCCor.R
index 6d70a6f..857147e 100644
--- a/R/FCCor.R
+++ b/R/FCCor.R
@@ -4,7 +4,7 @@
 #' @param y A list of function values corresponding to the second process.
 #' @param Lt A list of time points for both \code{x} and \code{y}.
 #' @param bw A numeric vector for bandwidth of length either 5 or 1, specifying the bandwidths for E(X), E(Y), var(X), var(Y), and cov(X, Y). If \code{bw} is a scalar then all five bandwidths are chosen to be the same. 
-#' @param kern Smoothing kernel for mu and covariance; "rect", "gauss", "epan", "gausvar", "quar" (default: "gauss")
+#' @param kern Smoothing kernel for mu and covariance; "rect", "gauss", "epan", "gausvar", "quar", "triangular", "triweight", "tricube", "cosine", "logistic", "sigmoid", "silverman" (default: "gauss")
 #' @param Tout Output time points. Default to the sorted unique time points. 
 #'
 #' @details \code{FCCor} calculate only the concurrent correlation corr(X(t), Y(t)) (note that the time points t are the same). It assumes no measurement error in the observed values.
diff --git a/R/FCReg.R b/R/FCReg.R
index 99bb0d9..6a80c95 100644
--- a/R/FCReg.R
+++ b/R/FCReg.R
@@ -6,7 +6,7 @@
 #' @param userBwMu A scalar with bandwidth used for smoothing the mean
 #' @param userBwCov A scalar with bandwidth used for smoothing the auto- and cross-covariances
 #' @param outGrid A vector with the output time points
-#' @param kern Smoothing kernel choice, common for mu and covariance; "rect", "gauss", "epan", "gausvar", "quar" (default: "gauss")
+#' @param kern Smoothing kernel choice, common for mu and covariance; "rect", "gauss", "epan", "gausvar", "quar", "triangular", "triweight", "tricube", "cosine", "logistic", "sigmoid", "silverman" (default: "gauss")
 #' @param measurementError Indicator measurement errors on the functional observations should be assumed. If TRUE the diagonal raw covariance will be removed when smoothing. (default: TRUE)
 #' @param diag1D  A string specifying whether to use 1D smoothing for the diagonal line of the covariance. 
 #' 'none': don't use 1D smoothing; 'cross': use 1D only for cross-covariances; 'all': use 1D for both auto- and cross-covariances. (default : 'none')
diff --git a/R/FPCA.R b/R/FPCA.R
index f2ee6ef..10004ee 100644
--- a/R/FPCA.R
+++ b/R/FPCA.R
@@ -25,7 +25,7 @@
 #' \item{fitEigenValues}{Whether also to obtain a regression fit of the eigenvalues - default: FALSE}
 #' \item{FVEthreshold}{Fraction-of-Variance-Explained threshold used during the SVD of the fitted covariance function; numeric (0,1] - default: 0.99}
 #' \item{FVEfittedCov}{Fraction-of-Variance explained by the components that are used to construct fittedCov; numeric (0,1] - default: NULL (all components available will be used)}
-#' \item{kernel}{Smoothing kernel choice, common for mu and covariance; "rect", "gauss", "epan", "gausvar", "quar" - default: "gauss"; dense data are assumed noise-less so no smoothing is performed. }
+#' \item{kernel}{Smoothing kernel choice, common for mu and covariance; "rect", "gauss", "epan", "gausvar", "quar", "triangular", "triweight", "tricube", "cosine", "logistic", "sigmoid", "silverman" - default: "gauss"; dense data are assumed noise-less so no smoothing is performed. }
 #' \item{kFoldMuCov}{The number of folds to be used for mean and covariance smoothing. Default: 10}
 #' \item{lean}{If TRUE the 'inputData' field in the output list is empty. Default: FALSE}
 #' \item{maxK}{The maximum number of principal components to consider - default: min(20, N-2,nRegGrid-2), N:# of curves, nRegGrid:# of support points in each direction of covariance surface}
diff --git a/R/FSVD.R b/R/FSVD.R
index 08aa045..26e0493 100644
--- a/R/FSVD.R
+++ b/R/FSVD.R
@@ -18,7 +18,7 @@
 #' \item{userMu1}{The user defined mean of sample 1 used to centre it prior to the cross-covariance estimation. - default: determine automatically based by the FPCA of sample 1}
 #' \item{userMu2}{The user defined mean of sample 2 used to centre it prior to the cross-covariance estimation. - default: determine automatically based by the FPCA of sample 2}
 #' \item{maxK}{The maximum number of singular components to consider; default: min(20, N-1), N:# of curves.}
-#' \item{kernel}{Smoothing kernel choice, common for mu and covariance; "rect", "gauss", "epan", "gausvar", "quar" - default: "gauss"; dense data are assumed noise-less so no smoothing is performed.}
+#' \item{kernel}{Smoothing kernel choice, common for mu and covariance; "rect", "gauss", "epan", "gausvar", "quar", "triangular", "triweight", "tricube", "cosine", "logistic", "sigmoid", "silverman" - default: "gauss"; dense data are assumed noise-less so no smoothing is performed.}
 #' \item{rmDiag}{Logical describing if the routine should remove diagonal raw cov for cross cov estimation (default: FALSE) }
 #' \item{noScores}{Logical describing if the routine should return functional singular scores or not (default: TRUE) }
 #' \item{regulRS}{String describing if the regularisation of the composite cross-covariance matrix should be done using 'sigma1' or 'rho' (see ?FPCA for details) (default: 'sigma2') }
diff --git a/R/GCVLwls1D1.R b/R/GCVLwls1D1.R
index 07be196..5a1ae75 100644
--- a/R/GCVLwls1D1.R
+++ b/R/GCVLwls1D1.R
@@ -40,7 +40,7 @@ GCVLwls1D1 <- function(yy,tt, kernel, npoly, nder, dataType, verbose=TRUE) {
     h0 = 1.5 * Minb(t,npoly+1);
   }  
   if ( is.nan(h0) ){
-    if ( kernel == "gauss" ){
+    if ( kernel in c("gauss", "logistic", "sigmoid", "silverman") ){
       h0 = 0.2 * r;
     }else{
       stop("The data is too sparse, no suitable bandwidth can be found! Try Gaussian kernel instead!\n")
@@ -56,7 +56,9 @@ GCVLwls1D1 <- function(yy,tt, kernel, npoly, nder, dataType, verbose=TRUE) {
   # I would write them in a function (equivalent of mykernel.m) if it is worth it 
   # Similarly there is no reason to repeat the FOR-loop twice; this too can go into a seperate function
   k0_candidates <- list('quar' = 0.9375,  'epan' = 0.7500, 'rect' = 0.5000, 
-                        'gausvar' = 0.498677, 'gausvar1' = 0.598413,  'gausvar2' = 0.298415, 'other' = 0.398942)
+                        'gausvar' = 0.498677, 'gausvar1' = 0.598413,  'gausvar2' = 0.298415, 'other' = 0.398942,
+                        'triangular' = 1.000, 'triweight' = 35/32, 'tricube' = 70/81,
+                        'cosine' = 0.7853982, 'logistic' = 0.25, 'sigmoid' = 1/pi, 'silverman' = 0.35355339059)
   if( any(names(k0_candidates) == kernel)){ 
     k0 = as.numeric(k0_candidates[kernel])
   } else { 
diff --git a/R/GCVLwls2DV2.R b/R/GCVLwls2DV2.R
index f4a895c..8185b3c 100644
--- a/R/GCVLwls2DV2.R
+++ b/R/GCVLwls2DV2.R
@@ -244,6 +244,20 @@ KernelAt0 <- function(kern) {
     k0 <- 0.498677850501791
   else if (kern == 'gauss')
     k0 <- 0.398942280401433
+  else if (kern == 'triangular')
+    k0 <- 1.0
+  else if (kern == 'triweight')
+    k0 <- 35/32
+  else if (kern == 'tricube')
+    k0 <- 70/81
+  else if (kern == 'cosine')
+    k0 <- 0.7853982
+  else if (kern == 'logistic')
+    k0 <- 0.25
+  else if (kern == 'sigmoid')
+    k0 <- 1/pi
+  else if (kern == 'silverman')
+    k0 <- 0.35355339059
   else
     stop('Unknown kernel')
   
diff --git a/R/GetCovSurface.R b/R/GetCovSurface.R
index 8e23584..ed8b97e 100644
--- a/R/GetCovSurface.R
+++ b/R/GetCovSurface.R
@@ -14,7 +14,7 @@
 #' \item{methodBwMu}{The bandwidth choice method for the mean function; 'GMeanAndGCV' (the geometric mean of the GCV bandwidth and the minimum bandwidth),'CV','GCV' - default: 5\% of the support} 
 #' \item{dataType}{The type of design we have (usually distinguishing between sparse or dense functional data); 'Sparse', 'Dense', 'DenseWithMV', 'p>>n' - default:  determine automatically based on 'IsRegular'}
 #' \item{error}{Assume measurement error in the dataset; logical - default: TRUE}
-#' \item{kernel}{Smoothing kernel choice, common for mu and covariance; "rect", "gauss", "epan", "gausvar", "quar" - default: "gauss"; dense data are assumed noise-less so no smoothing is performed. }
+#' \item{kernel}{Smoothing kernel choice, common for mu and covariance; "rect", "gauss", "epan", "gausvar", "quar", "triangular", "triweight", "tricube", "cosine", "logistic", "sigmoid", "silverman" - default: "gauss"; dense data are assumed noise-less so no smoothing is performed. }
 #' \item{kFoldMuCov}{The number of folds to be used for mean and covariance smoothing. Default: 10}
 #' \item{lean}{If TRUE the 'inputData' field in the output list is empty. Default: FALSE}
 #' \item{methodMuCovEst}{The method to estimate the mean and covariance in the case of dense functional data; 'cross-sectional', 'smooth' - default: 'cross-sectional'}
diff --git a/R/GetMeanCurve.R b/R/GetMeanCurve.R
index ca6cba5..9c3f37d 100644
--- a/R/GetMeanCurve.R
+++ b/R/GetMeanCurve.R
@@ -12,7 +12,7 @@
 #' \item{methodBwMu}{The bandwidth choice method for the mean function; 'GMeanAndGCV' (the geometric mean of the GCV bandwidth and the minimum bandwidth),'CV','GCV' - default: 5\% of the support} 
 #' \item{dataType}{The type of design we have (usually distinguishing between sparse or dense functional data); 'Sparse', 'Dense', 'DenseWithMV', 'p>>n' - default:  determine automatically based on 'IsRegular'}
 #' \item{plot}{Plot mean curve; logical - default: FALSE}
-#' \item{kernel}{Smoothing kernel choice, common for mu and covariance; "rect", "gauss", "epan", "gausvar", "quar" - default: "gauss"; dense data are assumed noise-less so no smoothing is performed. }
+#' \item{kernel}{Smoothing kernel choice, common for mu and covariance; "rect", "gauss", "epan", "gausvar", "quar", "triangular", "triweight", "tricube", "cosine", "logistic", "sigmoid", "silverman" - default: "gauss"; dense data are assumed noise-less so no smoothing is performed. }
 #' \item{kFoldMuCov}{The number of folds to be used for mean and covariance smoothing. Default: 10}
 #' \item{methodMuCovEst}{The method to estimate the mean and covariance in the case of dense functional data; 'cross-sectional', 'smooth' - default: 'cross-sectional'}
 #' \item{numBins}{The number of bins to bin the data into; positive integer > 10, default: NULL}
diff --git a/R/Lwls2D.R b/R/Lwls2D.R
index dd28e2f..abb9870 100644
--- a/R/Lwls2D.R
+++ b/R/Lwls2D.R
@@ -2,7 +2,7 @@
 #'
 #' Two dimensional local weighted least squares smoother. Only local linear smoother for estimating the original curve is available (no higher order, no derivative). 
 #' @param bw A scalar or a vector of length 2 specifying the bandwidth.
-#' @param kern Kernel used: 'gauss', 'rect', 'gausvar', 'epan' (default), 'quar'.
+#' @param kern Kernel used: 'gauss', 'rect', 'gausvar', 'epan' (default), 'quar', "triangular", "triweight", "tricube", "cosine", "logistic", "sigmoid", "silverman".
 #' @param xin An n by 2 data frame or matrix of x-coordinate.
 #' @param yin A vector of y-coordinate.
 #' @param win A vector of weights on the observations. 
diff --git a/R/Lwls2DDeriv.R b/R/Lwls2DDeriv.R
index 4cce5f2..0ec6954 100644
--- a/R/Lwls2DDeriv.R
+++ b/R/Lwls2DDeriv.R
@@ -2,7 +2,7 @@
 #'
 #' Two dimensional local weighted least squares smoother. Only a local linear smoother for estimating the original curve is available (no higher order)
 #' @param bw A scalar or a vector of length 2 specifying the bandwidth.
-#' @param kern Kernel used: 'gauss', 'rect', 'gausvar', 'epan' (default), 'quar'.
+#' @param kern Kernel used: 'gauss', 'rect', 'gausvar', 'epan' (default), 'quar', "triangular", "triweight", "tricube", "cosine", "logistic", "sigmoid", "silverman".
 #' @param xin An n by 2 data frame or matrix of x-coordinate.
 #' @param yin A vector of y-coordinate.
 #' @param win A vector of weights on the observations. 
diff --git a/R/SetOptions.R b/R/SetOptions.R
index 6759dd1..7411d96 100755
--- a/R/SetOptions.R
+++ b/R/SetOptions.R
@@ -186,7 +186,8 @@ SetOptions = function(y, t, optns){
       kernel = "gauss";  # kernel: Gaussian
     }
   }
-  kernNames = c("rect", "gauss", "epan", "gausvar", "quar");
+  kernNames = c("rect", "gauss", "epan", "gausvar", "quar", "triangular", 
+                "triweight", "tricube", "cosine", "logistic", "sigmoid", "silverman");
   if(!(kernel %in% kernNames)){ # Check suitability of kernel
     message(paste('kernel', kernel, 'is unrecognizable! Reset to automatic selection now!\n')); 
     kernel = NULL; 
diff --git a/src/CPPlwls1d.cpp b/src/CPPlwls1d.cpp
index 0970d9f..3bc9c6e 100644
--- a/src/CPPlwls1d.cpp
+++ b/src/CPPlwls1d.cpp
@@ -49,6 +49,13 @@ Eigen::VectorXd CPPlwls1d( const double & bw, const std::string kernel_type, con
   possibleKernels["epan"]    = 1;   possibleKernels["rect"]    = 2;
   possibleKernels["gauss"]   = 3;   possibleKernels["gausvar"] = 4;
   possibleKernels["quar"]    = 5;
+  possibleKernels["triangular"] = 6; // Added triangular kernel
+  possibleKernels["triweight"] = 7; // Added triweight kernel
+  possibleKernels["tricube"] = 8; // Added tricube kernel
+  possibleKernels["cosine"] = 9; // Added cosine kernel
+  possibleKernels["logistic"] = 10; // Added logistic kernel
+  possibleKernels["sigmoid"] = 11; // Added sigmoid kernel
+  possibleKernels["silverman"] = 12; // Added silverman kernel
 
   // If the kernel_type key exists set KernelName appropriately
   int KernelName = 0;
@@ -84,7 +91,7 @@ Eigen::VectorXd CPPlwls1d( const double & bw, const std::string kernel_type, con
     const double* upper ;
 
     //if the kernel is not Gaussian
-    if ( KernelName != 3 && KernelName != 4) {
+    if ( KernelName != 3 && KernelName != 4 && KernelName != 10 && KernelName != 11 && KernelName != 12) {
       //construct listX as vectors / size is unknown originally
       // for (unsigned int y = 0; y != nXGrid; ++y){ if ( std::fabs( xout(i) - xin(y) ) <= bw ) { indx.push_back(y); }  }
       // Get iterator pointing to the first element which is not less than xou(u)
@@ -141,6 +148,27 @@ Eigen::VectorXd CPPlwls1d( const double & bw, const std::string kernel_type, con
         temp = (lw.array()) *
                ((1.-llx.array().pow(2)).array().pow(2)).array() * (15./16.);
         break;
+    case 6 : // truangular
+      temp = (lw.array()) * (1.-llx.array().abs()).array();
+      break;
+    case 7 : // triweight
+      temp = (lw.array()) * (1.-llx.array().pow(2)).array().pow(3) * (35./32.);
+      break;
+    case 8 : // tricube
+      temp = (lw.array()) * (1.-llx.array().abs().pow(3)).array().pow(3) * (70./81.);
+      break;
+    case 9 : // cosine
+      temp = (lw.array()) * M_PI / 4.0 * (llx.array() * M_PI / 2.0).cos().array();
+      break;
+    case 10 : // logistic
+      temp = (lw.array()) / ((-llx.array()).exp() + 2. + llx.array().exp()).array();
+      break;
+    case 11 : // sigmoid
+      temp = (lw.array()) / ((-llx.array()).exp() + llx.array().exp()).array() * 2.0 / M_PI;
+      break;
+    case 12 : // silverman
+      temp = (lw.array()) * 0.5 * (-llx.array().abs() / sqrt(2)).exp() * (M_PI / 4.0 + llx.array().abs() * sqrt(2)).sin().array();
+      break;
     }
 
     if(nder >= indxSize){
diff --git a/src/Rmullwlsk.cpp b/src/Rmullwlsk.cpp
index bb4c5f5..6131024 100644
--- a/src/Rmullwlsk.cpp
+++ b/src/Rmullwlsk.cpp
@@ -23,6 +23,13 @@ Eigen::MatrixXd Rmullwlsk( const Eigen::Map<Eigen::VectorXd> & bw, const std::st
   possibleKernels["epan"]    = 1;   possibleKernels["rect"]    = 2;
   possibleKernels["gauss"]   = 3;   possibleKernels["gausvar"] = 4; 
   possibleKernels["quar"]    = 5; 
+  possibleKernels["triangular"] = 6; // Added triangular kernel
+  possibleKernels["triweight"] = 7; // Added triweight kernel
+  possibleKernels["tricube"] = 8; // Added tricube kernel
+  possibleKernels["cosine"] = 9; // Added cosine kernel
+  possibleKernels["logistic"] = 10; // Added logistic kernel
+  possibleKernels["sigmoid"] = 11; // Added sigmoid kernel
+  possibleKernels["silverman"] = 12; // Added silverman kernel
    
   // The following test is here for completeness, we mightwant to move it up a 
   // level (in the wrapper) in the future. 
@@ -59,7 +66,7 @@ Eigen::MatrixXd Rmullwlsk( const Eigen::Map<Eigen::VectorXd> & bw, const std::st
       //locating local window (LOL) (bad joke)
       std::vector <unsigned int> indx; 
       //if the kernel is not Gaussian
-      if ( KernelName != 3) { 
+      if ( KernelName != 3 && KernelName != 10 && KernelName != 11 && KernelName != 12) { 
         //construct listX as vectors / size is unknown originally
         for (unsigned int y = 0; y != tPairs.cols(); y++){ 
           if ( std::abs( tPairs(0,y) - xgrid(j) ) <= (bw(0)+ bufSmall) && std::abs( tPairs(1,y) - ygrid(i) ) <= (bw(1)+ bufSmall) ) {
@@ -132,6 +139,33 @@ Eigen::MatrixXd Rmullwlsk( const Eigen::Map<Eigen::VectorXd> & bw, const std::st
                      ((1.-llx.row(0).array().pow(2)).array().pow(2)).array() *
                      ((1.-llx.row(1).array().pow(2)).array().pow(2)).array() * (225./256.);
             break;
+        case 6 : // truangular
+          temp = (lw.transpose().array()) * (1.-llx.row(0).array().abs()).array() * (1.-llx.row(1).array().abs()).array();
+          break;
+        case 7 : // triweight
+          temp = (lw.transpose().array()) * (1.-llx.row(0).array().pow(2)).array().pow(3) * (35./32.) *
+            (1.-llx.row(1).array().pow(2)).array().pow(3) * (35./32.);
+          break;
+        case 8 : // tricube
+          temp = (lw.transpose().array()) * (1.-llx.row(0).array().abs().pow(3)).array().pow(3) * (70./81.) *
+            (1.-llx.row(1).array().abs().pow(3)).array().pow(3) * (70./81.);
+          break;
+        case 9 : // cosine
+          temp = (lw.transpose().array()) * M_PI / 4.0 * (llx.row(0).array() * M_PI / 2.0).cos().array() *
+            M_PI / 4.0 * (llx.row(1).array() * M_PI / 2.0).cos().array();
+          break;
+        case 10 : // logistic
+          temp = (lw.transpose().array()) / ((-llx.row(0).array()).exp() + 2. + llx.row(0).array().exp()).array() /
+            ((-llx.row(1).array()).exp() + 2. + llx.row(1).array().exp()).array();
+          break;
+        case 11 : // sigmoid
+          temp = (lw.transpose().array()) / ((-llx.row(0).array()).exp() + llx.row(0).array().exp()).array() * 2.0 / M_PI /
+            ((-llx.row(1).array()).exp() + llx.row(1).array().exp()).array() * 2.0 / M_PI;
+          break;
+        case 12 : // silverman
+          temp = (lw.transpose().array()) * 0.5 * (-llx.row(0).array().abs() / sqrt(2)).exp() * (M_PI / 4.0 + llx.row(0).array().abs() * sqrt(2)).sin().array() *
+            0.5 * (-llx.row(1).array().abs() / sqrt(2)).exp() * (M_PI / 4.0 + llx.row(1).array().abs() * sqrt(2)).sin().array();
+          break;
         } 
 
         // make the design matrix
diff --git a/src/RmullwlskCC.cpp b/src/RmullwlskCC.cpp
index 4e11e25..6fa7400 100644
--- a/src/RmullwlskCC.cpp
+++ b/src/RmullwlskCC.cpp
@@ -23,6 +23,13 @@ Eigen::MatrixXd RmullwlskCC( const Eigen::Map<Eigen::VectorXd> & bw, const std::
   possibleKernels["epan"]    = 1;   possibleKernels["rect"]    = 2;
   possibleKernels["gauss"]   = 3;   possibleKernels["gausvar"] = 4; 
   possibleKernels["quar"]    = 5; 
+  possibleKernels["triangular"] = 6; // Added triangular kernel
+  possibleKernels["triweight"] = 7; // Added triweight kernel
+  possibleKernels["tricube"] = 8; // Added tricube kernel
+  possibleKernels["cosine"] = 9; // Added cosine kernel
+  possibleKernels["logistic"] = 10; // Added logistic kernel
+  possibleKernels["sigmoid"] = 11; // Added sigmoid kernel
+  possibleKernels["silverman"] = 12; // Added silverman kernel
    
   // The following test is here for completeness, we mightwant to move it up a 
   // level (in the wrapper) in the future. 
@@ -59,7 +66,7 @@ Eigen::MatrixXd RmullwlskCC( const Eigen::Map<Eigen::VectorXd> & bw, const std::
       //locating local window (LOL) (bad joke)
       std::vector <unsigned int> indx; 
       //if the kernel is not Gaussian
-      if ( KernelName != 3) { 
+      if ( KernelName != 3 && KernelName != 10 && KernelName != 11 && KernelName != 12) { 
         //construct listX as vectors / size is unknown originally
         for (unsigned int y = 0; y != tPairs.cols(); y++){ 
           if (  (std::abs( tPairs(0,y) - xgrid(i) ) <= (bw(0)+ bufSmall ) && std::abs( tPairs(1,y) - ygrid(j) ) <= (bw(1)+ bufSmall)) ) {
@@ -135,6 +142,33 @@ Eigen::MatrixXd RmullwlskCC( const Eigen::Map<Eigen::VectorXd> & bw, const std::
                      ((1.-llx.row(0).array().pow(2)).array().pow(2)).array() *
                      ((1.-llx.row(1).array().pow(2)).array().pow(2)).array() * (225./256.);
             break;
+        case 6 : // truangular
+          temp = (lw.transpose().array()) * (1.-llx.row(0).array().abs()).array() * (1.-llx.row(1).array().abs()).array();
+          break;
+        case 7 : // triweight
+          temp = (lw.transpose().array()) * (1.-llx.row(0).array().pow(2)).array().pow(3) * (35./32.) *
+            (1.-llx.row(1).array().pow(2)).array().pow(3) * (35./32.);
+          break;
+        case 8 : // tricube
+          temp = (lw.transpose().array()) * (1.-llx.row(0).array().abs().pow(3)).array().pow(3) * (70./81.) *
+            (1.-llx.row(1).array().abs().pow(3)).array().pow(3) * (70./81.);
+          break;
+        case 9 : // cosine
+          temp = (lw.transpose().array()) * M_PI / 4.0 * (llx.row(0).array() * M_PI / 2.0).cos().array() *
+            M_PI / 4.0 * (llx.row(1).array() * M_PI / 2.0).cos().array();
+          break;
+        case 10 : // logistic
+          temp = (lw.transpose().array()) / ((-llx.row(0).array()).exp() + 2. + llx.row(0).array().exp()).array() /
+            ((-llx.row(1).array()).exp() + 2. + llx.row(1).array().exp()).array();
+          break;
+        case 11 : // sigmoid
+          temp = (lw.transpose().array()) / ((-llx.row(0).array()).exp() + llx.row(0).array().exp()).array() * 2.0 / M_PI /
+            ((-llx.row(1).array()).exp() + llx.row(1).array().exp()).array() * 2.0 / M_PI;
+          break;
+        case 12 : // silverman
+          temp = (lw.transpose().array()) * 0.5 * (-llx.row(0).array().abs() / sqrt(2)).exp() * (M_PI / 4.0 + llx.row(0).array().abs() * sqrt(2)).sin().array() *
+            0.5 * (-llx.row(1).array().abs() / sqrt(2)).exp() * (M_PI / 4.0 + llx.row(1).array().abs() * sqrt(2)).sin().array();
+          break;
         } 
 
         // make the design matrix
diff --git a/src/RmullwlskCCsort2.cpp b/src/RmullwlskCCsort2.cpp
index 7733844..5d20424 100644
--- a/src/RmullwlskCCsort2.cpp
+++ b/src/RmullwlskCCsort2.cpp
@@ -28,6 +28,13 @@ Eigen::MatrixXd RmullwlskCCsort2( const Eigen::Map<Eigen::VectorXd> & bw, const
   possibleKernels["epan"]    = 1;   possibleKernels["rect"]    = 2;
   possibleKernels["gauss"]   = 3;   possibleKernels["gausvar"] = 4; 
   possibleKernels["quar"]    = 5; 
+  possibleKernels["triangular"] = 6; // Added triangular kernel
+  possibleKernels["triweight"] = 7; // Added triweight kernel
+  possibleKernels["tricube"] = 8; // Added tricube kernel
+  possibleKernels["cosine"] = 9; // Added cosine kernel
+  possibleKernels["logistic"] = 10; // Added logistic kernel
+  possibleKernels["sigmoid"] = 11; // Added sigmoid kernel
+  possibleKernels["silverman"] = 12; // Added silverman kernel
    
   // The following test is here for completeness, we mightwant to move it up a 
   // level (in the wrapper) in the future. 
@@ -89,7 +96,7 @@ Eigen::MatrixXd RmullwlskCCsort2( const Eigen::Map<Eigen::VectorXd> & bw, const
       std::vector <unsigned int> indx; 
       
       //if the kernel is not Gaussian
-      if ( KernelName != 3) { 
+      if ( KernelName != 3 && KernelName != 10 && KernelName != 11 && KernelName != 12) { 
       // Search the lower and upper bounds increasingly.
         ylIt = std::lower_bound(ylIt, yval.end(), valIndPair(yl, 0), compPair);
         yuIt = std::upper_bound(yuIt, yval.end(), valIndPair(yu, 0), compPair);
@@ -171,6 +178,33 @@ Eigen::MatrixXd RmullwlskCCsort2( const Eigen::Map<Eigen::VectorXd> & bw, const
                      ((1.-llx.row(0).array().pow(2)).array().pow(2)).array() *
                      ((1.-llx.row(1).array().pow(2)).array().pow(2)).array() * (225./256.);
             break;
+        case 6 : // truangular
+          temp = (lw.transpose().array()) * (1.-llx.row(0).array().abs()).array() * (1.-llx.row(1).array().abs()).array();
+          break;
+        case 7 : // triweight
+          temp = (lw.transpose().array()) * (1.-llx.row(0).array().pow(2)).array().pow(3) * (35./32.) *
+            (1.-llx.row(1).array().pow(2)).array().pow(3) * (35./32.);
+          break;
+        case 8 : // tricube
+          temp = (lw.transpose().array()) * (1.-llx.row(0).array().abs().pow(3)).array().pow(3) * (70./81.) *
+            (1.-llx.row(1).array().abs().pow(3)).array().pow(3) * (70./81.);
+          break;
+        case 9 : // cosine
+          temp = (lw.transpose().array()) * M_PI / 4.0 * (llx.row(0).array() * M_PI / 2.0).cos().array() *
+            M_PI / 4.0 * (llx.row(1).array() * M_PI / 2.0).cos().array();
+          break;
+        case 10 : // logistic
+          temp = (lw.transpose().array()) / ((-llx.row(0).array()).exp() + 2. + llx.row(0).array().exp()).array() /
+            ((-llx.row(1).array()).exp() + 2. + llx.row(1).array().exp()).array();
+          break;
+        case 11 : // sigmoid
+          temp = (lw.transpose().array()) / ((-llx.row(0).array()).exp() + llx.row(0).array().exp()).array() * 2.0 / M_PI /
+            ((-llx.row(1).array()).exp() + llx.row(1).array().exp()).array() * 2.0 / M_PI;
+          break;
+        case 12 : // silverman
+          temp = (lw.transpose().array()) * 0.5 * (-llx.row(0).array().abs() / sqrt(2)).exp() * (M_PI / 4.0 + llx.row(0).array().abs() * sqrt(2)).sin().array() *
+            0.5 * (-llx.row(1).array().abs() / sqrt(2)).exp() * (M_PI / 4.0 + llx.row(1).array().abs() * sqrt(2)).sin().array();
+          break;
         } 
 
         // make the design matrix
diff --git a/src/RmullwlskUniversal.cpp b/src/RmullwlskUniversal.cpp
index 5000524..a4224e2 100644
--- a/src/RmullwlskUniversal.cpp
+++ b/src/RmullwlskUniversal.cpp
@@ -29,6 +29,13 @@ Eigen::MatrixXd RmullwlskUniversal( const Eigen::Map<Eigen::VectorXd> & bw, cons
   possibleKernels["epan"]    = 1;   possibleKernels["rect"]    = 2;
   possibleKernels["gauss"]   = 3;   possibleKernels["gausvar"] = 4; 
   possibleKernels["quar"]    = 5; 
+  possibleKernels["triangular"] = 6; // Added triangular kernel
+  possibleKernels["triweight"] = 7; // Added triweight kernel
+  possibleKernels["tricube"] = 8; // Added tricube kernel
+  possibleKernels["cosine"] = 9; // Added cosine kernel
+  possibleKernels["logistic"] = 10; // Added logistic kernel
+  possibleKernels["sigmoid"] = 11; // Added sigmoid kernel
+  possibleKernels["silverman"] = 12; // Added silverman kernel
    
   // The following test is here for completeness, we mightwant to move it up a 
   // level (in the wrapper) in the future. 
@@ -90,7 +97,7 @@ Eigen::MatrixXd RmullwlskUniversal( const Eigen::Map<Eigen::VectorXd> & bw, cons
       std::vector <unsigned int> indx; 
       
       //if the kernel is not Gaussian
-      if ( KernelName != 3) { 
+      if ( KernelName != 3 && KernelName != 10 && KernelName != 11 && KernelName != 12) { 
       // Search the lower and upper bounds increasingly.
         ylIt = std::lower_bound(ylIt, yval.end(), valIndPair(yl, 0), comparePair);
         yuIt = std::upper_bound(yuIt, yval.end(), valIndPair(yu, 0), comparePair);
@@ -169,6 +176,33 @@ Eigen::MatrixXd RmullwlskUniversal( const Eigen::Map<Eigen::VectorXd> & bw, cons
                      ((1.-llx.row(0).array().pow(2)).array().pow(2)).array() *
                      ((1.-llx.row(1).array().pow(2)).array().pow(2)).array() * (225./256.);
             break;
+        case 6 : // truangular
+          temp = (lw.transpose().array()) * (1.-llx.row(0).array().abs()).array() * (1.-llx.row(1).array().abs()).array();
+          break;
+        case 7 : // triweight
+          temp = (lw.transpose().array()) * (1.-llx.row(0).array().pow(2)).array().pow(3) * (35./32.) *
+            (1.-llx.row(1).array().pow(2)).array().pow(3) * (35./32.);
+          break;
+        case 8 : // tricube
+          temp = (lw.transpose().array()) * (1.-llx.row(0).array().abs().pow(3)).array().pow(3) * (70./81.) *
+            (1.-llx.row(1).array().abs().pow(3)).array().pow(3) * (70./81.);
+          break;
+        case 9 : // cosine
+          temp = (lw.transpose().array()) * M_PI / 4.0 * (llx.row(0).array() * M_PI / 2.0).cos().array() *
+            M_PI / 4.0 * (llx.row(1).array() * M_PI / 2.0).cos().array();
+          break;
+        case 10 : // logistic
+          temp = (lw.transpose().array()) / ((-llx.row(0).array()).exp() + 2. + llx.row(0).array().exp()).array() /
+            ((-llx.row(1).array()).exp() + 2. + llx.row(1).array().exp()).array();
+          break;
+        case 11 : // sigmoid
+          temp = (lw.transpose().array()) / ((-llx.row(0).array()).exp() + llx.row(0).array().exp()).array() * 2.0 / M_PI /
+            ((-llx.row(1).array()).exp() + llx.row(1).array().exp()).array() * 2.0 / M_PI;
+          break;
+        case 12 : // silverman
+          temp = (lw.transpose().array()) * 0.5 * (-llx.row(0).array().abs() / sqrt(2)).exp() * (M_PI / 4.0 + llx.row(0).array().abs() * sqrt(2)).sin().array() *
+            0.5 * (-llx.row(1).array().abs() / sqrt(2)).exp() * (M_PI / 4.0 + llx.row(1).array().abs() * sqrt(2)).sin().array();
+          break;
         } 
 
         // make the design matrix
diff --git a/src/RmullwlskUniversalDeriv.cpp b/src/RmullwlskUniversalDeriv.cpp
index 9a4f525..29fdde6 100644
--- a/src/RmullwlskUniversalDeriv.cpp
+++ b/src/RmullwlskUniversalDeriv.cpp
@@ -45,6 +45,13 @@ Eigen::MatrixXd RmullwlskUniversalDeriv(
   possibleKernels["epan"]    = 1;   possibleKernels["rect"]    = 2;
   possibleKernels["gauss"]   = 3;   possibleKernels["gausvar"] = 4; 
   possibleKernels["quar"]    = 5; 
+  possibleKernels["triangular"] = 6; // Added triangular kernel
+  possibleKernels["triweight"] = 7; // Added triweight kernel
+  possibleKernels["tricube"] = 8; // Added tricube kernel
+  possibleKernels["cosine"] = 9; // Added cosine kernel
+  possibleKernels["logistic"] = 10; // Added logistic kernel
+  possibleKernels["sigmoid"] = 11; // Added sigmoid kernel
+  possibleKernels["silverman"] = 12; // Added silverman kernel
    
   // The following test is here for completeness, we mightwant to move it up a 
   // level (in the wrapper) in the future. 
@@ -118,7 +125,7 @@ Eigen::MatrixXd RmullwlskUniversalDeriv(
       std::vector <unsigned int> indx; 
       
       //if the kernel is not Gaussian
-      if ( KernelName != 3) { 
+      if ( KernelName != 3 && KernelName != 10 && KernelName != 11 && KernelName != 12) { 
       // Search the lower and upper bounds increasingly.
         ylIt = std::lower_bound(ylIt, yval.end(), valIndPair(yl, 0), comparePair);
         yuIt = std::upper_bound(yuIt, yval.end(), valIndPair(yu, 0), comparePair);
@@ -197,6 +204,33 @@ Eigen::MatrixXd RmullwlskUniversalDeriv(
                      ((1.-llx.row(0).array().pow(2)).array().pow(2)).array() *
                      ((1.-llx.row(1).array().pow(2)).array().pow(2)).array() * (225./256.);
             break;
+        case 6 : // truangular
+          temp = (lw.transpose().array()) * (1.-llx.row(0).array().abs()).array() * (1.-llx.row(1).array().abs()).array();
+          break;
+        case 7 : // triweight
+          temp = (lw.transpose().array()) * (1.-llx.row(0).array().pow(2)).array().pow(3) * (35./32.) *
+            (1.-llx.row(1).array().pow(2)).array().pow(3) * (35./32.);
+          break;
+        case 8 : // tricube
+          temp = (lw.transpose().array()) * (1.-llx.row(0).array().abs().pow(3)).array().pow(3) * (70./81.) *
+            (1.-llx.row(1).array().abs().pow(3)).array().pow(3) * (70./81.);
+          break;
+        case 9 : // cosine
+          temp = (lw.transpose().array()) * M_PI / 4.0 * (llx.row(0).array() * M_PI / 2.0).cos().array() *
+            M_PI / 4.0 * (llx.row(1).array() * M_PI / 2.0).cos().array();
+          break;
+        case 10 : // logistic
+          temp = (lw.transpose().array()) / ((-llx.row(0).array()).exp() + 2. + llx.row(0).array().exp()).array() /
+            ((-llx.row(1).array()).exp() + 2. + llx.row(1).array().exp()).array();
+          break;
+        case 11 : // sigmoid
+          temp = (lw.transpose().array()) / ((-llx.row(0).array()).exp() + llx.row(0).array().exp()).array() * 2.0 / M_PI /
+            ((-llx.row(1).array()).exp() + llx.row(1).array().exp()).array() * 2.0 / M_PI;
+          break;
+        case 12 : // silverman
+          temp = (lw.transpose().array()) * 0.5 * (-llx.row(0).array().abs() / sqrt(2)).exp() * (M_PI / 4.0 + llx.row(0).array().abs() * sqrt(2)).sin().array() *
+            0.5 * (-llx.row(1).array().abs() / sqrt(2)).exp() * (M_PI / 4.0 + llx.row(1).array().abs() * sqrt(2)).sin().array();
+          break;
         } 
 
         // make the design matrix
diff --git a/src/RrotatedMullwlsk.cpp b/src/RrotatedMullwlsk.cpp
index 7b244f6..858d8f0 100644
--- a/src/RrotatedMullwlsk.cpp
+++ b/src/RrotatedMullwlsk.cpp
@@ -21,6 +21,13 @@ Eigen::VectorXd Rrotatedmullwlsk( const Eigen::Map<Eigen::VectorXd> & bw, const
   possibleKernels["epan"]    = 1;   possibleKernels["rect"]    = 2;
   possibleKernels["gauss"]   = 3;   possibleKernels["gausvar"] = 4; 
   possibleKernels["quar"]    = 5; 
+  possibleKernels["triangular"] = 6; // Added triangular kernel
+  possibleKernels["triweight"] = 7; // Added triweight kernel
+  possibleKernels["tricube"] = 8; // Added tricube kernel
+  possibleKernels["cosine"] = 9; // Added cosine kernel
+  possibleKernels["logistic"] = 10; // Added logistic kernel
+  possibleKernels["sigmoid"] = 11; // Added sigmoid kernel
+  possibleKernels["silverman"] = 12; // Added silverman kernel
    
   // The following test is here for completeness, we mightwant to move it up a 
   // level (in the wrapper) in the future. 
@@ -60,7 +67,7 @@ Eigen::VectorXd Rrotatedmullwlsk( const Eigen::Map<Eigen::VectorXd> & bw, const
     //locating local window (LOL) (bad joke)
     std::vector <unsigned int> indx; 
     //if the kernel is not Gaussian or Gaussian-like
-    if ( KernelName != 3 && KernelName != 4 ) { 
+    if ( KernelName != 3 && KernelName != 4 && KernelName != 10 && KernelName != 11 && KernelName != 12) { 
       //construct listX as vectors / size is unknown originally
       std::vector <unsigned int> list1, list2; 
       for (unsigned int y = 0; y != tPairs.cols(); y++){ 
@@ -126,6 +133,33 @@ Eigen::VectorXd Rrotatedmullwlsk( const Eigen::Map<Eigen::VectorXd> & bw, const
                  ((1.-llx.row(0).array().pow(2)).array().pow(2)).array() *
                  ((1.-llx.row(1).array().pow(2)).array().pow(2)).array() * (225./256.);
           break;
+      case 6 : // truangular
+        temp = (lw.transpose().array()) * (1.-llx.row(0).array().abs()).array() * (1.-llx.row(1).array().abs()).array();
+        break;
+      case 7 : // triweight
+        temp = (lw.transpose().array()) * (1.-llx.row(0).array().pow(2)).array().pow(3) * (35./32.) *
+          (1.-llx.row(1).array().pow(2)).array().pow(3) * (35./32.);
+        break;
+      case 8 : // tricube
+        temp = (lw.transpose().array()) * (1.-llx.row(0).array().abs().pow(3)).array().pow(3) * (70./81.) *
+          (1.-llx.row(1).array().abs().pow(3)).array().pow(3) * (70./81.);
+        break;
+      case 9 : // cosine
+        temp = (lw.transpose().array()) * M_PI / 4.0 * (llx.row(0).array() * M_PI / 2.0).cos().array() *
+          M_PI / 4.0 * (llx.row(1).array() * M_PI / 2.0).cos().array();
+        break;
+      case 10 : // logistic
+        temp = (lw.transpose().array()) / ((-llx.row(0).array()).exp() + 2. + llx.row(0).array().exp()).array() /
+          ((-llx.row(1).array()).exp() + 2. + llx.row(1).array().exp()).array();
+        break;
+      case 11 : // sigmoid
+        temp = (lw.transpose().array()) / ((-llx.row(0).array()).exp() + llx.row(0).array().exp()).array() * 2.0 / M_PI /
+          ((-llx.row(1).array()).exp() + llx.row(1).array().exp()).array() * 2.0 / M_PI;
+        break;
+      case 12 : // silverman
+        temp = (lw.transpose().array()) * 0.5 * (-llx.row(0).array().abs() / sqrt(2)).exp() * (M_PI / 4.0 + llx.row(0).array().abs() * sqrt(2)).sin().array() *
+          0.5 * (-llx.row(1).array().abs() / sqrt(2)).exp() * (M_PI / 4.0 + llx.row(1).array().abs() * sqrt(2)).sin().array();
+        break;
       } 
       
       // make the design matrix