---
title: "Introduction to the rice package"
output:
  html_vignette:
    toc: true
    toc_depth: 3
    number_sections: true
vignette: >
  %\VignetteIndexEntry{Introduction to the rice package}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

# Introduction

Radiocarbon dating requires a range of calculations, for example calibration[^1][^2][^3][^4], translations between pMC, F14C, C14 age and D14C, and assessing the impacts of contamination. This package provides functions to do so in R. 

# Installation

On first usage of the package, it has to be installed: 

```{r, eval=FALSE}
install.packages('rice')
``` 

The companion data package 'rintcal' which has the radiocarbon calibration curves will be installed if it isn't already. New versions of R packages appear regularly, so please re-issue the above command regularly to remain up-to-date, or use:

```{r, eval=FALSE}
update.packages()
``` 

To obtain access to the calibration curves and radiocarbon functions, first the package has to be loaded:

```{r}
library(rice)
```

# Calibration curves

The calibration curves can be plotted:

```{r, fig.width=4, fig.asp=.8}
draw.ccurve()
```

Or, comparing two calibration curves, on the BC/AD scale:

```{r, fig.width=4, fig.asp=.8}
draw.ccurve(1000, 2020, BCAD=TRUE, cc2='marine20', add.yaxis=TRUE)
```

Or zooming in to between AD 1600 and 2000:

```{r, fig.width=4, fig.asp=.8}
draw.ccurve(1600, 1950, BCAD=TRUE)
```

Note the wiggles, or recurring C-14 ages? Calibration curves such as IntCal20 come with 3 columns: cal BP ($\theta$), the corresponding C-14 ages ($\mu$), and their errors ($\sigma$). Whereas each cal BP or BC/AD year has a single $\mu$ and $\sigma$, owing to wiggles in the calibration curves, any single $\mu$ will often have multiple $\theta$'s:

```{r}
BCADtoC14(1694) 
C14toBCAD(129)
```

Let's add these values to the previous plot:

```{r, fig.width=4, fig.asp=.8}
draw.ccurve(1600, 1950, BCAD=TRUE)
abline(h=BCADtoC14(1694)[1], lty=3) 
abline(v=C14toBCAD(129), lty=3)
```


Atmospheric nuclear tests around AD 1950-1964 caused a huge human-made C-14 peak (resulting in highly negative C-14 ages):

```{r, fig.width=4, fig.asp=.8}
draw.ccurve(1600, 2020, BCAD=TRUE, cc2='nh1')
```

Since the postbomb curve dwarfs the IntCal20 curve, we could also plot both on separate vertical axes:

```{r, fig.width=4, fig.asp=.8}
draw.ccurve(1600, 2020, BCAD=TRUE, cc2='nh1', add.yaxis=TRUE)
```

The calibration curves can also be plotted in the 'realms' of F14C, pMC or D14C (see below), e.g.:

```{r, fig.width=4, fig.asp=.8}
draw.ccurve(50000, 35000, realm="D")
```


# Realms

This package provides functions to translate values between the radiocarbon-relevant 'realms' of cal BP (calendar years before AD 1950), BC/AD (here we mean cal BC/AD, but shorten this to BC/AD), C14, F14C, pMC and D14C. The following Table lists the available functions:

| from\\to  | calBP         | BCAD          | C14          | F14C         | pMC         | D14C         |
|-----------|---------------|---------------|--------------|--------------|-------------|--------------|
| **calBP** | -             | `calBPtoBCAD` | `calBPtoC14` | `calBPtoF14C`| `calBPtopMC`| `calBPtoD14C`|
| **BCAD**  | `BCADtocalBP` | -             | `BCADtoC14`  | `BCADtoF14C` | `BCADtopMC` | `BCADtoD14C` |
| **C14**   | `C14tocalBP`  | `C14toBCAD`   | -            | `C14toF14C`  | `C14topMC`  | `C14toD14C`  |
| **F14C**  | NA            | NA            | `F14CtoC14`  | -            | `F14CtopMC` | `F14CtoD14C` |
| **pMC**   | NA            | NA            | `pMCtoC14`   | `pMCtoF14C`  |  -          | `pMCtoD14C`  |
| **D14C**  | NA            | NA            | `D14CtoC14`  | `D14CtoF14C` | `D14CtopMC` |  -           |

