王彥東地震波作業(yè)

1、10 平面波在介質(zhì)分界面上的反射與透射 學(xué) 院：地球物理與空間信息學(xué)院 學(xué) 號(hào)：20121002799 班級(jí)序號(hào)：060121-17 姓 名：王彥東 指導(dǎo)教師：張學(xué)強(qiáng) 目錄一、 方向假設(shè)：3二、理論基礎(chǔ)3三、 圖像繪制及解釋4四、能量守恒6五、 總結(jié)6六、源代碼6平面波在介質(zhì)分界面上的反射與透射1、 方向假設(shè)：假設(shè)所有縱波及橫波的傳播方向及質(zhì)點(diǎn)振動(dòng)方向如圖二、理論基礎(chǔ) 佐普里茲（Zoppritz）方程 反射系數(shù)與透射系數(shù)定義 R 稱為反射系數(shù)，T 稱為透射系數(shù)；其下標(biāo)表示波的類型。反射或透射波與入射波同屬一個(gè)類型，稱為同類型波；否則，稱為轉(zhuǎn)換波。3、 圖像繪制及解釋利用書(shū)上數(shù)據(jù)：系列一對(duì)應(yīng)AB

2、橫波反射系數(shù)系列二對(duì)應(yīng)DE橫波透射系數(shù)系列三對(duì)應(yīng)GH縱波反射系數(shù)系列四對(duì)應(yīng)JK縱波透射系數(shù) 如圖所見(jiàn)，當(dāng)P波入射角不大時(shí)id<20-30。反射系數(shù)和透射系數(shù)變化不大。當(dāng)入射角較大時(shí)，轉(zhuǎn)換波反射系數(shù)和透射系數(shù)為極大值、其振幅最大。當(dāng)入射波垂直或平行入射到分界面時(shí)，無(wú)波的轉(zhuǎn)換現(xiàn)象發(fā)生。 還可以看到，縱波反射系數(shù)以及橫波透射系數(shù)為負(fù)，說(shuō)明其質(zhì)點(diǎn)振動(dòng)方向與假設(shè)方向相反，當(dāng)縱波在介質(zhì)表面發(fā)生發(fā)射時(shí)質(zhì)點(diǎn)振動(dòng)方向與波的傳播方向相反，此時(shí)為拉伸狀態(tài)，過(guò)后變?yōu)閴嚎s狀態(tài)并且兩種狀態(tài)不斷交替重復(fù)。更改數(shù)據(jù)后：系列一對(duì)應(yīng)AB橫波反射系數(shù)系列二對(duì)應(yīng)DE橫波透射系數(shù)系列三對(duì)應(yīng)GH縱波反射系數(shù)系列四對(duì)應(yīng)JK縱波透射系

3、數(shù) 將上下層參數(shù)對(duì)調(diào)后求出的圖像如圖所示，從所得的數(shù)據(jù)發(fā)現(xiàn)id>53之后所有的系數(shù)均為常值，此時(shí)透射角大于90度，產(chǎn)生了不均勻平面波，此時(shí)的反射系數(shù)與透射系數(shù)為復(fù)數(shù)，程序并未對(duì)這樣的系數(shù)進(jìn)行求解，所以顯示為常數(shù)。從圖中可以發(fā)現(xiàn)，當(dāng)id=48度左右時(shí)橫波透射系數(shù)為零，且曲線左右橫波透射系數(shù)的值發(fā)生異號(hào)，這說(shuō)明質(zhì)點(diǎn)的振動(dòng)方向發(fā)生了突變。當(dāng)id<30度時(shí)，透射系數(shù)變化不大，當(dāng)入射角增加時(shí)（小于53度），可以見(jiàn)到縱波反射系數(shù)的極小值以及橫波反射系數(shù)的極大值（絕對(duì)值），此時(shí)其振幅最大，反射能量最大。與第一種情況相比，此時(shí)的縱波反射系數(shù)為正，說(shuō)明其質(zhì)點(diǎn)振動(dòng)方向與假設(shè)方向相同，當(dāng)縱波在介質(zhì)表面發(fā)

4、生發(fā)射時(shí)質(zhì)點(diǎn)振動(dòng)方向與波的傳播方向相同，此時(shí)為壓縮狀態(tài)，過(guò)后變?yōu)槔鞝顟B(tài)并且兩種狀態(tài)不斷交替重復(fù)。四、能量守恒能量守恒公式： 其中，縱波與橫波速度，為介質(zhì)密度，id、is、td、ts分別問(wèn)縱波入射（反射）角度，橫波反射角度，縱波透射角度，橫波透射角度。 通過(guò)程序很容易可以驗(yàn)證上述兩種情況的能量均守恒。5、 總結(jié) 在老師的指點(diǎn)和同學(xué)的幫助下簡(jiǎn)單的畫(huà)出了反射與透射系數(shù)的圖像并對(duì)圖像進(jìn)行了一定程度的分析，但是由于代碼的局限性只能對(duì)有限的情況進(jìn)行圖像的繪制與分析，但對(duì)于反射與透射系數(shù)為復(fù)數(shù)情況還不能順利的完成。在這次實(shí)習(xí)中認(rèn)識(shí)到了自己對(duì)于c語(yǔ)言的不足，有的地方還不能搞懂，只是懂了一些皮毛，還有對(duì)于su

