MaynardSmith <- function(r1,
                         r2,
                         K,
                         C,
                         x0=NULL,
                         T=30                         
                         ){
  write("\nCall Maynard Smith predator-prey model","")
  if (!(r1>1) | !(r2>1) | !(K>0) | !(C>0)){
    write(paste(" Error: parameters insufficient (r1=",r1," r2=",r2,
                " K=",K," C=",C,"\n",sep=""),"")
    NULL}
  else{if (length(x0)<2){x0 <- K*c(1,0.001)}
       if (T < 1){T <- 1
                  write(" warning: time of simulation set to 1","")}

       write(paste("\n Prey growth rate            r1 = ",r1,sep=""),"")
       write(paste(" Carrying cappacity for prey  K = ",K,sep=""),"")
       write(paste(" Predator growth rate        r2 = ",r2,sep=""),"")
       write(paste(" Efficiency of predation      C = ",C,sep=""),"")

       Equil <- array(0,c(3,2))
       Equil[2,1] <- K
       Equil[3,1] <- K/r2; Equil[3,2] <- K*(r1-1)*(r2-1)/(C*r2)

       Stab <- c(FALSE,FALSE,((r2 < 2) & (r2 > 3*(r1-1)/(r1+3))))
       Osc <- (r2 < 0.5*(1+sqrt(r1))) & (r2 > 2*(1-1/r1))
       SStab <- ""
       if (Stab[3]){SStab <- "asymptotically stable"
                    if (Osc){
                      SStab <- paste(SStab,", (monotone convergence)",sep="")}}
       else{if ((r2 > 2) | (r2 < 3*(r1-1)/(r1+3))){SStab <- "unstable"}}

       write("\n Equilibria: (0 , 0) - unstable","")
       write(paste("             (",K," , 0) - unstable",sep=""),"")
       write(paste("             (",Equil[3,1]," , ",Equil[3,2],")\n",
                   "                     - ",SStab,sep=""),"")

       t <- seq(0,floor(T))
       x <- array(NA,c(2,T+1))
       x[,1] <- x0
       for (i in 1:T){x[1,i+1] <- (r1-((r1-1)*x[1,i]+C*x[2,i])/K)*x[1,i]
                      x[2,i+1] <- r2*x[1,i]*x[2,i]/K}
       MaynardSmith <- list(time=t,prey=x[1,],predator=x[2,],
                            equilibria=Equil,stability=Stab,
                            type="Maynard Smith")}
}


PlotsPP <- function(PP,xylims=NULL,t.xy=TRUE,phase=TRUE,
                    ch.prey=NULL,col.prey="green",
                    ch.predator=NULL,col.predator="red",
                    ch.phase=19,col.phase="blue",line.phase=TRUE){

  if (t.xy & phase){split.screen(c(2,1))
                    screen(1)}
  if (t.xy){
    B <- max(c(PP$prey,PP$predator))
    plot(c(PP$time[1],PP$time[length(PP$time)]),1.1*c(0,B),
         type="n",bty="l",
         main=paste(PP$type,"Predator-Prey Model"),
         xlab="time",ylab="prey, predator")
    if (!(length(ch.predator)==0)){
      points(PP$time,PP$predator,type="p",pch=ch.predator,col=col.predator)}
    points(PP$time,PP$predator,type="l",col=col.predator)
    if (!(length(ch.prey)==0)){
      points(PP$time,PP$prey,type="p",pch=ch.prey,col=col.prey)}
    points(PP$time,PP$prey,type="l",col=col.prey)
    legend(PP$time[1],1.1*B,c("prey","predator"),col=c(col.prey,col.predator),
           lty=1,pch=c(ch.prey,ch.predator),horiz=TRUE,yjust=0.5,bty="n")
  }

  if (t.xy & phase){screen(2)}
  if (phase){
    if (length(xylims) == 0){xlm <- c(min(c(PP$prey,PP$equilibria[,1])),
                                      max(c(PP$prey,PP$equilibria[,1])))
                             ylm <- c(min(c(PP$predator,PP$equilibria[,2])),
                                      max(c(PP$predator,PP$equilibria[,2])))}
    else{xlm <- c(xylims[1],xylims[3])
         ylm <- c(xylims[2],xylims[4])}
  
  if (t.xy & phase){screen(2)}
    
    plot(PP$prey,PP$predator,type="p",pch=ch.phase,col=col.phase,bty="l",
         xlab="prey",ylab="predator",xlim=xlm,ylim=ylm)
    if(line.phase){points(PP$prey,PP$predator,type="l",col=col.phase)}
    points(PP$equilibria[,1],PP$equilibria[,2],type="p",col="red")
  }

  if (t.xy & phase){close.screen(all=TRUE)}
}