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

1、地球物理重磁勘探实验 实验名称：直立六面体正演程序设计实验 学院名称：专业名称：姓名： 学生学号：指导老师：一、基本原理（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输入平面规则网文件

2、：Inputfilename重力异常输出文件：gravityfilename磁力异常输出文件：magnaticfilename化极磁力异常：RTPfilename地磁的倾角：ai，ax，ay点数线数：mpoint，nline平面规则网数据：xmin，xmax，ymin，ymax，zmin，zmax（2）类型曲面规则网：输出用GRD格式平面规则网：输出也用GRD格式曲面非规则网：输出采用DAT格式，形式如下：平面非规则网：开始四，总体设计 定义变量参数调用定义文件文名称的子程序 I 调用子程序读入所需的数据 调用子程序依次计算重力异常，磁力异常，化极异常 P 得到结果，输出结果 O结束重力异常，

3、磁力异常，化极磁力异常由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_so

4、urce(:,:),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 (

5、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_s

6、ource),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:mp

7、oint,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,y

8、min,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,huajifilenam

9、e,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,*) inputfilenameRE

10、AD(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

11、(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_sou

12、rce(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

13、(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,ymaxop

14、en(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

15、,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+(

16、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_SOU

17、RCE(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=-HUG

18、E(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)')'DSAA'write (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!*

19、（计算磁力异常）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),

20、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(k,l)-z)*2) a=cosd(AI_S(i)*cosd(AX_S(i) b=cosd(AI_S(i)*sind(AX_S(i) c=sind(A