Added Segall subsidence -- time varying not working.

DESCRIPTION

@@ -1,8 +1,8 @@
 Package: deform
 Type: Package
 Title: Crustal deformation tools in R
-Version: 0.0-1
-Date: 2014-10-20
+Version: 0.0-2
+Date: 2014-11-04
 Author: Andrew J. Barbour; some of these tools were
     adapted from Matlab code by Brendan Crowell (UW),
     and Francois Beauducel (IPGP).

NAMESPACE

@@ -1,10 +1,20 @@
 # Generated by roxygen2 (4.0.1): do not edit by hand
 
+export(.surface_g)
 export(D_project)
 export(agency.root)
 export(degreelen)
 export(ellipsoid)
+export(erf)
+export(erfc)
 export(getEPSG)
+export(hydraulic_diffusivity)
+export(ierf)
+export(ierfc)
+export(rigid_hydraulic_diffusivity)
+export(segall85)
+export(surface_displacement)
+export(timevarying_surface_displacement)
 exportPattern("^[[:alpha:]]+")
 importFrom(Rcpp,evalCpp)
 useDynLib(deform)

NOBUILD/Sandbox/segall.R

@@ -0,0 +1,19 @@
+library(deform)
+
+y. <- c(-10:-2, seq(-1.9,1.9,by=0.1), 2:10)
+su <- surface_displacement(y.*1e3, C.=1e13, z_src=1e3)
+#plot(uz ~ y, su, type="b", ylim=c(-2,1)*10)
+#lines(uy ~ y, su, type="b", col="red")
+#lines((uy+uz) ~ y, su, type="b", col="blue")
+
+# Fig 8
+#y. <- 0:10
+t. <- c(0.1,0.5,1:5)
+Vdot <- 2e6
+D. <- 1e3
+L. <- 10e3
+B. <- 0.6
+c. <- 0.1
+matplot(y.,t.,
+  timevarying_surface_displacement(y., t., Vdot, B., L., D., c., y_src.=1),
+  type="l")

R/deform-utilities.R

@@ -0,0 +1,24 @@
+#' Error functions
+#' @param x numeric;
+#' @export
+erf <- function(x){
+  2 * pnorm(x * sqrt(2)) - 1
+}
+
+#' @rdname erf
+#' @export
+erfc <- function(x){
+  2 * pnorm(x * sqrt(2), lower = FALSE)
+}
+
+#' @rdname erf
+#' @export
+ierf <- function (x){
+  qnorm((1 + x)/2)/sqrt(2)
+}
+
+#' @rdname erf
+#' @export
+ierfc <- function (x){
+  qnorm(x/2, lower = FALSE)/sqrt(2)
+}

R/degreelen.R

@@ -5,6 +5,7 @@
 #' @author A.J.Barbour, adapted from \code{degreelen.m} by B.W.Crowell
 #' @references \url{http://en.wikipedia.org/wiki/World_Geodetic_System}
 #' @examples
+#' \dontrun{
 #' lats <- 30:32
 #' # Lon should be about 111.4 km, and Lat 96 - 94 km
 #' degreelen(lats)
@@ -19,6 +20,7 @@
 #' options(digits=13)
 #' lats <- 30.1111111123
 #' rbind(degreelen(lats,"IRTF"),degreelen(lats,"WGS84"))