As an example of the above functions, the IntCal20 C14 age and error belonging to one or more cal BP or cal BCAD ages can be found (interpolating linearly where necessary):  

```{r}
calBPtoC14(10.5)
BCADtoC14(1940:1950)
```

To translate between cal BC/AD and cal BP ages, we can use (the last example avoids 0 BC/AD, since some calendars do not include zero):

```{r}
BCADtocalBP(2024)
BCADtocalBP(-1, zero=TRUE)
BCADtocalBP(-1, zero=FALSE)
```

D14C ($\Delta$^14^C, a proxy for atmospheric ^14^C concentration at *t* cal BP) can be transferred to F^14^C, and the other way around:

```{r}
D14CtoF14C(152, t=4000)
F14CtoD14C(0.71, t=4000)
```

This can also be done with C14 ages:

```{r}
C14toD14C(152, t=4000)
D14CtoC14(592, t=4000)
```

These functions can be used to investigate $\Delta^{14}C$ over time:

```{r, fig.width=4, fig.asp=.8}
x <- seq(0, 55e3, length=1e3)
cc <- calBPtoC14(x)
Dcc <- C14toD14C(cc[,1], cc[,2], x)

par(mar=c(4,3,1,3), bty="l")
plot(x/1e3, Dcc[,1]+Dcc[,2], type="l", xlab="kcal BP", ylab="")
mtext(expression(paste(Delta, ""^{14}, "C")), 2, 1.7)
lines(x/1e3, Dcc[,1]-Dcc[,2])

par(new=TRUE)
plot(x/1e3, (cc[,1]-cc[,2])/1e3, type="l", xaxt="n", yaxt="n", col=4, xlab="", ylab="")
lines(x/1e3, (cc[,1]+cc[,2])/1e3, col=4)
mtext(expression(paste(""^{14}, "C kBP")), 4, 2, col=4)
axis(4, col=4, col.axis=4, col.ticks=4)
```

# Pooling dates

Sometimes, the same material is measured using multiple radiocarbon dates. If we can be sure that the dated material stems from one single age in time (e.g., multiple dates on the same single bone, or perhaps the cereal grains within one bowl which could be assumed to all stem from the same season), then we can check to which degree the dates agree using a Chi^2^-test (Ward & Wilson 1978)[^5]. If they do agree, then a pooled mean can be calculated. For example, take the Shroud of Turin, which was dated multiple times in three different labs:

```{r}
data(shroud)
shroud
pool(shroud$y,shroud$er) 
Zu <- grep("ETH", shroud$ID) # Zurich lab only
pool(shroud$y[Zu],shroud$er[Zu])
```

Then the weighted mean can be calculated, with the error taken as the largest of the weighted uncertainty and the standard deviation:

```{r}
weighted_means(shroud$y[Zu],shroud$er[Zu])
```

If it can indeed be safely assumed that all dates stem from the same (unknown) calendar year, then the age distribution of that single year can be plotted together with the individual calibrated ages:

```{r, fig.width=5, fig.asp=.8}
as.one(shroud$y,shroud$er)
```
 
It would however often be much safer to assume that the multiple dates were deposited over not just one calendar year but rather over a period, e.g., over 50 years. To do so, a moving bin is made, and for each bin placement it is checked how much of the calibrated distribution of each date fits within that bin. Here is an example, using a bin width of 100 years, moving at 25 year steps (the top value indicates that around 660 cal BP, a total of around 4-5 dates fit within the 100-year bin):

```{r, fig.width=5, fig.asp=.8}
as.bin(shroud$y,shroud$er, 100, move.by=25)
```

The spread of multiple calibrated dates can be visualised and summarised:

