Nothing
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)
})
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.