# demo/RFOC.R In RFOC: Graphics for Spherical Distributions and Earthquake Focal Mechanisms

```############  to illustrate the graphical representation of
############   an earthquake focal mechanism (for teaching, for example)

print("Demonstration of vectors for Focal Mechanisms")

TEACHFOC(65, 32, -34, up=TRUE)

#############  animation of faults for teaching:sfanc
if(FALSE)
{

##############  these are animations...to be run outside of the demo
anim= seq(from=0, to=1, by=.1)

strikeslip.fault(anim=anim, Light=c(45,90) )

normal.fault(45, anim=anim, KAPPA=4, Light=c(-20, 80))

thrust.fault(anim=anim, KAPPA=4, Light=c(-20, 80))
}

###############  end animation part

print("Illustration of differnt kinds of faults")

###############   PLOT ALL three styles
par(mfrow=c(3,1))

anim=.5

############
###   graphics.off();  X11(w=8, h=10)

strikeslip.fault(anim=anim, Light=c(45,90) )

MFOC1 = SDRfoc(65,90,1, u=FALSE, ALIM=c(-1,-1, +1, +1), PLOT=FALSE)
Fcol1 = foc.color(foc.icolor(MFOC1\$rake1), pal=1)

justfocXY( MFOC1, fcol = Fcol1, .5, .7 , focsiz=.4 )
title("Strike-slip fault")

############
normal.fault(45, anim=anim, KAPPA=4, Light=c(-20, 80))

MFOC2 = SDRfoc(135,45,-90, u=FALSE, ALIM=c(-1,-1, +1, +1), PLOT=FALSE)
Fcol2 = foc.color(foc.icolor(MFOC2\$rake1), pal=1)

justfocXY( MFOC2, fcol = Fcol2, .5, 1 , focsiz=.4 )
title("Normal fault")

############

thrust.fault(anim=anim, KAPPA=4, Light=c(-20, 80))

MFOC3 = SDRfoc(135,45,90, u=FALSE, ALIM=c(-1,-1, +1, +1), PLOT=FALSE)
Fcol3 = foc.color(foc.icolor(MFOC3\$rake1), pal=1)

justfocXY( MFOC3, fcol = Fcol3, 0, -1 , focsiz=.4 )

title("Reverse (Thrust) fault")

###################################
############################# ############################# #############################
############################# ############################# #############################
############################# ############################# #############################

data(KAMCORN)
par(mfrow=c(1,1))

print("Plot Focal Mechanisms Geographically")

plot(KAMCORN\$LON, KAMCORN\$LAT, xlab="LON", ylab="LAT" , main="Kamchatka-Aleutian Inersection", asp=1)

####  set up these vectors for later use
Paz =vector()
Pdip =vector()
Taz =vector()
Tdip =vector()
h = vector()
v = vector()

IFcol = vector()
Fcol = vector()

for(i in 1:length(KAMCORN\$LON))
{
Msdr = CONVERTSDR(KAMCORN\$STRIKE[i], KAMCORN\$DIP[i], KAMCORN\$RAKE[i]   )
MEC = MRake(Msdr\$M)
MEC\$UP = FALSE
IFcol[i] = foc.icolor(MEC\$rake1)
Fcol[i] = foc.color(IFcol[i], 1)

az1 = Msdr\$M\$az1
dip1 = Msdr\$M\$d1
az2 = Msdr\$M\$az2
dip2 = Msdr\$M\$d2
BBB = Bfocvec(az1, dip1,  az2,  dip2)
V = ternfoc.point(BBB\$Bdip, Msdr\$M\$pd, Msdr\$M\$td )
Paz[i] = Msdr\$M\$paz
Pdip[i] = Msdr\$M\$pd
Taz[i] = Msdr\$M\$taz
Tdip[i] = Msdr\$M\$td
h[i] = V\$h
v[i] = V\$v

justfocXY( MEC, fcol = Fcol[i], KAMCORN\$LON[i], KAMCORN\$LAT[i] , focsiz=.2 )
}

print("Plot Focal Mechanisms showing only fault plane and slip vector point")

#############  plot only fault planes with slip vector:

plot(KAMCORN\$LON, KAMCORN\$LAT, xlab="LON", ylab="LAT" , main="Kamchatka-Aleutian Inersection", asp=1, type='n')

for(i in 1:length(KAMCORN\$LON) )
{
Msdr = CONVERTSDR(KAMCORN\$STRIKE[i], KAMCORN\$DIP[i], KAMCORN\$RAKE[i]   )
MEC = MRake(Msdr\$M)
MEC\$UP = FALSE
##  points(fxy\$x, fxy\$y)
nipXY( MEC,  KAMCORN\$LON[i], KAMCORN\$LAT[i] , fcol = Fcol[i], nipcol=Fcol[i],  focsiz=.2  )

}

print("show summary steronet plot of P and T-axes from Kamchatka Data")

####################################################
####################################################
####################################################  combined P-T axes
####################################################

net()
PZZ =  focpoint(Paz, Pdip, col='red',  pch=3, lab="", UP=FALSE)
TZZ =  focpoint(Taz, Tdip, col='blue',  pch=4, lab="", UP=FALSE)
text(0,1.04,labels="P&T-axes Projected", font=2, cex=1.2)
legend("topright",c("P", "T") , col=c('red','blue') , pch=c(3,4))

net()
PZZ =  focpoint(Paz, Pdip, col='red',  pch=3, lab="", UP=FALSE)

ALPH = alpha95(Paz, Pdip)

net()
TZZ =  focpoint(Taz, Tdip, col='blue',  pch=4, lab="", UP=FALSE)
ALPH = alpha95(Taz, Tdip)

print("Contour P and T-axis poles")

library(MASS)

KP = kde2d(PZZ\$x, PZZ\$y, n=50, lims=c(-1, 1, -1, 1))
KT = kde2d(TZZ\$x, TZZ\$y, n=50, lims=c(-1, 1, -1, 1) )

par(mfrow=c(1,3))
#################  plot 1
par(mai=c(0.2,0,.2,0))
CC = PLTcirc(PLOT=FALSE, add=FALSE,  ndiv=36,  angs=c(-pi, pi))

GEOmap::antipolygon(CC\$x,CC\$y,col="white")

focpoint(Paz, Pdip, col='red',  pch=3, lab="", UP=FALSE)
text(0,1.04,labels="P-axes 2D Density", font=2, cex=1.2)

#################  plot 2
CC = PLTcirc(PLOT=FALSE, add=FALSE,  ndiv=36,  angs=c(-pi, pi))

GEOmap::antipolygon(CC\$x,CC\$y,col="white")
focpoint(Taz, Tdip, col='blue',  pch=4, lab="", UP=FALSE)
text(0,1.04,labels="T-axes 2D Density", font=2, cex=1.2)

#################  plot 3
CC = PLTcirc(PLOT=FALSE, add=FALSE,  ndiv=36,  angs=c(-pi, pi))

##  focpoint(Paz, Pdip, col='red',  pch=3, lab="", UP=FALSE)

GEOmap::antipolygon(CC\$x,CC\$y,col="white")
text(0,1.04,labels="Combined P-T 2D Density", font=2, cex=1.2)

print("Show summary ternary plots of P and T-axes from Kamchatka Data")

####################################################
####################################################
####################################################  Ternary Plots
####################################################
graphics.off()
par(mfrow=c(1,1))

PlotTernfoc(h,v,x=0, y=0, siz=1, fcols=Fcol, add=FALSE, LAB=TRUE)

MFOC1 = SDRfoc(65,90,1, u=FALSE, ALIM=c(-1,-1, +1, +1), PLOT=FALSE)
Fcol1 = foc.color(foc.icolor(MFOC1\$rake1), pal=1)
MFOC2 = SDRfoc(135,45,-90, u=FALSE, ALIM=c(-1,-1, +1, +1), PLOT=FALSE)
Fcol2 = foc.color(foc.icolor(MFOC2\$rake1), pal=1)
MFOC3 = SDRfoc(135,45,90, u=FALSE, ALIM=c(-1,-1, +1, +1), PLOT=FALSE)
Fcol3 = foc.color(foc.icolor(MFOC3\$rake1), pal=1)

justfocXY( MFOC3, fcol = Fcol3, 1.2, -0.9,  focsiz=.4 )
justfocXY( MFOC2, fcol = Fcol2, -1.2, -0.9, focsiz=.4 )
justfocXY( MFOC1, fcol = Fcol1, 0, 1.414443+.2, focsiz=.4 )

mtext("Ternary Plot of focal mecahnisms", side=1, line = 1, font=2)

######################################
graphics.off()
plot(KAMCORN\$LON, KAMCORN\$LAT, xlab="LON", ylab="LAT" , main="Kamchatka-Aleutian Inersection", asp=1, type='n')

u = par("usr")

RI = RectDense(KAMCORN\$LON , KAMCORN\$LAT, icut=20, u=u, ndivs=6)

points(KAMCORN\$LON , KAMCORN\$LAT, pch='.', cex=2)

rect(RI\$icorns[,1],RI\$icorns[,2],RI\$icorns[,3],RI\$icorns[,4], col=NA, border='blue')

######################################

for(i in 1:length(RI\$ipass))
{
flag = KAMCORN\$LON>RI\$icorns[i,1]& KAMCORN\$LAT>RI\$icorns[i,2] & KAMCORN\$LON<RI\$icorns[i,3] & KAMCORN\$LAT<RI\$icorns[i,4]
jh =h[flag]
jv= v[flag]
}
legend("topright",focleg(1:7), col=foc.color(1:7, pal=1), pch=16)

####################################################
####################################################
####################################################  Ternary Plots
print("show geographics distribution of P-T axes around corner")

plot(KAMCORN\$LON, KAMCORN\$LAT, xlab="LON", ylab="LAT" , main="Kamchatka-Aleutian Inersection", asp=1, type='n')
points(KAMCORN\$LON, KAMCORN\$LAT, pch='.', cex=2)

KPspat =  matrix(NA, nrow=length(RI\$ipass), ncol=10)
KTspat =  matrix(NA, nrow=length(RI\$ipass), ncol=10)
colnames(KTspat) = c("x", "y", "n", "Ir",  "Dr", "R", "K", "S", "Alph95" , "Kappa")
colnames(KPspat) = c("x", "y", "n", "Ir",  "Dr", "R", "K", "S", "Alph95" , "Kappa")

for(i in 1:length(RI\$ipass))
{
flag = KAMCORN\$LON>RI\$icorns[i,1]& KAMCORN\$LAT>RI\$icorns[i,2] & KAMCORN\$LON<RI\$icorns[i,3] & KAMCORN\$LAT<RI\$icorns[i,4]
paz=Paz[flag]
pdip=Pdip[flag]
taz=Taz[flag]
tdip=Tdip[flag]
x=mean(RI\$icorns[i,c(1,3)])
y=mean(RI\$icorns[i,c(2,4)])
siz=(RI\$icorns[1,3]-RI\$icorns[1,1])/2.5

PlotPTsmooth(paz, pdip, x=x, y=y, siz=siz, border=NA, bcol='white' , LABS=FALSE, add=FALSE, IMAGE=TRUE, CONT=FALSE)
PlotPTsmooth(taz, tdip, x=x, y=y, siz=siz, border=NA, bcol='white' , LABS=FALSE, add=TRUE, IMAGE=FALSE, CONT=TRUE)

######dev.set(2)

ALPH = alpha95(paz, pdip)
n = length( paz)
KPspat[i,] =   c(x, y, n, ALPH\$Ir,  ALPH\$Dr, ALPH\$R, ALPH\$K, ALPH\$S, ALPH\$Alph95 , ALPH\$Kappa)

ALPH = alpha95(taz, tdip)
KTspat[i,] =    c(x, y, n, ALPH\$Ir,  ALPH\$Dr, ALPH\$R, ALPH\$K, ALPH\$S, ALPH\$Alph95 , ALPH\$Kappa)

}
```

## Try the RFOC package in your browser

Any scripts or data that you put into this service are public.

RFOC documentation built on Jan. 20, 2018, 9:46 a.m.