### R code from vignette source 'focalmec.Rnw'

###################################################
### code chunk number 1: focalmec.Rnw:66-72
###################################################
library(RFOC)
data(PKAM)
numk = length(PKAM$LATS)
payr = paste(collapse="-", range(PKAM$yr))




###################################################
### code chunk number 2: fig1plot
###################################################
load(file='./EYEVIEW.RData')

`emptyplot` <-
function(...)
  {

  plot(c(-1,1), c(-1,1), asp=1, type='n', axes=FALSE, ann=FALSE)
}
`rotsym` <-
function(sym , alpha, reflect=FALSE)
  {
    if(missing(reflect)) { reflect=FALSE }
    rphi = pi*alpha/180
    cphi = cos(rphi)
    sphi = sin(rphi)


    ###########  rotation
   for(i in 1:length(sym)  )
      {

         if(reflect==TRUE)
           {

             sym[[i]]$x = (-1)*sym[[i]]$x

           }

        ax = cphi*sym[[i]]$x  + sphi*sym[[i]]$y
        ay = -sphi*sym[[i]]$x  + cphi*sym[[i]]$y


        sym[[i]]$x = ax
          sym[[i]]$y =ay
      }
    return(sym)
  }
`recentersym` <-
function(sym)
  {

    gx = rep(NA, 2)
    gy = rep(NA, 2)

    for(i in 1:length(sym)  )
      {

        gx = range(c(gx, sym[[i]]$x), na.rm=TRUE)
        gy = range(c(gy, sym[[i]]$y), na.rm=TRUE)

      }
    rx = gx[2] - gx[1]
    ry = gy[2] - gy[1]
    for(i in 1:length(sym)  )
      {
        sym[[i]]$x = 2* (((sym[[i]]$x-gx[1])/rx)-0.5)
        sym[[i]]$y = 2* (((sym[[i]]$y-gy[1])/ry)-0.5)


      }


    attr(sym, "DX") = c(-1,1)
    attr(sym, "DY") = c(-1,1)
    attr(sym, "OUT") = TRUE

     return(sym)

  }

`draw.sym` <-
function(x,  y, sym, siz=siz, col=rgb(83/255, 157/255, 194/255), border="black" )
  {
    if(missing(col)) { col=rep(rgb(83/255, 157/255, 194/255),length=length(x)) }
    if(missing(siz)) { siz=rep(.2,length=length(x)) }
    if(missing(border)) { border=rep("black",length=length(x)) }

    if(length(siz)==1) { siz=rep(siz,length=length(x)) }
    if(length(col)==1) { col=rep(col,length=length(x)) }
    if(length(border)==1) { border=rep(border,length=length(x)) }

    up = par("usr")
    ui = par("pin")

    ratx = (up[2]-up[1])/ui[1]
    raty=  (up[4]-up[3])/ui[2]

    for(i in 1:length(x))
      {
        x0 = x[i]
        y0 = y[i]

        usizx = siz[i]*ratx
        usizy = siz[i]*raty
###  print(i)
        A = prep.sym(sym, X=c(x0-usizx, x0+usizx),  Y=c(y0-usizy, y0+usizy),  delx=attr(sym, "DX"), dely=attr(sym, "DY"))
        attr(A, "OUT") = FALSE
        JREdraw(A)
      }
  }

`prep.sym` <-
function(KALL, X=c(0,1),  Y=c(0,1),  delx=c(0,1), dely=c(0,1))
  {
    ########  remove the mean value and scale everything from
    ########  -1 to 1  or -.5 to .5
    K = 0
    NALL = list()
     NKALL = length(KALL)
    if(NKALL<1) return
    for(i in 1:NKALL )
      {
        if(is.null(KALL[[i]])) next
       K = K+1
        NALL[[K]] = KALL[[i]]
        NALL[[K]]$x = RESCALE(KALL[[i]]$x, X[1], X[2], delx[1], delx[2])
        NALL[[K]]$y = RESCALE(KALL[[i]]$y, Y[1], Y[2], dely[1], dely[2])

      }

    attr(NALL, "DX") =  X
    attr(NALL, "DY") =  Y
     attr(NALL, "OUT") = TRUE
    return(NALL)

  }
`RESCALE` <-
function(x, nx1, nx2, minx, maxx)
{
  #    rescale a vector
  nx = nx1+(nx2-nx1)*(x-minx)/(maxx-minx)
  return(nx)
}

