---
title: rice
output:
  html_vignette:
    toc: true
    toc_depth: 3
    number_sections: true
vignette: >
  %\VignetteIndexEntry{rice}
  %\VignetteEngine{knitr::rmarkdown}
  \usepackage[utf8]{inputenc}
---

# 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

## Plots

Calibration curves can be plotted:

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

Or, comparing two calibration curves:

```{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 (using the BCAD scale):

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

Interesting things happened after 1950, as can be seen by adding a postbomb curve:

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

The postbomb curve dwarfs the IntCal20 curve, so 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)
```

To find the IntCal20 radiocarbon year belonging to a certain calendar year, e.g., 900 cal BP:

```{r}
calBP.14C(900)
```

# Calculations

This package also provides functions related to radiocarbon 'realms'. First there are two functions to calculate radiocarbon ages from pMC values (in this case of a postbomb date):

```{r}
pMC.age(150, 1)
```

and the other way round:

```{r}
age.pMC(-2300, 40)
```

The same for calculations in the F^14^C realm:

```{r}
F14C.age(.150, .01)
```

and the other way round:

```{r}
age.F14C(-2300, 40)
```

To transfer $\Delta^{14}C$ (a proxy for atmospheric ^14^C concentration at *t* cal BP) to F^14^C, and the other way around:

```{r}
F14C.D14C(0.71, t=4000)
D14C.F14C(152, 4000)
```

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

```{r, fig.width=4, fig.asp=.8}
cc <- rintcal::ccurve()
cc.Fmin <- age.F14C(cc[,2]+cc[,3])
cc.Fmax <- age.F14C(cc[,2]-cc[,3])
cc.D14Cmin <- F14C.D14C(cc.Fmin, cc[,1])
cc.D14Cmax <- F14C.D14C(cc.Fmax, cc[,1])
par(mar=c(4,3,1,3), bty="l")
plot(cc[,1]/1e3, cc.D14Cmax, type="l", xlab="kcal BP", ylab="")
mtext(expression(paste(Delta, ""^{14}, "C")), 2, 1.7)
lines(cc[,1]/1e3, cc.D14Cmin)
par(new=TRUE)
plot(cc[,1]/1e3, (cc[,2]+cc[,3])/1e3, type="l", xaxt="n", yaxt="n", col=4, xlab="", ylab="")
lines(cc[,1]/1e3, (cc[,2]-cc[,3])/1e3, col=4)
axis(4, col=4, col.axis=4)
mtext(expression(paste(""^{14}, "C kBP")), 4, 2, col=4)
```

## Contamination

The above functions can be used 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 +- 20 ^14^C BP would be contaminated with 1% of modern carbon (F^14^C=1)?

```{r}
contaminate(5000, 20, .01, 1)
```

The effect of different levels of contamination can also be visualised:

```{r, fig.width=6, fig.asp=.8}
real.14C <- seq(0, 50e3, length=200)
contam <- seq(0, .1, 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))
  lines(real.14C, contaminate(real.14C, c(), contam[i], 1, decimals=5), col=contam.col[i])
contam.legend <- seq(0, .1, length=6)
contam.col <- rainbow(length(contam.legend)-1)
text(50e3, contaminate(50e3, c(), contam.legend, 1), 
  labels=contam.legend, col=contam.col, cex=.7, offset=0, adj=c(0,.8))
```

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

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

## 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)
```

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(130,10)
```

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)
```


## Marine offsets

Dates on marine material will often have to be calibrated with the Marine20 calibration curve[^5], and many coastal locations will have an additional regional reservoir offset (deltaR) (Reimer and Reimer 2006)[^6]. 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
```

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

You can also draw one or more 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(300*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[^7], Joy Division[^8] and the Hallstatt Plateau[^9]):

```{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), n:1,
  mirror=FALSE, draw.base=FALSE, draw.hpd=FALSE, col="white",
  threshold=1e-28, age.lim=c(2250, 2800), ex=.8)
```

[^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]: 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>

[^6]: 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>

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

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

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