moviefield<-function(lon,lat,array3d,start=1,end=10,delay=0.2, n=11,...){

# Produce a movie of a given three-dimensional array with p longitude points
# and q latitude points as the first two dimensions and n as the third time 
# dimension 
#
# Description:
#
#      Animates a sequence of maps 
#
# Usage:
#
#      moviefield(lon,lat,array3d,start=1,end=10,delay=0.2,n=11,...)
#
# Input:
#
#      array3d: a three-dimensional array with p longitude points and q latitude
#               points as the first two dimensions and n as the third time 
#               dimension 
#
#      lon: vector of longitudes
#
#      lat: vector of latitudes
#
#      start: first time step of the animation
#
#      end: last time step of the animation
#
#      delay: number of waiting seconds between maps
#
#      n: number of colours/intervals.
#      
#      ...: Additional arguments passed to `image'
#
# Authors:
#
#      David Stephenson <d.b.stephenson@reading.ac.uk> 16 June 2005
#      Christopher Ferro <c.a.t.ferro@reading.ac.uk>
#      Caio Coelho <c.a.d.s.coelho@reading.ac.uk> 
#      Dag Johan Steinskog <dag.johan.steinskog@nersc.no> 
#
# Examples:
#
#      x <- seq(-20, 20, 5)
#      y <- seq(30, 60, 5)
#      dim <- c(length(x), length(y), 1000)
#      data <- array(rnorm(prod(dim)), dim)
#      moviefield(x,y,data)
#      moviefield(x,y,data,start=1,end=50,delay=0.5,zlim=c(0.5,0.8))
#



rg<-range(array3d, na.rm=T)

lowerlim <- rg[1]
upperlim <- rg[2]

maximum <- max(abs(c(lowerlim,upperlim)))

if(lowerlim<0) breaks <- seq(-maximum,maximum,(maximum-(-maximum))/n)
else breaks <- seq(lowerlim,upperlim,(upperlim-(lowerlim))/n)


wait <- function(time) {
  now <- start <- proc.time()[3]
  while(now < (start + time)) now <- proc.time()[3]
}


# check if lon is from 0 to 360 or -180 to 180 to use appropriate world map
if (min(lon)<0){
coast <- map('world',interior = FALSE,resolution = 1) # Low resolution world map (lon -180 to 180)
}
else{
coast <- map('world2',interior = FALSE,resolution = 1) # Low resolution world map (lon 0 360)
}


# check if range includes negative values to use appropriate colour scale
if (min(array3d,na.rm=T) <0) {
colours = bluered(seq(0,n-1,by=1), "linear", yellow =TRUE)
}
else {
#colours = rev(heat.colors(n+1))
colours = rev(heat.colors(n))
}

par(mai=c(2, 1, 2, 0.5))
#par(mfrow = c(1, 1), mar = c(1, 1, 1, 1), pty = "m")
for(i in start:end) {
  image(lon, lat, array3d[,,i], col = colours, xlab='', ylab='', axes = FALSE,breaks=breaks,...)  
  map(coast, add = TRUE)
  title(paste("Frame",i,"of",end,sep=" "))
  box()
  wait(delay)
}


}