Estimate of the internal structure of the earth by the simulation of the earthquake wave (seismic ray)

• English
• 日本語

An earthquake wave to be transmitted through the inside of the earth includes the phenomenon called the shadow zone. This is that from 103 to 143 degrees at an angle of from the centrosphere dose not come the seismic wave. The reason consists of the liquids which the outer core in the centrosphere was able to dissolve iron in, and it is thought that therefore an earthquake wave is greatly refracted and reaches it.

I rode a programming cord of Python representing the shadow zone of the earthquake wave. I can show it by a class as animation。

I made an execute file. You download it from this, and please try it.　→　From this

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