频率域位场处理与转换(12.25-)

CALL Multi_sub(Freal,Fimage,Fcount_Real,Fcount_image,m0,m3,n0,n3)
CALL FFt2(Freal,Fimage,m3-m0+1,n3-n0+1,1)
CALL Putout_GRD(Freal,output_filename,m0,m1,m2,m3,n0,n1,n2,n3,eigval,Xmin,Xmax,Ymin,Ymax)
ALLOCATE(U(m0:m3),V(n0:n3),W(m0:m3,n0:n3),Fcount_Real(m0:m3,n0:n3),Fcount_image(m0:m3,n0:n3))
ALLOCATE(Field(m0:m3,n0:n3),Freal(m0:m3,n0:n3),Fimage(m0:m3,n0:n3))
end do
end do
do i=m0,m3
Field(i,n0)=(Field(i,n1)+Field(i,n2))/2
Field(i,n3)=Field(i,n0)
end do
do i=m0,m3
do j=n0+1,n1-1
Field(i,j)=Field(i,n0)+cos(pi/2.0*(n1-j)/(n1-n0))*(Field(i,n1)-Field(i,n0))
四．总体设计
五．测试结果
分析：将底高度向上延拓3.3m后结果与高高度原场值图形近似，但是场值变小了，且异常部位相对不明显了。
六．结论与建议
结论：向上延拓可以压制浅部的重磁异常特征，突出深部异常。向下延拓是由实测磁源异常得到。
建议：本次实验程序编写比较困难，以后要提高编程能力。最后得到了扩边后的数据，分析能力比较差，分析结果不够准确，多看资料，加强这方面能力。
CALL Read_mn(input_filename,mpoint,line,Xmin,Xmax,Ymin,Ymax)
CALL Calculate_m1m2_dx(mpoint,m0,m1,m2,m3,m,Xmin,Xmax,dx)
CALL Calculate_m1m2_dx(line,n0,n1,n2,n3,n,Ymin,Ymax,dy)
m0=1
factor_m=1.5
dx=(Xmax-Xmin)/(mpoint-1)
mtemp=mpoint
do while((mod(mtemp,2).eq.0).and.(mtemp.ne.0))
mtemp=mtemp/2
end do
if(mtemp.eq.1) then
m=mpoint*2
else
! call search_unit(10,nunit_in)
OPEN(12,file=input_filename)
READ(12,*)
READ(12,*) mpoint,line
READ(12,*) Xmin,Xmax
READ(12,*) Ymin,Ymax
! read(12,*)
close(12)
end do源自文库
do j=n2+1,n3-1
Field(i,j)=Field(i,n2)+cos(pi/2.0*(n3-j)/(n3-n2))*(Field(i,n3)-Field(i,n2))
end do
end do
END SUBROUTINE
!*******************输出扩边Expound_a20_2D.grd***********************
END SUBROUTINE
!……………………
!计算扩边点位
!……………………
SUBROUTINE Calculate_m1m2_dx(mpoint,m0,m1,m2,m3,m,Xmin,Xmax,dx)
integer mpoint,mtemp,m0,m1,m2,m3,m
real Xmin,Xmax,dx,factor_m,mu
! end do
!END SUBROUTINE
!…………………
!读入数据文件
!…………………
SUBROUTINE Read_CMD(cmd_filename,height,input_filename,Expound_a20_2D_filename,output_filename,eigval)
real height,eigval
七．源程序及代码
PROGRAM wavenumber_continuation
INTEGER mpoint,line,m0,m1,m2,m3,n0,n1,n2,n3,m,n
REAL height,eigval,Xmin,Xmax,Ymin,Ymax,dx,dy
CHARACTER*80 input_filename,output_filename,cmd_filename,Expound_a20_2D_filename
read(13,*)
read(13,*)
read(13,*) Xmin,Xmax
read(13,*) Ymin,Ymax
read(13,*)
read(13,*) ((Field(i,j),i=m1,m2),j=n1,n2)
Close(13)
END SUBROUTINE
!**********************2-D扩边****************************
设场源位于z=H平面以下（H>0），则磁场为在z=H平面以上对x、y、z连续的函数，若z=0观测平面的磁场T（x，y，0）为已知，则由外部的狄利克莱问题：
T（x，y，z）= d d
对其进行变量x，y进行傅里叶变换成
因此：
=
在频率域中，延拓变成了对观测数据的傅里叶变换乘以延拓因子。
三
1.
（1）建立文件：wave.par
SUBROUTINE Expound_2D(Field,m0,m1,m2,m3,n0,n1,n2,n3,line)
integer m0,m1,m2,m3,n0,n1,n2,n3
real Field(m0:m3,n0:n3)
real::pi=3.141592654
do i=n1,n2
Field(m0,i)=(Field(m1,i)+Field(m2,i))/2
（2）观测面高度之差：height=3.3m
（3）低高度网格化数据：A20_mag.grd
（4）特征值：
（6）输出文件名：expound_a20_2D.grd
2.
（1）输入文件：input_filename
(2) 输出文件：out_filename
(3) 点数：mpoint
(4)线数：nline
(5)异常值：filed
character*(*) input_filename
integer m0,m1,m2,m3,n0,n1,n2,n3
real eigval
real Field(m0:m3,n0:n3)
Field=eigval
! call search_unit(10,nunit_in)
open(13,file=input_filename,form='formatted')
Freal=Field
Fimage=0.
CALL FFt2(Freal,Fimage,m3-m0+1,n3-n0+1,2) !2说明是正变换
CALL WAVE2D_sub(U,V,W,m,n,dx,dy)
CALL Factor_conti(height,m0,m3,n0,n3,U,V,W,Fcount_Real,Fcount_image) !计算延拓因子
!读入源文件中的值
!………………………
SUBROUTINE Read_mn(input_filename,mpoint,line,Xmin,Xmax,Ymin,Ymax)
integer mpoint,line
real Xmin,Xmax,Ymin,Ymax
character*80 input_filename
REAL,ALLOCATABLE::Field(:,:),Freal(:,:),Fimage(:,:)
REAL,ALLOCATABLE::U(:),v(:),w(:,:),Fcount_Real(:,:),Fcount_image(:,:)
cmd_filename='wave.par'
CALL Read_CMD(cmd_filename,height,input_filename,Expound_a20_2D_filename,output_filename,eigval)
Field(m3,i)=Field(m0,i)
end do
do j=n1,n2
do i=m0+1,m1-1
Field(i,j)=Field(m0,j)+cos(pi/2.0*(m1-i)/(m1-m0))*(Field(m1,j)-Field(m0,j))
end do
do i=m2+1,m3-1
Field(i,j)=Field(m2,j)+cos(pi/2.0*(m3-i)/(m3-m2))*(Field(m3,j)-Field(m2,j))
重磁实验报告
实验名称：频率域位场处理与转换
*******
专业名称：勘查技术与工程
学生学号：************
指导老师：巨鹏 王万银
一
（1）实验内容：
对一个网格化数据进行向上延拓，得到延拓后的结果，并计算延拓后的均方误差。
（2）要求：
进行软件设计、代码编写、软件测试、编写实验报告。
二．基本原理
分别介绍空间域和频率域位场延拓的计算公式。
! logical unit_open
! integer nunit_in,nunit_start
! nunit_in=nunit_start
! unit_open=.TRUE.
! do while (unit_open)
! nunit_in=nunit_in+1
! inquire(unit=nunit_in,opened=unit_open)
CALL Input_GRD(input_filename,m0,m1,m2,m3,n0,n1,n2,n3,Field,eigval)
CALL Expound_2D(Field,m0,m1,m2,m3,n0,n1,n2,n3,line)
CALL Putout_GRD_Expound_A20_2D(Field,Expound_a20_2D_filename,m0,m3,n0,n3,m,n,eigval,Xmin,Xmax,Ymin,Ymax)
（
延拓分为向上延拓和向下延拓。
向上延拓是只将实测平面上的数据换到平面之上的平面上的计算。对于二度体（令z向下为正）：
U(x,z)=- （z<0）
其中U（ ，0）为剖面上各点的实测值。Z为延拓的高度，即转换后的平面和观测平面的距离，由于z坐标向下为正，因此z<0。空间域的延拓实际是积分的计算。
向下延拓是由实测磁场向磁源方向延拓。
character*80 input_filename,output_filename,cmd_filename,Expound_a20_2D_filename
! call search_unit(10,nunit_in)
OPEN(11,file=cmd_filename,status='old')
Deallocate(Field,Freal,Fimage,U,V,W,Fcount_Real,Fcount_image)
END PROGRAM
!……………………………………
!搜寻未打开的文件的通道号
!……………………………………
!SUBROUTINE SEARCH_UNIT(NUNIT_START,NUNIT_IN)
read(11,*) height
read(11,*) input_filename
read(11,*) Expound_a20_2D_filename
read(11,*) output_filename
read(11,*) eigval
close(11)
END SUBROUTINE
!………………………
mu=int(log(float(mpoint))/0.693147+factor_m)
m=2**mu
end if
m1=1+(m-mpoint)/2
m2=m1+mpoint-1
m3=m
END SUBROUTINE
!………………………………
!输入源文件的重磁异常值
!………………………………
SUBROUTINE Input_GRD(input_filename,m0,m1,m2,m3,n0,n1,n2,n3,Field,eigval)