# R/OptimizeDist.R In TauP.R: Earthquake Traveltime Calculations for 1-D Earth Models

#### Documented in OptimizeDist

```OptimizeDist <- function(alphalimit,deltalimit,phase,h,imodel){

### Define internal function:
ConvAng2pWrap=function(phase,h,angle,vp,vs,rp){

if( length(vp)==2){
if( angle<=90){
whichvelocity=2
}else{
whichvelocity=1
}
}else{
whichvelocity=1
}

p=ConvAng2p(phase,h,angle,model=NULL,vp=vp[whichvelocity],vs=vs[whichvelocity],rp=rp)
return(p)
}

QualifyExtremum=function(a,d){ # is the center of 3 points a local min, max, or neither
sorter=order(a)
a=a[sorter]
d=d[sorter]

if( (d[2] < d[1]) && (d[2] < d[3])){
qualifier=-1 # local minimum
}else{
if( (d[2] > d[1]) && (d[2] > d[3])){
qualifier=1 # local maximum
}else{
qualifier=0 # neither
}
}
return(qualifier)
}

### Done defining internal functions

tt=proc.time()
eps = .Machine\$double.eps
extremalpha=NULL
extremp=NULL
extremdelta=NULL

sorter=order(alphalimit)
alphalimit=alphalimit[sorter]
deltalimit=deltalimit[sorter]

focus=InterpModel(imodel,h,'simple')
indy=which(focus\$z==h)
vp=focus\$vp[indy]
vs=focus\$vs[indy]
identepsilon=1e-6 # changed from 1e-8 because sometimes diff(deltalimit)==0 while diff(alphalimit) > 1e-8

if( (diff(alphalimit) < identepsilon) && (diff(deltalimit) < identepsilon) ){
extremalpha=alphalimit[1]
extremdelta=deltalimit[1]
extremp=ConvAng2pWrap(phase, h, extremalpha, vp, vs, imodel\$rp)
return(list(extremalpha=extremalpha,extremp=extremp,extremdelta=extremdelta))
}
phi=(1+sqrt(5))/2
goldrat=phi-1
goldepsilon=eps^.25

goldmaxcnt=ceiling((log(alphalimit[length(alphalimit)]-log(alphalimit[1]))-log(goldepsilon))/log(phi))
goldmaxcnt=goldmaxcnt+1
alphalist=alphalimit
deltalist=deltalimit
alphalen=length(alphalist)
rayplist=alphalist * 0
for(i in 1:length(alphalist)){ # This for loop was added 5/16/11--previously it just had rayplist = CongAng2pWrap(phase, h, alphalist, vp, vs, imodel\$rp), which caused a warning.
rayplist[i] = ConvAng2pWrap(phase, h, alphalist[i], vp, vs, imodel\$rp)
}
a=NULL
a[1]=alphalimit[1]
a[3]=alphalimit[2]
a[2]=a[1]+goldrat*(a[3]-a[1])
p=NULL
p[1]=rayplist[1]
p[3]=rayplist[2]
p[2]=ConvAng2pWrap(phase,h,a[2],vp,vs,imodel\$rp)
d=NULL
d[1]=deltalimit[1]
d[3]=deltalimit[2]
d[2]=FindDist4p(phase,h,imodel,p=p[2],takeoff=a[2])[[1]]
extremumfound=0
done=FALSE
cnt=0
while(!done){

alphalist=c(alphalist, a[2])
deltalist=c(deltalist, d[2])
rayplist=c(rayplist, p[2])

if(sum(is.na(c(p,a,d)))){return(list(NaN,NaN,NaN))}

qualifier=QualifyExtremum(a,d)

if(qualifier == -1){
searchwhat=qualifier
done=TRUE
}else if(qualifier == 1){
searchwhat=qualifier
done=TRUE
}else if(qualifier == 0){
searchwhat=0
if (cnt==0){
a[2]=a[3]-goldrat*abs(a[3]-a[1])
p[2]=ConvAng2pWrap(phase,h,a[2],vp,vs,imodel\$rp)
d[2]=FindDist4p(phase,h,imodel,p=p[2],takeoff=a[2])[[1]]

alphalist=c(alphalist,a[2])
deltalist=c(deltalist,d[2])
rayplist=c(rayplist,p[2])
}else{

intervallength=sqrt((a[3]-a[1])^2+(d[3]-d[1])^2) # JFA can't imagine why one would use the 2-D distance here, but it shouldn't matter and if it aint broke, don't fix it...5/18/11

prightofa=ConvAng2pWrap(phase,h,rightofa,vp,vs,imodel\$rp)
pleftofc=ConvAng2pWrap(phase,h,leftofc,vp,vs,imodel\$rp)
aa=FindDist4p(phase,h,imodel,p=c(prightofa,pleftofc),takeoff=c(rightofa,leftofc))\$dist
frightofa=aa[1]
fleftofc=aa[2]
if(is.na(frightofa)|is.na(fleftofc)){return(list(NaN,NaN,NaN))}

alphalist=c(alphalist, rightofa, leftofc)
deltalist=c(deltalist, frightofa,fleftofc)
rayplist=c(rayplist, prightofa, pleftofc)
fdiscr=sign(d[1]-d[3])

if(fdiscr == -1){
searchwhat = 1
extremalpha=a[3]
extremp=p[3]
extremdist=d[3]
done=1
extremumfound=1
}else{
searchwhat=+1
a[2]=rightofa
p[2]=prightofa
d[2]=frightofa
}
}else{

searchwhat=-1
a[2]=rightofa
p[2]=prightofa
d[2]=frightofa
done=1
}else{
return(list(extremalpha=NaN,extremdelta=NaN,extremp=NaN))
}
}

}else if(fdiscr == 0){

searchwhat=-1
a[2]=rightofa
p[2]=prightofa
d[2]=frightofa
done=1
stop('OptimizeDist: insufficient information!')
searchwhat=+1
a[2]=rightofa
p[2]=prightofa
d[2]=frightofa
done=1
}else{
stop('OptimizeDist: unexpected gradient sign (1)!')
}
}else if(fdiscr==1){

searchwhat=-1
a[2]=leftofc
p[2]=pleftofc
d[2]=fleftofc
done=1
}else{
searchwhat=-1
extremalpha=a[3]
extremp=p[3]
extremdist=d[3]
done=1
extremumfound=1
}

searchwhat=-1
extremalpha=a[1]
extremp=p[1]
extremdist=d[1]
done=1
extremumfound=1
searchwhat= 1
a[2]=rightofa
p[2]=prightofa
d[2]=frightofa
done=1
}else{
stop('OptimizeDist: unexpected gradient sign (2)!')
}
}else{
warning(paste('OptimizeDist: unexpected exception. Gradients: A: ', (grada), ', C: ',(gradc),' (Delta=', (deltalimit[1]), '...', (deltalimit[2]), ')'))
}
}
}

cnt=cnt+1

done=(cnt>2) | done
extremumfound=(cnt>2) | extremumfound
}
sorter=order(a)
a=a[sorter]
p=p[sorter]
d=d[sorter]
if( extremumfound==0){
done=FALSE
cnt=0
while(!done){
intlen=abs(diff(a))
if( sum(intlen)>goldepsilon){
if( intlen[1]>=intlen[2]){
newa=a[1]+goldrat*intlen[1]
}else{
newa=a[3]-goldrat*intlen[2]
}

newp=ConvAng2pWrap(phase,h,newa,vp,vs,imodel\$rp)
newd=FindDist4p(phase,h,imodel,p=newp,takeoff=newa)[[1]]

alphalist=c(alphalist, newa)
rayplist=c(rayplist, newp)
deltalist=c(deltalist,newd)

a=c(a, newa)
p=c(p, newp)
d=c(d, newd)

sorter=order(a)
a=a[sorter]
p=p[sorter]
d=d[sorter]

if(searchwhat== -1){
extremindy=which(min(d)==d)[1]
}else if(searchwhat == 1){
extremindy=which(max(d)==d)[1]
}else{
stop(paste('OptimizeDist: unexpected SEARCHWHAT: ', (searchwhat)))
}
if(extremindy==1){
newpoints=1:3
}else if(extremindy==4){
newpoints=2:4
}else{
newpoints=extremindy + (-1):1
}
a=a[newpoints]
p=p[newpoints]
d=d[newpoints]

}else{
done=TRUE
}

cnt=cnt+1
if( cnt>goldmaxcnt){
done=TRUE
print('OptimizeDist: iteration exhaust!')
}
}
}
if(searchwhat==-1){
extremdelta=min(deltalist)
indy=(deltalist==extremdelta)
}else if(searchwhat==1){
extremdelta=max(deltalist)
indy=(deltalist==extremdelta)
}else if(searchwhat==0){
indy=NULL
}else{
stop('OptimizeDist: unexpected SEARCHWHAT value.')
}
return(list(extremalpha=(alphalist[indy])[1],extremp=(rayplist[indy])[1],extremdelta=extremdelta))
}
```

## Try the TauP.R package in your browser

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

TauP.R documentation built on May 2, 2019, 3:25 a.m.