##### plots of effects of intensity and prevalence on transmission

fs=0.28
par(mfrow=c(1,2),mai=c(2,0.85,2,0.85))
factor=0.0
#set up a vector of numbers on the range of 0 to 100.
num = (0:100)
# prelavenvce is the probability that a contact is infectious = I/N
prev=num/100
# c is the contact rate in a population, with density dependence c = k(N/A), where k is constant, and A is the area occupied.

Nlarge=10000
Nmedium=1000
Nsmall=100
A=10000
clarge=Nlarge/A
cmedium=Nmedium/A
csmall=Nsmall/A
# nu is the infectivity of a contact
nu = .1
# infectivity values in a vector 
nua=num/100

# Based on Begon et al 2002 and Antonovics et al 1995 (equations 1-5) for 
# density-dependent transmission 
beta1=1-exp(-(nu*clarge*Nsmall*prev))
beta2=1-exp(-(nu*cmedium*Nsmall*prev))
beta3=1-exp(-(nu*csmall*Nsmall*prev))
plot(prev,beta1,xlab = "Prevalence within populations",ylab="P(infection)",csi=fs,xlim=c(0,1),pch=" ")
title("Density-dependent: N = 100; c = 0.01, 0.1, 1; nu = 0.1",cex=(fs*2))
lines(prev,beta1,col=3)
lines(prev,beta2,col=2)
lines(prev,beta3,col=4)

# Anonovic et al (1995) for frequency-dependent transmission, where n is the contact rate is determined by the ratio of
# seartch time for hosts (T) to the handling time for infection (Th) (Antonovics et al Eq 6-7).  This rate (n) can then scale 1 (T = Th)
# to infinity (Th=0).  Using the same Poisson model for encounters, this becomes:   
# 
n=1
beta4=1-exp(-(nu*n*prev)/(1-prev))
plot(prev,beta4,xlab = "Prevalence within populations",ylab="",csi=fs,xlim=c(0,1),pch=" ")
title("Frequency-dependent, n=1, 10, 100; nu=0.1",cex=(fs*2))
lines(prev,beta4,col=6)

n=10
beta5=1-exp(-(nu*n*prev)/(1-prev))
lines(prev,beta5,col=7)

n=100
beta6=1-exp(-(nu*n*prev)/(1-prev))
lines(prev,beta6,col=1)

##############
# 3-D graphs #
##############

#  set up matrix of values to graph, C for contact rate, N for population size
C=num/100
N=num


threed <- data.frame(p=seq(from=min(prev), to=max(prev), length=25), 
   c1=seq(from=min(C), to=max(C), length=25),
	nu1=seq(from=min(nua), to=max(nua), length=25),
	n1=seq(from=min(N),to=max(N), length=25))


################################
#  set number 1: contact rates #
################################

par(mfrow=c(1,2),mai=c(2,0.85,2,0.85))
# Density dependent graphing function for prevlance and contact rate for nu = 0.1 and population size = 100
nu=0.1
Popsize=100 
beta.dd <- function(x, y){
	1-exp(-(x*y*nu*Popsize))
	 }
#Density dependent graph of effects of prevlance and contact rate on probability of S->I
persp(threed[,1],threed[,2], outer(threed[,1],threed[,2], 
      FUN =  beta.dd), 
      xlab = "Prevalence", ylab = "Contact rate", zlab = "P(infection)", axes = T,ar=1,zlim=0,  
      box = T) 
#eye = c(0, -500, 3))
title("Density-dependent, N = 100, nu=0.1",cex=(fs*2))


#Frequency-dependent 3-D graphs

# Frequency-dependent graphing function for prevalence and contact rate on probability of S->I, nu=0.1
nu=0.1 
beta.fd <- function(x, y){
	1-exp(-(x*y*nu/(1-x)))
	 }
	
persp(threed[,1],threed[,4], outer(threed[,1],threed[,4],
      FUN =  beta.fd), 
      xlab = "Prevlance", ylab = "Contact rate", zlab = "P(infection)", axes = T,ar=1,zlim=0,  
      box = T) 
#eye = c(0, -500, 3))
title("Frequency-dependent, nu = 0.1",cex=(fs*2))


################################
#  set number 2: infectivity   #
################################


par(mfrow=c(1,2),mai=c(2,0.85,2,0.85))

# Density dependent graphing function for prevlance and infectivity for Population size = 100, contact rate = 0.01 (csmall)
csmall=0.01
Popsize=100 
beta.dd <- function(x, y){
	1-exp(-(x*y*csmall*Popsize))
	 }


# Density dependent graph of effects of prevlance and infectivity on probability of S->I
persp(threed[,1],threed[,3], outer(threed[,1],threed[,3], 
      FUN =  beta.dd), 
      xlab = "Prevalence", ylab = "Infectivity", zlab = "P(infection)", axes = T,ar=1,zlim=0,  
      box = T) 
#eye = c(0, -500, 3))
title("Density-dependent, N = 100, c=0.01",cex=(fs*2))


#Frequency-dependent 3-D graph

# function to test effects of prevalence and infectivity probability of S->I, with contact rate n = 1
n=1   
beta.fd <- function(x, y){
	1-exp(-(x*y*n/(1-x)))
	 }
	
persp(threed[,1],threed[,3], outer(threed[,1],threed[,3],
      FUN =  beta.fd), 
      xlab = "Prevlance", ylab = "Infectivity", zlab = "P(infection)", axes = T,ar=1,zlim=0,  
      box = T) 
#eye = c(0, -500, 3))
title("Frequency-dependent, contact = 1",cex=(fs*2))