实验二Word格式.docx
《实验二Word格式.docx》由会员分享，可在线阅读，更多相关《实验二Word格式.docx（16页珍藏版）》请在冰豆网上搜索。
对其进行变量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（:
）
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:
n3）,Fimage（m0:
n3）
&
Fconti_real（m0:
n3）,Fconti_image（m0:
n3））
allocate（U（m0:
m3）,V（n0:
n3）,w（m0:
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='
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）
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
subroutinecalculate_n1n2（line,n0,n1,n2,n3,ymax,ymin）
integerline
integern0,n1,n2,n3,n,k
n1=1
n2=ymax-ymin+1
n=line
n0=1
n1=（2**k-n）/2
n2=n1+n-1
n3=2**k
subroutineinput_GRD（input_filename,m0,m1,m2,m3,n0,n1,n2,n3,Field
integerm0,m1,m2,m3,n0,n1,n2,n3
realField（m0:
realeigval
open（33,file=input_filename）
Field=eigval
read（33,*）
doj=n1,n2
read（33,*）（Field（i,j）,i=m1,m2）
enddo
close（33）
subroutineexpend_2D（Field,m0,m1,m2,m3,n0,n1,n2,n3）
realpi
pi=3.141592654
Field（m0,j）=（Field（m1,j）+Field（m2,j））*0.5
Field（m3,j）=Field（m0,j）
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）
doi=m0,m3
Field（i,n0）=（Field（i,n1）+Field（i,n2））*0.5
Field（i,n3）=Field（i,n0）
doj=n0,n1-1
Field（i,j）=（Field（i,n1）-Field（i,n0））*（cos（0.5*pi*（n1-j）/（n1-
n0）））+Field（i,n0）
doj=
升级会员 