#!/bin/sh

# 震源メカニズムを書くスクリプト．
# 一行めには，"#!"の後に起動するシェルを指定する．
# この場合はBourne Shellを起動している．

rm -f .gmtcommand .gmtdefaults

# .gmtcommand .gmtdefaults のファイルがあると，
# GMTの動きを邪魔されることがあるので，
# 消しておく方が無難かな？


OPTION="-Jm40.0 -R130.15/130.55/31.80/32.13 -B10mNseW"

# -J，-R，-Bのような共通のオプションを，
# シェル変数"OPTION"に代入している．
# シェル変数に代入する文字列に空白が含まれる場合には，
# クォーテーションでくくる必要がある．
#
# -Jm40.0     メルカトル図法で表示する．
#             "40.0"の値は，１度を40cmで書くように指定している．
#
# -R130.15/130.55/31.80/32.13
#             表示する範囲の指定．
#             -R/経度最小/経度最大/緯度最小/経度最大 という具合．
#             経度は130.15から130.55まで，緯度は31.80から32.13の意．
#
# -B10mNseW   枠に表示する数値や目盛の指定．
#             "10m"は，10minite（分）ごとに表示する意．
#             "NseW"は"東西南北"を表し，大文字になっている方角に
#             数値と目盛を表示する指定．
#             この場合は，N（上）W（左）に数値と目盛を表示するの意．


gmtset ANOT_FONT_SIZE 15 DEGREE_FORMAT 3

### gmtset: GMTのparameter値を変更 ###
#
# ANOT_FONT_SIZE 15  軸の数値を書くフォントの大きさを指定（ポイント）．
#
# DEGREE_FORMAT 3    緯度・経度（軸の数値）の書き方を，
#                    ０〜５の数値での指定．
#                    ３の場合には，経度が0-180までの値で，
#                    緯度が-90〜90までの値．
#
#
#

psbasemap ${OPTION} -X3 -Y5 -P -K >| map.ps

### psbasemap: 枠や目盛・数値・タイトルなどを表示する ###
#${OPTION}については，"set OPTION = ...."の所を参照．
#
# -X3 -Y5   図の原点の位置を，紙の左下からの距離で指定．
#           X方向に3cm，Y方向に5cmの位置を原点にしている．
#
# -P        紙の方向を縦置き(Portrait)に指定．
#           このオプションをつけない時には，横置き(Landscape)になる．
#
# -K        継続オプション．
#           更に続けて絵を書く場合に指定する．
#           このオプションをつけないと，次の命令で絵を重ね書き出来ない．
#
# >| map.ps 出力リダイレクトの強制上書きモード．
#           既にmap.psというファイルがある場合でも，
#           強制的に上書きして保存する． 


pscoast -Jm -R -B -Lf130.45/31.85/31.85/10 -S200 -Df -W1 -P -O -K >>map.ps

### pscoast: 海岸線などを表示 ###
# -Jm -R -B は，前の指定をそのまま使う．
#
# -Lf130.45/31.85/31.85/10
#          スケールの表示をさせる．
#          fは"fancyな"（白黒のきれいなやつ）スケールを表示の意．
#          スケールを表示する位置が(130.45deg,33deg)，
#          スケールの基準を緯度31.85degで指定している．（２つ目の31.85の意）
#          長さ10kmのスケールを表示するように指定．
#
# -S200    海(sea)の色を指定．
#          この場合，グレースケールで指定してあり，
#          0で真っ黒，255で真っ白の意なので，黒っぽいグレーな色．
#
# -Df      海岸線のデータの解像度を指定．
#          この場合は"f"(full)なので，最も高い解像度で海岸線を表示している．
#
# -W1      海岸線をフチドリさせる．線の太さは１．数値が大きいほど太い．
#
# -P       縦置き(Portrait)オプション．
#
# -O       上書き(over write)オプション．
#          他のGMTの命令の後に絵を上書きする場合に付ける必要あり．
#
# -K       継続オプション．
#
# >> map.ps 出力リダイレクトの追加書きモード．
#           map.psファイルの後に追加して保存させている．


# plot aftershocks
psxy -Jm -R -B -Sc0.1 -G180 -P -O -K <lon_lat.txt >>map.ps

### psxy: 余震の位置をプロット ###
# 位置データは，"lon_lat.txt"に格納されている．
#
# -Jm -R -B オプションは，前の指定を使う．
#
# -Sc0.1  位置をプロットするデータのシンボルを指定．
#         "c"は円(circle)でデータの位置をプロットするの意．
#         直径を0.1cmで表示させている．
#
# -G180   シンボルの中を塗りつぶす色の指定．
#         この場合はグレースケールで指定してある．
#         0で真っ黒，255で真っ白．
#
# -P      縦置き(Portrait)オプション．
#
# -O      上書き(over write)オプション
#
# -K      継続オプション．
#
# <lon_lat.txt  入力リダイレクト．
#              位置データを，ファイル"lon_lat.txt"から入力している．
#
# >> map.ps     出力リダイレクトの追加書きモード．


# plot star
psxy -Jm -R -B -Sa -G0 -P -O -K <<END >>map.ps
130.37997   31.96959   0.5
END

