In EarthSystemDiagnostics/spectraluncertainty:

knitr::opts_chunk$set(echo = TRUE, cache = TRUE)

library(spectraluncertainty)
library(knitr)
library(tidyverse)

ModelSpectrum

ms1 <- ModelSpectrum(latitude = 45)
plot(ms1$freq, ms1$spec, log = "xy", type = "l")

Proxy Error Spectrum

Default parameters for proxy error spectrum

pars.1 <- GetSpecPars(proxy.type = "Mg_Ca")

kable(t(data.frame(pars.1)))

PES.1 <- do.call(ProxyErrorSpectrum, pars.1)

PlotSpecError(PES.1)

pars.2 <- GetSpecPars(proxy.type = "Mg_Ca", clim.spec.fun = "ModelSpectrum",
                    clim.spec.fun.args = list("latitude" = 1), delta_t = 2, tau_r = 2, T = 1e05)
PES.2 <- do.call(ProxyErrorSpectrum, pars.2)
PlotSpecError(PES.2)

Get error variance after smoothing to a given timescale.

spec.integral <- IntegrateErrorSpectra(PES.1, method = "sincfilter")
#spec.integral <- IntegrateErrorSpectra(PES.1, method = "cumsum")
#PlotTSDVariance(spec.integral)

kable(spec.integral)

GetProxyError(spec.integral, timescale = 100)

PlotTSDVariance(spec.integral)
PlotTSDVariance(spec.integral,
                include.constant.errors = TRUE)

Expected correlation

PES.MgCa <- RenameSpecComponents(do.call(ProxyErrorSpectrum, GetSpecPars(proxy.type = "Mg_Ca")))
PES.MgCa.300 <- RenameSpecComponents(
  do.call(
    ProxyErrorSpectrum, GetSpecPars(proxy.type = "Mg_Ca", N = 300)))

PES.Uk37 <- RenameSpecComponents(do.call(ProxyErrorSpectrum, GetSpecPars(proxy.type = "Uk37")))

ExpectedCorrelation(IntegrateErrorSpectra(PES.MgCa, method = "sincfilter"))
ExpectedCorrelation(IntegrateErrorSpectra(PES.MgCa.300, method = "sincfilter"))
ExpectedCorrelation(IntegrateErrorSpectra(PES.Uk37, method = "sincfilter"))

GetExpectedCorrelation <- function(proxy.type, lon, lat, tau_b){

   orbital.pars <- RelativeAmplitudeModulation(latitude = lat,
                                                  maxTimeKYear = 23,
                                                  minTimeKYear = 1,
                                                  bPlot = FALSE)

   amps <- AmpFromLocation(longitude = lon, latitude = lat,
                          depth.upr = 0, depth.lwr = -50,
                          proxy.type = proxy.type)

   seas.amp.pars <- list(
     mean.amp = signif(amps$mean.amp, 3), sig.sq_c = signif(amps$sig.sq_c, 3),
     sig.sq_a = round(orbital.pars$sig.sq_a, 3))


  pars <- GetSpecPars(proxy.type = proxy.type, 
                      tau_b = tau_b,
                      sig.sq_a = seas.amp.pars$sig.sq_a,
                      sig.sq_c = seas.amp.pars$sig.sq_c,
                      clim.spec.fun = "ModelSpectrum",
                      clim.spec.fun.args = list("latitude" = lat))

  PES <- do.call(ProxyErrorSpectrum, pars)
  PES <- RenameSpecComponents(PES)
  ExpectedCorrelation(IntegrateErrorSpectra(PES, method = "sincfilter", tau_smooth = 100))
}

GetExpectedCorrelation("Mg_Ca", 1.5, 1.5, tau_b = 1e03 * 10 / 5)

tmp <- breitkreuz.amp %>%
  select(longitude, latitude) %>% 
  distinct() %>% 
  filter(ceiling(latitude) %% 4 == 0,
         ceiling(longitude) %% 4 == 0) %>% 
  group_by(longitude, latitude) %>% 
  do({
  data.frame(corr = GetExpectedCorrelation("Mg_Ca", lat = .$latitude, lon = .$longitude,
                                           tau_b = 1e03 * 10/25))
  })


p <- tmp %>% 
  ggplot(aes(x = longitude, y = latitude, fill = corr)) +
  geom_raster()
p

hist(tmp$corr)
plot(corr~latitude, data = tmp)
plot(corr~longitude, data = tmp)

tmp %>% 
  ggplot(aes(x = latitude, y = corr, group = longitude, color = longitude)) +
  geom_line()