實驗四位場邊緣識別程序設(shè)計實驗匯總

1、實驗四位場邊沿鑒識程序設(shè)計實驗匯總實驗四位場邊沿鑒識程序設(shè)計實驗匯總20/20實驗四位場邊沿鑒識程序設(shè)計實驗匯總重磁資料辦理與解說實驗四位場邊沿鑒識程序設(shè)計實驗專業(yè)名稱：地球物理學(xué)學(xué)生姓名：學(xué)生學(xué)號：指導(dǎo)老師：王萬銀、紀(jì)新林、紀(jì)曉琳、邱世燦提交日期：2016-1-31、基本源理地質(zhì)目標(biāo)體邊沿時指斷裂結(jié)構(gòu)線、不一樣樣地質(zhì)體的界限限，其實是擁有必然密度或磁性差其余地質(zhì)體的界限限，在地質(zhì)體的邊沿周邊，重、磁異樣變化率較大，故全部的邊沿鑒識方法均利用這一特色進(jìn)行設(shè)計。此刻有重、磁位場邊沿鑒識方法分為數(shù)理統(tǒng)計、數(shù)值計算以及其余三大類。數(shù)值類邊沿鑒識方法均利用極大值地點或零值地點確立地質(zhì)體的邊沿地點，其依

2、照的理論基礎(chǔ)是二度體鉛垂臺階模型的重力異樣特色。在該模型邊沿處重力異樣總水平導(dǎo)數(shù)和分析信號振幅達(dá)到極大值、垂導(dǎo)游數(shù)達(dá)到零值。故能夠利用這些特色地點來確立二度體鉛垂臺階的邊沿地點，確立傾斜二度體、不規(guī)則二度體及三度體邊沿地點的理論均為二度體鉛垂臺階模型理論的推行，但確立的邊沿地點和真切的地點有必然的誤差。該誤差跟著地質(zhì)體界限形狀、埋深、水平尺寸及物性差別等的變化而變化。所以，邊沿鑒識結(jié)果是一種定性或半定量解說結(jié)果，與定量解說結(jié)果有必然差別，鑒識結(jié)果可作為邊沿地點定量反演的初值。(1)垂導(dǎo)游數(shù)：垂導(dǎo)游數(shù)方法利用零值地點確立地質(zhì)體的邊沿地點，重力異樣能夠直接使用，對磁力異樣必然轉(zhuǎn)變?yōu)榇旁粗亓Ξ悩踊蚧?/p>

3、極磁力異樣才能夠使用VDR(,xy,)zg(x,y,z)z(1.1)（2）分析信號振幅：分析信號振幅也是利用極大值地點來確立地質(zhì)體的邊沿地點合用于重、磁力異樣ASMTHDR2VDR2(1.2)（3）總水平導(dǎo)數(shù)（THDR）g(x,y,z)2g(x,y,z)2THDR(x,y,z)y(1.3)x2、輸入/輸出數(shù)據(jù)格式設(shè)計依照上述原理，此刻對上述各樣邊沿鑒識方法進(jìn)行程序設(shè)計。2.1輸入數(shù)據(jù)格式設(shè)計本次實驗給了正演的重力異樣數(shù)據(jù)，為.GRD格式，均為實型變量。比方：DSAA201201-1000.0000001000.000000-1000.0000001000.0000005.549671E-012

4、3.5398465.549671E-015.634658E-015.721339E-015.808522E-015.897312E-015.987253E-016.078691E-016.171604E-012.2輸出數(shù)據(jù)格式設(shè)計計算結(jié)果輸出數(shù)據(jù)格式與輸入格式對應(yīng)，格式為.GRD格式，均為實型變量。比方：DSAA201201-1000.0001000.000-1000.0001000.000-0.14650840.3190881-3.6523044E-02-3.3485338E-02-3.3061244E-02-3.2748722E-02-3.2688729E-02-3.2654848E-02

