R/classes.R

In pvermees/IsoplotR: Statistical Toolbox for Radiometric Geochronology

#' @title
#' Geochronological data classes
#'
#' @description S3 classes to store geochronological data generated by
#'     \code{\link{read.data}} or \code{\link{diseq}}.
#'
#' @details \code{IsoplotR} uses the following S3 classes to store
#'     geochronological data: \code{UPb}, \code{PbPb}, \code{ThPb},
#'     \code{KCa}, \code{UThHe}, \code{fissiontracks},
#'     \code{detritals} and \code{PD}, where the latter is the parent
#'     class of the simple parent-daughter chronometers \code{RbSr},
#'     \code{SmNd}, \code{LuHf} and \code{ReOs}. All these classes
#'     have overloaded versions of the generic \code{length()}
#'     function and \code{`[`} subsetting method.
#'
#'     Additional functions for each class include \code{as.X(x)},
#'     which converts the data table \code{x} to an object of class
#'     \code{X}; and \code{is.X(x)}, which checks if \code{x} has
#'     class \code{X}.
#'
#' \itemize{
#' 
#' \item\code{UPb}: a list containing:
#'
#' \describe{
#'
#' \item{\code{x}}{a matrix containing the isotopic measurements}
#'
#' \item{\code{format}}{a number between 1 and 8}
#'
#' \item{\code{d}}{an object of class \code{diseq}, i.e. the output of
#' the \code{\link{diseq}} function}
#'
#' }
#'
#' \item\code{ArAr}: a list containing
#'
#' \describe{
#'
#' \item{\code{x}}{a matrix containing the isotopic measurements}
#'
#' \item{\code{J}}{a two-element vector with the J-factor and its
#' uncertainty}
#' 
#' \item{\code{format}}{a number between 1 and 3}
#'
#' }
#'
#' \item\code{ThU}: a list containing
#'
#' \describe{
#'
#' \item{\code{x}}{a matrix containing the isotopic measurements}
#'
#' \item{\code{format}}{a number between 1 and 4}
#'
#' \item{\code{Th02}}{a two element vector with the assumed initial
#'     \eqn{^{230}}Th/\eqn{^{232}}Th-ratio of Th-bearing
#'     detritus. Only aplicable to formats 1 and 2.}
#' 
#' \item{\code{Th02U48}}{9-element vector with the measured
#'     composition of Th-bearing detritus}
#'
#' \item{\code{U8Th2}}{the measured \eqn{^{238}}U/\eqn{^{232}}Th
#'       activity ratio of the whole rock. Only applicable to formats
#'       3 and 4}
#'
#' }
#' 
#' \item\code{PbPb}, \code{ThPb}, \code{KCa}, \code{PD}, \code{RbSr},
#' \code{SmNd}, \code{LuHf}, or \code{ReOs}: a list containing
#'
#' \describe{
#'
#' \item{\code{x}}{a matrix containing the isotopic measurements}
#'
#' \item{\code{format}}{a number between 1 and 3}
#'
#' }
#'
#' \item\code{UThHe}: a matrix of He, U, Th (and Sm) measurements
#'
#' \item\code{fissiontracks}: a list containing
#'
#' \describe{
#'
#' \item{\code{format}}{a number between 1 and 3}
#'
#' \item{\code{x}}{a matrix of spontaneous and induced fission track
#' counts (only included if \code{format=1})}
#'
#' \item{\code{rhoD}}{the track density of the dosimeter glass,
#' extracted from the input data (only included if \code{format=1})}
#'
#' \item{\code{zeta}}{the zeta calibration constant extracted from the
#' input data (only included if \code{format<3})}
#'
#' \item{\code{Ns}}{a list containing the spontaneous fission track
#' counts (only included if \code{format>1})}
#'
#' \item{\code{U}}{a list of lists containing the U-concentration or
#' U/Ca-ratio measurements for each of the analysed grains (only
#' included if \code{format>1})}
#'
#' \item{\code{sU}}{a list of lists containing the standard errors of
#' the U-concentration or U/Ca-ratio measurements for each of the
#' analysed grains (only include if \code{format>1})}
#'
#' \item{\code{spotSize}}{the laser ablation spot size (only included
#' if \code{format>1})}
#'
#' }
#'
#' \item\code{detritals}: a list of named vectors, one for each
#' detrital sample.
#' 
#' \item\code{diseq}: is a class that contains the output of the
#' \code{\link{diseq}} function, which stores initial disequilibrium
#' data for U--Pb geochronology.
#'
#' }
#'
#' @param x a data object returned by \code{\link{read.data}} or
#'     \code{\link{diseq}}.
#' 
#' @param format data format. See \code{\link{read.data}} for details.
#' 
#' @param ierr input error. See \code{\link{read.data}} for details.
#'
#' @param d an object of class \code{\link{diseq}}.
#' 
#' @param Th02i 2-element vector with the assumed initial
#'     \eqn{^{230}}Th/\eqn{^{232}}Th-ratio of the detritus (for
#'     Th-U formats 1 and 2) and its standard error.
#' 
#' @param Th02U48 9-element vector with the measured composition of
#'     the detritus, containing \code{X=0/8}, \code{sX}, \code{Y=2/8},
#'     \code{sY}, \code{Z=4/8}, \code{sZ}, \code{rXY}, \code{rXZ},
#'     \code{rYZ}.
#' 
#' @param U8Th2 \eqn{^{238}}U/\eqn{^{232}}Th activity-ratio of the
#'     whole rock. Used to estimate the initial
#'     \eqn{^{230}}Th/\eqn{^{238}}U disequilibrium (for Th-U formats 3
#'     and 4).
#' 
#' @return \code{is.X(x)} returns a logical value.
#'
#' \code{as.X(x)} returns an object of class \code{X}.
#' 
#' @examples
#' attach(examples)
#' ns <- length(UPb)
#' concordia(UPb[-ns,])
#' if (is.PD(RbSr)) print('RbSr has class PD')
#' 
#' @seealso read.data diseq
#' @name classes
#' @aliases UPb PbPb ThPb ThU ArAr KCa PD RbSr SmNd LuHf ReOs UThHe
#'     fissiontracks detritals is.UPb is.PbPb is.ThPb is.ThU is.ArAr
#'     is.KCa is.PD is.RbSr is.SmNd is.LuHf is.ReOs is.UThHe
#'     is.fissiontracks is.detritals is.other as.UPb as.PbPb as.ThPb
#'     as.ThU as.ArAr as.KCa as.PD as.RbSr as.SmNd as.LuHf as.ReOs
#'     as.UThHe as.fissiontracks as.detritals as.other
NULL

