oce: Analysis of Oceanographic Data

library(oce)

## See page 2 of Foreman 1977 for the format of tide3.dat,
## which is provided in t-tide as (what seems to be) an exact
## copy of the appendix in Foreman (1977).

## Notes
## 1. most of the tests here check against the T_TIDE values.
## 2. variable names are chosen to match T_TIDE, except that e.g. "const" is "tideconst"

file <- file("tide3.dat.gz", "r")

############################
## Constituents [const]   ##
############################

nc <- 146
name <- kmpr <- vector("character", nc)
freq <- ikmpr <- df <-
    d1 <- d2 <- d3 <- d4 <- d5 <- d6 <-
    semi <- isat <- nsat <- ishallow <- nshallow <- doodsonamp <- doodsonspecies <-
    vector("numeric", nc)
doodsonamp <- rep(NA, nc)
doodsonspecies <- rep(NA, nc)

ishallow <- NA
ic <- 1
while (TRUE) {
    items <- scan(file, "character", nlines=1, quiet=TRUE)
    nitems <- length(items)
    if (nitems == 0)
        break
    if (nitems != 2 && nitems != 3) stop("wrong number of entries on line", ic)
    name[ic] <- gsub(" ", "", items[1])
    freq[ic] <- as.numeric(items[2])
    kmpr[ic] <- if (nitems == 2) ""  else gsub(" ", "", items[3])
    ic <- ic + 1
}
for (ic in 1:nc) {
    if (kmpr[ic] != "")  {
        ikmpr[ic] <- which(name == kmpr[ic])
        df[ic]    <- freq[ic] - freq[ikmpr[ic]]
    }
}
df[1] <- 0                              # Z0 is special (???) BUG: fixme

############################
## Satellites [sat]       ##
############################
get.satellite <- function(x, o)
{
    ldel  <- as.numeric(substr(x, o+01, o+03))            # I3
    mdel  <- as.numeric(substr(x, o+04, o+06))            # I3
    ndel  <- as.numeric(substr(x, o+07, o+09))            # I3
    ph    <- as.numeric(substr(x, o+10, o+13))            # F4.2
    ee    <- as.numeric(substr(x, o+14, o+20))            # F7.4
    ir    <- as.numeric(substr(x, o+22, o+22))            # 1x, I1
    c(ldel, mdel, ndel, ph, ee, ir)
}

scan(file, "character", nlines=3, quiet=TRUE) # skip 3 lines

ns <- 162
deldood <- matrix(NA, ns, 3)
phcorr  <- matrix(NA, ns, 1)
amprat  <- matrix(NA, ns, 1)
ilatfac <- matrix(NA, ns, 1)
iconst  <- matrix(NA, ns, 1)

