Added Vasco 1988 formulation, and tested for integration stability. T…

…his produces a more understandable framework than segall 1985

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 <[email protected]>
 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

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)

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

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)
+}

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}
 

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}}
+}
+