# Estimating the mean of angles (Assuncao, 1994), Ap - 26.07.2007. Translated to R on 19.09.2012. Updated on 13.02.2013.

rm(list = ls())										

#
# Calculates Euclidean distance (with periodic boundary conditions).
#
getEuclideanDistance <- function(xmax, ymax, x1, y1, x2, y2) {
  dx <- abs(x1 - x2)
  dy <- abs(y1 - y2)
  dx <- min(dx, xmax - dx)
  dy <- min(dy, ymax - dy)
  dz <- sqrt(dx * dx + dy * dy)
  return(dz)
}

#
# Calculates Euclidean distance (without periodic boundary conditions).
#
getSimpleEuclideanDistance <- function(x1, y1, x2, y2) {
  dx <- abs(x2 - x1)
  dy <- abs(y2 - y1)
  dz <- dx^2 + dy^2
  return(dz^0.5)
}

#
# Calculates the number of grid points.
#
calcGridPoints <- function(xmax, ymax, gridwidth)
  return((xmax - gridwidth) * (ymax - gridwidth) * (1 / (gridwidth * gridwidth)))

#
# Calculates the grid point coordinates. xrow -> x direction, ycolumn -> y 
# direction, numbers start in the bottom left corner and progress in x direction.
#
 calcGridPointCoordinates <- function(xmax, gridwidth, number) {
  offset <- 3.5
  number.per.row <- floor((xmax / gridwidth) - gridwidth)
  xnumber <- number.per.row
  xrow <- 1
  if(number > number.per.row)
  #repeat {
  while(number > xnumber) {
  	xnumber <- xnumber + number.per.row
  	xrow <- xrow + 1
  	# if(number <= xnumber) break
  } 
  else
    xrow <- 1
  ycolumn <- number - (number.per.row * (xrow - 1))
  x <- offset + ycolumn * gridwidth
  y <- offset + xrow * gridwidth
  new.point <- data.frame(cbind(x, y))
  return(new.point) 	
}

#
# This function identifies the nerest tree neighbour of a grid point.
#
findNearestNeighbourOfXY <- function(TreeList, x, y, ind, xmax, ymax) {
  min.distance <- 1E20
  index <- 0
  for (i in 1 : length(TreeList$Number)) {
    if(i != ind) {
      actualDistance <- getEuclideanDistance(xmax, ymax, x, y, TreeList$x[i], TreeList$y[i])
      if(actualDistance < min.distance) {
      	min.distance <- actualDistance
      	index <- i
      }	
    }
  }
  return(index)
}

#
# Calculates angle between one of the trees and the test location in relation to 
# zero/the x axis. Includes periodic boundary conditions.
#
calcAngle <- function(TreeList, index, x, y, xmax, ymax){
  dx <- TreeList$x[index] - x
  dy <- TreeList$y[index] - y
  r <- dx^2 + dy^2
  xx <- TreeList$x[index]
  yy <- TreeList$y[index]
  if(dx > xmax * 0.5) {
  	w <- (dx - xmax)^2
  	if(w + dy^2 < r)
  	  xx <- TreeList$x[index] - xmax  
  }
  if(dx > -xmax * 0.5) {
  	w <- (dx + xmax)^2
  	if(w + dy^2 < r)
  	  xx <- TreeList$x[index] + xmax  
  }
  if(dy > ymax * 0.5) {
  	w <- (dy - ymax)^2
  	if(w + dx^2 < r)
  	  yy <- TreeList$y[index] - ymax  
  }
  if(dy > -ymax * 0.5) {
  	w <- (dy + ymax)^2
  	if(w + dx^2 < r)
  	  yy <- TreeList$y[index] + ymax  
  }
  dx <- xx - x
  dy <- yy - y
  i <- 0
  if((dx > 0) && (dy > 0))
    i <- 0
  if((dx <= 0) && (dy <= 0))
    i <- 2
  if((dx <= 0) && (dy > 0))
    i <- 1
  if((dx > 0) && (dy <= 0))
    i <- 3
  if(abs(dx) < 0.00000001)
    dx <- 0.00000001  
  if(abs(dy) < 0.00000001)
    dy <- 0.00000001  
  if ((i == 0) || (i == 2))
    angle <- atan(dy / dx)
  else
    angle <- atan(-dx / dy)
  angle <- angle + i * pi * 0.5
  angle <- angle * 180 / pi
  return(angle)
}

#
# Calculates the angle between two trees in relation to the test location.
#
calcAngleTwoPoints <- function(angle1, angle2) {
  if(angle1 > angle2) {
  	finalAngle <- angle1
  	angle1 <- angle2
  	angle2 <- finalAngle
  }
  finalAngle <- angle2 - angle1
  if(finalAngle > 180) {
  	if(angle1 > 180)
  	  finalAngle <- 360 - angle1 + angle2
  	if(angle2 > 180)
  	  finalAngle <- 360 - angle2 + angle1
  }
  return(finalAngle)
}

#
# MeanOfAngles. Index calculation.
#
calcMeanOfAngles <- function(TreeList, xmax, ymax, gridwidth) {
  offset <- 7
  number.of.points <- round(calcGridPoints(xmax - offset, ymax - offset, gridwidth))
  cat("Number of points: ", number.of.points, "\n")
  xmean <- 0
  points <- 0
  for (i in 1 : number.of.points) {
    xpoint <- calcGridPointCoordinates(xmax - offset, gridwidth, i)
    index1 <- findNearestNeighbourOfXY(TreeList, xpoint$x, xpoint$y, 0, xmax, ymax)
    index2 <- findNearestNeighbourOfXY(TreeList, xpoint$x, xpoint$y, index1, xmax, ymax)
    angle1 <- calcAngle(TreeList, index1, xpoint$x, xpoint$y, xmax, ymax)
    angle2 <- calcAngle(TreeList, index2, xpoint$x, xpoint$y, xmax, ymax)
    xpoint$angle <- calcAngleTwoPoints(angle1, angle2)
    xmean <- xmean + xpoint$angle
    if(i == 1)
      points <- xpoint[c("x","y","angle")]
    else 
      points <- rbind(points, xpoint)
  }
  xmean <- xmean / number.of.points
  cat("Mean of angles: ", xmean, "\n")
  return(points)
}

# Implementation

# Settings
# xpath <- "/Users/afs80b/Dropbox/Rcpp/"
xpath <- "/Users/arnepommerening/Dropbox/Rcpp/"
# xpath <- "/Rcpp/"

xFile <- paste(xpath, "Clg6.txt", sep = "")
xmax <- 98 
ymax <- 107

# Load data
sdata <- read.table(xFile, header = T)
dim(sdata)										
names(sdata)
length(sdata)

pointList <- calcMeanOfAngles(sdata, xmax, ymax, 0.25)			
pointList$weight <- 1 / length(pointList$angle) 
d <- density(pointList$angle, weights = pointList$weight, kernel = "epanechnikov", bw = 15, n = length(pointList$angle)) #, from = 0, to = 180)
par(mar = c(2, 5, 1, 1))
# par(new = T)
plot(d, xlab = "", ylab = "", main = "", las = 1, axes = FALSE, lwd = 2, col = "red", xlim = c(0, 210), ylim = c(0, 0.010)) # col = "blue", ylim = c(0, 2.0), xlim = c(0, 1),
axis(side = 1, lwd = 2, las = 1, cex.axis = 1.7)
axis(side = 2, lwd = 2, las = 1, cex.axis = 1.7)
box(lwd = 2)

max(pointList$angle)
sum(pointList$weight)