`JREdraw` <-
function(KALL)
  {

    if(attr(KALL,"OUT")==TRUE)
      {
        DX = attr(KALL, "DX")
        DY = attr(KALL, "DY")
        rect(DX[1], DY[1], DX[2] , DY[2], col=NA, border='black' )
      }

    NKALL = length(KALL)
    if(NKALL<1) return
    for(i in 1:NKALL )
      {if(is.null(KALL[[i]])) next

        if(KALL[[i]]$type=="lines")
          { lines(KALL[[i]]$x, KALL[[i]]$y, col=KALL[[i]]$col, lwd=KALL[[i]]$lwd) }

        if(KALL[[i]]$type=="points")
          { points(KALL[[i]]$x, KALL[[i]]$y, col=KALL[[i]]$col) }

        if(KALL[[i]]$type=="polygon")
          { polygon(KALL[[i]]$x, KALL[[i]]$y, col=KALL[[i]]$col, border=KALL[[i]]$border) }

        if(KALL[[i]]$type=="circles")
          { makecircles(KALL[[i]], col=KALL[[i]]$col, border=KALL[[i]]$border) }

        if(KALL[[i]]$type=="rects")
          { makerects(KALL[[i]], col=KALL[[i]]$col, border=KALL[[i]]$border,  lwd=KALL[[i]]$lwd) }

        if(KALL[[i]]$type=="arrows")
          { makearrows(KALL[[i]], col=KALL[[i]]$col,  lwd=KALL[[i]]$lwd) }

      }

  }



C = circle()

cex1 = 1.5

i=seq(from=0, to=-pi, length=180)
cx = cos(i);
cy = sin(i);
emptyplot()
title("Lower Hemisphere Projection View from Side")
####axis(1)
####axis(2)
lines(cx, cy)

i=seq(from=0, to=-pi/4, length=180)

cx = cos(i);
cy = sin(i);

polygon(c(0, cx, 0), c(0, cy, 0) , col=grey(0.85), border='black')

i=seq(from=-3*pi/4, to=-pi, length=180)

cx = cos(i);
cy = sin(i);

polygon(c(0, cx, 0), c(0, cy, 0) , col=grey(0.85), border='black')

segments(-0.25*cx[1], -0.25*cy[1], cx[1], cy[1], lwd=2, lty=2)
text(cx[1], cy[1], labels="Auxilliary Plane", srt=45, adj=c(1,0), font=2, xpd=TRUE, cex=cex1)

i=seq(from=0, to=-pi/4, length=180)

cx = cos(i);
cy = sin(i);

segments(-0.25*cx[180], -0.25*cy[180], cx[180], cy[180], lwd=2, lty=2)

segments(-0.25*cx[180], -0.25*cy[180], 0.25*cx[180], 0.25*cy[180], lwd=4, lty=1)

segments(-0.1*cx[180]+.01, -0.1*cy[180]+.01, 0.1*cx[180]+.01, 0.1*cy[180]+.01, lwd=2, lty=1)
segments(0.1*cx[180]+.01, 0.1*cy[180]+.01,0.1*cx[180]+.01, 0.1*cy[180]+.01+.05, lwd=2, lty=1)


segments(-0.1*cx[180]-0.01, -0.1*cy[180]-0.01, 0.1*cx[180]-0.01, 0.1*cy[180]-0.01, lwd=2, lty=1)

segments(-0.1*cx[180]-0.01, -0.1*cy[180]-0.01, -0.1*cx[180]-0.01, -0.1*cy[180]-0.01-.05, lwd=2, lty=1)


text(-0.25*cx[180], -0.25*cy[180] , labels="Fault Plane",  adj=c(1,0), font=2, xpd=TRUE, cex=cex1)
text(cx[180], cy[180] , labels="Fault Plane Projection",  adj=c(0,0), srt=-45, font=2, xpd=TRUE, cex=cex1)




abline(h=.35)
text(0,0.35, "Surface of Earth", pos=3, cex=cex1)

arrows(-1.05, 0, -1.05, -0.5)
text(-1.05, -0.25, "Depth", srt=-90, adj=c(1,-0.1), cex=cex1)

points(c(1, 0), c(0, -1), pch=20, cex=2, col=c("blue", "red") )

text(1,0, labels="T", pos=4, xpd=TRUE, font=2, col="blue", cex=cex1)

text(0,-1, labels="P", pos=1, xpd=TRUE, font=2, col="red", cex=cex1)

############  eye part
###  get.syms()

### data(EYEVIEW)
EYEVIEW2 = rotsym(EYEVIEW, -90, reflect=FALSE)

EYEVIEW2 = recentersym(EYEVIEW2)

draw.sym(0,0.8, EYEVIEW2, siz=.4)

MFOC2 = SDRfoc(0,45,-90, u=FALSE, ALIM=c(-1,-1, +1, +1), PLOT=FALSE)
    ###  Fcol2 = foc.color(foc.icolor(MFOC2$rake1), pal=1)
  Fcol2 =  grey(0.85)
    justfocXY( MFOC2, fcol = Fcol2, .9, .9 , size = c(.2,.2) )





###################################################
### code chunk number 3: focalmec.Rnw:382-383
###################################################
load(file='./EYEVIEW.RData')

