# redpuma.nawk modified from earlier programs by Jim Lawson of the Oklahoma
# Geological Survey on 2001 MAR29. Oklahoma Geological Survey Observatory
# programs are Public Domain.
BEGIN{
# Change ~/jims/nawk/redpuma.nawk to your own path to redpuma.nawk
# redpuma_ami_list.html is what you typically get from
# seismo.ethz.ch/redpuma/redpuma_ami_list.html
# if you use "save as source" on your browser.
# you may change output file redpuma to another name, and include a
# path if you want it in other than your current directory
# the NAWK interpreter is assumed to be in your path. enter "which nawk"
# to check this. If it is not in your path, add its location, eg.
# /usr/local/bin/nawk. Or you may put the full path in the usage statement
# in place of "nawk". GNU "gawk" may also be used.
print "usage: nawk -f ~/jims/nawk/redpuma.nawk redpuma_ami_list.html > redpuma","\n"
#print "usage: nawk|gawk -f yourpath/redpuma.nawk redpuma_ami_list.html > redpum#a","\n" Run from the directory which includes redpuma_ami_list.html to
#write in the same directory the file "redpuma"
#You could use different directories by writing something like:
#nawk|gawk path1/redpuma.nawk path2/redpuma_ami_list.html>
#path3/redpuma
#START CONSTANT DECLARATIONS
# The next three variables must be changed if any station other than
# TUL is used
STA = "TUL"
geographic_latitude_STA = 35.910556 #use +,- not N,S
geographic_longitude_STA = -95.7925 #use +,- not E,W
colatSTA=deg2rad(90-g2c(geographic_latitude_STA))
if(geographic_longitude_STA >= 0){elonSTA = deg2rad(geographic_longitude_STA)}
if(geographic_longitude_STA < 0){elonSTA = deg2rad(360 + geographic_longitude_STA)}
# colat is colatitude, elon is east longitude
# function deg2rad converts degrees to radians
# function g2c converts geographic latitude to geocentric latitude
#P=0 phase is Pb depth 33km
#P=1 to P=13 phase is Pn depth 33km
#P=14 to P=98 phase is P depth 33km
#P=99 to P=113 phase is Pdiff depth 33km
#P=114 to P=180 phase is PKPdf depth 33km
P[0]=5.45 ;P[1]=17.69 ;P[2]=31.45 ;P[3]=45.20 ;P[4]=58.95
P[5]=72.69 ;P[6]=86.43 ;P[7]=100.16 ;P[8]=113.88 ;P[9]=127.60
P[10]=141.30 ;P[11]=154.99 ;P[12]=168.67 ;P[13]=182.33 ;P[14]=203.47
P[15]=214.59 ;P[16]=225.69 ;P[17]=236.77 ;P[18]=247.81 ;P[19]=258.80
P[20]=269.73 ;P[21]=280.56 ;P[22]=291.29 ;P[23]=301.90 ;P[24]=311.59
P[25]=320.70 ;P[26]=329.77 ;P[27]=338.79 ;P[28]=347.75 ;P[29]=356.65
P[30]=365.51 ;P[31]=374.33 ;P[32]=383.11 ;P[33]=391.85 ;P[34]=400.54
P[35]=409.17 ;P[36]=417.75 ;P[37]=426.27 ;P[38]=434.73 ;P[39]=443.12
P[40]=451.45 ;P[41]=459.71 ;P[42]=467.90 ;P[43]=476.03 ;P[44]=484.08
P[45]=492.07 ;P[46]=499.98 ;P[47]=507.83 ;P[48]=515.60 ;P[49]=523.30
P[50]=530.93 ;P[51]=538.49 ;P[52]=545.97 ;P[53]=553.38 ;P[54]=560.72
P[55]=567.99 ;P[56]=575.19 ;P[57]=582.31 ;P[58]=589.36 ;P[59]=596.34
P[60]=603.24 ;P[61]=610.08 ;P[62]=616.84 ;P[63]=623.52 ;P[64]=630.14
P[65]=636.68 ;P[66]=643.15 ;P[67]=649.55 ;P[68]=655.87 ;P[69]=662.12
P[70]=668.30 ;P[71]=674.41 ;P[72]=680.44 ;P[73]=686.40 ;P[74]=692.28
P[75]=698.09 ;P[76]=703.83 ;P[77]=709.49 ;P[78]=715.08 ;P[79]=720.59
P[80]=726.02 ;P[81]=731.38 ;P[82]=736.67 ;P[83]=741.87 ;P[84]=747.00
P[85]=752.04 ;P[86]=757.01 ;P[87]=761.91 ;P[88]=766.69 ;P[89]=771.40
P[90]=776.08 ;P[91]=780.73 ;P[92]=785.35 ;P[93]=789.95 ;P[94]=794.54
P[95]=799.10 ;P[96]=803.63 ;P[97]=808.13 ;P[98]=812.60 ;P[99]=817.04
P[100]=821.48 ;P[101]=825.92 ;P[102]=830.36 ;P[103]=834.80 ;P[104]=839.24
P[105]=843.68 ;P[106]=848.11 ;P[107]=852.55 ;P[108]=856.99 ;P[109]=861.43
P[110]=865.87 ;P[111]=870.31 ;P[112]=874.75 ;P[113]=879.19 ;P[114]=1115.42
P[115]=1117.34 ;P[116]=1119.25 ;P[117]=1121.16 ;P[118]=1123.08 ;P[119]=1124.99
P[120]=1126.90 ;P[121]=1128.82 ;P[122]=1130.72 ;P[123]=1132.63 ;P[124]=1134.54
P[125]=1136.44 ;P[126]=1138.34 ;P[127]=1140.24 ;P[128]=1142.13 ;P[129]=1144.03
P[130]=1145.91 ;P[131]=1147.79 ;P[132]=1149.67 ;P[133]=1151.54 ;P[134]=1153.40
P[135]=1155.25 ;P[136]=1157.10 ;P[137]=1158.94 ;P[138]=1160.76 ;P[139]=1162.58
P[140]=1164.38 ;P[141]=1166.17 ;P[142]=1167.94 ;P[143]=1169.70 ;P[144]=1171.43
P[145]=1173.15 ;P[146]=1174.84 ;P[147]=1176.50 ;P[148]=1178.14 ;P[149]=1179.76
P[150]=1181.34 ;P[151]=1182.88 ;P[152]=1184.40 ;P[153]=1185.87 ;P[154]=1187.31
P[155]=1188.70 ;P[156]=1190.05 ;P[157]=1191.36 ;P[158]=1192.62 ;P[159]=1193.83
P[160]=1195.00 ;P[161]=1196.11 ;P[162]=1197.17 ;P[163]=1198.18 ;P[164]=1199.13
P[165]=1200.04 ;P[166]=1200.88 ;P[167]=1201.67 ;P[168]=1202.41 ;P[169]=1203.08
P[170]=1203.70 ;P[171]=1204.26 ;P[172]=1204.77 ;P[173]=1205.21 ;P[174]=1205.60
P[175]=1205.92 ;P[176]=1206.19 ;P[177]=1206.40 ;P[178]=1206.55 ;P[179]=1206.64
P[180]=1206.67
#END CONSTANT DECLARATIONS
}
#START FUNCTIONS
#functions delta and azimuth are based on: Bullen, K. E., 1965. An
#introduction to the theory of seismology, Cambridge University Press, Third
#Edition, pages 154-155.
function deg2rad(degrees) {
return (degrees/180)*3.141592635
}
function rad2deg(radians) {
return (radians/3.141592635)*180
}
function arccos(x) {
return atan2(sqrt(1-(x*x)),x)
}
function g2c(geoglatdeg, geoclatrad)
{
geoglatrad=(geoglatdeg/180)*3.141592635
geoclatrad= atan2(.99327*sin(geoglatrad),cos(geoglatrad))
return (geoclatrad/3.141592635)*180
}
function delta(colatsta,elonsta,colatepi,elonepi, Asta,Bsta,Csta,Aepi,Bepi,Cepi) {
Asta=sin(colatsta)*cos(elonsta)
Bsta=sin(colatsta)*sin(elonsta)
Csta=cos(colatsta)
Aepi=sin(colatepi)*cos(elonepi)
Bepi=sin(colatepi)*sin(elonepi)
Cepi=cos(colatepi)
return arccos((Asta*Aepi)+(Bsta*Bepi)+(Csta*Cepi))
}
function azimuth(colatsta,elonsta,colatepi,elonepi,drad, Aepi,Bepi,Cepi,Dsta,Esta,Gsta,Hsta,Ksta,COSZ,leftside1,leftside2,SINZ) {
Aepi=sin(colatepi)*cos(elonepi)
Bepi=sin(colatepi)*sin(elonepi)
Cepi=cos(colatepi)
Dsta=sin(elonsta)
Esta=-cos(elonsta)
Gsta=cos(colatsta)*cos(elonsta)
Hsta=cos(colatsta)*sin(elonsta)
Ksta=-sin(colatsta)
leftside1=(((Aepi-Gsta)^2)+((Bepi-Hsta)^2)+((Cepi-Ksta)^2))
COSZ=(((leftside1)-2)/(2*sin(drad)))
leftside2=(((Aepi-Dsta)^2)+((Bepi-Esta)^2)+((Cepi)^2))
SINZ=(((leftside2)-2)/(2*sin(drad)))
return (atan2(SINZ,COSZ))
}
function decsec(HH,MM,SS) {
return ((3600*HH)+(60*MM)+SS)
}
# function dectimetosextime changes decimal time (in this program the
# number of seconds after 000000) to sexagesimal time: hhmmss
function dectime2sextime(decimalsecond, hour,minute,second) {
hour=int(decimalsecond/3600)
decimalsecond=decimalsecond - (3600*hour)
minute=int(decimalsecond/60)
second=decimalsecond - (60*minute)
return((10000*hour)+(100*minute)+second)
}
function PPKPTT(distdeg, ip,fp) {
ip=int(distdeg)
fp=distdeg-ip
if (ip==113) {return 879.19+fp*4.44} #Pdiff for ip=113
else {return P[ip]+(fp*(P[ip+1]-P[ip]))}
}
function LRTT(distancedeg) {
return ((distancedeg*111.199)/4.0)
}
#END FUNCTIONS
/^[0-3][0-9](Jan|Feb|Mar|Apr|May|Jun|Jul|Aug|Sep|Oct|Nov|Dec)/ {
print $0
hh=substr($2,1,2)
mm=substr($2,4,2)
ss=substr($2,7,2)
if(match($0,/[0-9]+\.[0-9]+[NS]/)){
latnumber=substr($0,RSTART, RLENGTH-1)
latsign=substr($0,RSTART+RLENGTH-1,1)
lat=latnumber
if(latsign=="S") {lat=-lat}
lat=g2c(lat)}
#print latnumber,latsign,lat
if(match($0,/[0-9]+\.[0-9]+[EW]/)){
lonnumber=substr($0,RSTART, RLENGTH-1)
lonsign=substr($0,RSTART+RLENGTH-1,1)
lon=lonnumber
if(lonsign=="W") {lon=-lon}}
#print lonnumber,lonsign,lon
colat=90-lat #colatitude 0deg at North Pole, 180deg at South Pole
elon=lon
if(lon<0) {elon=360+lon}
DELRAD=delta(colatSTA,elonSTA,deg2rad(colat),deg2rad(elon))
DELDEG=rad2deg(DELRAD)
ZDEG=rad2deg(azimuth(colatSTA,elonSTA,deg2rad(colat),deg2rad(elon),DELRAD))
if (ZDEG < 0) {ZDEG=ZDEG+360}
lrarrival=dectime2sextime(decsec(hh,mm,ss)+LRTT(DELDEG))
parrival =dectime2sextime(decsec(hh,mm,ss)+PPKPTT(DELDEG))
if (DELDEG < 14){pname="Pn"}
else if (DELDEG >= 14 && DELDEG < 99){pname="P"}
else if (DELDEG >= 99 && DELDEG < 114){pname="PDIFF"}
else if (DELDEG >= 114){pname="PKP"}
printf("D%7.3f Z%05.1f %s %s%06.0f LR%06.0f\n\n",DELDEG,ZDEG,STA,pname,parrival,lrarrival)
}