重磁实验二Word文件下载.docx

1、半空间狄利克莱问题解析解：其中：称为延拓因子，为计算面Z坐标，Z轴向下为正方向，为计算面频率域位场，为延拓面的频率域位场。2 输入/输出数据格式设计2.1 输入/输出数据文件名输入数据和输出数据文件名均保存在“parameter.txt”中。第一行为输入的低高度观测面数据文件；第二行为输出的高高度观测面数据文件；第三行第四行依次为输入的扩边比例因子和延拓高度。A20_mag.grdA53_mag.grd1.53.32.2重要变量名filename_Field:低高度观测面数据文件filename_Conti:高高度观测面数据文件Field（m,n）: 低高度观测面数据Conti（m,n）: 高

2、高度观测面数据error: 延拓后的均方误差factor_x: 扩边比例因子（1.0）height: 延拓高度（0：向上延拓，向下延拓）Factor_Conti: 延拓因子point: 点数line: 线数m: 扩边以后总点数（满足2的幂次方）n: 扩边以后总线数（满足2的幂次方）3 总体设计I输入有关参数：filename_Field,filename_Conti,factor_x,height从文件输入有关参数：point,line计算有关参数：m,m1,m2,n,n1,n2从文件输入有关数据：xmin,xmax,ymin,ymax,dx,dy,FieldP进行扩边处理傅里叶正变换计算延拓

3、因子进行乘积运算傅里叶反变换O输出计算结果：Conti,error4 测试结果5 结论及建议（1） 由高高度观测面位场等值线图（图2）可以看出，向上延拓的效果较好，能清晰地反映磁异常位置，等值线光滑，均方误差仅为1.660418。（2） 与低高度观测（图1）相比，向上延拓3.3m的高高度观测（图2）得到的磁异常幅度有所降低，位场等值线图明显对1、2号这两个局部异常体进行了压制，突出了3号区域异常体。（3） 本次程序设计过程中的一些认识：扩边参数m,n计算完成后再开辟数组，否则会导致在扩边过程中数组空间不够用；二维扩边可利用一维扩边子程序，将大大减轻二维的繁琐程度；对老师提供的FFT和计算圆频率

4、的子程序要理解输入/输出参数的含义并灵活运用；理解延拓因子的计算公式中各个参量的含义才能进行程序实现；最大的感受就是在一开始虽然对整体的程序流程有了大致的把握，但落实到具体的程序编写时却发现了对各种细节理解的偏颇之处，正所谓纸上谈兵就是如此吧，只有在实践中才能真正地理解和提升。源程序代码!*c 功能：频率域位场延拓 文件参数说明： filename_Field: filename_Conti: 输入参数说明： factor_x: height: point: line: xmin,xmax: 点坐标的最大值，最小值 ymin,ymax: 线坐标的最大值，最小值 dx: x 方向点距 dy: y

5、 方向线距 扩边参数说明： m1,m2: 实际数据起点位置和终点位置（m2-m1+1=point） m: n1,n2: 实际数据起点位置和终点位置（n2-n1+1=line） n: 输入数据说明： Field（m,n）: Freal（m,n）: 低高度观测面数据数据实部 Fimage（m,n）: 低高度观测面数据数据虚部 延拓参数说明： U（m）: x 方向对应的圆频率 V（n）: y 方向对应的圆频率 W（m,n）: 径向圆频率 Factor_Conti: 输出数据说明： Conti（m,n）: Creal（m,n）: 高高度观测面数据数据实部 Cimage（m,n）: 高高度观测面数据数据

6、虚部 error:-cProgram PotenConti implicit none character*（80） filename_Field,filename_Conti real factor_x,height integer point,line real xmin,xmax,ymin,ymax,dx,dy integer m,m1,m2,n,n1,n2 real error real,allocatable:Field（:,:）,Freal（:）,Fimage（:）Conti（:）,Creal（:）,Cimage（:U（:）,V（:）,W（:Factor_Conti（: !读取文件

7、名及参数 call Get_parameter（filename_Field,filename_Conti,factor_x,height）读取点线数 call Get_point_line（filename_Field,point,line）计算扩边数据 call GET_M1M2m（point,factor_x,m,m1,m2） call GET_M1M2m（line,factor_x,n,n1,n2）开辟数组 allocate（Field（1:m,1:n）,Freal（1:n）,Fimage（1:n） allocate（Conti（1:n）,Creal（1:n）,Cimage（1: al

8、locate（U（1:m）,V（1:n）,W（1: allocate（Factor_Conti（1:读取低高度观测面数据 call Input_Field（filename_Field,m,m1,m2,n,n1,n2,xmin,xmax,ymin,ymax,dx,dy,Field）二维扩边 call Expend_cos_2D（m,m1,m2,n,n1,n2,Field）2D FFT正变换 Freal=Field Fimage=0.0 call FFT2（Freal,Fimage,m,n,2） 计算二维圆频率 call WAVE2（m,n,dx,dy,U,V,W） call Calculate

9、_Factor_Conti（m,n,W,height,Factor_Conti）乘以延拓因子 call Product_Factor_Conti（m,n,Freal,Fimage,Factor_Conti,Creal,Cimage）2D FFT反变换 call FFT2（Creal,Cimage,m,n,1）输出高高度观测面数据 Conti=Creal call Output_Conti（filename_Conti,point,line,m,m1,m2,n,n1,n2,xmin,xmax,ymin,ymax,Conti）计算延拓后均方误差 call Calculate_error（m,m1,

10、m2,n,n1,n2,Field,Conti,error）End Program PotenConti读取文件名及参数子程序Subroutine Get_parameter（filename_Field,filename_Conti,factor_x,height） character*（*） filename_Field,filename_Conti open（10,file=parameter.txt,status=old read（10,*） filename_Field read（10,*） filename_Conti read（10,*） factor_x read（10,*） h

11、eight close（10）End Subroutine Get_parameter读取点线数子程序Subroutine Get_point_line（filename_Field,point,line） character*（*） filename_Field open（20,file=filename_Field,status= read（20,*） read（20,*） point,line close（20）End Subroutine Get_point_line扩边数据位置确定子程序 mpoint: 实际数据个数 输出参数说明： 实际数据起点位置和终点位置（m2-m1+1=mpoint）SUBROUTINE GET_M1M2m（mpoint,factor_x,m,m1,m2） INTEGER mpoint,m,m1,m2 REAL factor_x integer mtemp,mu 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 mu=int（log（float（mpoint）/0.693147+factor_x） m=2*mu end if m