this.sat <- 1
while (TRUE) {
    x <- readLines(file, n=1)
    nx <- nchar(x)
    if (this.sat > ns || nx < 10) break
    ## Line format and content
    ## 6X,A5,1X,6I3,F5.2,I4
    ## kon, ii,jj,kk,ll,mm,nn semi nj
    kon <- gsub(" ", "", substr(x, 7, 11))
    which.constituent <- which(name == kon)

    d1[which.constituent] <- ii <- as.numeric(substr(x, 13, 15))
    d2[which.constituent] <- jj <- as.numeric(substr(x, 16, 18))
    d3[which.constituent] <- kk <- as.numeric(substr(x, 19, 21))
    d4[which.constituent] <- ll <- as.numeric(substr(x, 22, 24))
    d5[which.constituent] <- mm <- as.numeric(substr(x, 25, 27))
    d6[which.constituent] <- nn <- as.numeric(substr(x, 28, 30))

    semi[which.constituent] <- as.numeric(substr(x, 31, 35))
    nj <- as.numeric(substr(x, 36, 39)) # number of satellites
    nsat[which.constituent] <- nj

    if (nj > 0) {
        ## ALP1    1 -4  2  1  0  0 -.25   2
        ## ALP1  -1  0  0 .75 0.0360R1   0 -1  0 .00 0.1906
        is <- 1
        while (is <= nj) {
            xs <- readLines(file, n=1)
            nxs <- nchar(xs)
            if (nxs != 31 && nxs != 33 && nxs != 39 && nxs != 54 && nxs != 56 && nxs != 77 && nxs != 79) {
                cat("GOT BAD LINE AS FOLLOWS:\n12345678901234567890123456789012345678901234567890\n", xs, "\n", sep="")
                stop("need 31, 33, 39, 54, 56, 77 or 79 chars, but got ", nxs)
            }
            s <- get.satellite(xs, 11)
            deldood[this.sat, 1:3] <- s[1:3]
            phcorr[this.sat]  <- s[4]
            amprat[this.sat]  <- s[5]
            ilatfac[this.sat] <- s[6]
            iconst[this.sat]  <- which.constituent # constituent to which this satellite is attached
            this.sat <- this.sat + 1
            is <- is + 1
            if (nxs > 50) {
                s <- get.satellite(xs, 34)
                deldood[this.sat, 1:3] <- s[1:3]
                phcorr[this.sat]  <- s[4]
                amprat[this.sat]  <- s[5]
                ilatfac[this.sat] <- s[6]
                iconst[this.sat]  <- which.constituent # constituent to which this satellite is attached
                this.sat <- this.sat + 1
                is <- is + 1
            }
            if (nxs > 70) {
                s <- get.satellite(xs, 57)
                deldood[this.sat, 1:3] <- s[1:3]
                phcorr[this.sat]  <- s[4]
                amprat[this.sat]  <- s[5]
                ilatfac[this.sat] <- s[6]
                iconst[this.sat]  <- which.constituent # constituent to which this satellite is attached
                this.sat <- this.sat + 1
                is <- is + 1
            }
        }
    }
}
if ((this.sat - 1) != ns) stop("failed to read all ", ns, " satellite entries.  Only got ", this.sat)


sat <- list(deldood=deldood, phcorr=phcorr, amprat=amprat, ilatfac=ilatfac, iconst=iconst)


## This portion of the file ends as follows
##      M3      3  0  0  0  0  0 -.50   1
##      M3     0 -1  0 .50  .0564
## and see Foreman (1977) page 2 for how to read this.  Here,
## we'll check every aspect of the last satellite.

##
############################
## Shallow [shallow]      ##
############################

## (6X,A5,I1,2X,4(F5.2,A5,5X))
## KON = name of the shallow water constituent;
## NJ = number of main constituents from which it is derived;
## COEF,KONCO = combination number and name of these main constituents.

num.shallow <- 251
iconst <- vector("numeric", num.shallow)   # names as T_TIDE
coef   <- vector("numeric", num.shallow)
iname  <- vector("numeric", num.shallow)
this.shallow <- 1
while (TRUE) {
    x <- readLines(file, n=1)
    nx <- nchar(x)
    if (nx < 10) break
    kon <- gsub(" ", "", substr(x, 7, 11))
    which.constituent <- which(name == kon)
    nj <- as.numeric(substr(x, 12, 12))
    nshallow[which.constituent] <- nj
    if (nj > 0) {
        for (j in 1:nj) {
            o <- 15 + (j-1)*15
            iconst[this.shallow] <- which.constituent
            coef[this.shallow]  <- as.numeric(substr(x, o, o+4))
            konco <- gsub(" ", "", substr(x, o+5, o+9))
            iname[this.shallow] <- which(name == konco)
            ishallow[which.constituent] <- this.shallow - j + 1 # BUG: broken
            this.shallow <- this.shallow + 1
        }
    }
}
close(file)

shallow <- data.frame(iconst=iconst, coef=coef, iname=iname)


#################
## equilibrium ##
##################

