9月はじめの阿蘇火山地震観測実習で、本堂のトンネルに 広帯域地震計を設置しました。実習の後も、その地震計はずっと 地震動をとり続けています。データは、刻々と、 京都のＰＣに送られています。
9月29日15時19分32.0秒（世界標準時）（日本時間では30日00時19分32.0秒）に、 ニュージーランドのKermadec諸島で、マグニチュードMw7.0の地震が 発生しました。今年6月に起きた岩手・宮城内陸地震（Mw6.8）より、 少し大きな地震です。
皆さんが本堂に設置した広帯域地震計は、この地震からのＰ波をきれいに 記録しました。
今回の４回の演習では、そのＰ波の波形（"かたち"）をモデリングする（計算機上で 復元する）ことで、 地震の起こり方、断層の向き、断層のずれの継続時間、震源の深さ を推測します。
震源があるところでは、太平洋プレートが西に向かって オーストラリアプレートの下に沈み込もうとしています。 Kermadec地震は、太平洋プレートとオーストラリアプレートのプレート境界 の地震と思われます。
本堂の震央距離は、 79° です。 震央距離が 30-100°の範囲を、地震学では「遠地（teleseismic distance）」と 呼んでいます。 この範囲のＰ波の波形は、上部マントルや遷移層の複雑な構造の影響を受けにくく、 震源で起こったことを推測するときによく使用します。
地球の内部の速度構造を仮定すると、Ｐ波などの地震波が 伝播するにかかる時間を推測することができます。 "iasp91"という世界標準１次元速度モデルを用いると、 Ｐ波は本堂に719.6秒かかって到達し、 Ｓ波は本堂に1315.8秒かかって到達すると予想されます。 また、あわせて、 本堂に伝播するＰ波が、震源から鉛直下向きから上方に 23°の角度で射出したＰ波であることも推測されます。 （この角度を射出角take-off angleと呼んでいます）
地震の震源からみた本堂の方向（方位 azimuth）は、 -42°、本堂からみた地震の震源の方向（逆方位 back-azimuth）は、136°です。 いずれも、北をゼロとして時計回りを正に測られています。
★が地震、太い矢印が波の伝播方向、ihが射出角、φsが方位。
地動を見るためには、取得されたデータから、地震計の応答関数を 取り除く必要があります。 今回は、この作業を、SACというソフトウェアを使って行います。
阿蘇のデータがあるサーバにログインして、
作業を行ってください。
作業の詳細はココ
また、サーバへのログインやwin形式のデータを扱うコマンド等は、
第１回目の演習のメモにあります。
ココ
【演習1】断層運動から出るＰ波の振幅特性（Ｐ波放射特性）を計算する
断層の走向、傾斜角、すべり角（すべり方向）、Ｐ波の射出角と方位の５つの 変数を与えたときに、Ｐ波放射特性（Rp）を計算するプログラムをつくる。
この関係式はココ
プログラムの動作確認として、
(1) 横ずれ断層（走向90°、傾斜角90°、すべり角0°）で射出角が45°のとき、 地震波の伝播する方位を0から360°まで増やしていった場合の Ｐ波放射特性（Rp）をgnuplotで図にしてみる。
(2) 逆断層（走向0°、傾斜角45°、すべり角90°）で 地震波の伝播する方位が90°のとき、 射出角を0から180°まで増やしていった場合のＰ波放射特性（Rp）を gnuplotで図にしてみる。
Ｐ波が出来たら、ＳＶ波放射特性も同じように確認して。
【演習2】均質無限媒質中のＰ波（直達Ｐ波）の波形をつくる
この原理は上の放射特性のところに含まれている （ココ）。
ここでは、 畳み込み積分（convolution）を使いながら波形をつくってみる。 説明と手順は （ココ）。
【演習3】地表の効果として地表での反射波（pP）と変換波（sP）を加える
この説明はココ
【演習4】伝播の効果として非弾性減衰の効果を加える
この説明はココ
◇時間に余裕のある人（３回目の演習や４回目開始時にこの項目に達した人）
【演習3】の出来上がりの波形にフーリエ変換を行って、 周波数領域の中で、非弾性減衰の効果を掛け、 結果を逆フーリエ変換で時間領域に戻すように、 【演習3】で作ったプログラムを改良してください。 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 メインプログラムのファイル名 fcoorl.f
です。
逆フーリエ変換した後、実数部にだけデータがある（虚数部がゼロ）ためには、 周波数領域の複素データが、折り返し点をはさんで、 後半部と前半部で共役になっていることが必要です。 非弾性減衰を前半部の複素データに作用させ、 後半部をその共役な値で埋めることにより、これを実現します。
例えば、以下のような周波数領域で１６の複素データをもつ場合、 ９つ目のデータが折り返し点になります。 この折り返し点をはさんで前後のデータを共役にします。
上のFFTのプログラムでは、周波数領域でi番目の複素データ（前半部）は、
(i-1)/データ長
の周波数に対応します。
（ここで、データ長は、時間ステップ（△t）かける２のｎ乗です。)
◇時間に余裕のない人
非弾性減衰の応答関数として、デルタ関数に対する t*=1秒の場合にえられる時間関数（時間ステップ△t=0.2秒）を ココ におきました。このデータをプログラムで読み込み、 時間領域で畳み込み積分して、非弾性減衰の効果をいれてください。
(2) その結果をもとにしながら、断層の走向、傾斜角、すべり角、ずれの継続時間、 深さをどのように変えると、観測波形と計算波形が似るようになるかを調べる。
(3) もし同じ位置・深さで断層運動ではなく爆発が起こった場合、 どのような波形が予想されるかを調べる。
2009年1月22日（木）締め切り
Wordファイルかpdfファイルにまとめたものを久家宛にメールで送るか
（メールの件名（サブジェクト）を「ＤＣレポート」にすること）、
あるいは理2号館514号室前の封筒へ提出。
レポートには、名前を忘れずにいれること。