`emptyplot` <-
function(...)
  {

  plot(c(-1,1), c(-1,1), asp=1, type='n', axes=FALSE, ann=FALSE)
}
`rotsym` <-
function(sym , alpha, reflect=FALSE)
  {
    if(missing(reflect)) { reflect=FALSE }
    rphi = pi*alpha/180
    cphi = cos(rphi)
    sphi = sin(rphi)


    ###########  rotation
   for(i in 1:length(sym)  )
      {

         if(reflect==TRUE)
           {

             sym[[i]]$x = (-1)*sym[[i]]$x

           }

        ax = cphi*sym[[i]]$x  + sphi*sym[[i]]$y
        ay = -sphi*sym[[i]]$x  + cphi*sym[[i]]$y


        sym[[i]]$x = ax
          sym[[i]]$y =ay
      }
    return(sym)
  }
`recentersym` <-
function(sym)
  {

    gx = rep(NA, 2)
    gy = rep(NA, 2)

    for(i in 1:length(sym)  )
      {

        gx = range(c(gx, sym[[i]]$x), na.rm=TRUE)
        gy = range(c(gy, sym[[i]]$y), na.rm=TRUE)

      }
    rx = gx[2] - gx[1]
    ry = gy[2] - gy[1]
    for(i in 1:length(sym)  )
      {
        sym[[i]]$x = 2* (((sym[[i]]$x-gx[1])/rx)-0.5)
        sym[[i]]$y = 2* (((sym[[i]]$y-gy[1])/ry)-0.5)


      }


    attr(sym, "DX") = c(-1,1)
    attr(sym, "DY") = c(-1,1)
    attr(sym, "OUT") = TRUE

     return(sym)

  }

`draw.sym` <-
function(x,  y, sym, siz=siz, col=rgb(83/255, 157/255, 194/255), border="black" )
  {
    if(missing(col)) { col=rep(rgb(83/255, 157/255, 194/255),length=length(x)) }
    if(missing(siz)) { siz=rep(.2,length=length(x)) }
    if(missing(border)) { border=rep("black",length=length(x)) }

    if(length(siz)==1) { siz=rep(siz,length=length(x)) }
    if(length(col)==1) { col=rep(col,length=length(x)) }
    if(length(border)==1) { border=rep(border,length=length(x)) }

    up = par("usr")
    ui = par("pin")

    ratx = (up[2]-up[1])/ui[1]
    raty=  (up[4]-up[3])/ui[2]

    for(i in 1:length(x))
      {
        x0 = x[i]
        y0 = y[i]

        usizx = siz[i]*ratx
        usizy = siz[i]*raty
###  print(i)
        A = prep.sym(sym, X=c(x0-usizx, x0+usizx),  Y=c(y0-usizy, y0+usizy),  delx=attr(sym, "DX"), dely=attr(sym, "DY"))
        attr(A, "OUT") = FALSE
        JREdraw(A)
      }
  }

`prep.sym` <-
function(KALL, X=c(0,1),  Y=c(0,1),  delx=c(0,1), dely=c(0,1))
  {
    ########  remove the mean value and scale everything from
    ########  -1 to 1  or -.5 to .5
    K = 0
    NALL = list()
     NKALL = length(KALL)
    if(NKALL<1) return
    for(i in 1:NKALL )
      {
        if(is.null(KALL[[i]])) next
       K = K+1
        NALL[[K]] = KALL[[i]]
        NALL[[K]]$x = RESCALE(KALL[[i]]$x, X[1], X[2], delx[1], delx[2])
        NALL[[K]]$y = RESCALE(KALL[[i]]$y, Y[1], Y[2], dely[1], dely[2])

      }

    attr(NALL, "DX") =  X
    attr(NALL, "DY") =  Y
     attr(NALL, "OUT") = TRUE
    return(NALL)

  }
`RESCALE` <-
function(x, nx1, nx2, minx, maxx)
{
  #    rescale a vector
  nx = nx1+(nx2-nx1)*(x-minx)/(maxx-minx)
  return(nx)
}

`JREdraw` <-
function(KALL)
  {

    if(attr(KALL,"OUT")==TRUE)
      {
        DX = attr(KALL, "DX")
        DY = attr(KALL, "DY")
        rect(DX[1], DY[1], DX[2] , DY[2], col=NA, border='black' )
      }

    NKALL = length(KALL)
    if(NKALL<1) return
    for(i in 1:NKALL )
      {if(is.null(KALL[[i]])) next

        if(KALL[[i]]$type=="lines")
          { lines(KALL[[i]]$x, KALL[[i]]$y, col=KALL[[i]]$col, lwd=KALL[[i]]$lwd) }

        if(KALL[[i]]$type=="points")
          { points(KALL[[i]]$x, KALL[[i]]$y, col=KALL[[i]]$col) }

        if(KALL[[i]]$type=="polygon")
          { polygon(KALL[[i]]$x, KALL[[i]]$y, col=KALL[[i]]$col, border=KALL[[i]]$border) }

        if(KALL[[i]]$type=="circles")
          { makecircles(KALL[[i]], col=KALL[[i]]$col, border=KALL[[i]]$border) }

        if(KALL[[i]]$type=="rects")
          { makerects(KALL[[i]], col=KALL[[i]]$col, border=KALL[[i]]$border,  lwd=KALL[[i]]$lwd) }

        if(KALL[[i]]$type=="arrows")
          { makearrows(KALL[[i]], col=KALL[[i]]$col,  lwd=KALL[[i]]$lwd) }

      }

  }



