# Ideal pendulum: unlimited oscillations (pend_id_un_o_1)
# trajectory of libration motion are plotted;
# nullclines and velocity field are plotted;
# equilibrium points are classified;
# x and y variables are plotted as a function of the time.

library(phaseR)
pend<-function(t,y,parameters){
x <- y[1]
y <- y[2]
a<-parameters[1]
b<-parameters[2]
c<-parameters[3]
d<-parameters[4]
dy<-numeric(2)
dy[1]<- a*x+b*y
dy[2]<- c*sin(x)+d*y
list(dy)
}

# LIBRATION

omega2<- 2.5
parms<-c(0,1,-omega2,0) #parameter values
t.end<- 100
delta<-0.05
x0<-   0
y0<-   2
xy0<-cbind(x0,y0) #initial conditions
plot(0,0,type="n" ,xlim=c(-2.5*pi,2.5*pi),ylim=c(-10,10),
cex.axis=1.2,cex.lab=1.5,font.lab=3,xlab="x",ylab="y")
trajectory<-trajectory(pend,y0=xy0,t.end=t.end,t.step=delta,parameters=parms,
lwd=2,lty=1,colour="black",add=TRUE)

arrows<-flowField(pend,lwd=2,x.lim=c(-2.5*pi,2.5*pi),y.lim=c(-10,10),
parameters=parms,points=19,arrow.head=0.05,frac=1,colour="gray60",add=TRUE)
null<-nullclines(pend,parameters=parms,lwd=2,points=151,
x.lim=c(-2.5*pi,2.5*pi),y.lim=c(-10,10),colour=c("red","blue") )

stability.1 <- stability(pend, y.star=c(0,0),parameters=parms)
stability.2 <- stability(pend, y.star=c(pi/2,0),parameters=parms)

windows()
numericalSolution<- numericalSolution(pend,y0=c(1,1),t.end=20,type="one",t.step=delta,
parameters=parms,colour=c("green","orange"),lty=c(1,2),
ylab="x, y",ylim=c(-2.5,2.5),lwd=2,cex.axis=1.2,cex.lab=1.5)
legend("bottomright",col=c("green","orange"),legend=c("x","y"),lty=c(1,2))