\(\alpha\) which determines the decay
\(\beta\) determines the frequency
and \(p\) is the phase
We generate a matrix which we've called h. Each element corresponds to an x and a y. We can convert that h matrix to an image and then cycle through values of phase to animate it. We will then have an image such as the one above.
Here is the R code that generates the h matrix:
# generate an image of a ripple
rippleMat <- function( numPx, decay, phase, freq){
x <- seq( -1, 1, 2/ (numPx-1))
y <- x
# We have a matrix r[i,j] = sqrt( x[i] * x[i] + y[j] * y[j])
r <- sqrt( outer( x*x, y*y, "+"))
h <- exp( - decay * r ) * cos( phase + freq * r )
return (h)
}
(The full code used to generate that animated gif is show below)
Alternatively, we could also plot a spiral. The basic equation is: \[r = exp( \phi \alpha + p ) \] Where r is the distance from the origin, \(\phi\) is the angle to the x-axis, p is the phase and \(\alpha\) determines how tightly coilded it is.
We could of course plot such a path using a line chart in R. But it might be a bit more interesting to color it in. Again we could iterate through values of the phase p to generate an animated GIF file:
We could also build an RGB color triangle. Each corner is set to one of { red, green, blue }. That color then fades across the triangle until it is completely absent on the opposite edge.
Just in case you're interested, here's the full source code for the triangle:
# Author: Philip Kinlen, Aug 2014
library(grid) # required for grid.raster(.)
# numPix should be an integer,
# and gapX < (numPix / 2)
# where gapX is the minimum number of horizontal pixels from the triangle to the edge of the image
RGBTriangle <- function(numPix = 100, gapX = 10, doSmoothing=FALSE) {
# the verticle gap between the triangle and edge of image
gapY <- round(numPix * ( 0.5 - sqrt(3)/4) + gapX * sqrt(3)/2)
xArr <- 1:numPix
yArr <- numPix:1 # the function call grid.raster(..) below will put the element with
# highest row number at the bottom.
# For our array, the last shall be at the bottom.
# The x's a constant in the matrix columns, y's are constant along the matrix rows.
xMat <- matrix(rep(xArr, numPix), numPix, numPix, F)
yMat <- matrix(rep(yArr, numPix), numPix, numPix, T)
m1 <- sqrt(3) # slope
c1 <- gapY - m1 * gapX # intercept
m2 <- -sqrt(3) # slope
c2 <- gapY - m2 * (numPix - gapX) # intercept
height <- numPix - 2 * gapY # height of triangle in pixels
red <- matrix(mapply(getRed, xMat, yMat, m1, c1, m2, c2, gapY, height), numPix, numPix, T)
green <- matrix(mapply(getGreen, xMat, yMat, m1, c1, m2, c2, gapY, height), numPix, numPix, T)
blue <- matrix(mapply(getBlue, xMat, yMat, m1, c1, m2, c2, gapY, height), numPix, numPix, T)
col <- rgb(red, green, blue)
dim(col) <- dim(red)
grid.raster(col, interpolate=doSmoothing)
# One possible source of minor confusion is that when we write a cartesian co-ordinate (x,y)
# the x comes first.
# but in a matrix, we specify the row first, which corresponds to y
}
#################################################################
getRed <- function( x, y, m1, c1, m2, c2, gapY, height){
if ( isInsideTriangle( x, y, m1, c1, m2, c2, gapY)){
res <- (y - gapY) / height # A rounding error may cause res to be outside the range [0,1]
return (max(0, min(1, res))) # If left untreated that could cause an error in grid.raster(.)
} else
return (0)
}
#################################################################
getBlue <- function( x, y, m1, c1, m2, c2, gapY, height){
if ( isInsideTriangle( x, y, m1, c1, m2, c2, gapY)){
res <- sqrt(3) * ((y - c2) / m2 - x) / ( 2 * height)
return (max(0, min(1, res)))
} else
return (0)
}
#################################################################
getGreen <- function( x, y, m1, c1, m2, c2, gapY, height){
if ( isInsideTriangle( x, y, m1, c1, m2, c2, gapY)){
res <- 0.5 * (m1 * x + c1 - y) / height
return (max(0, min(1, res)))
} else
return (0)
}
#################################################################
isInsideTriangle <- function(x,y, m1, c1, m2 ,c2, gapY){
if ( ( y > gapY)
& ( y < ( m1 * x + c1 ) )
& ( y < ( m2 * x + c2 ) ) )
return(T)
else
return (F)
}
# Author: Philip Kinlen, Aug 2014
###################################################################
library( caTools) # library required for write.gif(..)
###################################################################
animateRipples <- function(){
filePath <- paste(getwd(), "/animatedRipples.gif", sep="")
numImages <- 25
numSecForCycle <- 1
numPx <- getNumPixels()
phases <- -2 * pi * ( 0:(numImages-1)) / numImages
# The following line does the main work generating the data.
# The call to sapply(..) does the work that would have been done by a for-loop in
# a language such as a Java.
zeds <- sapply(phases, rippleMat, simplify="array")
delayCS <- 100 * numSecForCycle / numImages # in hundreths of a second
write.gif( zeds, filePath, col= getCustomColors(), delay = delayCS )
}
###################################################################
animateSpiral <- function(){
# for MS Windows might need to set the separator argument to a back-slash
filePath <- paste(getwd(), "animatedSpiral.gif", sep="/")
numImagesPerSec <- 20
numSecForCycle <- 5
numImages <- numImagesPerSec * numSecForCycle
numPx <- getNumPixels()
phases <- 2 * pi * ( 0:(numImages-1)) / (numImages * getSpiralBeta())
# The following line does the main work generating the data.
# The call to sapply(..) does the work that would have been done by a for-loop in
# a language such as a Java.
zeds <- sapply( phases, spiralMat, simplify="array")
delayCS <- 100 * numSecForCycle / numImages # in hundreths of a second
write.gif( zeds, filePath, col= getCustomColors(), delay = delayCS )
}
###################################################################
getSpiralBeta <- function(){
return ( 0.2)
}
###################################################################
getNumPixels <- function(){
return ( 200)
}
###################################################################
getCustomColors <- function(){
maxCol <- 255
rising <- 0:maxCol
falling <- maxCol - rising
myRed <- c(falling, rising) * 0.4 + 0.2
myGreen <- c(falling, rising) * 0.4 + 0.2
myBlue <- c(falling, rising) * 0.8 + 0.2
return (rgb( myRed, myGreen, myBlue, maxCol, maxColorValue=maxCol ))
}
###################################################################
# generate a matrix representation of a ripple
rippleMat <- function( phase, numPx = getNumPixels(), decay = 2.5, freq = 12 * pi){
x <- seq( -1, 1, 2/ (numPx-1))
y <- x
# We have a matrix r[i,j] = sqrt( x[i] * x[i] + y[j] * y[j])
r <- sqrt( outer( x*x, y*y, "+"))
h <- exp( - decay * r ) * cos( phase + freq * r )
return (h)
}
###################################################################
spiralMat <- function(phase, beta = getSpiralBeta() , numPx= getNumPixels(), generateImage= FALSE){
alpha <- 0.1 # the main spiral follows: r = exp( alpha * phi)
x <- seq(-1, 1, 2 / (numPx - 1))
y <- x
matR <- sqrt(outer(x*x, y*y, "+"))
matX <- matrix(rep(x, numPx), numPx, numPx, F)
matY <- matrix(rep(y, numPx), numPx, numPx, T)
theta <- ( atan2(matY, matX) + pi) + phase
# Note that:
# r >= exp( (2 * pi * (rev-1) + theta) * alpha)
# and r <= exp( (2 * pi * rev + theta) * alpha)
# where rev is the number of revolutions
rev <- ceiling((log(matR)/ alpha - theta) /( 2 * pi))
phi <- rev * 2 * pi + theta
z <- sin(phi * beta)
if ( generateImage)
image( z, col = getCustomColors())
else
return (z)
}
###################################################################