5、-3.2723978E-02-3.2787599E-02-3.2927759E-02-3.3138681E-022.3參數(shù)文件數(shù)據(jù)格式設(shè)計將以上部重量保留在一個文件中，該文件名變量為cmdfile，字符串變量，長度不超出80，全路徑名。在該文件中保留的參數(shù)以下：輸入數(shù)據(jù)文件名input_file，字符串變量，長度不超出80；輸出vdr數(shù)據(jù)文件名output_file_vdr，字符串變量，長度不超出80；輸出thdr數(shù)據(jù)文件名output_file_thdr，字符串變量，長度不超出80;輸出asm數(shù)據(jù)文件名output_file_asm，字符串變量，長度不超出80factor_m：擴邊比率因子，

6、實型變量（1）；3.整體設(shè)計此次程序采納IPO結(jié)構(gòu)設(shè)計，第一經(jīng)過讀取cmd文件，獲得有關(guān)輸入?yún)?shù)：輸入數(shù)據(jù)文件名gravity.grd、輸出vdr文件名field_vrd.grd、輸出thdr文件名field_thdr.grd、輸出asm文件名field_asm.grd、擴邊比率因子factor_m；此后確立確立擴邊網(wǎng)格的大小，擴邊數(shù)據(jù)點號地點；再從觀察面位場數(shù)據(jù)文件中讀取數(shù)據(jù)。下一步，進(jìn)行二維余弦擴邊，將擴完邊的數(shù)據(jù)進(jìn)行迅速二維傅里葉變換，變換到頻次域；接下來在頻次域求出在x,y,z方向的導(dǎo)數(shù)并反變換；最后求出VDR、THDR、ASM數(shù)據(jù)。最后去除擴邊部分后輸出。整體設(shè)計見表1。輸入?yún)?shù)：輸

7、入數(shù)據(jù)文件名gravity.grd、輸出vdr文件名gravity_vrd.grd、輸出thdr文件名gravity_thdr.grd、輸出asm文件名gravity_asm.grd、擴邊比率因子factor_m。確立擴邊網(wǎng)格的大小m*n(m,n均為2的冪次方)從輸入數(shù)據(jù)文件名中讀取數(shù)據(jù)對原始數(shù)據(jù)進(jìn)行二維余弦擴邊對擴邊后的數(shù)據(jù)進(jìn)行迅速二維傅里葉正變換將傅里葉變換后的數(shù)據(jù)與導(dǎo)數(shù)因子相乘求出重力數(shù)據(jù)在x,y,z方向的導(dǎo)數(shù)對導(dǎo)數(shù)進(jìn)行迅速二維傅里葉反變換分別求出VDR、THDR、ASM值。去除擴邊部分后輸出結(jié)果圖4.1整體設(shè)計N-S圖4.測試結(jié)果圖4.2重力異樣原始圖圖4.3垂導(dǎo)游數(shù)（VDR）圖4.4

8、分析信號振幅（ASM）圖4.5總水平導(dǎo)數(shù)（THDR）分析：由圖4.3可看出，VDR方法的零值線可較正確鑒識模型體的界限；由圖4.4可看出，ASM的極大值點界限可大概鑒識模型體界限，但精度不是很高。比較圖4.2到4.5能夠看出，THDR方法對模型界限的鑒識見效是最好的。結(jié)論及建議經(jīng)測試，VDR與THDR對模型體的界限地點鑒識見效較好，而ASM對模型體界限鑒識見效較差。三種方法中，THDR見效最好。附錄：邊沿鑒識程序源代碼!*!程序功能：實現(xiàn)頻次域位場導(dǎo)數(shù)運算進(jìn)行邊沿鑒識cmd文件參數(shù):cmdfile:寄存有關(guān)參數(shù)的文件名變量input_file:觀察面位場數(shù)據(jù)文件output_file_vdr：

9、場值的水平導(dǎo)數(shù)數(shù)據(jù)文件output_file_thdr:場值垂導(dǎo)游數(shù)數(shù)據(jù)文件output_file_asm:場值的分析信號振幅數(shù)據(jù)文件!factor_m：擴邊因子.grd文件參數(shù)：!N_point,N_line：點數(shù)(x方向)、線數(shù)(y方向)!x_min,x_max：x的最小值和最大值!y_min,y_max：y的最小值和最大值Ur:初始觀察面場值擴邊參數(shù)：!m1,m2：x方向?qū)嵸|(zhì)數(shù)據(jù)起點和終點點號地點!1,m：x方向擴邊后數(shù)據(jù)起點和終點點號地點!n1,n2：y方向?qū)嵸|(zhì)數(shù)據(jù)起點和終點點號地點!1,n3：y方向擴邊后數(shù)據(jù)起點和終點點號地點求導(dǎo)參數(shù)：field_re:初始觀察面信號的實部field

