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')
（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
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
重磁实验报告
实验名称：频率域位场处理与转换
*******
专业名称：勘查技术与工程
学生学号：************
指导老师：巨鹏 王万银
一
（1）实验内容：
对一个网格化数据进行向上延拓，得到延拓后的结果，并计算延拓后的均方误差。
（2）要求：
进行软件设计、代码编写、软件测试、编写实验报告。
二．基本原理
分别介绍空间域和频率域位场延拓的计算公式。
七．源程序及代码
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
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')
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)
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 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))
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)
! 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)
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
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))
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)
! 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)
Deallocate(Field,Freal,Fimage,U,V,W,Fcount_Real,Fcount_image)
END PROGRAM
!……………………………………
!搜寻未打开的文件的通道号
!……………………………………
!SUBROUTINE SEARCH_UNIT(NUNIT_START,NUNIT_IN)
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) !计算延拓因子
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)
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 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)
read(11,*) height
read(11,*) input_filename
read(11,*) Expound_a20_2D_filename
read(11,*) output_filename
read(11,*) eigval
close(11)
END SUBROUTINE