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波没有临界角
文章评论