R : Copyright 2005, The R Foundation for Statistical Computing
Version 2.1.0  (2005-04-18), ISBN 3-900051-07-0

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for a HTML browser interface to help.
Type 'q()' to quit R.

> ### * <HEADER>
> ###
> attach(NULL, name = "CheckExEnv")
> assign(".CheckExEnv", as.environment(2), pos = length(search())) # base
> ## add some hooks to label plot pages for base and grid graphics
> setHook("plot.new", ".newplot.hook")
> setHook("persp", ".newplot.hook")
> setHook("grid.newpage", ".gridplot.hook")
> 
> assign("cleanEx",
+        function(env = .GlobalEnv) {
+ 	   rm(list = ls(envir = env, all.names = TRUE), envir = env)
+            RNGkind("default", "default")
+ 	   set.seed(1)
+    	   options(warn = 1)
+ 	   delayedAssign("T", stop("T used instead of TRUE"),
+ 		  assign.env = .CheckExEnv)
+ 	   delayedAssign("F", stop("F used instead of FALSE"),
+ 		  assign.env = .CheckExEnv)
+ 	   sch <- search()
+ 	   newitems <- sch[! sch %in% .oldSearch]
+ 	   for(item in rev(newitems))
+                eval(substitute(detach(item), list(item=item)))
+ 	   missitems <- .oldSearch[! .oldSearch %in% sch]
+ 	   if(length(missitems))
+ 	       warning("items ", paste(missitems, collapse=", "),
+ 		       " have been removed from the search path")
+        },
+        env = .CheckExEnv)
> assign("..nameEx", "__{must remake R-ex/*.R}__", env = .CheckExEnv) # for now
> assign("ptime", proc.time(), env = .CheckExEnv)
> grDevices::postscript("SoPhy-Examples.ps")
> assign("par.postscript", graphics::par(no.readonly = TRUE), env = .CheckExEnv)
> options(contrasts = c(unordered = "contr.treatment", ordered = "contr.poly"))
> options(warn = 1)    
> library('SoPhy')
> 
> assign(".oldSearch", search(), env = .CheckExEnv)
> assign(".oldNS", loadedNamespaces(), env = .CheckExEnv)
> cleanEx(); ..nameEx <- "CDE"
> 
> ### * CDE
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: ADE
> ### Title: Advection-dispersion equation
> ### Aliases: ADE CDE
> ### Keywords: spatial
> 
> ### ** Examples
> 
> par(cex=1.5, mar=c(4.2,4,0.2,0.4))
> z <- ADE(z=1:100, t=10^(0:3), C0=1, d=1, v=1)
> matplot(z, outer(-1:-100, rep(1,4)), type='l',
+         xlab='concentration', ylab='z', 
+         lwd=3, cex.lab=1.5, cex.axis=1.5, col=1)
> legend(1.0, -100, legend=paste('t=', 10^(0:3), sep=''), xj=1, yj=0,
+         lwd=3, lty=1:4, cex=1.5, col = 1)
> 
> 
> 
> graphics::par(get("par.postscript", env = .CheckExEnv))
> cleanEx(); ..nameEx <- "read.swms2d.table"
> 
> ### * read.swms2d.table
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: read.swms2d.table
> ### Title: Standard SWMS2D input files
> ### Aliases: read.swms2d.table
> ### Keywords: spatial
> 
> ### ** Examples
> 
> ##################################################################
> ## global printing paramters and functions
> ##################################################################
> par(cex=1.5, mar=c(4.2,4,0.2,0.2))
> my.legend <- function(lu.x, lu.y, zlim, col, cex=1.5) {
+     ## uses already the legend code of R-1.3.0
+     cn <- length(col)
+     filler <- vector("character", length=(cn-3)/2)
+     legend(lu.x, lu.y, y.i=0.04, x.i=0.1, yj=0,
+            legend=c(format(zlim[2], dig=2), filler,
+              format(mean(zlim), dig=2), filler,
+              format(zlim[1], dig=2)),
+            lty=1, col=rev(col), cex=cex)
+   }
> 
> ##################################################################
> ##  SWMS2D, Example 1
> ##################################################################
> 
> path <- paste(system.file(package='SoPhy'), 'swms2d', sep="/")
> x <- read.swms2d.table(path)
> x$TPrint <- seq(10, 5400, 10)
> 
> z <- swms2d(x)$hQ[3, , ]
> i <- ((length(x$nCodeM$z) / 2):1) * 2 - 1
> image(x=c(0, x$TPrint) / 60, y=x$nCodeM$z[i] - max(x$nCodeM$z),
+       z=t(z[i, ]), xlab='time [min]', ylab='z [cm]',
+       col=grey(seq(1, 0.15, -0.01)), cex.lab=1.5, cex.axis=1.5)
> 
> my.legend(0, -max(x$nC$z[i]), zlim=range(z),
+           col=grey(seq(1,0.15,-0.01)))
> 
> 
> 
> graphics::par(get("par.postscript", env = .CheckExEnv))
> cleanEx(); ..nameEx <- "swms2d"
> 
> ### * swms2d
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: swms2d
> ### Title: SWMS2D function call
> ### Aliases: swms2d SWMS2D
> ### Keywords: spatial
> 
> ### ** Examples
> 
> ##################################################################
> ##  SWMS2D, Example 1, stochastically modified material properties
> ##################################################################
> # E = exp(mu + 1/2 * sd^2)
> # Var = exp(2*mu + sd^2)*(exp(sd^2) - 1)
> 
> path <- paste(system.file(package='SoPhy'), 'swms2d', sep="/")
> x <- read.swms2d.table(path)
> x$TPrint <- 20 * 60
> sd <- 0.3 ## then sd of log.gauss is 0.306878
> mu <- log(x$nCodeM[, 10]) - 0.5 * sd^2
> n <- 100 # better 10000
> front <- integer(n)
> 
> for (j in 1:n) {
+   cat(j, "\n")
+   x$nCodeM[, 10] <- rlnorm(length(x$nCodeM$z), m = mu, sd = sd)
+   x$nCodeM[, 11] <- x$nCodeM[, 10]^(-2)
+   z <- swms2d(x)$hQ[3, , ]
+   front[j] <- (min(which(z[-1:-2, 2] / z[-1:-2, 1] < 1.005)) %/% 2) + 1
+ }
1 
2 
3 
4 
5 
6 
7 
8 
9 
10 
11 
12 
13 
14 
15 
16 
17 
18 
19 
20 
21 
22 
23 
24 
25 
26 
27 
28 
29 
30 
31 
32 
33 
34 
35 
36 
37 
38 
39 
40 
41 
42 
43 
44 
45 
46 
47 
48 
49 
50 
51 
52 
53 
54 
55 
56 
57 
58 
59 
60 
61 
62 
63 
64 
65 
66 
67 
68 
69 
70 
71 
72 
73 
74 
75 
76 
77 
78 
79 
80 
81 
82 
83 
84 
85 
86 
87 
88 
89 
90 
91 
92 
93 
94 
95 
96 
97 
98 
99 
100 
> 
> idx <- 24 
> cat('The probability that the front advances at most', 
+     max(x$nCodeM$z) - x$nCodeM$z[2 * idx],'cm is about',
+     format(mean(front <= idx), dig=2),'\n')
The probability that the front advances at most 19 cm is about 0.06 
> 
> par(cex=1, mar=c(4.2,4.2,1.2,0.4))
> distance <- max(x$nCodeM$z) - x$nCodeM$z[2 * front]
> h <- sort(unique(distance))
> dh <- min(diff(h) / 2)
> hist(distance, freq=FALSE,
+      breaks=seq(h[1]-dh, h[length(h)]+dh, 2 * dh),
+      main="Histogramm for the depth of the water front",
+      xlab="depth", cex.axis=1.5, cex.main=1.2, cex.lab=1.5)
> 
> 
> 
> graphics::par(get("par.postscript", env = .CheckExEnv))
> ### * <FOOTER>
> ###
> cat("Time elapsed: ", proc.time() - get("ptime", env = .CheckExEnv),"\n")
Time elapsed:  26.72 0.16 26.87 0 0 
> grDevices::dev.off()
null device 
          1 
> ###
> ### Local variables: ***
> ### mode: outline-minor ***
> ### outline-regexp: "\\(> \\)?### [*]+" ***
> ### End: ***
> quit('no')