直立六面體正演程序設(shè)計實驗

1、地球物理重磁勘探實驗 實驗名稱：直立六面體正演程序設(shè)計實驗 學(xué)院名稱：專業(yè)名稱：姓名： 學(xué)生學(xué)號：指導(dǎo)老師：一、基本原理（1）模型示意及有關(guān)變量的意義，并分別給出重力異常和磁力異常的計算公式，化極。（2）方向描述及磁化強度分量計算（3）計算點輸入坐標格式設(shè)計(4)計算點輸入坐標格式設(shè)計三，輸入/輸出數(shù)據(jù)格式設(shè)計（1）參數(shù)場源系數(shù)設(shè)計：density，Jmay，ai_source,ax_source,ay_source,x1_source,x2_source,y1_source,y2_source,z1_source,z2_source場源參數(shù)文件：sourcefilename輸入平面規(guī)則網(wǎng)文件

2、：Inputfilename重力異常輸出文件：gravityfilename磁力異常輸出文件：magnaticfilename化極磁力異常：RTPfilename地磁的傾角：ai，ax，ay點數(shù)線數(shù)：mpoint，nline平面規(guī)則網(wǎng)數(shù)據(jù)：xmin，xmax，ymin，ymax，zmin，zmax（2）類型曲面規(guī)則網(wǎng)：輸出用GRD格式平面規(guī)則網(wǎng)：輸出也用GRD格式曲面非規(guī)則網(wǎng)：輸出采用DAT格式，形式如下：平面非規(guī)則網(wǎng)：開始四，總體設(shè)計 定義變量參數(shù)調(diào)用定義文件文名稱的子程序 I 調(diào)用子程序讀入所需的數(shù)據(jù) 調(diào)用子程序依次計算重力異常，磁力異常，化極異常 P 得到結(jié)果，輸出結(jié)果 O結(jié)束重力異常，

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!*（讀取場源個數(shù)）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!*（讀取平面地形數(shù)據(jù)）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!*（輸出結(jié)果文件）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