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 = 1/12,
             delta_t = 100, T = 101*100,
             sig.sq_c = VarSine(4), phi_c = 0, delta_phi_c = pi, nu_c = 1,
             sig.sq_a = 0.4, phi_a = pi/2, nu_a = 1/23e03,
             N = 1)

tmp.long <- tmp[["freq"]] %>%
  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")


tmp.long <- tmp[["time"]] %>%
  gather(component, sigma, -tax)

tmp.long %>%
  #filter(nu > 0) %>%
  ggplot(aes(x = tax, y=sigma, colour = component)) +
  geom_line() +
  facet_wrap(~component, scales = "free_y")
  #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




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")