###################################################################################################################
# Title: r-script to compute population trends for different species at a common diameter as a basis for comparison 
# Notes: eliminates differences in tree size distribution as a confounding variable in species comparisons  
# Author: Christian O. Marks (The Nature Conservancy)
# Date: January 21, 2015
####################################################################################################################

# read in the data:
setwd("/Users/cmarks/Documents/Data/demography")
sppnames <- read.csv(file="tblSpecies.csv",header=T)
spplist <- read.csv(file="spplist.csv",header=T)
spptable <- merge(spplist,sppnames,by.x="Species", by.y="Species_code")#includes common species names, maximum measured size, and replication
#read in demographic function parameters:
growth <- read.csv(file="GrowthCoef.csv",header=T)
mort <- read.csv(file="MortCoef.csv",header=T)
density <- read.csv(file="DensCoef.csv",header=T)

###################################################################################################################
# define demographic functions:

neg.exponential.function <- function(diam,a,b) { a*exp(-b*diam) }

neg.exponential.function.prime <- function(diam,a,b) { -b*a*exp(-b*diam) }

neg.power.function <- function(diam,a,b) { a*(diam^(-b)) }

neg.power.function.prime <- function(diam,a,b) { -b*a*(diam^(-b-1)) }

power.function <- function(diam,a,b) { a/10*(diam^b) } # growth function

power.function.prime <- function(diam,a,b) { b*a/10*(diam^(b-1)) }

double.exponential.function <- function(diam,a,b,c,d)
{   1 - ((1-a*exp(b*diam/100)) * (exp(c*((diam/100)^d))))}

weibull.function <- function(diam,a,b,c)
{    (a*(diam)^(c-1))*(exp(-b*(diam)^c))  }

weibull.function.prime <- function(diam,a,b,c)
{    (a*(c-1)*(diam)^(c-2))*(exp(-b*(diam)^c)) + (a*(diam)^(c-1))*(exp(-b*(diam)^c))*(-b*c*(diam)^(c-1))  }

basemort.negexp <- function(diam,a1,b1,a2,b2) 
{    - (power.function(diam,a1,b1)) / (neg.exponential.function(diam,a2,b2)) * (neg.exponential.function.prime(diam,a2,b2)) - (power.function.prime(diam,a1,b1))  } 

basemort.negpow <- function(diam,a1,b1,a2,b2) 
{    - (power.function(diam,a1,b1)) / (neg.power.function(diam,a2,b2)) * (neg.power.function.prime(diam,a2,b2)) - (power.function.prime(diam,a1,b1))  } 

basemort.weibull <- function(diam,a1,b1,a2,b2,c2) 
{    - (power.function(diam,a1,b1)) / (weibull.function(diam,a2,b2,c2)) * (weibull.function.prime(diam,a2,b2,c2)) - (power.function.prime(diam,a1,b1))  } 

delta.term2.negexp <- function(diam,a1,b1,a2,b2) 
{    - (power.function(diam,a1,b1)) * (neg.exponential.function.prime(diam,a2,b2)) - (neg.exponential.function(diam,a2,b2)) * (power.function.prime(diam,a1,b1))  } 

delta.term2.negpow <- function(diam,a1,b1,a2,b2) 
{    - (power.function(diam,a1,b1)) * (neg.power.function.prime(diam,a2,b2)) - (neg.power.function(diam,a2,b2)) * (power.function.prime(diam,a1,b1))  } 

delta.term2.weibull <- function(diam,a1,b1,a2,b2,c2) 
{    - (power.function(diam,a1,b1)) * (weibull.function.prime(diam,a2,b2,c2)) - (weibull.function(diam,a2,b2,c2)) * (power.function.prime(diam,a1,b1))  } 

#####################################################################################################
# compute population trends for different common species at a given dbh:
dbh = 30 # dbh in cm to be used as the basis of comparison

#take data for common species individually:
common.species <- subset(spptable,spptable$N>=50 & spptable$maxdiam>=dbh)
sppnum<-dim(common.species)[1]

res.table <- data.frame(common.species, trend=0)

for (i in 1:sppnum) {

sp1 <- as.character(common.species[i,"Species"])
Dtype <- subset(density$Function,density$Species==sp1)
x <- dbh # dbh to be used as the basis of comparison
a1 <- subset(growth$a,growth$Species==sp1)
b1 <- subset(growth$b,growth$Species==sp1)
a2 <- subset(density$a,density$Species==sp1)
b2 <- subset(density$b,density$Species==sp1)
a <- subset(mort$a, mort$Species==sp1)
b <- subset(mort$b, mort$Species==sp1)
c <- subset(mort$c, mort$Species==sp1)
d <- subset(mort$d, mort$Species==sp1)

# modeled trend in population structure as is only:
if(Dtype=="exponential"){ 
	y0 <- -double.exponential.function(x,a,b,c,d) *  neg.exponential.function(x,a2,b2) + delta.term2.negexp(x,a1,b1,a2,b2) 
} else { 
	y0 <- -double.exponential.function(x,a,b,c,d) *  neg.power.function(x,a2,b2) + delta.term2.negpow(x,a1,b1,a2,b2) 
}
res.table[i,8]=y0
}

# save species results table:
write.csv(res.table, file="Population trends at 30 cm dbh.csv")