BMat=rbind(c(0.333,0,0,0.333,-0.333,0),c(0.167,-0.289,0.289,0.167,-0.0830,0.144),c(0.167,0.289,-0.289,0.167,0.083,0.144),c(0.333,0,0,0.333,0.333,0))

# Initial conditions:

x=0

y=0

plot(0,0,xlim=c(-0.5,0.5),ylim=c(0,1),col="white")

COLOR=c("green","red","blue","yellow")

for(j in 1:10000)

{

i=sample(1:4,1) # ,prob=c(0.25,0.25,0.25,0.25)

x3=x

x=BMat[i,1]*x+BMat[i,2]*y+BMat[i,5]

y=BMat[i,3]*x3+BMat[i,4]*y+BMat[i,6]

points(x,y,pch=".",cex=1, col=COLOR[i])

}

## Thursday, March 29, 2007

## Wednesday, March 21, 2007

### A flower?

BMat=rbind(c(2/3,1/3,0,2/3,0,2/3),c(-1/3,1/3,0,1/3,-1/3,2/3),c(2/3,1/3,1/3,1/3,1/3,2/3),c(1/3,0,2/3,-2/3,1/3,0))

# Initial conditions:

x=0

y=0

X11()

plot(0,0,xlim=c(-1,3),ylim=c(-1,3),col="white")

COLOR=c("green","red","green","red")

for(j in 1:30000)

{

i=sample(1:4,1) # ,prob=c(0.25,0.25,0.25,0.25)

x=BMat[i,1]*x+BMat[i,2]*y+BMat[i,5]

y=BMat[i,3]*x+BMat[i,4]*y+BMat[i,6]

points(x,y,pch=".",cex=1, col=COLOR[i])

}

# Initial conditions:

x=0

y=0

X11()

plot(0,0,xlim=c(-1,3),ylim=c(-1,3),col="white")

COLOR=c("green","red","green","red")

for(j in 1:30000)

{

i=sample(1:4,1) # ,prob=c(0.25,0.25,0.25,0.25)

x=BMat[i,1]*x+BMat[i,2]*y+BMat[i,5]

y=BMat[i,3]*x+BMat[i,4]*y+BMat[i,6]

points(x,y,pch=".",cex=1, col=COLOR[i])

}

### Mandelbrot set

Here you can find the code for the Mandelbrot set:

http://tolstoy.newcastle.edu.au/R/help/03b/3377.html

http://tolstoy.newcastle.edu.au/R/help/03b/3377.html

### Barnsley fern

BMat=rbind(c(0.849,0.037,-0.037,0.849,0.075,0.183),c(0.197,-0.226,0.226,0.197,0.400,0.049),c(-0.150,0.283,0.260,0.237,0.575,-0.084),c(0,0,0,0.160,0.500,0))

# Initial conditions:

x=0

y=0

X11()

plot(0,0,xlim=c(0,1),ylim=c(0,1),col="white")

COLOR=c("green","red","blue","yellow")

for(j in 1:10000)

{

i=sample(1:4,1) # ,prob=c(0.25,0.25,0.25,0.25)

x3=x

x=BMat[i,1]*x+BMat[i,2]*y+BMat[i,5]

y=BMat[i,3]*x3+BMat[i,4]*y+BMat[i,6]

points(x,y,pch=".",cex=1, col=COLOR[i])

}

#--------------------------------------#

# Alternative code by Stefano Guazzetti

#--------------------------------------#

# The probabilities used here give more dense fern.

niter<-150000

kind<-sample(1:4, niter, repl=T, prob=c(.01, .07, .07, .85))

x<-numeric(niter+1)

y<-numeric(niter+1)

x[1]<-0

y[1]<-0

for (i in 1:niter) {

x[i+1]<- 0*(kind[i]==1)+(0.2*x[i]-0.26*y[i])*(kind[i]==2)+

(-0.15*x[i]+0.28*y[i])*(kind[i]==3) +

(0.85*x[i]+0.04*y[i])*(kind[i]==4)

y[i+1]<- 0.16*y[i]*(kind[i]==1)+(0.23*x[i]+0.22*y[i]+1.6)*(kind[i]==2)+

(0.26*x[i]+0.24*y[i]+.44)*(kind[i]==3) +

(-0.04*x[i]+0.85*y[i]+1.6)*(kind[i]==4)

}

par(mar=c(0.1,0.1,0.1,0.1))

plot(x, y, pch=17, cex=.3, col="darkgreen", axes=F, ann=F)

# Initial conditions:

x=0

y=0

X11()

plot(0,0,xlim=c(0,1),ylim=c(0,1),col="white")

COLOR=c("green","red","blue","yellow")

for(j in 1:10000)

{

i=sample(1:4,1) # ,prob=c(0.25,0.25,0.25,0.25)

x3=x

x=BMat[i,1]*x+BMat[i,2]*y+BMat[i,5]

y=BMat[i,3]*x3+BMat[i,4]*y+BMat[i,6]

points(x,y,pch=".",cex=1, col=COLOR[i])

}

#--------------------------------------#

# Alternative code by Stefano Guazzetti

#--------------------------------------#

# The probabilities used here give more dense fern.

niter<-150000

kind<-sample(1:4, niter, repl=T, prob=c(.01, .07, .07, .85))

x<-numeric(niter+1)

y<-numeric(niter+1)

x[1]<-0

y[1]<-0

for (i in 1:niter) {

x[i+1]<- 0*(kind[i]==1)+(0.2*x[i]-0.26*y[i])*(kind[i]==2)+

(-0.15*x[i]+0.28*y[i])*(kind[i]==3) +

(0.85*x[i]+0.04*y[i])*(kind[i]==4)

y[i+1]<- 0.16*y[i]*(kind[i]==1)+(0.23*x[i]+0.22*y[i]+1.6)*(kind[i]==2)+

(0.26*x[i]+0.24*y[i]+.44)*(kind[i]==3) +

(-0.04*x[i]+0.85*y[i]+1.6)*(kind[i]==4)

}

par(mar=c(0.1,0.1,0.1,0.1))

plot(x, y, pch=17, cex=.3, col="darkgreen", axes=F, ann=F)

### Sierpinsky gasked

# Sierpinsky gasked, Chaos game -algorithm:

X11(width=6,height=6)

plot(0,0,xlim=c(0,1),ylim=c(0,1), col="white")

A=c(0,0)

B=c(1,0)

C=c(0.5,0.8)

D=sample(0:1000,2)/1000

for(i in 1:10000)

{

Corner=sample(1:3,1)

if(Corner==1){CornerS=A}

if(Corner==2){CornerS=B}

if(Corner==3){CornerS=C}

D=(D+CornerS)/2

points(D[1],D[2],pch=".",cex=1)

}

X11(width=6,height=6)

plot(0,0,xlim=c(0,1),ylim=c(0,1), col="white")

A=c(0,0)

B=c(1,0)

C=c(0.5,0.8)

D=sample(0:1000,2)/1000

for(i in 1:10000)

{

Corner=sample(1:3,1)

if(Corner==1){CornerS=A}

if(Corner==2){CornerS=B}

if(Corner==3){CornerS=C}

D=(D+CornerS)/2

points(D[1],D[2],pch=".",cex=1)

}

Subscribe to:
Posts (Atom)