setwd("")
library(caTools)

TDRdata <- read.table(file="TDRdata.txt", sep="\t", header=T, colClasses=c("POSIXct","numeric"))
# TDRdata.txt contains 2 columns: 
# "Date.Time" (POSIXct format) and "Depth" (in meters)

n <- nrow(TDRdata)

# calculation of the vertical speed
    numerator <- diff(TDRdata$Depth)
    denominator <- diff(as.numeric(TDRdata$Date.Time))
    Vertical.Speed <- c(0, -numerator / denominator)

# smoothing of the vertical speed over a 10 seconds window 
# (data temporal resolution = 2 seconds)
    TDRdata$Smooth.Vert.Speed <- runmean(Vertical.Speed,5)

# creation of the drift vectors: 
# one for the negative drift phases and the other for the positive ones
    TDRdata$drift.N <- TDRdata$drift.P <- logical(n)

# parameters of drift detection (to be adjusted according to the data resolution)

# "alpha" corresponds to the drift detection threshold (in m/s):
# under this vertical speed, the animal is thought to swim passively
    alpha <- 0.6 

# "beta" is the minimal duration for a passive phase to be considered as a drift phase:
# an animal is drifting only if its passive swimming lasted more than beta minutes
    beta <- 90 	# 3min = 180 sec = 90 data points
    
# "gamma" corresponds to the maximal duration tolerated outside the threshold:
# if the animal swims over the threshold during more than gamma seconds,
# the drift phase is considered to be over
    gamma <- 4	# 4 data points = 8 seconds
	
# "delta" is the homogeneity parameter:
# the animal is considered drifting when the vertical speed variance
# over the passive phase is inferior to delta
    delta <- 0.05
    
# Detection of the Negative Drift Phases
    Neg <- vector("numeric",length=n)
    sel <- which((TDRdata$Smooth.Vert.Speed < 0) & (TDRdata$Smooth.Vert.Speed > -alpha))
    Neg[sel] <- TRUE
    diff.N <- c(0,diff(Neg))
    # location of the beginning of the negative passive phases:
    start.N <- which(diff.N > 0) 
    if(start.N[length(start.N)]==n)){
    start.N <- start.N[-length(start.N)] }
    # location of the end of the negative passive phases:
    end.N <- which(diff.N < 0)

    # passive phases separeted by less than gamma seconds are pooled together:
    dum <- start.N[2:length(start.N)] - end.N[1:(length(start.N)-1)]
    sel.short.seq <- which(dum < (gamma + 1))
    start.N <- start.N[-(sel.short.seq+1)]
    end.N <- end.N[-sel.short.seq]

    # passive phases lasting less than beta minutes are discarded
    start.drift.N <- start.N[which((end.N - start.N) > beta)]
    end.drift.N <- end.N[which((end.N - start.N) > beta)]

    for (j in 1:length(start.drift.N)){
     TDRdata$drift.N[start.drift.N[j]:end.drift.N[j]] <- TRUE }

# Detection of the Positive Drift Phases  
    Pos <- vector("numeric",length=n)
    sel <- which((TDRdata$Smooth.Vert.Speed > 0) & (TDRdata$Smooth.Vert.Speed < alpha))
    Pos[sel] <- TRUE
    diff.P <- c(0,diff(Pos))
    # location of the beginning of the positive passive phases:
    start.P <- which(diff.P > 0) 
    if(start.P[length(start.P)] == n)){ start.P <- start.P[-length(start.P)] }
    # location of the end of the positive passive phases:
    end.P <- which(diff.P < 0)

    # passive phases separeted by less than gamma seconds are pooled together:
    dum <- start.P[2:length(start.P)] - end.P[1:(length(start.P)-1)]
    sel.short.seq <- which(dum < (gamma + 1))
    start.P <- start.P[-(sel.short.seq+1)]
    end.P <- end.P[-sel.short.seq]

    # passive phases lasting less than beta minutes are discarded
    start.drift.P <- start.P[which((end.P - start.P) > beta)]
    end.drift.P <- end.P[which((end.P - start.P) > beta)]

    for (j in 1:length(start.drift.P)){
     	TDRdata$drift.P[start.drift.P[j]:end.drift.P[j]] <- TRUE 
    }