### psxy: 震源の位置をプロット ###
# -Sa     位置をプロットするデータのシンボルを指定． 
#         "a"は星型でデータの位置をプロットするの意．
#         シンボルの外接円の半径は，入力するデータの第３フィールドで指定．
#         この場合は0.5cmに指定．
#
# -G0     シンボルを塗りつぶす色の指定． 
#         0なので真っ黒．
#
# <<END   ヒアドキュメントの指定．
#         ENDという文字列だけのデータの行までを，
#         psxyコマンドに標準入力させる．
#         この場合，130.37997，31.96959，0.5というデータが入力される．
#         データの意味は，経度，緯度，シンボルの半径．



# plot diamond
psxy -Jm -R -B -Sd -G0 -P -O -K <<END >>map.ps
130.33528   31.93846   0.4
END

### psxy: 震源の位置をプロット ###
# -Sd     位置をプロットするシンボルの指定．
#         "d"は菱形(diamond)でデータの位置をプロットするの意．
#         シンボルの辺の長さは，入力データの第３フィールドで指定する．
#         この場合は0.4cmに指定．


# draw line 
psxy -Jm -R -B -W3 -M -P -O -K <<END >>map.ps
130.37033   32.02033
130.37997   31.96959
>
130.28167   31.92833
130.33528   31.93846
END

### psxy: データを線で結ぶ ###
# psxyでシンボルの指定(-S)をしない場合，データの位置を線で結ぶ．
#
# -W3   データとデータを結ぶ線の太さを指定．
#       この場合は太さ"3"に指定．数値が大きいほど太い．
#
# -M    複数のデータセットを結ぶ指定．(Multi data)
#       データの区切りは">"のみの行になる．
#       この場合は区切りが１つなので，データセットが２つ． 
#       よって２本の線を引いている．



# plot station
psxy -Jm -R -B -St -G0 -P -O -K <<END >>map.ps
130.3524    31.9677   0.4
END

### psxy: 観測点のプロット ###
# -St   プロットするデータのシンボルを指定．
#       "t"は三角形(triangle)でデータをプロットするの意．
#       シンボルの辺の長さはデータの第３フィールドで指定(0.4cm)


# plot station name
pstext -Jm -R -P -O <<END >>map.ps
130.3524    31.9627   13   90.0    1   10   Station
END

### pstext: 観測点名を書く ###
# 入力データは，
# X座標（経度），Y座標（緯度），フォントサイズ， 回転角（反時計周り），
# フォント番号，位置の基準，表示する文字列
# の７項目．よって，データは第７フィールドまで必要．
#
# この場合は，
# 経度130.3524，緯度31.9627の位置に，フォントサイズ13で，回転角は0.0，
# フォント番号１(Helvetica)，位置の基準は10
# 表示する文字列は"Station"
# となっている．


# draw mechanism
# For data format of GMT3.0
#function psvelomeca(){ psmeca -o ;}
#psvelomeca -Jm -R -Sa3.0 -G150 -L -P -O <<END>>map.ps
#130.36997   32.06059    11   88  -169   6.5  0  0  Mechanism1 
#130.25511   31.91137   192   87   156   6.2  0  0  Mechanism2
#END

### psvelomeca: 震源のメカニズムを表示 (GMT3.0-3.3)###
# GMTの古いバージョンでは，震源メカニズムの表示は"psvelomeca"を使っていた．
# 新しいバージョンのGMTになって，震源メカニズムを表示するコマンドが
# "psmeca"に変わり，入力するデータのフォーマットも変更された．
#
# psvelomecaのデータフォーマットを使う場合は，
# psmecaを"-o"オプションで実行する．
# シェル関数を使って，psvelomecaコマンドで"psmeca -o"コマンドを実行している．
#
# 入力データのフォーマットは以下のとおり．必要なデータは９項目．
# よって，データは第９フィールドまで必要．
# offset_xとoffset_yは，longitudeとlatitudeからの距離．
# psvelomecaでは，深さのデータが無いのに注意！
#
# longitude latitude strike dip rake magnitude offset_x offset_y comment
#
# -Sa3.0  データフォーマットの指定．
#         aオプションのとき，strike，dip，rake は Aki & Richard の定義に従う．
#         3.0の値は，Magnitudeが５の時の直径．
#         この場合，メカニズムの図の直径は Magnigude * 3 / 5 (cm) になる．
#         （マニュアルには半径とあるが，直径が正解のようだ．）
#
# -G150   押し引き図の塗りつぶしの色指定．
#         0から255までの数値を指定する．0が黒，255が白．
#         この場合は150なので，中間的な灰色になる．
#
# -L      メカニズムの図を縁取りする指定．

# For GMT3.3.1 User
psmeca -Jm -R -Sa3.0 -G150 -L -P -O <<END>>map.ps
130.36997   32.06059  8.0   11   88  -169   6.5  0  0  Mechanism1 
130.25511   31.91137  8.0  192   87   156   6.2  0  0  Mechanism2
END

### psmeca: 震源のメカニズムを表示 (GMT3.3.1-)###
#
# 新しいバージョンのGMTでは，震源のメカニズム表示に"psmeca"コマンドを使う．
# psvelomecaとの違いは入力データのフォーマットで，深さが必要．
# コマンドのオプションなどについては，psvelomecaと同じ．
#
# 以下に，データフォーマットを示す．必要なデータは１０項目．
# よって，データは第１０フィールドまで必要．
# 第３フィールドに「深さ」のデータがあることに注意．
# （ただし，「深さ」の値は作図には直接関係ないようです．）
#
# longitude latitude 'depth' strike dip rake magnitude offset_x offset_y comment