library(ggplot2)
library(latex2exp)

# Rs Seawater is the Standard for stable Isotope analysis of water
# The ratios are given with their errors
R18O = 2005.2E-6
R18Oe= 0.45E-6
R2H  = 155.76E-6
R2He= 0.05E-6
d18Oe = 0.1
d2He  = 1.0

# Majoube 18O/16O, equilibrium coefficients for water - vapour system
D18O = 1.137
E18O = -0.42
F18O = -2.07

# 2H/H
D2H =  24.844
E2H = -76.25
F2H =  52.61

# Temperature C & K - define a temperature range
TC  = seq(0,35,5)

# Calculate equilibrium fractionation
equilibrium <- function(TC){
  # convert TC to TK
  TK <- TC*0
  TK  = TC+273.15
  lnta18O = D18O*10^6/TK^2+E18O*10^3/TK+F18O 
  lnta2H  = D2H *10^6/TK^2+ E2H*10^3/TK+F2H 
  lntas   = data.frame(TC,TK,lnta18O,lnta2H)

  # ln(alpha)*1000 -> alpha
  alpha18O    = exp(lnta18O/1000)
  alpha2H     = exp(lnta2H/1000)
  lntas$alpha18O <- alpha18O   
  lntas$alpha2H  <- alpha2H

  # R(Start) -> R(Sample)
  S18O = R18O * alpha18O
  S2H  = R2H  * alpha2H
  lntas["S18O"] <- S18O                 # Add new column to data
  lntas["S2H"]  <- S2H                  # Add new column to data

  # Rs -> delta
  d18O = (R18O- S18O) /R18O*1000
  d2H  = (R2H - S2H)  / R2H*1000

  lnta  <- lntas                        # Replicate example data
  lntas <- cbind(lnta, d18O = d18O)     # Add new column to data
  lnta  <- lntas                        # Replicate example data
  lntas <- cbind(lnta, d2H  = d2H)     # Add new column to data
  return(lntas)
}

lntas <- equilibrium(TC)
attach(lntas)

ggplot(data=lntas, mapping=aes(x=TC, y=lnta18O)) +
  # ggtitle("Equilibrium fractionation of water") +
  geom_line(linetype="dashed",  size = 0.5, color="blue") +
  geom_point(shape=21, fill="white", color="blue", size=3) +
  xlab("Temperature (Kelvin)") +
  ylab(TeX("$1000 \\cdot ln(\\alpha)$")) 

ggplot(data=lntas, mapping=aes(x=d18O, y=d2H)) +
  # ggtitle("Equilibrium fractionation of water") +
  geom_line(linetype="dashed",  size = 0.5, color="blue") +
  geom_point(shape=21, fill="white", color="blue", size=3) +
  geom_errorbar(aes(ymin=(d2H-d2He), ymax=(d2H+d2He)), color="grey") +
  geom_errorbarh(aes(xmin = d18O-d18Oe, xmax = d18O+d18Oe), color="grey") +
  xlab(TeX("\\delta 18O $")) +
  ylab(TeX("\\delta 2H  $")) +
  geom_text(label=TC, nudge_x=0.15, nudge_y=0.1, check_overlap=T)