model <- function(t, state, parms) {
  with(as.list(c(state,parms)), {
    H <- h/(1-0.5*gamma)
    f <- 2*R/(H + R + sqrt((H+R)^2 - 4*gamma*H*R));   
    dR <- r*R*(1 - R/K) - a*N*f
    dN <- c*a*N*f - delta*N
    return(list(c(dR, dN)))  
  }) 
}  

reverse <- function(t, state, parms) {
  with(as.list(c(state,parms)), {
    H <- h/(1-0.5*gamma)
    f <- 2*R/(H + R + sqrt((H+R)^2 - 4*gamma*H*R));   
    dR <- r*R*(1 - R/K) - a*N*f
    dN <- c*a*N*f - delta*N
    return(list(c(-dR, -dN)))  
  }) 
} 

H <- 0.1; gamma <- 0
curve(2*R/(H+R+sqrt((H+R)^2-4*gamma*H*R)),from=0,to=1,ylim=c(0,1),xname="R",col="red",lwd=2)
for (gamma in c(0.5,0.75,0.99,1)) 
  curve(2*R/(H+R+sqrt((H+R)^2-4*gamma*H*R)),xname="R",col="red",lwd=2,add=TRUE)

h <- 0.1; gamma <- 0; H <- h/(1-0.5*gamma)
curve(2*R/(H+R+sqrt((H+R)^2-4*gamma*H*R)),from=0,to=1,ylim=c(0,1),xname="R",col="red",lwd=2)
for (gamma in c(0.5,0.75,0.99,1)) {
  H <- h/(1-0.5*gamma)
  curve(2*R/(H+R+sqrt((H+R)^2-4*gamma*H*R)),xname="R",col="red",lwd=2,add=TRUE)
}

s <- c(R=1,N=0.01)
p <- c(r=1,K=1,h=0.1,a=0.5,c=1,delta=0.3,gamma=0)
p["gamma"] <- 0.95;plane();run(tstep=0.1,traject=TRUE)
p["gamma"] <- 0.99;plane();run(tstep=0.1,traject=TRUE)



# Make figures
# size <- 5 #inch
# 
# pdf("hyperA.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")
# h <- 0.1; gamma <- 0; H <- h/(1-0.5*gamma)
#   curve(2*R/(H+R+sqrt((H+R)^2-4*gamma*H*R)),from=0,to=1,ylim=c(0,1.1),xname="R",col="red",lwd=2)
# for (gamma in c(0.5,0.75,0.99,1)) {
#   H <- h/(1-0.5*gamma)
#   curve(2*R/(H+R+sqrt((H+R)^2-4*gamma*H*R)),xname="R",col="red",lwd=2,add=TRUE)
# }
# dev.off()
# 
# pdf("hyperB.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")
# s <- c(R=1,N=0.01)
# p <- c(r=1,K=1,h=0.1,a=0.5,c=1,delta=0.3,gamma=0)
# p["gamma"] <- 0.95;plane();run(tstep=0.1,traject=TRUE)
# dev.off()
# 
# pdf("hyperC.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")
# s <- c(R=1,N=0.01)
# p <- c(r=1,K=1,h=0.1,a=0.5,c=1,delta=0.3,gamma=0)
# p["gamma"] <- 0.99;plane();run(tstep=0.1,traject=TRUE)
# dev.off()
# 
# pdf("hyperD.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(tstep=0.1,ymax=0.8,portrait=TRUE,grid=3,legend=FALSE)
# dev.off()
# 
# s <- c(R=0.15,N=0.2)
# s <- run(odes=reverse)
# pdf("hyperE.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(tstep=0.1,ymax=0.8,portrait=TRUE,grid=3,legend=FALSE)
# run(tmax=20,tstep=0.1,odes=reverse,traject=T,col="darkgreen")
# dev.off()