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

delta_t = 100
T = 101*delta_t


tmp <- OrbitalError(tau_s = 1e03*1/5, tau_b = 1e03*10/5, tau_p = 6/12,
             delta_t = 100, T = 1001*100,
             sig.sq_c = VarSine(4), phi_c = pi, delta_phi_c = 0, nu_c = 1,
             sig.sq_a = 0, phi_a = pi/2, nu_a = 1/23e03,
             N = 1)

tmp.long <- tmp %>% 
  gather(component, spec, -nu)

tmp.long %>% 
  #filter(nu > 0) %>% 
  ggplot(aes(x = nu, y=spec, colour = component)) +
  geom_line() +
  scale_y_continuous(trans = "log10") +
  scale_x_continuous(trans = "log10")

sqrt(sum(tmp$SB * diff(tmp$nu)[1]))
sqrt(sum(tmp$SU4 * diff(tmp$nu)[1]))
sqrt(VarSine(4))

tmp %>% filter(nu == 0)

sqrt(162276*diff(tmp$nu)[1])

delta_t = 100
T = 100*pi*delta_t

n.dt <- round(T / delta_t)
n.dt[n.dt %% 2 == 0] <- n.dt + 1
T <- n.dt * delta_t

OrbitalError <- function(nu = NULL,
                         tau_s, tau_b, tau_p, delta_t, T,
                         sig.sq_c, phi_c, delta_phi_c, nu_c,
                         sig.sq_a, phi_a, nu_a, N){
  
  if (is.null(nu))  nu <- GetNu(T = T, delta_t = delta_t)
  
  fax <- nu
  ### time scale parameters (all the same units - in the following examples: years)
  taus <- tau_s
  taub <- tau_b
  taup <- tau_p
  
  Deltat <- delta_t
  T <- T     ## T must be an odd multiple of Deltat
  
  ### amplitudes and phases of the seasonal cycle and the orbital modulation (nus in 1/time, here 1/years; phases in radians)
  
  sigc <- sqrt(sig.sq_c)
  phic.ev <- phi_c
  Deltaphic <- delta_phi_c
  nuc <- nu_c
  siga <- sqrt(sig.sq_a) 
  phia <- phi_a
  nua <- nu_a
  
  
  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)
  
  
  out <- data.frame(nu = nu,
                    SB = SB,
                    SU3 = SU3,
                    SU4=SU4)
  
  return(out)
  
}




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)


#png("paper.png", 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")

#dev.off()