diff --git a/R/GloSpheGeoReg.R b/R/GloSpheGeoReg.R
index 04025ec..bcde908 100644
--- a/R/GloSpheGeoReg.R
+++ b/R/GloSpheGeoReg.R
@@ -3,8 +3,8 @@
 #' @import trust trust
 
 GloSpheGeoReg <- function(xin, yin, xout) {
-  k = length(xout)
-  n = length(xin)
+  k = nrow(xout)
+  n = nrow(xin)
   m = ncol(yin)
   
   xbar <- colMeans(xin)
diff --git a/R/GloSpheReg.R b/R/GloSpheReg.R
index b4b78ab..538d757 100644
--- a/R/GloSpheReg.R
+++ b/R/GloSpheReg.R
@@ -46,7 +46,7 @@ GloSpheReg <- function(xin=NULL, yin=NULL, xout=NULL){
   if(is.vector(xout)){
     xout <- as.matrix(xout)
   }
-  if (length(xin)!=nrow(yin))
+  if (nrow(xin)!=nrow(yin))
     stop("The length of xin should be the same as the number of rows in yin.")
   if (sum(abs(rowSums(yin^2) - rep(1,nrow(yin))) > 1e-6)){
     yin = yin / sqrt(rowSums(yin^2))
diff --git a/tests/testthat/test_GloSpheReg.R b/tests/testthat/test_GloSpheReg.R
index 6ce9414..e3a9db1 100644
--- a/tests/testthat/test_GloSpheReg.R
+++ b/tests/testthat/test_GloSpheReg.R
@@ -22,6 +22,29 @@ test_that("truth: great circle (no noise)", {
   expect_lt( mean(sapply(seq_len(n), function(i) SpheGeoDist(res$yout[i,], ytrue[,i]))), 5e-3 )
 })
 
+test_that("truth: great circle (no noise) with multiple covariates", {
+  n <- 101
+  xin <- seq(-1,1,length.out = n)
+  theta_true <- rep(pi/2,n)
+  phi_true <- (xin + 1) * pi / 4
+  err_sd <- 0
+  ytrue <- apply( cbind( 1, phi_true, theta_true ), 1, pol2car )
+  if (err_sd > 0) {
+    basis <- apply( ytrue, 2, frameSphere ) # [1:3,] - 1st basis vector; [4:6,] - 2nd basis vector
+    set.seed(1)
+    yin_tg <- basis[1:3,] * rnorm(n, mean = 0, sd = err_sd) + 
+      basis[4:6,] * rnorm(n, mean = 0, sd = err_sd)
+    yin <- t( sapply( seq_len(n), function(i) expSphere( base = ytrue[,i], tg = yin_tg[,i] ) ) )
+  } else {
+    yin <- t( ytrue )
+  }
+  
+  xin <- cbind(xin, c(0, rep(c(-0.95, 0.95), 50)))
+  xout <- xin
+  res <- GloSpheReg(xin=xin, yin=yin, xout=xout)
+  expect_lt( mean(sapply(seq_len(n), function(i) SpheGeoDist(res$yout[i,], ytrue[,i]))), 5e-3 )
+})
+
 # How to generate data according to a global model?
 test_that("truth: great circle + noise", {
   n <- 101