############################
##LORENZ strange attractor##
############################
#Modified From:
#http://fractalswithr.blogspot.com/2007/04/lorenz-attractor.html


# install.packages("rgl") # Install if needed.
library(rgl)

### Settings
# Gaussian +/- noise.sd
add.noise = T
noise.sd  =.1


### Parameters
a  = 10
r  = 28
b  = 8/3
dt = 0.01
n  = 5000



### Initial Conditions
# Initial Condition Blue
Xa = 0.01
Ya = 0.01
Za = 0.01

# Initial Condition Red (Slight Difference)
Xb = 0.01 + .0001
Yb = 0.01
Zb = 0.01

# Initial Condition White (Large Difference)
Xc = 20
Yc = 20
Zc = .01




### Initialize
XYZa = array(0, dim = c(n, 3))
XYZb = array(0, dim = c(n, 3))
XYZc = array(0, dim = c(n, 3))

par3d(windowRect = c(20, 30, 800, 800),
      bg3d(color   = c("Black", "Black"),
           fogtype = "exp2",
       sphere  = TRUE,
       back    = "fill"))


### Run
for(i in 1:n){
  X1a = Xa
  Y1a = Ya
  Z1a = Za

  Xa = X1a + (-a*X1a   + a*Y1a)*dt
  Ya = Y1a + (-X1a*Z1a + r*X1a - Y1a)*dt
  Za = Z1a + (X1a*Y1a  - b*Z1a)*dt


  X1b = Xb
  Y1b = Yb
  Z1b = Zb

  Xb = X1b + (-a*X1b   + a*Y1b)*dt
  Yb = Y1b + (-X1b*Z1b + r*X1b - Y1b)*dt
  Zb = Z1b + (X1b*Y1b  - b*Z1b)*dt


  X1c = Xc
  Y1c = Yc
  Z1c = Zc

  Xc = X1c + (-a*X1c   + a*Y1c)*dt
  Yc = Y1c + (-X1c*Z1c + r*X1c - Y1c)*dt
  Zc = Z1c + (X1c*Y1c  - b*Z1c)*dt


  if(add.noise){
    Xa = rnorm(1, Xa, noise.sd)
    Xb = rnorm(1, Xb, noise.sd)
    Xc = rnorm(1, Xc, noise.sd)

    Ya = rnorm(1, Ya, noise.sd)
    Yb = rnorm(1, Yb, noise.sd)
    Yc = rnorm(1, Yc, noise.sd)

    Za = rnorm(1, Za, noise.sd)
    Zb = rnorm(1, Zb, noise.sd)
    Zc = rnorm(1, Zc, noise.sd)
  }

  XYZa[i, ] = c(Xa, Ya, Za)
  XYZb[i, ] = c(Xb, Yb, Zb)
  XYZc[i, ] = c(Xc, Yc, Zc)


  points3d(XYZa[i, 1], XYZa[i, 2], XYZa[i, 3], col = "Blue",  alpha = .7, add = T)
  points3d(XYZb[i, 1], XYZb[i, 2], XYZb[i, 3], col = "Red",   alpha = .7, add = T)
  points3d(XYZc[i, 1], XYZc[i, 2], XYZc[i, 3], col = "White", alpha = .7, add = T)

}