# Ideal pendulum: unlimited oscillations (pend_id_un_o)
library(phaseR)
#Set parameters for integration
t.end<- 20  #integrate from t=0 to t.end
delta<-0.05

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)
}

# ROTATION

omega2<- 1.
parms<-c(0,1,-omega2,0) #parameter values
x0<-c(-4*pi,4*pi)
y0<-c(4,-4)
xy0<-cbind(x0,y0) #initial conditions

plot(0,0,type="n" ,xlim=c(-4.5*pi,4.5*pi),ylim=c(-5,5),
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=3,lty=2,colour=rep("blue",2),add=TRUE)

omega2<- 1.01
parms<-c(0,1,-omega2,0) #parameter values
delta<-0.05
x0<-c(-4*pi,4*pi)
y0<-c(2.2,-2.2)
xy0<-cbind(x0,y0) #initial conditions
trajectory<-trajectory(pend,y0=xy0,t.end=t.end,t.step=delta,parameters=parms,
lwd=3,lty=2,colour=rep("blue",2),add=TRUE)


# LIBRATION

omega2<- 2.5
parms<-c(0,1,-omega2,0) #parameter values
t.end<- 100
delta<-0.05
x0<-c(1-4*pi,1-2*pi,1,1+2*pi,1+4*pi)
y0<-c(1,1,1,1,1)
xy0<-cbind(x0,y0) #initial conditions
trajectory<-trajectory(pend,y0=xy0,t.end=t.end,t.step=delta,parameters=parms,
lwd=2,lty=1,colour=rep("black",5),add=TRUE)
omega2<- 0.5
parms<-c(0,1,-omega2,0) #parameter values
t.end<- 20
delta<-0.05
x0<-c(1-4*pi,1-2*pi,1,1+2*pi,1+4*pi)
y0<-c(1,1,1,1,1)
xy0<-cbind(x0,y0) #initial conditions
trajectory<-trajectory(pend,y0=xy0,t.end=t.end,t.step=delta,parameters=parms,
lwd=2,lty=1,colour=rep("black",5),add=TRUE)

# SEPARATRIX

omega2<- 1
parms<-c(0,1,-omega2,0) #parameter values
t.end<- 100
delta<-0.05
x0<-c(-5*pi,-3*pi,-pi,pi,3*pi,5*pi)
y0<-c(0.001,0.001,0.001,-0.001,-0.001,-0.001)
xy0<-cbind(x0,y0) #initial conditions
trajectory<-trajectory(pend,y0=xy0,t.end=t.end,t.step=delta,parameters=parms,
lwd=2,lty=3,colour=rep("red",2),add=TRUE)

seq_pie<-seq(-4*pi,4*pi,by=2*pi)             # even
seq_pie
points(seq_pie,rep(0,5),col="black",pch=19,cex=0.8)
seq_pio<-seq(-5*pi,5*pi,by=2*pi)             # odd
seq_pio
segments(seq_pio,-5,seq_pio,5,col="green3",lwd=2,lty=4)