C = circle()

cex1 = 1.5

i=seq(from=0, to=-pi, length=180)
cx = cos(i);
cy = sin(i);
emptyplot()
title("Lower Hemisphere Projection View from Side")
####axis(1)
####axis(2)
lines(cx, cy)

i=seq(from=0, to=-pi/4, length=180)

cx = cos(i);
cy = sin(i);

polygon(c(0, cx, 0), c(0, cy, 0) , col=grey(0.85), border='black')

i=seq(from=-3*pi/4, to=-pi, length=180)

cx = cos(i);
cy = sin(i);

polygon(c(0, cx, 0), c(0, cy, 0) , col=grey(0.85), border='black')

segments(-0.25*cx[1], -0.25*cy[1], cx[1], cy[1], lwd=2, lty=2)
text(cx[1], cy[1], labels="Auxilliary Plane", srt=45, adj=c(1,0), font=2, xpd=TRUE, cex=cex1)

i=seq(from=0, to=-pi/4, length=180)

cx = cos(i);
cy = sin(i);

segments(-0.25*cx[180], -0.25*cy[180], cx[180], cy[180], lwd=2, lty=2)

segments(-0.25*cx[180], -0.25*cy[180], 0.25*cx[180], 0.25*cy[180], lwd=4, lty=1)

segments(-0.1*cx[180]+.01, -0.1*cy[180]+.01, 0.1*cx[180]+.01, 0.1*cy[180]+.01, lwd=2, lty=1)
segments(0.1*cx[180]+.01, 0.1*cy[180]+.01,0.1*cx[180]+.01, 0.1*cy[180]+.01+.05, lwd=2, lty=1)


segments(-0.1*cx[180]-0.01, -0.1*cy[180]-0.01, 0.1*cx[180]-0.01, 0.1*cy[180]-0.01, lwd=2, lty=1)

segments(-0.1*cx[180]-0.01, -0.1*cy[180]-0.01, -0.1*cx[180]-0.01, -0.1*cy[180]-0.01-.05, lwd=2, lty=1)


text(-0.25*cx[180], -0.25*cy[180] , labels="Fault Plane",  adj=c(1,0), font=2, xpd=TRUE, cex=cex1)
text(cx[180], cy[180] , labels="Fault Plane Projection",  adj=c(0,0), srt=-45, font=2, xpd=TRUE, cex=cex1)




abline(h=.35)
text(0,0.35, "Surface of Earth", pos=3, cex=cex1)

arrows(-1.05, 0, -1.05, -0.5)
text(-1.05, -0.25, "Depth", srt=-90, adj=c(1,-0.1), cex=cex1)

points(c(1, 0), c(0, -1), pch=20, cex=2, col=c("blue", "red") )

text(1,0, labels="T", pos=4, xpd=TRUE, font=2, col="blue", cex=cex1)

text(0,-1, labels="P", pos=1, xpd=TRUE, font=2, col="red", cex=cex1)

############  eye part
###  get.syms()

### data(EYEVIEW)
EYEVIEW2 = rotsym(EYEVIEW, -90, reflect=FALSE)

EYEVIEW2 = recentersym(EYEVIEW2)

draw.sym(0,0.8, EYEVIEW2, siz=.4)

MFOC2 = SDRfoc(0,45,-90, u=FALSE, ALIM=c(-1,-1, +1, +1), PLOT=FALSE)
    ###  Fcol2 = foc.color(foc.icolor(MFOC2$rake1), pal=1)
  Fcol2 =  grey(0.85)
    justfocXY( MFOC2, fcol = Fcol2, .9, .9 , size = c(.2,.2) )





###################################################
### code chunk number 4: focalmec.Rnw:395-399
###################################################
### par(asp=1)

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



###################################################
### code chunk number 5: focalmec.Rnw:438-443
###################################################

 mc = CONVERTSDR(65, 32, -34 )

printMEC(mc)



###################################################
### code chunk number 6: focalmec.Rnw:447-448
###################################################
plotMEC(mc, detail=2, up=FALSE)


