# Duffing attractor: Complete
library(tseriesChaos)
parms<- c(0.5,1,0.42)  # parameters: c, w, F 
# F=0.325 period 1; F=0.34875 period 2; F=0.42 chaos 
tinit<- 0
tfin<- 650   
step<- 0.1  
times<- seq(tinit,tfin,by=step)
funct<- function(t,integ,parms){
x<-integ[1]
y<-integ[2]
z<-integ[3]
c<- parms[1]
w<- parms[2]
F<- parms[3]
              
dx <- y
	                   
dy <- x*(1-x^2) - c*y + F*cos(z)
 
dz <- w


list(c(dx,dy,dz))
                          }   # end of funct
require(deSolve)
cinit<-c(0,0,0)
xyz<-lsoda(cinit,times,funct,parms)
#xyz # comment if you do not wish the xyz values printed
par(mfrow=c(2,1))
par(mar = c(6.3, 4.8, 1., 3))
par(cex.lab=2,cex.axis=1.6,lwd=1,lty=1)
plot(xyz[,1],xyz[,2],type="l",xlab="t",ylab="x(t)", # x(t) vs t
xlim=c(tinit,tfin),ylim=c(-1.5,1.5))
plot(xyz[,1],xyz[,3],type="l",xlab="t",ylab="y(t)", # y(t) vs t
xlim=c(tinit,tfin),ylim=c(-1,1))
#plot(xyz[,1],xyz[,4],type="l",xlab="t",ylab="z(t)", # z(t) vs t
#xlim=c(tinit,tfin),ylim=c(0,50))
windows()
# phase space portrait
require(scatterplot3d)
scatterplot3d(xyz[,2],xyz[,3],xyz[,4],type="l",#xlim=c(-10,20),
cex.lab=1.4,cex.axis=1.2)

windows()
par(mai=c(1.02,1.,0.82,0.42)+0.1)
plot(xyz[,2],xyz[,3],type="p",pch=20,cex=0.7,xlab="x(t)",ylab="y(t)",
cex.lab=1.5,cex.axis=1.2,lwd=3,lty=1,xlim=c(-1.5,1.5),ylim=c(-1,1))

trans<- 2000 # integration time step considered as transient
# discard initial transient:
x <- window(xyz[,2],trans)
y <- window(xyz[,3],trans)
z <- window(xyz[,4],trans)
t_start<- trans*step - step # new initial time
t_time<-seq(t_start,tfin,by=step) # new time interval
windows()
par(mfrow=c(2,1))
par(mar = c(6.3, 4.8, 1., 3))
par(cex.lab=2,cex.axis=1.6,lwd=1,lty=1)
# x(t), y(t), z(t) after the transient:
plot(t_time,x,type="l",xlab="t",ylab="x(t)",
xlim=c(t_start,tfin),ylim=c(-1.5,1.5))
plot(t_time,y,type="l",xlab="t",ylab="y(t)",
xlim=c(t_start,tfin),ylim=c(-1,1))
#plot(t_time,z,type="l",xlab="t",ylab="z(t)",
#xlim=c(t_start,tfin),ylim=c(0 ,50))

windows()
par(mai=c(1.02,1.,0.82,0.42)+0.1)
plot(x,y,type="p",pch=20,cex=0.7,xlab="x(t)",ylab="y(t)",
cex.lab=1.5,cex.axis=1.2,lwd=2,lty=1,xlim=c(-1.5,1.5),ylim=c(-1,1))



m.max<- 6 # embedding dimensions: from 1 to m_max
d<- 40 # tentative time delay (see below)
tw<- 100 # Theiler window
rt<- 10 # escape factor
eps<- sd(x)/10 # neighbourhood diameter
fn <- false.nearest(x,m.max,d,tw,rt,eps)
fn
windows()
plot(fn)

windows()
lm<- 60                # largest lag
mutual(x,lag.max = lm) # average mutual information to suggest d

