# Estimating the pair and mark correlation functions (Illian et al., 2008), Ap - 24.10.2013. Last updated on 24.10.2013.

rm(list = ls())										

library(Rcpp)

cppFile <- "/Users/pommeren/Dropbox/Rcpp/"

sourceCpp(paste(cppFile, "PairMarkCorrelation.cpp", sep = ""))


# Load data
dFile <- paste(cppFile, "PlotB.txt", sep = "")
TreeList <- read.table(dFile, header = TRUE)

# Settings for mingling and mark differentiation functions
m <- 200 # number of steps
dd <- 0.25 # step width
h <- 2.0 # bandwidth
m * dd
kernel <- 0 # 0: Epanechnikov, 1: box kernel

# Check coordinates
min(TreeList$x)
min(TreeList$y)
max(TreeList$x)
max(TreeList$y)

# Observation window
xmax <- 100 # width
ymax <- 100 # length

mdbh <- mean(TreeList$dbh)

# Calculate pair and mark correlation functions for observed data
gkl <- calcPairMark(TreeList, xmax, ymax, m, dd, h, mdbh, kernel)


# Calculate minimum distance and cut functions at that point
(dm <- calcMinDistance(TreeList$x, TreeList$y))
gkl <- gkl[gkl$r > dm, ]
gkl$r[1]


par(mfrow=c(2,2))
# pdf(file = paste(cppFile, "PairCorrelation.pdf", sep = ""))
par(mar = c(2.5, 4, 1, 1))
plot(gkl$r, gkl$g, las = 1, main = "",  xlab = "", ylab = "", xlim = c(0, 50), ylim = c(0, 2), cex.axis = 1.7, type = "l") 
curve(1 + x * 0, from = gkl$r[1], to = 50, lty = 2, add = TRUE)
box(lwd = 2)
# dev.off()

# pdf(file = paste(cppFile, "L.pdf", sep = ""))
par(mar = c(2.5, 3, 1, 1))
plot(gkl$r, gkl$L, las = 1, main = "",  xlab = "", ylab = "", xlim = c(0, 50), ylim = c(0, 50), lty = 1, type = "l",
cex.axis = 1.7) 
curve(1 * x, from = gkl$r[1], to = 50, lty = 2, add = TRUE)
box(lwd = 2)
# dev.off()

# pdf(file = paste(cppFile, "MarkCorrelation.pdf", sep = ""))
par(mar = c(2.5, 4, 1, 1))
plot(gkl$r, gkl$k, las = 1, main = "",  xlab = "", ylab = "", xlim = c(0, 50), ylim = c(0, 2), cex.axis = 1.7, type = "l")
curve(1 + x * 0, from = gkl$r[1], to = 50, lty = 2, add = TRUE)
box(lwd = 2)
# dev.off()

# pdf(file = paste(cppFile, "L_k.pdf", sep = ""))
par(mar = c(2.5, 3, 1, 1))
plot(gkl$r, gkl$Lk, las = 1, main = "",  xlab = "", ylab = "", xlim = c(0, 50), ylim = c(0, 50), lty = 1, type = "l",
cex.axis = 1.7) 
curve(1 * x, from = gkl$r[1], to = 50, lty = 2, add = TRUE)
box(lwd = 2)
# dev.off()

rm(TreeList)