###################################################
### code chunk number 7: focalmec.Rnw:457-472
###################################################
K = 10
rakes = runif(K, 30, 60)
dips = runif(K, 40, 55)
strikes = runif(K, 25,75)

 net(add=FALSE)

for(i in 1:K)
{
 mc = CONVERTSDR(floor(rakes[i]), floor( dips[i]),   floor(strikes[i]))
LP1 = lowplane( mc$M$az1-90, mc$M$d1, col='blue', UP=TRUE)
## LP2 = lowplane( mc$M$az2, mc$M$d2, col='red', UP=TRUE)
  U = focpoint(mc$V$az, mc$V$dip, col=6,  lab="", UP=TRUE)
}



###################################################
### code chunk number 8: faults1
###################################################
opar <- par(no.readonly = TRUE)

  par(MFA=c(1,3))
    par(mas=c(0.1, .1, 0.2, 0.1) )
    anim=0
############
###   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 , size = c(.4,.4) )
      mtext("Strike-slip fault", side = 3, line =-1.5)
#### box()


###################################################
### code chunk number 9: faults2
###################################################
############
    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 , size = c(.45,.45) )
       mtext("Normal fault", side = 3, line =-1.5)
#### box()


###################################################
### code chunk number 10: faults3
###################################################
############
    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 , size = c(.45,.45) )
    mtext("Reverse (Thrust) fault", side = 3, line =-1.5 )
####   box()
par(opar)



###################################################
### code chunk number 11: focArad
###################################################
`tomo.colors` <-
function(n, alpha = 1)
  {
    if ((n <- as.integer(n[1])) > 0)
      {
        k <- n%/%2
        h <- c(0/12, 2/12, 8/12)
        s <- c(1, 0, 1)
        v <- c(0.9, 0.9, 0.95)

        c(hsv(h = seq.int(h[1], h[2], length = k),
              s = seq.int(s[1], s[2], length = k),
              v = seq.int(v[1], v[2], length = k),
              alpha = alpha),
          hsv(h = seq.int(h[2], h[3], length = n - k + 1)[-1],
              s = seq.int(s[2], s[3], length = n -  k + 1)[-1],
              v = seq.int(v[2], v[3], length = n - k + 1)[-1],
              alpha = alpha))
      }

    else character(0)
  }
`Gcols` <-
function(plow=10, phi=10,  N=100, pal="rainbow", mingray=0.5)
{
  ###   get a palette with the upper or lower parts replaced
  if(missing(plow)) { plow = 10 }
  if(missing(phi)) { phi = 10 }
  if(missing(N)) { N = 100 }
  if(missing(pal)) { pal = "rainbow" }
  if(missing(mingray)) { mingray=0.5 }


  nlow = floor(plow*N/100)
  nhi = floor(phi*N/100)
  LOW = grey(seq(from=mingray, to =1, length=nlow))
  HI  = grey(seq(from=mingray, to =1, length=nhi))
  K = N-nlow-nhi

     FUN = match.fun(pal)
     Z = FUN(K)
  #####  Z = rainbow(K)
  return(c(LOW  , Z, HI) )
}



 mc  = CONVERTSDR(65, 32, -34 )
   
    MEC = MRake(mc$M)
    

     MEC$strike = mc$strike
    MEC$dipdir = mc$dipdir
    MEC$dip = mc$dip
    MEC$rake = mc$rake


    MEC$UP = FALSE
## 


mycol = Gcols(plow=0, phi=0,  N=100, pal="tomo.colors")

#######   graphics.off(); X11( width=15, height=7)

par(mfrow=c(1,3))
par(mai=c(.5,0.1,.5,0.1) )


radiateP(MEC, col=mycol )
text(0,1, labels="P-wave  radiation, Lower Hemisphere", pos=3, cex=1.5)

# dev.set(which = dev.next())
radiateSV(MEC, col=mycol )
text(0,1, labels="SV-wave  radiation, Lower Hemisphere", pos=3, cex=1.5)

# dev.set(which = dev.next())
radiateSH(MEC, col=mycol )
text(0,1, labels="SH-wave  radiation, Lower Hemisphere", pos=3, cex=1.5)




###################################################
### code chunk number 12: focalmec.Rnw:631-632
###################################################
`tomo.colors` <-
function(n, alpha = 1)
  {
    if ((n <- as.integer(n[1])) > 0)
      {
        k <- n%/%2
        h <- c(0/12, 2/12, 8/12)
        s <- c(1, 0, 1)
        v <- c(0.9, 0.9, 0.95)

        c(hsv(h = seq.int(h[1], h[2], length = k),
              s = seq.int(s[1], s[2], length = k),
              v = seq.int(v[1], v[2], length = k),
              alpha = alpha),
          hsv(h = seq.int(h[2], h[3], length = n - k + 1)[-1],
              s = seq.int(s[2], s[3], length = n -  k + 1)[-1],
              v = seq.int(v[2], v[3], length = n - k + 1)[-1],
              alpha = alpha))
      }

    else character(0)
  }
