oce: Analysis of Oceanographic Data

Documented in read.adv.nortek

# vim:textwidth=80:expandtab:shiftwidth=4:softtabstop=4

# abbreviations:
#   SIG     = System Integrator Guide
#   SIG2014 = system-integrator-manual_Dec2014_jan.pdf
#   IMU     = http://files.microstrain.com/3DM-GX3-35-Data-Communications-Protocol.pdf

#' @template readAdvTemplate
#'
#' @template encodingIgnoredTemplate
#'
#' @param haveAnalog1 A logical value indicating whether the data file has 'analog1' data.
#'
#' @param haveAnalog2 A logical value indicating whether the data file has 'analog2' data.
#'
#' @param type A string indicating which type of Nortek device produced the data file, `vector`
#' or `aquadopp`.
#'
#' @param header A logical value indicating whether the file starts with a header.
#' (This will not be the case for files that are created by data loggers that
#' chop the raw data up into a series of sub-files, e.g. once per hour.)
#'
#' @template adReadingMethodTemplate
read.adv.nortek <- function(
    file, from = 1, to, by = 1, tz = getOption("oceTz"),
    header = TRUE, longitude = NA, latitude = NA, encoding = NA, type = c("vector", "aquadopp"),
    haveAnalog1 = FALSE, haveAnalog2 = FALSE, debug = getOption("oceDebug"), monitor = FALSE, processingLog = NULL) {
    if (missing(file)) {
        stop("must supply 'file'")
    }
    if (is.character(file)) {
        if (!file.exists(file)) {
            stop("cannot find file \"", file, "\"")
        }
        if (0L == file.info(file)$size) {
            stop("empty file \"", file, "\"")
        }
    }
    if (!interactive()) {
        monitor <- FALSE
    }
    #   vvd=vector velocity data [p35 SIG], containing the data: pressure, vel, amp, corr (plus sensemble counter etc)
    #   vsd=velocity system data [p36 SIG], containing times, temperatures, angles, etc
    # NOTE: we interpolate from vsd to vvd, to get the final data$time, etc.
    type <- match.arg(type)
    oceDebug(debug, "read.adv.nortek(file=\"", file, "\", from=", format(from),
        ", to=", if (!missing(to)) format(to) else "(missing)",
        ", by=", by, ", tz=\"", tz, "\", header=", header, ",
        longitude=", longitude, ", latitude=", latitude, ", type=\"",
        type, "\", debug=", debug, ", monitor=", monitor, ",
        processingLog=(not shown)) {\n",
        sep = "", unindent = 1
    )
    if (is.numeric(by) && by < 1) {
        stop("cannot handle negative 'by' values")
    }
    if (by != 1) {
        warning("'by' argument only applies to fast-scale items such as velocity; temperature, etc. are not affected")
    }
    if (is.numeric(from) && from < 1) {
        stop("argument \"from\" must be 1 or larger")
    }
    if (!missing(to) && is.numeric(to) && to < 1) {
        stop("argument \"to\" must be 1 or larger")
    }
    if (is.character(file)) {
        filename <- fullFilename(file)
        file <- file(file, "rb")
        on.exit(close(file))
    }
    if (!inherits(file, "connection")) {
        stop("argument `file' must be a character string or connection")
    }
    if (!isOpen(file)) {
        filename <- "(connection)"
        open(file, "rb")
        on.exit(close(file))
    }
    if (!header) {
        stop("header must be TRUE")
    }
    oceDebug(debug, "  read.adv.nortek() about to read header\n")
    oceDebug(debug, "  read.adv.nortek() finished reading header\n")
    # find file length
    seek(file, 0, "end")
    fileSize <- seek(file, 0, "start")
    oceDebug(debug, "  fileSize=", fileSize, "\n")
    buf <- readBin(file, "raw", fileSize)
    header <- decodeHeaderNortek(buf, type = "vector", debug = debug - 1)
    if (debug > 1) {
        # Note: need high debugging to get this
        cat("\nheader is as follows:\n")
        str(header)
    }
    res <- new("adv")
    res@metadata$manufacturer <- "nortek"
    res@metadata$instrumentType <- "vector"
    res@metadata$filename <- filename
    res@metadata$longitude <- longitude
    res@metadata$latitude <- latitude
    res@metadata$numberOfSamples <- NA # filled in later
    res@metadata$numberOfBeams <- header$head$numberOfBeams # FIXME: check that this is correct
    res@metadata$measurementStart <- NA # FIXME
    res@metadata$measurementEnd <- NA # FIXME
    res@metadata$samplingRate <- 512 / header$user$averagingInterval # FIXME: why 512?
    res@metadata$serialNumber <- header$hardware$serialNumber
    res@metadata$internalCodeVersion <- header$hardware$picVersion
    res@metadata$softwareVersion <- header$user$swVersion
    res@metadata$hardwareRevision <- header$hardware$hwRevision
    res@metadata$recSize <- header$hardware$recSize
    res@metadata$velocityRange <- header$hardware$velocityRange
    res@metadata$firmwareVersion <- header$hardware$fwVersion
    res@metadata$hardwareConfiguration <- header$hardware$config
    res@metadata$configPressureSensor <- header$head$configPressureSensor
    res@metadata$configMagnetometerSensor <- header$head$configMagnetometerSensor
    res@metadata$configTiltSensor <- header$head$configTiltSensor
    res@metadata$beamAngle <- 25 # FIXME: should read from file
    res@metadata$tiltSensorOrientation <- header$head$tiltSensorOrientation
    res@metadata$frequency <- header$head$frequency
    res@metadata$headSerialNumber <- header$head$headSerialNumber
    res@metadata$bin1Distance <- header$user$blankingDistance # FIXME: is this right?
    res@metadata$blankingDistance <- header$user$blankingDistance
    res@metadata$measurementInterval <- header$user$measurementInterval
    res@metadata$transformationMatrix <- header$head$transformationMatrix
    res@metadata$deploymentName <- header$user$deploymentName
    res@metadata$cellSize <- header$user$cellSize
    res@metadata$velocityScale <- header$user$velocityScale
    res@metadata$originalCoordinate <- header$user$originalCoordinate
    res@metadata$oceCoordinate <- header$user$originalCoordinate
    res@metadata$oceBeamUnspreaded <- FALSE
    res@metadata$comments <- header$user$comments
    if (is.null(processingLog)) {
        processingLog <- paste(deparse(match.call()), sep = "", collapse = "")
    }
    hitem <- processingLogItem(processingLog)
    # Find the focus time by bisection, based on "sd" (system data, containing a time).
    vsdStart <- NULL # prevent scope warning from rstudio; defined later anyway
    bisectNortekVectorSd <- function(tFind, add = 0, debug = 0) {
        # tFind=time add=offset debug=debug
        oceDebug(debug, "\n")
        oceDebug(debug, "bisectNortekVectorSd(tFind=", format(tFind), ", add=", add, ", debug=", debug, ")\n")
        vsdLen <- length(vsdStart)
        lower <- 1
        upper <- vsdLen
        passes <- floor(10 + log(vsdLen, 2)) # won't need this many; only do this to catch coding errors
        for (pass in 1:passes) {
            middle <- floor((upper + lower) / 2) # nolint (no space before opening parenthesis)
            t <- ISOdatetime(2000 + bcdToInteger(buf[vsdStart[middle] + 8]), # year
                bcdToInteger(buf[vsdStart[middle] + 9]), # month
                bcdToInteger(buf[vsdStart[middle] + 6]), # day
                bcdToInteger(buf[vsdStart[middle] + 7]), # hour
                bcdToInteger(buf[vsdStart[middle] + 4]), # min
                bcdToInteger(buf[vsdStart[middle] + 5]), # sec
                tz = tz
            )
            if (tFind < t) {
                upper <- middle
            } else {
                lower <- middle
            }
            if (upper - lower < 2) {
                break
            }
            oceDebug(
                debug, "examine: t=", format(t), " middle=", middle,
                " lower=", lower, " upper=", upper, " pass=", pass, " of max=", passes, "\n"
            )
        }
        middle <- middle + add
        if (middle < 1) {
            middle <- 1
        }
        if (middle > vsdLen) {
            middle <- vsdLen
        }
        t <- ISOdatetime(2000 + bcdToInteger(buf[vsdStart[middle] + 8]), # year
            bcdToInteger(buf[vsdStart[middle] + 9]), # month
            bcdToInteger(buf[vsdStart[middle] + 6]), # day
            bcdToInteger(buf[vsdStart[middle] + 7]), # hour
            bcdToInteger(buf[vsdStart[middle] + 4]), # min
            bcdToInteger(buf[vsdStart[middle] + 5]), # sec
            tz = tz
        )
        oceDebug(debug, "result: t=", format(t), " at vsdStart[", middle, "]=", vsdStart[middle], "\n")
        return(list(index = middle, time = t)) # index is within vsd
    }
    # NOTE: system.time() indicates 0.2s to scan a 100Meg file [macpro desktop, circa 2009]
    # "vvd" stands for "Vector Velocity Data" [bottom of p35 of SIG]
    vvdStart <- .Call("locate_byte_sequences_old", buf, c(0xa5, 0x10), 24, c(0xb5, 0x8c), 0)
    vvdStartNew <- locateByteSequences(buf, c(0xa5, 0x10), 24, c(0xb5, 0x8c), 0)
    # FIXME <issue 2201> keep this test for a while
    if (!identical(vvdStart, as.integer(vvdStartNew))) {
        message("IMPORTANT: read.adv.nortek/locateByteSequences/vvdStart problem -- please report at github.com/dankelley/oce/issues")
        warning("IMPORTANT: read.adv.nortek/locateByteSequences/vvdStart problem -- please report at github.com/dankelley/oce/issues")
    }
    # "vsd" stands for "Vector System Data" [p36 of SIG]
    vsdStart <- .Call("locate_byte_sequences_old", buf, c(0xa5, 0x11), 28, c(0xb5, 0x8c), 0)
    vsdStartNew <- locateByteSequences(buf, c(0xa5, 0x11), 28, c(0xb5, 0x8c), 0)
    if (!identical(vsdStart, as.integer(vsdStartNew))) {
        message("IMPORTANT: read.adv.nortek/locateByteSequences/vsdStart problem -- please report at github.com/dankelley/oce/issues")
        warning("IMPORTANT: read.adv.nortek/locateByteSequences/vsdStart problem -- please report at github.com/dankelley/oce/issues")
    }
    # "vvdh" stands for "Vector Velocity Data Header" [p35 of SIG]
    vvdhStart <- .Call("locate_byte_sequences_old", buf, c(0xa5, 0x12), 42, c(0xb5, 0x8c), 0)
    vvdhStartNew <- locateByteSequences(buf, c(0xa5, 0x12), 42, c(0xb5, 0x8c), 0)
    if (!identical(vvdhStart, as.integer(vvdhStartNew))) {
        message("IMPORTANT: read.adv.nortek/locateByteSequences/vvdhStart problem -- please report at github.com/dankelley/oce/issues")
        warning("IMPORTANT: read.adv.nortek/locateByteSequences/vvdhStart problem -- please report at github.com/dankelley/oce/issues")
    }
    # "imu" stands for 'inertial motion unit' [p30 SIG2014]
    #<2201 FIXME FIXME CODE CODE> # FIXME <issue 2201> keep this test for a while
    imuStart <- .Call("locate_vector_imu_sequences_old", buf)
    imuStartNew <- locateVectorImuSequences(buf)
    if (!identical(imuStart, as.integer(imuStartNew))) {
        message("IMPORTANT: read.adv.nortek/locatorVectorImuSequences/imuStart problem -- please report at github.com/dankelley/oce/issues")
        warning("IMPORTANT: read.adv.nortek/locatorVectorImuSequences/imuStart problem -- please report at github.com/dankelley/oce/issues")
    }
    haveIMU <- length(imuStart) > 0
    if (haveIMU) {
        IMUtype <- "unknown"
        if (buf[imuStart[1] + 5] == 0xc3) {
            IMUtype <- "c3"
        } else if (buf[imuStart[1] + 5] == 0xcc) {
            IMUtype <- "cc"
        } else if (buf[imuStart[1] + 5] == 0xd2) {
            IMUtype <- "d2"
        } else if (buf[imuStart[1] + 5] == 0xd3) {
            IMUtype <- "d3"
        } else {
            warning(
                "unknown IMU type, with 5th byte 0x", buf[imuStart[1] + 5],
                "; only 0xc3, 0xcc, 0xd2 and 0xd3 are recognized"
            )
        }
        IMUlength <- length(imuStart)
        B4 <- sort(c(imuStart, imuStart + 1, imuStart + 2, imuStart + 3))
        # Note: a "tick" of the internal timestamp clock is 16 microseconds [IMU p 78]
        if (IMUtype == "c3") {
            # desribed in [1C] of the refernces of ?read.adv
            res@data$IMUdeltaAngleX <- readBin(buf[B4 + 6], "numeric", size = 4, n = IMUlength, endian = "little")
            res@data$IMUdeltaAngleY <- readBin(buf[B4 + 10], "numeric", size = 4, n = IMUlength, endian = "little")
            res@data$IMUdeltaAngleZ <- readBin(buf[B4 + 14], "numeric", size = 4, n = IMUlength, endian = "little")
            res@data$IMUdeltaVelocityX <- readBin(buf[B4 + 18], "numeric", size = 4, n = IMUlength, endian = "little")
            res@data$IMUdeltaVelocityY <- readBin(buf[B4 + 22], "numeric", size = 4, n = IMUlength, endian = "little")
            res@data$IMUdeltaVelocityZ <- readBin(buf[B4 + 26], "numeric", size = 4, n = IMUlength, endian = "little")
            res@data$IMUrotation <- array(NA, dim = c(3, 3, IMUlength))
            res@data$IMUrotation[1, 1, ] <- readBin(buf[B4 + 30], "numeric", size = 4, n = IMUlength, endian = "little")
            res@data$IMUrotation[1, 2, ] <- readBin(buf[B4 + 34], "numeric", size = 4, n = IMUlength, endian = "little")
            res@data$IMUrotation[1, 3, ] <- readBin(buf[B4 + 38], "numeric", size = 4, n = IMUlength, endian = "little")
            res@data$IMUrotation[2, 1, ] <- readBin(buf[B4 + 42], "numeric", size = 4, n = IMUlength, endian = "little")
            res@data$IMUrotation[2, 2, ] <- readBin(buf[B4 + 46], "numeric", size = 4, n = IMUlength, endian = "little")
            res@data$IMUrotation[2, 3, ] <- readBin(buf[B4 + 50], "numeric", size = 4, n = IMUlength, endian = "little")
            res@data$IMUrotation[3, 1, ] <- readBin(buf[B4 + 54], "numeric", size = 4, n = IMUlength, endian = "little")
            res@data$IMUrotation[3, 2, ] <- readBin(buf[B4 + 58], "numeric", size = 4, n = IMUlength, endian = "little")
            res@data$IMUrotation[3, 3, ] <- readBin(buf[B4 + 62], "numeric", size = 4, n = IMUlength, endian = "little")
            res@data$IMUtime <- readBin(buf[B4 + 66], "integer", size = 4, n = IMUlength, endian = "little") / 62500
            # test to show nortek the byte codes {
            #> for (ii in 1:2) {
            #>     message("IMU entry number: ", ii)
            #>     B <- imuStart[ii]
            #>     message("    starting byte (B, say) in file: ", B)
            #>     message("    byte[B+0  =", B+0, "]: 0x", buf[B+0], " ... check: should be 0xa5")
            #>     message("    byte[B+1  =", B+1, "]: 0x", buf[B+1], " ... check: should be 0x71")
            #>     message("    byte[B+5  =", B+5, "]: 0x", buf[B+5], " ... check: should be 0xc3")
            #>     message("    byte[B+66 =", B+66,"]: 0x", buf[B+66])
            #>     message("    byte[B+67 =", B+67,"]: 0x", buf[B+67])
            #>     message("    byte[B+68 =", B+68,"]: 0x", buf[B+68])
            #>     message("    byte[B+69 =", B+69,"]: 0x", buf[B+69])
            #>     lit <- readBin(buf[B+66:69], "integer", size=4, n=IMUlength, endian="little")/62500
            #>     big <- readBin(buf[B+66:69], "integer", size=4, n=IMUlength, endian="big")/62500
            #>     message("timestamp (in seconds) if little endian: ", lit)
            #>     message("timestamp (in seconds) if big endian:    ", big)
            #> }
            # } test to show nortek the byte codes
            res@metadata$IMUtype <- IMUtype
            res@metadata$units$IMUdeltaAngleX <- list(unit = expression(degree), scale = "")
            res@metadata$units$IMUdeltaAngleY <- list(unit = expression(degree), scale = "")
            res@metadata$units$IMUdeltaAngleZ <- list(unit = expression(degree), scale = "")
            res@metadata$units$IMUdeltaVelocityX <- list(unit = expression(m / s), scale = "")
            res@metadata$units$IMUdeltaVelocityY <- list(unit = expression(m / s), scale = "")
            res@metadata$units$IMUdeltaVelocityZ <- list(unit = expression(m / s), scale = "")
            res@metadata$units$IMUrotation <- list(unit = expression(), scale = "")
            res@metadata$units$IMUtime <- list(unit = expression(s), scale = "")
        } else if (IMUtype == "cc") {
            # described in [1B] of the references of ?read.adv
            # a "tick" of the internal timestamp clock is 16 microseconds [IMU p 78]
            res@data$IMUaccelX <- readBin(buf[B4 + 6], "numeric", size = 4, n = IMUlength, endian = "little")
            res@data$IMUaccelY <- readBin(buf[B4 + 10], "numeric", size = 4, n = IMUlength, endian = "little")
            res@data$IMUaccelZ <- readBin(buf[B4 + 14], "numeric", size = 4, n = IMUlength, endian = "little")
            res@data$IMUangrtX <- readBin(buf[B4 + 18], "numeric", size = 4, n = IMUlength, endian = "little")
            res@data$IMUangrtY <- readBin(buf[B4 + 22], "numeric", size = 4, n = IMUlength, endian = "little")
            res@data$IMUangrtZ <- readBin(buf[B4 + 26], "numeric", size = 4, n = IMUlength, endian = "little")
            res@data$IMUmagrtX <- readBin(buf[B4 + 30], "numeric", size = 4, n = IMUlength, endian = "little")
            res@data$IMUmagrtY <- readBin(buf[B4 + 34], "numeric", size = 4, n = IMUlength, endian = "little")
            res@data$IMUmagrtZ <- readBin(buf[B4 + 38], "numeric", size = 4, n = IMUlength, endian = "little")
            res@data$IMUrotation <- array(NA, dim = c(3, 3, IMUlength))
            res@data$IMUrotation[1, 1, ] <- readBin(buf[B4 + 42], "numeric", size = 4, n = IMUlength, endian = "little")
            res@data$IMUrotation[1, 2, ] <- readBin(buf[B4 + 46], "numeric", size = 4, n = IMUlength, endian = "little")
            res@data$IMUrotation[1, 3, ] <- readBin(buf[B4 + 50], "numeric", size = 4, n = IMUlength, endian = "little")
            res@data$IMUrotation[2, 1, ] <- readBin(buf[B4 + 54], "numeric", size = 4, n = IMUlength, endian = "little")
            res@data$IMUrotation[2, 2, ] <- readBin(buf[B4 + 58], "numeric", size = 4, n = IMUlength, endian = "little")
            res@data$IMUrotation[2, 3, ] <- readBin(buf[B4 + 62], "numeric", size = 4, n = IMUlength, endian = "little")
            res@data$IMUrotation[3, 1, ] <- readBin(buf[B4 + 66], "numeric", size = 4, n = IMUlength, endian = "little")
            res@data$IMUrotation[3, 2, ] <- readBin(buf[B4 + 70], "numeric", size = 4, n = IMUlength, endian = "little")
            res@data$IMUrotation[3, 3, ] <- readBin(buf[B4 + 74], "numeric", size = 4, n = IMUlength, endian = "little")
            res@data$IMUtime <- readBin(buf[B4 + 78], "integer", size = 4, n = IMUlength, endian = "little") / 62500
            res@metadata$IMUtype <- IMUtype
            res@metadata$units$IMUaccelX <- list(unit = expression(m / s^2), scale = "")
            res@metadata$units$IMUaccelY <- list(unit = expression(m / s^2), scale = "")
            res@metadata$units$IMUaccelZ <- list(unit = expression(m / s^2), scale = "")
            res@metadata$units$IMUangrtX <- list(unit = expression(degree / s), scale = "")
            res@metadata$units$IMUangrtY <- list(unit = expression(degree / s), scale = "")
            res@metadata$units$IMUangrtZ <- list(unit = expression(degree / s), scale = "")
            res@metadata$units$IMUmagrtX <- list(unit = expression(gauss), scale = "")
            res@metadata$units$IMUmagrtY <- list(unit = expression(gauss), scale = "")
            res@metadata$units$IMUmagrtZ <- list(unit = expression(gauss), scale = "")
            res@metadata$units$IMUrotation <- list(unit = expression(), scale = "")
            res@metadata$units$IMUtime <- list(unit = expression(s), scale = "")
        } else if (IMUtype == "d2") {
            # described in [1B] of the references of ?read.adv
            # a "tick" of the internal timestamp clock is 16 microseconds [IMU p 78]
            res@data$IMUaccelX <- readBin(buf[B4 + 6], "numeric", size = 4, n = IMUlength, endian = "little")
            res@data$IMUaccelY <- readBin(buf[B4 + 10], "numeric", size = 4, n = IMUlength, endian = "little")
            res@data$IMUaccelZ <- readBin(buf[B4 + 14], "numeric", size = 4, n = IMUlength, endian = "little")
            res@data$IMUangrtX <- readBin(buf[B4 + 18], "numeric", size = 4, n = IMUlength, endian = "little")
            res@data$IMUangrtY <- readBin(buf[B4 + 22], "numeric", size = 4, n = IMUlength, endian = "little")
            res@data$IMUangrtZ <- readBin(buf[B4 + 26], "numeric", size = 4, n = IMUlength, endian = "little")
            res@data$IMUmagrtX <- readBin(buf[B4 + 30], "numeric", size = 4, n = IMUlength, endian = "little")
            res@data$IMUmagrtY <- readBin(buf[B4 + 34], "numeric", size = 4, n = IMUlength, endian = "little")
            res@data$IMUmagrtZ <- readBin(buf[B4 + 38], "numeric", size = 4, n = IMUlength, endian = "little")
            res@data$IMUtime <- readBin(buf[B4 + 42], "integer", size = 4, n = IMUlength, endian = "little") / 62500
            res@metadata$IMUtype <- IMUtype
            res@metadata$units$IMUaccelX <- list(unit = expression(m / s^2), scale = "")
            res@metadata$units$IMUaccelY <- list(unit = expression(m / s^2), scale = "")
            res@metadata$units$IMUaccelZ <- list(unit = expression(m / s^2), scale = "")
            res@metadata$units$IMUangrtX <- list(unit = expression(degree / s), scale = "")
            res@metadata$units$IMUangrtY <- list(unit = expression(degree / s), scale = "")
            res@metadata$units$IMUangrtZ <- list(unit = expression(degree / s), scale = "")
            res@metadata$units$IMUmagrtX <- list(unit = expression(gauss), scale = "")
            res@metadata$units$IMUmagrtY <- list(unit = expression(gauss), scale = "")
            res@metadata$units$IMUmagrtZ <- list(unit = expression(gauss), scale = "")
            res@metadata$units$IMUrotation <- list(unit = expression(), scale = "")
            res@metadata$units$IMUtime <- list(unit = expression(s), scale = "")
        } else if (IMUtype == "d3") {
            # described in [1B] of the references of ?read.adv
            res@data$IMUdeltaAngleX <- readBin(buf[B4 + 6], "numeric", size = 4, n = IMUlength, endian = "little")
            res@data$IMUdeltaAngleY <- readBin(buf[B4 + 10], "numeric", size = 4, n = IMUlength, endian = "little")
            res@data$IMUdeltaAngleZ <- readBin(buf[B4 + 14], "numeric", size = 4, n = IMUlength, endian = "little")
            res@data$IMUdeltaVelocityX <- readBin(buf[B4 + 18], "numeric", size = 4, n = IMUlength, endian = "little")
            res@data$IMUdeltaVelocityY <- readBin(buf[B4 + 22], "numeric", size = 4, n = IMUlength, endian = "little")
            res@data$IMUdeltaVelocityZ <- readBin(buf[B4 + 26], "numeric", size = 4, n = IMUlength, endian = "little")
            res@data$IMUdeltaMagVectorX <- readBin(buf[B4 + 30], "numeric", size = 4, n = IMUlength, endian = "little")
            res@data$IMUdeltaMagVectorY <- readBin(buf[B4 + 34], "numeric", size = 4, n = IMUlength, endian = "little")
            res@data$IMUdeltaMagVectorZ <- readBin(buf[B4 + 38], "numeric", size = 4, n = IMUlength, endian = "little")
            res@data$IMUtime <- readBin(buf[B4 + 42], "integer", size = 4, n = IMUlength, endian = "little") / 62500
            res@metadata$IMUtype <- IMUtype
            res@metadata$units$IMUdeltaAngleX <- list(unit = expression(degree), scale = "")
            res@metadata$units$IMUdeltaAngleY <- list(unit = expression(degree), scale = "")
            res@metadata$units$IMUdeltaAngleZ <- list(unit = expression(degree), scale = "")
            res@metadata$units$IMUdeltaVelocityX <- list(unit = expression(m / s), scale = "")
            res@metadata$units$IMUdeltaVelocityY <- list(unit = expression(m / s), scale = "")
            res@metadata$units$IMUdeltaVelocityZ <- list(unit = expression(m / s), scale = "")
            res@metadata$units$IMUdeltaMagVectorRateX <- list(unit = expression(gauss), scale = "")
            res@metadata$units$IMUdeltaMagVectorRateY <- list(unit = expression(gauss), scale = "")
            res@metadata$units$IMUdeltaMagVectorRateZ <- list(unit = expression(gauss), scale = "")
            res@metadata$units$IMUtime <- list(unit = expression(s), scale = "")
        } else {
            warning("unsupported IMU type '", IMUtype, "'; only c3, cc, d2 and d3 are allowed")
        }
    }
    vvdhTime <- ISOdatetime(2000 + bcdToInteger(buf[vvdhStart + 8]), buf[vvdhStart + 9],
        buf[vvdhStart + 6], buf[vvdhStart + 7], buf[vvdhStart + 4], buf[vvdhStart + 5],
        tz = tz
    )
    vvdhRecords <- readBin(buf[sort(c(vvdhStart, vvdhStart + 1)) + 10], "integer", size = 2, n = length(vvdhStart), signed = FALSE, endian = "little")
    # Velocity scale.  Nortek's System Integrator Guide (p36) says
    # the velocity scale is in bit 1 of "status" byte (at offset 23)
    # in the Vector System Data header.  However, they seem to count
    # bits in the opposite way as oce does, so their bit 1 (starting
    # from 0) corresponds to our bit 7 (ending at 8).
    #
    # NOTE: the S.I.G. is confusing on the velocity scale, and this
    # confusion resulted in a change to the present code on
    # 2010-09-13.  Page 35 of S.I.G. clearly states that velocities
    # are in mm/s, which was used in the code previously.  However,
    # p44 contradicts this, saying that there are two possible scale
    # factors, namely 1mm/s and 0.1mm/s.  Starting on 2010-09-13, the
    # present function started using this possibility of two scale
    # factors, as determined in the next code line, following p36.
    res@metadata$velocityScale <- if ("0" == substr(byteToBinary(buf[vsdStart[1] + 23], endian = "big"), 7, 7)) 1e-3 else 0.1e-3
    oceDebug(
        debug, "velocityScale=", res@metadata$velocityScale, "m/s (from VSD header byte 24, 0x",
        as.raw(buf[vsdStart[1] + 23]), "(bit 7 of",
        byteToBinary(buf[vsdStart[1] + 23], endian = "big"), ")\n"
    )
    # Measurement start and end times.
    vsdLen <- length(vsdStart)
    res@metadata$measurementStart <- ISOdatetime(2000 + bcdToInteger(buf[vsdStart[1] + 8]), # year
        bcdToInteger(buf[vsdStart[1] + 9]), # month
        bcdToInteger(buf[vsdStart[1] + 6]), # day
        bcdToInteger(buf[vsdStart[1] + 7]), # hour
        bcdToInteger(buf[vsdStart[1] + 4]), # min
        bcdToInteger(buf[vsdStart[1] + 5]), # sec
        tz = tz
    )
    res@metadata$measurementEnd <- ISOdatetime(2000 + bcdToInteger(buf[vsdStart[vsdLen] + 8]), # year
        bcdToInteger(buf[vsdStart[vsdLen] + 9]), # month
        bcdToInteger(buf[vsdStart[vsdLen] + 6]), # day
        bcdToInteger(buf[vsdStart[vsdLen] + 7]), # hour
        bcdToInteger(buf[vsdStart[vsdLen] + 4]), # min
        bcdToInteger(buf[vsdStart[vsdLen] + 5]), # sec
        tz = tz
    )
    vvdLen <- length(vvdStart)
    res@metadata$measurementDeltat <- (as.numeric(res@metadata$measurementEnd) - as.numeric(res@metadata$measurementStart)) / (vvdLen - 1)
    toGiven <- !missing(to)
    if (!toGiven) {
        to <- length(vvdStart)
        oceDebug(debug, "No 'to' given, so using whole dataset: to=", to, "\n")
    }
    # Window data buffer, using bisection in case of a variable number of vd between sd pairs.
    if (inherits(from, "POSIXt")) {
        if (!inherits(to, "POSIXt")) {
            stop("if 'from' is POSIXt, then 'to' must be, also")
        }
        fromPair <- bisectNortekVectorSd(from, -1, debug - 1)
        from <- fromIndex <- fromPair$index
        toPair <- bisectNortekVectorSd(to, 1, debug - 1)
        to <- toIndex <- toPair$index
        byTime <- ctimeToSeconds(by)
        oceDebug(
            debug,
            "  from=", format(fromPair$t), " yields vsdStart[", fromIndex, "]\n",
            "  to  =", format(toPair$t), " yields vsdStart[", toIndex, "]\n",
            "  by=", by, "byTime=", byTime, "s\n",
            "vsdStart[", fromPair$index, "]=", vsdStart[fromPair$index], "at time", format(fromPair$t), "\n",
            "vsdStart[", toPair$index, "]=", vsdStart[toPair$index], "at time", format(toPair$t), "\n"
        )
        twoTimes <- ISOdatetime(2000 + bcdToInteger(buf[vsdStart[1:2] + 8]), # year
            bcdToInteger(buf[vsdStart[1:2] + 9]), # month
            bcdToInteger(buf[vsdStart[1:2] + 6]), # day
            bcdToInteger(buf[vsdStart[1:2] + 7]), # hour
            bcdToInteger(buf[vsdStart[1:2] + 4]), # min
            bcdToInteger(buf[vsdStart[1:2] + 5]), # sec  NOTE: nortek files lack fractional seconds
            tz = tz
        )
        vsd.dt <- as.numeric(twoTimes[2]) - as.numeric(twoTimes[1]) # FIXME: need # samplesPerBurst here
        # Next two lines suggest that readBin() can be used instead of bcdToInteger ... I imagine it would be faster
        # cat("month=", readBin(buf[vsdStart[1]+9], "integer", n=1, size=1, endian="little"), "(as readBin)\n")
        # cat("month=", bcdToInteger(buf[vsdStart[1]+9]), "(as bcd)\n")
        oceDebug(debug, "nrecords=", readBin(buf[vsdStart[1] + 10:11], "integer", n = 1, size = 2, endian = "little"), "\n")
        oceDebug(debug, "vsd.dt=", vsd.dt, "(from twoTimes)\n")
        vvdStart <- vvdStart[vsdStart[fromIndex] <= vvdStart & vvdStart <= vsdStart[toIndex]]
        # vvdDt <- vsd.dt * (toIndex - fromIndex) / length(vvdStart)
        # find vvd region that lies inside the vsd [from, to] region.
        # vvdStartFrom <- max(1, vvdStart[vvdStart < fromPair$index])
        # vvdStartTo   <- min(length(vvdStart), vvdStart[vvdStart > toPair$index])
    } else {
        # Window data buffer, using bisection in case of a variable number of vd between sd pairs.
        oceDebug(debug, "numeric values for args from=", from, "to=", to, "by=", by, "\n")
        fromIndex <- from
        toIndex <- to
        if (toIndex < 1 + fromIndex) {
            stop("need more separation between from and to")
        }
        oceDebug(debug, "fromIndex=", fromIndex, "toIndex=", toIndex, "\n")
        oceDebug(debug, vectorShow(vvdStart, "before subset, vvdStart is"))
        oceDebug(debug, vectorShow(vsdStart, "before subset, vsdStart is"))
        vvdStart <- vvdStart[fromIndex:toIndex]
        oceDebug(debug, vectorShow(vvdStart, "after  subset, vvdStart is"))
        # ensure that vvdStart pointers are bracketed by vsdStart pointers
        # FIXME: but will this invalidate fromIndex and toIndex?
        vsdStartFrom <- which(vvdStart[1] >= vsdStart)[1]
        vsdStartTo <- which(vsdStart >= vvdStart[length(vvdStart)])[1]
        oceDebug(debug, "vsdStartFrom=", vsdStartFrom, "and vsdStartTo=", vsdStartTo, "(before NA check)\n")
        if (is.na(vsdStartTo)) {
            vsdStartTo <- length(vsdStart)
        }
        oceDebug(debug, "vsdStartFrom=", vsdStartFrom, "and vsdStartTo=", vsdStartTo, "(after NA check)\n")
        vsdStart <- vsdStart[seq(vsdStartFrom, vsdStartTo)]
        oceDebug(debug, vectorShow(vsdStart, "after  subset, vsdStart is"))
    }
    oceDebug(
        debug, "about to trim vsdStart, based on vvdStart[1]=", vvdStart[1],
        " and vvdStart[length(vvdStart)]=", vvdStart[length(vvdStart)], "\n"
    )
    oceDebug(debug, vectorShow(vsdStart, "before trimming, vsdStart:"))
    oceDebug(debug, "from=", from, "to=", to, "\n")
    # Find spanning subset, expanded a little for now
    subsetStart <- head(which(vvdStart[1] < vsdStart), 1)
    if (subsetStart > 1) {
        subsetStart <- subsetStart - 1
    } # extend a bit (for now)
    subsetEnd <- tail(which(vsdStart < vvdStart[length(vvdStart)]), 1)
    oceDebug(debug, "first guess: subsetEnd=", subsetEnd, "\n")
    if (subsetEnd < length(vsdStart)) {
        subsetEnd <- subsetEnd + 1
    }
    oceDebug(debug, "try start vsdStart[subsetStart=", subsetStart, "] = ", vsdStart[subsetStart], "\n")
    oceDebug(debug, "try end   vsdStart[subsetEnd=  ", subsetEnd, "] = ", vsdStart[subsetEnd], "\n")
    oceDebug(debug, vectorShow(vsdStart, "before trimming, vsdStart:"))
    vsdStart <- vsdStart[seq(subsetStart, subsetEnd - 1, 1)]
    oceDebug(debug, vectorShow(vsdStart, "after  trimming, vsdStart:"))
    if (2 > length(vsdStart)) {
        stop("need at least 2 velocity-system-data chunks to determine the timing; try increasing the difference between 'from' and 'to'")
    }
    if (toIndex <= fromIndex) {
        stop("no data in specified range from=", format(from), " to=", format(to))
    }
    # we make the times *after* trimming, because this is a slow operation
    # NOTE: the ISOdatetime() call takes 60% of the entire time for this function.
    vsdTime <- ISOdatetime(2000 + bcdToInteger(buf[vsdStart + 8]), # year
        bcdToInteger(buf[vsdStart + 9]), # month
        bcdToInteger(buf[vsdStart + 6]), # day
        bcdToInteger(buf[vsdStart + 7]), # hour
        bcdToInteger(buf[vsdStart + 4]), # min
        bcdToInteger(buf[vsdStart + 5]), # sec
        tz = tz
    )
    oceDebug(debug, "reading Nortek Vector, and using timezone: ", tz, "\n")
    # update res@metadata$measurementDeltat
    res@metadata$measurementDeltat <- mean(diff(as.numeric(vsdTime)), na.rm = TRUE) * length(vsdStart) / length(vvdStart) # FIXME
    vsdLen <- length(vsdStart)
    vsdStart2 <- sort(c(vsdStart, 1 + vsdStart))
    voltage <- 0.1 * readBin(buf[vsdStart2 + 10], "integer", size = 2, n = vsdLen, signed = FALSE, endian = "little")
    heading <- 0.1 * readBin(buf[vsdStart2 + 14], "integer", size = 2, n = vsdLen, signed = TRUE, endian = "little")
    oceDebug(debug, vectorShow(heading, "heading"))
    pitch <- 0.1 * readBin(buf[vsdStart2 + 16], "integer", size = 2, n = vsdLen, signed = TRUE, endian = "little")
    oceDebug(debug, vectorShow(pitch, "pitch"))
    roll <- 0.1 * readBin(buf[vsdStart2 + 18], "integer", size = 2, n = vsdLen, signed = TRUE, endian = "little")
    oceDebug(debug, vectorShow(roll, "roll"))
    temperature <- 0.01 * readBin(buf[vsdStart2 + 20], "integer", size = 2, n = vsdLen, signed = TRUE, endian = "little")
    oceDebug(debug, vectorShow(temperature, "temperature"))
    salinity <- header$user$salinity
    oceDebug(debug, "salinity (in res@metadata):", salinity, "\n")
    # byte 22 is an error code
    # byte 23 is status, with bit 0 being orientation (p36 of Nortek's System Integrator Guide)
    status <- buf[vsdStart[floor(0.5 * length(vsdStart))] + 23]
    res@metadata$orientation <- if ("0" == substr(byteToBinary(status, endian = "big"), 1, 1)) "upward" else "downward"
    # FIXME: should read roll and pitch "out of range" or "OK" here, in bites 3 and 2
    # FIXME was wrong# res@metadata$burstLength <- round(length(vvdStart) / length(vsdStart), 0) # FIXME: surely this is in the header (?!?)
    # FIXME was wrong# oceDebug(debug, vectorShow(res@metadata$burstLength, "burstLength"))
    vvdStart2 <- sort(c(vvdStart, 1 + vvdStart))
    vvdLen <- length(vvdStart) # FIXME: should be subsampled with 'by' ... but how???
    if (haveAnalog1) {
        # FIXME: shouldn't this be auto-detected from 'USER' header?
        analog1 <- readBin(buf[sort(c(vvdStart + 8, vvdStart + 9))],
            "integer",
            n = vvdLen, size = 2, endian = "little", signed = FALSE
        )
    }
    if (haveAnalog2) {
        # FIXME: shouldn't this be auto-detected from 'USER' header?
        analog2 <- readBin(buf[sort(c(vvdStart + 2, vvdStart + 5))],
            "integer",
            n = vvdLen, size = 2, endian = "little", signed = FALSE
        )
    }
    p.MSB <- as.numeric(buf[vvdStart + 4])
    p.LSW <- readBin(buf[vvdStart2 + 6], "integer", size = 2, n = vvdLen, signed = FALSE, endian = "little")
    pressure <- (65536 * p.MSB + p.LSW) / 1000
    oceDebug(debug, vectorShow(pressure, "pressure"))
    v <- array(double(), dim = c(vvdLen, 3))
    v[, 1] <- res@metadata$velocityScale * readBin(buf[vvdStart2 + 10], "integer", size = 2, n = vvdLen, signed = TRUE, endian = "little")
    v[, 2] <- res@metadata$velocityScale * readBin(buf[vvdStart2 + 12], "integer", size = 2, n = vvdLen, signed = TRUE, endian = "little")
    v[, 3] <- res@metadata$velocityScale * readBin(buf[vvdStart2 + 14], "integer", size = 2, n = vvdLen, signed = TRUE, endian = "little")
    if (debug > 0.9) {
        oceDebug(debug, "v[", dim(v), "] begins...\n")
        print(matrix(as.numeric(v[1:min(3, vvdLen), ]), ncol = 3))
    }
    a <- array(raw(), dim = c(vvdLen, 3))
    a[, 1] <- buf[vvdStart + 16]
    a[, 2] <- buf[vvdStart + 17]
    a[, 3] <- buf[vvdStart + 18]
    if (debug > 0.9) {
        oceDebug(debug, "a[", dim(a), "] begins...\n")
        print(matrix(as.numeric(a[1:min(3, vvdLen), ]), ncol = 3))
    }
    q <- array(raw(), dim = c(vvdLen, 3))
    q[, 1] <- buf[vvdStart + 19]
    q[, 2] <- buf[vvdStart + 20]
    q[, 3] <- buf[vvdStart + 21]
    if (debug > 0.9) {
        cat("q[", dim(q), "] begins...\n")
        print(matrix(as.numeric(q[1:min(3, vvdLen), ]), ncol = 3))
    }
    # sec <- as.numeric(vsdTime) - as.numeric(vsdTime[1])
    # vds <- var(diff(sec))
    # BAD: vvdSec <- .Call("stutter_time", sec, 8)
    # vvdSec <- approx(seq(0, 1, length.out=length(vsdTime)), vsdTime, seq(0, 1, length.out=length(vvdStart)))$y
    # oceDebug(debug, vectorShow(vvdSec, "vvdSec"))
    oceDebug(debug, vectorShow(vsdStart, "vsdStart"))
    oceDebug(debug, vectorShow(vvdStart, "vvdStart"))
    rm(buf)
    gc()
    # subset using 'by'
    # by.orig <- by
    if (is.character(by)) {
        oceDebug(debug, "by=\"", by, "\" given as argument to read.adv.nortek()\n", sep = "")
        oceDebug(debug, " ... infer to be", ctimeToSeconds(by), "s\n")
        by <- ctimeToSeconds(by) / res@metadata$measurementDeltat
        oceDebug(debug, " ... so step by", by, "through the data\n")
    }
    len <- length(vvdStart)
    look <- seq(1, len, by = by)
    oceDebug(debug, "length(vvdStart)=", length(vvdStart), "\n")
    # vvdStart.orig <- vvdStart
    vvdStart <- vvdStart[look]
    oceDebug(debug, "length(vvdStart)=", length(vvdStart), "(after 'look'ing) with by=", by, "\n")
    # vvdSec <- vvdSec[look]
    pressure <- pressure[look] # only output at burst headers, not with velo (FIXME: huh??)
    v <- v[look, ]
    a <- a[look, ]
    q <- q[look, ]
    if (0 < sum(vvdhRecords)) {
        res@metadata$samplingMode <- "burst"
        # Note: if we knew that all bursts were of the same length, we could use the same method
        # as for the continuous case, specifying e.g. vvdhRecords[1] instead of 0.  But do we know that?
        # Also, what I'm doing here is probably fine, since bursts last an hour and so looping
        # also them won't be expensive.
        sss <- NULL
        for (b in seq_along(vvdhRecords)) {
            sss <- c(sss, as.numeric(vvdhTime[b]) + seq(0, by = 1 / res@metadata$samplingRate, length.out = vvdhRecords[b]))
        }
        time <- sss[look] + (vsdTime[1] - as.numeric(vsdTime[1]))
        delayForWarmup <- 2 + 1 / (res@metadata$samplingRate * 2) # FIXME: this is from a forum posting, not an official doc.
        time <- time + delayForWarmup
    } else {
        res@metadata$samplingMode <- "continuous"
        time <- numberAsPOSIXct(do_adv_vector_time(vvdStart, vsdStart, vsdTime, vvdhStart, vvdhTime, 0, res@metadata$samplingRate))
    }
    res@metadata$numberOfSamples <- dim(v)[1]
    res@metadata$numberOfBeams <- dim(v)[2]
    res@metadata$velocityResolution <- res@metadata$velocityScale / 2^15
    res@metadata$salinity <- salinity
    # FIXME: guess-based kludge to infer whether continuous or burst-mode sample
    res@data$v <- v
    res@data$a <- a
    res@data$q <- q
    res@data$time <- time
    res@data$pressure <- pressure
    res@data$timeBurst <- vvdhTime
    res@data$recordsBurst <- vvdhRecords
    res@data$voltageSlow <- voltage
    res@data$timeSlow <- vsdTime
    res@data$headingSlow <- heading
    res@data$pitchSlow <- pitch
    res@data$rollSlow <- roll
    res@data$temperatureSlow <- temperature
    if (haveAnalog1) {
        res@data$analog1 <- analog1
    }
    if (haveAnalog2) {
        res@data$analog2 <- analog2
    }
    res@metadata$velocityResolution <- res@metadata$velocityScale
    res@metadata$velocityMaximum <- res@metadata$velocityScale * 2^15
    res@metadata$units$v <- list(unit = expression(m / s), scale = "")
    res@processingLog <- unclass(hitem)
    res@metadata$units$v <- list(unit = expression(m / s), scale = "")
    res@metadata$units$pressure <- list(unit = expression(dbar), scale = "")
    res@metadata$units$headingSlow <- list(unit = expression(degree), scale = "")
    res@metadata$units$pitchSlow <- list(unit = expression(degree), scale = "")
    res@metadata$units$rollSlow <- list(unit = expression(degree), scale = "")
    res@metadata$units$temperatureSlow <- list(unit = expression(degree * C), scale = "")
    oceDebug(debug, "} # read.adv.nortek(file=\"", filename, "\", ...)\n", sep = "", unindent = 1)
    res
}