基本地震学理论（射线路径绘制）

C-4. 编写程序：给定介质两边速度，用C-2的结果生成如图2.6-11 所示的射线；用该程序计算核幔边界可能入射波的射线路径（问题C-3）

下地幔： α = 13.7 km/s, β = 7.2 km/s, ρ = 5.5 g/cm³

地核： α₂ = 8.0 km/s, β₂ = 0.0 km/s, ρ₂ = 9.9 g/cm³

一、 实验原理

运用斯奈尔定律求出特定地层模型下，P波或者S波不同入射角下的反射波和透射波的角度，根据波射线的角度绘制一个入射点所有波射线的图像。根据不同入射角划分x轴为多个点，分别绘制不同入射角的波产生的所有射线。

二、 程序源代码

1.主函数

i=input('请输入第几幅图 i');%图号

vp1=input('input vp1:');%上层p波速度

vp2=input('input vp2:');%下层p波速度

vs1=input('input vs1:');%上层s波速度

vs2=input('input vs2:');%下层s波速度

if i==1;%第一幅图p波从上方入射

x=[0 0];

y=[2 -1];

plot(x,y,'k')

hold on

for a=10:10:90

[ a1,a2,b1,b2 ]=snel( i,a,vp1,vp2,vs1,vs2 );

k=-tand(90-a);

k1=tand(90-a1);

k2=tand(90-a2);

k3=-tand(90-b1);

k4=-tand(90-b2);%计算斜率

x1=-2/k;%与x轴交点

x=[0 x1];

y1=k*x+2;

plot(x,y1,'k')%入射波

hold on

x=[x1 (x1+abs(1/k1))];

y2=k1*x-k1*x1;

plot(x,y2,'k')%反射波

hold on

x=[x1 (x1+abs(1/k2))];

y3=k2*x-k2*x1;

plot(x,y3,'--k')%反射波

x=[x1 (x1+abs(1/k3))];

y4=k3*x-k3*x1;

plot(x,y4,'k')%透射波

x=[x1 (x1+abs(1/k4))];

y5=k4*x-k4*x1;

plot(x,y5,'--k');%透射波

x=[0 18];y=[0 0];

plot(x,y,'k');

end

end

if i==2%第二幅图p波从下方入射

x=[0 0];

y=[2 -1];

plot(x,y,'--k')

hold on

for a=10:10:90

[ a1,a2,b1,b2 ]=snel( i,a,vp1,vp2,vs1,vs2 );

k=-tand(90-a);

k1=tand(90-a1);

k2=tand(90-a2);

k3=-tand(90-b1);

k4=-tand(90-b2);%计算斜率

x1=-2/k;

x=[0 x1];

y1=k*x+2;

plot(x,y1,'--k')

hold on

x=[x1 (x1+abs(1/k1))];

y2=k1*x-k1*x1;

plot(x,y2,'--k')

hold on

x=[x1 (x1+abs(1/k2))];

y3=k2*x-k2*x1;

plot(x,y3,'k')

x=[x1 (x1+abs(1/k3))];

y4=k3*x-k3*x1;

plot(x,y4,'--k')

x=[x1 (x1+abs(1/k4))];

y5=k4*x-k4*x1;

plot(x,y5,'k')

x=[0 18];y=[0 0];

plot(x,y,'k');

end

end

if i==3%第三幅图s波从上方入射

x=[0 0];

y=[-2 1];

plot(x,y,'k')

hold on

t=i-2;

for a=10:10:90

[ a1,a2,b1,b2 ]=snel( t,a,vp2,vp1,vs2,vs1 );

k=tand(90-a);

k1=-tand(90-a1);

k2=-tand(90-a2);

k3=tand(90-b1);

k4=tand(90-b2);%计算斜率

x1=2/k;

x=[0 x1];

y1=k*x-2;

plot(x,y1,'k')

hold on

x=[x1 (x1+abs(1/k1))];

y2=k1*x-k1*x1;