```{r, fig.width=5, fig.asp=.8}
spread(shroud$y,shroud$er)
```

We can assess the overlap between the calibrated distributions of multiple dates (as the minimum calibrated height for each calendar year in a range). For example, imagine two dates from the same archaeological layer, a twig (3820 &pm; 40 C-14 BP, requiring calibration with IntCal20) and a marine shell (4430 &pm; 40 C-14 BP, requiring Marine20, with a regional marine offset of 90 &pm; 25). To what degree do their calibrated distributions overlap?

```{r, fig.width=5, fig.asp=.8}
y <- c(3820, 4590-90) 
er <- c(40, sqrt(40^2 + 25^2)) 
cc <- c(1,2)
overlap(y, er, cc=cc)
```

# Contamination

To calculate the effect of contamination on radiocarbon ages, e.g. what age would be observed if material with a "true" radiocarbon age of 5000 &pm; 20 ^14^C BP would be contaminated with 10\% (&pm; 0\%) of modern carbon (F^14^C=1)?

```{r, fig.width=5, fig.asp=.8}
contaminate(5000, 20, 10, 0, 1)
```

Or imagine that you measured a dinosaur bone, dating to far beyond the limit of radiocarbon dating, and the sample is very clean as it contains only 0.5\% &pm; 0.1\% modern contamination:

```{r, fig.width=5, fig.asp=.8}
contaminate(66e6, 1e6, 0.5, 0.1, 1)
```

The other way round, e.g., inferring what would happen to an observed age if its assumed 10% modern contamination were to be removed:

```{r, fig.width=5, fig.asp=.8}
clean(9000, 100, 10)
```

We can also calculate the amount of contamination, or muck, required to 'explain away' certain ages. For example, one of the measurements of the Shroud of Turin dates to 591 &pm; 30 C14 BP. How much modern contamination would have to be inferred for the material to really date to, say, AD 40?

```{r, fig.width=5, fig.asp=.8}
muck(591, 30, BCADtoC14(40)[,1], 0, 1)
```

So we'd require the sample to have been contaminated by 67% modern carbon to still date to around AD 40. But what if the sample had been repaired in, say, AD 1400, thus adding material of an not-entirely-modern F14C value (i.e., taking into account both the atmospheric C-14 concentrations in AD 1400 and the fact that some of the C14 will have decayed since then)?

```{r, fig.width=5, fig.asp=.8}
perFaith <- BCADtoC14(40)
repairF <- BCADtoF14C(1400)
muck(591, 30, perFaith[,1], perFaith[,2], repairF[,1], repairF[,2])
```

This means that the dated sample would have to consist almost entirely of Medieval age material - which is exactly what was found by Damon et al. (1989)[^6].

The effect of different levels of contamination can be visualised:

```{r, fig.width=6, fig.asp=.8}
real.14C <- seq(0, 50e3, length=200)
contam <- seq(0, 10, length=101) # 0 to 10% contamination
contam.col <- rainbow(length(contam))
plot(0, type="n", xlim=c(0, 55e3), xlab="real 14C age", ylim=range(real.14C), ylab="observed 14C age")
for (i in 1:length(contam)) {
  observed <- contaminate(real.14C, 0, contam[i], 0, 1, decimals=5, talk=FALSE, MC=FALSE)
  lines(real.14C, observed[,1], col = contam.col[i])
}
```

If that is too much code for you, try this function instead:

```{r, fig.width=6, fig.asp=.8}
draw.contamination()
```

# Fractions

Sometimes, one needs to estimate a missing radiocarbon age from a sample which has C14 dates on both the entire sample and on fractions, but where one of the samples was too small to be dated. This can be used in for example soils separated into size fractions, or samples dated using both bulk and humic/humin fractions, where one of the samples turns out to be too small to be dated. This equation requires the bulk age, the ages of the dated fractions, and the carbon contents and weights of all fractions.

