diff --git a/DESCRIPTION b/DESCRIPTION
index 4966bc0..31eaf9c 100755
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -4,9 +4,13 @@ Version: 0.8.38
 Author: Thomas Wutzler
 Maintainer: Thomas Wutzler <twutz@bgc-jena.mpg.de>
 Description: Density, distribution, quantile and random generation function for the logitnormal distribution. Estimation of the mode and the first two moments. Estimation of distribution parameters.
-Depends: 
-Suggests: RUnit, knitr, markdown, ggplot2, reshape2
+Suggests: 
+    knitr,
+    markdown,
+    ggplot2,
+    testthat (>= 3.0.0)
 VignetteBuilder: knitr
 License: GPL-2
 LazyData: true
 RoxygenNote: 7.1.1
+Config/testthat/edition: 3
diff --git a/NEWS.md b/NEWS.md
index 6c3b038..5c4c64d 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,3 +1,7 @@
+# logitnorm 0.8.38
+- Reworked twCoefLogitnormMLEFlat
+- Migrated from RUnit to testthat
+
 # logitnorm 0.8.37
 
 - Density returns 0 instead of NAN outside (0,1)
diff --git a/R/logitnorm.R b/R/logitnorm.R
index c1bc52d..5ccd96a 100755
--- a/R/logitnorm.R
+++ b/R/logitnorm.R
@@ -97,6 +97,8 @@ qlogitnorm <- function(
   plogis(qn)
 }
 
