`ssmMeta` <- function (loc.list, model = "DCRWSmeta_cov", n.covar=1,...) {
# 
# Function ssmMeta - Takes output (loc.list) from dat4bugsCOV, sets up initial values, 
#   calls WinBUGS, and aggregates results. Implements 'DCRWSmeta_cov1oneway', 
#   a Bayesian state-space switching model with covariates influencing the 
#   behavioural switch from 'directed' into 'resident' state

    if (!model %in% "DCRWSmeta_cov") stop("model not implemented.")
    
    nt <- length(loc.list)
    y <- NULL
    itau2 <- NULL
    nu <- NULL
	  covar <-  NULL  
    j <- NULL
    RegN <- NULL
    id <- NULL
    first.date <- NULL
    tstep <- NULL
    for(i in 1:nt){
    	y <- rbind(y, loc.list[[i]]$y)
    	itau2 <- rbind(itau2, loc.list[[i]]$itau2)
    	nu <- rbind(nu, loc.list[[i]]$nu)
    	if (n.covar>0) {covar <- rbind(covar, loc.list[[i]]$covar)}
    	j <- c(j, loc.list[[i]]$j)
    	RegN <- c(RegN, loc.list[[i]]$RegN)
    	id <- c(id, loc.list[[i]]$id)
    	first.date <- c(first.date, loc.list[[i]]$first.date)
    	tstep <- c(tstep, loc.list[[i]]$tstep)
    	}
    tmp <- list(NULL)	
    for(i in 1:nt){
		  tmp[[i]] <- loc.list[[i]]$idx
   	}
 	  for(i in 2:nt){
	    tmp[[i]] <- max(tmp[[i-1]]) + tmp[[i]][-1]-1
 		}
 	
 	  extract <- function(k) nrow(k$y)
 		
 	  idx <- unlist(tmp)
	  Xidx <- cumsum(c(1, RegN-1)) 
	  Yidx <- cumsum(c(1, unlist(lapply(loc.list, extract))))
        
    x <- rbind(y[idx[-length(idx)], ], y[max(idx) - 1, ])    
    row.na <- which(is.na(x[, 1]))
    x[row.na, 1] <- approx(seq(nrow(x)), x[, 1], xout = row.na, 
            rule = 2)$y
    row.na <- which(is.na(x[, 2]))
    x[row.na, 2] <- approx(seq(nrow(x)), x[, 2], xout = row.na, 
            rule = 2)$y
    x <- x[-nrow(x),]
    row.na <- which(is.na(y[, 1]) | is.na(y[, 2]))
    if (length(row.na) == 0) {
        y.init <- y
        y.init[!is.na(y.init)] <- NA
    	} else {
       y.init <- y
       y.init[!is.na(y.init)] <- 9999
       yinit.xna <- approx(seq(nrow(y)), y[, 1], xout = row.na, 
                rule = 2)$y
       yinit.yna <- approx(seq(nrow(y)), y[, 2], xout = row.na, 
                rule = 2)$y
       y.init[row.na, ] <- cbind(yinit.xna, yinit.yna)
       y.init[y.init == 9999] <- NA
       }     
    tstep.sec <- tstep[1] * 86400 
    lidx.fn <- function(k) length(k$idx) - 1
    lidx <- sapply(loc.list, lidx.fn)
    steplims <- unlist(sapply(1:nt, function(i){
    				seq(first.date[i], by = tstep.sec, length = lidx[i])}))
    				
    #---------------------------------------------------------------------------
    # Setup data list for BUGS plus initial values 
    bugs.data <- list(y = y, itau2 = itau2, nu = nu, idx = idx, j = j, Xidx=Xidx, Yidx=Yidx, N=nt)         
    iSigma <- matrix(c(1, 0, 0, 1), 2, 2)
    bmode <- rep(1, max(Xidx)-1) 
    parameters <- c("Sigma", "x", "theta", "gamma", "bmode", "psi") # consistent across all DCRWS models

    inits1 <- list(iSigma = iSigma, gamma = c(NA, 0.2), 
    		dev=0.6, tmp = c(0.5, 0.5), lambda = c(0.5, NA), 
    				bmode = bmode, x = x, logpsi = rep(0, nt), y = y.init)
    
    inits2 <- list(iSigma = iSigma, gamma = c(NA, 0.1), 
    		dev=0.7, tmp = c(0.5, 0.5), lambda = c(0.5, NA), 
    		bmode = bmode, x = x, logpsi = rep(0, nt), y = y.init)
    
    # Parameters specific to covariate model
   	inits1$beta <- c(0.1,-0.1); inits2$beta <- c(0,0);  	      
    m <- matrix(0,2,ncol(covar)); m[2,] <- NA
   	bugs.data$covar <- covar
    inits1$m <- m ; inits2$m <- m 	
    parameters <- c(parameters, "beta", "m", "phi")   
	  inits <- list(inits1,inits2)
    	  
    # Call to BUGS
	  tmp <- bugs(data = bugs.data, inits = inits, parameters.to.save = parameters, ...)
    save("tmp",file="bugs.out.RData")    # save in case directly need sims.array 
    
    #---------------------------------------------------------------------------            
    ###Create output object called 'out'. This object contains the full McMC
    ###sample, convergence diagnostics, and a data frame with summary statistics
    ###for the estimated states. 
    out <- tmp$sims.list
    out$steplims <- steplims   
    out$id <- id
  	out$model <- model	
  	out$DIC <- tmp$DIC; out$pD <- tmp$pD    
  	out$'convdiag' <- tmp$summary[,c('Rhat','n.eff')]
  	
    lon.mean <- apply(out$x[, , 1], 2, mean)
  	lat.mean <- apply(out$x[, , 2], 2, mean)
  	lon.quants <- t(apply(out$x[, , 1], 2, quantile, c(0.025, 0.5, 0.975)))
  	lat.quants <- t(apply(out$x[, , 2], 2, quantile, c(0.025, 0.5, 0.975)))
  	newid <- list(NULL)
  	newid <- unlist(lapply(1:nt, function(i){rep(id[i], RegN[i]-1)}))		
  	# location estimate information
    out$summary <- data.frame(id=newid, 
      date=as.POSIXct(as.numeric(steplims), origin="1970-01-01 00:00:00", tz="GMT"),
  		lon=lon.mean, lat=lat.mean,
  		lon.025=lon.quants[, 1],
  		lat.025=lat.quants[, 1],
  		lon.5=lon.quants[, 2],
  		lat.5=lat.quants[, 2],
  		lon.975=lon.quants[, 3],
  		lat.975=lat.quants[, 3])
 		# behavioural state information
 		out$summary <- data.frame(out$summary, bmode=tmp$mean$bmode, bmode.5=tmp$median$bmode)   	
  	# covariate information
    out$summary$covar <- covar[-cumsum(RegN),]
   		# covar needs to be trimmed to RegN-length(N) timesteps to match ssm output	
  		phi1 <- tmp$median$phi[,1]; out$summary$phi.1 <- NA
  		phi2 <- tmp$median$phi[,2]; out$summary$phi.2 <- NA
      nobs <- diff(Xidx)-1 # phi estimates per animal each lack 1st obs
      for (i in 1:nt) { # pad each start with an NA  
  		  out$summary$phi.1[Xidx[i]:c(Xidx[i+1]-1)] <- c(NA,phi1[c(Xidx[i]-i+1):c(Xidx[i]-i+nobs[i])]) 
  		  out$summary$phi.2[Xidx[i]:c(Xidx[i+1]-1)] <- c(NA,phi2[c(Xidx[i]-i+1):c(Xidx[i]-i+nobs[i])]) 
      }
	out		
} # end function