#### x.new<-A*x.old*(1-x.old) # Logistic Function

# Settings
iters<-200  # Number of iterations
step<-.1 # Increase A by this step size
burn.in.steps<-150 # Throw Away first n steps when plotting A vs x.new
sleep.time<-.0000001 # Control speed (in seconds)

x.old.init<-seq(0.1,.9,by=.9) # Initial x values
A.max<-4 # Max A value
A<- 1 # Initial value for A parameter



#Initialize Variables and Chart
x.old<-x.old.init[1] # Set initial x value to first value determined above
out=cbind(A,x.old.init[1],x.old) # Initialize output matrix

dev.new(); plot(0,0,type="n", ylim=c(0,1), xlim=c(1,iters))

#Run function for increasing A values until equal to A max
while(A<=A.max){
 
  # Set initial x value to next one in vector determined above
  for(j in x.old.init){
    x.old<-j
   
    # Run for number of iterations set by iters
    for(i in 1:iters){
     
      x.new<-A*x.old*(1-x.old) # Logistic Function (update x.new)
     
      out<-rbind(out,cbind(A,j,x.new)) # Update output
     
      # Plot current run
      plot(c(j,out[which(out[,1]==A & out[,2]==j),3]), type="n", lwd=2,
           col="White",
           ylim=c(0,1), xlim=c(0,iters), pch=16,
           xlab="Iteration", ylab="X.new",
           main=c("X.new = A*X.old*( 1 - X.old )",
                  paste("A=",A,  "   X.initial=",j),
                  paste("X.new=", round(x.new,3))))
      for(k in x.old.init){
        color<-rainbow(length(x.old.init))[which(x.old.init==k)]
        points(c(j,out[which(out[,1]==A & out[,2]==k),3]), type="b", lwd=2,
               col=color,
               ylim=c(0,1), xlim=c(0,iters), pch=16,
               xlab="Iteration", ylab="X.new",
               main=c("X.new = A*X.old*( 1 - X.old )",
                      paste("A=",A,  "   X.initial=",j),
                      paste("X.new=", round(x.new,3))))
      }
      ##
     
     
      x.old<-x.new # x.old becomes x.new
    }
    Sys.sleep(sleep.time)
   
  }
  A<-round(A+step,2) # update A parameter
 
}

#Plot "A" parameter vs Results
dev.new()
plot(out[,1],out[,3], type="n", ylab="X.new", xlab="A",
     main="X.new = X.old*( 1 - X.old )"
)
for(i in 1:length(unique(out[,1]))){
  points(out[which(out[,1]==unique(out[,1])[i])[-(1:burn.in.steps)],1],
         out[which(out[,1]==unique(out[,1])[i])[-(1:burn.in.steps)],3])
}

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


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

#####Settings
add.noise=F  #Gaussian +/- noise.sd
noise.sd=.01
live.plot=T
##

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

n=5000   #Iterations
##


####Initial Conditions
Xa=0.01; Ya=0.01; Za=0.01    #Initial Condition Blue
Xb=0.01+.0001; Yb=0.01; Zb=0.01    #Initial Condition Red (Slight Difference)
Xc=20; Yc=20; Zc=.01    #Initial Condition Black (Large Difference)
##


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

par3d(font=2, family="serif",
      bg3d(color=c("darkslategray3","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==T){
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)

if(live.plot==T){
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="Black", alpha=.7, add=T)

points3d(XYZa[i,1],XYZa[i,2],XYZa[i,3], col="Blue", alpha=1, size=10, add=T)
points3d(XYZb[i,1],XYZb[i,2],XYZb[i,3], col="Red", alpha=1, size=10, add=T)
points3d(XYZc[i,1],XYZc[i,2],XYZc[i,3], col="Black", alpha=1, size=10, add=T)

rgl.pop()
rgl.pop()
rgl.pop()
}

}
##

if(!live.plot==T){
points3d(XYZa[,1],XYZa[,2],XYZa[,3], col="Blue", alpha=.5, add=T)
points3d(XYZb[,1],XYZb[,2],XYZb[,3], col="Red", alpha=.5, add=T)
points3d(XYZc[,1],XYZc[,2],XYZc[,3], col="Black", alpha=.5, add=T)
}