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

#### Documented in FindPrange

FindPrange <-function(phase,imodel,h,dangle){
angles=NULL

phasev=strsplit(phase,'')[[1]]

if(!is.na(imodel\$cmb) && h > imodel\$cmb){stop('Focal depth must be above the CMB')}
if(!is.na(imodel\$cmb) && h == imodel\$cmb && phasev[2] %in% c('c', 'K')){stop("Focal depth must be strictly above CMB for this phase")}

if( phasev[1] %in% c('p','s') ) {

if(h == 0){
stop('Depth phases do not exist for h = 0')
}
depthleg=phasev[1]
surfaceleg=phasev[2:length(phasev)]

aa=FindPrange(paste(surfaceleg,collapse=''),imodel,0,10)
angles=aa[[1]]
minangle=aa[[2]]
maxangle=aa[[3]]

raypminangle=ConvAng2p(paste(surfaceleg, collapse = ''),0,minangle,model = NULL, imodel\$vp[1],imodel\$vs[1],imodel\$rp)
raypmaxangle=ConvAng2p(paste(surfaceleg, collapse = ''),0,maxangle,model = NULL, imodel\$vp[1],imodel\$vs[1],imodel\$rp)

focus=InterpModel(imodel,h,'simple')

minangle=ConvP2Ang(depthleg,h,raypminangle,vp=focus\$vp[1],vs=focus\$vs[1],rp=focus\$rp)
maxangle=ConvP2Ang(depthleg,h,raypmaxangle,vp=focus\$vp[1],vs=focus\$vs[1],rp=focus\$rp)

minangle=minangle[!(minangle<90)]
maxangle=maxangle[!(maxangle<90)]

if( sum(is.na(minangle)) ){minangle=90}
if( sum(is.na(maxangle)) ){maxangle=90}

if( minangle>maxangle ){
swapper=minangle
minangle=maxangle
maxangle=swapper
}

if( is.na(minangle+dangle+maxangle) ){
angles=NULL
}else{
angles=seq(from = minangle, by = dangle, to = maxangle)
if( angles[length(angles)]!=maxangle ){angles=c(angles,maxangle)}
}

return(list(angles=angles,minangle=minangle,maxangle=maxangle))

}

if(!('criticalrays' %in% names(imodel)) ){

imodel=ImproveModel(imodel)[[1]]
}

imodel\$criticalrays\$z = round(imodel\$criticalrays\$z, 3)

focus=InterpModel(imodel,h,'simple')

if( prod(focus\$vs == 0) ){ #5/16/11: changed from focus\$vs == 0 to prod(focus\$vs) == 0 because it produced warnings on discontinuities.  used prod instead of sum to permit sources on the CMB but not inside the outer core.
angles=NULL
return(list(angles=angles,minangle=minangle,maxangle=maxangle))
}

psurf=imodel\$rp/imodel\$vp[1]
ssurf=imodel\$rp/imodel\$vs[1]
psurf=psurf * pi/180
ssurf=ssurf * pi/180
psurfangle=ConvP2Ang('P',0,psurf,vp=imodel\$vp[1],vs=imodel\$vs[1],rp=imodel\$rp)
ssurfangle=ConvP2Ang('S',0,ssurf,vp=imodel\$vp[1],vs=imodel\$vs[1],rp=imodel\$rp)

psteepangle=0.1
ssteepangle=0.1

aa=ConvVdepth2p(imodel,imodel\$z[length(imodel\$z)-1])
psteep=aa[[1]]
ssteep=aa[[2]]
psteepangle=min(ConvP2Ang('P',h,psteep,vp=focus\$vp[1],vs=focus\$vs[1],rp=imodel\$rp));
ssteepangle=min(ConvP2Ang('S',h,ssteep,vp=focus\$vp[1],vs=focus\$vs[1],rp=imodel\$rp));

if( !is.na(imodel\$cmb) ){

indies=which(0 == round(imodel\$criticalrays\$z - imodel\$cmb, 3))
pcmb=imodel\$criticalrays\$p[indies]
scmb=imodel\$criticalrays\$s[indies]

indies=which(!is.na(pcmb))
pcmb=pcmb[indies[1]]
indies=which(!is.na(scmb))
scmb=scmb[indies[1]]

if( h<=imodel\$cmb ){

vp=focus\$vp[length(focus\$vp)]
vs=focus\$vs[length(focus\$vs)]
}else{

vp=focus\$vp[1]
vs=focus\$vs[1]
}

pcmbangle=ConvP2Ang('P',h,pcmb,vp=vp,vs=vs,rp=imodel\$rp);
scmbangle=ConvP2Ang('S',h,scmb,vp=vp,vs=vs,rp=imodel\$rp);

if( h<=imodel\$cmb ){

pcmbangle=pcmbangle[length(pcmbangle)]
scmbangle=scmbangle[length(scmbangle)]
}else{

pcmbangle=pcmbangle[1]
scmbangle=scmbangle[1]
}

}else{

pcmb=psteep
scmb=ssteep
pcmbangle=psteepangle
scmbangle=ssteepangle
}

if( !is.na(imodel\$icb) ){

indies=which(0 == round(imodel\$criticalrays\$z - imodel\$icb, 3))
picb=imodel\$criticalrays\$p[indies]
sicb=imodel\$criticalrays\$s[indies]

indies=which(!is.na(picb))
picb=picb[indies[1]]
indies=which(!is.na(sicb))
sicb=sicb[indies[1]]

if( h<=imodel\$icb ){

vp=focus\$vp[length(focus\$vp)]
vs=focus\$vs[length(focus\$vs)]
}else{

vp=focus\$vp[1]
vs=focus\$vs[1]
}

picbangle=ConvP2Ang('P',h,picb,vp=vp,vs=vs,rp=imodel\$rp)
sicbangle=ConvP2Ang('S',h,sicb,vp=vp,vs=vs,rp=imodel\$rp)
skicbangle=ConvP2Ang('S',h,picb,vp=vp,vs=vs,rp=imodel\$rp)

if( h<=imodel\$cmb ){

picbangle=picbangle[length(picbangle)]
sicbangle=sicbangle[length(sicbangle)]
skicbangle=skicbangle[length(skicbangle)]
}else{

picbangle=picbangle[1]
sicbangle=sicbangle[1]
skicbangle=skicbangle[1]
}

}else{

picb=psteep
sicb=ssteep
picbangle=psteepangle
sicbangle=ssteepangle
skicbangle=psteepangle
}

if( !is.na(imodel\$cmb) ){

cmbindies=which(imodel\$z==imodel\$cmb)

vpcmb=imodel\$vp[cmbindies]

cmbreflectp=(imodel\$rp-imodel\$cmb)/vpcmb[2]

cmbreflectp=cmbreflectp * pi/180

pcmbreflectangle=ConvP2Ang('P',h,cmbreflectp,vp=vp,vs=vs,rp=imodel\$rp)
scmbreflectangle=ConvP2Ang('S',h,cmbreflectp,vp=vp,vs=vs,rp=imodel\$rp)

if( h<=imodel\$cmb ){

pcmbreflectangle=pcmbreflectangle[length(pcmbreflectangle)]
scmbreflectangle=scmbreflectangle[length(scmbreflectangle)];
}else{

pcmbreflectangle=pcmbreflectangle[1]
scmbreflectangle=scmbreflectangle[1]
}

}else{

scmbreflectangle=NaN
pcmbreflectangle=NaN
cmbreflectp=NaN
}

if( !is.na(imodel\$icb) ){

icbindies=which(imodel\$z==imodel\$icb)

vpicb=imodel\$vp[icbindies[2]]
vsicb=imodel\$vs[icbindies[2]]

picbshallowp=(imodel\$rp-imodel\$icb)/vpicb
sicbshallowp=(imodel\$rp-imodel\$icb)/vsicb

picbshallowp=picbshallowp * pi/180
sicbshallowp=sicbshallowp * pi/180

picbshallowangle=ConvP2Ang('P',h,picbshallowp,vp=vp,vs=vs,rp=imodel\$rp)
sicbshallowangle=ConvP2Ang('S',h,sicbshallowp,vp=vp,vs=vs,rp=imodel\$rp)
skicbshallowangle=ConvP2Ang('S',h,picbshallowp,vp=vp,vs=vs,rp=imodel\$rp)

if( h<=imodel\$icb ){

picbshallowangle=picbshallowangle[length(picbshallowangle)];
sicbshallowangle=sicbshallowangle[length(sicbshallowangle)]
skicbshallowangle=skicbshallowangle[length(skicbshallowangle)]
}else{

picbshallowangle=picbshallowangle[1]
sicbshallowangle=sicbshallowangle[1]
skicbshallowangle=skicbshallowangle[1]
}

}else{

sicbshallowangle=NaN
picbshallowangle=NaN
skicbshallowangle=NaN
skicbshallowp=NaN
picbshallowp=NaN
sicbshallowp=NaN
}

aa=StripRepetitions(phase)
phase=aa[[1]]
repetitions=aa[[2]]

if(phase == 'P'){
minangle=pcmbangle
if(is.na(minangle)){minangle=psteepangle}
if(h>0){
maxangle=180
}else{
maxangle=90
}
}else if(phase %in% c('PP','PPP')){
minangle=pcmbangle
if(is.na(minangle)){
minangle=psteepangle
}
maxangle=90
}else if(phase %in% c('PcP','PcPPcP')){
minangle=psteepangle
maxangle=pcmbangle
}else if(phase %in% c('PKP','PKKP','PKPPKP')){
minangle=psteepangle
maxangle=min(c(pcmbreflectangle,pcmbangle))
}else if(phase == 'PKiKP'){
minangle=psteepangle
maxangle=picbangle
}else if(phase %in% c('PKIKP','PKIKPPKIKP')){
minangle=psteepangle
maxangle=picbshallowangle
}else if(phase == 'S'){
minangle=scmbangle
if(is.na(minangle)){
minangle=ssteepangle
}
if( h>0){
maxangle=180
}else{
maxangle=90
}
}else if(phase %in% c('SS','SSS')){
minangle=scmbangle
if( is.na(minangle)){

minangle=ssteepangle
}
maxangle=90
}else if(phase %in% c('ScS','ScSScS')){
minangle=ssteepangle
if(is.na(imodel\$icb)){

minangle=psteepangle
}else{
minangle=ssteepangle
}
maxangle=scmbangle
}else if(phase %in% c('SKS','SKKS','SKSSKS')){
if( imodel\$vs[length(imodel\$vs)-1]==0){

minangle=psteepangle
}else{

minangle=ssteepangle
}
maxangle=min(c(scmbreflectangle,scmbangle))
}else if(phase == 'SKiKS'){
minangle=ssteepangle
maxangle=skicbangle
}else if(phase == 'SKIKS'){
minangle=ssteepangle
maxangle=skicbshallowangle
}else if(phase == 'PS'){
if(!is.na(imodel\$cmb)){

minangle=ConvP2Ang('P',h,scmb,vp=vp,vs=vs,rp=imodel\$rp)
}else{

minangle=ConvP2Ang('P',h,ssteep,vp=focus\$vp[length(focus\$vp)],vs=focus\$vs[length(focus\$vs)],rp=imodel\$rp[length(imodel\$rp)])
}
minangle=minangle[length(minangle)]
if( is.na(minangle)){
minangle=psteepangle
}
maxangle=min(c(psurfangle, 90))
}else if(phase == 'SP'){
minangle=scmbangle
if(is.na(minangle)){
minangle=ssteepangle
}
maxangle=min(c(psurfangle, 90))
}else if(phase %in% c('PcS','PcSPcS')){
minangle=psteepangle
maxangle=pcmbangle
}else if(phase %in% c('ScP','ScPScP')){
minangle=psteepangle
maxangle=scmbangle
}else{
print(paste('MKSMARTTAKEOFF: unknown phase ', phase, ', using default.'))
minangle=psteepangle
if( h>0){
maxangle=180
}else{
maxangle=90
}
}

if( is.na(minangle+dangle+maxangle)){
angles=NULL
}else{
angles=seq(from=minangle,by=dangle,to=maxangle)
if( angles[length(angles)]!=maxangle){
angles=c(angles, maxangle)
}
}

return(list(angles=angles,minangle=minangle,maxangle=maxangle))
}

## 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.