`Gcols` <-
function(plow=10, phi=10,  N=100, pal="rainbow", mingray=0.5)
{
  ###   get a palette with the upper or lower parts replaced
  if(missing(plow)) { plow = 10 }
  if(missing(phi)) { phi = 10 }
  if(missing(N)) { N = 100 }
  if(missing(pal)) { pal = "rainbow" }
  if(missing(mingray)) { mingray=0.5 }


  nlow = floor(plow*N/100)
  nhi = floor(phi*N/100)
  LOW = grey(seq(from=mingray, to =1, length=nlow))
  HI  = grey(seq(from=mingray, to =1, length=nhi))
  K = N-nlow-nhi

     FUN = match.fun(pal)
     Z = FUN(K)
  #####  Z = rainbow(K)
  return(c(LOW  , Z, HI) )
}



 mc  = CONVERTSDR(65, 32, -34 )
   
    MEC = MRake(mc$M)
    

     MEC$strike = mc$strike
    MEC$dipdir = mc$dipdir
    MEC$dip = mc$dip
    MEC$rake = mc$rake


    MEC$UP = FALSE
## 


mycol = Gcols(plow=0, phi=0,  N=100, pal="tomo.colors")

#######   graphics.off(); X11( width=15, height=7)

par(mfrow=c(1,3))
par(mai=c(.5,0.1,.5,0.1) )


radiateP(MEC, col=mycol )
text(0,1, labels="P-wave  radiation, Lower Hemisphere", pos=3, cex=1.5)

# dev.set(which = dev.next())
radiateSV(MEC, col=mycol )
text(0,1, labels="SV-wave  radiation, Lower Hemisphere", pos=3, cex=1.5)

# dev.set(which = dev.next())
radiateSH(MEC, col=mycol )
text(0,1, labels="SH-wave  radiation, Lower Hemisphere", pos=3, cex=1.5)




###################################################
### code chunk number 13: focmap1
###################################################

data(KAMCORN)

plot(range(KAMCORN$LON), range(KAMCORN$LAT), type='n', xlab="LON", ylab="LAT", asp=1) 
 
for(i in 1:length(KAMCORN$LAT) )
{
  mc  = CONVERTSDR(KAMCORN$STRIKE[i], KAMCORN$DIP[i], KAMCORN$RAKE[i] )
 
  MEC <- MRake(mc$M)
  MEC$UP = FALSE
  Fcol <- foc.color(foc.icolor(MEC$rake1), pal=1)

  justfocXY( MEC,  x=KAMCORN$LON[i],  y=KAMCORN$LAT[i], size = c(.5,.5), fcol = Fcol, xpd=FALSE )


}



###################################################
### code chunk number 14: focalmec.Rnw:681-682
###################################################

data(KAMCORN)

plot(range(KAMCORN$LON), range(KAMCORN$LAT), type='n', xlab="LON", ylab="LAT", asp=1) 
 
for(i in 1:length(KAMCORN$LAT) )
{
  mc  = CONVERTSDR(KAMCORN$STRIKE[i], KAMCORN$DIP[i], KAMCORN$RAKE[i] )
 
  MEC <- MRake(mc$M)
  MEC$UP = FALSE
  Fcol <- foc.color(foc.icolor(MEC$rake1), pal=1)

  justfocXY( MEC,  x=KAMCORN$LON[i],  y=KAMCORN$LAT[i], size = c(.5,.5), fcol = Fcol, xpd=FALSE )


}



###################################################
### code chunk number 15: focNIP
###################################################
plot(range(KAMCORN$LON), range(KAMCORN$LAT), type='n', xlab="LON", ylab="LAT", asp=1) 
 
for(i in 1:length(KAMCORN$LAT) )
{
  mc  = CONVERTSDR(KAMCORN$STRIKE[i], KAMCORN$DIP[i], KAMCORN$RAKE[i] )
 
  MEC <- MRake(mc$M)
  MEC$UP = FALSE
  Fcol <- foc.color(foc.icolor(MEC$rake1), pal=1)

  nipXY( MEC,  x=KAMCORN$LON[i],  y=KAMCORN$LAT[i], size = c(.5,.5), fcol = Fcol )


}





###################################################
### code chunk number 16: focalmec.Rnw:718-719
###################################################
plot(range(KAMCORN$LON), range(KAMCORN$LAT), type='n', xlab="LON", ylab="LAT", asp=1) 
 
for(i in 1:length(KAMCORN$LAT) )
{
  mc  = CONVERTSDR(KAMCORN$STRIKE[i], KAMCORN$DIP[i], KAMCORN$RAKE[i] )
 
  MEC <- MRake(mc$M)
  MEC$UP = FALSE
  Fcol <- foc.color(foc.icolor(MEC$rake1), pal=1)

  nipXY( MEC,  x=KAMCORN$LON[i],  y=KAMCORN$LAT[i], size = c(.5,.5), fcol = Fcol )


}





