# Source grind.R and cube.R first!

model <- function(t, state, parms) {
  with(as.list(c(state,parms)), {
    dR1 <- r*R1*(1 - R1 - a12*R2) - a1*R1*N
    dR2 <- r*R2*(1 - R2 - a21*R1) - a2*R2*N
    dN  <- N*(c*a1*R1 + c*a2*R2 - 1)
    return(list(c(dR1, dR2, dN)))  
  }) 
}  

par(mar=c(3.6,2.6,1.6,0.2),mgp=c(1.5,0.5,0)) 
p <- c(r=1,a1=6,a12=1,a2=1,a21=1.5,c=0.5,delta=0.5)
s <- c(R1=0.1,R2=0.1,N=0.01)

plane()
run(state=c(R1=0.1,R2=0.1,N=0),traject=T)
cube(zmax=0.2)
run3d(zmax=0.2,add=TRUE)

# Make figures
# size <- 5
# pdf("rrna.pdf",width=size,height=size)
# par(mar=c(0.1,0.1,0.1,0.1),xaxt="n",yaxt="n",ann=F,xaxs="i",yaxs="i")
# p["a1"] <- 6;run(100,0.1,main="(a)")
# dev.off()
# 
# pdf("rrnb.pdf",width=size,height=size)
# par(mar=c(0.1,0.1,0.1,0.1),xaxt="n",yaxt="n",ann=F,xaxs="i",yaxs="i")
# p["a1"] <- 8;run(200,0.1,main="(b)",legend=FALSE)
# dev.off()
# 
# pdf("rrnc.pdf",width=size,height=size)
# par(mar=c(0.1,0.1,0.1,0.1),xaxt="n",yaxt="n",ann=F,xaxs="i",yaxs="i")
# p["a1"] <- 10;run(500,0.2,main="(c)",legend=FALSE)
# dev.off()
# 
# f <- s
# 
# pdf("rrn3a.pdf",width=size,height=size)
# par(mar=c(0,0,0,0)) # narrow margins
# p["a1"] <- 6
# data <- run(1000,table=TRUE,timeplot=FALSE)
# xmax=max(data$R1); ymax=max(data$R2); zmax=max(data$N)
# #cube(x=3,y=1,z=2,theta=130,xmax=0.5,ymax=1.5,zmax=1.5)
# f[1:3] <- as.numeric(data[1000,2:4])
# run3d(xmax=xmax,ymax=ymax,zmax=zmax,theta=20,phi=25,state=f,tmax=1000,tstep=0.1,col="red")
# dev.off()
# 
# pdf("rrn3b.pdf",width=size,height=size)
# par(mar=c(0,0,0,0)) # narrow margins
# p["a1"] <- 8
# data <- run(1000,table=TRUE,timeplot=FALSE)
# xmax=max(data$R1); ymax=max(data$R2); zmax=max(data$N)
# #cube(x=3,y=1,z=2,theta=130,xmax=0.5,ymax=1.5,zmax=1.5)
# f[1:3] <- as.numeric(data[1000,2:4])
# run3d(xmax=xmax,ymax=ymax,zmax=zmax,theta=20,phi=25,state=f,tmax=1000,tstep=0.1,col="red")
# dev.off()
# 
# pdf("rrn3c.pdf",width=size,height=size)
# par(mar=c(0,0,0,0)) # narrow margins
# p["a1"] <- 10
# data <- run(1000,table=TRUE,timeplot=FALSE)
# xmax=max(data$R1); ymax=max(data$R2); zmax=max(data$N)
# #cube(x=3,y=1,z=2,theta=130,xmax=0.5,ymax=1.5,zmax=1.5)
# f[1:3] <- as.numeric(data[1000,2:4])
# run3d(xmax=xmax,ymax=ymax,zmax=zmax,theta=20,phi=25,state=f,tmax=5000,tstep=0.1,col="red")
# dev.off()

# Use rgl
# install.packages("rgl")
# library(rgl)
# p["a1"] <- 10;data <- run(5000,tstep=0.1,table=T)
# plot3d(x=data$R1,y=data$R2,z=data$N,xlab="R1",ylab="R2",zlab="N",type="l")