斯奈尔定律和Zoeppritz方程

斯奈尔定律和Zoeppritz方程姓名：学号：专业：地球物理勘察技术2012级一、实验目的1.利用密度、上下界面的纵横波速度通过斯奈尔定律求出该界面的0°~90°入射角的反射角度和透射角度。PVVVVVspspp222211111sinsinsinsinsin2.利用Zoeppritz方程绘制反射系数和透射系数曲线。二、实验步骤1、模型：Vp1=3300m/sVs1=1585m/sρ1=2.4g/cm3Vp1=3100m/sVs1=1989m/sρ1=2.24g/cm3计算从0～90度的反射透射系数曲线如图1-1-4三、实验结果P波入射的反射透射角度：图1P波入射的反射透射系数曲线：图2SV波入射的反射透射角度正弦值：图3SV波入射的反射透射系数曲线：图4SH波入射的反射透射角度正弦值：图5SH波入射的反射透射系数曲线：图6四、实验分析因为各波反射和透射波振幅系数与其能量成正比，由此可以看出其能量的变化以及入射波能量的分配。当上层介质为密介质，下层介质为疏介质时：由图2知，P波从上层密介质入射到界面时：随着入射角增大，P波反射系数、P波透射系数在减小，即随着入射角的增大，反射P波、透射P波的能量在减小，而S波反射系数、S波透射系数在增大，即随着入射角的增大，反射S波、透射S波的能量在增大。说明当入射角发生变化时，入射波的能量分配在改变。由图4可知，SV波从上层密介质入射到界面时：入射角小于29度时，随着入射角增大，SV波反射系数在减小，P波反射系数、SV波透射系数、P波透射系数在增大，SV波透射系数最大，即反射SV波能量在减小，反射P波、透射SV波、透射P波能量在增大，透射SV波的能量最大；入射角大于等于29度且小于31度时，SV波反射系数、SV波透射系数减小，P波反射系数、P波透射系数在增大，SV波透射系数最大，即随着入射角的增大，反射SV波、透射SV波能量在减小，反射P波、透射P波能量在增大，透射SV波能量最大；入射角大于等于31度且小于53度时，SV波透射系数增大，SV波反射系数、P波反射系数、P波透射系数先减小后增大，SV波透射系数仍最大，即随着入射角的增大，透射SV波能量增大，反射SV波、反射P波、透射P波能量先减小后增大，透射SV波能量最大；入射角大于53度时，SV波反射系数为1，SV波透射系数、P波反射系数、P波透射系数减小，即随着入射角的增大，反射SV波能量不变，透射SV波、反射P波、透射P波能量减小。由图6可知，SH波从上层密介质入射到界面时：不产生转换波，入射角小于临界角时，SH波的反射系数、透射系数均随着入射角的增大而增大，即反射、透射SH波能量增大；入射角大于等于临界角时，随着入射角的增大，SH波的反射系数几乎不变，透射系数减小，即反射SH波的能量减小，透射SH波的能量几乎不变。五、附：源程序代码P波入射时，程序：#includestdio.h#includemath.h#include6GAUS.C#definePI3.1415926voidmain(){FILE*fp1,*fp2;inti,n[91];doubleipp,x1,x2,y1,y2,pr,sr,pt,st,a1[95],a2[95],b1[95],b2[95],vp1=3300,vp2=3100,vs1=1585,vs2=1989,den1=2.4,den2=2.24,k=den2/den1;staticdoublea[4][4]={0.0},b[4]={0.0};fp1=fopen(snell.csv,w);fp2=fopen(P波入射反射透射系数.csv,w);for(i=0;i=90;i++){n[i]=i;ipp=i*PI/180;x1=ipp;a1[i]=x1*180/PI;x2=sin(ipp)*vs1/vp1;a2[i]=asin(x2)*180/PI;y1=sin(ipp)*vp2/vp1;b1[i]=asin(y1)*180/PI;y2=sin(ipp)*vs2/vp1;b2[i]=asin(y2)*180/PI;fprintf(fp1,%d,%f,%f,%f,%f\n,n[i],a1[i],a2[i],b1[i],b2[i]);//*输出P波入射反射投射角度*//pr=a1[i]*PI/180;sr=a2[i]*PI/180;pt=b1[i]*PI/180;st=b2[i]*PI/180;a[0][0]=sin(pr);a[0][1]=cos(sr);a[0][2]=-sin(pt);a[0][3]=-cos(st);a[1][0]=cos(pr);a[1][1]=-sin(sr);a[1][2]=cos(pt);a[1][3]=-sin(st);a[2][0]=sin(2*pr);a[2][1]=cos(2*sr)*vp1/vs1;a[2][2]=sin(2*pt)*k*vp1*pow(vs2/vs1,2)/vp2;a[2][3]=cos(2*st)*k*vp1*vs2/pow(vs1,2);a[3][0]=cos(2*sr);a[3][1]=-sin(2*pr)*vs1/vp1;a[3][2]=-cos(2*st)*k*vp2/vp1;a[3][3]=sin(2*pt)*k*vs2/vp1;//*输入系数矩阵*//b[0]=-sin(pr);b[1]=cos(pr);b[2]=sin(2*pr);b[3]=-cos(2*sr);//*输入常数矩阵*//if(gaus(a,b,4)!=0)fprintf(fp2,%d,%f,%f,%f,%f\n,n[i],b[0],b[1],b[2],b[3]);}fclose(fp1);fclose(fp2);}SV波入射时，程序：#includestdio.h#includemath.h#include4CINV.C#include4TCMUL.C#definePI3.1415926voidmain(){FILE*fp1,*fp2;inti,n[91];doublein,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=3300,vp2=3100,vs1=1585,vs2=1989,den1=2.4,den2=2.24,k=den2/den1;staticdoublear[4][4],ai[4][4],br[4],bi[4],cr[4][1],ci[4][1];fp1=fopen(snell.csv,w);fp2=fopen(SV波入射反射透射系数.csv,w);ipp1=asin(vs1/vp1);ipp2=asin(vs1/vs2);ipp3=asin(vs1/vp2);//*临界角*//for(i=0;i=90;i++){n[i]=i;in=i*PI/180;p=sin(in)/vs1;x1=p*vs1;x2=p*vp1;y1=p*vs2;y2=p*vp2;//*snell定律a1[i]=x1;a2[i]=x2;b1[i]=y1;b2[i]=y2;fprintf(fp1,%d,%f,%f,%f,%f\n,n[i],a1[i],a2[i],b1[i],b2[i]);//*输出S波入射反射透射角度正弦值*//if(inipp1){sr=asin(x1);pr=asin(x2);st=asin(y1);pt=asin(y2);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]=2*cos(pr)*sin(pr)*vs1/vp1;ar[2][1]=-(1-2*sin(sr)*sin(sr));ar[2][2]=2*cos(pt)*sin(pt)*k*vs2*vs2/(vs1*vp1);ar[2][3]=(1-2*sin(st)*sin(st))*k*vs2/vs1;ar[3][0]=(1-2*sin(sr)*sin(sr))*vp1/vs1;ar[3][1]=2*cos(sr)*sin(sr);ar[3][2]=-(1-2*sin(st)*sin(st))*k*vp2/vp1;ar[3][3]=2*cos(st)*sin(st)*k*vs2/vp1;//*输入系数矩阵br[0]=cos(sr);br[1]=sin(sr);br[2]=1-2*sin(sr)*sin(sr);br[3]=2*cos(sr)*sin(sr);if(cinv(ar,ai,4)!=0)//*求系数矩阵的逆矩阵*//{tcmul(ar,ai,br,bi,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]);//*求S波入射反射透射系数*//fprintf(fp2,%d,%f,%f,%f,%f\n,n[i],rsp,rss,tsp,tss);//*输出S波入射反射透射系数*//}}elseif(in=ipp1&&inipp3){sr=asin(x1);st=asin(y1);pt=asin(y2);ar[0][0]=-p*vp1,ai[0][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);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]=0,ai[2][0]=2*p*vp1*sqrt(p*p*vp1*vp1-1)*vs1/vp1;ar[2][1]=-(1-2*sin(sr)*sin(sr)),ai[2][1]=0;ar[2][2]=2*cos(pt)*sin(pt)*k*vs2*vs2/(vs1*vp1),ai[2][2]=0;ar[2][3]=(1-2*sin(st)*sin(st))*k*vs2/vs1,ai[2][3]=0;ar[3][0]=(1-2*sin(sr)*sin(sr))*vp1/vs1,ai[3][0]=0;ar[3][1]=2*cos(sr)*sin(sr),ai[3][1]=0;ar[3][2]=-(1-2*sin(st)*sin(st))*k*vp2/vp1,ai[3][2]=0;ar[3][3]=2*cos(st)*sin(st)*k*vs2/vp1,ai[3][3]=0;//*输入系数矩阵br[0]=cos(sr),bi[0]=0;br[1]=sin(sr),bi[1]=0;br[2]=1-2*sin(sr)*sin(sr),bi[2]=0;br[3]=2*cos(sr)*sin(sr),bi[3]=0;//*输入常数矩阵if(cinv(ar,ai,4)!=0)//*求系数矩阵的逆矩阵*//{tcmul(ar,ai,br,bi,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]*c