10、_im：初始觀察面信號的虛部Px_re:x方導(dǎo)游數(shù)信號的實部Px_im:x方導(dǎo)游數(shù)信號的虛部Py_re:y方導(dǎo)游數(shù)信號的實部Py_im:y方導(dǎo)游數(shù)信號的虛部Pz_re:z方導(dǎo)游數(shù)信號的實部Pz_im:z方導(dǎo)游數(shù)信號的虛部!W(m,n)：徑向圓頻次field_vdr:對場值作水平導(dǎo)數(shù)的結(jié)果field_thdr:對場值作垂導(dǎo)游數(shù)的結(jié)果field_asm:場值的分析信號振幅的結(jié)果!*programdeviationparameter(eigval=3.701411*1e5)character*(80)cmdfilecharacter*80input_file,output_file_vdr,outp

11、ut_file_thdr,output_file_asmreal,allocatable:field_re(:,:),field_im(:,:)real,allocatable:Px_re(:,:),Px_im(:,:),Py_re(:,:),Py_im(:,:),Pz_re(:,:),Pz_im(:,:)real,allocatable:field_vdr(:,:),field_thdr(:,:),field_asm(:,:)real,allocatable:U(:),V(:),W(:,:)integerN_point,N_lineintegerm,n,m1,m2,n1,n2realfact

12、or_mrealxmin,xmax,ymin,ymax,dx,dycmdfile=deviation.cmdcallread_cmd(cmdfile,factor_m,input_file,output_file_vdr,output_file_thdr,output_file_asm)callread_grd(input_file,N_point,N_line,Xmin,Xmax,Ymin,Ymax)callcalculate_mn(N_point,N_line,m,n,m1,m2,n1,n2,factor_m)allocate(field_re(1:m,1:n),field_im(1:m,

13、1:n)allocate(Px_re(1:m,1:n),Px_im(1:m,1:n),Py_re(1:m,1:n),Py_im(1:m,1:n),Pz_re(1:m,1:n),Pz_im(1:m,1:n)allocate(field_vdr(1:m,1:n),field_thdr(1:m,1:n),field_asm(1:m,1:n)allocate(U(1:m),V(1:n),W(1:m,1:n)callinput_grd(field_re,input_file,m1,m2,n1,n2,m,n)callexpand_cos_2D(m1,m2,m,n1,n2,n,field_re,field_

14、im)callFFT2(field_re,field_im,m,n,2)CALLcal_dxdy(xmin,xmax,ymin,ymax,N_POINT,N_LINE,dx,dy)callWAVE2D(m,n,dx,dy,U,V,W)callfactor(m,n,field_re,field_im,u,v,w,px_re,px_im,py_re,py_im,pz_re,pz_im)callFFT2(px_re,px_im,m,n,1)callFFT2(py_re,py_im,m,n,1)callFFT2(pz_re,pz_im,m,n,1)calldeviration(m,n,px_re,py

15、_re,pz_re,field_vdr,field_thdr,field_asm)callOUTPUT_GRD(field_vdr,output_file_vdr,m1,m2,n1,n2,m,n,eigval,xmin,xmax,ymin,ymax)callOUTPUT_GRD(field_thdr,output_file_thdr,m1,m2,n1,n2,m,n,eigval,xmin,xmax,ymin,ymax)callOUTPUT_GRD(field_asm,output_file_asm,m1,m2,n1,n2,m,n,eigval,xmin,xmax,ymin,ymax)deall

16、ocate(field_re,field_im,px_re,px_im,py_re,py_im,pz_re,pz_im,field_vdr,field_thdr,field_asm,u,v,w)endprogram!*!子程序：read_cmd功能：讀取參數(shù)文件輸入?yún)?shù)說明：!cmdfile：參數(shù)文件名輸出參數(shù)說明：input_file:觀察面位場數(shù)據(jù)文件output_file_vdr：對場值作水平導(dǎo)數(shù)辦理后的數(shù)據(jù)文件output_file_thdr:對場值作垂導(dǎo)游數(shù)辦理后的數(shù)據(jù)文件output_file_asm:對場值作總導(dǎo)數(shù)辦理后的數(shù)據(jù)文件!factor_m：擴邊因子!*Subroutin

17、eread_cmd(cmdfile,factor_m,input_file,output_file_vdr,output_file_thdr,output_file_asm)implicitnonecharacter*80strcharacter*(*)cmdfilecharacter*(*)input_file,output_file_vdr,output_file_thdr,output_file_asmrealfactor_mopen(10,file=cmdfile,status=old)read(10,*)str,input_fileread(10,*)str,output_file_

18、vdrread(10,*)str,output_file_thdrread(10,*)str,output_file_asmread(10,*)str,factor_mclose(10)endSubroutineread_cmd!*!子程序：read_grd功能：從原始觀察.grd文件中讀取有關(guān)參數(shù)輸入?yún)?shù)說明：!filename_obser：輸入文件名輸出參數(shù)說明：!N_point,N_line：點數(shù)、線數(shù)!x_min,x_max：x的最小值和最大值!y_min,y_max：y的最小值和最大值!*subroutineread_grd(input_file,N_point,N_line,Xmi

19、n,Xmax,Ymin,Ymax)implicitnonecharacter*(*)input_fileintegerN_point,N_linerealXmin,Xmax,Ymin,Ymaxopen(10,file=input_file,status=old)Read(10,*)Read(10,*)N_line,N_pointRead(10,*)Xmin,XmaxRead(10,*)Ymin,YmaxClose(10)endsubroutineread_grd!*!子程序：calculate_mn功能：確立擴邊數(shù)據(jù)點號地點輸入?yún)?shù)說明：!factor_m:擴邊比率因子（1.0）!a,b：點數(shù)

