#
#
################ STRANGE ATTRACTORS ################
#

##############################################
#
# The Henon attractor
# 

gc()

# First, establish some original values for x & y.

x = 0.4
y = 0.2

# To have even smaller 'dots' in the plot, insert: pch = "." into both plot 
# statements (directly below & inside the function).
# For example: plot(x, y, xlim = c(-1.75, 1.75), ylim = c(-.6, .6), pch = ".", col = "blue")

plot(x, y, xlim = c(-1.75, 1.75), ylim = c(-.6, .6), col = "blue")

for (i in 1:10000){
 x.new = y + 1 - (1.4 * (x^2))
 y.new = x * 0.3
 par(new = T)
 plot(x.new, y.new, xlim = c(-1.75, 1.75), ylim = c(-.6, .6), col = "blue")
 y = y.new
 x = x.new
 }


###############################################
#
# Lorenz attractor 
#

library(scatterplot3d)

gc()

a = 10
r = 28
b = 8/3
dt = 0.01
X2 = 0.01
Y2 = 0.01
Z2 = 0.01
x <- as.vector(0)
y <- as.vector(0)
z <- as.vector(0)
n = 15000

for(i in 1:n){
 X1 = X2
 Y1 = Y2
 Z1 = Z2
 X2 = X1 + (-a * X1 + a * Y1) * dt
 Y2 = Y1 + (-X1 * Z1 + r * X1 - Y1) * dt
 Z2 = Z1 + (X1 * Y1 - b * Z1) * dt
 x[i] = X2
 y[i] = Y2
 z[i] = Z2
 scatterplot3d(x, y, z, highlight.3d = TRUE, xlab = "x", ylab = "y", zlab = "z", type = "p", 
   col.axis = "blue", angle = 55, col.grid = "lightblue", scale.y = 0.7, pch = ".", 
   xlim = c(-30, 30), ylim = c(-30, 30), zlim = c(0, 60))
}




# End Jan. 2011