-
Notifications
You must be signed in to change notification settings - Fork 1
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
- Loading branch information
1 parent
5803177
commit fbf7058
Showing
49 changed files
with
1,212 additions
and
205 deletions.
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -1,13 +1,16 @@ | ||
| Package: geophys | ||
| Type: Package | ||
| Title: Geophysics, Continuum Mechanics, Mogi Model | ||
| Version: 1.2-1 | ||
| Date: 2010-05-12 | ||
| Depends: R (>= 2.13), RPMG, RFOC, MASS, cluster, stats | ||
| Version: 1.3-6 | ||
| Date: 2013-01-16 | ||
| Depends: R (>= 2.15) | ||
| Imports: RPMG, RSEIS, RFOC, GEOmap, cluster | ||
| Suggests: stats | ||
| Author: Jonathan M. Lees | ||
| Maintainer: Jonathan M. Lees<jonathan.lees@unc.edu> | ||
| Description: Geophysics, Continuum Mechanics, Mogi Model, Okada Model | ||
| License: GPL | ||
| Packaged: 2011-10-17 19:46:11 UTC; lees | ||
| Maintainer: Jonathan M. Lees <jonathan.lees@unc.edu> | ||
| Description: Geophysics, Continuum Mechanics, Mogi Model, Okada Model, Gravity | ||
| License: GPL (>= 2) | ||
| Packaged: 2014-01-06 17:33:59 UTC; lees | ||
| NeedsCompilation: no | ||
| Repository: CRAN | ||
| Date/Publication: 2011-10-19 07:39:56 | ||
| Date/Publication: 2014-01-13 08:25:19 |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -1,2 +1,7 @@ | ||
| # Export all names | ||
| exportPattern(".") | ||
| import(RPMG) | ||
| import(RSEIS) | ||
| import(RFOC) | ||
| import(GEOmap) | ||
| importFrom(cluster, ellipsoidhull, predict.ellipsoid) |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,110 @@ | ||
| BMOD<-function(bill, nstn=100, PLOT=TRUE, obs=NULL ) | ||
| { | ||
| ####### two parts - model and plotting | ||
| if(missing(nstn)) nstn = 100 | ||
| if(missing(PLOT)) PLOT=TRUE | ||
|
|
||
| ZCOLS = RPMG::pastel.colors(12, seed=1 ) | ||
|
|
||
| DZ = abs( bill$zmax - bill$zmin) | ||
|
|
||
| PZ = DZ*0.2 | ||
|
|
||
| ZMAX = bill$zmin+PZ | ||
| ZMIN = bill$zmin+DZ*0.02 | ||
|
|
||
| xstart = bill$xmin | ||
| xend = bill$xmax | ||
|
|
||
| xs = seq(from=xstart, by=(xend-xstart)/nstn , length=nstn) | ||
| zs = rep(0, length=length(xs)) | ||
|
|
||
| den = 0.2 | ||
| n = length(bill$mod) | ||
|
|
||
| ALLgravZ = rep(0, length(xs)) | ||
|
|
||
| for(i in 1:bill$n) | ||
| { | ||
| den = bill$mod[[i]]$rho | ||
| if(is.na(den)) den=0.2 | ||
| if(is.null(den)) den=0.2 | ||
|
|
||
|
|
||
| pol = flipZEE(bill$mod[[i]]) | ||
|
|
||
| dc = dircheck(pol) | ||
| ## print(dc) | ||
| L1 = length(which(dc>0)) | ||
| if(L1>=(length(pol$x)/2)) | ||
| { | ||
| cat(paste("Reversing", i, sep=":") , sep="\n") | ||
| pol = rev2RH(pol) | ||
| } | ||
|
|
||
| Ngrav = DGzx(xs, zs, pol$x, pol$y, den) | ||
| ALLgravZ = ALLgravZ + Ngrav$Gz | ||
| } | ||
| ## plot(xs, ALLgravZ ) | ||
| ######################### | ||
| if(PLOT==TRUE) | ||
| { | ||
| plot(c(bill$xmin, bill$xmax), c( bill$zmax, ZMAX), type='n' , axes=FALSE, ann=FALSE ) | ||
| grid() | ||
| title(xlab="X-m", ylab="Depth, m") | ||
| rect(bill$xmin, bill$zmin, bill$xmax, bill$zmax, col=NA, border="black", xpd=TRUE ) | ||
|
|
||
| for(i in 1:bill$n) | ||
| { | ||
|
|
||
| den = bill$mod[[i]]$rho | ||
| pol = (bill$mod[[i]]) | ||
| polygon(pol$x, pol$y, col=ZCOLS[i], border="black") | ||
|
|
||
| ### need to check to make sure the polygon is right handed | ||
|
|
||
| text(pol$x, pol$y, labels=1:length(pol$x)) | ||
| cenP = centroid(pol) | ||
| text( cenP[1], cenP[2], paste(i,":", den, sep="" )) | ||
| } | ||
|
|
||
| g1 = min(ALLgravZ) | ||
| g2 = max(ALLgravZ) | ||
| RGz = RESCALE(ALLgravZ, ZMIN, ZMAX, g1, g2) | ||
|
|
||
| rect( bill$xmin, ZMIN, bill$xmax, ZMAX, col='white', border=grey(.9) ) | ||
|
|
||
| axis(1) | ||
| paz = pretty(c( bill$zmin, bill$zmax)) | ||
| axis(2, at=paz) | ||
|
|
||
| ## gravaxis | ||
| gzmin = min(ALLgravZ) | ||
| gzmax = max(ALLgravZ) | ||
|
|
||
| text(bill$xmax, ZMIN, labels=format(gzmin, digits=4), pos=4, xpd=TRUE ) | ||
| text(bill$xmax, ZMAX, labels=format(gzmax, digits=4), pos=4, xpd=TRUE ) | ||
|
|
||
| pgz = pretty(range(ALLgravZ)) | ||
| pgz = pgz[ pgz>gzmin & pgz<gzmax ] | ||
| rpgz = RESCALE(pgz, ZMIN, ZMAX, g1, g2) | ||
|
|
||
| segments(bill$xmin, rpgz, bill$xmax, rpgz, lty=2, col=grey(.9) ) | ||
| lines(xs, RGz, col='blue') | ||
|
|
||
|
|
||
| if(!is.null(obs)) | ||
| { | ||
| ## obsGz = RESCALE(obs$g, ZMIN, ZMAX, g1, g2) | ||
| obsGz = RESCALE(obs$g, ZMIN, ZMAX, min(obs$g), max(obs$g) ) | ||
| points(obs$x, obsGz, pch=3, col='red') | ||
| text(bill$xmin, ZMIN, labels=format(min(obs$g), digits=4), pos=2, xpd=TRUE, col='red' ) | ||
| text(bill$xmin, ZMAX, labels=format( max(obs$g), digits=4), pos=2, xpd=TRUE, col='red' ) | ||
|
|
||
| } | ||
|
|
||
|
|
||
| } | ||
|
|
||
| invisible(ALLgravZ) | ||
| } |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,81 @@ | ||
| DGzx<-function(xs, zs, xv, zv, den) | ||
| { | ||
| ### Compute the gravity anomaly following Won and Bevis | ||
| twopi = 2*pi | ||
| con=13.3464E-03 | ||
|
|
||
| nvert = length(xv) | ||
|
|
||
| #### make sure the polygon is closed | ||
|
|
||
| if(xv[1]!=xv[nvert] & zv[1]!=zv[nvert]) | ||
| { | ||
| xv = c(xv, xv[1]) | ||
| zv = c(zv, zv[1]) | ||
|
|
||
| } | ||
| nvert = length(xv) | ||
| ### create output vector | ||
| gravz = rep(NA, length(xs)) | ||
| gravx = rep(NA, length(xs)) | ||
| ### loop over the stations | ||
| for(i in 1:length(xs)) | ||
| { | ||
|
|
||
| xst = xs[i] | ||
| zst = zs[i] | ||
|
|
||
| x1 = xv[1:(nvert-1)]-xst; | ||
| z1 = zv[1:(nvert-1)]-zst; | ||
| x2 = xv[2:(nvert)]-xst; | ||
| z2 = zv[2:(nvert)]-zst; | ||
|
|
||
| ### calculate the angles | ||
| theta1 = atan2(z1, x1); | ||
| theta2 = atan2(z2, x2); | ||
| ### need to get rid of pathology here. | ||
|
|
||
| dsign = sign(z1) != sign(z2) | ||
| if(any(dsign)) | ||
| { | ||
| theta1[ dsign & (x1*z2<x2*z1) & z2>=0] = theta1[ dsign & (x1*z2<x2*z1) & z2>=0]+twopi | ||
| theta2[ dsign & (x1*z2>x2*z1) & z1>=0] = theta2[ dsign & (x1*z2>x2*z1) & z1>=0 ]+twopi | ||
| } | ||
|
|
||
| dx = x2-x1; | ||
| dz = z2-z1; | ||
|
|
||
| r1 = sqrt(x1*x1 + z1*z1); | ||
| r2 = sqrt(x2*x2 + z2*z2); | ||
|
|
||
|
|
||
| dxz2 = (dx*dx + dz*dz) | ||
| A = dx*( x1*z2 - x2*z1 )/(dx*dx + dz*dz); | ||
|
|
||
| B = dz/dx; | ||
|
|
||
|
|
||
| ZEE = A*( (theta1 - theta2) + B*log(r2/r1)) | ||
| EX = A*(-((theta1 - theta2) )*B + log(r2/r1)) | ||
|
|
||
| ZEE[x1*z2 == x2*z1] = 0 | ||
| EX[x1*z2 == x2*z1] = 0 | ||
|
|
||
| ZEE[ (x1==0 & z1==0) | (x2==0 & z2==0) ] = 0 | ||
| EX[ (x1==0 & z1==0) | (x2==0 & z2==0) ] = 0 | ||
|
|
||
| ZEE[x1==x2] = x1[x1==x2] * log(r2[x1==x2]/r1[x1==x2]) | ||
|
|
||
| EX[x1==x2] = -1*x1[x1==x2] *(theta1[x1==x2] - theta2[x1==x2]) | ||
|
|
||
| Z = sum( ZEE ); | ||
| X = sum( EX ); | ||
|
|
||
| gravz[i] = con*den*Z | ||
| gravx[i] = con*den*X | ||
|
|
||
| } | ||
|
|
||
| invisible(list(Gz=gravz, Gx=gravx)) | ||
|
|
||
| } |
Oops, something went wrong.