//SV波自上层介质左侧下行入射下层介质的Zopperitz求解 反射透射系数
//SV波入射时,存在临界角,当入射角大于临界角时,反射的P波不存在,代之以沿界面传播的非均匀波,此时出现全反射现象
//出现全反射时,波速满足的不等式存在多种情况,作为情形之一,仅讨论
//(1)不涉及全反射的sv波入射,类似于p波入射情况
//(2)vp1>vp2的情况
#include<stdio.h>
#include<math.h>
#include"矩阵求逆(子函数).C"
#include"矩阵相乘(子函数).C"
#define PI 3.1415926
//*在sv波入射时会出现全反射的现象,反射sv波的角度总是小于反射p波的角度,当反射p波角度大于PI/2时,反射p波沿x1方向传播,临界反射
void main()
{
FILE *fp1,*fp2;
int i,n[91];
float in,ipp1,ipp2,ipp3,x1,x2,y1,y2,pr,sr,pt,st,p,rsp,rss,tsp,tss,a1[91],a2[91],b1[91],b2[91],vp1,vp2,vs1,vs2,den1,den2,k;
printf("出现全反射时,波速满足的不等式存在多种情况,作为情形之一,仅讨论\n") ;
printf("(1)不涉及全反射的sv波入射,类似于p波入射情况(入射角度较小时)\n") ;
printf("(2)vp1>vp2>vs2>vs1的情况(入射角度大时存在全反射)\n\n ") ;
printf("请输出上层介质纵波速度(单位:m/s)\n");
scanf("%f", &vp1);
printf("请输出下层介质纵波速度(单位:m/s)\n");
scanf("%f", &vp2);
printf("请输出上层介质横波速度(单位:m/s)\n");
scanf("%f", &vs1);
printf("请输出下层介质横波速度(单位:m/s)\n");
scanf("%f", &vs2);
printf("请输出上层介质密度(单位:g/cm^3)\n");
scanf("%f", &den1);
printf("请输出下层介质密度(单位:g/cm^3)\n");
scanf("%f", &den2);
k=den2/den1;
static double ar[4][4],ai[4][4],br[4],bi[4],cr[4][1],ci[4][1];
fp1=fopen("SV波入射的反射透射角度正弦值.xls","w");
fp2=fopen("SV波入射的反射透射系数.xls","w");
//令反射或透射角度临界,求出入射SV波角度的最大值(临界值),方便后续讨论
ipp1=asin(vs1/vp1);//利用反射P波求出的入射角临界值 (单位PI)
ipp2=asin(vs1/vs2);//利用透射SV波求出的入射角临界值 (单位PI)
ipp3=asin(vs1/vp2);//利用透射P波求出的入射角临界值 (单位PI)
for(i=0;i<=90;i++)
{
n[i]=i;
in=i*PI/180;//入射角转化为单位PI格式
p=sin(in)/vs1;//求出该界面各个波所共同的射线参数P(斯奈尔定律)
x1=p*vs1; //求出反射SV波的正弦值
x2=p*vp1; //求出反射P波的正弦值
y1=p*vs2; //求出透射SV波的正弦值
y2=p*vp2; //求出透射P波的正弦值
a1[i]=x1;
a2[i]=x2;
b1[i]=y1;
b2[i]=y2;//遍历角度0-90°分别将各个波角度正弦值返回数组 (可能存在大于1的正弦值)
fprintf(fp1,"%d,%f,%f,%f,%f\n",n[i],a1[i],a2[i],b1[i],b2[i]);
//利用文件指针,分别输出整形格式的角度和浮点格式的反射SV波、反射P波、透射SV波、透射P波的正弦值到EXCEL文件
//讨论1:不存在全反射
//sinα、sinβ、sinα^、sinβ^ 值均小于1
if(in<ipp1)
{
sr=asin(x1);//求出反射SV波角度(单位PI)
pr=asin(x2);//求出反射P波角度
st=asin(y1);//求出透射SV波角度
pt=asin(y2);//求出透射P波角度
ar[0][0]=sin(pr);
ar[0][1]=cos(sr);
ar[0][2]=-sin(pt);
ar[0][3]=cos(st);//矩阵第一行
ar[1][0]=cos(pr);
ar[1][1]=-sin(sr);
ar[1][2]=cos(pt);
ar[1][3]=sin(st);//矩阵第二行
ar[2][0]=-vp1/vs1*cos(2*sr);
ar[2][1]=sin(2*sr);
ar[2][2]=k*vp2/vs1*cos(2*st);
ar[2][3]=k*vs2/vs1*sin(2*st);//矩阵第三行
ar[3][0]=-vs1/vp1*sin(2*pr);
ar[3][1]=-cos(2*sr);
ar[3][2]=-k*pow(vs2,2)/vs1/vp2*sin(2*pt);
ar[3][3]=k*vs2/vs1*cos(2*st);//矩阵第四行
br[0]=cos(sr);
br[1]=sin(sr);
br[2]=sin(2*sr);
br[3]=cos(2*sr);//右侧4*1矩阵
if(cinv(ar,ai,4)!=0)//调用子函数求系数矩阵的逆矩阵
{
tcmul(ar,ai,br,bi,4,4,1,cr,ci);//调用子函数复矩阵相乘,系数4*4矩阵乘右侧4*1矩阵后得到各个反射透射系数的复述形式并存在cr和ci行列式中
rsp=sqrt(cr[0][0]*cr[0][0]+ci[0][0]*ci[0][0]);
rss=sqrt(cr[1][0]*cr[1][0]+ci[1][0]*ci[1][0]);
tsp=sqrt(cr[2][0]*cr[2][0]+ci[2][0]*ci[2][0]);
tss=sqrt(cr[3][0]*cr[3][0]+ci[3][0]*ci[3][0]);//求sv波入射反射透射系数(复数)的大小
fprintf(fp2,"%d,%f,%f,%f,%f\n",n[i],rsp,rss,tsp,tss);
//利用文件指针,分别输出整形格式的角度和浮点格式的反射p波、反射sv波、透射p波、透射sv波的反射和透射系数
}
}
//讨论2:入射SV波角度大于反射P波求出的入射SV波的临界值,小于用透射P波求出的入射角临界值
//sinβ、sinα^、sinβ^ 值均小于1 ,sinα值均大于1
else if(in>=ipp1&&in<ipp3) //此时入射角度大于反射P波临界,但小于透射p波临界,反射p波角度为PI/2,沿x1轴正方向
{
sr=asin(x1);//求出反射SV波角度(单位PI)
//反射p波角度正弦值大于1 sinα=sin(PI/2-i*β) shβ^2=chβ^2-1=sqrt(cp1^2/
st=asin(y1);//求出透射SV波角度
pt=asin(y2);//求出透射P波角度
//反射p波正弦值大于1,对Zopperitz方程进行修改,其中p为射线参数
ar[0][0]=p*vp1, ai[0][0]=0;//根据snell定律,sin(in)/vs1=sin(反射纵波)/vp1,此处反射纵波正弦值大于1,得ar[0][0]
ar[0][1]=cos(sr), ai[0][1]=0;
ar[0][2]=-sin(pt), ai[0][2]=0;
ar[0][3]=cos(st), ai[0][3]=0;//矩阵第一行
ar[1][0]=0, ai[1][0]=-sqrt(p*p*vp1*vp1-1);//根据公式sin^2(α)+ cos^2(β)=1,结合ar[0][0] ,得ai[1][0]
ar[1][1]=-sin(sr), ai[1][1]=0;
ar[1][2]=cos(pt), ai[1][2]=0;
ar[1][3]=sin(st), ai[1][3]=0;//矩阵第二行
ar[2][0]=-cos(2*sr)*vp1/vs1, ai[2][0]=0;
ar[2][1]=sin(2*sr), ai[2][1]=0;
ar[2][2]=cos(2*st)*k*vp2/vs1, ai[2][2]=0;
ar[2][3]=sin(2*st)*k*vs2/vs1, ai[2][3]=0;//矩阵第三行
ar[3][0]=0, ai[3][0]=2*p*vp1*sqrt(p*p*vp1*vp1-1)*vs1/vp1;//第一行矩阵和第二行矩阵2*sin(α)*cos(α)
ar[3][1]=-cos(2*sr), ai[3][1]=0;
ar[3][2]=-sin(2*pt)*k*vs2*vs2/(vs1*vp2), ai[3][2]=0;
ar[3][3]=cos(2*st)*k*vs2/vs1, ai[3][3]=0;//矩阵第四行
br[0]=cos(sr), bi[0]=0;
br[1]=sin(sr), bi[1]=0;
br[2]=sin(2*sr), bi[2]=0;
br[3]=cos(2*sr), bi[3]=0;//右侧4*1矩阵
if(cinv(ar,ai,4)!=0)//调用子函数求系数矩阵的逆矩阵
{
tcmul(ar,ai,br,bi,4,4,1,cr,ci);//调用子函数复矩阵相乘,系数4*4矩阵乘右侧4*1矩阵后得到各个反射透射系数的复述形式并存在cr和ci行列式中
rsp=sqrt(cr[0][0]*cr[0][0]+ci[0][0]*ci[0][0]);
rss=sqrt(cr[1][0]*cr[1][0]+ci[1][0]*ci[1][0]);
tsp=sqrt(cr[2][0]*cr[2][0]+ci[2][0]*ci[2][0]);
tss=sqrt(cr[3][0]*cr[3][0]+ci[3][0]*ci[3][0]);//求sv波入射反射透射系数(复数)的大小
fprintf(fp2,"%d,%f,%f,%f,%f\n",n[i],rsp,rss,tsp,tss);
//利用文件指针,分别输出整形格式的角度和浮点格式的反射p波、反射sv波、透射p波、透射sv波的反射和透射系数
}
}
//讨论3:入射SV波角度大于透射P波求出的入射SV波的临界值,小于用透射SV波求出的入射角临界值
//sinβ、sinβ^ 值均小于1,sinα、sinα^值均大于1
else if(in>=ipp3&&in<ipp2)
{
sr=asin(x1);//求出反射SV波角度(单位PI)
st=asin(y1);//求出透射SV波角度(单位PI)
ar[0][0]=p*vp1, ai[0][0]=0;
ar[0][1]=cos(sr), ai[0][1]=0;
ar[0][2]=-p*vp2, ai[0][2]=0;
ar[0][3]=cos(st), ai[0][3]=0;//矩阵第一行
ar[1][0]=0, ai[1][0]=-sqrt(p*p*vp1*vp1-1);
ar[1][1]=-sin(sr), ai[1][1]=0;
ar[1][2]=0, ai[1][2]=-sqrt(p*p*vp2*vp2-1);
ar[1][3]=sin(st), ai[1][3]=0;//矩阵第二行
ar[2][0]=-cos(2*sr)*vp1/vs1, ai[2][0]=0;
ar[2][1]=sin(2*sr), ai[2][1]=0;
ar[2][2]=cos(2*st)*k*vp2/vs1, ai[2][2]=0;
ar[2][3]=sin(2*st)*k*vs2/vs1, ai[2][3]=0;//矩阵第三行
ar[3][0]=0, ai[3][0]=2*p*vp1*sqrt(p*p*vp1*vp1-1)*vs1/vp1;//第一行矩阵和第二行矩阵2*sin(α)*cos(α)
ar[3][1]=-cos(2*sr), ai[3][1]=0;
ar[3][2]=0, ai[3][2]=2*p*vp2*sqrt(p*p*vp2*vp2-1)*k*vs2*vs2/(vs1*vp2);//第一行矩阵和第二行矩阵2*sin(α)*cos(α)
ar[3][3]=cos(2*st)*k*vs2/vs1, ai[3][3]=0;//矩阵第四行
br[0]=cos(sr), bi[0]=0;
br[1]=sin(sr), bi[1]=0;
br[2]=sin(2*sr), bi[2]=0;
br[3]=cos(2*sr), bi[3
评论4