# embedding Procedure
m<- 3                # choose embedding dimension
d<- 24               # choose time delay 
xyz <- embedd(x,m,d) # embed the 'observed' series
windows()
par(mai=c(1.02,1.,0.82,0.42)+0.1)
plot(xyz[,2],xyz[,3],type="p",pch=20,cex=0.7,xlab="x(t)",ylab="y(t)",
cex.lab=1.5,cex.axis=1.2,lwd=3,lty=1,xlim=c(-1.5,1.5),ylim=c(-2,2))

windows()
n<- 1000
par(mai=c(1.02,1.,0.82,0.42)+0.1) 
plot(0,0,type="n",xlab="i",ylab="j",xlim=c(1,n),ylim=c(1,n),cex.lab=1.6,cex.axis=1.2,
family="Times",font.lab=3)
distx<- matrix(,n,n)
eps<- 1
xx<- xyz[,1]
yy<- xyz[,2]
zz<- xyz[,3]
for(i in 1:n){                   # construction of the thresholded recurrence plot
for(j in 1:n){
distx[i,j]<- sqrt( (xx[i]-xx[j])^2 + (yy[i]-yy[j])^2 + (zz[i]-zz[j])^2 ) # Euclidean distance
if(distx[i,j]<eps) points(i,j,pch=19,cex=0.01)   
     }
     }

                                 # unthresholded recurrence plot
  
time<- seq(1:n)
windows()
par(mai=c(1.02,1.,0.82,0.42)+0.1,cex.lab=1.6,cex.axis=1.2)  
require(gplots)
distx<- distx/max(distx)
distx<- 1-distx
filled.contour(time,time,distx,col=gray.colors(10,start=1,end=0),nlevels=10,
xlab = "i", ylab = "j", main = "",xlim=c(0,n),ylim=c(0,n),las=0,
key.axes = axis(4, las=1) )

# max exp di Lyapunov (Kantz)   	

S_nu <- lyap_k(x,m=3,d=24,t=50,k=5,ref=4200,s=400,eps=5) 
windows()
par(mai=c(1.02,1.,0.82,0.42)+0.1)
plot(S_nu,xlab = expression(paste(nu)),ylab=expression(paste("S",(nu))),
lwd=2,lty=1,cex=1.5)
lmin<- 60
lmax<- 160
abline(v=lmin,lty=4,lwd=2,col="black")   # to add the vertical lines delimiting the linear region
abline(v=lmax,lty=4,lwd=2,col="black")
lr<- S_nu[lmin:lmax]     # output of lyap_k(.) in [lmin, lmax]
l_lr<- length(lr)-1
lt<- seq(lmin*step,lmax*step,by=step)
lm(lr~lt)                # regression analysis in R
abline(lyap(S_nu,lmin,lmax),lty=2,lwd=2,col="black")


# estimate of D2
C.m <- d2(x,m=6,d=30,t=50,eps.min=0.005,neps=100)      # correlation integral, m = 1,...,6
C.m <- data.frame(unclass(C.m))         # class attribute removed                        
C.3 <- subset(C.m,eps > 0.08 & eps < 0.8, select=c(eps,m3)) # eps in [0.08, 0.8]
lm(log(m3) ~ log(eps), data = C.3)      # D2 estimate with m=3
windows()                               # if other plots are made before
par(mai=c(1.02,1.,0.82,0.42)+0.2)       	
plot(C.m[,1],C.m[,4], type="l",log="xy",main="",  
xlab=expression(paste(epsilon)),ylab=expression(paste(widehat(C),(epsilon))),
lwd=2,lty=1,cex.axis=1.3,cex.lab=2.0, xlim=c(0.004,100),ylim=c(0.000001,1))

# stop here to plot only C(m=3) 
##########################################


# add the lines below to plot C(m=1,...,6)
n.m<- ncol(C.m)

for(i in (n.m-0):2) lines(C.m[,c(1,i)],lwd=2,lty=1)
abline(v=0.08,lty=4,lwd=2,col="black")  # non messe
abline(v=0.8,lty=4,lwd=2,col="black")