efile <- file("t_equilib.dat.gz", "r")
edat <- readLines(efile)
ne <- length(edat)
for (i in 10:ne) {                      # 9 lines of header
    # Name Species A B
    kon <- gsub(" ", "", substr(edat[i], 1, 4))
    which.constituent <- which(name == kon)
    if (length(which.constituent) < 1) stop("cannot understand equilibirum constituent", kon)
    species <- as.numeric(substr(edat[i], 8, 8))
    A <- as.numeric(substr(edat[i],  9, 15))
    B <- as.numeric(substr(edat[i], 16, 21))
    if (A != 0) {
        doodsonamp[which.constituent] <- A / 1e5
        doodsonspecies[which.constituent] <- species
    } else {
        doodsonamp[which.constituent] <- B / 1e5
        doodsonspecies[which.constituent] <- -species
    }
}

close(efile)

const <- data.frame(name=name,
                    freq=freq,
                    kmpr=kmpr,
                    ikmpr=ikmpr,
                    df=df,
                    d1=d1, d2=d2, d3=d3, d4=d4, d5=d5, d6=d6,  # T_TIDE has these as matrix 'doodson'
                    semi=semi,
                    nsat=nsat,
                    ishallow=ishallow,
                    nshallow=nshallow,
                    doodsonamp=doodsonamp,
                    doodsonspecies=doodsonspecies,
                    stringsAsFactors=FALSE)
tidedata <- list(const=const, sat=sat, shallow=shallow)

test_that("deldood", {
          expect_equal(sum(nsat), ns)
          expect_equal(sat$deldood[1, ], c(-1,  0, 0))
          expect_equal(sat$deldood[2, ], c(0,  -1, 0))
          expect_equal(sat$deldood[3, ], c(-2, -2, 0))
})

test_that("constituents were read correctly", {
        expect_equal(name[1], "Z0")
        expect_equal(name[2], "SA")
        expect_equal(kmpr[2], "SSA")
        i <- which(name == "M2")
        expect_equal(freq[i], 1/12.4206011981605)
})

test_that("shallow constitutents", {
          shallow <- data.frame(iconst=iconst, coef=coef, iname=iname)
          expect_equal(tidedata$shallow$iconst[1:5], c(26, 26, 27, 27, 30))
          expect_equal(tidedata$shallow$coef[1:5], c(2, -1, 1, -1, 2))
          expect_equal(tidedata$shallow$iname[1:5], c(19, 13, 57, 13, 48))
          expect_equal(df[48], 0.08051140070000)
          expect_equal(ishallow[143:146], c(242, 245, 246, 248))
          expect_equal(nshallow[143:146], c(3, 1, 2, 4))
})

test_that("doodson", {
          expect_equal(tidedata$const$doodsonspecies[c(2, 3, 10)], c(0, 0, -1))
          expect_equal(tidedata$const$doodsonamp[c(2, 3, 10)], c(0.01160000000000, 0.07299000000000, 0.01153000000000))
})

test_that("ancillary code", {
    ## Test against matlab t_astron, with a randomly-picked calibration time
    t <- as.POSIXct("2008-01-22 18:50:24", tz="GMT")
    a <- tidemAstron(t)
    expect_equal(a$astro, c(1.28861316428, 0.33390620851, 0.83751937277, 0.14234854462, 0.08559663825, 0.78633079279), tolerance=1e-8)
    expect_equal(a$ader,  c(0.96613680803, 0.03660110127, 0.00273790931, 0.00030945407, 0.00014709388, 0.00000013082), tolerance=1e-8)
    vuf <- tidemVuf(t, 48)
    expect_equal(vuf$v, c(0.57722632857477), tolerance=1e-8)
    expect_equal(vuf$u, c(0), tolerance=1e-8)
    expect_equal(vuf$f, c(1), tolerance=1e-8)
    vuf <- tidemVuf(t, c(48, 49), 45)
    expect_equal(vuf$v, c(0.57722632857477, 0.62841490855698), tolerance=1e-8)
    expect_equal(vuf$u, c(0.00295677805220, 0.00180270946435), tolerance=1e-8)
    expect_equal(vuf$f, c(0.96893771510868, 0.98142639461951), tolerance=1e-8)
})