实验二.docx

1、实验二 一、基本原理1.1空间域延拓原理延拓分为：向上延拓和向下延拓。向上延拓是只将实测平面上的数据换到平面之上的平面上的计算。对于二度体（令z向下为正）： U(x,z)=- （z0） （1）其中U（，0）为剖面上各点的实测值。Z为延拓的高度，即转换后的平面和观测平面的距离，由于z坐标向下为正，因此z0），则磁场为在z=H平面以上对x、y、z连续的函数，若z=0观测平面的磁场T（x，y，0）为已知，则由外部的狄利克莱问题： T（x，y，z）=dd （2） 对其进行变量x，y进行傅里叶变换成 （3） 因此：= （4）在波数域中，延拓变成了对观测数据的傅里叶变换乘以延拓因子。二、输入/输出数据格式

2、设计1、输入数据 （1）控制变量名：cmd.txt （2）观测面高度之差：height=3.3m或height=-3.3m （3）低高度网格化数据：从文件“A20_mag.grd”输入 （4）高高度网格化数据：从文件“A53_mag.grd”输入 （5）特征值： （6）输出文件名：out.grd2、输出数据通过修改参数，运行两次程序，向上延拓结果及向下延拓结果均输入“output.grd”先做向上延拓，将结果利用surfer画图后，再做向下延拓三、总体设计1、程序流程图图3.1 程序流程框图3.1 程序流程框图2、参数说明mpoint,line 点数，线数height 观测面高度差xmin,x

3、max,ymax,ymin x方向最小值、最大值，y方向最大值、最小值eigval 特征值m0,m1,m2,m3 x方向扩边点位n0,n1,n2,n3 y方向扩边点位Field 场值Field_real,Field_image 傅里叶变换中场值的实部，虚部Fconti_real Fconti_image 延拓因子实部、虚部四、测试结果1、原场值图：图4.1 低高度场值图图4.2 高高度场值2、延拓结果图：图4.3 向上延拓结果图4.4 向下延拓结果分析：将底高度向上延拓3.3m后结果与高高度原场值图形近似，但是场值变小了，且异常部位相对不明显了，将高高度向下延拓3.3m后，异常部位没变，但是场

4、值绝对值变大了近5倍，且在边缘部位有几个异常点。五、结论及建议 结论：在本次实验中，通过编写位场延拓程序代码，对位场延拓的原理及步骤有了进一步的了解。但是在程序编写时还是遇到了很多的问题，整个程序写完，调试后出来的结果根本不是所需要的，检查后发现是延拓公式写错了，改正之后还是不对，又请同学帮忙修改后才作出结果。这反映出我编写程序是对原理的理解不清晰，编写公式时的不细心，这都是需要改正的。 建议:本次实验只做了一组向上延拓及向下延拓，对比性不是很强，我认为应该做三到四组不同范围的范围，这样能更好的对比延拓在什么范围内延拓结果最好。附录源程序代码： Program Wavenumber_conti

5、nuation Integer mpoint,line Real height,xmin,xmax,ymin,ymax,eigval Integer m0,m1,m2,m3,n0,n1,n2,n3 character*80 input_filename,output_filename,cmdfile real,allocatable :Field(:,:),Freal(:,:),Fimage(:,:), &Fconti_real(:,:),Fconti_image(:,:) real,allocatable :U(:),V(:),W(:,:) cmdfile=cmd.txt call read

