地震計でとられたデータで使えるものは、 地震波が到達した時刻だけではありません。 地震波の波形（”かたち”）にはたくさんの情報が含まれています。 理論をもとに、観測された地震波の波形を計算機でモデリングする（復元する） ことによって、 地震波が伝わってきた地下構造や、地震波を放射した地震そのものの断層運動など の特性を知ることができます。
今回の４回の演習では、実際に、簡単な理論から Ｐ波の波形（"かたち"）をモデリングする（計算機でつくる）ことで、 地震の断層・ずれの向き、断層のずれの継続時間、震源の深さ などの一部を推測してみましょう。
9月の阿蘇火山地震観測実習で、本堂のトンネルに 広帯域地震計を1台設置しました。 今回の演習では、この地震計が取得したデータを用いて、 9月28日にインドネシアのスラウェシ島で起こった、 マグニチュードMw7.5の地震を調べます。 この地震では、津波で大きな被害がでました。 津波の被害のようすは、発生当初、テレビなどで何度も報道されていました。 今でも、メディア等のページや東北大学災害科学国際研究所のインドネシアパル地震の特別ページなどで、どのような津波の被害だったかを見ることができます。
アメリカ地質調査所（USGS）によって決められている地震の震央の位置は、
南緯 0.256° 東経 119.846°
です。
地震の発生時刻は、9月28日10時02分45.25秒（世界標準時）（日本時間では28日19時02分45.25秒）と推定されています。 地震の規模（モーメントマグニチュード）は Mw7.5と推定されています。
（左図）これまでの地震の震央分布（IRIS https://www.iris.edu/hq/retm/event/5647 から引用）
（右図）本震から10日以内に起こった地震の震央分布（IRIS https://www.iris.edu/hq/retm/event/5647 から引用）
（上図）インドネシアからパプアニューギニア周辺のテクトニクス
（USGS から引用）。
本堂までの震央距離は、 34.6° です。 距離は、地球中心を見込む角度で表しています。 震央距離が 30-100°の範囲を、地震学では「遠地（teleseismic distance）」と 呼んでいます。 この範囲のＰ波の波形は、 地球内部を伝播するときに、地震波速度が急変する 上部マントルやマントル遷移層の影響を受けにくいために、 地震の震源域で起こったことを推測するときによく使用します。
地球の内部の速度構造を仮定すると、Ｐ波などの地震波が 伝播するにかかる時間を推測することができます。 "ak135"と呼ばれる世界標準１次元地震波速度モデルを用いると、 Ｐ波は本堂に408秒かかって到達し、 Ｓ波は本堂に736秒かかって到達すると予想されます。 また、あわせて、 本堂に伝播するＰ波が、震源から鉛直下向きから上方に 30.6°の角度で射出したＰ波であることも推測されます。 （この角度を射出角take-off angleと呼んでいます）
なお、地震の震央からみた本堂の方向（方位 azimuth）は、 16.8°、本堂からみた地震の震央の方向（逆方位 back-azimuth）は、-160.0°です。 いずれも、北をゼロとして時計回りを正に測られています。
★が地震、太い矢印が波の伝播方向、ihが射出角、φsが方位。
地動を見るためには、取得されたデータから、地震計の応答関数を 取り除く必要があります。 今回は、この作業を、SACというソフトウェアを使って行います。
阿蘇のデータがあるサーバにログインして、
作業を行ってください。
作業の詳細はココ
※地震計の応答関数を与える変数の情報を修正しました（18/12/18）
遠地P波の変位波形（変位の時間関数）は、波線理論をもとに、
で、近似的に計算できます。
ここで、tは時間で、
S(t)は地震の断層運動によって決まる時間関数、
E(t)は地下構造によって生じる種々の波の到着を与える時間関数
（ここでは、直達P波とともに、震源そばの地表での
反射波や変換波を与える時間関数）、
P(t)は震源から観測点まで伝わる時の非弾性の効果を与える時間関数。
演習では、これらを計算するプログラムを作成します。
内容を説明するpdfノートはココ
※修正版に変更しました（18/12/12）
【1-1】 断層運動から出るＰ波とＳ波の放射特性を計算する
断層の走向、傾斜角、すべり角（すべり方向）、Ｐ波の射出角と方位の５つの 変数を与えたときに、Ｐ波放射特性（Rp）を計算するプログラムをつくります。
プログラムの動作確認として、
(1) 横ずれ断層（走向0°、傾斜角90°、すべり角0°）で射出角が30°のとき、 地震波の伝播する方位を0から360°まで増やしていった場合の Ｐ波放射特性（Rp）をgnuplotで図にしてみます。
(2) 逆断層（走向0°、傾斜角45°、すべり角90°）で 地震波の伝播する方位が90°のとき、 射出角を0から180°まで増やしていった場合のＰ波放射特性（Rp）を gnuplotで図にしてみます。
Ｐ波が出来たら、ＳＶ波放射特性(Rsv)も作成して、同じように確認してください。
【1-2】地表での反射波（pP）と変換波（sP）を考慮して、E(t)を作成する
[1-1]のRpとRsvを使って、 深さが与えられたときのE(t)を作成するプログラムを作ります。 できたE(t)をgnuplotで図にして、出てくる結果を確認してください。
【2-1】震源の時間関数S(t)を計算するプログラムをつくり、動作を確認
継続時間が与えられたときの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) (1)で選んだケースにおいて、（1）で仮定した値から、 断層の走向、傾斜角、すべり角、ずれの継続時間、深さを変えると、 本堂の理論変位波形がどのように変化するかを調べてください。
(3)（2）の結果と、本堂の観測変位波形との比較から、
断層の走向、傾斜角、すべり角、ずれの継続時間、深さがどのぐらいの値になるか、
推測してください。
どのパラメタが決まりやすいのか、決まりにくいのかも、考察してください。
2019年1月31日（木）締め切り
Wordファイルかpdfファイルにまとめたものを久家宛にメールで送ってください。
メールの件名（サブジェクト）を「ＤＣレポート」にしてください。
レポート内には、名前を忘れずにいれてください。