Hi-netの地震波形の記録から震源決定の教材

#必要なファイル
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
```
完成図
![スクリーンショット 2021-10-25 22.57.04.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1101426/150fc48f-b472-2ed1-e542-25e2cd70926f.png)

![スクリーンショット 2021-10-25 22.56.53.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1101426/3f930381-c571-56b2-0cef-e11532422dcb.png)

#観測点地図
`$ ./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進数)
を同じディレクトリに入れておく

![スクリーンショット 2021-10-25 23.01.14.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1101426/f591694e-09fa-c500-f300-6c12ad3c083c.png)