20、、線數(shù)輸出參數(shù)說明：!m1,m2:x方向?qū)嵸|(zhì)數(shù)據(jù)起點地點和終點地點點號!m:擴邊后數(shù)據(jù)終點地點點號(起點地點點號為1)!n1,n2:y方向?qū)嵸|(zhì)數(shù)據(jù)起點地點和終點地點點號!n:擴邊后數(shù)據(jù)終點地點點號(起點地點點號為1)!*subroutinecalculate_mn(a,b,m,n,m1,m2,n1,n2,factor_m)implicitnoneintegera,b,m,n,m1,m2,n1,n2integermtemp,mu,nurealfactor_mmtemp=aDOWHILE(mod(mtemp,2).eq.0).and.(mtemp.ne.0)mtemp=mtemp/2EnddoIF

21、(mtemp.eq.1)THENm=a*2ELSEmu=int(log(float(a)/0.693147+factor_m)m=2*muENDIFm1=1+(m-a)/2m2=m1+a-1write(*,*)m,apausemtemp=bDOWHILE(mod(mtemp,2).eq.0).and.(mtemp.ne.0)mtemp=mtemp/2EnddoIF(mtemp.eq.1)THENn=b*2ELSEnu=int(log(float(b)/0.693147+factor_m)n=2*nuENDIFn1=1+(n-b)/2n2=n1+b-1write(*,*)m1,m2,n1,n2,

22、m,npauseendsubroutinecalculate_mn!*!子程序：INPUT_GRD功能：讀取grd文件中的數(shù)據(jù)輸入?yún)?shù)說明：!filename_obser：輸入文件名!m1,m2:x方向?qū)嵸|(zhì)數(shù)據(jù)起點地點和終點地點點號!m:擴邊后數(shù)據(jù)終點地點點號(起點地點點號為1)!n1,n2:y方向?qū)嵸|(zhì)數(shù)據(jù)起點地點和終點地點點號!n:擴邊后數(shù)據(jù)終點地點點號(起點地點點號為1)輸出參數(shù)說明：!A：寄存輸出數(shù)據(jù)的二維數(shù)組名!*SUBROUTINEINPUT_GRD(A,input_file,m1,m2,n1,n2,m,n)character*(*)input_fileintegerm1,m2,n1

