实验二.docx

实验二.docx

《实验二.docx》由会员分享，可在线阅读，更多相关《实验二.docx（16页珍藏版）》请在冰豆网上搜索。

实验二

一、基本原理

1.1空间域延拓原理

延拓分为：

向上延拓和向下延拓。

向上延拓是只将实测平面上的数据换到平面之上的平面上的计算。

对于二度体（令z向下为正）：

U（x,z）=-（z<0）

（1）

其中U（，0）为剖面上各点的实测值。

Z为延拓的高度，即转换后的平面和观测平面的距离，由于z坐标向下为正，因此z<0。

空间域的延拓实际是积分的计算。

向下延拓是由实测磁场向磁源方向延拓。

其计算公式和

（1）相同。

1.2波数域延拓原理

设场源位于z=H平面以下（H>0），则磁场为在z=H平面以上对x、y、z连续的函数，若z=0观测平面的磁场T（x，y，0）为已知，则由外部的狄利克莱问题：

T（x，y，z）=dd

（2）

对其进行变量x，y进行傅里叶变换成

（3）

因此：

=（4）

在波数域中，延拓变成了对观测数据的傅里叶变换乘以延拓因子。

二、输入/输出数据格式设计

1、输入数据

（1）控制变量名：

’cmd.txt’

（2）观测面高度之差：

height=3.3m或height=-3.3m

（3）低高度网格化数据：

从文件“A20_mag.grd”输入

（4）高高度网格化数据：

从文件“A53_mag.grd”输入

（5）特征值：

（6）输出文件名：

out.grd

2、输出数据

通过修改参数，运行两次程序，向上延拓结果及向下延拓结果均输入“output.grd”

先做向上延拓，将结果利用surfer画图后，再做向下延拓

三、总体设计

1、程序流程图

图3.1程序流程框

图3.1程序流程框图

2、参数说明

mpoint,line点数，线数

height观测面高度差

xmin,xmax,ymax,yminx方向最小值、最大值，y方向最大值、最小值

eigval特征值

m0,m1,m2,m3x方向扩边点位

n0,n1,n2,n3y方向扩边点位

Field场值

Field_real,Field_image傅里叶变换中场值的实部，虚部

Fconti_realFconti_image延拓因子实部、虚部

四、测试结果

1、原场值图：

图4.1低高度场值图

图4.2高高度场值

2、延拓结果图：

图4.3向上延拓结果

图4.4向下延拓结果

分析：

将底高度向上延拓3.3m后结果与高高度原场值图形近似，但是场值变小了，且异常部位相对不明显了，将高高度向下延拓3.3m后，异常部位没变，但是场值绝对值变大了近5倍，且在边缘部位有几个异常点。

五、结论及建议

结论：

在本次实验中，通过编写位场延拓程序代码，对位场延拓的原理及步骤有了进一步的了解。

但是在程序编写时还是遇到了很多的问题，整个程序写完，调试后出来的结果根本不是所需要的，检查后发现是延拓公式写错了，改正之后还是不对，又请同学帮忙修改后才作出结果。

这反映出我编写程序是对原理的理解不清晰，编写公式时的不细心，这都是需要改正的。

建议:

本次实验只做了一组向上延拓及向下延拓，对比性不是很强，我认为应该做三到四组不同范围的范围，这样能更好的对比延拓在什么范围内延拓结果最好。

附录

源程序代码：

ProgramWavenumber_continuation

Integermpoint,line

Realheight,xmin,xmax,ymin,ymax,eigval

Integerm0,m1,m2,m3,n0,n1,n2,n3

character*80input_filename,output_filename,cmdfile

real,allocatable:

:

Field（:

:

）,Freal（:

:

）,Fimage（:

:

）,

&Fconti_real（:

:

）,Fconti_image（:

:

）

real,allocatable:

:

U（:

）,V（:

）,W（:

:

）

cmdfile='cmd.txt'

callreadcmd（cmdfile,height,input_filename,output_filename,eigval）

callreadmn（input_filename,mpoint,line,xmin,xmax,ymin,ymax,dx,dy）

callcalculate_m1m2（mpoint,m0,m1,m2,m3,xmax,xmin）

