实验一 基本地震学理论（斯奈尔定律）

C-2. 编写子函数：对于固-固介质界面，已知P波和S波速度以及入射角，求出P波和S波反射及透射角。（注意：该子函数可以计算并列出入射波的临界角）

一、 实验原理

二、程序源代码

function [ a1,a2,b1,b2] = snel( i,a,vp1,vp2,vs1,vs2 )%本程序一些地方有中文，在某些版本的matlab中可能无法正确显示，但不影响使用。

i=input('请输入波形，p波请输入1，S波请输入2，SH波请输入3:');%........................%

a=input('请输入入射角 a:');

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

vp2=input('上层S波速度 vp2:');%如果作为子函数使用最好把这一块输入去掉，把输入调整到主函数中%

vs1=input('下层P波速度 vs1:');

vs2=input('下层S博速度 vs2:');%..........................................................%

if i==1 %P波

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; %S波

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 %SH波

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

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

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

end

end

三、 结果及分析

这个小问题并没有什么结果，特采用P波以30度角测试一下，选取模型，VP1=6.8km/s，VP2=8.0 km/s，VS1=3.9 km/s,VS2=4.8km/s。

30.000000 16.664362 36.031879 20.667316

p波临界角为：58.211669

sv波没有临界角