23、,n2,m,nrealxmin,xmax,ymin,ymaxrealA(1:m,1:n)reali,j,kdoj=1,n,1doi=1,m,1A(i,j)=3.701411*1e10enddoenddoOpen(20,file=input_file,status=old)read(20,*)read(20,*)read(20,*)xmin,xmaxread(20,*)ymin,ymaxread(20,*)read(20,*)(A(i,j),i=m1,m2),j=n1,n2)Close(20)ENDSUBROUTINEINPUT_GRD!*!子程序：expand_cos_2D功能：二維擴邊子程序

24、并為信號虛部賦值輸入?yún)?shù)說明：!m1,m2:x方向?qū)嵸|(zhì)數(shù)據(jù)起點地點和終點地點點號!m:擴邊后數(shù)據(jù)終點地點點號(起點地點點號為1)!n1,n2:y方向?qū)嵸|(zhì)數(shù)據(jù)起點地點和終點地點點號!n:擴邊后數(shù)據(jù)終點地點點號(起點地點點號為1)Ur:初始觀察面信號的實部!Ui：初始觀察面信號的虛部輸出參數(shù)說明：Ur:初始觀察面信號的實部!Ui：初始觀察面信號的虛部!*subroutineexpand_cos_2D(m1,m2,m,n1,n2,n,Ur,Ui)implicitnoneintegerm,n,m1,m2,n1,n2realUr(1:m,1:n),Ui(1:m,1:n)real,allocatable:

25、u(:),r(:)integerj,i,kallocate(u(1:m)doj=n1,n2,1doi=1,m,1u(i)=Ur(i,j)enddocallexpand_cos_1d(1,m1,m2,m,u(1)doi=1,m,1Ur(i,j)=u(i)enddoenddodeallocate(u)allocate(r(1:n)doi=1,m,1doj=1,n,1r(j)=Ur(i,j)enddocallexpand_cos_1d(1,n1,n2,n,r(1)doj=1,n,1Ur(i,j)=r(j)enddoenddodeallocate(r)doi=1,mdoj=1,nUi(i,j)=0en

26、ddoenddoendsubroutineexpand_cos_2D!*!子程序：expand_cos_1d功能：一維擴邊子程序輸入?yún)?shù)說明：!n0,n3：擴邊后數(shù)據(jù)起點地點和終點地點!n1,n2：實質(zhì)數(shù)據(jù)起點地點和終點地點!feild(i),(i=n1,n1+1,.,n2)：實質(zhì)數(shù)據(jù)輸出參數(shù)說明：!field(i),(i=n0,.,n1-1)：擴邊后左側(cè)的數(shù)據(jù)!field(i),(i=n2+1,.,n3)：擴邊后右側(cè)的數(shù)據(jù)!*Subroutineexpand_cos_1d(n0,n1,n2,n3,Field)RealField(n0:n3)pi=3.141592654Field(n0)=(F

27、ield(n1)+Field(n2)/2.0Field(n3)=Field(n0)i1=n0i2=n1DOi=i1,i2-1,1Field(i)=Field(i1)+cos(pi/2.0*(i2-i)/(i2-i1)*(Field(i2)-Field(i1)Enddoi1=n2i2=n3DOi=i1+1,i2,1Field(i)=Field(i1)+cos(pi/2.0*(i2-i)/(i2-i1)*(Field(i2)-Field(i1)EnddoEndsubroutineexpand_cos_1d!*!功能：FFT2功能：復(fù)數(shù)組2-D迅速Fourier變換輸入?yún)?shù)說明：!m0,m3：x方向

28、的起點和終點!n0,n3：y方向的起點和終點!field：輸入信號(需要賦值給Freal，實部)!m,n：x,y方向擴邊后數(shù)據(jù)終點點號地點（初步點號為1）!NF:正、反變換標(biāo)記量(1:反變換;2:正變換)輸出參數(shù)說明：!Freal：信號的實部!Fimage：信號的虛部(關(guān)于實信號而言，賦值為0)對應(yīng)頻次散布說明：!數(shù)據(jù)Freal(m,n)和Fimage(m,n)對應(yīng)的頻次散布地點為:!m方向：0,1,m/2-1,m/2,-(m/2-1),-1!n方向：0,1,n/2-1,n/2,-(n/2-1),-1!*SUBROUTINEFFT2(Freal,Fimage,m,n,nf)implicitno

29、neINTEGERm,n,nfREALFreal(1:m,1:n),Fimage(1:m,1:n)real,ALLOCATABLE:Treal(:),Timage(:)integernmmax,ierr,i,jnmmax=max(m,n)allocate(Treal(1:nmmax),Timage(1:nmmax),STAT=ierr)if(ierr.ne.0)STOPDOi=1,m,1IF(n.ne.1)THENdoj=1,n,1Treal(j)=Freal(i,j)Timage(j)=Fimage(i,j)enddocallFFT(Treal,Timage,n,nf)doj=1,n,1Fr

30、eal(i,j)=Treal(j)Fimage(i,j)=Timage(j)enddoENDIFENDDODOj=1,n,1IF(m.ne.1)THENdoi=1,m,1Treal(i)=Freal(i,j)Timage(i)=Fimage(i,j)enddocallFFT(Treal,Timage,m,nf)doi=1,m,1Freal(i,j)=Treal(i)Fimage(i,j)=Timage(i)enddoENDIFENDDOdeallocate(Treal,Timage,STAT=ierr)ENDSUBROUTINEFFT2!*!子程序：FFT功能：復(fù)數(shù)組1-D迅速Fourier變

31、換輸入?yún)?shù)說明：!Xreal(n):輸入數(shù)據(jù)實部!Ximage(n):輸入數(shù)據(jù)虛部!N:點數(shù)（N必然為2的冪次方）!NF:正、反變換標(biāo)記量(1:反變換;2:正變換)輸出參數(shù)說明：!Xreal(n):變換后頻譜實部!Ximage(n):變換后頻譜虛部對應(yīng)頻次散布說明：!數(shù)據(jù)Xreal(n)和Ximage(n)對應(yīng)的頻次散布地點為!0,1,n/2-1,n/2,-(n/2-1),-1:!*SUBROUTINEFFT(Xreal,Ximage,n,nf)implicitnoneINTEGERn,nfREALXreal(1:n),Ximage(1:n)integernu,n2,nu1,k,k1,k1n2

32、,l,i,ibitrrealf,p,arg,c,s,treal,timagenu=int(log(float(n)/0.693147+0.001)n2=n/2nu1=nu-1f=float(-1)*nf)k=0DOl=1,nu,1DOwhile(k.lt.n)doi=1,n2,1p=ibitr(k/2*nu1,nu)arg=6.2831853*p*f/float(n)c=cos(arg)s=sin(arg)k1=k+1k1n2=k1+n2treal=Xreal(k1n2)*c+Ximage(k1n2)*stimage=Ximage(k1n2)*c-Xreal(k1n2)*sXreal(k1n2

33、)=Xreal(k1)-trealXimage(k1n2)=Ximage(k1)-timageXreal(k1)=Xreal(k1)+trealXimage(k1)=Ximage(k1)+timagek=k+1enddok=k+n2ENDDOk=0nu1=nu1-1n2=n2/2ENDDODOk=1,n,1i=ibitr(k-1,nu)+1if(i.gt.k)thentreal=Xreal(k)timage=Ximage(k)Xreal(k)=Xreal(i)Ximage(k)=Ximage(i)Xreal(i)=trealXimage(i)=timageendifENDDOIF(nf.ne.

34、1)THENdoi=1,n,1Xreal(i)=Xreal(i)/float(n)Ximage(i)=Ximage(i)/float(n)enddoENDIFENDSUBROUTINEFFTFUNCTIONIBITR(J,NU)implicitnoneintegeribitr,j,nuintegerj1,itt,i,j2j1=jitt=0doi=1,nu,1j2=j1/2itt=itt*2+(j1-2*j2)ibitr=ittj1=j2enddoENDFUNCTIONIBITR!*!子程序：cal_dxdy功能：計算x,y方向的點距輸入?yún)?shù)說明：!x_min,x_max：x的最小值和最大值!y

35、_min,y_max：y的最小值和最大值!N_point,N_line：點數(shù)(x方向)、線數(shù)(y方向)輸出參數(shù)說明：!dx,dy：x,y方向的點距!*subroutinecal_dxdy(xmin,xmax,ymin,ymax,N_POINT,N_LINE,dx,dy)implicitnonerealxmin,xmax,ymin,ymaxintegerN_POINT,N_LINErealdx,dydx=(xmax-xmin)/N_POINTdy=(ymax-ymin)/N_LINEendsubroutinecal_dxdy!*!子程序：WAVE2D功能：計算2D徑向圓頻次W輸入?yún)?shù)說明：!dx