###################################################
### code chunk number 17: ALLPT
###################################################

 net()
  PZZ =  focpoint(PKAM$Paz, PKAM$Pdip, col='red',  pch=3, lab="", UP=FALSE)
 TZZ =  focpoint(PKAM$Taz, PKAM$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))



###################################################
### code chunk number 18: focalmec.Rnw:756-757
###################################################

 net()
  PZZ =  focpoint(PKAM$Paz, PKAM$Pdip, col='red',  pch=3, lab="", UP=FALSE)
 TZZ =  focpoint(PKAM$Taz, PKAM$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))



###################################################
### code chunk number 19: TERN1
###################################################

fcols = foc.color(PKAM$IFcol,1)

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



 MFOC3 = SDRfoc(135,45,90, u=FALSE, ALIM=c(-1,-1, +1, +1), PLOT=FALSE)
   
 Fcol3 = foc.color(foc.icolor(MFOC3$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)
  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( MFOC3, fcol = Fcol3, 1.2, -0.9, size = c(.1,.1) )
justfocXY( MFOC2, fcol = Fcol2, -1.2, -0.9, size = c(.1,.1) )
justfocXY( MFOC1, fcol = Fcol1, 0, 1.414443+.2, size = c(.1,.1) )
## mtext("Ternary Plot of focal mecahnisms", side = 1, line = 0 , font=2,  xpd=TRUE )
 


###################################################
### code chunk number 20: focalmec.Rnw:804-805
###################################################

fcols = foc.color(PKAM$IFcol,1)

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



 MFOC3 = SDRfoc(135,45,90, u=FALSE, ALIM=c(-1,-1, +1, +1), PLOT=FALSE)
   
 Fcol3 = foc.color(foc.icolor(MFOC3$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)
  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( MFOC3, fcol = Fcol3, 1.2, -0.9, size = c(.1,.1) )
justfocXY( MFOC2, fcol = Fcol2, -1.2, -0.9, size = c(.1,.1) )
justfocXY( MFOC1, fcol = Fcol1, 0, 1.414443+.2, size = c(.1,.1) )
## mtext("Ternary Plot of focal mecahnisms", side = 1, line = 0 , font=2,  xpd=TRUE )
 


###################################################
### code chunk number 21: CONTRDPT
###################################################


 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) )
opar <- par(no.readonly = TRUE)

 par(mfrow=c(1,3))

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

    image(KP$x, KP$y, KP$z, add=TRUE, col=terrain.colors(100))

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

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

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


    image(KT$x, KT$y, KT$z, add=TRUE, col=terrain.colors(100))
    
    antipolygon(CC$x,CC$y,col="white")
    net(add=1)
    focpoint(PKAM$Taz, PKAM$Tdip, col='blue',  pch=4, lab="", UP=FALSE)
    text(0,1.04,labels="T-axes 2D Density", font=2, cex=1.2)

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

    image(KP$x, KP$y, KP$z, add=TRUE, col=terrain.colors(100))


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

    net(add=1)


    contour(KT$x, KT$y, KT$z, add=TRUE, lwd=1.2)
    

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



###################################################
### code chunk number 22: focalmec.Rnw:866-867
###################################################


 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) )
opar <- par(no.readonly = TRUE)

 par(mfrow=c(1,3))

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

    image(KP$x, KP$y, KP$z, add=TRUE, col=terrain.colors(100))

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

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

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


    image(KT$x, KT$y, KT$z, add=TRUE, col=terrain.colors(100))
    
    antipolygon(CC$x,CC$y,col="white")
    net(add=1)
    focpoint(PKAM$Taz, PKAM$Tdip, col='blue',  pch=4, lab="", UP=FALSE)
    text(0,1.04,labels="T-axes 2D Density", font=2, cex=1.2)

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

    image(KP$x, KP$y, KP$z, add=TRUE, col=terrain.colors(100))


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

    net(add=1)


    contour(KT$x, KT$y, KT$z, add=TRUE, lwd=1.2)
    

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



###################################################
### code chunk number 23: ternmapkam1
###################################################
x = fmod(PKAM$LONS, 360)
y = PKAM$LATS
plot(x, y, asp=1, type="p", xlab="LON", ylab="LAT" )

 u = par("usr")

    RI = RectDense( x, y, icut=5, u=u, ndivs=10)
   rect(RI$icorns[,1],RI$icorns[,2],RI$icorns[,3],RI$icorns[,4], col=NA, border='blue')






