地球内部を伝わる地震波には、シャドーゾーンと呼ばれる現象があります。震源から地球の中心からの角度で103°から143°には内部から地震が伝わってこない範囲があることをいうのですが、その理由は地球の中心にある外核が鉄が融けた液体からなり、そこで地震波が大きく屈折して伝わるためと考えられています。地学の教科書には、ほぼ上のような図とともに解説がありますが、そもそも地震波がなぜこのように弧を描くように伝わるのか、という理解も難しい面があります。波の進み方かたの光の屈折と同じなのですが、その根本まで説明して納得してもらうのは、大変です。そこで、というわけでもないのですが、簡単なプログラミングでシミュレーションすることができるのを示すのがよいのではないかと考えました。最近注目されているプログラミング言語のPythonをかじって、参考になる文献を見つけ 昔のNECパソコンのN88Basicで書かれたコードを移植してみました。 https://www.toray-sf.or.jp/awards/education/pdf/s61_13.pdf
あまり上手なプログラミングではありませんが、参考まで。アニメーションとして授業で見せることができます。
#数学関数とタートル（亀）グラフィックスのインポート
import math
import turtle
t = turtle.Turtle() # t;タートルオブジェクト（変数）
wn = turtle.Screen()
wn.mode("standard")
wn.bgcolor("#97e8ff")
wn.title("地震波線のシミュレーション")
#深さと地震波速度の関係式のパラメータ
dt=5.0 ;T=0.0
L=6400.0 ; L2=3500 ; L3=1300 #dtが計算単位時間、Tは走時、Lは球殻の半径
vs1=7.8 ; vs2=5.0 ; vs3=9.8 ; K1=0.0022 ; K2=0.001 ;K3=0.0003 #初速度と係数
#リストの作成
r=[0.0]*1000;a=[0.0]*1000;v=[0.0]*1000;q=[0.0]*1000;x=[0.0]*1000;y=[0.0]*1000
#地球を描く
t.speed(0)
t.pu() ;t.setpos(0,L*0.05);t.pd()
t.begin_fill();t.color("#000000","#f8f6cc");t.circle(-L*0.05)
t.end_fill(); t.pu()
t.setpos(0,L2*0.05);t.pd()
t.begin_fill();t.color("#000000","#d2d9f1")
t.circle(-L2*0.05) ;t.end_fill(); t.pu()
t.setpos(0,L3*0.05);t.pd()
t.begin_fill();t.color("#000000","#d7c8d4")
t.circle(-L3*0.05) ;t.end_fill();t.pu()
t.setpos(-L*0.05,0) ; t.pd()
t.hideturtle()
#深さR(半径）と中心角Aの計算
def calcRA():
r[i] = math.sqrt((r[i-1]**2)+((v[i-1]*dt)**2)-((-1)**flag)*2*r[i-1]*v[i-1]*dt*math.cos(q[i-1]))
c1=v[i-1]*dt*math.sin(q[i-1])/r[i]
a[i]=math.asin(c1)
return
#速さVによる屈折角Qの計算
def calcVQ():
global flag
c2=v[i]*math.sin(q[i-1]+((-1)**flag)*a[i])/v[i-1]
if c2 >= .999: #全反射するかの判定
q[i]=q[i-1]+a[i]
flag=1
return
else:
q[i]=math.asin(c2)
return
#描画（座標の計算）
def display():
global A
A += a[i]
y[i]=r[i]*math.sin(A)*0.05 ; x[i]=-r[i]*math.cos(A)*0.05
t.goto(x[i],y[i])
return
try:
#地震波線のループ
Q0=8.0 # 初めの入射角（伏角）
for n in range(0,40):
Q0 +=2.0 #増やす角度
q[0] = math.radians(90.0-Q0)
t.setpos(-L*0.05,0) ; t.pd() #原点に戻す
flag = 0 ; A = 0.0 ; T = 00 ; r[0] = L ; a[0] = 0.0
x[0]=0.0 ; y[0]=0.0 ; c1=0.0 ; c2=0.0
i = 0
while r[i] <= r[0]:
i +=1
T += dt
calcRA()
if r[i]<L3:
v[i]=K3*(r[0]-r[i])+vs3 ;v[0]=10.9 #内核のとき
elif r[i]<L2:
v[i]=K2*(r[0]-r[i])+vs2 ;v[0]=7.8 #外核のとき
else:
v[i]=K1*(r[0]-r[i])+vs1 ;v[0]=vs1 #マントルの内部
if v[i] !=0.0: #走時の計算
T += dt/v[i]
calcVQ()
display()
print("入射角（伏角）＝",Q0,"走時＝",round(T,1)," 到達角距離＝",round(math.degrees(A),1),"°")
t.pu()
wn.exitonclick()
except:
exit()
行77～84のif文は、地震波速度を深さによる一次式（直線）で近似していることを示します。不明な点がありましたら、お問い合わせください。