```{r}
Cs <- c(.02, .05, .03, .04) # carbon contents of each fraction
wghts <- c(5, 4, 2, .5) # weights for all fractions, e.g., in mg
ages <- c(130, 130, 130, NA) # ages of all fractions. The unmeasured one is NA
errors <- c(10, 12, 10, NA) # errors, unmeasured is NA
fractions(150, 20, Cs, wghts, ages, errors) # assuming a bulk age of 150 +- 20 C14 BP 
```

# Calibration

Now on to calibration of radiocarbon dates. We can obtain the calibrated probability distributions from radiocarbon dates, e.g., one of 130 &pm; 10 C14 BP:

```{r, fig.width=4, fig.asp=.8}
calib.130 <- caldist(130, 10, BCAD=TRUE)
plot(calib.130, type="l")
```

It is also possible to find the likelihood of a single calendar year for our radiocarbon age, e.g., 145 cal BP:

```{r}
l.calib(145, 130, 10)
```

To sample n=10 random values from a calibrated distribution (where the probability of a year being sampled is proportional to its calibrated height):

```{r, fig.width=4, fig.asp=.8}
dice <- r.calib(100, 130, 10)
plot(density(dice))
rug(dice)
```

For reporting purposes, calibrated dates are often reduced to their 95% highest posterior density (hpd) ranges (please report all, not just your favourite one!):

```{r}
hpd(calib.130)
```

Additionally, calibrated dates are often reduced to single point estimates. Note however how poor representations they are of the entire calibrated distribution!

```{r, fig.width=4, fig.asp=.8}
calib.2450 <- caldist(2450, 20)
plot(calib.2450, type="l")
points.2450 <- point.estimates(calib.2450)
points.2450
abline(v=points.2450, col=1:4, lty=2)
```

Want a plot of the radiocarbon and calibrated distributions, together with their hpd ranges?
```{r, fig.width=5, fig.asp=1}
calibrate(2450, 40)
```

Sometimes one would want to smooth a calibration curve to take into account the fact that material has accumulated over a certain time (e.g., a tree, or peat). To do so, a calibration curve can be smoothed to produce a tailor-made calibration curve, after which this one is used to calibrate the date:

```{r, fig.width=5, fig.asp=1}
mycurve <- smooth.ccurve(smooth=50)
calibrate(2450, 40, thiscurve=mycurve)
```

Calibrating 'young' radiocarbon dates (close to 0 C14 BP) can cause an error, because a bomb curve might be required to capture the youngest ages. Do not worry, there is an option to avoid that error:

```{r, fig.width=5, fig.asp=1}
try(calibrate(130,30))
calibrate(130, 30, bombalert=FALSE)
```

It is also possible to analyse the calibrated probability distributions, e.g. what is the probability (between 0 and 1) that the date stems from material that is of the age of 150 cal BP or younger? Or that it is older than that age?

```{r}
younger(150, 130, 10)
older(150, 130, 10)
```

Or if you wish to calculate the probability covered between two age points (e.g., here between 300 and 150 cal BP):

```{r}
p.range(300, 150, 130, 10)
```

