diff --git a/DESCRIPTION b/DESCRIPTION
index bb9dc670022f18a2782294d69fb89dca6f9ce543..fbdecc7d31d17abdbb0526ac101d8c922e7d714a 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -14,6 +14,7 @@ Imports:
     Rcpp (>= 0.12.18),
     R6 (>= 2.2.2),
     splines (>= 3.1.1),
+    pbs,
     digest,
 LinkingTo: Rcpp, RcppArmadillo
 Suggests:
diff --git a/R/as.data.list.R b/R/as.data.list.R
index 5645d9f0032b64db173a92bbd02ee7d67843de4a..fc685a0610e9ce69a821f6b6a64837e01fbaba4b 100644
--- a/R/as.data.list.R
+++ b/R/as.data.list.R
@@ -48,7 +48,7 @@ as.data.list <- function(object){
 #' as.data.list(X)
 #'
 #' # Convert a dataframe with time, forecast and an observed variable
-#' X <- data.frame(t=1:10, x.k1=1:10, x.k2=10:1, y.obs=1:10, y.k1=1:10, y.k2=1:10)
+#' X <- data.frame(t=1:10, x.k1=1:10, x.k2=10:1, yobs=1:10, y.k1=1:10, y.k2=1:10)
 #' as.data.list(X)
 #'
 #' # Can be converted back and forth
diff --git a/R/bspline.R b/R/bspline.R
index d3529410bd6f6716a61f888910e94e6469605726..39103ed6961806577a5714b16f33bb4b5818dc8f 100644
--- a/R/bspline.R
+++ b/R/bspline.R
@@ -17,44 +17,56 @@
 #' 
 #' @param X data.frame (as part of data.list) with horizons as columns named \code{kxx} (i.e. one for each horizon)
 #' @param Boundary.knots The value is NA: then the boundaries are set to the range of each horizons (columns in X). See \code{?splines::bs}
-#' @param intercept Default value is TRUE: in an onlineforecast model there is no intercept per default (set by \code{one()}). It's default to FALSE in \code{?splines::bs}.
+#' @param intercept See \code{?splines::bs}.
 #' @param df See \code{?splines::bs}
 #' @param knots See \code{?splines::bs}
 #' @param degree See \code{?splines::bs}
 #' @param bknots Is just a short for Boundary.knots and replace Boundary.knots (if Boundary.knots is not given)
+#' @param periodic Default FALSE. If TRUE, then \code{pbs::pbs} is called and periodic splines are generated.
 #' @return List of data frames with the computed base splines, each with columns for the same horizons as in X
 #' @rdname bs
 #' @examples
 #'
 #' # How to make a diurnal curve using splines
-#' # Select first 54 hours from the load data
-#' D <- subset(Dbuilding, 1:54, kseq=1:4)
-#' # Make the hour of the day as a forecast input
-#' D$tday <- make_tday(D$t, kseq=1:4)
-#' D$tday
+#'  # Select first 54 hours from the load data
+#'  D <- subset(Dbuilding, 1:76, kseq=1:4)
+#'  # Make the hour of the day as a forecast input
+#'  D$tday <- make_tday(D$t, kseq=1:4)
+#'  D$tday
 #' 
-#' # Calculate the base splines for each column in tday
-#' L <- bspline(D$tday)
+#'  # Calculate the base splines for each column in tday
+#'  L <- bspline(D$tday)
 #'
-#' # Now L holds a data.frame for each base spline
-#' str(L)
-#' # Hence this will result in four inputs for the regression model
+#'  # Now L holds a data.frame for each base spline
+#'  str(L)
+#'  # Hence this will result in four inputs for the regression model
 #'
-#' # Plot (note that the splines period starts at tday=0)
+#'  # Plot (note that the splines period starts at tday=0)
+#'  plot(D$t, L$bs1$k1, type="s")
+#'  for(i in 2:length(L)){
+#'    lines(D$t, L[[i]]$k1, col=i, type="s")
+#'  }
+#'
+#'  # In a model formulation it will be:
+#'  model <- forecastmodel$new()
+#'  model$add_inputs(mutday = "bspline(tday)")
+#'  # Such that at the transform stage will give the same as above
+#'  model$transform_data(D)
+#'
+#' # Periodic splines are useful for modelling a diurnal harmonical functions
+#' L <- bspline(D$tday, bknots=c(0,24), df=4, periodic=TRUE)
+#' # or
+#' L <- pbspline(D$tday, bknots=c(0,24), df=4)
+#' # Note, how it has to have high enough df, else it generates an error
+#'
+#' # Plot
 #' plot(D$t, L$bs1$k1, type="s")
 #' for(i in 2:length(L)){
-#'   lines(D$t, L[[i]]$k1, col=i, type="s")
+#'     lines(D$t, L[[i]]$k1, col=i, type="s")
 #' }
 #'
-#' # In a model formulation it will be:
-#' model <- forecastmodel$new()
-#' model$add_inputs(mutday = "bspline(tday)")
-#' # Such that at the transform stage will give the same as above
-#' model$transform_data(D)
-#'
-#'
 #' @export
-bspline <- function(X, Boundary.knots = NA, intercept = TRUE, df = NULL, knots = NULL, degree = 3, bknots = NA) {
+bspline <- function(X, Boundary.knots = NA, intercept = FALSE, df = NULL, knots = NULL, degree = 3, bknots = NA, periodic = FALSE) {
     # bknots is just a short for Boundary.knots and replace if Boundary.knots are not given.
     if(is.na(Boundary.knots)){
         Boundary.knots <- bknots
@@ -64,7 +76,7 @@ bspline <- function(X, Boundary.knots = NA, intercept = TRUE, df = NULL, knots =
         # Call again for each element
         val <- lapply(1:length(X), function(i) {
             bspline(X[[i]], df = df, knots = knots, degree = degree, intercept = intercept, 
-                    Boundary.knots = Boundary.knots)
+                    Boundary.knots = Boundary.knots, periodic = periodic)
         })
         nams(val) <- nams(X)
         return(val)
@@ -77,8 +89,21 @@ bspline <- function(X, Boundary.knots = NA, intercept = TRUE, df = NULL, knots =
       if (is.na(Boundary.knots[1])) {
           Boundary.knots <- range(X[, i], na.rm=TRUE)
       }
-      spls <- splines::bs(X[, i], Boundary.knots = Boundary.knots, degree = degree, df = df,
-                          knots = knots, intercept = intercept)
+      if(periodic){
+          # Periodic splines (annoyingly it use 'missing()' to conclude if the argument was given, not 'is.null()' as bs does)
+          if(is.null(knots)){
+              # Call without knots
+              spls <- pbs::pbs(X[, i], Boundary.knots = Boundary.knots, degree = degree, df = df,
+                               intercept = intercept)
+          }else{
+              # Call without df
+              spls <- pbs::pbs(X[, i], Boundary.knots = Boundary.knots, degree = degree,
+                               knots = knots, intercept = intercept)
+          }
+      }else{
+          spls <- splines::bs(X[, i], Boundary.knots = Boundary.knots, degree = degree, df = df,
+                              knots = knots, intercept = intercept)
+      }
       return(spls)
     })
     # Now we have a bs value in L for each horizon
@@ -95,3 +120,24 @@ bspline <- function(X, Boundary.knots = NA, intercept = TRUE, df = NULL, knots =
     nams(L) <- pst("bs", 1:length(L))
     return(L)
 }
+
+#' Wrapper for \code{bspline} with \code{periodic=TRUE}
+#'
+#' Simply a wrapper.
+#' 
+#' @title Wrapper for \code{bspline} with \code{periodic=TRUE}
+#' @param X see \code{?bspline}
+#' @param Boundary.knots see \code{?bspline}
+#' @param intercept see \code{?bspline}
+#' @param df see \code{?bspline}
+#' @param knots see \code{?bspline}
+#' @param degree see \code{?bspline}
+#' @param bknots see \code{?bspline}
+#' @family Transform stage functions
+#'
+#' @export
+pbspline <- function(X, Boundary.knots = NA, intercept = FALSE, df = NULL, knots = NULL, degree = 3, bknots = NA){
+    bspline(X, Boundary.knots = Boundary.knots, degree = degree, df = df,
+            knots = knots, intercept = intercept, bknots = bknots, periodic = TRUE)
+}
+    
diff --git a/R/data.R b/R/data.R
index 06409f452787c8daae6a05ddc5567d12f4ff6cd9..ec202c0a3c7b10c3788bfd4b0580521775943986 100644
--- a/R/data.R
+++ b/R/data.R
@@ -11,8 +11,8 @@
 #'   \item{t}{Time in GMT as POSIXct}
 #'   \item{heatload}{The heatload of a single family building in W}
 #'   \item{heatloadtotal}{The average heatload of a 16 single family buildings in W}
-#'   \item{Ta.obs}{Observed ambient temperature at the weather station in Celcius}
-#'   \item{I.obs}{Observed global radiation at the weather station in W/m^2}
+#'   \item{Taobs}{Observed ambient temperature at the weather station in Celcius}
+#'   \item{Iobs}{Observed global radiation at the weather station in W/m^2}
 #'   \item{Ta}{Weather forecasts of ambient temperature up to 36 hours ahead from DMI in Celcius}
 #'   \item{Ta}{Weather forecasts of global radiation up to 36 hours ahead from DMI in W/m^2}
 #' }
diff --git a/R/data.list.R b/R/data.list.R
index d09d42e1e9797d2c8ceb67a98df0da353383d3e4..9987b1eb8efcd79fd8d087092f39a91b8f53df3e 100644
--- a/R/data.list.R
+++ b/R/data.list.R
@@ -25,12 +25,12 @@
 #' # The time vector
 #' time <- seq(ct("2019-01-01"),ct("2019-01-02"),by=3600)
 #' # Observations time series (as vector)
-#' x.obs <- rnorm(length(time))
+#' xobs <- rnorm(length(time))
 #' # Forecast input as data.frame
 #' X <- data.frame(matrix(rnorm(length(time)*3), ncol=3))
 #' names(X) <- pst("k",1:3)
 #' 
-#' D <- data.list(t=time, x.obs=x.obs, X=X)
+#' D <- data.list(t=time, xobs=xobs, X=X)
 #'
 #' # Check it
 #' check(D)
@@ -58,7 +58,7 @@ data.list <- function(...) {
 #' # Use the data.list with building heat load 
 #' D <- Dbuilding
 #' # Take a subset for the example
-#' D <- subset(D, 1:10, nms=c("t","Ta.obs","Ta","I.obs","I"), kseq=1:3)
+#' D <- subset(D, 1:10, nms=c("t","Taobs","Ta","Iobs","I"), kseq=1:3)
 #' 
 #' # Take subset index 2:4
 #' subset(D, 2:4)
@@ -74,7 +74,7 @@ data.list <- function(...) {
 #' subset(D, nms=c("I","Ta"), kseq = 1)
 #' 
 #' # Lag the forecasts such that they are aligned in time with observations
-#' subset(D, nms=c("Ta.obs","Ta"), kseq = 2:3, lagforecasts = TRUE)
+#' subset(D, nms=c("Taobs","Ta"), kseq = 2:3, lagforecasts = TRUE)
 #' 
 #' # The order follows the order in nms
 #' subset(D, nms=c("Ta","I"), kseq = 2)
@@ -87,10 +87,10 @@ data.list <- function(...) {
 #' 
 #' # A scatter plot between the forecast and the observations
 #' # (try lagforecasts = FALSE and see the difference)
-#' plot(X$Ta$k10, X$Ta.obs)
+#' plot(X$Ta$k10, X$Taobs)
 #'
 #' # Fit a model for the 10-step horizon
-#' abline(lm(Ta.obs ~ Ta.k10, X), col=2)
+#' abline(lm(Taobs ~ Ta.k10, X), col=2)
 #'
 #' @export
 subset.data.list <- function(x, subset = NA, nms = NA, kseq = NA, lagforecasts = FALSE, pattern = NA, ...) {
@@ -200,7 +200,7 @@ subset.data.list <- function(x, subset = NA, nms = NA, kseq = NA, lagforecasts =
 #' #' # Use the data.list with building heat load 
 #' D <- Dbuilding
 #' # Take a subset
-#' D <- subset(D, 1:5, nms=c("t","Ta.obs","Ta","I.obs","I"), kseq=1:3)
+#' D <- subset(D, 1:5, nms=c("t","Taobs","Ta","Iobs","I"), kseq=1:3)
 #'
 #' # Convert to a data.frame, note the names of the forecasts are appended .kxx (i.e. for Ta and I)
 #' as.data.frame(D)
diff --git a/R/plot_ts.R b/R/plot_ts.R
index c2cb95392d298883882ac1f6440ae35477bfabbb..00103dacef60acaaae782b13696c63add7c70fdd 100644
--- a/R/plot_ts.R
+++ b/R/plot_ts.R
@@ -44,7 +44,7 @@
 #' # Make two plots (and set the space for the legend)
 #' plot_ts(D, c("heatload","Ta"), kseq=c(1,24), legendspace=11)
 #' # Only the Ta observations 
-#' plot_ts(D, c("heatload","Ta.obs$"), kseq=c(1,24), legendspace=11)
+#' plot_ts(D, c("heatload","Taobs$"), kseq=c(1,24), legendspace=11)
 #'
 #' # Give labels
 #' plot_ts(D, c("heatload","Ta"), kseq=c(1,24), xlab="Time", ylabs=c("Heat (kW)","Temperature (C)"))
diff --git a/R/plotly_ts.R b/R/plotly_ts.R
index 97c04458cca0c4ac5eee70b4ecafb813b12f4558..5f34a218b61e7b25be3eca6bf82744cb22777f9f 100644
--- a/R/plotly_ts.R
+++ b/R/plotly_ts.R
@@ -20,7 +20,7 @@
 #' D <- Dbuilding
 #' plotly_ts(D, c("heatload","Ta"), kseq=c(1,24))
 #' plotly_ts(D, c("heatload","Ta"), kseq=c(1,24))
-#' plotly_ts(D, c("heatload","Ta$|Ta.obs$"), kseq=c(1,24))
+#' plotly_ts(D, c("heatload","Ta$|Taobs$"), kseq=c(1,24))
 #' }
 #'
 #' @export
diff --git a/R/rls_fit.R b/R/rls_fit.R
index f872fe2a44eba7405866e4cc819830b73582ab62..72db60e4000120921bb32b46769177861274c7d9 100644
--- a/R/rls_fit.R
+++ b/R/rls_fit.R
@@ -202,7 +202,7 @@ rls_fit <- function(prm=NA, model, data, scorefun = NA, returnanalysis = TRUE,
         # The estimated coefficients
         Lfitval <- getse(Lresult, "Theta", fun=as.data.frame)
         # Return the model validation data
-        invisible(structure(list(Yhat = Yhat, model = model$clone_deep(), data = data, Lfitval = Lfitval, scoreval = scoreval), class = c("forecastmodel_fit","rls_fit")))
+        invisible(structure(list(Yhat = Yhat, model = model$clone_deep(), data = data, datatr = datatr, Lfitval = Lfitval, scoreval = scoreval), class = c("forecastmodel_fit","rls_fit")))
     }else{
         # Only the summed score returned
         val <- sum(scoreval, na.rm = TRUE)
diff --git a/R/rls_update.R b/R/rls_update.R
index 6b1124fc706dd2466cb2aa56ed4fbfc5e4811241..72efe74289247f8106741dbdf388788b77a61ffe 100644
--- a/R/rls_update.R
+++ b/R/rls_update.R
@@ -133,7 +133,7 @@ rls_update <- function(model, datatr = NA, y = NA, runcpp=TRUE) {
             }
             #
             # Give names to the matrices
-            colnames(P) <- names(datatr)
+            colnames(R) <- names(datatr)
             colnames(Theta) <- names(datatr)
             # Return the fit and result
             return(list(fit=list(k=k, theta=theta, R=R, yhat=yhat[(n-kmax+1):n]),
diff --git a/R/rmse.R b/R/rmse.R
index afe1dd4123e32ac9a6871365698d4b70cd85896e..12de721a1433446a57d6cf8be810c1a57c0b98e9 100644
--- a/R/rmse.R
+++ b/R/rmse.R
@@ -13,7 +13,7 @@
 #' @title Computes the RMSE score.
 #' @param x a numerical vector of residuals.
 #' @return The RMSE score.
-#' @seealso \code{\link{score_for_k}()} for calculation of a score for the k'th horizon, and \code{\link{score_fit}()} which takes a forecastmodel fit and returns score taking scoreperiod etc. into account.
+#' @seealso \code{\link{score}()} for calculation of a score for the k'th horizon, and \code{\link{score_fit}()} which takes a forecastmodel fit and returns score taking scoreperiod etc. into account.
 #' @name rmse
 #' @examples
 #'
diff --git a/R/score_fit.R b/R/score_fit.R
index 5752624cc37fde712fdabcc77a522019760b0d90..5ddb7a0f41ad39d84ba63ce0ec02ad5dce8dddac 100644
--- a/R/score_fit.R
+++ b/R/score_fit.R
@@ -8,7 +8,7 @@
 #' @param scoreperiod as an index (logical or integer) defining which points to inlude in the score calculation. If NA, the \code{scoreperiod} from the \code{fit$model} object is used.
 #' @param usecomplete Only use points where 
 #' @param scorefun The score function applied, per default \code{\link{rmse}}.
-#' @seealso \code{\link{rmse}} and \code{\link{score_for_k}} which are used in this function.
+#' @seealso \code{\link{rmse}} and \code{\link{score}} which are used in this function.
 #' @return A list with:
 #'   - \code{scoreval} is the score value
 #'   - \code{scoreperiod} is the score period used (can be different that the one in arguments \code{fit} or \code{scoreperiod})
@@ -41,7 +41,7 @@ score_fit <- function(fit, scoreperiod = NA, usecomplete = TRUE, scorefun = rmse
     }
 
     # Calculate the score
-    tmp <- score_for_k(Residuals, scoreperiod, usecomplete)
+    tmp <- score(Residuals, scoreperiod, usecomplete)
     scoreval <- tmp$scoreval
     scoreperiod <- tmp$scoreperiod
     
diff --git a/R/score_for_k.R b/R/score_for_k.R
index 45bffeb874bd4588128d6a2617d558f369b0c0d6..6b9cea980ce759697539a464534f982781699727 100644
--- a/R/score_for_k.R
+++ b/R/score_for_k.R
@@ -2,7 +2,7 @@
 #library(devtools)
 #document()
 #load_all(as.package("../../onlineforecast"))
-#?score_for_k
+#?score
 
 #' Calculates the score for each horizon for a matrix with residuals for each horizon.
 #'
@@ -24,15 +24,15 @@
 #' Resid <- residuals(Yhat, y)
 #'
 #' # Calculate the score for the k1 horizon
-#' score_for_k(Resid)$scoreval
+#' score(Resid)$val
 #'
 #' # The first values were excluded, since there are NAs
 #' head(Resid)
-#' score_for_k(Resid)$scoreperiod
+#' score(Resid)$scoreperiod
 #'
 #' @importFrom stats complete.cases
 #' @export
-score_for_k <- function(Residuals, scoreperiod = NA, usecomplete = TRUE, scorefun = rmse){
+score <- function(Residuals, scoreperiod = NA, usecomplete = TRUE, scorefun = rmse){
     # If no scoreperiod is given, then use all
     if(is.na(scoreperiod[1])){
         scoreperiod <- rep(TRUE,nrow(Residuals))
@@ -50,10 +50,10 @@ score_for_k <- function(Residuals, scoreperiod = NA, usecomplete = TRUE, scorefu
         scoreperiod <- scoreperiod & complete.cases(Residuals)
     }
     # Calculate the objective function for each horizon
-    scoreval <- sapply(1:ncol(Residuals), function(i){
+    val <- sapply(1:ncol(Residuals), function(i){
         scorefun(Residuals[scoreperiod,i])
     })
-    nams(scoreval) <- gsub("h","k",nams(Residuals))
+    nams(val) <- gsub("h","k",nams(Residuals))
     # 
-    return(list(scoreval=scoreval,scoreperiod=scoreperiod))
+    return(list(val=val,scoreperiod=scoreperiod))
 }
diff --git a/data/Dbuilding.rda b/data/Dbuilding.rda
index 236ca3811d9dc73d7fd42dddd107493fd48a17a7..6c952cd91c28e38a7504d126f5d2332c69c3962a 100644
Binary files a/data/Dbuilding.rda and b/data/Dbuilding.rda differ
diff --git a/data/all/make.R b/data/all/make.R
index 0fd52226b088e892e4478efce074066828c19e6d..279da6f6144ad537986c0751ca17d7c6751f871d 100644
--- a/data/all/make.R
+++ b/data/all/make.R
@@ -53,6 +53,8 @@ names(data[["Heatload"]]) <- pst("house", 1:16)
 data[["cosAoi"]] <- data_or[, "cosAoi.obs"]
 data[["sunElevation"]] <- data_or[, "sunElevation.obs"]
 
+# .obs to obs
+names(data) <- gsub("\\.","",names(data))
 
 # # The time of day
 # ncol <- ncol(data$Ta)
@@ -78,9 +80,9 @@ data$totaln <- sapply(1:nrow(data$Heatload), function(i){
 plot(data$t, data$totaln)
 
 # Write for building heat load forecasting
-#Dbuilding <- subset(data, c("2010-12-15","2011-03-01"), nms=c("t","Heatload","Ta","I","Ws","Wd","Ta.obs","I.obs","Wd.obs","Ws.obs","cosAoi","sunElevation","tday"))
+#Dbuilding <- subset(data, c("2010-12-15","2011-03-01"), nms=c("t","Heatload","Ta","I","Ws","Wd","Taobs","Iobs","Wdobs","Wsobs","cosAoi","sunElevation","tday"))
 data$heatload <- data$Heatload$house9
-Dbuilding <- subset(data, c("2010-12-15","2011-03-01"), nms=c("t","heatload","heatloadtotal","Ta.obs","I.obs","Ta","I"))
+Dbuilding <- subset(data, c("2010-12-15","2011-03-01"), nms=c("t","heatload","heatloadtotal","Taobs","Iobs","Ta","I"))
 rownames(Dbuilding$Ta) <- NULL
 Dbuilding$Ta <- Dbuilding$Ta[ ,1:36]
 rownames(Dbuilding$I) <- NULL
diff --git a/make.R b/make.R
index c53cacfda18b432306ae2f6490d55ca91a9c4eb1..0e706057e8d8f5248d40b78be2801d9b1bb32411 100644
--- a/make.R
+++ b/make.R
@@ -16,7 +16,6 @@
 ## install.packages("plotly")
 
 
-
 #----------------------------------------------------------------
 # Go
 library(devtools)
@@ -67,6 +66,14 @@ build(".", vignettes=TRUE)
 install.packages("../onlineforecast_0.9.1.tar.gz")
 
 library(onlineforecast)
+# ----------------------------------------------------------------
+
+
+# ----------------------------------------------------------------
+# Build binary package
+#system("R CMD INSTALL --build ../onlineforecast_0.9.1.tar.gz")
+# ----------------------------------------------------------------
+
 
 
 # ----------------------------------------------------------------
diff --git a/vignettes/forecast-evaluation.Rmd b/vignettes/forecast-evaluation.Rmd
index 516a1a15a8ab1ddda2eb8d8aff4b6c9e66d1cb48..8da010ea71d6d8345fed6a5012e57d46b90271be 100644
--- a/vignettes/forecast-evaluation.Rmd
+++ b/vignettes/forecast-evaluation.Rmd
@@ -326,14 +326,14 @@ Now the residuals can be calculated and the score:
 # Use the residuals function
 R <- residuals(D$Yhat1, D$y)
 # And the score as a function of the horizon
-score_for_k(R, scoreperiod=ok)$scoreval
+score(R, scoreperiod=ok)$val
 ```
 
 
 Calculated the score (default is RMSE) for all models:
 ```{r}
 RMSE <- sapply(nms, function(nm){
-    score_for_k(residuals(D[[nm]],D$y), ok)$scoreval
+    score(residuals(D[[nm]],D$y), ok)$val
 })
 ```
     
@@ -386,7 +386,7 @@ fittmp <- rls_fit(model$prm, model, D)
 
 Finally, the score can be calculated on the period following the train period by:
 ```{r scorefit}
-score_fit(fittmp, !D$trainperiod)$scoreval
+score_fit(fittmp, !D$trainperiod)$val
 ```
 
 In this way it's rather easy to set up different schemes, like optimizing the
diff --git a/vignettes/make.R b/vignettes/make.R
index 6180d210cc85bab537087e0f5c13d44883e0d7c2..89532033717af48b131aa7e02106dcc506c2ce51 100644
--- a/vignettes/make.R
+++ b/vignettes/make.R
@@ -19,7 +19,7 @@ makeit <- function(nam, openit=FALSE, clean=TRUE){
 
 #
 unlink(paste0(dirnam,"tmp-setup-data/"), recursive=TRUE)
-makeit("setup-data", openit=FALSE)
+makeit("setup-data", openit=TRUE)
 
 #
 unlink(paste0(dirnam,"tmp-setup-and-use-model/"), recursive=TRUE)
diff --git a/vignettes/setup-and-use-model.Rmd b/vignettes/setup-and-use-model.Rmd
index bae49a76878adc620ded08c14954274380e31f69..45725b198a849c479b1423678f4fe2797250d085 100644
--- a/vignettes/setup-and-use-model.Rmd
+++ b/vignettes/setup-and-use-model.Rmd
@@ -441,7 +441,7 @@ See the help `?make_tday` for more details.
 
 If we want to use observations in inputs to a model, we can use e.g.:
 ```{r}
-D$Tao <- make_input(D$Ta.obs, kseq=1:36)
+D$Tao <- make_input(D$Taobs, kseq=1:36)
 model$add_inputs(Tao = "lp(Tao, a1=0.99)")
 ```
 
diff --git a/vignettes/setup-data.Rmd b/vignettes/setup-data.Rmd
index c2371e949183d5f2a5a02e5b5aa8e629a3b4e68a..7c646169a942aad3b53353c2f7e6a86652c74793 100644
--- a/vignettes/setup-data.Rmd
+++ b/vignettes/setup-data.Rmd
@@ -284,7 +284,7 @@ observations) can be generated by:
 i <- 200:296
 plot(D$t[i], D$I$k8[i], type="l", col=2, xlab="Time", ylab="Global radiation (W/m²)")
 ## Add the observations
-lines(D$t[i], D$I.obs[i])
+lines(D$t[i], D$Iobs[i])
 legend("topright", c("8-step forecasts","Observations"), bg="white", lty=1, col=2:1)
 ```
 
@@ -292,7 +292,7 @@ Notice how the are not aligned, since the forecasts are 8 hours ahead. To align
 them the forecasts must be lagged 8 steps by:
 ```{r}
 plot(D$t[i], lagvec(D$I$k8[i], 8), type="l", col=2, xlab="Time", ylab="Global radiation (W/m²)")
-lines(D$t[i], D$I.obs[i])
+lines(D$t[i], D$Iobs[i])
 legend("topright", c("8-step forecasts lagged","Observations"), bg="white", lty=1, col=2:1)
 ```
 
@@ -341,7 +341,7 @@ So lagging (thus aligning in time) makes less slightly less scatter.
 A wrapper for the `pairs` function is provided for a `data.list`, which can
 generate very useful explorative plots:
 ```{r, fig.height=figwidth}
-pairs(D, nms=c("heatload","Ta.obs","Ta","t"), kseq=c(1,8,24))
+pairs(D, nms=c("heatload","Taobs","Ta","t"), kseq=c(1,8,24))
 ```
 Note how the sequence of included horizons are specified in the `kseq` argument,
 and note that the forecasts are lagged to be aligned in time. See `?pairs.data.list` for more details.
@@ -354,7 +354,7 @@ x <- lagvec(D$I$k8, 8)
 ## Take a smaller range
 x <- x[i]
 ## Take the observations
-y <- D$I.obs[i]
+y <- D$Iobs[i]
 ## Fit a linear regression model
 fit <- lm(y ~ x)
 ## Plot the result
@@ -429,6 +429,3 @@ Dsub2 <- as.data.list(Df)
 summary(Dsub2)
 summary(Dsub)
 ```
-
-
-## Literature