# This program finds the best value of b for the Holling type 2 model.
#	r is the intrinsic aphid growth rate
#	q is the maximum predator consumption rate
#	P[t] is the number of predators for day t of the experiment
#	Nexpt is the vector of aphid populations
#		Nexpt must have one more element than P
#	Ntheory is the vector of populations from the model
#	b is the semisaturation constant times the area

# DATA

r <- 0.317
q <- 51.6
P <- c(2,2,2,2,5,5)
Nexpt <- c(894,1099,1303,1499,1654,1191,1015)
DEsteps <- 5

# INITIALIZATION

c <- q*P/r

# FUNCTIONS

Nprime <- function(N,b,r,c)
	{
	deriv <- r*N*(N+b-c)/(N+b)
	return(deriv)
	}

desolve <- function(N0,b,r,c,DEsteps)
	{
	dt <- 1/DEsteps
	N <- N0
	for (i in 1:DEsteps)
		{
		k1 <- Nprime(N,b,r,c)
		k2 <- Nprime(N+.5*k1*dt,b,r,c)
		k3 <- Nprime(N+.5*k2*dt,b,r,c)
		k4 <- Nprime(N+k3*dt,b,r,c)
		N <- N+(k1+2*k2+2*k3+k4)*dt/6
		}
	return(N)
	}

model <- function(N0,b,r,c,DEsteps,T)
	{
	result <- rep(0,T)
	result[1] <- N0
	for (t in 1:T)
		result[t+1] <- desolve(result[t],b,r,c[t],DEsteps)
	return(result)
	}

RSS <- function(Nexpt,b,r,c,DEsteps)
	{
	T <- length(Nexpt)-1
	N0 <- Nexpt[1]
	Ntheory <- model(N0,b,r,c,DEsteps,T)
	error <- Nexpt-Ntheory
	result <- c(sum(error^2),Ntheory)
	return(result)
	}

# MAIN PROGRAM

T <- length(Nexpt)-1
bvalues <- 10*(0:50)
F <- rep(0,length(bvalues))
for (j in 1:length(bvalues))
	{
	b <- bvalues[j]
	result <- RSS(Nexpt,b,r,c,DEsteps)
	F[j] <- result[1]
	}
j0 <- which.min(F)
b0 <- bvalues[j0]
result <- RSS(Nexpt,b0,r,c,DEsteps)
F0 <- result[1]
Ntheory <- result[-1]

# OUTPUT

#par(mfrow=c(1,2))

ymax <- 20*F0
xrange <- which(F<ymax)
plot(bvalues,F,ylim=c(0,ymax),xlim=c(bvalues[min(xrange)],bvalues[max(xrange)]),xlab="b",ylab="RSS")
plot(0:T,Nexpt,ylim=c(0,max(max(Ntheory),max(Nexpt))),col="green3",pch=15,xlab="t",ylab="N")
points(0:T,Ntheory,col="blue3",pch=18)
legend(0.5,3300,c("experiment","model"),col=c("green3","blue3"),pch=c(15,18))
lines(c(0,6),c(400,400),lwd=0.5)
for (xx in c(0,4,6))
	lines(c(xx,xx),c(400,500),lwd=0.5)
text(2,250,"2 coccinellids")
text(5,250,"5 coccinellids")

b0
F0