plot(x,y2,'k')

hold on

x=[x1 (x1+abs(1/k2))];

y3=k2*x-k2*x1;

plot(x,y3,'--k')

x=[x1 (x1+abs(1/k3))];

y4=k3*x-k3*x1;

plot(x,y4,'k')

x=[x1 (x1+abs(1/k4))];

y5=k4*x-k4*x1;

plot(x,y5,'--k')

x=[0 18];y=[0 0];

plot(x,y,'k');

end

end

if i==4%第四幅图s波从下方入射

x=[0 0];

y=[-2 1];

plot(x,y,'--k')

hold on

t=i-2;

for a=10:10:90

[ a1,a2,b1,b2 ]=snel( t,a,vp2,vp1,vs2,vs1 );

k=tand(90-a);

k1=-tand(90-a1);

k2=-tand(90-a2);

k3=tand(90-b1);

k4=tand(90-b2);%计算斜率

x1=2/k;

x=[0 x1];

y1=k*x-2;

plot(x,y1,'--k')

hold on

x=[x1 (x1+abs(1/k1))];

y2=k1*x-k1*x1;

plot(x,y2,'--k')

hold on

x=[x1 (x1+abs(1/k2))];

y3=k2*x-k2*x1;

plot(x,y3,'k')

x=[x1 (x1+abs(1/k3))];

y4=k3*x-k3*x1;

plot(x,y4,'--k')

x=[x1 (x1+abs(1/k4))];

y5=k4*x-k4*x1;

plot(x,y5,'k')

x=[0 18];y=[0 0];

plot(x,y,'k');

end

end

hold off

2.子函数

function [ a1,a2,b1,b2] = snel( i,a,vp1,vp2,vs1,vs2 )

%i=input('input i:');

%a=input('input a:');

%vp1=input('input vp1:');

%vp2=input('input vp2:');

%vs1=input('input vs1:');

%vs2=input('input vs2:');

if i==1

a1=a;

a2=asind(sind(a)*vs1/vp1);

b1=asind(sind(a)*vp2/vp1);

b2=asind(sind(a)*vs2/vp1);

fprintf('%f\t%f\t%f\t%f\n',a1,a2,b1,b2);

if vp2>vp1

r1=asind(vp1/vp2);

fprintf('p波临界角为：%f\n',r1);

else fprintf('p波没有临界角\n');

end

if vs2>vp1

r2=asind(vp1/vs2);

fprintf('sv波临界角为：%f\n',r2);

else fprintf('sv波没有临界角\n');

end

end

if i==2;

a2=asind(sind(a)*vp1/vs1);

a1=asind(sind(a)*vs1/vs1);

b2=asind(sind(a)*vp2/vs1);

b1=asind(sind(a)*vs2/vs1);

fprintf('%f\t%f\t%f\t%f\n',a1,a2,b1,b2);

if vp2>vs1

r2=asind(vs1/vp2);

fprintf('p波临界角为：%f\n',r2);

else fprintf('p波没有临界角\n');

end

if vs2>vs1

r1=asind(vs1/vs2);

fprintf('sv波临界角为：%f\n',r1);

else fprintf('sv波没有临界角\n');

end

end

if i==3

a2=asind(sind(a)*vs1/vs1);

b2=asind(sind(a)*vs2/vs1);

fprintf('%f\n%f\n',a2,b2);

end

end

1. 结果及分析

备注：所有图片中实线为P波，虚线为S波

1.P波从下地幔入射到外核

2.S波从下地幔入射到外核

3.P波从外核入射到下地慢

结果分析：因为外核为液体，所以外核没有S波。故没有从外核入射到下地慢的S波。因为上层速度大于下层速度，所以当P波从下层入射时，当角度大于临界角的时候，就不会再有透射P波，因为下地慢S波速度小于外核，故透射S波不会消失。因为上地幔P波速度大于下地慢P波速度，故无论P波还是S波从下地慢入射到地核，都不会产生折射波。外核P波总是存在。