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

```ReadND <-
function(filename,verbose=TRUE){

# turn off warnings--coercion creates NAs, which is normal here
options(warn=-1)

vmod=-1
dz=NULL
dname=NULL
model=-1

pubyear=NULL
comment=FALSE

fid=file(filename,'r')
vmodini=rep(NaN,6)
vmod=NULL

parmcnt=1
namecnt=1
linecnt=1
done=0
while( done==0){
if( length(line)==0 ){
done=1
}else if(comment==TRUE & !('/' %in% strsplit(line,'*',fixed=TRUE)[[1]][2:length(strsplit(line,'*',fixed=TRUE)[[1]])]) ){
}else if(comment==TRUE){
comment=FALSE
}else{
if(length(strsplit(line,'/*',fixed=TRUE)[[1]])>1){
comment=TRUE # ignore lines until */ is reached
}
line=strsplit(line,'/*',fixed=TRUE)[[1]][1]
## strip // and # comments
line=strsplit(line,'//',fixed=TRUE)[[1]][1]
line=strsplit(line,'#',fixed=TRUE)[[1]][1]

if(nchar(line)>0 & !is.na(line)){
vec = strsplit(line,' ')[[1]]
if(is.na(sum(as.numeric(vec[nzchar(vec)])))){
# line includes letters
nws=paste(strsplit(line,' ')[[1]],collapse='')
if(substr(nws,1,1)=='!'){

vec = strsplit(line,' ')[[1]]
keyword = tolower(vec[which(nzchar(vec))[1]])
par = vec[which(nzchar(vec))[2]]
}else if(keyword=='!year'){
pubyear=as.numeric(par)
}else if(keyword=='!name'){
mname=par
}else{
if(verbose){
}
}
}else{
vec = strsplit(line,' ')[[1]]
w = which(nzchar(vec))
name = paste(vec[min(w):max(w)],collapse = ' ')
dz[namecnt]=vmod[parmcnt-1,1]
dname=c(dname,name)
if(verbose){
}
namecnt=namecnt+1
}
}else{

pcnt=1
done2=0
vmod=rbind(vmod,vmodini,deparse.level=0)

vec = strsplit(line,' ')[[1]]
vmod[parmcnt,1:sum(nzchar(vec))] = as.numeric(vec[nzchar(vec)])
parmcnt=parmcnt+1
}
}
linecnt=linecnt+1
}
}

close(fid)

##### replace all -1 by NaN
indies=which(vmod==-1)
vmod[indies]=NaN
indies=which(dz==-1)
dz[indies]=-1

}

model=list()
model\$z=vmod[,1]
model\$vp=vmod[,2]
model\$vs=vmod[,3]
model\$rho=vmod[,4]
model\$qp=vmod[,5]
model\$qs=vmod[,6]
model\$name=mname
model\$year=pubyear

model\$conr=NaN
model\$moho=NaN
model\$d410=NaN
model\$d520=NaN
model\$d660=NaN
model\$cmb=NaN
model\$icb=NaN
model\$dz=NULL
model\$dname=NULL
anz=length(dz)
cnt=1

for( indy in 1:anz){
vec = strsplit(dname[indy],' ')[[1]]
w = which(nzchar(vec))
p=tolower(paste(vec[min(w):max(w)],collapse=' '))
model\$conr=dz[indy]
}else if(p %in% c('moho','mantle')){
model\$moho=dz[indy]
}else if(p %in% c('olivine alpha beta','transition zone')){
model\$d410=dz[indy]
}else if(p == 'olivine beta gamma'){
model\$d520=dz[indy]
}else if(p %in% c('olivine gamma perovskite','lower mantle')){
model\$d660=dz[indy]
}else if(p %in% c('cmb','outer core','outer-core')){
model\$cmb=dz[indy]
}else if(p %in% c('icb','inner core','inner-core')){
model\$icb=dz[indy]
}else{
if(verbose){
}
model\$dz=c(model\$dz, dz[indy])
model\$dname=c(model\$dname,p)
}
}

if(verbose){