R/egg.geodesic.R

In egg: Distances and areas on the surface of an egg

#' Estimate the geodesic
#'
#' Estimates the geodesic between the two input points
#'
#' @param egg an object of class \code{egg.fit}
#' @param coords the spherical coordinates of two points on the egg surface, as obtained from \code{egg.coords}
#' @param n the number of iterative rotations between \code{p} and \code{q}
#' @param plot logical, whether or not to plot the geodesic in the current plot window
#' @param return.points logical, should a vector of points along the geodesic be returned
#' @param ... additional parameters passed to the \code{\link[graphics]{lines}} function
#'
#' @return The geodesic between two points
#'
#' @examples
#' stop()
#'
#' @export
egg.geodesic = function(egg, coords, n=1000, plot=FALSE, return.points=FALSE, ...) {
  a = egg$a
  b = egg$b
  c = egg$c
  h = egg$length
  p = coords[1,]
  q = coords[2,]
  p.cart = norm(sph2cart(p))
  q.cart = norm(sph2cart(q))
  alpha = angle.between(p.cart, q.cart)
  points = matrix(sph2cart(p), nrow=1, ncol=3)
  dist = 0
  if(alpha > 1e-6) {
    points = rbind(points, matrix(0, nrow=n, ncol=3))
    for(i in 2:(n+1)) {
      points[i,] = rotate.towards(p.cart, q.cart, (alpha/n)*(i-1))
      phi = acos(points[i,3]/1)
      r = get.r(phi, egg)
      points[i,] = points[i,]*r
      dist = dist + magnitude(points[i,] - points[(i-1),])
    }
  } else {
    points = matrix(rbind(points, sph2cart(q)), ncol=3)
  }
  if(plot) {
    lines(points[,1]+egg$mid.point[1], points[,3]+egg$mid.point[2], ...)
  }
  if(return.points) {
    for(i in 1:n) {
      points[i,] = cart2sph(points[i,])
    }
    return(list(dist=dist, points=points))
  } else {
    return(dist)
  }
}