diff --git a/DESCRIPTION b/DESCRIPTION index 16e9244..c8f258c 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,16 +1,16 @@ Package: deform Type: Package Title: Crustal deformation modeling tools in R -Version: 0.0-7.9999 -Date: 2015-08-05 +Version: 0.0-8.9999 +Date: 2015-08-19 Author: Andrew J. Barbour; some of these tools were adapted from Matlab(tm) codes by Brendan Crowell (UW), and Francois Beauducel (IPGP). Maintainer: A.J.Barbour Description: Crustal deformation modeling tools. Currently this code supports calculation of surface deformation associated with - (1) a depleting layer using a Segall-type model, and - (2) a volume source using a Mogi-type model. + (1) a depleting reservoir using Segall's models, and + (2) a spherical volume source using a Mogi's model. License: GPL Depends: rgdal (>= 0.9-2) @@ -23,5 +23,6 @@ Suggests: TeachingDemos, maps, sp, - plyr + plyr, + pracma LinkingTo: Rcpp diff --git a/NAMESPACE b/NAMESPACE index f985e08..62bb3f0 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -2,6 +2,7 @@ export(.setleft) export(.surface_g) +export(.surface_g_pt) export(.surface_g_xi) export(D_project) export(Tilt) @@ -25,10 +26,12 @@ export(segall85) export(segall89) export(static) export(surface_displacement) +export(surface_displacement_point) export(surface_displacement_reservoir) export(timevarying_fluidmass) export(timevarying_porepressure) export(timevarying_surface_displacement) +export(vasco88) exportPattern("^[[:alpha:]]+") importFrom(RCurl,getURL) importFrom(RCurl,url.exists) diff --git a/R/segall85.R b/R/segall85.R index 9708c30..014c298 100644 --- a/R/segall85.R +++ b/R/segall85.R @@ -115,13 +115,13 @@ segall85 <- function(help=FALSE){ #' @param x numeric; spatial coordinate relative to extraction point #' @param C. numeric; units of force, proportional to source strength (e.g., extraction rate) #' @param mu. numeric; the shear modulus in Pascals -#' @param ... additional arguments to \code{\link{.surface_g}} +#' @param ... additional arguments passed to \code{\link{.surface_g}} #' @rdname segall85 #' @export surface_displacement <- function(x, C.=1, mu.=1e9, ...){ sg <- .surface_g(x, ...) c. <- C./mu. - sg <- mutate(sg, + sg <- plyr::mutate(sg, ux = gx * c., uz = gz * c., uxz.mag = sqrt(ux^2 + uz^2), diff --git a/R/vasco88.R b/R/vasco88.R new file mode 100644 index 0000000..16df53c --- /dev/null +++ b/R/vasco88.R @@ -0,0 +1,64 @@ +#' Surface deformation associated with fluid withdrawl: an alternative formulation +#' @name vasco88 +#' @aliases vasco1988 +#' @export +#' @param help logical; load documentation for \code{\link{vasco88}} +#' @seealso \code{\link{segall85}}, or \code{\link{segall89}}; \code{\link{Simple-deformation}} +#' @export +#' @examples +#' # Reproduce... +#' r <- seq(-10,10,by=0.2) # km +#' r.m <- r*1e3 +#' xx <- surface_displacement_point(r.m, depth=2e3, delV.=1e7) +#' plot(uz ~ x, xx, type='l') +vasco88 <- function(){ + cat("\nThis function is simply a placeholder. See the documentation ( `?vasco88` ).\n") + if (help) ?vasco88 +} + +#' @rdname vasco88 +#' @export +.surface_g_pt <- function(x=0, x_src=0, z_src=0, nuu=1/3){ + # vasco88 11 + sc <- (1 + nuu)/(3*pi) + xrel <- (x - x_src) + g <- sc * xrel / (xrel^2 + z_src^2)**(3/2) + data.frame(x, g, xz = x/z_src) +} + +#' @rdname vasco88 +#' @export +#' @inheritParams surface_displacement_reservoir +#' @param depth numeric; the depth below the surface to the source +#' @param delV. numeric; uniform change in pore-fluid volume (positive = loss) +#' @param tol numeric; the numerical tolerance in the integration; if +#' supplied, this should be much smaller than the smallest difference between +#' any values in \code{x} divided by the depth +#' @param ... additional arguments passed to \code{\link{.surface_g_pt}} +surface_displacement_point <- function(x, depth, delV., B.=1, rho_f.=1000, tol, ...){ + x <- as.vector(x) + depth <- depth[1] + if (depth <= 0) stop('source must be below the surface (positive depth)') + delV. <- delV.[1] + # Setup some vectors/info for interpolation + if (missing(tol)) tol <- 0.01 + xr <- range(pretty(x)) + mx <- ceiling(max(abs(x), na.rm = TRUE)/depth) + mxsc <- 15 # this many times the x/z distance + xcalc <- depth * seq.int(from = -1* mxsc * mx, to = mxsc * mx, by = tol) + # Calculate point source response + sgcalc <- .surface_g_pt(xcalc, z_src = depth, ...) + # Calculate response function + ix <- sgcalc$x + iy <- sgcalc$g + n <- length(ix) + ih <- ix[2] - ix[1] + ca <- (iy[2] - iy[1]) / ih + cb <- (iy[n] - iy[n - 1]) / ih + Gfun <- approxfun(ix, pracma::cumtrapz(ix, iy) - (cb - ca) * ih^2/12) + # Use interpolation function to resample integrated response function + g <- Gfun(x) + C. <- B. * delV. / rho_f. + sg <- data.frame(x, g, xz = x/depth, uz = g * C.) + return(sg) +} diff --git a/man/segall85.Rd b/man/segall85.Rd index 38f28b5..246e9a8 100644 --- a/man/segall85.Rd +++ b/man/segall85.Rd @@ -34,7 +34,7 @@ timevarying_porepressure(x, z, Time, Vdot., B., L., D., HD., t., mu.gpa., \item{mu.}{numeric; the shear modulus in Pascals} -\item{...}{additional arguments to \code{\link{.surface_g}}} +\item{...}{additional arguments passed to \code{\link{.surface_g}}} \item{x_src}{numeric; the horizontal distance from the source} diff --git a/man/vasco88.Rd b/man/vasco88.Rd new file mode 100644 index 0000000..c21664e --- /dev/null +++ b/man/vasco88.Rd @@ -0,0 +1,48 @@ +% Generated by roxygen2 (4.1.1): do not edit by hand +% Please edit documentation in R/vasco88.R +\name{vasco88} +\alias{.surface_g_pt} +\alias{surface_displacement_point} +\alias{vasco1988} +\alias{vasco88} +\title{Surface deformation associated with fluid withdrawl: an alternative formulation} +\usage{ +vasco88() + +.surface_g_pt(x = 0, x_src = 0, z_src = 0, nuu = 1/3) + +surface_displacement_point(x, depth, delV., B. = 1, rho_f. = 1000, tol, ...) +} +\arguments{ +\item{x}{numeric; the horizontal (or radial) distance from the source} + +\item{depth}{numeric; the depth below the surface to the source} + +\item{delV.}{numeric; uniform change in pore-fluid volume (positive = loss)} + +\item{B.}{numeric; Skempton's ratio for the reservoir material} + +\item{rho_f.}{numeric; the density of the fluid in the reservoir material} + +\item{tol}{numeric; the numerical tolerance in the integration; if +supplied, this should be much smaller than the smallest difference between +any values in \code{x} divided by the depth} + +\item{...}{additional arguments passed to \code{\link{.surface_g_pt}}} + +\item{help}{logical; load documentation for \code{\link{vasco88}}} +} +\description{ +Surface deformation associated with fluid withdrawl: an alternative formulation +} +\examples{ +# Reproduce... +r <- seq(-10,10,by=0.2) # km +r.m <- r*1e3 +xx <- surface_displacement_point(r.m, depth=2e3, delV.=1e7) +plot(uz ~ x, xx, type='l') +} +\seealso{ +\code{\link{segall85}}, or \code{\link{segall89}}; \code{\link{Simple-deformation}} +} +