《重磁资料处理与解释》实验二频率域位场处理和转换实验学院：地测学院专业名称：勘查技术与工程学生姓名：学生学号：指导老师：提交日期：2018年1月9日二0一八年一月目录1 基本原理 (2)1.1位场的方程 (2)1.2二维傅里叶变换及卷积性质 (2)（1）傅里叶变换 (2)（2）卷积性质 (2)1.3频率域位场延拓原理 (3)2 输入/输出数据格式设计 (3)2.1 输入数据格式设计 (3)2.2 输出数据格式设计 (3)2.3 参数文件数据格式设计 (3)3 总体设计 (4)3.1频率域位场处理与转换的一般步骤 (4)3.2软件总体设计结果流程图 (4)4 测试结果 (5)4.1 测试参数 (5)（1）向上延拓 (5)（2）向下延拓 (5)4.2 测试结果 (6)5 结论及建议 (7)附录：源程序代码 (8)1 基本原理1.1位场的方程由场论知识可知，位场方程分为两大类：有源的Possion 方程()02≠∇U ，以及无源的Laplace 方程()02=∇U 。

Laplace 方程的第一边值问题()1|f U S =通常为Dirichlet 问题，第二边值问题⎪⎭⎫⎝⎛=∂∂2|f n U s 通常称为Nueman 问题。

若P 点在S 平面内称为内部问题，反之称为外部问题。

由唯一性定理可知，Dirichlet 的内部和外部问题的解是唯一的，而Nueman 内部问题的解不是唯一的，有一常数差，但其外部问题解是唯一的。

外部问题的解的唯一性的原因：。

