#! /bin/csh

#  〈作者の久保篤規さんによる解説〉
#
#  xyzの位置データを入力として、任意の方位の断面図を描く。主に用いるコマンドは
#  projectである。
#
#  断面の指定方法は３通りがあるが（-C , -A; -C -E; -C -T),ここでは、直交する
#  二つの断面をとることを考えて投影の基準位置と投影面の方位を指定する方法を用い
#  る(-Clon/lat -Aazimuth)。
#
#  -Lで投影面に沿った選択範囲を指定、-Wで投影面に直交方向の選択範囲を指定する。
#
#  -Fxyzpqで指定した形式の投影データファイルが作られる。ここでxyz はもとのデータ
#  p, q がそれぞれ投影断面に平行、及び直交方向の座標値である。選択範囲内のデータ
#  だけが処理される。
#
#  -Qはp,qの値をkmで返すために指定する。
#
#  p,z平面をpsxyで描けば断面図の完成。
#
#  例では簡単に選択領域を示すために水平面図に選択されたデータを上書きしてみまし
#  た。


# plot cross section for arbitrary azimuth
# input file: kanto.pde written by USGSscreen file format
# written by Astsuki Kubo
# prepare xyz data
awk '{ print $7,$6,$8 }' kanto.pde >! xyz.dat
#
# draw coast line
pscoast -K -X6. -Y20 -P -JM10 -R135/145/35/40 -Dh -B2/2::WESN -W1 -G220 >! hypo_section.ps
# plot  all epicenters
psxy -K -O -JM -R -Sc0.1 -G0 -B2/1::WeSn xyz.dat >> hypo_section.ps
# Set parameters for two perpendicular cross sections
set lat1=37.0
set lon1=141.5
set AZ=105
set AZ2=15
#
# projection for first cross section
project -Q -C$lon1/$lat1 -A$AZ xyz.dat -Fxyzpqrs -L-600/600 -W-50/50  >! xyzpqrs
awk '{ print $1,$2 }' xyzpqrs >! plot2.xy
#
# plot epicenters  selected for first cross section (overwrite)
psxy -K -O -JM -R -Sc0.1 -G200/0/0 plot2.xy >> hypo_section.ps
#
# projection for second cross section
project -Q -C$lon1/$lat1 -A$AZ2 xyz.dat -Fxyzpqrs -L-450/450 -W-50/50  >! xyzpqrs.2
awk '{ print $1,$2 }' xyzpqrs.2 >! plot3.xy
# plot epicenters selected for second cross section (overwrite)
psxy -K -O -JM -R -Sc0.1 -G0/0/200 plot3.xy >> hypo_section.ps
#
# plot center location of coordinate
echo $lon1 $lat1 >! tmp
psxy -K -O -Sx.5/.5 -JM -R -W10::0 tmp >> hypo_section.ps
# draw cross section 1
psbasemap -K -O -Y-8.5 -JX10/-6 -R-450/450/0/150 -B200/50 >> hypo_section.ps
awk '{ print $4,$3 }' xyzpqrs >! plot4.xy
psxy -K -O -JX -R -Sc0.1 -G200/0/0 plot4.xy >> hypo_section.ps
# draw cross section 2
psbasemap -K -O -Y-8.5 -JX10/-6 -R-450/450/0/150 -B200/50 >> hypo_section.ps
awk '{ print $4,$3 }' xyzpqrs.2 >! plot5.xy
psxy -O -JX -R -Sc0.1 -G0/0/200 plot5.xy >> hypo_section.ps