###################################################
### code chunk number 24: focalmec.Rnw:902-903
###################################################
x = fmod(PKAM$LONS, 360)
y = PKAM$LATS
plot(x, y, asp=1, type="p", xlab="LON", ylab="LAT" )

 u = par("usr")

    RI = RectDense( x, y, icut=5, u=u, ndivs=10)
   rect(RI$icorns[,1],RI$icorns[,2],RI$icorns[,3],RI$icorns[,4], col=NA, border='blue')






###################################################
### code chunk number 25: ternmapkam2
###################################################


Fcol = foc.color(PKAM$IFcol, pal=1)
i = 1
sizy = RI$icorns[i,4]-RI$icorns[i,2]
sizx = RI$icorns[i,3]-RI$icorns[i,1]
siz = .5*min(c(sizy, sizx))


plot(x, y, asp=1, type="p" , xlab="LON", ylab="LAT" )

 u = par("usr")

    RI = RectDense( x, y, icut=5, u=u, ndivs=10)
   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 = x>RI$icorns[i,1]& y>RI$icorns[i,2] & x<RI$icorns[i,3] & y<RI$icorns[i,4]
        jh =PKAM$h[flag]
        jv= PKAM$v[flag]
        PlotTernfoc(jh,jv,x=mean(RI$icorns[i,c(1,3)]), y=mean(RI$icorns[i,c(2,4)]), siz=siz, fcols=Fcol[flag], add=TRUE)
      }




###################################################
### code chunk number 26: focalmec.Rnw:943-944
###################################################


Fcol = foc.color(PKAM$IFcol, pal=1)
i = 1
sizy = RI$icorns[i,4]-RI$icorns[i,2]
sizx = RI$icorns[i,3]-RI$icorns[i,1]
siz = .5*min(c(sizy, sizx))


plot(x, y, asp=1, type="p" , xlab="LON", ylab="LAT" )

 u = par("usr")

    RI = RectDense( x, y, icut=5, u=u, ndivs=10)
   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 = x>RI$icorns[i,1]& y>RI$icorns[i,2] & x<RI$icorns[i,3] & y<RI$icorns[i,4]
        jh =PKAM$h[flag]
        jv= PKAM$v[flag]
        PlotTernfoc(jh,jv,x=mean(RI$icorns[i,c(1,3)]), y=mean(RI$icorns[i,c(2,4)]), siz=siz, fcols=Fcol[flag], add=TRUE)
      }




###################################################
### code chunk number 27: smomapkam2
###################################################
##   plot(x, y, asp=1, type="p" , xlab="LON", ylab="LAT" )
x = fmod(PKAM$LONS, 360)
y = PKAM$LATS
plot(x, y, asp=1, type="p", xlab="LON", ylab="LAT" )

 u = par("usr")
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 = x>RI$icorns[i,1]& y>RI$icorns[i,2] & x<RI$icorns[i,3] & y<RI$icorns[i,4]
        paz=PKAM$Paz[flag]
        pdip=PKAM$Pdip[flag]
         taz=PKAM$Taz[flag]
        tdip=PKAM$Tdip[flag]
        ax=mean(RI$icorns[i,c(1,3)])
        ay=mean(RI$icorns[i,c(2,4)])
        siz=(RI$icorns[1,3]-RI$icorns[1,1])/2.5

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

######dev.set(2)
######pnet(MN, add=FALSE)
  
    ALPH = alpha95(paz, pdip)
        n = length( paz)
      KPspat[i,] =   c(ax, ay, n, ALPH$Ir,  ALPH$Dr, ALPH$R, ALPH$K, ALPH$S, ALPH$Alph95 , ALPH$Kappa)

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

        

      }

   
   



###################################################
### code chunk number 28: focalmec.Rnw:1001-1002
###################################################
##   plot(x, y, asp=1, type="p" , xlab="LON", ylab="LAT" )
x = fmod(PKAM$LONS, 360)
y = PKAM$LATS
plot(x, y, asp=1, type="p", xlab="LON", ylab="LAT" )

 u = par("usr")
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 = x>RI$icorns[i,1]& y>RI$icorns[i,2] & x<RI$icorns[i,3] & y<RI$icorns[i,4]
        paz=PKAM$Paz[flag]
        pdip=PKAM$Pdip[flag]
         taz=PKAM$Taz[flag]
        tdip=PKAM$Tdip[flag]
        ax=mean(RI$icorns[i,c(1,3)])
        ay=mean(RI$icorns[i,c(2,4)])
        siz=(RI$icorns[1,3]-RI$icorns[1,1])/2.5

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

######dev.set(2)
######pnet(MN, add=FALSE)
  
    ALPH = alpha95(paz, pdip)
        n = length( paz)
      KPspat[i,] =   c(ax, ay, n, ALPH$Ir,  ALPH$Dr, ALPH$R, ALPH$K, ALPH$S, ALPH$Alph95 , ALPH$Kappa)

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

        

      }

   
   



