#!/bin/csh

＜作者の関口さん＠京大防災研による解説＞
作業は、地震（余震等）の震源を任意の鉛直断面とそれに共役な面とに投
影し、投影した鉛直断面図を描くというものです。

入力ファイルは、震源の経度、緯度、深さ（ｋｍ）、マグニチュードの４
列からなるデータです。
そのほか与えるべき情報は、当該領域での緯度・経度１度に対応する長さ
（ｋｍ），投影面上の基準点座標、投影面の伸びる方向（角度）等です。

このスクリプトのキーポイントは、緯度経度座標からメートル法座標への
変換、及び、任意の鉛直断面への投影に必要なデータの２次元面内での回
転を、awk を使って、全部シェルスクリプト内でやってしまう点です。
awk の文章中にシェル中の変数を直接入力することができないので、echo 
を介して awk のコマンド１行からなるサブシェルスクリプトを作り、これ
を動かす、という方法を取りました。サイン、コサインの計算をシェルで
行う方法が無いと思われるので、このためにごくごく短いフォートランの
プログラム（sine.fとcosine.f）を使います。

#^^^^^^ plotting is done with measure unit = cm
gmtset MEASURE_UNIT	cm
gmtset DOTS_PR_INCH	300

set post = aftershock.section.ps

#^^^^^^ input information
# input hypocenter data file 
set fin = hypo.in
# length for 1deg in latitude and in longitude
set ulat = 110.922
set ulon = 91.7216
# latitude and logintude of origin 
set olat = 34.598
set olon = 135.044
# direction of plane extension to project data points
set theta = 50
set theta2 = `echo "$theta + 90" |bc`
# amplification of magnitude for plotting circles with varying size due to magnitude
set magamp = 0.0125

\rm $post out.tmp1 out.tmp2

#^^^^^^ transform lon/lat cordinate to metric cordinate
\rm bat1
echo "awk '{print ("'$1'"-$olon)*$ulon,("'$2'"-$olat)*$ulat,"'$3'","'$4'"*$magamp}' $fin > out.tmp1" >bat1
chmod +x bat1
bat1

#^^^^^^ project data points onto the plane
\rm sine.in
echo $theta > sine.in
set sn = `sine < sine.in`
set cs = `cosine < sine.in`
\rm bat2
echo "awk '{print ("'$1'"*$sn+"'$2'"*$cs),("'$1'"*$cs-"'$2'"*$sn),"'$3'","'$4'"}' out.tmp1 > out.tmp2" > bat2
chmod +x bat2
bat2

#^^^^^^ plotting section view along "theta" degree direction
awk '{print $1,$3,$4}' out.tmp2 | psxy -R-30/50/-30/0 -Jx0.2 -Sc -K -Ba5f10 -X6. -Y-1. > $post
pstext -R -Jx -N -K -O <<END>> $post
-30  8.  20 0 7 1 Section View of Aftershock Hypocenters
37  3.3  15 0 7 1 N${theta}E
10  2.5  10 0 0 2 km
END
psxy -R-30/50/-30/30 -Jx -W6 -N -K -O <<END>> $post
35 3
50 3
47.5  5
END

#^^^^^^ plotting section view along "theta2" degree direction

awk '{print $2,$3,$4}' out.tmp2 | psxy -R-10/10/-30/0 -Jx0.2 -Sc -K -O -B -X9. -Y-9. >> $post
pstext -R -Jx -N -K -O <<END>> $post
-3  3.3  15 0 7 1  N${theta2}E
END
psxy -R-10/70/-30/30 -Jx -W6 -N -O <<END>> $post
-5 3
10 3
7.5  5
END

gmtset MEASURE_UNIT	inch
gmtset DOTS_PR_INCH	300