spectraluncertainty:

#!/usr/bin/Rscript

### time scale parameters (all the same units - in the following examples: years)

taus <- 100
taub <- 2000
taup <- 1/3

Deltat <- 100
T <- 45100     ## T must be an odd multiple of Deltat

### numbers (no units)

N <- 1

### amplitudes and phases of the seasonal cycle and the orbital modulation (nus in 1/time, here 1/years; phases in radians)

sigc <- sqrt(1/2)
phic.ev <- 0
Deltaphic <- 2*pi/6
nuc <- 1
siga <- sqrt(1/2) * 1/2
phia <- pi/2
nua <- 1 / (23000)

######### From here on: computation of bias and uncertainty components

nh <- (T/Deltat - 1) / 2     ## defined 4 lines after Eq. (94)

tax <- seq(from=-(nh * Deltat), to=+(nh * Deltat), by=Deltat)     ## time axis

fax <- seq(-1/2+1/(2*T/Deltat), 1/2-1/(2*T/Deltat), length.out=T/Deltat) / Deltat     ## frequency axis

sinc <- function(x){     ## the sinc function
    s <- ifelse(x != 0, sin(pi*x) / (pi*x), 1)
    return(s)
}

asinc <- function(x, Tt, Dt){     ## the aliased sinc function, defined 3 lines before Eq. (95)
    s <- ifelse(x != 0, sin(pi*x*Tt) / (sin(pi*x*Dt) * Tt/Dt), 1)
    return(s)
}

phib1 <- 2*pi*nua*taub - atan(2*pi*nua*taub)     ## defined by Eq. (77)
phib2 <- 4*pi*nua*taub - atan(4*pi*nua*taub)     ## defined by Eq. (A17)

Mb1 <- 1 / sqrt(1 + (2*pi*nua*taub)^2)     ## defined by Eq. (78)
Mb2 <- 1 / sqrt(1 + (4*pi*nua*taub)^2)     ## defined by Eq. (A18)

snc.cp  <- sinc(nuc*taup)     ## some frequently used variables to make the lines shorter below
snc.2cp <- sinc(2*nuc*taup)
snc.as  <- sinc(nua*taus)
snc.2as <- sinc(2*nua*taus)
snc.D   <- sinc(Deltaphic / (2*pi))
snc.2D  <- sinc(2*Deltaphic / (2*pi))

cos.c <- cos(phic.ev)     ## ...some more...
cos.2c <- cos(2*phic.ev)

gamma1 <- cos.c * snc.D * snc.cp     ## defined by Eq. (88)
gamma2 <- 1 - snc.D^2 + cos.2c * (snc.2D - snc.D^2)     ## defined by Eq. (90)

orb1   <- cos(2*pi*nua*tax +   phia+phib1)     ## the time dependent cosines that appear in Eq. (A12); this first one appears also in Eq. (80)
orb2d  <- cos(4*pi*nua*tax + 2*phia+phib2)
orb2dd <- cos(4*pi*nua*tax + 2*phia+phib1)

And <- 1 + siga*sqrt(2) * Mb1 * snc.as * orb1     ## defined by Eq. (80)

VB0   <- sigc^2 * (1 - snc.cp^2 + cos.2c * snc.2D * (snc.2cp - snc.cp^2)) + sigc^2*siga^2 * (1 - Mb1^2 * snc.as^2 * snc.cp^2 + cos.2c * snc.2D * (snc.2cp - Mb1^2 * snc.as^2 * snc.cp^2))     ## the VB terms, defined by Eqs. (A13-A16)
VB1   <- sigc^2*siga*sqrt(8) * (Mb1 * snc.as * (1 - snc.cp^2) + cos.2c * snc.2D * Mb1 * snc.as * (snc.2cp - snc.cp^2))
VB2d  <- sigc^2*siga^2 * (Mb2 * snc.2as * (1 + cos.2c * snc.2D * snc.2cp))
VB2dd <- -sigc^2*siga^2 * (Mb1^2 * snc.as^2 * snc.cp^2 * (1 + cos.2c * snc.2D))

B <- sigc*sqrt(2) * gamma1 * And     ## the bias, defined by Eq. (87)

U3.sq <- sigc^2 * gamma2 * And^2     ## third squared uncertainty component, defined by Eq. (89)

var.Bnj.1   <- VB1   * orb1     ## preparing for the next step...
var.Bnj.2d  <- VB2d  * orb2d
var.Bnj.2dd <- VB2dd * orb2dd

var.Bnj <- VB0 + var.Bnj.1 + var.Bnj.2d + var.Bnj.2dd     ## defined by Eq. (A12)

var.Bnj.mean <- VB0 + var.Bnj.1[nh+1] * sinc(nua*T) + var.Bnj.2d[nh+1] * sinc(2*nua*T) + var.Bnj.2dd[nh+1] * sinc(2*nua*T)     ## defined by Eq. (A19)

U4.sq <- var.Bnj / N     ## fourth squared uncertainty component, defined by Eq. (92)

dm <- rep(0, T/Deltat)     ## the single-argument Kronecker delta, defined one line after Eq. (99)
dm[nh+1] <- 1

xi.p <- asinc(fax + nua, T, Deltat)     ## defined by Eq. (99)
xi.m <- asinc(fax - nua, T, Deltat)

Sc <- dm     ## defined by Eq. (97), left
Sca <- dm * siga*sqrt(2) * Mb1 * snc.as * 2 * xi.p[nh+1] * cos(phia+phib1)     ## defined by Eq. (97), right
Sa <- (siga^2/2) * Mb1^2 * snc.as^2 * (xi.p^2 + xi.m^2 + 2 * xi.p * xi.m * cos(2 * (phia+phib1)))     ## defined by Eq. (98)
S0 <- T * (Sc + Sca + Sa)     ## defined by Eq. (96)

SB <- 2 * sigc^2 * gamma1^2 * S0     ## defined by Eq. (100)

SU3 <- sigc^2 * gamma2 * S0     ## defined by Eq. (101)

SU4 <- rep(Deltat/N * var.Bnj.mean, T/Deltat)     ## defined by Eq. (102)


postscript("paper.ps", horizontal=TRUE)

par(mfcol=c(1, 2))

### time series

plot(tax, B, type="l", lwd=1, col="blue", ylim=c(-2*sigc*sqrt(2), 2*sigc*sqrt(2)))     ## bias
lines(tax, sqrt(U3.sq), type="l", lwd=1, col="cyan")     ## bias uncertainty
lines(tax, sqrt(U4.sq), type="l", lwd=1, col="green")     ## white noise
lines(tax, sqrt(U3.sq+U4.sq), type="l", lwd=2, lty=2, col="black")     ## total scatter
lines(tax, sqrt(B^2+U3.sq+U4.sq), type="l", lwd=2, lty=2, col="gray")     ## squared error, defined by Eq. (23), but without the X-contributions to U

### power spectra

#plot(fax, SB, type="l", lwd=1, col="blue", log="y", ylim=c(1e-3, 1e5))     ## linear frequency axis
plot(fax, SB, type="l", lwd=1, col="blue", log="xy", xlim=c(1/T, 1/(2*Deltat)), ylim=c(1e-3, 1e5))     ## logarithmic frequency axis
lines(fax, SU3, type="l", lwd=1, col="cyan")
lines(fax, SU4, type="l", lwd=1, col="green")
lines(fax, SU3+SU4, type="l", lwd=2, lty=2, pch=19, col="black")
lines(fax, SB+SU3+SU4, type="l", lwd=2, lty=2, pch=19, col="gray")