5、rfer的使用也只是皮毛。只有好好鍛煉過(guò)才能完美的應(yīng)用。六、源代碼#include "stdio.h"#include "stdlib.h"#include "math.h"short AGAUS(double *A,double *B,short N,double *X,short *JS);void main(void)FILE *fp;double dsta,vp1,vs1,vp2,vs2,pi,sid,sis,cis,cid,ctd,sts,std,cts,vsis,vsid,vsts,vstd,p1,p2,z1,z2,w1,

6、w2;double *rpp,*rps,*tpp,*tps,*angl,*sneng,rvps1,rvps2;double *a,*b,*x,s;short i,m,n,*js,L; char cmrk3="TX"p1=2700.0;p2=900.0;vp1=4500.0; vs1=2200.0; vp2=3600.0;vs2=1700.0;rvps1=vs1/vp1;rvps2=vs2/vp2;z1=p1*vp1;z2=p2*vp2;w1=p1*vs1;w2=p2*vs2;m=181; n=4;dsta=90.0/(m-1);pi=4.0*atan(1.0);/申請(qǐng)內(nèi)存空

7、間 rpp=new doublem; rps=new doublem; tpp=new doublem; tps=new doublem; angl=new doublem; sneng=new doublem; js=new shortn; x =new doublen; b=new doublen; a=new double*n; for(i=0;i<n;i+) ai=new doublen;for(i=0;i<m;i+)angli=i*dsta;/縱波入射角度vsid=angli/180.0*pi; /縱波入射角及反射角弧度vsis=asin(vs1*sin(vsid)/vp

8、1);/橫波反射角弧度vsts=asin(vs2*sin(vsid)/vp1);/橫波折射角弧度vstd=asin(vp2*sin(vsid)/vp1);/縱波折射角弧度sid=sin(2.0*vsid); sis=sin(2.0*vsis); cis=cos(2.0*vsis);cid=cos(2.0*vsid);sts=sin(2.0*vsts);ctd=cos(2.0*vstd);std=sin(2.0*vstd);cts=cos(2.0*vsts);/系數(shù)設(shè)置a00=sin(vsid); a01=cos(vsis); a02=-sin(vstd); a03=cos(vsts); b0=

9、-sin(vsid); a10=cos(vsid); a11=-sin(vsis); a12=cos(vstd); a13=sin(vsts); b1=cos(vsid);a20=z1*cis; a21=-w1*sis; a22=-z2*cts; a23=-w2*sts; b2=-z1*cis;a30=rvps1*w1*sid; a31=w1*cis; a32=rvps2*w2*std; a33=-w2*cts; b3=rvps1*w1*sid;L=AGAUS(a,b,n,x,js);printf("%d , L=%dn",i,L);rppi=x0;rpsi=x1;tppi

10、=x2;tpsi=x3;s=cos(vsid);if(s<1.0e-30)s=1.0e-30;elses=cos(vsid); /防止分母為0snengi=1.0-(rppi*rppi+rpsi*rpsi*vs1*cos(vsis)/(s*vp1)+p2*vp2*tppi*tppi*cos(vstd)/(s*p1*vp1)+p2*vs2*tpsi*tpsi*cos(vsts)/(s*p1*vp1);/能量守恒fp=fopen("縱波反射系數(shù).bln","wt");fprintf(fp,"%d,1n",m);for(i=0;i&l

11、t;m;i+) fprintf(fp,"%f %fn",angli,rppi);fclose(fp);fp=fopen("橫波反射系數(shù).bln","wt");fprintf(fp,"%d,1n",m);for(i=0;i<m;i+) fprintf(fp,"%f %fn",angli,rpsi);fclose(fp);fp=fopen("縱波透射系數(shù).bln","wt");fprintf(fp,"%d,1n",m);for(i=0

12、;i<m;i+) fprintf(fp,"%f %fn",angli,tppi);fclose(fp);fp=fopen("橫波透射系數(shù).bln","wt");fprintf(fp,"%d,1n",m);for(i=0;i<m;i+) fprintf(fp,"%f %fn",angli,tpsi);fclose(fp);fp=fopen("能量守恒.bln","wt");fprintf(fp,"%d,1n",m);for(i

13、=0;i<m;i+) fprintf(fp,"%f %fn",angli,snengi*100.0);fclose(fp);short AGAUS(double *A,double *B,short N,double *X,short *JS)short I,J,K,L,IS; double D,T;L=1;for (K=0; K<N-1; K+) D=0.0; for (I=K;I<N;I+) for (J=K;J<N;J+) if (fabs(AIJ)>D) D=fabs(AIJ); JSK=J; IS=I; if (D<1.0e-30) L=0; else if (JSK!=K) for (I=0;I<N;I+) T=AIK; AIK=AIJSK; AIJSK=T; if (IS!=K) for(J=K;J<N;J+) T=AKJ; AKJ=AISJ; AISJ=T; T=BK; BK=BIS; BIS=T; if (L=0) return L; for(J=K+1;J<N;J+) AKJ=AKJ/AKK; BK=BK/AKK; for(I=K+1;I<N;I+) for(J=K+1;J<N;J+) AIJ=AIJ-AIK*AKJ; BI=BI-