+#' }
 degreelen <- function(lat, model){
   #
   pic <- pi/180

R/diffusivity.R

@@ -0,0 +1,15 @@
+#' Calculate hydraulic diffusivity of a fluid saturated medium
+#' @export
+hydraulic_diffusivity <- function(k., eta., mu., B., nu.=1/4, nuu.=1/3){
+  # segall85 eq 11
+  n1 <- 1 - nu.
+  n1u <- 1 - nuu.
+  n2 <- 1 - 2*nu.
+  k. * B.^2 * (2*mu.*n1/n2) * (n2 * (1 + nuu.)^2 / n1u / (nuu. - nu.)) / eta. / 9
+}
+#' @rdname hydraulic_diffusivity
+#' @export
+rigid_hydraulic_diffusivity <- function(k., eta., phi., Beta.){
+  # segall85 eq 12
+  k./eta./phi./Beta.
+}

R/segall85.R

@@ -0,0 +1,64 @@
+#' Surface deformation associated with fluid withdrawl
+#' @name segall85
+#' @export
+#' @examples
+#' \dontrun{
+#' y. <- c(-10:-2, seq(-1.9,1.9,by=0.1), 2:10)
+#' su <- surface_displacement(y.*1e3, z_src=1e3)
+#' plot(uz ~ y, su, type="b")
+#' }
+segall85 <- function(){
+}
+
+#  spatial coordinates in segall
+#   x is positive downwards
+#   y is positive away from source
+#   z is along line source
+#  source coordinates:
+#   z_src is zeta (depth)
+#   y_src is xi (distance away)
+#  parameters:
+#   mu shear modulus
+#   C units of force, proportional to the source strength
+#' @rdname segall85
+#' @export
+surface_displacement <- function(y, C.=1, mu.=1e9, ...){
+  sg <- .surface_g(y, ...)
+  c. <- C./mu.
+  sg$uz <- sg$gz * c.
+  sg$uy <- sg$gy * c.
+  return(sg)
+}
+
+#' @rdname segall85
+#' @export
+.surface_g <- function(y=0, z_src=0, y_src=0, nuu=1/3){
+  # segall85 C9
+  gz <- -2*(1 - nuu)*z_src/(z_src^2 + (y - y_src)^2)
+  gy <-  2*(1 - nuu)*(y - y_src)/(z_src^2 + (y - y_src)^2)
+  data.frame(y, gz, gy)
+}
+
+# Time varying deformation associatedw with fluid extraction
+
+#' @rdname segall85
+#' @export
+timevarying_surface_displacement <- function(y, Time, Vdot., B., L., D., HD., nuu.=1/3, y_src.=0){
+  # segall85 eq 26
+  # at each time slice, calculate a profile
+  .tvsd <- function(yi, ti, .Vdot, .B, .L, .D, .HD, .nuu, .y_src){
+    message(paste(.Vdot, .B, .L, .D, .HD, .nuu, .y_src))
+    # source
+    mod <- 2*.B*(1 + .nuu)*.Vdot*.D*sqrt(ti/.HD)/(3*pi*.L)
+    # line source correction
+    sc <- if (.y_src==0){
+      1
+    } else {
+      Time.src <- sqrt(.y_src^2 / 4 / .HD / ti)
+      err <- qnorm(Time.src/2, lower = FALSE)/sqrt(2)  #ierfc
+      sum(err/(.D^2 + (yi - .y_src)^2), na.rm=TRUE)
+    }
+    return(mod*sc)
+  }
+  outer(X=y, Y=Time, FUN=.tvsd, .Vdot=Vdot., .B=B., .L=L., .D=D., .HD=HD., .nuu=nuu., .y_src=y_src.)
+}

man/degreelen.Rd

@@ -17,6 +17,7 @@ ellipsoid(model = c("IRTF", "WGS84"))
 Compute the length of 1-degree of latitude and longitude
 }
 \examples{
+\dontrun{
 lats <- 30:32
 # Lon should be about 111.4 km, and Lat 96 - 94 km
 degreelen(lats)
@@ -32,6 +33,7 @@ options(digits=13)
 lats <- 30.1111111123
 rbind(degreelen(lats,"IRTF"),degreelen(lats,"WGS84"))
 }
+}
 \author{
 A.J.Barbour, adapted from \code{degreelen.m} by B.W.Crowell
 }

man/erf.Rd

@@ -0,0 +1,23 @@
+% Generated by roxygen2 (4.0.1): do not edit by hand
+\name{erf}
+\alias{erf}
+\alias{erfc}
+\alias{ierf}
+\alias{ierfc}
+\title{Error functions}
+\usage{
+erf(x)
+
+erfc(x)
+
+ierf(x)
+
+ierfc(x)
+}
+\arguments{
+\item{x}{numeric;}
+}
+\description{
+Error functions
+}
+

man/hydraulic_diffusivity.Rd

@@ -0,0 +1,14 @@
+% Generated by roxygen2 (4.0.1): do not edit by hand
+\name{hydraulic_diffusivity}
+\alias{hydraulic_diffusivity}
+\alias{rigid_hydraulic_diffusivity}
+\title{Calculate hydraulic diffusivity of a fluid saturated medium}
+\usage{
+hydraulic_diffusivity(k., eta., mu., B., nu. = 1/4, nuu. = 1/3)
+
+rigid_hydraulic_diffusivity(k., eta., phi., Beta.)
+}
+\description{
+Calculate hydraulic diffusivity of a fluid saturated medium
+}
+

man/segall85.Rd

@@ -0,0 +1,28 @@
+% Generated by roxygen2 (4.0.1): do not edit by hand
+\name{segall85}
+\alias{.surface_g}
+\alias{segall85}
+\alias{surface_displacement}
+\alias{timevarying_surface_displacement}
+\title{Surface deformation associated with fluid withdrawl}
+\usage{
+segall85()
+
+surface_displacement(y, C. = 1, mu. = 1e+09, ...)
+
+.surface_g(y = 0, z_src = 0, y_src = 0, nuu = 1/3)
+
+timevarying_surface_displacement(y, Time, Vdot., B., L., D., HD., nuu. = 1/3,
+  y_src. = 0)
+}
+\description{
+Surface deformation associated with fluid withdrawl
+}
+\examples{
+\dontrun{
+y. <- c(-10:-2, seq(-1.9,1.9,by=0.1), 2:10)
+su <- surface_displacement(y.*1e3, z_src=1e3)
+plot(uz ~ y, su, type="b")
+}
+}
+