install.packages("palinsol")

# install.packages("devtools")
devtools::install_github("special-uor/palinsol")

# Load the pipe operator from `magrittr`
`%>%` <- magrittr::`%>%`

# Create sequence (vector) with years
years <- seq(from = -10^5, to = 0, by = 1000)

# Find orbital parameters
orb_param <- years %>%
  purrr::map_df(palinsol::astro, solution = palinsol::ber78, degree = TRUE)
#> load BER78data

# Append the years to the table with the orbital parameters
orb_param <- orb_param %>%
  dplyr::mutate(year = years, .before = 1)

# Print table with the first 10 rows
orb_param %>%
  dplyr::slice(1:10) %>%
  knitr::kable()

orb_param %>%
  tidyr::pivot_longer(cols = c(eps, ecc, varpi, epsp),
                      names_to = "var") %>%
  ggplot2::ggplot(ggplot2::aes(year, value)) +
  ggplot2::geom_point() +
  ggplot2::geom_line() +
  ggplot2::facet_wrap(facets = ~var, 
                      scales = "free",
                      labeller = ggplot2::labeller(var = c("ecc" = "Eccentricity (ecc)",
                                                           "eps" = "Obliquity (eps)",
                                                           "varpi" = "True Solar Longitude of the Perihelion (varpi)",
                                                           "epsp" = "(epsp)"))) +
  ggplot2::theme_bw()

# Load the pipe operator from `magrittr`
`%>%` <- magrittr::`%>%`

# Generate a sequence (vector) with mid-month values
mid_month <- seq(from = 15.5, to = 345.5, by = 30)

# Calculate the True Solar Longitudes (TSL) for present-day mid-month values
tt_present <- 0.0
orbit_present <- palinsol::astro(t = tt_present, 
                                 solution = palinsol::ber78, 
                                 degree = FALSE)
mid_month_tsl_present <- palinsol::day2l(orbit = orbit_present, 
                                     day = mid_month)

# Calculate orbital parameters for the last 100 ka years
## Create sequence (vector) with years
years <- seq(from = -10^5, to = 0, by = 1000)

## Find orbital parameters
orb_param <- years %>%
  purrr::map_df(palinsol::astro, solution = palinsol::ber78, degree = FALSE)

## Append the years to the table with the orbital parameters
orb_param <- orb_param %>%
  dplyr::mutate(year = years, .before = 1)


# Calculate the insolation values at the location of interest
## Set the latitude of interest
lat <- 45.469951

insol_tbl <- mid_month_tsl_present %>%
  purrr::map(palinsol::Insol, orbit = orb_param, lat = lat * pi / 180, S0 = 1365) %>%
  magrittr::set_names(c("Jan", "Feb", "Mar", "Apr", "May", "Jun", "Jul", "Aug", "Sep", "Oct", "Nov", "Dec")) %>%
  tibble::as_tibble()

## Append the years to the table with the insolation values
insol_tbl <- insol_tbl %>%
  dplyr::mutate(year = years, .before = 1)

## Print table with the first 10 rows
insol_tbl %>%
  dplyr::slice(1:10) %>%
  knitr::kable()

insol_tbl %>%
  tidyr::pivot_longer(cols = 2:13, names_to = "month") %>%
  dplyr::mutate(year = year / 1000) %>%
  dplyr::filter(month %in% c("Jan")) %>%
  ggplot2::ggplot(ggplot2::aes(year, value)) +
  ggplot2::geom_point(size = 0.5) +
  ggplot2::geom_line(colour = "navyblue") +
  ggplot2::labs(x = "Year [kyr]", 
                y = "Insolation [W m-2]",
                title = "January's insolation for the last 100ka yrs") +
  ggplot2::scale_x_continuous(breaks = scales::pretty_breaks(10)) +
  ggplot2::theme_bw()

insol_tbl %>%
  tidyr::pivot_longer(cols = 2:13, names_to = "month") %>%
  dplyr::mutate(month = factor(month,
                               levels = c("Jan", "Feb", "Mar", "Apr", "May", "Jun", "Jul", "Aug", "Sep", "Oct", "Nov", "Dec")),
                season = ifelse(month %in% c("Mar", "Apr", "May"),
                                "Spring",
                                ifelse(month %in% c("Jun", "Jul", "Aug"),
                                       "Summer",
                                       ifelse(month %in% c("Sep", "Oct", "Nov"),
                                              "Autumn",
                                              "Winter"))),
                season = factor(season,
                                levels = c("Spring", "Summer", "Autumn", "Winter"))) %>%
  ggplot2::ggplot(ggplot2::aes(year, value, colour = month)) +
  ggplot2::geom_point(size = 0.5) +
  ggplot2::geom_line() +
  ggplot2::labs(x = "Year", y = "Insolation [W m-2]") + 
  ggplot2::facet_wrap(facets = ~season, nrow = 2) +
  ggplot2::theme_bw()

insol_tbl %>%
  dplyr::group_by(year) %>% 
  dplyr::summarise(value = sum(Jan:Dec)) %>%
  dplyr::mutate(value = value / 1000,
                year = year / 1000) %>%
  ggplot2::ggplot(ggplot2::aes(year, value)) +
  ggplot2::geom_point(size = 0.5) +
  ggplot2::geom_line(colour = "navyblue") +
  ggplot2::labs(x = "Year [kyr]", 
                y = "Insolation [kW m-2]",
                title = "Annual insolation for the last 100ka yrs") +
  ggplot2::scale_x_continuous(breaks = scales::pretty_breaks(10)) +
  ggplot2::scale_y_continuous(breaks = scales::pretty_breaks(10)) +
  ggplot2::theme_bw()