callcalculate_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_real（m0:

m3,n0:

n3）,Fconti_image（m0:

m3,n0:

n3））

allocate（U（m0:

m3）,V（n0:

n3）,w（m0:

m3,n0:

n3））

callinput_GRD（input_filename,m0,m1,m2,m3,n0,n1,n2,n3,Field,

&eigval）

callexpend_2D（Field,m0,m1,m2,m3,n0,n1,n2,n3）

Freal=Field

Fimage=0.0

callFFT2（Freal,Fimage,m3-m0+1,n3-n0+1,2）

callwave2（m0,m3,dx,U,n0,n3,dy,V,W）

callfactor_conti（height,W,Freal,Fimage,m0,m3,n0,n3,Fconti_real,

&Fconti_image）

callmulti_sub（Freal,Fimage,Fconti_real,Fconti_image,m0,m3,n0,n3）

callFFT2（Freal,Fimage,m3-m0+1,n3-n0+1,1）

calloutput_GRD（output_filename,m0,m1,m2,m3,n0,n1,n2,n3,xmin,xmax,

&ymin,ymax,Freal,eigval）

deallocate（Field,U,V,W,Freal,Fimage,Fconti_real,Fconti_image）

endprogram

subroutinereadcmd（cmdfile,height,input_filename,output_filename

&,eigval）

character*80cmdfile,input_filename,output_filename

realheigt,eigval

open（11,file=cmdfile,status='old'）

read（11,*）height

read（11,*）input_filename,output_filename

read（11,*）eigval

close（11）

endsubroutine

subroutinereadmn（input_filename,mpoint,line,xmin,xmax,ymin,ymax,

&dx,dy）

character*80input_filename

integermpoint,line

realxmin,xmax,ymin,ymax,dx,dy

open（22,file=input_filename,status='old'）

read（22,*）

read（22,*）mpoint,line

read（22,*）xmin,xmax

read（22,*）ymin,ymax

close（22）

dx=（xmax-xmin）/（mpoint-1）

dy=（ymax-ymin）/（line-1）

endsubroutine

subroutinecalculate_m1m2（mpoint,m0,m1,m2,m3,xmin,xmax）

integermpoint,n,k

integerm0,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

endsubroutine

subroutinecalculate_n1n2（line,n0,n1,n2,n3,ymax,ymin）

integerline

integern0,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

endsubroutine

subroutineinput_GRD（input_filename,m0,m1,m2,m3,n0,n1,n2,n3,Field

&,eigval）

integerm0,m1,m2,m3,n0,n1,n2,n3

character*80input_filename

realField（m0:

m3,n0:

n3）

realeigval

open（33,file=input_filename）

Field=eigval

read（33,*）

read（33,*）

read（33,*）

read（33,*）

read（33,*）

doj=n1,n2

read（33,*）（Field（i,j）,i=m1,m2）

enddo

close（33）

endsubroutine

subroutineexpend_2D（Field,m0,m1,m2,m3,n0,n1,n2,n3）

integerm0,m1,m2,m3,n0,n1,n2,n3

realpi

realField（m0:

m3,n0:

n3）

pi=3.141592654

doj=n1,n2

Field（m0,j）=（Field（m1,j）+Field（m2,j））*0.5

Field（m3,j）=Field（m0,j）

enddo

doj=n1,n2

doi=m0+1,m1-1

Field（i,j）=（Field（m1,j）-Field（m0,j））*（cos（0.5*pi*（m1-i）/（m1-m0）

&））+Field（m1,j）

enddo

doi=m2+1,m3-1

Field（i,j）=（Field（m3,j）-Field（m2,j））*（cos（0.5*pi*（m3-i）/（m3-m2

&）））+Field（m2,j）

enddo

enddo

doi=m0,m3

Field（i,n0）=（Field（i,n1）+Field（i,n2））*0.5

Field（i,n3）=Field（i,n0）

enddo

doi=m0,m3

doj=n0,n1-1

Field（i,j）=（Field（i,n1）-Field（i,n0））*（cos（0.5*pi*（n1-j）/（n1-

&n0）））+Field（i,n0）

enddo

doj=