6、cmd(cmdfile,height,input_filename,output_filename,eigval) call readmn(input_filename,mpoint,line,xmin,xmax,ymin,ymax,dx,dy) call calculate_m1m2(mpoint,m0,m1,m2,m3,xmax,xmin) call calculate_n1n2(mpoint,n0,n1,n2,n3,ymax,ymin) allocate(Field(m0:m3,n0:n3),Freal(m0:m3,n0:n3),Fimage(m0:m3,n0:n3) &,Fconti_

7、real(m0:m3,n0:n3),Fconti_image(m0:m3,n0:n3) allocate(U(m0:m3),V(n0:n3),w(m0:m3,n0:n3) call input_GRD(input_filename,m0,m1,m2,m3,n0,n1,n2,n3,Field, &eigval) call expend_2D(Field,m0,m1,m2,m3,n0,n1,n2,n3) Freal=Field Fimage=0.0 call FFT2(Freal,Fimage,m3-m0+1,n3-n0+1,2) call wave2(m0,m3,dx,U,n0,n3,dy,V,

8、W) call factor_conti(height,W,Freal,Fimage,m0,m3,n0,n3,Fconti_real, &Fconti_image) call multi_sub(Freal,Fimage,Fconti_real,Fconti_image,m0,m3,n0,n3) call FFT2(Freal,Fimage,m3-m0+1,n3-n0+1,1) call output_GRD(output_filename,m0,m1,m2,m3,n0,n1,n2,n3,xmin,xmax, &ymin,ymax,Freal,eigval) deallocate(Field,

9、U,V,W,Freal,Fimage,Fconti_real,Fconti_image) end program subroutine readcmd(cmdfile,height,input_filename,output_filename &,eigval) character*80 cmdfile,input_filename,output_filename real heigt,eigval open(11,file=cmdfile,status=old) read(11,*) height read(11,*) input_filename,output_filename read(

10、11,*) eigval close(11) end subroutine subroutine readmn(input_filename,mpoint,line,xmin,xmax,ymin,ymax, &dx,dy) character*80 input_filename integer mpoint,line real xmin,xmax,ymin,ymax,dx,dy open(22,file=input_filename,status=old) read(22,*) read(22,*) mpoint,line read(22,*) xmin,xmax read(22,*) ymi

11、n,ymax close(22) dx=(xmax-xmin)/(mpoint-1) dy=(ymax-ymin)/(line-1) end subroutine subroutine calculate_m1m2(mpoint,m0,m1,m2,m3,xmin,xmax) integer mpoint,n,k integer m0,m1,m2,m3 m1=1.0 m2=xmax-xmin+1 n=mpoint k=int(log(n*1.0)/log(2.0)+0.5)+1 m0=1 m1=(2*k-n)/2 m2=m1+n-1 m3=2*k end subroutine subroutin

12、e calculate_n1n2(line,n0,n1,n2,n3,ymax,ymin) integer line integer n0,n1,n2,n3,n,k n1=1 n2=ymax-ymin+1 n=line k=int(log(n*1.0)/log(2.0)+0.5)+1 n0=1 n1=(2*k-n)/2 n2=n1+n-1 n3=2*k end subroutine subroutine input_GRD(input_filename,m0,m1,m2,m3,n0,n1,n2,n3,Field &,eigval) integer m0,m1,m2,m3,n0,n1,n2,n3

13、character*80 input_filename real Field(m0:m3,n0:n3) real eigval open(33,file=input_filename) Field=eigval read(33,*) read(33,*) read(33,*) read(33,*) read(33,*) do j=n1,n2 read(33,*)(Field(i,j),i=m1,m2) end do close(33) end subroutine subroutine expend_2D(Field,m0,m1,m2,m3,n0,n1,n2,n3) integer m0,m1

14、,m2,m3,n0,n1,n2,n3 real pi real Field(m0:m3,n0:n3) pi=3.141592654 do j=n1,n2 Field(m0,j)=(Field(m1,j)+Field(m2,j)*0.5 Field(m3,j)=Field(m0,j) end do do j=n1,n2 do i=m0+1,m1-1 Field(i,j)=(Field(m1,j)-Field(m0,j)*(cos(0.5*pi*(m1-i)/(m1-m0) &)+Field(m1,j) end do do i=m2+1,m3-1 Field(i,j)=(Field(m3,j)-Field(m2,j)*(cos(0.5*pi*(m3-i)/(m3-m2 &)+Field(m2,j) end do end do do i=m0,m3 Field(i,n0)=(Field(i,n1)+Field(i,n2)*0.5 Field(i,n3)=Field(i,n0) end do do i=m0,m3 do j=n0,n1-1 Field(i,j)=(Field(i,n1)-Field(i,n0)*(cos(0.5*pi*(n1-j)/(n1- &n0)+Field(i,n0) end do do j=