dat4bugsDN <- function(data, dir = 'lat', usek = F)
{
#
#	Code to convert typical ARGOS data into proper format for Winbugs
#
#	Created by Ian Jonsen
#
#		created on: 04/01/2022
#	  last modified on: 01/04/2022
#Note data file must include following named columns in the first 5 positions: date, time, lc, lat, lon
# 		
#date must be a numeric object of class "dates times", this can be accomplished as follows:
#  data$date <- chron(as.character(data$date), format='d/m/y', out.format='m/d/y')
#  where format (eg, ='d/m/y') must be set to format the date is in, and out.format is set to 'm/d/y'
#
#time must be converted to "minutes after midnight" format, this can be accomplished as follows:
#  data$time <- as.numeric(times(as.character(data$time))) * 1440				
#
#Try to convert date and time to correct format if not already so
require(chron)
if(is.factor(data$date) || is.character(data$date)) data$date <- chron(as.character(data$date), format = 'd/m/y', out.format='m/d/y')
else if(class(data$date)[1] != "dates") {stop('Error - date is not in correct format, see comments in dat4bugs')}
if(!is.numeric(data$time)) data$time <- as.numeric(times(as.character(data$time))) * 1440
else if(max(data$time) > 1440 || min(data$time) < 0){stop('Error - time is not in correct format, see comments in dat4bugs')}
#fill in missing dates with NA's
newdata <- fill.missing.dates(data)
time <- newdata$time
sunset <- newdata$sunset
sunrise <- newdata$sunrise
day <- newdata$day
lon <- newdata$lon
lat <- newdata$lat
date <- newdata$date
date.double <- rep(unique(date), each = 2)
row.out<-list(NULL)
i <- 1
j <- 1
# sort locations into day and night timesteps
if(day[1] == 2){
	while(i <= length(date.double) & j <= (length(date.double))){
		if(i == 1){
			if(time[1] > sunset[1]){
				row.out[[1]] <- row.names(newdata)[day == 2 & (date == date.double[1] & time > sunset | date == date.double[3] & time < sunrise)]
				row.out[[2]] <- row.names(newdata)[day == 1 & date == date.double[3]]
				j <- j + 2
				}
			else{ 
				row.out[[1]] <- row.names(newdata)[day == 2 & date == date.double[1] & time < sunrise]
				row.out[[2]] <- row.names(newdata)[day == 1 & date == date.double[1]]
				}
			}
		else{
			row.out[[i]] <- row.names(newdata)[day == 2 & (date == date.double[j-1] & time > sunset | date == date.double[j] & time < sunrise) & !is.na(day)]
			row.out[[i+1]] <- row.names(newdata)[day == 1 & date == date.double[j] & !is.na(day)]
			}
		j <- j + 2
		i <- i + 2
		}
	}
else{
	while(i <= length(date.double) - 1){
		if(i == 1){
			row.out[[1]] <- row.names(newdata)[day == 1 & date == date.double[1]]
			row.out[[2]] <- row.names(newdata)[day == 2 & (date == date.double[1] | (date == date.double[3] & time < sunrise))]
			}
		else if(i > 1 && i < length(date.double) - 1){
			row.out[[i]] <- row.names(newdata)[day == 1 & date == date.double[i]]
			row.out[[i+1]] <- row.names(newdata)[day ==2 & (date == date.double[i+1] & time > sunset | date == date.double[i+2] & time < sunrise)]
			}
		else if (i > 1 && i == length(date.double) - 1){
			row.out[[i]] <- row.names(newdata)[day == 1 & date == date.double[i]]
			row.out[[i+1]] <- row.names(newdata)[day ==2 & date == date.double[i+1] & time < sunset]
			}
		i <- i + 2
		}
	}	
tmp <- sapply(1:length(row.out),function(i){if(length(row.out[[i]]) == 0) row.out[[i]] <- NA else row.out[[i]]})
idx <- cumsum(c(1,sapply(tmp,length)))
nd.len <- c()
day.out <- c()
# Calculate night and day lengths in decimal hours
for(i in 1:length(tmp)){
	if(length(unique(newdata[tmp[[i]],]$date)) > 1){
		if(newdata[tmp[[i]][1],]$day == 2 && !is.na(newdata[tmp[[i]][1],]$day[1])){ nd.len[i] <- (1440 - newdata[tmp[[i]][1],]$sunset + newdata[tmp[[i]][nrow(newdata[tmp[[i]],])],]$sunrise) / 60 }
		else { nd.len[i] <- (newdata[tmp[[i]][1],]$sunset - newdata[tmp[[i]][nrow(newdata[tmp[[i]],])],]$sunrise) / 60 }
		}
	else{
		## determine sunset time for previous day to determine start of current night period
		sunset1 <- sunriseset(newdata[tmp[[i]][1],]$date - 1, newdata[tmp[[i]][1],]$lat, newdata[tmp[[i]][1],]$lon, 0)$sunset/60
		if(newdata[tmp[[i]][1],]$day == 2 && !is.na(newdata[tmp[[i]][1],]$day[1])){ nd.len[i] <- (1440 - sunset1 + newdata[tmp[[i]][nrow(newdata[tmp[[i]],])],]$sunrise)/60 }
		else { nd.len[i] <- (sunset1 - newdata[tmp[[i]][nrow(newdata[tmp[[i]],])],]$sunrise) / 60 }
		}	
	# Create index for 'regular' day-night intervals
	day.out[i] <- newdata[tmp[[i]][1],]$day
	}
# fill in NA's in day.out with correct day or night designation 
if(day.out[1] == 1) day.out.ind <- rep(c(1,2), length = length(day.out))
else {day.out.ind <- rep(c(2,1), length = length(day.out))}
day.out[is.na(day.out)] <- day.out.ind[is.na(day.out)]
newdata <- newdata[unlist(tmp),]
if(usek == T){
	# determine fraction of regular timestep at which each observation is made, k
	k <- c()
	for(i in 1:(length(idx) - 1)){
		if(day.out[i] == 2){
			tmp <- ((1440 - newdata$sunset[idx[i]:(idx[i+1]-1)]) + newdata$time[idx[i]:(idx[i+1]-1)])
			k1 <- ifelse(tmp > 1440, tmp - 1440, tmp) / (nd.len[i] * 60)
			}
		else {
			k1 <- (newdata$time[idx[i]:(idx[i+1]-1)] - newdata$sunrise[idx[i]:(idx[i+1]-1)]) / (nd.len[i] * 60)
			}
		k <- c(k, k1)
		k[k > 1 & k < 1.1] <- 1
		}
	}
# set up t-distribution parameters and create look-up values for the different location classess
sigma.lon <- c()
sigma.lat <- c()
nu.lon <- c()
nu.lat <- c()
# parameters for t-distributions (obtained by MLE); estimated from Vincent et al. 2002 ARGOS data
#   these values are published in Jonsen et al (submitted to Ecology 2004), appendix A, Table A1
sigma.lon[1] <- 0.2898660
sigma.lon[2] <- 0.3119293
sigma.lon[3] <- 0.9020423
sigma.lon[4] <- 2.1625936
sigma.lon[5] <- 0.5072920
sigma.lon[6] <- 4.2050261
sigma.lat[1] <- 0.1220553
sigma.lat[2] <- 0.2605126
sigma.lat[3] <- 0.4603374
sigma.lat[4] <- 1.607056
sigma.lat[5] <- 0.5105468
sigma.lat[6] <- 3.041276
nu.lon[1] <- 3.070609
nu.lon[2] <- 1.220822
nu.lon[3] <- 2.298819
nu.lon[4] <- 0.9136517
nu.lon[5] <- 0.786954
nu.lon[6] <- 1.079216
nu.lat[1] <- 2.075642
nu.lat[2] <- 6.314726
nu.lat[3] <- 3.896554
nu.lat[4] <- 1.010729
nu.lat[5] <- 1.057779
nu.lat[6] <- 1.331283
# add following code to restrict nu to be >= 2; required for Winbugs
nu.lon[nu.lon < 2] <- 2
nu.lat[nu.lat < 2] <- 2
# convert sigma's from km back to degrees
sigma.lon <- (sigma.lon/6367 * 180)/pi
sigma.lat <- (sigma.lat/6367 * 180)/pi
# convert standard errors to precisions (se^-2); required for Winbugs
itau1 <- sigma.lon[newdata$lc] ^ -2
itau2 <- sigma.lat[newdata$lc] ^ -2
# convert NA's in itau to poorest location class standard errors
itau1[is.na(itau1)] <- sigma.lon[6]
itau2[is.na(itau2)] <- sigma.lat[6]
# converts NA's in nu to smallest nu allowed by Winbugs
nu1 <- nu.lon[newdata$lc]
nu2 <- nu.lat[newdata$lc]
nu1[is.na(nu1)] <- 2
nu2[is.na(nu2)] <- 2
nd.len[is.na(nd.len)] <- 12
if(usek == T){
	# converts NA's in k to 0.5
	k[is.na(k)] <- 0.5
	}
if(dir == 'lon' && usek == F) return(list(y = newdata$lon, itau = itau1, nu = nu1, day = day.out, idx = idx, h = nd.len, RegN = length(idx) - 1))
else if(dir == 'lat' && usek == F) return(list(y = newdata$lat, itau = itau2, nu = nu2, day = day.out, idx = idx, h = nd.len, RegN = length(idx) - 1))
else if(dir == '2d' && usek == T) return(list(y = cbind(newdata$lon, newdata$lat), itau = cbind(itau1, itau2), nu = cbind(nu1, nu2), day = day.out, idx = idx, h = nd.len, k = k, RegN = length(idx) - 1))
else if(dir == '2d' && usek == F) return(list(y = cbind(newdata$lon, newdata$lat), itau = cbind(itau1, itau2), nu = cbind(nu1, nu2), day = day.out, idx = idx, h = nd.len, RegN = length(idx) - 1))
}