必要なファイル
channels.tbl
win.prm
win32tools(ディレクトリ)
win.prmの中身
win.prm
../trg /* default directory for data file */
channels.tbl /* channel table file */
../etc/zones.tbl /* file for grouping of stations */
../pick /* directory for pick files */
../bin/hypomh /* hypomh program */
../etc/struct.tbl /* velocity structure for hypomh */
../etc/map.japan.km /* map data file */
. /* output directory for cut-out wave data */
B4 /* data format (B4/B2/L4/L2) */
../etc/filter.prm /* filter setting file */
lp /* printer */
upper /* upper or lower hemisphere projection */
../etc/labels.tbl /* NOISE, FAR etc. */
../etc/final.man /* directory for hypocenter data */
100.0 /* printer's DPI for hardcopy */
/var/tmp /* temporary working directory */
cntをファイルをwinファイルに変換
$ w32tow1 2021100722410101VM.cnt 20211007.2241
$ w32tow1 2021100722420101VM.cnt 20211007.2242
$ w32tow1 2021100722430101VM.cnt 20211007.2243
$ w32tow1 2021100722440101VM.cnt 20211007.2244
$ w32tow1 2021100722450101VM.cnt 20211007.2245
結合
一分毎のwinファイルを五分間分のwinファイルにまとめる
$ cat 0211007.224? > 0211007.2241_5
winファイルをsacファイルに変換
$ win2sac_32 0211007.2241_5 7773 s > gomi.txt
※数字は観測点コード
これで3成分を作成
sacを起動
$ sac
> r N.KOTH.X.s
> rmean
> w KOTH_EW.s
これをやらないと後になぜかエラーが起きる
テキストファイルに変換（3成分）
$ gfortran -o sac2txt sac2txt.f90
$ ./sac2txt KOTH_EW.s > KOTH_EW.txt
これを3成分作成する
sac2txt.f90
program rensyumain
implicit none
integer i, a(40), maxx, maxy, maxz, p, maxnum, npts
real b(70),c(2),d,wa,maxtime,dt,ap,as,maxv
real, allocatable :: x(:), y(:), z(:), h(:), o(:),zz(:),xx(:),yy(:)
character :: file1*80, file2*80,file3*80, arg*10, khd(48)*4
integer :: j
real :: hfl(20), gl, eps
integer :: mfl, nnfc,number
real :: fp, fs, fl, fh
if(iargc()>0) then
call getarg (1, file1) !コマンドラインで引数!
! call getarg (2, file2) !コマンドラインで引数
! call getarg (3, file3) !コマンドラインで引数
! call getarg (2, arg) !コマンドラインで引数所得
! read(arg,*) number
endif
maxx=10000000
maxy=10000000
maxz=10000000
allocate (x(maxx),y(maxy),z(maxz),xx(maxx),yy(maxy),zz(maxz))
call arsac(file1,b,a,khd,x,maxx)
npts = a(10)
dt = b(1)
! call arsac(file2,b,a,khd,y,maxy)
! call arsac(file3,b,a,khd,z,maxz)
do i = 1, npts
xx(i) = abs(x(i))
! yy(i) = abs(y(i))
! zz(i) = abs(z(i))
enddo
maxv = maxval(xx)
do i=1,npts
x(i) = x(i)/maxv
enddo
!maxv = maxval(yy)
!do i=1,npts
! y(i) = y(i)/maxv
!enddo
!maxv = maxval(zz)
!do i=1,npts
! z(i) = z(i)/maxv - 2
!enddo
do i=1,npts
write(*,*)i*dt, x(i)
enddo
end program rensyumain
subroutine arsac(fn,fhd,nhd,khd,y,maxy)
real fhd(70),y(maxy)
integer nhd(40),i,i1,i2,i3,npts
character khd(48)*4,fn*80
open(1,file=fn,access='direct',recl=632,status='old')
read(1,rec=1)(fhd(i1),i1=1,70),(nhd(i2),i2=1,40),(khd(i3),i3=1,48)
npts=nhd(10)
close(1)
open(1,file=fn,access='direct',recl=4,status='old')
do i=1,npts
! read(1,rec=i+158,end=100)y(i) Mar.1 H.Aoyama
read(1,rec=i+158)y(i)
enddo
! 100 continue
close(1)
return
end subroutine arsac
波形を作成
$ ./wave.sh
wave.sh
# !/bin/bash
gmt set MAP_TICK_LENGTH_PRIMARY -0.50c
gmt set MAP_ANNOT_OFFSET_PRIMARY 0.6c
# 切り抜く範囲
range=0/100
# 拡大範囲
en_range=30/50
# 3つの観測点
st1="CBAH"
st2="KOTH"
st3="SHMH"
gmt psxy ${st1}_UD.txt -R${range}/-1/1 -Bf1a5:"Time [s]":/f1:"U-D Velocity":WeSn -W0.5 -JX25/5 -K > ${st1}.eps
gmt psxy ${st1}_NS.txt -R${range}/-1/1 -Bf1/f1:"N-S Velocity":WeSn -W0.5 -JX25/5 -O -K -Y5 >> ${st1}.eps
gmt psxy ${st1}_EW.txt -R${range}/-1/1 -Bf1/f1:"E-W Velosity":WeSn:.${st1}: -W0.5 -JX25/5 -O -Y5 >> ${st1}.eps
gmt psxy ${st2}_UD.txt -R${range}/-1/1 -Bf1a5:"Time [s]":/f1:"U-D Velocity":WeSn -W0.5 -JX25/5 -K > ${st2}.eps
gmt psxy ${st2}_NS.txt -R${range}/-1/1 -Bf1/f1:"N-S Velocity":WeSn -W0.5 -JX25/5 -O -K -Y5 >> ${st2}.eps
gmt psxy ${st2}_EW.txt -R${range}/-1/1 -Bf1/f1:"E-W Velosity":WeSn:.${st2}: -W0.5 -JX25/5 -O -Y5 >> ${st2}.eps
gmt psxy ${st3}_UD.txt -R${range}/-1/1 -Bf1a5:"Time [s]":/f1:"U-D Velocity":WeSn -W0.5 -JX25/5 -K > ${st3}.eps
gmt psxy ${st3}_NS.txt -R${range}/-1/1 -Bf1/f1:"N-S Velocity":WeSn -W0.5 -JX25/5 -O -K -Y5 >> ${st3}.eps
gmt psxy ${st3}_EW.txt -R${range}/-1/1 -Bf1/f1:"E-W Velosity":WeSn:.${st3}: -W0.5 -JX25/5 -O -Y5 >> ${st3}.eps
# 拡大波形
gmt psxy ${st1}_UD.txt -R${en_range}/-1/1 -Bf1a5:"Time [s]":/f1:"U-D Velocity":WeSn -W0.5 -JX25/5 -K > ${st1}_en.eps
gmt psxy ${st1}_NS.txt -R${en_range}/-1/1 -Bf1/f1:"N-S Velocity":WeSn -W0.5 -JX25/5 -O -K -Y5 >> ${st1}_en.eps
gmt psxy ${st1}_EW.txt -R${en_range}/-1/1 -Bf1/f1:"E-W Velosity":WeSn:.${st1}_Enlarged: -W0.5 -JX25/5 -O -Y5 >> ${st1}_en.eps
gmt psxy ${st2}_UD.txt -R${en_range}/-1/1 -Bf1a5:"Time [s]":/f1:"U-D Velocity":WeSn -W0.5 -JX25/5 -K > ${st2}_en.eps
gmt psxy ${st2}_NS.txt -R${en_range}/-1/1 -Bf1/f1:"N-S Velocity":WeSn -W0.5 -JX25/5 -O -K -Y5 >> ${st2}_en.eps
gmt psxy ${st2}_EW.txt -R${en_range}/-1/1 -Bf1/f1:"E-W Velosity":WeSn:.${st2}_Enlarged: -W0.5 -JX25/5 -O -Y5 >> ${st2}_en.eps
gmt psxy ${st3}_UD.txt -R${en_range}/-1/1 -Bf1a5:"Time [s]":/f1:"U-D Velocity":WeSn -W0.5 -JX25/5 -K > ${st3}_en.eps
gmt psxy ${st3}_NS.txt -R${en_range}/-1/1 -Bf1/f1:"N-S Velocity":WeSn -W0.5 -JX25/5 -O -K -Y5 >> ${st3}_en.eps
gmt psxy ${st3}_EW.txt -R${en_range}/-1/1 -Bf1/f1:"E-W Velosity":WeSn:.${st3}_Enlarged: -W0.5 -JX25/5 -O -Y5 >> ${st3}_en.eps
完成図
観測点地図
$ ./map.sh
map.sh
# !/bin/bash
gmt gmtset FONT_ANNOT_PRIMARY 16p
# 範囲
REGION='-R139.5/140.5/35/36'
# 投影
PROJECTION='-JM6.5i -P'
# 軸
AXIS='-Ba1f0.5'
# 海の色
COLOR_SEA='-S51/204/255'
# 枠を表示
gmt psbasemap -K ${AXIS} ${REGION} ${PROJECTION} > map.eps
gmt pscoast -K -O -W1 -Dh $COLOR_SEA ${MAP_SCALE} ${REGION} ${PROJECTION} -Bg0.5 >> map.eps #grd
# 観測点
gmt psxy station.txt -Gred -W0.5 -K -O ${REGION} ${PROJECTION} -St0.5 >> map.eps
gmt psxy station.txt -Gred -W0.5 -K -O ${REGION} ${PROJECTION} -Sx0.5 >> map.eps
# 震源
# gmt psxy hypocenter.txt -G255/165/0 -W0.5 -K -O ${REGION} ${PROJECTION} -Sa0.5 >> map_nagano.eps
# 県境界を表示
gmt psxy GMT_JP_kenkai.dat -K -O -W1 ${REGION} ${PROJECTION} >>map.eps
# 最後のおまじない
gmt psbasemap -O ${AXIS} ${REGION} ${PROJECTION} >> map.eps
GMT_JP_kenkai.dat 県境界用
station.txt 観測点経度緯度(10進数)
を同じディレクトリに入れておく
