diff --git a/DESCRIPTION b/DESCRIPTION
index 15e660b..5ab94ba 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -5,7 +5,9 @@ Version: 0.6.0
 Date: 2023-01-10
 Authors@R: c(
     person("Simon", "Garnier", email = "garnier@njit.edu", role = c("aut", "cre"),
-      comment = c(ORCID = "0000-0002-3886-3974"))
+      comment = c(ORCID = "0000-0002-3886-3974")),
+    person("Marina", "Papadopoulou", email = "m.papadopoulou.rug@gmail.com", role = c("ctb"),
+      comment = c(ORCID = "0000-0002-6478-8365"))
   )
 Maintainer: Simon Garnier <garnier@njit.edu>
 Description: Function library for processing collective movement data (e.g. fish
diff --git a/NAMESPACE b/NAMESPACE
index ade0951..eacc8db 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -21,6 +21,7 @@ export(linear_acc)
 export(linear_dist)
 export(linear_speed)
 export(nn)
+export(nnba)
 export(nnd)
 export(nsd)
 export(pdist)
diff --git a/R/nn.R b/R/nn.R
index 2e4bdba..1964d2f 100644
--- a/R/nn.R
+++ b/R/nn.R
@@ -1,35 +1,35 @@
 #' @title Pairwise Distance Matrix
-#' 
-#' @description Given a set of locations, this function computes the distances 
+#'
+#' @description Given a set of locations, this function computes the distances
 #'  between each possible pair of locations.
-#'  
-#' @param x A vector of x (or longitude) coordinates. 
-#' 
+#'
+#' @param x A vector of x (or longitude) coordinates.
+#'
 #' @param y A vector of y (or latitude) coordinates.
-#' 
-#' @param geo A logical value indicating whether the locations are defined by 
-#'  geographic coordinates (pairs of longitude/latitude values). Default: FALSE. 
-#'  
-#' @return A square matrix representing pairwise distances between each possible 
+#'
+#' @param geo A logical value indicating whether the locations are defined by
+#'  geographic coordinates (pairs of longitude/latitude values). Default: FALSE.
+#'
+#' @return A square matrix representing pairwise distances between each possible
 #'  pair of locations.
-#' 
+#'
 #' @author Simon Garnier, \email{garnier@@njit.edu}
-#' 
+#'
 #' @seealso \code{\link{nn}}, \code{\link{nnd}}
-#' 
+#'
 #' @examples
 #' x <- rnorm(25)
 #' y <- rnorm(25, sd = 3)
 #' pdist(x, y)
-#' 
+#'
 #' @export
 pdist <- function(x, y, geo = FALSE) {
-  if (length(x) != length(y)) 
+  if (length(x) != length(y))
     stop("x and y should have the same length.")
-  
+
   if (!is.numeric(x) | !is.numeric(y))
     stop("x and y should be numeric.")
-  
+
   if (geo) {
     l <- length(x)
     idx <- expand.grid(row = 1:l, col = 1:l)
@@ -44,45 +44,45 @@ pdist <- function(x, y, geo = FALSE) {
 
 #' @title Nearest Neighbor
 #'
-#' @description Given the locations of different objects, this function 
-#'  determines the identity of the nearest neighboring object to each object. 
-#' 
-#' @param x A vector of x (or longitude) coordinates. 
-#' 
+#' @description Given the locations of different objects, this function
+#'  determines the identity of the nearest neighboring object to each object.
+#'
+#' @param x A vector of x (or longitude) coordinates.
+#'
 #' @param y A vector of y (or latitude) coordinates.
-#' 
+#'
 #' @param id A vector corresponding to the unique identities of each track.
-#' 
-#' @param geo A logical value indicating whether the locations are defined by 
-#'  geographic coordinates (pairs of longitude/latitude values). Default: FALSE. 
-#'  
-#' @return A vector of the same length as x and y representing the identity of 
-#'  the nearest neighboring object to each object. 
-#' 
+#'
+#' @param geo A logical value indicating whether the locations are defined by
+#'  geographic coordinates (pairs of longitude/latitude values). Default: FALSE.
+#'
+#' @return A vector of the same length as x and y representing the identity of
+#'  the nearest neighboring object to each object.
+#'
 #' @author Simon Garnier, \email{garnier@@njit.edu}
-#' 
+#'
 #' @seealso \code{\link{nnd}}
-#' 
+#'
 #' @examples
 #' x <- rnorm(25)
 #' y <- rnorm(25, sd = 3)
 #' id <- 1:25
 #' nn(x, y, id)
-#' 
+#'
 #' @export
 nn <- function(x, y, id, geo = FALSE) {
   if (!all(length(x) == c(length(y), length(id))))
     stop("x, y and id should have the same length.")
-  
+
   if (!is.numeric(x) | !is.numeric(y))
     stop("x and y should be numeric.")
-  
+
   d <- pdist(x, y, geo = geo)
   diag(d) <- NA
   d[is.na(x) | is.na(y), ] <- NA
   d[, is.na(x) | is.na(y)] <- NA
-  
-  idx <- apply(d, 2, 
+
+  idx <- apply(d, 2,
                function(x) {
                  if (sum(is.na(x)) != length(x)) {
                    which(x == min(x, na.rm = TRUE))[1]
@@ -96,43 +96,43 @@ nn <- function(x, y, id, geo = FALSE) {
 
 #' @title Nearest Neihgbor Distance
 #'
-#' @description Given the locations of different objects, this function 
-#'  determines the distance of the nearest neighboring object to each object. 
-#' 
-#' @param x A vector of x (or longitude) coordinates. 
-#' 
+#' @description Given the locations of different objects, this function
+#'  determines the distance of the nearest neighboring object to each object.
+#'
+#' @param x A vector of x (or longitude) coordinates.
+#'
 #' @param y A vector of y (or latitude) coordinates.
-#' 
-#' @param geo A logical value indicating whether the locations are defined by 
-#'  geographic coordinates (pairs of longitude/latitude values). Default: FALSE. 
-#'  
-#' @return A vector of the same length as x and y representing the distance to 
-#'  the nearest neighboring object for each object. 
-#' 
+#'
+#' @param geo A logical value indicating whether the locations are defined by
+#'  geographic coordinates (pairs of longitude/latitude values). Default: FALSE.
+#'
+#' @return A vector of the same length as x and y representing the distance to
+#'  the nearest neighboring object for each object.
+#'
 #' @author Simon Garnier, \email{garnier@@njit.edu}
-#' 
+#'
 #' @seealso \code{\link{nn}}
-#' 
+#'
 #' @examples
 #' x <- rnorm(25)
 #' y <- rnorm(25, sd = 3)
 #' id <- 1:25
 #' nnd(x, y)
-#' 
+#'
 #' @export
 nnd <- function(x, y, geo = FALSE) {
-  if (length(x) != length(y)) 
+  if (length(x) != length(y))
     stop("x and y should have the same length.")
-  
+
   if (!is.numeric(x) | !is.numeric(y))
     stop("x and y should be numeric.")
-  
+
   d <- pdist(x, y, geo = geo)
   diag(d) <- NA
   d[is.na(x) | is.na(y), ] <- NA
   d[, is.na(x) | is.na(y)] <- NA
-  
-  apply(d, 2, 
+
+  apply(d, 2,
         function(x) {
           if (sum(is.na(x)) != length(x)) {
             min(x, na.rm = TRUE)
@@ -141,3 +141,70 @@ nnd <- function(x, y, geo = FALSE) {
           }
         })
 }
+
+
+#' @title Nearest Neighbor Bearing Angle
+#'
+#' @description Given the locations and headings of different objects,
+#'  this function determines the angle between the heading of each object
+#'  and the position to the nearest neighboring object (bearing angle).
+#'
+#' @param x A vector of x (or longitude) coordinates.
+#'
+#' @param y A vector of y (or latitude) coordinates.
+#'
+#' @param hs A vector of headings (angle in rads).
+#'
+#' @param geo A logical value indicating whether the locations are defined by
+#'  geographic coordinates (pairs of longitude/latitude values). Default: FALSE.
+#'
+#' @return A vector of the same length as x, y and hs representing the bearing
+#' angle to the nearest neighboring object for each object.
+#'
+#' @author Simon Garnier, \email{garnier@@njit.edu},
+#'  Marina Papadopoulou, \email{m.papadopoulou.rug@@gmail.com}
+#'
+#' @seealso \code{\link{pdist}}
+#'
+#' @examples
+#' x <- rnorm(25)
+#' y <- rnorm(25, sd = 3)
+#' hs <- rnorm(25, sd = 1)
+#' nnba(x, y, hs)
+#'
+#' @export
+nnba <- function(x, y, hs, geo = FALSE) {
+  if (!all(length(x) == c(length(y), length(hs))))
+    stop("x, y and hs should have the same length.")
+
+  if (!is.numeric(x) || !is.numeric(y) || !is.numeric(hs))
+    stop("x, y and hs should be numeric.")
+
+  d <- swaRm::pdist(x, y, geo = geo)
+  diag(d) <- NA
+  d[is.na(x) | is.na(y), ] <- NA
+  d[, is.na(x) | is.na(y)] <- NA
+  idx <- apply(d, 2, function(x) {
+    if (sum(is.na(x)) != length(x)) {
+      which(x == min(x, na.rm = TRUE))[1]
+    } else {
+      as.numeric(NA)
+    }
+  })
+
+  if (geo) {
+    m1 <- cbind(x, y)
+    m2 <- cbind(x[idx], y[idx])
+    br <- geosphere::bearing(m1, m2) * pi / 180
+  } else {
+    dy <- y[idx] - y
+    dx <- x[idx] - x
+    br <- atan2(y = dy, x = dx)
+  }
+  db <- hs - br
+
+  db[db <= (-pi) & !is.na(db)] <- 2 * pi + db[db <= (-pi) & !is.na(db)]
+  db[db > pi & !is.na(db)] <- db[db > pi & !is.na(db)] - 2 * pi
+
+  db
+}
diff --git a/man/nn.Rd b/man/nn.Rd
index 8d95f98..b1f58a3 100644
--- a/man/nn.Rd
+++ b/man/nn.Rd
@@ -13,15 +13,15 @@ nn(x, y, id, geo = FALSE)
 
 \item{id}{A vector corresponding to the unique identities of each track.}
 
-\item{geo}{A logical value indicating whether the locations are defined by 
+\item{geo}{A logical value indicating whether the locations are defined by
 geographic coordinates (pairs of longitude/latitude values). Default: FALSE.}
 }
 \value{
-A vector of the same length as x and y representing the identity of 
+A vector of the same length as x and y representing the identity of
  the nearest neighboring object to each object.
 }
 \description{
-Given the locations of different objects, this function 
+Given the locations of different objects, this function
  determines the identity of the nearest neighboring object to each object.
 }
 \examples{
diff --git a/man/nnba.Rd b/man/nnba.Rd
new file mode 100644
index 0000000..72b0458
--- /dev/null
+++ b/man/nnba.Rd
@@ -0,0 +1,41 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/nn.R
+\name{nnba}
+\alias{nnba}
+\title{Nearest Neighbor Bearing Angle}
+\usage{
+nnba(x, y, hs, geo = FALSE)
+}
+\arguments{
+\item{x}{A vector of x (or longitude) coordinates.}
+
+\item{y}{A vector of y (or latitude) coordinates.}
+
+\item{hs}{A vector of headings (angle in rads).}
+
+\item{geo}{A logical value indicating whether the locations are defined by
+geographic coordinates (pairs of longitude/latitude values). Default: FALSE.}
+}
+\value{
+A vector of the same length as x, y and hs representing the bearing
+angle to the nearest neighboring object for each object.
+}
+\description{
+Given the locations and headings of different objects,
+ this function determines the angle between the heading of each object
+ and the position to the nearest neighboring object (bearing angle).
+}
+\examples{
+x <- rnorm(25)
+y <- rnorm(25, sd = 3)
+hs <- rnorm(25, sd = 1)
+nnba(x, y, hs)
+
+}
+\seealso{
+\code{\link{pdist}}
+}
+\author{
+Simon Garnier, \email{garnier@njit.edu},
+ Marina Papadopoulou, \email{m.papadopoulou.rug@gmail.com}
+}
diff --git a/man/nnd.Rd b/man/nnd.Rd
index 2ad6c82..c64c955 100644
--- a/man/nnd.Rd
+++ b/man/nnd.Rd
@@ -11,15 +11,15 @@ nnd(x, y, geo = FALSE)
 
 \item{y}{A vector of y (or latitude) coordinates.}
 
-\item{geo}{A logical value indicating whether the locations are defined by 
+\item{geo}{A logical value indicating whether the locations are defined by
 geographic coordinates (pairs of longitude/latitude values). Default: FALSE.}
 }
 \value{
-A vector of the same length as x and y representing the distance to 
+A vector of the same length as x and y representing the distance to
  the nearest neighboring object for each object.
 }
 \description{
-Given the locations of different objects, this function 
+Given the locations of different objects, this function
  determines the distance of the nearest neighboring object to each object.
 }
 \examples{
diff --git a/man/pdist.Rd b/man/pdist.Rd
index 9c1d955..2d728bf 100644
--- a/man/pdist.Rd
+++ b/man/pdist.Rd
@@ -11,15 +11,15 @@ pdist(x, y, geo = FALSE)
 
 \item{y}{A vector of y (or latitude) coordinates.}
 
-\item{geo}{A logical value indicating whether the locations are defined by 
+\item{geo}{A logical value indicating whether the locations are defined by
 geographic coordinates (pairs of longitude/latitude values). Default: FALSE.}
 }
 \value{
-A square matrix representing pairwise distances between each possible 
+A square matrix representing pairwise distances between each possible
  pair of locations.
 }
 \description{
-Given a set of locations, this function computes the distances 
+Given a set of locations, this function computes the distances
  between each possible pair of locations.
 }
 \examples{
diff --git a/tests/testthat/test-nnba.R b/tests/testthat/test-nnba.R
new file mode 100644
index 0000000..93962ff
--- /dev/null
+++ b/tests/testthat/test-nnba.R
@@ -0,0 +1,13 @@
+test_that("bearing angle works", {
+  expect_equal(nnba(x = rep(1,4),
+                    y = c(1, 2.1, 3, 4.1), 
+                    hs = rep(pi/2, 4)
+                    ), 
+               c(0, 0, pi, pi))
+  
+  expect_equal(nnba(y = rep(1,4),
+                    x = c(1, 3, 4.1, 5), 
+                    hs = rep(pi/4, 4)
+                    ), 
+               c(pi/4, pi/4, pi/4, -3*pi/4))
+})