#######################################################################################
#
#   Program: CalcHR.R
#   Author:  John Fieberg 
#       
#
#   Description:    Function for calculating kernel HR estimates for a particular contour  
#
#######################################################################################
 
 
#-------------------- Function to calculate HR and 95% and 50% contours --------

calcHR<-function(est, p=c(0.5, 0.95)){

#  p = quantiles for estimation and display purposes

 
 x<-est$x1
 y<-est$x2
 z<-est$fhat
 
# Calculate HR size
#  Get size of grid
dx<-x[2]-x[1]
dy<-y[2]-y[1]

# Estimate total volume and make sure =1 or close to
totp<-sum(z)*dx*dy

if(round(totp,4)<0.95){
    print("ERROR, total probability not equal to 1")
    print(paste("totp = ",totp))}
else z<-z/totp # to ensure integrates to 1

#  Get minimum number of squares to encompass pHR% probability
nps<-length(p)
zsort <- sort(as.vector(z)) 
pt <- cumsum(zsort)/sum(zsort)


#  p's for plotting and number of cells for calculating HR
pplot<-matrix(0,nps,1)
ncells<-matrix(0,nps,1)
for(i in 1:nps){
  pplot[i]<-min(zsort[pt>=(1-p[i])])
  ncells[i]<-length(pt[pt>=(1-p[i])])
}

# Home range area
HR<-dx*dy*ncells
out<-list(HR, pplot, totp)
names(out)<-c("HR", "ps", "totp")
return(out)
}