#' @rdname classes
#' @export
is.UPb <- function(x) inherits(x,"UPb")
#' @rdname classes
#' @export
is.PbPb <- function(x) inherits(x,"PbPb")
#' @rdname classes
#' @export
is.ThPb <- function(x) inherits(x,"ThPb")
#' @rdname classes
#' @export
is.ArAr <- function(x) inherits(x,"ArAr")
#' @rdname classes
#' @export
is.KCa <- function(x) inherits(x,"KCa")
#' @rdname classes
#' @export
is.PD <- function(x) inherits(x,"PD")
#' @rdname classes
#' @export
is.RbSr <- function(x) inherits(x,"RbSr")
#' @rdname classes
#' @export
is.SmNd <- function(x) inherits(x,"SmNd")
#' @export
#' @rdname classes
is.LuHf <- function(x) inherits(x,"LuHf")
#' @rdname classes
#' @export
is.ReOs <- function(x) inherits(x,"ReOs")
#' @rdname classes
#' @export
is.ThU <- function(x) inherits(x,"ThU")
#' @rdname classes
#' @export
is.UThHe <- function(x) inherits(x,"UThHe")
#' @rdname classes
#' @export
is.fissiontracks <- function(x) inherits(x,"fissiontracks")
#' @rdname classes
#' @export
is.detritals <- function(x) inherits(x,"detritals")
#' @rdname classes
#' @export
is.other <- function(x) inherits(x,"other")
#' @rdname classes
#' @export
is.diseq <- function(x) inherits(x,"diseq")

