地震計でとられたデータで使えるものは、 地震波が到達した時刻だけではありません。 地震波の波形（”かたち”）にはたくさんの情報が含まれています。 理論をもとに、観測された地震波の波形を計算機でモデリングする（復元する） ことによって、 地震波が伝わってきた地下構造や、地震波を放射した地震そのものの断層運動など の特性を知ることができます。
今回の４回の演習では、実際に、簡単な理論から Ｐ波の波形（"かたち"）をモデリングする（計算機でつくる）ことで、 地震の断層・ずれの向き、断層のずれの継続時間、震源の深さ などの一部を推測してみましょう。
9月の阿蘇火山地震観測実習で、本堂のトンネルに 広帯域地震計を1台設置しました。実習の後も、その地震計はずっと 地震動をとり続けています。今回の演習では、このデータを用いて、 10月16日にパプアニューギニアの海域で起こった、 マグニチュードMw6.8の地震を調べます。
調べる地震は、 10月16日10時30分58秒（世界標準時）（日本時間では16日19時30分58秒）に 発生しました。 地震の規模（マグニチュード）は、 1995年兵庫県南部地震（Mw6.9）より少し小さいです。
パプアニューギニアのテクトニクス
（USGS NEIC
http://comcat.cr.usgs.gov/earthquakes/eventpage/usb000kemb#summary
から引用）。
本堂の震央距離は、 45.2° です。 距離は、地球中心からの角度で表しています。 震央距離が 30-100°の範囲を、地震学では「遠地（teleseismic distance）」と 呼んでいます。 この範囲のＰ波の波形は、 地球内部を伝播するときに、地震波速度が急変する 上部マントルやマントル遷移層の影響を受けにくいために、 地震の震源域で起こったことを推測するときによく使用します。
地球の内部の速度構造を仮定すると、Ｐ波などの地震波が 伝播するにかかる時間を推測することができます。 "ak135"と呼ばれる世界標準１次元地震波速度モデルを用いると、 Ｐ波は本堂に493.6秒かかって到達し、 Ｓ波は本堂に891.1秒かかって到達すると予想されます。 また、あわせて、 本堂に伝播するＰ波が、震源から鉛直下向きから上方に 35°の角度で射出したＰ波であることも推測されます。 （この角度を射出角take-off angleと呼んでいます）
なお、地震の震央からみた本堂の方向（方位 azimuth）は、 -29°、本堂からみた地震の震央の方向（逆方位 back-azimuth）は、146°です。 いずれも、北をゼロとして時計回りを正に測られています。
★が地震、太い矢印が波の伝播方向、ihが射出角、φsが方位。
地動を見るためには、取得されたデータから、地震計の応答関数を 取り除く必要があります。 今回は、この作業を、SACというソフトウェアを使って行います。
阿蘇のデータがあるサーバにログインして、
作業を行ってください。
作業の詳細はココ
サーバへのログインやwin形式のデータを扱うコマンド等は、
第１回目の演習のメモにあります。
ココ
遠地P波の変位波形（変位の時間関数）は、波線理論をもとに、
で、近似的に計算できます。
ここで、
S(t)は地震の断層運動によって決まる時間関数、
E(t)は地下構造によって生じる種々の波の到着を与える時間関数
（ここでは、直達P波とともに、震源そばの地表での
反射波や変換波を与える時間関数）、
P(t)は震源から観測点まで伝わる時の非弾性の効果を与える時間関数。
演習では、これらを計算するプログラムを作成します。
内容を説明するpdfノートはココ
【1-1】 断層運動から出るＰ波とＳ波の放射特性を計算する
断層の走向、傾斜角、すべり角（すべり方向）、Ｐ波の射出角と方位の５つの 変数を与えたときに、Ｐ波放射特性（Rp）を計算するプログラムをつくる。
プログラムの動作確認として、
(1) 横ずれ断層（走向90°、傾斜角90°、すべり角0°）で射出角が45°のとき、 地震波の伝播する方位を0から360°まで増やしていった場合の Ｐ波放射特性（Rp）をgnuplotで図にしてみる。
(2) 逆断層（走向0°、傾斜角45°、すべり角90°）で 地震波の伝播する方位が90°のとき、 射出角を0から180°まで増やしていった場合のＰ波放射特性（Rp）を gnuplotで図にしてみる。
Ｐ波が出来たら、ＳＶ波放射特性(Rsv)も作成して、同じように確認して。
【1-2】地表での反射波（pP）と変換波（sP）を考慮して、E(t)を作成する
E(t)をgnuplotで図にして、出てくる結果を確認。
【2-1】震源の時間関数S(t)を計算するプログラムをつくり、動作を確認
S(t)をgnuplotで図にして、出てくる結果を確認。
【2-2】出来たS(t)を計算するプログラムを使って、E(t)*S(t)を計算させる
E(t)*S(t)をgnuplotで図にして、出てくる結果を確認。
時間の余裕にあわせて、以下２つから、どちらかを選択して行う。
◇時間に余裕のある人（３回目の演習や４回目開始時にこの項目に達した人）
【ステップ2】の出来上がりの波形にフーリエ変換を行って、 周波数領域の中で、非弾性減衰の効果を掛け、 結果を逆フーリエ変換で時間領域に戻すように、 【ステップ2】で作ったプログラムを改良してください。 t*は１秒を使ってください。
フーリエ変換のプログラムはfcoolr.f
メインプログラムから
integer:: n
complex:: f(2048)
real:: rind
call fcoolr( n, f, rind )
のようにして呼び出すことができます。
ここで、配列fに含まれるデータは(2のn乗)個。 フーリエ変換ではrind=-1、逆フーリエ変換ではrind=1となります。 フーリエ変換のとき、変換された結果は、(2のn乗)倍になっていますので、 注意してください。
コンパイルの仕方は
gfortran メインプログラムのファイル名 fcoolr.f
です。
逆フーリエ変換した後、実数部にだけデータがある（虚数部がゼロ）ためには、 周波数領域の複素データが、折り返し点をはさんで、 後半部と前半部で共役になっていることが必要です。 非弾性減衰を前半部の複素データに作用させ、 後半部をその共役な値で埋めることにより、これを実現します。
例えば、以下のような周波数領域で１６の複素データをもつ場合、 ９つ目のデータが折り返し点になります。 この折り返し点をはさんで前後のデータを共役にします。
上のFFTのプログラムでは、周波数領域でi番目の複素データ（前半部）は、
(i-1)/データ長
の周波数に対応します。
（ここで、データ長は、時間ステップ（△t）かける２のｎ乗です。)
折り返し点が、ナイキスト周波数(振動数）になります。
◇時間に余裕のない人
非弾性減衰の応答関数として、デルタ関数に対する t*=1秒の場合にえられる時間関数（時間ステップ△t=0.2秒）を ココ におきました。このデータをプログラムで読み込み、 時間領域で畳み込み積分して、非弾性減衰の効果をいれてください。
(2) その結果をもとにしながら、パプアニューギニアの地震において、 断層の走向、傾斜角、すべり角、ずれの継続時間、 深さがどのぐらいの値であるかを、 計算波形と観測波形の比較から推測する。 どのパラメタが決まりやすいのか、決まりにくいのかも考える。
2014年1月17日（金）締め切り
Wordファイルかpdfファイルにまとめたものを久家宛にメールで送ってください。
メールの件名（サブジェクト）を「ＤＣレポート」にしてください。
レポート内には、名前を忘れずにいれてください。