The probabilities reported by `p.range` should be similar to those reported by the `hpd` function mentioned earlier. Any discrepancies can be explained by rounding errors in the calculations (e.g., hpd's settings of `every`, `age.round` and `prob.round`).

At times there could be a lag in the date, for example if material has accumulated over decades. We can calibrate a date and then add a lag, by adding or subtracting either a normal distribution or a gamma one:

```{r, fig.width=5, fig.asp=1}
push.normal(2450,40, 400,20)
```

# Marine offsets

Dates on marine material will often have to be calibrated with the Marine20 calibration curve[^7], and many coastal locations will have an additional regional reservoir offset (deltaR) (Reimer and Reimer 2006)[^8]. The on-line database at <http://calib.org/marine/> is very useful for this; it features the radiocarbon ages and deltaR of many shells of known collection date. The data from this database were downloaded (in August 2024) and can be queried. For example, a map can be drawn with all shell data within certain coordinates:

```{r, fig.width=5, fig.asp=1}
myshells <- map.shells(S=54, W=-8, N=61, E=0) # the northern part of the UK
```

If you don't have the R package 'rnaturalearth' installed, you will be suggested to install it if possible, because the resulting maps will look much nicer. 

The output can also be queried:

```{r, fig.width=5, fig.asp=1}
head(myshells)
shells.mean(myshells)
```

You can also extract say the 20 shells closest to a coordinate, e.g., 120 East and 10 North:
```{r, fig.width=5, fig.asp=1}
myshells <- find.shells(120, 10, 20)
shells.mean(myshells, distance=TRUE)
```

# Multiple calibrations

It is also possible to draw multiple calibrated distributions:

```{r, fig.width=4, fig.asp=1}
set.seed(123)
dates <- sort(sample(500:2500,5))
errors <- .05*dates
depths <- 1:length(dates)
my.labels <- c("my", "very", "own", "simulated", "dates")
draw.dates(dates, errors, depths, BCAD=TRUE, labels=my.labels, age.lim=c(0, 1800))
```

or add them to an existing plot:

```{r, fig.width=4, fig.asp=1}
plot(250*1:5, 5:1, xlim=c(0, 1800), ylim=c(5,0), xlab="AD", ylab="dates")
draw.dates(dates, errors, depths, BCAD=TRUE, add=TRUE, labels=my.labels, mirror=FALSE)
```

or get creative (inspired by Jocelyn Bell Burnell[^9], Joy Division[^10] and the Hallstatt Plateau[^11]):

```{r, fig.width=4, fig.asp=1}
par(bg="black", mar=rep(1, 4))
n <- 50; set.seed(1)
draw.dates(rnorm(n, 2450, 30), rep(25, n), 1:n,
  mirror=FALSE, draw.base=FALSE, draw.hpd=FALSE, col="white",
  threshold=1e-28, age.lim=c(2250, 2800), up=TRUE, ex=-20)
```


[^1]: Stuiver, R., Polach, H.A., 1977. Discussion: reporting of 14C data.
*Radiocarbon* 19, 355-363 <http://dx.doi.org/10.1017/S0033822200003672>

[^2]: Reimer, P.J., Brown, T.A., Reimer, R.W., 2004. Discussion: reporting
and calibration of post-bomb 14C Data. *Radiocarbon* 46, 1299-1304 <http://dx.doi.org/10.1017/S0033822200033154>

[^3]: Millard, R., 2014. Conventions for reporting radiocarbon determinations. *Radiocarbon* 56, 555-559 <http://dx.doi.org/10.2458/56.17455>

[^4]: Reimer, P.J., et al., 2020. The IntCal20 Northern Hemisphere radiocarbon age calibration curve (0-55 cal kBP). *Radiocarbon* 62, 725-757 <http://dx.doi.org/10.1017/S0033822200032999>

[^5]: Ward, G.K., Wilson, S., 1978. Procedures for comparing and combining radiocarbon age determinations: A critique. *Archaeometry* 20, 19-31 <http://dx.doi.org/10.1111/j.1475-4754.1978.tb00208.x>

[^6]: Damon, P., et al., 1989. Radiocarbon dating of the Shroud of Turin. *Nature* 337, 611–615. <https://doi.org/10.1038/337611a0>

[^7]: Heaton, T.J., et al., 2020. Marine20-the marine radiocarbon age calibration curve (0-55,000 cal BP). *Radiocarbon* 62, 779-820 <http://dx.doi.org/10.1017/RDC.2020.68>

[^8]: Reimer, P.J., Reimer, R.W., 2006. A marine reservoir correction database and on-line interface. *Radiocarbon* 43, 461-463. <http://dx.doi.org/10.1017/S0033822200038339>

[^9]: https://www.cam.ac.uk/stories/journeysofdiscovery-pulsars

[^10]: https://www.radiox.co.uk/artists/joy-division/cover-joy-division-unknown-pleasures-meaning/

[^11]: https://en.wikipedia.org/wiki/Hallstatt_plateau
