#     R script to create the plots of medians in the paper
#     "Human hair from the wreck of HMS Erebus of the Franklin
#      Expedition, 1845: elemental chemistry revealed by double-
#      ablation LA-ICP-MS"

library(ggplot2)

# Read the data

dmedall  <- read.table("medians.txt", header = TRUE)
# Set the levels of factors to control the order of groups when plotted
dmedall$source <- factor(dmedall$source,
         levels = c("Erebus-human", "Erebus-hog", "Beechey", "reference"))
dmedall$element <- factor(dmedall$element, 
         levels = c("23Na", "24Mg", "27Al", "31P", "34S", "39K", "44Ca",
                    "51V", "52Cr", "55Mn", "57Fe", "59Co", "60Ni", 
                    "63Cu", "66Zn", "75As", "77Se", "88Sr", "111Cd",
                    "137Ba", "202Hg", "208Pb"))

# High Erebus (hog and human) values

el_list <- c("202Hg", "27Al", "51V", "57Fe", "59Co", "60Ni", "63Cu", "75As")
dmed  <- subset(dmedall,    element %in% el_list)
# Add a fake point to adjust the scales to avoid label over-printing
dmed <- rbind(c("51V", "380", "reference", 0.007), dmed)
dmed$median <- as.numeric(dmed$median)
dmed$fake_alpha = c(0, rep(1, nrow(dmed) - 1))
ggplot(dmed, aes(x = median, y = source)) +
   geom_jitter(size = 0.5, height = 0.1, alpha = 0) +
   geom_jitter(size = 0.5, height = 0.1, alpha = dmed$fake_alpha) +
   scale_x_log10() +
   theme(axis.text.x = element_text(hjust = 0.5)) +
   facet_wrap(~ element, ncol = 4, scales = "free_x")

# All groups elevated in comparison with reference values

el_list <- c("24Mg", "55Mn", "88Sr", "208Pb")
dmed <- subset(dmedall,    as.character(element) %in% el_list)
ggplot(dmed, aes(x = median, y = source)) +
   geom_jitter(size = 0.5, height = 0.1, alpha = 0) +
   geom_jitter(size = 0.5, height = 0.1) +
   scale_x_log10() +
   theme(axis.text.x = element_text(hjust = 0.7)) +
   facet_wrap(~ element, ncol = 4, scales = "free_x")

# The remaining elements

el_list <- c("111Cd", "137Ba", "23Na", "31P", "34S", "44Ca", "52Cr", "66Zn", "77Se", "39K")
dmed <- subset(dmedall, as.character(element) %in% el_list)
ggplot(dmed, aes(x = median, y = source)) +
   geom_jitter(size = 0.5, height = 0.1, alpha = 0) +
   geom_jitter(size = 0.5, height = 0.1) +
   scale_x_log10() +
   facet_wrap(~ element, ncol = 4, scales = "free_x")