**by Allan Roberts**

It’s December. Perhaps it’s dark, cold, rainy, and not very snowy at the Bamfield Marine Sciences Centre. But that can’t stop you from making paper snowflakes. In fact, you don’t need paper. That’s right, you can make snowflakes with R. (See Figure 1.)

When you run the script provided below, a six-sided flake will appear in the graphics window, with fold lines for reference. Click on the graphics window. You will see a “+”. Click 3 times to cut a triangle out of the snowflake. You have 3 triangles to design your snowflake. (See Figure 2.) After you cut out the 3rd triangle, you get a snowflake! (At any time you can end your snowflake making session by closing the graphics window.)

How does it work? Matrix algebra. Matrices also have biological applications. See, for example, Caswell (2001).

For more about snow & ice see the following: sea ice data from Cambridge Bay; NOAA arctic web cam videos; the ecological and cultural importance of sea ice at Arctic Eider Society.

**Reading
**

Caswell, Hall. 2001. Matrix population models: construction, analysis and interpretation, 2nd ed. Sunderland, Mass.: Sinauer Associates.

**R Script**

#Script written by Allan Roberts, 2012.

snowflake <- function(k=3){

Rotate <- function(a){

#returns a matrix for a rotation of a.

A <- matrix(0,2,2);

A[1,1] <- cos(a);

A[1,2] <- -sin(a);

A[2,1] <- sin(a);

A[2,2] <- cos(a);

return(A);

}

Reflect <- function(a){

#a reflection across the line determined by angle a.

A <- matrix(0,2,2);

A[1,1] <- 1;

A[2,2] <- -1;

B <- Rotate(a) %*% A %*% Rotate(-a);

return(B);

}

locator.with.point <- function(n){

X <- numeric(n);

Y <- numeric(n);

for (i in 1:n){P<-locator(1);X[i]<-P$x;Y[i]<-P$y; points(X[i],Y[i],pch=19,col=4)}

return(list(X,Y))

}

par( bg = 1 )

X1 <- numeric(6);

Y1 <- numeric(6);

for (i in 1:6) {X1[i] <- sin(i*2*pi/6)}

for (i in 1:6) {Y1[i] <- cos(i*2*pi/6)}

plot.new();

par(usr=c(-2,2,-2,2), bg=1);

polygon(X1,Y1,col=rgb(1,1,1));

abline(0,tan(pi/6),lty=3);

abline(0,-tan(pi/6),lty=3);

lines(c(0,0),c(-1,1),lty=3);

C <- list(matrix(0,2,3),matrix(0,2,3),matrix(0,2,3));

for (i in 1:k){

X <- locator.with.point(3);

X <- rbind(X[[1]],X[[2]]);

C[[i]] <- X;

polygon(t(C[[i]]),col=rgb(1,1,1,0),border=1,lty=2,cex=0.2);

}

#Draw the Snowflake

plot.new();

par(usr=c(-2,2,-2,2), bg=1);

polygon(X1,Y1,col=rgb(1,1,1));

for (j in 1:k){

for (i in 1:6){ R <- Rotate(i*pi/3) %*% Reflect(pi/3) %*% C[[j]]; polygon(t(R), col=1)}

for (i in 1:6){R <- Rotate(i*pi/3) %*% C[[j]]; polygon(t(R),col=1);}

}

}

snowflake( )

Pingback: Let it snow! « Probability and statistics blog