library(deSolve)

################################################################################
SLODE <- function(t,state,parameters){
    with(as.list(c(state,parameters)),{
        dx <- A[1,1]*x+A[1,2]*y
        dy <- A[2,1]*x+A[2,2]*y
        list(c(dx,dy))})
}

################################################################################
# Plot solution

PlotSol <- function(outLSode,Tlm=NULL,Xlm=NULL){
    if (length(Tlm) == 0){
        Tlm <- c(min(outLSode[,"time"]),max(outLSode[,"time"]))}
    if (length(Xlm) == 0){
        xmax <- max(abs(outLSode[,"x"]),abs(outLSode[,"y"]));
        Xlm <- c(-xmax,xmax)}

    plot(Tlm,Xlm,type="n",bty="l",xlab="t",ylab="x, y",main="Solution")
    lines(Tlm,c(0,0))
    lines(outLSode[,"time"],outLSode[,"x"],col="red",lwd=2)
    lines(outLSode[,"time"],outLSode[,"y"],col="blue",lwd=2)
}

################################################################################
# Plot phase portrait

PlotPhP <- function(A,outLSode,xlm=NULL,ylm=NULL){
    if (length(xlm) == 0){xe <- max(abs(outLSode[,"x"]))
        xlm <- c(-xe,xe)}
    if (length(ylm) == 0){ye <- max(abs(outLSode[,"y"]))
        ylm <- c(-ye,ye)}

    plot(outLSode[,"x"],outLSode[,"y"],type="l",lwd="2",xlab="x",ylab="y",
         xlim=xlm,ylim=ylm,main="Phase portrait")

    colEq="gray"
    detA <- det(A); trA <- A[1,1]+A[2,2]
    if (detA < 0){colEq <- "cyan"}
    else {if (detA > 0){colEq <- "brown"
                        if (trA < 0){colEq <- "black"}
                        if (trA > 0){colEq <- "green"}
          }}
    points(0,0,pch=19,col=colEq)              
}

################################################################################
# Visualize analysis

PlotAn <- function(A,X=NULL,Y=NULL,dens=20){
    detA <- det(A)
    trA <- A[1,1]+A[2,2]
    if (length(X) == 0){X <- 1.5*max(1,abs(trA))}
    if (length(Y) == 0){Y <- 1.5*abs(detA)}


    x <- seq(-X,X,length.out=201)
    x1 <- seq(-0.8*X,0,length.out=201)
    x2 <- seq(0,0.8*X,length.out=201)
    y <- 0.25*x^2
    y1 <- 0.25*x1^2
    y2 <- 0.25*x2^2
    Xm <- max(x2[which(x2<6/7*X)])
    Ym <- max(y2[which(x1<6/7*X)])
    Yp <- 0.9*Y

    plot(X*6/7*c(-1,1),Y*c(-1,1),type="n",bty="l",xlab="trA",ylab="detA")
    polygon(Xm*c(-1,1,1,-1),Yp*c(-1,-1,0,0),angle=90,density=dens,col="cyan",
            border=NA)
    if (Yp > Ym){
        polygon(c(x1,-Xm,-Xm),c(y1,0,max(y1)),angle=-45,density=dens,border=NA,
                col="green")
        polygon(c(x2,Xm,0),c(y2,0,0),angle=45,density=dens,border=NA,
                col="green")
        polygon(c(x1,0,-Xm),c(y1,Yp,Yp),angle=-45,density=dens,border=NA,
                col="gold")
        polygon(c(x2,Xm,0),c(y2,Yp,Yp),angle=45,density=dens,border=NA,
                col="gold")
        Yg <- Ym
    }
    else {
        polygon(c(x1[which(y1<=Yp)],-Xm,-Xm),c(y1[which(y1<=Yp)],0,Yp),
                angle=-45,density=dens,border=NA,col="green")
        polygon(c(x2[which(y2<=Yp)],Xm,Xm),c(y2[which(y2<=Yp)],Yp,0),
                angle=45,density=dens,border=NA,col="green")
        polygon(c(x1[which(y1<=Yp)],0),c(y1[which(y1<=Yp)],Yp),angle=-45,
                density=dens,border=NA,col="gold")
        polygon(c(x2[which(y2<=Yp)],0),c(y2[which(y2<=Yp)],Yp),angle=45,
                density=dens,border=NA,col="gold")
        Yg <- Yp
    }
    
    lines(X*c(-1,1),c(0,0))
    lines(c(0,0),Y*c(-1,1))
    lines(x,y)
    lines(x1[which(y1<Yg)],y1[which(y1<Yg)],lwd=3,col="green")
    lines(x2[which(y2<Yg)],y2[which(y2<Yg)],lwd=3,lty=2,col="green")
    lines(c(0,0),Yp*c(0,1),lwd=3,col="brown")
    lines(Xm*c(-1,1),c(0,0),col="white")
    points(0,0,pch=19,col="white")

    points(trA,detA,pch=23,col="red",bg="red")
}
################################################################################
# Plot results

PlotRes <- function(A,outLSode,
                    Tlm=NULL,Xlm=NULL,xlm=NULL,ylm=NULL,
                    X=NULL,Y=NULL){
    split.screen(c(2,1))
    screen(1); PlotSol(outLSode,Tlm,Xlm)
    screen(2); split.screen(c(1,2)); screen(3); PlotPhP(A,outLSode,xlm,ylm)
                                     screen(4); PlotAn(A,X,Y)
    close.screen(all=TRUE)
}

################################################################################
################################################################################

S2Lode <- function(A,tmax=50,InitVal=c(x=0,y=1),
                   Tlm=NULL,Xlm=NULL,   # limits for graph
                   xlm=NULL,ylm=NULL,   # limits for phase portrait
                   X=NULL,Y=NULL        # limits for parameter plot
                   ){
    out <- ode(InitVal,seq(0,tmax,by=tmax/500),func=SLODE,parms=A)
    PlotRes(A,out,Tlm,Xlm,xlm,ylm,X,Y)

    detA <- det(A); trA <- A[1,1]+A[2,2]
    write("\n** System of 2 linear ODE's **\n","")
    write("Matrix of the system:",""); print(A)
    if (detA == 0){write("\nEquilibrium: not unique","")}
    else {if (detA < 0){write("\nEquilibrium: saddle point","")}
          else {if (trA == 0){write("\nEquilibrium: center","")}
                else {if (trA^2 >4*detA){tp <- "node"}
                      else{tp <- "focus"}}}}
    if (detA>0 & trA>0){write(paste("\nEquilibrium: source,",tp,sep=" "),"")}
    if (detA>0 & trA<0){write(paste("\nEquilibrium: sink,",tp,sep=" "),"")}
}


################################################################################
################################################################################


tmax <- 15
A <- matrix(c(-0.1,-0.5,0.5,0.2),c(2,2))

S2Lode(A)