model <- function(t, state, parms) {
  with(as.list(c(state,parms)), {
    dR1 <- s1 - d1*R1 - c11*N1*R1 - c21*N2*R1
    dR2 <- s2 - d2*R2 - c12*N1*R2 - c22*N2*R2
    dN1 <- (a11*c11*R1 + a12*c12*R2 - delta1)*N1
    dN2 <- (a21*c21*R1 + a22*c22*R2 - delta2)*N2
    dN3 <- (a31*c31*R1 + a32*c32*R2 - delta3)*N3
    return(list(c(dR1, dR2, dN1, dN2, dN3)))  
  }) 
}  

qss <- function(t, state, parms) {
  with(as.list(c(state,parms)), {
    R1 <- s1/(d1 + c11*N1 + c21*N2)
    R2 <- s2/(d2 + c12*N1 + c22*N2)
    dN1 <- (a11*c11*R1 + a12*c12*R2 - delta1)*N1
    dN2 <- (a21*c21*R1 + a22*c22*R2 - delta2)*N2
    dN3 <- (a31*c31*R1 + a32*c32*R2 - delta3)*N3
    return(list(c(dN1, dN2, dN3)))  
  }) 
}  

vectors <- function(parms,size1=1,size2=0.5) {
  with(as.list(parms), {
    R2 <- as.numeric((a21*c21*delta1/(a11*c11)-delta2)/(a12*c12*a21*c21/(a11*c11) - a22*c22))
    R1 <- as.numeric((delta1 - a12*c12*R2)/(a11*c11))
    # v1 = (c11R1 c12R2) & v2 = (c21R1 c22R2) & v3 = (s1-d1R1 s2-d2R2)
    arrows(R1,R2,R1-size1*c11*R1,R2-size1*c12*R2,col=colors[3],lwd=2,length=0.1)
    arrows(R1,R2,R1-size1*c21*R1,R2-size1*c22*R2,col=colors[4],lwd=2,length=0.1)
    dx <- model(0,c(R1=R1,R2=R2,N1=0,N2=0,N3=0),parms)[[1]]
    arrows(R1,R2,R1+size2*dx[1],R2+size2*dx[2],lwd=2,length=0.1)
    return(c(R1,R2))
  })
}

s <- c(R1=1,R2=1,N1=0.1,N2=0.1,N3=0.1)

pR <- c(s1=1.5,s2=1.7,d1=1,d2=1)
pN <- c(delta1=0.5,delta2=0.5,delta3=0.5)
pC <- c(c11=0.5,c22=0.5,c12=0.25,c21=0.25,c31=0.3,c32=0.3)
pA <- c(a11=1,a22=1,a12=1,a21=1,a31=1,a32=1)
#pC <- c(c11=0.5,c22=0.5,c12=0.5,c21=0.5)  # Same diet
p <- c(pR,pN,pC,pA)
initR1 <- as.numeric(p["s1"]/p["d1"]) 
initR2 <- as.numeric(p["s2"]/p["d2"])
plane(show=c("N1","N2"),xmax=2,ymax=2,zero=FALSE)
vectors(p)
newton(c(R1=initR1,R2=initR2,N1=0,N2=0,N3=0),plot=TRUE)
run(state=c(R1=initR1,R2=initR2,N1=0.1,N2=0.1,N3=0.1),traject=TRUE)
f <- newton(run(200,state=c(R1=initR1,R2=initR2,N1=0.1,N2=0.1,N3=0.1),timeplot=FALSE),plot=TRUE)

plane(odes=qss,state=c(N1=0,N2=0,N3=0),xmax=5,ymax=5)


# Make Figures
# size<-5#inch
# 
# pdf("tilmanLVA.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")
# plane(show=c("N1","N2"),zero=FALSE,xmax=2,ymax=2)
# f <- newton(run(200,state=c(R1=initR1,R2=initR2,N1=0.1,N2=0.1,N3=0.1),timeplot=FALSE),plot=TRUE)
# dev.off()
# 
# pdf("tilmanLVB.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")
# plane(show=c("N1","N2","N3"),zero=FALSE,xmax=2,ymax=2)
# f <- newton(run(200,state=c(R1=initR1,R2=initR2,N1=0.1,N2=0.1,N3=0.1),timeplot=FALSE),plot=TRUE)
# dev.off()
# 
# pdf("tilmanLVC.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")
# pC <- c(c11=0.5,c22=0.5,c12=0.25,c21=0.25,c31=0.3,c32=0.3)
# p <- c(pR,pN,pC,pA)
# plane(show=c("N1","N2"),zero=FALSE,xmax=2,ymax=2)
# vectors(p)
# newton(c(R1=initR1,R2=initR2,N1=0,N2=0,N3=0),plot=TRUE)
# f <- newton(run(200,state=c(R1=initR1,R2=initR2,N1=0.1,N2=0.1,N3=0.1),timeplot=FALSE),plot=TRUE)
# dev.off()
# 
# pdf("tilmanLVD.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")
# pC <- c(c11=0.4,c22=0.5,c12=0.25,c21=0.25,c31=0.3,c32=0.3)
# p <- c(pR,pN,pC,pA)
# plane(show=c("N1","N2"),zero=FALSE,xmax=2,ymax=2)
# vectors(p)
# newton(c(R1=initR1,R2=initR2,N1=0,N2=0,N3=0),plot=TRUE)
# dev.off()
# 
# pdf("tilmanLVE.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")
# pC <- c(c11=0.5,c22=0.5,c12=0.25,c21=0.25,c31=0.3,c32=0.3)
# pA <- c(a11=0.5,a22=0.5,a12=1.25,a21=1.25,a31=1,a32=1)
# p <- c(pR,pN,pC,pA)
# plane(show=c("N1","N2"),zero=FALSE,xmax=2,ymax=2)
# vectors(p)
# f <- newton(run(10,state=c(R1=initR1,R2=initR2,N1=0.1,N2=0.1,N3=0.1),timeplot=FALSE),plot=TRUE)
# newton(c(R1=initR1,R2=initR2,N1=0,N2=0,N3=0),plot=TRUE)
# dev.off()
# 
# pdf("tilmanLVF.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")
# pC <- c(c11=0.5,c22=0.5,c12=0.25,c21=0.25,c31=0,c32=0)
# pA <- c(a11=0.5,a22=0.5,a12=1.25,a21=1.25,a31=1,a32=1)
# p <- c(pR,pN,pC,pA)
# plane(xmax=2.5,ymax=2.3,odes=qss,state=c(N1=1,N2=1,N3=0),portrait=TRUE,tmax=200)
# newton(state=c(N1=2,N2=0,N3=0),odes=qss,plot=TRUE)
# newton(state=c(N1=0,N2=2,N3=0),odes=qss,plot=TRUE)
# newton(state=c(N1=1,N2=1,N3=0),odes=qss,plot=TRUE)
# dev.off()