36、:x方向點距!dy:y方向線距!m:點數(shù)（M必然為2的冪次方）!n:線數(shù)（N必然為2的冪次方）輸出參數(shù)說明：!W(m,n):徑向圓頻次!*SUBROUTINEWAVE2D(m,n,dx,dy,U,V,W)implicitnoneINTEGERm,nREALdx,dyREALW(1:m,1:n),U(1:m),V(1:n)realpi,delx,delyintegermidm,midn,i,j,xx,yymidm=m/2+1midn=n/2+1delx=float(m)/dxdely=float(n)/dydoj=1,n,1yy=jif(yy.gt.midn)yy=yy-nv(j)=dely*(

37、yy-1)enddodoi=1,m,1xx=iif(xx.gt.midm)xx=xx-mu(i)=delx*(xx-1)enddodoj=1,n,1doi=1,m,1w(i,j)=sqrt(u(i)*u(i)+v(j)*v(j)enddoenddoENDSUBROUTINEWAVE2D!*!子程序：factor功能：計算x,y,z方導(dǎo)游數(shù)的實部和虛部!輸入?yún)?shù)說明：!m:點數(shù)（M必然為2的冪次方）n:線數(shù)（N必然為2的冪次方）field_re:初始觀察面信號的實部field_im：初始觀察面信號的虛部!W(m,n)：徑向圓頻次輸出參數(shù)說明：Px_re:x方導(dǎo)游數(shù)信號的實部Px_im:x方導(dǎo)游數(shù)

38、信號的虛部Py_re:y方導(dǎo)游數(shù)信號的實部Py_im:y方導(dǎo)游數(shù)信號的虛部Pz_re:z方導(dǎo)游數(shù)信號的實部Pz_im:z方導(dǎo)游數(shù)信號的虛部!*subroutinefactor(m,n,field_re,field_im,u,v,w,px_re,px_im,py_re,py_im,pz_re,pz_im)implicitnoneintegerm,nrealfield_re(1:m,1:n),field_im(1:m,1:n)realpx_re(1:m,1:n),px_im(1:m,1:n),py_re(1:m,1:n),py_im(1:m,1:n),pz_re(1:m,1:n),pz_im(1:

39、m,1:n)realU(1:m),V(1:n),W(1:m,1:n)integeri,jrealpipi=3.1415926doi=1,m,1doj=1,n,1px_re(i,j)=field_im(i,j)*u(i)*(-1)*pi*2.0px_im(i,j)=field_re(i,j)*u(i)*pi*2.0py_re(i,j)=field_im(i,j)*v(j)*(-1)*pi*2.0py_im(i,j)=field_re(i,j)*v(j)*pi*2.0pz_re(i,j)=field_re(i,j)*w(i,j)*pi*2.0pz_im(i,j)=field_im(i,j)*w(i,j)*pi*2.0enddoenddoendsubroutinefactor!*!子程序：deviration功能：計算異樣的水平導(dǎo)數(shù)、垂導(dǎo)游數(shù)及總導(dǎo)數(shù)!輸入?yún)?shù)說明：!m:點數(shù)（M必然為2的冪次方）n:線數(shù)（N必然為2的冪次方）Px_re:x方導(dǎo)游數(shù)信號的實部Px_im:x方導(dǎo)游數(shù)信號的虛部Py_re:y方導(dǎo)游數(shù)信號的實部Py_im:y方導(dǎo)游數(shù)信號的虛部Pz_re: