必要なファイル

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進数)
を同じディレクトリに入れておく