#' @export
length.UPb  <- function(x){ nrow(x$x) }
#' @export
length.PbPb <- function(x){ nrow(x$x) }
#' @export
length.ArAr <- function(x){ nrow(x$x) }
#' @export
length.PD <- function(x){ nrow(x$x) }
#' @export
length.ThPb <- function(x){ nrow(x$x) }
#' @export
length.KCa <- function(x){ nrow(x$x) }
#' @export
length.ThU <- function(x){ nrow(x$x) }
#' @export
length.UThHe <- function(x){ nrow(x) }
#' @export
length.KDE <- function(x){ length(x$ages) }
#' @export
length.KDEs <- function(x){ length(x$kdes) }
#' @export
length.fissiontracks <- function(x){
    if (x$format==1) return(nrow(x$x))
    else return(length(x$Ns))
}
#' @export
length.other <- function(x){
    if (x$format==1) return(length(x$x))
    else if (x$format==6) return(nrow(x$x)/2)
    else return(nrow(x$x))
}

#' @export
`[.UPb` <- function(x,...){
    out <- x
    out$x <- x$x[...]
    if ('x.raw' %in% names(x)){
        out$x.raw <- x$x.raw[...]
    }
    out
}
#' @export
`[.ArAr` <- function(x,...){ `[helper`(x,...) }
#' @export
`[.PbPb` <- function(x,...){ `[helper`(x,...) }
#' @export
`[.ThPb` <- function(x,...){ `[helper`(x,...) }
#' @export
`[.KCa` <- function(x,...){ `[helper`(x,...) }
#' @export
`[.PD` <- function(x,...){ `[helper`(x,...) }
#' @export
`[.ThU` <- function(x,...){ `[helper`(x,...) }
#' @export
`[.fissiontracks` <- function(x,...){
    if (x$format==1){
        out <- `[helper`(x,...)
    } else {
        out <- x
        out$Ns <- x$Ns[...]
        out$A <- x$A[...]
        out$U <- x$U[...]
        out$sU <- x$sU[...]
    }
    out
}
#' @export
`[.diseq` <- function(x,i){
    out <- x
    for (ratio in c('U48','ThU','RaU','PaU')){
        j <- min(length(out[[ratio]]$x),i)
        out[[ratio]]$x <- out[[ratio]]$x[j]
        k <- min(length(out[[ratio]]$sx),i)
        out[[ratio]]$sx <- out[[ratio]]$sx[k]
    }
    out
}
`[helper` <- function(x,...){
    out <- x
    out$x <- x$x[...]
    out
}

#' @export
subset.UPb  <- function(x,...){
    out <- x
    out$x <- subset(x$x,...)
    if ('x.raw' %in% names(x)){
        out$x.raw <- subset(x$x.raw,...)
    }
    out
}
#' @export
subset.PbPb <- function(x,...){ subset_helper(x,...) }
#' @export
subset.ArAr <- function(x,...){ subset_helper(x,...) }
#' @export
subset.ThPb <- function(x,...){ subset_helper(x,...) }
#' @export
subset.KCa <- function(x,...){ subset_helper(x,...) }
#' @export
subset.PD <- function(x,...){ subset_helper(x,...) }
#' @export
subset.ThU <- function(x,...){ subset_helper(x,...) }
#' @export
subset.UThHe <- function(x,...){
    out <- subset.matrix(x,...)
    class(out) <- class(x)
    out
}
#' @export
subset.fissiontracks <- function(x,...){
    if (x$format==1){
        out <- subset_helper(x,...)
    } else {
        out <- x
        out$Ns <- subset(x$Ns,...)
        out$A <- subset(x$A,...)
        out$U <- subset(x$U,...)
        out$sU <- subset(x$sU,...)
    }
    out
}
#' @export
subset.detritals <- function(x,subset,...){
    out <- x[subset]
    class(out) <- class(x)
    out
}
#' @export
subset.other <- function(x,subset,...){
    out <- x
    if (x$format==1){
        out$x <- x$x[,subset]
    } else if (x$format==6){
        out$x <- subset_ogls(x$x,subset)
    } else {
        out$x <- subset.matrix(out$x,subset=subset,...)
    }
    out
}
subset_helper <- function(x,...){
    out <- x
    out$x <- subset.matrix(x$x,...)
    out
}
subset_ogls <- function(x,subset){
    ns <- nrow(x)/2
    iX <- (1:ns)[subset]
    iY <- ((ns+1):(2*ns))[subset]
    i <- c(iX,iY)
    x[i,c(1,i+1)]
}

mediand <- function(d){
    out <- d
    for (ratio in c('U48','ThU','RaU','PaU')){
        out[[ratio]]$x <- stats::median(d[[ratio]]$x)
        out[[ratio]]$sx <- stats::median(d[[ratio]]$sx)
    }
    out    
}