# Source grind.R and cube.R first!

model <- function(t, state, parms) {
  with(as.list(c(state,parms)), {
    dR <- r*R*(1-R) + s - d*R - c1*N1*R - c2*N2*R 
    dN1 <- N1*(b1*c1*R/(h1+c1*R) - d1)
    dN2 <- N2*(b2*c2*R/(h2+c2*R) - d2)
    return(list(c(dR, dN1, dN2)))  
  }) 
}
  
s <- c(R=1, N1=0.01,N2=0.01)

size<-5#inch

pdf("compC.pdf",width=size,height=size)
par(mar=c(0,0,0,0))
p <- c(r=0,s=1,d=1,b1=2,b2=2,d1=1,d2=1,c1=1,c2=2,h1=0.2,h2=0.8)
f1 <- newton(c(R=0.2,N1=4,N2=0))
f2 <- newton(c(R=0.4,N1=0,N2=0.75))
f3 <- newton(c(R=1,N1=0,N2=0))
cube(xmax=1,ymax=5,zmax=3,theta=30)
cube(shade="yz",show="R",add=TRUE)
points3D(x=f1[1],y=f1[2],z=f1[3],add=TRUE,pch=19)
points3D(x=f2[1],y=f2[2],z=f2[3],add=TRUE,pch=1)
points3D(x=f3[1],y=f3[2],z=f3[3],add=TRUE,pch=1)
dev.off()

pdf("compD.pdf",width=size,height=size)
par(mar=c(0,0,0,0))
p <- c(r=1,s=0,d=0,b1=2,b2=2,d1=1,d2=1,c1=4,c2=4,h1=0.8,h2=1.4)
f1 <- newton(c(R=0.2, N1=0.2,N2=0))
f2 <- newton(c(R=0.35, N1=0,N2=0.16))
f3 <- newton(c(R=1,N1=0,N2=0))
cube(xmax=1,ymax=0.25,zmax=0.25,theta=30)
cube(shade="xyz",show="R",add=TRUE)
points3D(x=f1[1],y=f1[2],z=f1[3],add=TRUE,pch=19)
points3D(x=f2[1],y=f2[2],z=f2[3],add=TRUE,pch=1)
points3D(x=f3[1],y=f3[2],z=f3[3],add=TRUE,pch=1)
points3D(x=0,y=0,z=0,add=TRUE,pch=1)
dev.off()