地学教材に「地震波経路のシミュレーションによる地球の内部構造の推定」を追加しました。

• English
• 日本語

地球内部を伝わる地震波には、シャドーゾーンと呼ばれる現象があります。震源から地球の中心核で103°から143°には内部から地震が伝わってこない範囲があることをいうのですが、この地震波の伝わるようすを、Pythonのプログラミングでシミュレーションしてみました。ページは→こちら

#数学関数とタートル（亀）グラフィックスのインポート
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()