+# derivative of logit(x)
+.dlogit <- function(x) 1/x + 1/(1-x)
 
 #------------------ estimating parameters from fitting to distribution 
 .ofLogitnorm <- function(
@@ -307,8 +309,6 @@ twCoefLogitnormMLE <- function(
 attr(twCoefLogitnormMLE,"ex") <- function(){
   # estimate the parameters, with mode 0.7 and upper quantile 0.9
   mode = 0.7; upper = 0.9
-  mode = 0.2; upper = 0.7
-  #mode = 0.5; upper = 0.9
   (theta <- twCoefLogitnormMLE(mode,upper))
   x <- seq(0,1,length.out = 41)[-c(1,41)]	# plotting grid
   px <- plogitnorm(x,mu = theta[1],sigma = theta[2])	#percentiles function
@@ -321,18 +321,56 @@ attr(twCoefLogitnormMLE,"ex") <- function(){
   (theta <- twCoefLogitnormMLEFlat(mode))
 }
 
+.tmp.f <- function(){
+  (theta = twCoefLogitnormMLEFlat(0.9))
+  (theta = twCoefLogitnormMLEFlat(0.1))
+  (theta = twCoefLogitnormMLEFlat(0.5))
+  mle = c(0.7,0.2)
+  mle = seq(0.55,0.95,by=0.05)
+  mle = seq(0.95,0.99,by=0.01)
+  theta = twCoefLogitnormMLEFlat(mle)
+  sigma = theta[,2]
+  plot(sigma ~ mle)
+  (sigmac = 1/(2*mle*(1-mle)))
+  sigmac/sigma
+  plot(sigmac ~ mle)
+}
+
 twCoefLogitnormMLEFlat <- function(
-  ### Estimating coefficients of a maximally flat unimodal logitnormal distribution from mode 	
-  mle						##<< numeric vector: the mode of the density function
+  ### Estimating coefficients of a maximally flat unimodal logitnormal distribution given the mode 	
+  mle		##<< numeric vector: the mode of the density function
 ){
-  ##details<<
-  ## When increasing the sigma parameter, the distribution becomes
-  ## eventually becomes bi-model, i.e. has two maxima.	
-  ## This function estimates parameters for given mode, so that the distribution assigns high  
-  ## density to a maximum range, i.e. is maximally flat, but still is unimodal.
-  twCoefLogitnormMLE(mle,0)
+  vectorList <- lapply(mle, .twCoefLogitnormMLEFlat_scalar)
+  do.call(rbind, vectorList)
+}
+  
+  
+.twCoefLogitnormMLEFlat_scalar <- function(mle){
+  if (mle == 0.5) {
+    #sigma =  sqrt(dlogit(mle)/2) # 1.414214
+    return(c(mu = 0.0, sigma = sqrt(2)))
+  } 
+  is_right <- mle > 0.5
+  mler = ifelse(is_right, mle, 1-mle)
+  xt <- optimize(.ofModeFlat, interval = c(0,0.5), m=mler, logitm=logit(mler))$minimum
+  sigma2 <- (1/xt + 1/(1-xt))/2
+  mur <- logit(mler) - sigma2*(2*mler - 1)
+  mu <- ifelse(is_right, mur, -mur)
+  c(mu = mu, sigma = sqrt(sigma2))
 }
 
+.ofModeFlat <- function(x, m, logitm = logit(m)){
+  # lhs = 1/x + 1/(1-x)
+  # rhs = (logitm - logit(x))/(m-x)
+  # lhs = (m-x)/x + (m-x)/(1-x)
+  # rhs = (logitm - logit(x))
+  lhs = m/x + (m-1)/(1-x)
+  rhs = logitm - logit(x)
+  d = lhs - rhs
+  d*d
+}
+
+
 .ofLogitnormE <- function(
   ### Objective function used by \code{\link{coefLogitnormE}}. 
   theta				##<< theta[1] is the mu, theta[2] is the sigma
diff --git a/man/twCoefLogitnormMLE.Rd b/man/twCoefLogitnormMLE.Rd
index 0147aba..5e5b3b1 100755
--- a/man/twCoefLogitnormMLE.Rd
+++ b/man/twCoefLogitnormMLE.Rd
@@ -21,8 +21,6 @@ rows correspond to rows in \code{mle}, \code{quant}, and \code{perc}}
 \examples{
 # estimate the parameters, with mode 0.7 and upper quantile 0.9
 mode = 0.7; upper = 0.9
-mode = 0.2; upper = 0.7
-#mode = 0.5; upper = 0.9
 (theta <- twCoefLogitnormMLE(mode,upper))
 x <- seq(0,1,length.out = 41)[-c(1,41)]	# plotting grid
 px <- plogitnorm(x,mu = theta[1],sigma = theta[2])	#percentiles function
diff --git a/man/twCoefLogitnormMLEFlat.Rd b/man/twCoefLogitnormMLEFlat.Rd
index 6f246bb..5362514 100644
--- a/man/twCoefLogitnormMLEFlat.Rd
+++ b/man/twCoefLogitnormMLEFlat.Rd
@@ -1,15 +1,12 @@
 \name{twCoefLogitnormMLEFlat}
 \alias{twCoefLogitnormMLEFlat}
 \title{twCoefLogitnormMLEFlat}
-\description{Estimating coefficients of a maximally flat unimodal logitnormal distribution from mode 	}
+\description{Estimating coefficients of a maximally flat unimodal logitnormal distribution given the mode 	}
 \usage{twCoefLogitnormMLEFlat(mle)}
 \arguments{
   \item{mle}{numeric vector: the mode of the density function}
 }
-\details{When increasing the sigma parameter, the distribution becomes
-eventually becomes bi-model, i.e. has two maxima.	
-This function estimates parameters for given mode, so that the distribution assigns high  
-density to a maximum range, i.e. is maximally flat, but still is unimodal.}
+
 
 
 \author{Thomas Wutzler}
diff --git a/tests/doRUnit.R b/tests/doRUnit.R
deleted file mode 100644
index 38d478c..0000000
--- a/tests/doRUnit.R
+++ /dev/null
@@ -1,63 +0,0 @@
-## unit tests will not be done if RUnit is not available
-if (require("RUnit", quietly = TRUE)) {
-	
-	## --- Setup ---
-	
-	#pkg <- "PKG" # <-- Change to package name!
-	pkg <- "logitnorm" # <-- Change to package name!
-	if (Sys.getenv("RCMDCHECK") == "FALSE") {
-		## Path to unit tests for standalone running under Makefile (not R CMD check)
-		## PKG/tests/../inst/unitTests
-		path <- file.path(getwd(), "..", "inst", "unitTests")
-	} else {
-		## Path to unit tests for R CMD check
-		## PKG.Rcheck/tests/../PKG/unitTests
-		path <- system.file(package = pkg, "unitTests")
-	}
-	cat("\nRunning unit tests\n")
-	print(list(pkg = pkg, getwd = getwd(), pathToUnitTests = path))
-	
-	library(package = pkg, character.only = TRUE)
-	
-	## If desired, load the name space to allow testing of private functions
-	## if (is.element(pkg, loadedNamespaces()))
-	##     attach(loadNamespace(pkg), name = paste("namespace", pkg, sep = ":"), pos = 3)
-	##
-	## or simply call PKG:::myPrivateFunction() in tests
-	
-	## --- Testing ---
-	
-	## Define tests
-	testSuite <- defineTestSuite(name = paste(pkg, "unit testing"),
-			dirs = path)
-	## Run
-	tests <- runTestSuite(testSuite)
-	
-	## Return stop() to cause R CMD check stop in case of
-	##  - failures i.e. FALSE to unit tests or
-	##  - errors i.e. R errors
-	tmp <- getErrors(tests)
-	if (tmp$nFail > 0 | tmp$nErr > 0) {
-		stop(paste("\n\nunit testing failed (#test failures: ", tmp$nFail,
-						", #R errors: ",  tmp$nErr, ")\n\n", sep = ""))
-	}
-} else {
-	warning("cannot run unit tests -- package RUnit is not available")
-}
-
-.tmp.f <- function(){
-  #twutz 1712 removed writing test protocols, as its against CRAN policy
-  ## Default report name
-  pathReport <- file.path(path, "report")
-  
-  ## Report to stdout and text files
-  cat("------------------- UNIT TEST SUMMARY ---------------------\n\n")
-  printTextProtocol(tests, showDetails = FALSE)
-  printTextProtocol(tests, showDetails = FALSE,
-                    fileName = paste(pathReport, "Summary.txt", sep = ""))
-  printTextProtocol(tests, showDetails = TRUE,
-                    fileName = paste(pathReport, ".txt", sep = ""))
-  
-  ## Report to HTML file
-  printHTMLProtocol(tests, fileName = paste(pathReport, ".html", sep = ""))  
-}
diff --git a/tests/testthat.R b/tests/testthat.R
new file mode 100644
index 0000000..41c8c29
--- /dev/null
+++ b/tests/testthat.R
@@ -0,0 +1,4 @@
+library(testthat)
+library(logitnorm)
+
+test_check("logitnorm")
diff --git a/inst/unitTests/runitdlogitnormLog.R b/tests/testthat/test_dlogitnormLog.R
similarity index 52%
rename from inst/unitTests/runitdlogitnormLog.R
rename to tests/testthat/test_dlogitnormLog.R
index e5ec086..f61508f 100644
--- a/inst/unitTests/runitdlogitnormLog.R
+++ b/tests/testthat/test_dlogitnormLog.R
@@ -1,22 +1,22 @@
-#require(logitnorm)
-
-test.dlogitnorm <- function() {
-	# Check dlogitnorm with arbitrary parameters at arbitrary points.
-	m <- 1.4
-	s <- 0.8
-	x <- c(0.1, 0.8);
-	# taken from https://en.wikipedia.org/w/index.php?title=Logit-normal_distribution&oldid=708557844
-	x.density <- 1/s/sqrt(2*pi) * exp(-(log(x/(1-x)) - m)^2/2/s^2) / x / (1-x)
-	stopifnot(all.equal(x.density, dlogitnorm(x, m, s)))
-	stopifnot(all.equal(log(x.density), dlogitnorm(x, m, s, log=TRUE)))
-}
-
-test.plogitnorm <- function() {
-	# Check plogitnorm with arbitrary parameters at arbitrary points.
-	m <- 1.4
-	s <- 0.8
-	x <- c(0.1, 0.8);
-	stopifnot(all.equal(pnorm(log(x/(1-x)), m, s), plogitnorm(x, m, s)))
-	stopifnot(all.equal(pnorm(log(x/(1-x)), m, s, log=TRUE), log(plogitnorm(x, m, s))))
-}
-
+#require(logitnorm)
+
+test_that("dlogitnorm",  {
+	# Check dlogitnorm with arbitrary parameters at arbitrary points.
+	m <- 1.4
+	s <- 0.8
+	x <- c(0.1, 0.8);
+	# taken from https://en.wikipedia.org/w/index.php?title=Logit-normal_distribution&oldid=708557844
+	x.density <- 1/s/sqrt(2*pi) * exp(-(log(x/(1-x)) - m)^2/2/s^2) / x / (1-x)
+	expect_equal(x.density, dlogitnorm(x, m, s))
+	expect_equal(log(x.density), dlogitnorm(x, m, s, log=TRUE))
+})
+
+test_that("plogitnorm",  {
+	# Check plogitnorm with arbitrary parameters at arbitrary points.
+	m <- 1.4
+	s <- 0.8
+	x <- c(0.1, 0.8);
+	expect_equal(pnorm(log(x/(1-x)), m, s), plogitnorm(x, m, s))
+	expect_equal(pnorm(log(x/(1-x)), m, s, log=TRUE), log(plogitnorm(x, m, s)))
+})
+
diff --git a/inst/unitTests/runitLogitnorm.R b/tests/testthat/test_logitnorm.R
old mode 100755
new mode 100644
similarity index 67%
rename from inst/unitTests/runitLogitnorm.R
rename to tests/testthat/test_logitnorm.R
index 68c0a83..9a9177e
--- a/inst/unitTests/runitLogitnorm.R
+++ b/tests/testthat/test_logitnorm.R
@@ -1,86 +1,105 @@
 .tmp.f <- function(){
-  library(RUnit)
+  library(testthat)
 }
 
-.setUp <- function() {
-	#library(MASS)
-	.setUpDf <- within( list(),{
-			x <- seq(0,1,length.out = 41)[-c(1,41)]; 
-			#x[1] = x[1] + .Machine$double.eps; 
-			#x[length(x)] <- x[length(x)]- .Machine$double.eps
-			lx <- logit(x)			
-		})
-	attach(.setUpDf)
+# .setUp <- function() {
+# 	#library(MASS)
+# 	.setUpDf <- within( list(),{
+# 			x <- seq(0,1,length.out = 41)[-c(1,41)]; 
+# 			#x[1] = x[1] + .Machine$double.eps; 
+# 			#x[length(x)] <- x[length(x)]- .Machine$double.eps
+# 			lx <- logit(x)			
+# 		})
+# 	attach(.setUpDf)
+# }
+# 
+# .tearDown <- function() {
+# 	#detach(.setUpDf)
+# 	detach()
+# }
+
+local_x_l <- function(env = parent.frame()) {
+	env$x <- x <- seq(0,1,length.out = 41)[-c(1,41)];
+	env$lx <- lx <- logit(x)
 }
 
-.tearDown <- function() {
-	#detach(.setUpDf)
-	detach()
+.tmp.f <- function(){
+  local_x_l()
+  print(x) # replaced x
 }
+.tmp.f()
 
-
-
-test.inverseSame <- function(){
+test_that("inverseSame", {
+  local_x_l()
+  x
 	xnorm <- logit(x)
 	xinv <- invlogit(xnorm)
-	checkEquals(x, xinv)
-}
+	expect_equal(x, xinv)
+})
 
-test.plogitnorm <- function(){
-	px <- plogitnorm(x)	#percentiles
-	checkEquals( pnorm(lx), px)
+test_that("plogitnorm", {
+  local_x_l()
+  px <- plogitnorm(x)	#percentiles
+	expect_equal( pnorm(lx), px)
 	px2 <- plogitnorm(x,mu = 2,sigma = 1)
-	checkEquals( pnorm(lx,mean = 2,sd = 1), px2)
+	expect_equal( pnorm(lx,mean = 2,sd = 1), px2)
 	#plot( px ~ x)
 	#plot( px ~ logit(x))
 	#lines( pnorm( logit(x)) ~ logit(x) )
-}
+})
 
-test.dlogitnorm <- function(){
+test_that("dlogitnorm", {
+  local_x_l()
   q <- c(-1,0,0.5,1,2)
   set.seed(0815)
   ans <- suppressWarnings(dlogitnorm(q))
-  checkEquals(c(0,0,1.595769,0,0), ans, tolerance = 1e-7)
-}
+  expect_equal(c(0,0,1.595769,0,0), ans, tolerance = 1e-7)
+})
 
 
-test.twCoefLogitnorm <- function(){
-	theta <- twCoefLogitnorm(0.7,0.9,perc = 0.999)
+test_that("twCoefLogitnorm", {
+  local_x_l()
+  theta <- twCoefLogitnorm(0.7,0.9,perc = 0.999)
 	px <- plogitnorm(x,mu = theta[1],sigma = theta[2])	#percentiles function
 	dx <- dlogitnorm(x,mu = theta[1],sigma = theta[2])	#density function
 	#plot(px~x); abline(v = c(0.7,0.9)); abline(h = c(0.5,0.975))
 	#plot(dx~x); abline(v = c(0.7,0.9))
 	# upper percentile at 0.9
-	checkEquals(which.min(abs(px - 0.999)), which(x == 0.9) )
-	checkEquals(which.min(abs(px - 0.5)), which.min(abs(x - 0.7)) )
+	expect_equal(which.min(abs(px - 0.999)), which(x == 0.9) )
+	expect_equal(which.min(abs(px - 0.5)), which.min(abs(x - 0.7)) )
 	# mode at 0.7
-	#checkEquals(which(abs(x-0.7)<.Machine$double.eps), which.max(dx) )
-}
+	#expect_equal(which(abs(x-0.7)<.Machine$double.eps), which.max(dx) )
+})
 
-test.twCoefLogitnormN <- function(){
-	quant = c(0.7,0.8,0.9)
+test_that("twCoefLogitnormN", {
+  local_x_l()
+  quant = c(0.7,0.8,0.9)
 	perc = c(0.5,0.75,0.975)
 	(theta <- twCoefLogitnormN( quant = quant, perc = perc ))
+	# regression to previous results
+	expect_equal(theta, c(mu=0.86, sigma = 0.76), tolerance = 0.01)
 	#px <- plogitnorm(x,mu = theta[1],sigma = theta[2])	#percentiles function
 	#dx <- dlogitnorm(x,mu = theta[1],sigma = theta[2])	#density function
 	#plot(px~x); abline(v = quant,col = "gray"); abline(h = perc,col = "gray")
-}
+})
 
 
-test.twCoefLogitnormMLE <- function(){
-	theta <- twCoefLogitnormMLE(0.7,0.9,perc = 0.975)
+test_that("twCoefLogitnormMLE", {
+  local_x_l()
+  theta <- twCoefLogitnormMLE(0.7,0.9,perc = 0.975)
 	px <- plogitnorm(x,mu = theta[1],sigma = theta[2])	#percentiles function
 	dx <- dlogitnorm(x,mu = theta[1],sigma = theta[2])	#density function
 	#plot(px~x); abline(v = c(0.7,0.9)); abline(h = c(0.5,0.975))
 	#plot(dx~x); abline(v = c(0.7,0.9))
 	# upper percentile at 0.9
-	checkEquals(which.min(abs(px - 0.975)), which(x == 0.9) )
+	expect_equal(which.min(abs(px - 0.975)), which(x == 0.9) )
 	# mode at 0.7
-	checkEquals(which.min(abs(x - 0.7)), which.max(dx) )
-}
+	expect_equal(which.min(abs(x - 0.7)), which.max(dx) )
+})
 
 
-test.twCoefLogitnormE <- function(){
+test_that("twCoefLogitnormE", {
+  local_x_l()
   set.seed(0815)
   theta <- twCoefLogitnormE(0.7,0.9)
 	px <- plogitnorm(x,mu = theta[1],sigma = theta[2])	#percentiles function
@@ -88,16 +107,17 @@ test.twCoefLogitnormE <- function(){
 	#plot(px~x); abline(v = c(0.7,0.9)); abline(h = c(0.5,0.975))
 	#plot(dx~x); abline(v = c(0.7,0.9))
 	# upper percentile at 0.9
-	checkEquals(which.min(abs(px - 0.975)), which(x == 0.9) )
+	expect_equal(which.min(abs(px - 0.975)), which(x == 0.9) )
 	# mean at 0.7
-	checkEqualsNumeric( momentsLogitnorm(
-	  mu = theta[1],sigma = theta[2])["mean"], 0.7, tolerance = 1e-3)
+	expect_equal( unname(momentsLogitnorm(
+	  mu = theta[1],sigma = theta[2])["mean"]), 0.7, tolerance = 1e-3)
 	z <- rlogitnorm(1e5, mu = theta[1],sigma = theta[2])
-	checkEqualsNumeric(0.7, mean(z), tolerance = 5e-3 )
-}
+	expect_equal(unname(mean(z)), 0.7, tolerance = 5e-3 )
+})
 
 .tmp.f <- function(){
-	px <- plogitnorm(x)	#percentiles
+  local_x_l()
+  px <- plogitnorm(x)	#percentiles
 	plot( px ~ x )	
 	plot( qlogitnorm(px) ~ x ) 	#one to one line
 	plot( dlogitnorm(x,mu = 0.9) ~ x, type = "l" )
diff --git a/tests/testthat/test_logitnorm_flat.R b/tests/testthat/test_logitnorm_flat.R
new file mode 100644
index 0000000..b544ea9
--- /dev/null
+++ b/tests/testthat/test_logitnorm_flat.R
@@ -0,0 +1,45 @@
+expect_monotone_logit_slope <- function(theta, upper=0.5, lower=0.0, decreasing = FALSE) {
+  x <- seq(lower,upper,length.out = 41)[-c(1,41)]	# plotting grid
+  dx <- dlogitnorm(x,mu = theta[1],sigma = theta[2])	#density function
+  #plot(dx~x); abline(v = c(mode,upper),col = "gray")
+  if (decreasing)
+    expect_true(all(diff(dx) <= 0))
+  else
+    expect_true(all(diff(dx) >= 0))
+}
+
+test_that("twCoefLogitnormMLEFlat", {
+  # standard case
+  theta9 <- theta <- twCoefLogitnormMLEFlat(0.9)
+  xmode <- modeLogitnorm(theta[,"mu"],theta[,"sigma"])
+  expect_equal(xmode, 0.9, tolerance = 1e-3)
+  expect_monotone_logit_slope(theta, xmode)
+  expect_monotone_logit_slope(theta, 1, xmode, decreasing = TRUE)
+  #
+  # mle < 0.5
+  theta1 <- theta <- twCoefLogitnormMLEFlat(0.1)
+  expect_equal(theta1[,"sigma"], theta9[,"sigma"])
+  expect_equal(-theta1[,"mu"], theta9[,"mu"])
+  xmode <- modeLogitnorm(theta[,"mu"],theta[,"sigma"])
+  expect_equal(xmode, 0.1, tolerance = 1e-3)
+  expect_monotone_logit_slope(theta, xmode)
+  expect_monotone_logit_slope(theta, 1, xmode, decreasing = TRUE)
+  #
+  # mle == 0.5
+  theta5 <- theta <- twCoefLogitnormMLEFlat(0.5)
+  expect_equal(unname(theta5[,"mu"]), 0.0)
+  expect_equal(unname(theta5[,"sigma"]), sqrt(2)) # 1.414214)
+  xmode <- modeLogitnorm(theta[,"mu"],theta[,"sigma"])
+  expect_equal(xmode, 0.5, tolerance = 1e-3)
+  expect_monotone_logit_slope(theta, xmode)
+  expect_monotone_logit_slope(theta, 1, xmode, decreasing = TRUE)
+  # really maximum sigma?
+  theta[,"sigma"] <- theta[,"sigma"] +0.01
+  expect_error(expect_monotone_logit_slope(theta, xmode))
+  #
+  # vectorized
+  theta <- twCoefLogitnormMLEFlat((c(0.9,0.1,0.5)))
+  expect_equal(theta[1,], theta9[1,])
+  expect_equal(theta[2,], theta1[1,])
+  expect_equal(theta[3,], theta5[1,])
+})
diff --git a/inst/unitTests/runitmodeLogitnorm.R b/tests/testthat/test_modeLogitnorm.R
old mode 100755
new mode 100644
similarity index 69%
rename from inst/unitTests/runitmodeLogitnorm.R
rename to tests/testthat/test_modeLogitnorm.R
index f3b743e..4e8c3dc
--- a/inst/unitTests/runitmodeLogitnorm.R
+++ b/tests/testthat/test_modeLogitnorm.R
@@ -1,12 +1,4 @@
-#TODO
-
-.setUp <-function () {
-}
-
-.tearDown <- function () {
-}
-
-test.rightMode <- function(){
+test_that("rightMode", {
 	theta0 <- c(mu=1.5, sigma=0.8)
 	#plot the true and the rediscovered distributions
 	xGrid = seq(0,1, length.out=81)[-c(1,81)]
@@ -20,11 +12,11 @@ test.rightMode <- function(){
 	#check by monte carlo integration
 	#z <- rlogitnorm(1e6, mu=theta0[1], sigma=theta0[2]);	var(z)
 	#dz <- density(z)
-	#checkEqualsNumeric( dz$x[which.max(dz$y)], mle, tolerance=5e-2)
-	checkEqualsNumeric( 0.88, mle, tolerance=1e-2)
-}
+	#expect_equal( dz$x[which.max(dz$y)], mle, tolerance=5e-2)
+	expect_equal( 0.88, mle, tolerance=1e-2)
+})
 
-test.leftMode <- function(){
+test_that("leftMode", {
 	theta0 <- c(mu=-1.5, sigma=0.8)
 	#plot the true and the rediscovered distributions
 	xGrid = seq(0,1, length.out=81)[-c(1,81)]
@@ -39,7 +31,7 @@ test.leftMode <- function(){
 	# deprecated: did not run on Windows
 	#z <- rlogitnorm(1e6, mu=theta0[1], sigma=theta0[2]);	var(z)
 	#dz <- density(z)
-	#checkEqualsNumeric( dz$x[which.max(dz$y)], mle, tolerance=5e-2)
-	checkEqualsNumeric( 0.12, mle, tolerance=1e-2)		# regression 0.12 calculated previously
-}
+	#expect_equal( dz$x[which.max(dz$y)], mle, tolerance=5e-2)
+	expect_equal( 0.12, mle, tolerance=1e-2)		# regression 0.12 calculated previously
+})
 
diff --git a/inst/unitTests/runitmomentsLogitnorm.R b/tests/testthat/test_momentsLogitnorm.R
old mode 100755
new mode 100644
similarity index 52%
rename from inst/unitTests/runitmomentsLogitnorm.R
rename to tests/testthat/test_momentsLogitnorm.R
index 98dbcc3..3956b6e
--- a/inst/unitTests/runitmomentsLogitnorm.R
+++ b/tests/testthat/test_momentsLogitnorm.R
@@ -1,12 +1,4 @@
-#TODO
-
-.setUp <-function () {
-}
-
-.tearDown <- function () {
-}
-
-test.1 <- function(){
+test_that("1", {
   set.seed(0815)
   theta0 <- c(mu=1.5, sigma=0.8)
 	#plot the true and the rediscovered distributions
@@ -17,17 +9,17 @@ test.1 <- function(){
 	moments <- momentsLogitnorm(mu=theta0[1], sigma=theta0[2] )
 	#check by monte carlo integration
 	z <- rlogitnorm(1e6, mu=theta0[1], sigma=theta0[2]);	var(z)
-	checkEqualsNumeric( mean(z), moments["mean"], tolerance=1e-3)
-	checkEqualsNumeric( var(z), moments["var"], tolerance=1e-2)
-}
+	expect_equal(moments["mean"], c(mean=mean(z)), tolerance=1e-3)
+	expect_equal(moments["var"], c(var=var(z)), tolerance=1e-2)
+})
 
-test.momentsLogitnorm41 <- function(){
+test_that("momentsLogitnorm41", {
 	(res <- momentsLogitnorm(4,1))
-	checkEqualsNumeric( c(0.97189602, 0.00101663), res)
-}
+	expect_equal( unname(res), c(0.97189602, 0.00101663))
+})
 
-test.momentsLogitnorm501 <- function(){
+test_that("momentsLogitnorm501", {
   set.seed(0815)
   (res <- momentsLogitnorm(5,0.1))
-	checkEqualsNumeric( c(9.932743e-01, 4.484069e-07), res, tolerance=1e-7)
-}
+	expect_equal( unname(res), c(9.932743e-01, 4.484069e-07), tolerance=1e-7)
+})
diff --git a/vignettes/logitnorm.Rmd b/vignettes/logitnorm.Rmd
index 389e67f..1fcce6b 100644
--- a/vignettes/logitnorm.Rmd
+++ b/vignettes/logitnorm.Rmd
@@ -25,7 +25,7 @@ if( !require(ggplot2) ){
 	print("To generate this vignette, ggplo2 is required.")
 	exit(0)
 }
-library(reshape2)
+#library(reshape2)
 # genVigs("logitnorm")
 ```
 
@@ -49,21 +49,24 @@ sigma=10^seq(-0.5,0.5,length.out=5))
 n <- nrow(theta0)
 
 .calcDensityGrid <- function(
-		### Calculate lognormal desnity for given combinations
+		### Calculate logitnormal density for given combinations
 		theta0	##<< matrix with columns mu and sigma
 		,xGrid = seq(0,1, length.out=81)[-c(1,81)]
 ){ 
 	dx <- apply( theta0, 1, function(theta0i){
 				dx <- dlogitnorm(xGrid, mu=theta0i[1], sigma=theta0i[2])
 			})
-	dimnames(dx) <- list(iX=NULL,iTheta=NULL)
-	ds <- melt(dx)	# two variables over time
-	ds[1:10,]
-	ds$mu <- rep(as.factor(round(theta0[,1],2)), each=length(xGrid))
-	ds$sigma <- rep(as.factor(round(theta0[,2],2)), each=length(xGrid))
-	ds$x <- rep(xGrid, nrow(theta0))
-	ds
-	### data frame with columns value,x,mu and sigma 
+	# dimnames(dx) <- list(iX=NULL,iTheta=NULL)
+	# ds <- melt(dx)	# two variables over time
+	# ds <- gather(as.data.frame(dx))
+	# ds[1:10,]
+	ds <- data.frame(
+	  mu = rep(as.factor(round(theta0[,1],2)), each=length(xGrid)),
+	  sigma = rep(as.factor(round(theta0[,2],2)), each=length(xGrid)),
+	  x = rep(as.vector(xGrid), nrow(theta0)),
+	  value = as.vector(dx)
+	)
+	### data frame with columns value,x,mu and sigma, value (density at x) 
 }
 
 ds <- .calcDensityGrid(theta0,xGrid=xGrid)