[R] Creating Movies with R

TEMPL Matthias Matthias.Templ at statistik.gv.at
Fri Sep 22 13:00:21 CEST 2006


There is also a write.gif function in package caTools.
See the example there.

Best, 
Matthias

> 
> 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
> 
> ______________________________________________
> R-help at stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide 
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>



More information about the R-help mailing list