test_that("means are equal", {
  x <- rnorm(500, sd = 2) # angle in radians
  expect_equal(circ_mean_rad(x), circ_descr(x)[['mu']])
  expect_equal(circ_mean_rad(x), as.vector(circular::mean.circular(circular::circular(x))))
  expect_equal(circ_mean_rad(x), circ_mean_360(x/pi*180)/180*pi)
  expect_equal(circ_mean_rad(x), weighted_circ_mean(x, rep(1, length(x))))
  expect_equal(circ_mean_rad(x), weighted_circ_mean2(x, rep(1, length(x))))
})


test_that("weighted means are equal", {
  x <- rnorm(500, sd = 2) # angle in radians
  w <- runif(500) # weights

  expect_equal(weighted_circ_mean(x, w), weighted_circ_mean2(x, w))

})

test_that("SDs are equal", {
  x <- rnorm(500, sd = 2) # angle in radians

  expect_equal(circ_sd_rad(x), circ_descr(x)[['sigma']])
  expect_equal(circ_sd_rad(x), circ_sd_360(x/pi*180)/180*pi)
  expect_equal(circ_sd_rad(x), weighted_circ_sd(x, rep(1, length(x))))
  expect_equal(circ_sd_rad(x), as.vector(circular::sd.circular(circular::circular(x))))


})

test_that("Circular correlation works properly for uniform and conditional von Mises", {
  # see Jammalamadaka & SenGupta, pp. 181-182
  kappa <- runif(1, 5, 200) # NB: for lower kappa, the results are not always within tolerance limits.
  n <- 1000000
  x <- runif(n, -pi, pi)
  y <- x + as.vector(circular::rvonmises(n, mu = circular::circular(0), kappa = kappa))
  exp_r <- a_fun(kappa)
  obs_r <- circ_corr(x, y, ill_defined = T)
  tolerance <- 1e-3
  print(sprintf('Difference between expected and observed correlation is %.6f', exp_r-obs_r))
  expect_lt(abs(exp_r-obs_r), tolerance)
})

test_that("Circular correlation close to Pearson for narrow cases", {
  # see Jammalamadaka & SenGupta
  data <- mgcv::rmvn(10000, c(0,0), V = matrix(c(1,0.5,0.5,1), ncol = 2))
  data <- data/8
  pearson_r <- cor(data[,1], data[,2])
  obs_r <- circ_corr(data[,1], data[,2])
  tolerance <- 1e-3
  print(sprintf('Difference between expected and observed correlation is %.6f', pearson_r-obs_r))
  expect_lt(abs(pearson_r-obs_r), tolerance)
})

test_that("Circular correlation matches BAMBI::circ_cor", {
  data <- mgcv::rmvn(10000, c(0,0), V = matrix(c(1,0.5,0.5,1), ncol = 2))
  data <- data*2
  obs_r <- circ_corr(data[,1], data[,2])
  exp_r <- BAMBI::circ_cor(data)[[1]]
  tolerance <- 1e-3
  print(sprintf('Difference between expected and observed correlation is %.6f', exp_r-obs_r))
  expect_lt(abs(exp_r-obs_r), tolerance)
})

test_that("conversion from circular SD to kappa works both ways", {
  test_sd_deg <- runif(1,0,100) # SD in degrees
  test_sd_rad <- test_sd_deg/180*pi

  kappa_from_deg <- vm_circ_sd_deg_to_kappa(test_sd_deg)
  kappa_from_rad <- vm_circ_sd_to_kappa(test_sd_rad)
  expect_equal(kappa_from_deg, kappa_from_rad)

  tolerance <- 1e-3
  expect_lt(abs(test_sd_deg-vm_kappa_to_circ_sd_deg(kappa_from_deg)), tolerance)
  expect_lt(abs(test_sd_rad-vm_kappa_to_circ_sd(kappa_from_deg)), tolerance)
})