[R] Creating Movies with R

[Lorenzo Isella]

Dear All,

I'd like to know if it is possible to create animations with R.
To be specific, I attach a code I am using for my research to plot
some analytical results in 3D using the lattice package. It is not
necessary to go through the code.
Simply, it plots some 3D density profiles at two different times
selected by the user.
I wonder if it is possible to use the data generated for different
times to create something like an .avi file.

Here is the script:

rm(list=ls())
library(lattice)

# I start defining the analytical functions needed to get the density
as a function of time

expect_position <- function(t,lam1,lam2,pos_ini,vel_ini)
{1/(lam1-lam2)*(lam1*exp(lam2*t)-lam2*exp(lam1*t))*pos_ini+
1/(lam1-lam2)*(exp(lam1*t)-exp(lam2*t))*vel_ini
}

sigma_pos<-function(t,q,lam1,lam2)
{
q/(lam1-lam2)^2*(
(exp(2*lam1*t)-1)/(2*lam1)-2/(lam1+lam2)*(exp(lam1*t+lam2*t)-1) +
(exp(2*lam2*t)-1)/(2*lam2) )
}

rho_x<-function(x,expect_position,sigma_pos)
{
1/sqrt(2*pi*sigma_pos)*exp(-1/2*(x-expect_position)^2/sigma_pos)
}

#### Now the physical parameters
tau<-0.1
beta<-1/tau
St<-tau ### since I am in dimensionless units and tau is already in
units of 1/|alpha|
D=2e-2
q<-2*beta^2*D
############### Now the grid in space and time
time<-5  # time extent
tsteps<-501 # time steps
newtime<-seq(0,time,len=tsteps)
#### Now the things specific for the dynamics along x
lam1<- -beta/2*(1+sqrt(1+4*St))
lam2<- -beta/2*(1-sqrt(1+4*St))
xmin<- -0.5
xmax<-0.5
x0<-0.1
vx0<-x0
nx<-101 ## grid intervals along x
newx<-seq(xmin,xmax,len=nx) # grid along x

# M1 <- do.call("g", c(list(x = newx), mypar))


mypar<-c(q,lam1,lam2)
sig_xx<-do.call("sigma_pos",c(list(t=newtime),mypar))
mypar<-c(lam1,lam2,x0,vx0)
exp_x<-do.call("expect_position",c(list(t=newtime),mypar))

#rho_x<-function(x,expect_position,sigma_pos)

#NB: at t=0, the density blows up, since I have a delta as the initial state!
# At any t>0, instead, the result is finite.
#for this reason I now redefine time by getting rid of the istant t=0
to work out
# the density


rho_x_t<-matrix(ncol=nx,nrow=tsteps-1)
for (i in 2:tsteps)
{mypar<-c(exp_x[i],sig_xx[i])
myrho_x<-do.call("rho_x",c(list(x=newx),mypar))
rho_x_t[ i-1, ]<-myrho_x
}

### Now I also define a scaled density

rho_x_t_scaled<-matrix(ncol=nx,nrow=tsteps-1)
for (i in 2:tsteps)
{mypar<-c(exp_x[i],sig_xx[i])
myrho_x<-do.call("rho_x",c(list(x=newx),mypar))
rho_x_t_scaled[ i-1, ]<-myrho_x/max(myrho_x)
}

###########Now I deal with the dynamics along y

lam1<- -beta/2*(1+sqrt(1-4*St))
lam2<- -beta/2*(1-sqrt(1-4*St))
ymin<- 0
ymax<- 1
y0<-ymax
vy0<- -y0

mypar<-c(q,lam1,lam2)
sig_yy<-do.call("sigma_pos",c(list(t=newtime),mypar))
mypar<-c(lam1,lam2,y0,vy0)
exp_y<-do.call("expect_position",c(list(t=newtime),mypar))


# now I introduce the function giving the density along y: this has to
include the BC of zero
# density at wall

rho_y<-function(y,expect_position,sigma_pos)
{
1/sqrt(2*pi*sigma_pos)*exp(-1/2*(y-expect_position)^2/sigma_pos)-
1/sqrt(2*pi*sigma_pos)*exp(-1/2*(y+expect_position)^2/sigma_pos)
}

newy<-seq(ymin,ymax,len=nx) # grid along y with the same # of points
as the one along x


rho_y_t<-matrix(ncol=nx,nrow=tsteps-1)
for (i in 2:tsteps)
{mypar<-c(exp_y[i],sig_yy[i])
myrho_y<-do.call("rho_y",c(list(y=newy),mypar))
rho_y_t[ i-1, ]<-myrho_y
}

rho_y_t_scaled<-matrix(ncol=nx,nrow=tsteps-1)
for (i in 2:tsteps)
{mypar<-c(exp_y[i],sig_yy[i])
myrho_y<-do.call("rho_y",c(list(y=newy),mypar))
rho_y_t_scaled[ i-1, ]<-myrho_y/max(myrho_y)
}


# The following 2 plots are an example of the plots I'd like to use to
make an animation


g <- expand.grid(x = newx, y = newy)

instant<-100
mydens<-rho_x_t[ instant, ]%o%rho_y_t[ instant, ]/(max(rho_x_t[
instant, ]%o%rho_y_t[ instant, ]))


lentot<-nx^2
dim(mydens)<-c(lentot,1)

g$z<-mydens
jpeg("dens-t-3.jpeg")
print(wireframe(z ~ x * y, g, drape = TRUE,shade=TRUE,
scales = list(arrows = FALSE),pretty=FALSE, aspect = c(1,1), colorkey = TRUE
,zoom=0.8, main=expression("Density at t=2"), zlab =
list(expression("density"),rot = 90),distance=0.0,
perspective=TRUE,#screen = list(z = 150, x = -55,y= 0)
,zlim=range(c(0,1))))
dev.off()


instant<-300
mydens<-rho_x_t[ instant, ]%o%rho_y_t[ instant, ]/(max(rho_x_t[
instant, ]%o%rho_y_t[ instant, ]))


lentot<-nx^2
dim(mydens)<-c(lentot,1)

g$z<-mydens
jpeg("dens-t-3.jpeg")
print(wireframe(z ~ x * y, g, drape = TRUE,shade=TRUE,
scales = list(arrows = FALSE),pretty=FALSE, aspect = c(1,1), colorkey = TRUE
,zoom=0.8, main=expression("Density at t=3"), zlab =
list(expression("density"),rot = 90),distance=0.0,
perspective=TRUE,#screen = list(z = 150, x = -55,y= 0)
,zlim=range(c(0,1))))
dev.off()




Kind Regards

Lorenzo