直立六面体正演程序设计实验.doc

地球物理重磁勘探实验 实验名称：直立六面体正演程序设计实验 学院名称：专业名称：姓名： 学生学号：指导老师：目录一，实验内容：2二，基本原理3（1）模型示意及有关变量的意义，并分别给出重力异常和磁力异常的计算公式，化极。3（2）方向描述及磁化强度分量计算4（3）计算点输入坐标格式设计5 (4) 计算点输入坐标格式设计5三，输入/输出数据格式设计6（1）参数6（2）类型6四，总体设计6五，测试结果7(1)命名文件：filename.txt8(2)场源参数文件：source.txt8(3)平面地形数据文件：bxyz53.txt9(4)输出的重力异常文件：gravity.dat9(5)输出的磁力异常文件：magnetic.dat10(6)化极磁力异常文件：huaji.dat10(7)重力异常，磁力异常，化极磁力异常由surfer得到的图形11六，结论及建议13附录：源程序代码（本程序按照自由格式编写，对于粘贴中难免有一些问题）13一，实验内容：对给定的密度体和磁性体计算平面规则网上的重力异常、磁力异常以及化极磁力异常。二，基本原理（1）模型示意及有关变量的意义，并分别给出重力异常和磁力异常的计算公式，化极。（2）方向描述及磁化强度分量计算（3）计算点输入坐标格式设计(4)计算点输入坐标格式设计三，输入/输出数据格式设计（1）参数场源系数设计：density，Jmay，ai_source,ax_source,ay_source,x1_source,x2_source,y1_source,y2_source,z1_source,z2_source场源参数文件：sourcefilename输入平面规则网文件：Inputfilename重力异常输出文件：gravityfilename磁力异常输出文件：magnaticfilename化极磁力异常：RTPfilename地磁的倾角：ai，ax，ay点数线数：mpoint，nline平面规则网数据：xmin，xmax，ymin，ymax，zmin，zmax（2）类型曲面规则网：输出用GRD格式平面规则网：输出也用GRD格式曲面非规则网：输出采用DAT格式，形式如下：平面非规则网：开始四，总体设计 定义变量参数调用定义文件文名称的子程序 I 调用子程序读入所需的数据 调用子程序依次计算重力异常，磁力异常，化极异常 P 得到结果，输出结果 O结束五，测试结果 测试参数如下：地磁场倾角为50，与x轴的夹角为85，与y轴夹角为5。平面坐标范围与曲面坐标范围相同，但其z坐标值为-5.3km。 (1)命名文件：filename.txt （2）场源参数文件：source.txt (3)平面地形数据文件：bxyz53.txt (4)输出的重力异常文件：gravity.dat (5)输出的磁力异常文件：magnetic.dat (6)化极磁力异常文件：huaji.dat (7)重力异常，磁力异常，化极磁力异常由surfer得到的图形重力异常： 磁力异常： 化极磁力异常：六，结论及建议 本次实验了解了在场源背景下球的重磁异常的方法，绘制出了地下异常的图像。加深了对地下异常分布的宏观认识，对于在课本上的一些抽象的理解在结果图像中有了直观的认识，对于异常的分布趋势以及为什么是这样分布的问题作了讨论。加深了对知道的了解。对于在实验原理以及在程序编写中老师给与的帮助表示感谢。附录：源程序代码（本程序按照自由格式编写，对于粘贴中难免有一些问题）!*(主程序)program Forward_Recfanglecharacter*(80) sourcefilename,Inputfilename,gravityfilename,magneticfilename,huajifilenamereal xmin,xmax,ymin,ymax,zmin,zmaxinteger number_source,mpoint,nlinereal,allocatable:density(:),X_source(:,:),Y_source(:,:),Z_source(:,:),Jmag(:),ai_source(:),ax_source(:),ay_source(:),AI_S(:),AX_S(:),AY_S(:),AI_SS(:),AX_SS(:),AY_SS(:),JMAY(:)real,allocatable:Dg(:,:),mg(:,:),HG(:,:)call Read_cmd(Inputfilename,sourcefilename,gravityfilename,magneticfilename,huajifilename)call get_number(sourcefilename,number_source)allocate (density(1:number_source),x_source(1:2,1:number_source),Y_source(1:2,1:number_source),Z_source(1:2,1:number_source),jmag(1:number_source),ai_source(1:number_source),ax_source(1:number_source),ay_source(1:number_source),AI_S(1:number_source),AX_S(1:number_source),AY_S(1:number_source),AI_SS(1:number_source),AX_SS(1:number_source),AY_SS(1:number_source),JMAY(1:number_source)call Read_source(sourcefilename,density,Jmag,ai_source,ax_source,ay_source,number_source,x_source,Y_source,Z_source)call Input_XYZ(Inputfilename,mpoint,nline,xmin,xmax,ymin,ymax,zmin,zmax)allocate (Dg(1:mpoint,1:nline),mg(1:mpoint,1:nline),HG(1:mpoint,1:nline)call gravity_sub(number_source,density,x_source,Y_source,z_source,Dg,mpoint,nline,xmin,xmax,ymin,ymax,zmin,zmax)call output_grd(Dg,gravityfilename,mpoint,nline,xmin,xmax,ymin,ymax)call magnetic(number_source,JMAg,X_SOURCE,Y_SOURCE,Z_SOURCE,AI_S,AX_S,AY_S,xmin,xmax,ymin,ymax,zmin,zmax,mpoint,nline,MG)call output_grd(Mg,magneticfilename,mpoint,nline,xmin,xmax,ymin,ymax)call RTP(AI_SS,AX_SS,AY_SS,number_source)call magnetic(number_source,JMAg,X_SOURCE,Y_SOURCE,Z_SOURCE,AI_SS,AX_SS,AY_SS,xmin,xmax,ymin,ymax,zmin,zmax,mpoint,nline,HG)call output_grd(HG,huajifilename,mpoint,nline,xmin,xmax,ymin,ymax)end program!*（命名文件）SUBROUTINE read_cmd(Inputfilename,sourcefilename,gravityfilename,magneticfilename,huajifilename)CHARACTER*80 Inputfilename,sourcefilename,gravityfilename,magneticfilename,huajifilenameOPEN(11,file=filename.txt)READ(11,*) inputfilenameREAD(11,*) sourcefilenameREAD(11,*) gravityfilenameREAD(11,*) magneticfilenameREAD(11,*) huajifilenameCLOSE(11)END SUBROUTINE read_cmd!*（读取场源个数）subroutine get_number(sourcefilename,number)character*80 sourcefilenameinteger numberOpen(1,file=sourcefilename,status=old)Number=0Do while(.not.eof(1)Read (1,*,end=100,ERR=100)x,y,zNumber=number+1100 End doCLOSE(1)end subroutine get_number!*subroutine Read_source(sourcefilename,density,Jmag,ai_source,ax_source,ay_source,number_source,x_source,Y_source,Z_source)character*(*) sourcefilenameinteger number_sourcereal density(1:number_source),x_source(1:2,1:number_source),Y_source(1:2,1:number_source),Z_source(1:2,1:number_source),jmag(1:number_source),ai_source(1:number_source),ax_source(1:number_source),ay_source(1:number_source)open(21,file=sourcefilename)do i=1,number_sourceread(21,*) density(i),jmag(i),ai_source(i),ax_source(i),ay_source(i),x_source(1,i),x_source(2,i),y_source(1,i),y_source(2,i),z_source(1,i),z_source(2,i)end doclose(21)end subroutine Read_source!*（读取平面地形数据）subroutine input_xyz(inputfilename,mpoint,nline,xmin,xmax,ymin,ymax,zmin,zmax)character*80 inputfilenameinteger mpoint,nlinereal zmin,zmax,xmin,xmax,ymin,ymaxopen(30,file=INPUTFILENAME)read(30,*)read(30,*) mpoint,nlineread(30,*) xmin,xmaxread(30,*) ymin,ymaxread(30,*) zmin,zmaxclose(30)end subroutine input_xyz!*（计算重力异常）subroutine gravity_sub(number_source,density,x_source,y_source,z_source,Dg,mpoint,nline,xmin,xmax,ymin,ymax,zmin,zmax)integer number_source,mpoint,nlinereal xmin,xmax,ax,zmin,zmax,sum,x,y,rreal Dg(1:mpoint,1:nline)real density(1:number_source),x_source(1:2,1:number_source),y_source(1:2,1:number_source),z_source(1:2,1:number_source)z=0.do n=1,nline,1dy=(ymax-ymin)/(nline-1)y=ymin+(n-1)*dy do m=1,mpoint,1dx=(xmax-xmin)/(mpoint-1)x=xmin+(m-1)*dxDg(m,n)=0.0 do l=1,number_source,1sum=0.0 do i=1,2,1 do j=1,2,1 do k=1,2,1 R=sqrt(x_source(i,l)-x)*2+(y_source(j,l)-y)*2+(z_source(k,l)-z)*2) sum=sum+(-1)*(i+j+k)*(X_SOURCE(i,l)-x)*LOG(Y_SOURCE(j,l)-y)+r)+(Y_SOURCE(j,l)-y)*LOG(X_SOURCE(i,l)-x)+r)-(Z_SOURCE(k,l)-z)*ATAN(X_SOURCE(i,l)-x)*(Y_SOURCE(j,l)-y)/(Z_SOURCE(k,l)-z)/r) end doend do end doDg(m,n)=Dg(m,n)+density(l)*sumend do end doend doend subroutine gravity_sub!*（输出结果文件）SUBROUTINE OUTPUT_GRD(Dg,gravity,m,n,xmin,xmax,ymin,ymax)integer m,nreal Dg(1:m,1:n)character*(*) gravityreal amin,amax,xmin,xmax,ymin,ymaxamin=HUGE(amin)amax=-HUGE(amax)DO j=1,n,1 do i=1,m,1 amin=MIN(amin,Dg(I,j) amax=MAX(amax,Dg(I,j) end doEND DOOPEN(20,file=gravity)write (20,(A)DSAAwrite (20,*)m,nwrite (20,*)xmin,xmaxwrite (20,*)ymin,ymaxwrite(20,*)amin,amaxDo j=1,n,1write(20,*) (Dg(i,j),i=1,m)End doClose(20)END subroutine OUTPUT_GRD!*（计算磁力异常）SUBROUTINE magnetic(number_source,JMAY,X_SOURCE,Y_SOURCE,Z_SOURCE,AI_S,AX_S,AY_S,xmin,xmax,ymin,ymax,zmin,zmax,mpoint,nline,MG)INTEGER number_source,mpoint,nlineREAL xmin,xmax,ymin,ymax,zmin,zmax,k1,k2,k3,k4,k5,k6,r,a,b,cREAL JMAY(1:number_source),AI_S(1:number_source),AX_S(1:number_source),AY_S(1:number_source),X_SOURCE(1:2,1:number_source),Y_SOURCE(1:2,1:number_source),Z_SOURCE(1:2,1:number_source)REAL MG(1:mpoint,1:nline)dx=(xmax-xmin)/(mpoint-1)dy=(ymax-ymin)/(nline-1)z=0.DO n=1,nline,1y=ymin+(n-1)*dy DO m=1,mpoint,1 x=xmin+(m-1)*dx MG(m,n)=0.0 DO l=1,number_source,1sum=0. DO i=1,2,1 DO j=1,2,1 DO k=1,2,1 r=SQRT(X_SOURCE(i,l)-x)*2+(Y_SOURCE(j,l)-y)*2+(Z_SOURCE