0;0=∂∂=∞→∞→r r nUU 无源区域位场可以表示为：ds n G W n W G p W ⎰⎥⎦⎤⎢⎣⎡∂∂-∂∂=π41)( （1-1）()()()()()[]()()z y x h W d d z y x W z z y W -=-+-+--=⎰⎰+∞∞-+∞∞-ξξηεηεξηεξηεπξ,,*,,,,2,,x 23222（1-2）1.2二维傅里叶变换及卷积性质（1）傅里叶变换[]⎰⎰+∞∞-+∞∞-+-==dxdy y x g y x g F v u G e vy ux i )(2),(),(),(π （1-3）[]⎰⎰+∞∞-+∞∞-+-==dudv v u G v u G y x g eF vy ux i )(21),(),(),(π （1-4）（2）卷积性质()()[]()()v u P v u G y x p y x F ,*,,*,g = （1-5）()()[]()()y x p y x v u P v u G F ,*,g ,*,1=- （1-6）1.3频率域位场延拓原理当已知实测平面的异常时，换算场源以外的异常称之为延拓，分为向上延拓和向下延拓。

半空间狄利克莱问题解析解：()[]()[]()[]()[]()z v u ey x W F z y x h F y x W F z y x W F -+-=-=ζπζζζ222,,,,,,,, （1-7）其中：ez v u )(222-+-ζπ称为延拓因子，ζ为计算面Z 坐标 ,Z 轴向下为正方向，()[]ζ,,y x W F 为计算面频率域位场，()[]z y x W F ,,为延拓面的频率域位场。

2 输入/输出数据格式设计2.1 输入数据格式设计观测面位场数据保存在filename_obser 文件中，为.grd 格式。

计算延拓误差时的精确场值文件保存在filename_obser2中，为.grd 格式。

2.2 输出数据格式设计实际计算得到的数据保存在filename_output1和文件filename_output2中，为.grd 格式。

2.3 参数文件数据格式设计将所要读取的参数保存在一个文件中，该文件名变量为cmdfile ，字符串变量，长度不超过80，全路径名。

在该文件中保存的参数如下：filename_obser:低高度观测面位场数据文件 filename_output1：向上延拓后位场数据文件 filename_output2：向下延拓后位场数据文件 filename_obser2:高高度观测面数据文件factor_m ：扩边因子distance1:向上延拓的高度（m/z 轴向下为正方向） distance2:向下延拓的高度（m/z 轴向下为正方向）3 总体设计3.1频率域位场处理与转换的一般步骤()[][][][][]即可。

去掉扩边部分，输出步第；进行反变换得到步第；计算步第；以及方向转换因子、延拓因子计算导数因子步第；对扩边后的数据进行步第进行扩边处理；对网格化的数据步第),,(:6),,(),,(:5),,(),,(F :4),,,,( )Y(),,,,D(:3,,W FFT :2),,(:1121z y x Q z y x Q F F z y x Q f Y D y x W F z y x Q M M t t u,v f z u,v,t t t u,v y x F y x W n -=⋅⋅⋅=''-⇒ζζζζ3.2软件总体设计结果流程图此次程序采用IPO 结构设计，首先通过读取cmd 文件，得到相关输入参数：观测面位场数据文件名变量、延拓后位场数据文件名变量、延拓后准确位场数据文件名变量、扩边因子、延拓的高度（m/z 轴向下为正方向）；然后确定确定扩边网格的大小，扩边数据点号位置；再从观测面位场数据文件中读取数据。

下一步，进行二维余弦扩边，将扩完边的数据进行快速二维傅里叶变换，转换到频率域；接下来计算延拓因子并且将扩完边的数据进行快速二维傅里叶变换后在频率域与延拓因子相乘；最后进行快速二维傅里叶反变换并且去除扩边部分后输出。

总体设计见表1。

4 测试结果4.1 测试参数（1）向上延拓原始场值数据保存在’’a20_mag.grd”中，向上延拓3.3m，延拓后理论数据保存在“a53_mag.grd”中,延拓后的数据保存在output1.grd中。

网格大小为27*27（m），Xmin=-26m, Xmax=26m, Ymin=-26m, Ymax=26m.（2）向下延拓原始场值数据保存在’’a53_mag.grd”中，向下延拓3.3m，延拓后理论数据保存在“a20_mag.grd”中延拓后的数据保存在output2.grd中。

网格大小为27*27（m），Xmin=-26m, Xmax=26m, Ymin=-26m, Ymax=26m.有关参数保存在filename.cmd文件中，如下:A20_mag.grdA53_mag.grdoutput1.grdoutput2.grd3.3 -3.34.2 测试结果图1 低高度观测面(a20_mag)位场等值线图图2 高高度观测面(a53_mag)位场等值线图图3 高高度观测面(a53_mag)向下延拓位场等值线图图4 低高度观测面(a20_mag)向上延拓位场等值线图图5 运行结果图5 结论及建议由测试结果可知，向上延拓可以对浅部异常进行压制，且误差较小；向下延拓可突出浅部异常，且延拓结果误差较大。

通过本次频率域位场处理与转换设计实验，我收获颇丰，同时也感触很深！首先，刚开始我对空间域转频率域以及位场延拓的概念与处理是有不理解的，但是经过老师耐心细致的讲解和同学们之间的帮助，最终克服了这些困难。

其次，加深了我对于位场延拓的作用和具体延拓步骤的理解。

感谢老师孜孜不倦的教导与同学不厌其烦的讲解。

谢谢！附录：源程序代码!******************************************************************** **********! 程序功能：实现频率域位场延拓! cmd文件参数:! cmdfile:存放有关参数的文件名变量! filename_obser:低高度观测面位场数据文件! filename_output1：向上延拓后位场数据文件! filename_output2：向下延拓后位场数据文件! filename_obser2:高高度观测面数据文件! factor_m：扩边因子! distance1:向上延拓的高度（m/z轴向下为正方向）! distance2:向下延拓的高度（m/z轴向下为正方向）! .grd文件参数：! N_point,N_line：点数(x方向)、线数(y方向)! x_min,x_max：x的最小值和最大值! y_min,y_max：y的最小值和最大值! Ur:初始观测面场值! 扩边参数：! m1,m2：x方向实际数据起点和终点点号位置! 1,m：x方向扩边后数据起点和终点点号位置! n1,n2：y方向实际数据起点和终点点号位置! 1,n3：y方向扩边后数据起点和终点点号位置! 延拓参数：! Ur:初始观测面信号的实部! Ui：初始观测面信号的虚部! Fr：延拓观测面的信号的实部! Fi：延拓观测面的信号的虚部! U：延拓后准确场值! W(m,n)：径向圆频率! Y(m,n)：延拓因子! 误差参数：! error:延拓后场值与真实场值的均方误差!******************************************************************** **********Program plyytimplicit nonereal eigvalparameter(eigval=3.701411*1e5)character*(80)cmdfile,filename_obser,filename_output1,filename_output2,filename _obser2real,allocatable::Ur(:,:),Ui(:,:),y(:,:),Fr(:,:),Fi(:,:),U(:,:),W(:,:)integer N_point,N_lineinteger m,n,m1,m2,n1,n2integer factor_mreal xmin,xmax,ymin,ymax,dx,dyreal distance1,distance2,errorcmdfile='filename.cmd'callread_cmd(cmdfile,filename_obser,filename_output1,filename_output2,factor_m,& filename_obser2,distance1,distance2)call read_grd(filename_obser,N_point,N_line,Xmin,Xmax,Ymin,Ymax)call calculate_mn(N_point,N_line,m,n,m1,m2,n1,n2,factor_m)allocate(Ur(1:m,1:n),Ui(1:m,1:n),y(1:m,1:n),Fr(1:m,1:n),Fi(1:m,1:n),U(1:m,1:n),W(1:m ,1:n))call input_grd(Ur,filename_obser,m1,m2,n1,n2,m,n)call input_grd(U,filename_obser2,m1,m2,n1,n2,m,n)!低高度向上延拓call expand_cos_2D(m1,m2,m,n1,n2,n,Ur,Ui)call FFT2(Ur,Ui,m,n,2)CALL cal_dxdy(xmin,xmax,ymin,ymax,N_POINT,N_LINE,dx,dy)call WAVE2D(m,n,dx,dy,W)call calculate_Y(m,n,W,y,distance1)call calculate_FmulY(m,n,Ur,Ui,y,Fr,Fi)call FFT2(Fr,Fi,m,n,1)calloutput_grd(filename_output1,N_point,N_line,m,n,m1,m2,n1,n2,eigval,Fr,Xmin,Xmax ,Ymin,Ymax)call calculate_error(distance1,Fr,U,m1,m2,n1,n2,m,n,error)!高高度向下延拓call expand_cos_2D(m1,m2,m,n1,n2,n,U,Ui)call FFT2(U,Ui,m,n,2)CALL cal_dxdy(xmin,xmax,ymin,ymax,N_POINT,N_LINE,dx,dy)call WAVE2D(m,n,dx,dy,W)call calculate_Y(m,n,W,y,distance2)call calculate_FmulY(m,n,U,Ui,y,Fr,Fi)call FFT2(Fr,Fi,m,n,1)calloutput_grd(filename_output2,N_point,N_line,m,n,m1,m2,n1,n2,eigval,Fr,Xmin,Xmax ,Ymin,Ymax)call calculate_error(distance2,Fr,Ur,m1,m2,n1,n2,m,n,error)deallocate(Ur,Ui,y,Fr,Fi,u,w)end program!******************************************************************** ********! 子程序：read_cmd! 功能：读取参数文件! 输入参数说明：! cmdfile：参数文件名! 输出参数说明：! filename_obser：存放低高度观测面的数据文件! filename_output1：存放向上延拓观测面的数据文件! filename_output2：存放向下延拓观测面的数据文件! filename_obser2：存放高高度观测面的数据文件! distance1:向上延拓的高度（m/z轴向下为正方向）! distance2:向下延拓的高度（m/z轴向下为正方向）! factor_m：扩边因子!******************************************************************** ********subroutineread_cmd(cmdfile,filename_obser,filename_output1,filename_output2,factor_m,& filename_obser2,distance1,distance2)implicit nonecharacter*(*)cmdfile,filename_obser,filename_output1,filename_output2,filename_ obser2integer factor_mreal distance1,distance2character*80 stropen(10,file=cmdfile,status='old')read(10,*)filename_obserread(10,*)filename_obser2read(10,*)filename_output1read(10,*)filename_output2read(10,*)factor_mread(10,*)distance1read(10,*)distance2close(10)end subroutine read_cmd!******************************************************************** *******! 子程序：read_grd! 功能：从原始观测.grd文件中读取相关参数! 输入参数说明：! filename_obser：输入文件名! 输出参数说明：! N_point,N_line：点数、线数! x_min,x_max：x的最小值和最大值! y_min,y_max：y的最小值和最大值!******************************************************************** *******subroutine read_grd(filename_obser,N_point,N_line,Xmin,Xmax,Ymin,Ymax) implicit nonecharacter*(*)filename_obserinteger N_point,N_linereal Xmin,Xmax,Ymin,Ymaxopen(10,file=filename_obser,status='old')Read(10,*)Read(10,*)N_line,N_pointRead(10,*)Xmin,XmaxRead(10,*)Ymin,YmaxClose(10)end subroutine read_grd!******************************************************************** ******! 子程序：calculate_mn! 功能：确定扩边数据点号位置! 输入参数说明：! factor_m: 扩边比例因子（>1.0）! a,b：点数、线数! 输出参数说明：! m1,m2: x方向实际数据起点位置和终点位置点号! m:扩边后数据终点位置点号(起点位置点号为1)! n1,n2: y方向实际数据起点位置和终点位置点号! n:扩边后数据终点位置点号(起点位置点号为1)!******************************************************************** ******subroutine calculate_mn(a,b,m,n,m1,m2,n1,n2,factor_m)implicit noneinteger a,b,m,n,m1,m2,n1,n2integer mtemp,mu,nuinteger factor_mmtemp=aDO WHILE ((mod(mtemp,2).eq.0).and.(mtemp.ne.0))mtemp=mtemp/2End doIF (mtemp.eq.1) THENm=a*2ELSEmu=int(log(float(a))/0.693147+factor_m)m=2**muEND IFm1=1+(m-a)/2m2=m1+a-1pausemtemp=bDO WHILE ((mod(mtemp,2).eq.0).and.(mtemp.ne.0))mtemp=mtemp/2End doIF (mtemp.eq.1) THENn=b*2ELSEnu=int(log(float(b))/0.693147+factor_m)n=2**nuEND IFn1=1+(n-b)/2n2=n1+b-1pauseend subroutine calculate_mn!******************************************************************** *****! 子程序：INPUT_GRD! 功能：读取grd文件中的数据! 输入参数说明：! filename_obser：输入文件名! m1,m2: x方向实际数据起点位置和终点位置点号! m:扩边后数据终点位置点号(起点位置点号为1)! n1,n2: y方向实际数据起点位置和终点位置点号! n:扩边后数据终点位置点号(起点位置点号为1)! 输出参数说明：! A：存放输出数据的二维数组名!******************************************************************** *****SUBROUTINE INPUT_GRD(A,filename_obser,m1,m2,n1,n2,m,n)character*(*)filename_obserinteger m1,m2,n1,n2,m,nreal xmin,xmax,ymin,ymaxreal A(1:m,1:n)real i,j,kdo j=1,n,1do i=1,m,1A(i,j)=3.701411*1e10enddoenddoOpen(20,file=filename_obser,status='old')read(20,*)read(20,*)read(20,*)xmin,xmaxread(20,*)ymin,ymaxread(20,*)read(20,*) ((A(i,j),i=m1,m2),j=n1,n2)Close(20)END SUBROUTINE INPUT_GRD!******************************************************************** *****! 子程序：expand_cos_2D! 功能：二维扩边子程序并为信号虚部赋值! 输入参数说明：! m1,m2: x方向实际数据起点位置和终点位置点号! m:扩边后数据终点位置点号(起点位置点号为1)! n1,n2: y方向实际数据起点位置和终点位置点号! n:扩边后数据终点位置点号(起点位置点号为1)! Ur:初始观测面信号的实部! Ui：初始观测面信号的虚部! 输出参数说明：! Ur:初始观测面信号的实部! Ui：初始观测面信号的虚部!******************************************************************** *****subroutine expand_cos_2D(m1,m2,m,n1,n2,n,Ur,Ui)implicit noneinteger m,n,m1,m2,n1,n2real Ur(1:m,1:n),Ui(1:m,1:n)real,allocatable::u(:),r(:)integer j,i,kallocate(u(1:m))do j=n1,n2,1do i=1,m,1u(i)=Ur(i,j)enddocall expand_cos_1d(1,m1,m2,m,u(1)) do i=1,m,1Ur(i,j)=u(i)enddoenddodeallocate(u)allocate(r(1:n))do i=1,m,1do j=1,n,1r(j)=Ur(i,j)enddocall expand_cos_1d(1,n1,n2,n,r(1)) do j=1,n,1Ur(i,j)=r(j)enddoenddodeallocate(r)do i=1,mdo j=1,nUi(i,j)=0enddoenddoend subroutine expand_cos_2D!******************************************************************** ******! 子程序：expand_cos_1d! 功能：一维扩边子程序! 输入参数说明：! n0,n3：扩边后数据起点位置和终点位置! n1,n2：实际数据起点位置和终点位置! feild(i),(i=n1,n1+1,...,n2)：实际数据! 输出参数说明：! field(i),(i=n0,...,n1-1)：扩边后左边的数据! field(i),(i=n2+1,...,n3)：扩边后右边的数据!******************************************************************** ******Subroutine expand_cos_1d(n0,n1,n2,n3,Field)Real Field(n0:n3)pi=3.141592654Field (n0)=(Field (n1)+Field (n2))/2.0Field (n3)=Field (n0)i1=n0i2=n1DO i=i1,i2-1,1Field(i)=Field(i1)+cos(pi/2.0*(i2-i)/(i2-i1))*(Field(i2)-Field(i1)) End doi1=n2i2=n3DO i=i1+1,i2,1Field(i)=Field(i1)+cos(pi/2.0*(i2-i)/(i2-i1))*(Field(i2)-Field(i1)) End doEnd subroutine expand_cos_1d!******************************************************************** ****! 功能：FFT2! 功能：复数组2-D快速Fourier变换! 输入参数说明：! m0,m3：x方向的起点和终点! n0,n3：y方向的起点和终点! field：输入信号(需要赋值给Freal，实部)! m,n：x,y方向扩边后数据终点点号位置（起始点号为1）! NF: 正、反变换标志量(1:反变换;2:正变换)! 输出参数说明：! Freal：信号的实部! Fimage：信号的虚部(对于实信号而言，赋值为0)! 对应频率分布说明：! 数据Freal(m,n)和Fimage(m,n)对应的频率分布位置为:! m 方向：0,1,.......,m/2-1,m/2,-(m/2-1),......,-1! n 方向：0,1,.......,n/2-1,n/2,-(n/2-1),......,-1!******************************************************************** ****SUBROUTINE FFT2(Freal,Fimage,m,n,nf)implicit noneINTEGER m,n,nfREAL Freal(1:m,1:n),Fimage(1:m,1:n)real,ALLOCATABLE::Treal(:),Timage(:)integer nmmax,ierr,i,jnmmax=max(m,n)allocate(Treal(1:nmmax),Timage(1:nmmax),STAT=ierr)if(ierr.ne.0) STOPDO i=1,m,1IF (n.ne.1) THENdo j=1,n,1Treal(j)=Freal(i,j)Timage(j)=Fimage(i,j)end docall FFT(Treal,Timage,n,nf)do j=1,n,1Freal(i,j)=Treal(j)Fimage(i,j)=Timage(j)end doEND IFEND DODO j=1,n,1IF(m.ne.1) THENdo i=1,m,1Treal(i)=Freal(i,j)Timage(i)=Fimage(i,j)end docall FFT(Treal,Timage,m,nf)do i=1,m,1Freal(i,j)=Treal(i)Fimage(i,j)=Timage(i)end doEND IFEND DOdeallocate(Treal,Timage,STAT=ierr)END SUBROUTINE FFT2!****************************************************************** ! 子程序：FFT! 功能：复数组1-D快速Fourier变换! 输入参数说明：! Xreal(n): 输入数据实部! Ximage(n): 输入数据虚部! N: 点数（N 必须为2 的幂次方）! NF: 正、反变换标志量(1:反变换;2:正变换)! 输出参数说明：! Xreal(n): 变换后频谱实部! Ximage(n): 变换后频谱虚部! 对应频率分布说明：! 数据Xreal(n)和Ximage(n)对应的频率分布位置为:! 0,1,.......,n/2-1,n/2,-(n/2-1),......,-1!***************************************************************** SUBROUTINE FFT(Xreal,Ximage,n,nf)implicit noneINTEGER n,nfREAL Xreal(1:n),Ximage(1:n)integer nu,n2,nu1,k,k1,k1n2,l,i,ibitrreal f,p,arg,c,s,treal,timagenu=int(log(float(n))/0.693147+0.001)n2=n/2nu1=nu-1f=float((-1)**nf)k=0DO l=1,nu,1DO while (k.lt.n)do i=1,n2,1p=ibitr(k/2**nu1,nu)arg=6.2831853*p*f/float(n)c=cos(arg)s=sin(arg)k1=k+1k1n2=k1+n2treal=Xreal(k1n2)*c+Ximage(k1n2)*stimage=Ximage(k1n2)*c-Xreal(k1n2)*sXreal(k1n2)=Xreal(k1)-trealXimage(k1n2)=Ximage(k1)-timageXreal(k1)=Xreal(k1)+trealXimage(k1)=Ximage(k1)+timagek=k+1end dok=k+n2END DOk=0nu1=nu1-1n2=n2/2END DODO k=1,n,1i=ibitr(k-1,nu)+1if(i.gt.k) thentreal=Xreal(k)timage=Ximage(k)Xreal(k)=Xreal(i)Ximage(k)=Ximage(i)Xreal(i)=trealXimage(i)=timageend ifEND DOIF(nf.ne.1) THENdo i=1,n,1Xreal(i)=Xreal(i)/float(n)Ximage(i)=Ximage(i)/float(n)end doEND IFEND SUBROUTINE FFTFUNCTION IBITR(J,NU)implicit noneinteger ibitr,j,nuinteger j1,itt,i,j2j1=jitt=0do i=1,nu,1j2=j1/2itt=itt*2+(j1-2*j2)ibitr=ittj1=j2end doEND FUNCTION IBITR!******************************************************************** *******! 子程序：cal_dxdy! 功能：计算x,y方向的点距! 输入参数说明：! x_min,x_max：x的最小值和最大值! y_min,y_max：y的最小值和最大值! N_point,N_line：点数(x方向)、线数(y方向)! 输出参数说明：! dx,dy：x,y方向的点距!******************************************************************** *******subroutine cal_dxdy(xmin,xmax,ymin,ymax,N_POINT,N_LINE,dx,dy)implicit nonereal xmin,xmax,ymin,ymaxinteger N_POINT,N_LINEreal dx,dydx=(xmax-xmin)/N_POINTdy=(ymax-ymin)/N_LINEend subroutine cal_dxdy!****************************************************************** ! 子程序：WAVE2D! 功能：计算2D径向圆频率W! 输入参数说明：! dx: x 方向点距! dy: y 方向线距! m: 点数（M 必须为2 的幂次方）! n: 线数（N 必须为2 的幂次方）! 输出参数说明：! W(m,n): 径向圆频率!****************************************************************** SUBROUTINE WAVE2D(m,n,dx,dy,W)implicit noneINTEGER m,nREAL dx,dyREAL W(1:m,1:n),U(1:m),V(1:n)real pi,delx,delyinteger midm,midn,i,j,xx,yypi=3.1415926midm=m/2+1midn=n/2+1delx=pi*2.0/float(m)/dxdely=pi*2.0/float(n)/dydo j=1,n,1yy=jif(yy.gt.midn) yy=yy-nv(j)=dely*(yy-1)end dodo i=1,m,1xx=iif(xx.gt.midm) xx=xx-mu(i)=delx*(xx-1)end dodo j=1,n,1do i=1,m,1w(i,j)=sqrt(u(i)*u(i)+v(j)*v(j))end doend doEND SUBROUTINE WAVE2D!******************************************************************* ! 子程序：calculate_Y! 功能：计算延拓因子Y! 输入参数说明：! distance:延拓的高度（m/z轴向下为正方向）! m,n：x,y方向扩边后数据终点点号位置（起始点号为1）! W：径向圆频率! 输出参数说明：! Y：延拓因子!******************************************************************* subroutine calculate_Y(m,n,W,y,distance)implicit noneinteger m,nreal pi,distancereal y(1:m,1:n),W(1:m,1:n)real i,jpi=3.1415926do i=1,m,1do j=1,n,1Y(I,J)=1*exp(-1*w(i,j)*distance)enddoenddoend subroutine calculate_Y!******************************************************************** *******! 子程序：calculate_FmulY! 功能：将延拓因子与原信号做乘积! 输入参数说明：! m,n：x,y方向扩边后数据终点点号位置（起始点号为1）! Y：延拓因子! Ur：初始信号的实部! Ui：初始信号的虚部! 输出参数说明：! Fr：延拓信号的实部! Fi：延拓信号的虚部!******************************************************************** *******subroutine calculate_FmulY(m,n,Ur,Ui,y,Fr,Fi)implicit noneinteger m,nreal Ur(1:m,1:n),Ui(1:m,1:n),y(1:m,1:n),Fr(1:m,1:n),Fi(1:m,1:n)real i,jdo i=1,mdo j=1,nFr(i,j)=Ur(i,j)*y(i,j)Fi(i,j)=Ui(i,j)*y(i,j)enddoenddoend subroutine calculate_FmulY!******************************************************************** **************************! 子程序：output_grd! 功能：按照grd格式输出延拓后的场值! 输入参数说明：! N_point,N_line：点数(x方向)、线数(y方向)! filename_output：输出文件的文件名! m1,m2：x方向实际数据起点和终点点号位置! 1,m：x方向扩边后数据起点和终点点号位置! n1,n2：y方向实际数据起点和终点点号位置! 1,n3：y方向扩边后数据起点和终点点号位置! x_min,x_max,y_min,y_max：x/y的最小值、最大值! A：输出数组!******************************************************************** **************************subroutineoutput_grd(filename_output,N_point,N_line,m,n,m1,m2,n1,n2,eigval,A,Xmin,Xmax,Y min,Ymax)character*(*)filename_outputinteger N_point,N_lineinteger m,n,m1,m2,n1,n2real Xmin,Xmax,Ymin,Ymaxreal A(1:m,1:n)real amin,amax,i,jDO j=n1,n2,1do i=m1,m2,1if(ABS(A(i,j)).lt.eigval) thenamin=MIN(amin,A(I,j))amax=MAX(amax,A(I,j))endifEND DOOPEN(20,file=filename_output,status='unknown',form='formatted')write (20,'(A)')'DSAA'write (20,*)m2-m1+1,n2-n1+1write (20,*)Xmin,Xmaxwrite (20,*)Ymin,Ymaxwrite(20,*)amin,amaxDo j=n1,n2,1write(20,*) (A(i,j),i=m1,m2)End doClose(20)end subroutine output_grd!******************************************************************** **************************! 子程序：calculate_error! 功能：计算均方误差! 输入参数说明：! distance:延拓的方向! m1,m2：x方向实际数据起点和终点点号位置! 1,m：x方向扩边后数据起点和终点点号位置! n1,n2：y方向实际数据起点和终点点号位置! 1,n3：y方向扩边后数据起点和终点点号位置! Fr,U：待求均方误差的两数组! 输出参数说明：! error:延拓后场值与真实场值的均方误差!******************************************************************** **************************subroutine calculate_error(distance,Fr,U,m1,m2,n1,n2,m,n,error)implicit noneinteger m1,m2,n1,n2,m,nreal Fr(1:m,1:n),U(1:m,1:n)real distancereal i,jerror=0.0do j=n1,n2do i=m1,m2error=error+(U(i,j)-Fr(i,j))**2enddoenddoerror=sqrt(error)if(distance>0.0)thenwrite(*,*)'请输入向上延拓的均方误差' elsewrite(*,*)'请输入向下延拓的均方误差' end ifwrite(*,*)errorend subroutine。