model <- function(t, state, parms) {
  with(as.list(c(state,parms)), {
    dtA <- rA*A*(1 - (A+alpha*C)/kA)
    dtC <- rC*C*(1 - (C+beta*A)/kC)
    return(list(c(dtA,dtC)))  
  }) 
} 

p <- c(rA=1.1,rC=0.8,kA=105,kC=64,alpha=0,beta=0)
s <- c(A=2,C=2)

aurelia1 <- read.table("data/aurelia1.txt", header=TRUE)
aurelia2 <- read.table("data/aurelia2.txt", header=TRUE)
caudatum1 <- read.table("data/caudatum1.txt", header=TRUE)
caudatum2 <- read.table("data/caudatum2.txt", header=TRUE)

timePlot(aurelia1,draw=points)
timePlot(caudatum1,draw=points,add=TRUE,colMap=2)
run(18,add=TRUE)
                     
# Keep alpha=beta=0 and fit the two separate expts:

free <- c("rA","kA","rC","kC")
f1 <- fit(list(aurelia1,caudatum1),free=free,add=TRUE)
p[free] <- f1$par;p

# Next use these estimates to fit the two data sets from the mixed experiment:

timePlot(aurelia2,draw=points)
timePlot(caudatum2,draw=points,add=TRUE,colMap=2)
p["alpha"] <-1 ; p["beta"] <- 1
free <- c("alpha","beta")
f2 <- fit(list(aurelia2,caudatum2),free=free,add=TRUE)
p[free] <- f2$par; p

plane(xmax=110,ymax=100,eps=-0.01)
run(20,0.1,traject=TRUE)

# Finally fit all data together by fixing the initial conditions.
# Aurelia and Caudatum are absent from the first and second expt, respectively:

fixed <- list(A=c(2,0,2,2),C=c(0,2,2,2))
free <- c("rA","kA","rC","kC","alpha","beta")
f3 <- fit(list(aurelia1,caudatum1,aurelia2,caudatum2),free=free,fixed=fixed,add=TRUE)

f3 <- fit(list(aurelia1,caudatum1,aurelia2,caudatum2),free=free,fixed=fixed,add=TRUE,bootstrap=100)