重磁异常空间域正演

1、!对给定的密度体和磁性体计算起伏观测面以及平面上的重力异常、磁力异常以及化极 磁力异常。PROGRAM Complex_body_forwardIMPLICIT NONECHARACTER*80 cmdfile,sourcefile,stationfile,output_grafile,&output_magfile,output_polfileINTEGER component,order,plane场的分量标识，导数的阶数，观测面的形 态标识INTEGER n_source,n_station,i,j,mpoint,line,cs,cxyzREAL xmin,xmax,ymin,ymaxR

2、EAL,ALLOCATABLE:SOURCE(:,:),XYZ(:,:)REAL,ALLOCATABLE:ANOM_GRA(:),ANOM_MAG(:),ANOM_POL(:)cmdfile=cmd.datcs=11cxyz=3CALL INPUT_CMDFILE(component,order,plane,cmdfile,sourcefile,&stationfile,output_grafile,output_magfile,output_polfile)CALL CHECK_CMD(component,order,plane,stationfile)CALL GET_NUMBER_SOU

3、RCE(sourcefile,n_source)ALLOCATE(SOURCE(1:n_source,1:cs)CALL INPUT_SOURCE(sourcefile,SOURCE,n_source,cs)CALL GET_n_station(stationfile,plane,xmin,xmax,ymin,ymax,&mpoint,line,n_station)ALLOCATE(XYZ(1:n_station,1:cxyz)CALL INPUT_STATION(stationfile,plane,xmin,xmax,ymin,ymax,&mpoint,line,XYZ,n_station)

4、ALLOCATE(ANOM_GRA(1:n_station),ANOM_MAG(1:n_station)ALLOCATE(ANOM_POL(1:n_station)CALL CAL_ANOM(SOURCE,n_source,cs,XYZ,n_station,cxyz,&ANOM_GRA,ANOM_MAG,ANOM_POL,component)CALL OUTPUT(plane,XYZ,ANOM_GRA,ANOM_MAG,ANOM_POL,n_station,xmin&,xmax,ymin,ymax,mpoint,line,output_grafile,output_magfile&,outpu

5、t_polfile)ENDPROGRAM!读取场的分量标识、导数的阶数、观测面的形态标识!读取场源参数文件名、测点文件名、输出文件名SUBROUTINE INPUT_CMDFILE(component,order,plane,cmdfile,sourcefile,&stationfile,output_grafile,output_magfile,output_polfile)IMPLICIT NONEINTEGER component,order,planeCHARACTER*80 cmdfile,sourcefile,stationfile,output_grafile,&output_m

6、agfile,output_polfileOPEN(12,FILE=cmdfile,ACTION=READ)READ(12,*) output_polfileREAD(12,*) output_polfileREAD(12,*) output_polfileREAD(12,*) componentREAD(12,*) output_polfileREAD(12,*) orderREAD(12,*) output_polfileREAD(12,*) output_polfileREAD(12,*) output_polfileREAD(12,*) output_polfileREAD(12,*)

7、 output_polfileREAD(12,*) planeREAD(12,*) output_polfileREAD(12,*) sourcefileREAD(12,*) output_polfileREAD(12,*) stationfileREAD(12,*) output_polfileREAD(12,*) output_grafileREAD(12,*) output_polfileREAD(12,*) output_magfileREAD(12,*) output_polfileCLOSE(12)ENDSUBROUTINE!检查输入参数是否正确SUBROUTINE CHECK_C

8、MD(component,order,plane,stationfile)IMPLICIT NONEEXTERNAL UPPERINTEGER component,order,plane,ccCHARACTER*80 stationfile,UPPERCHARACTER signPRINT*,请确认以下信息是否正确：SELECT CASE(component)CASE(1)PRINT*,1计算重力或磁力异常CASE(2)PRINT*,1计算重力或磁力异常的导数CASE DEFAULTPRINT*,错误：场分量类型不合法，请修改参数文件！STOPENDSELECTzpnN&d(寸)冗0zpnN&

9、d(mmsvoEniis件Bzk zpnN&dQmsvo zpnN&d (Imsvo (兽2_dmsoID 山2SLL.IaNLU dons 二efin扇而苣妄蓄-啪鲤E&d NLUKL (。voqnz H Euouod eou)ll.iLL.K1N 山 dons j f esn、现WK 胡熟e仅亚、K-pnN&d n 山 kl (clxrdu . (8专8)挡 JU.Q1US&山 dd n)义 .aNV.(寸 Hu2_d&0.ZHu2_d)LILL.K1N 山 donsj f esn、现WK 胡熟e仅亚、K-pnN&d N 山工_1(0&0-:兽.(8毕8_巨0柘切)&山 ddn)义 .aN

10、.(mHHou2_d&oIHHou2_d)LL.I(挡 JU0Qss)WRLLN%H8ID 山_1山 SC1N 山dons 二efin扇而K胡熟匡K-啪鲤MN&dnnvLL.捋山 s=97.AND.ASC = 122) THENASC=ASC-32string(i:i)=CHAR(ASC)ENDIFENDDOUPPER=stringRETURN!读取场源的个数SUBROUTINE GET_NUMBER_SOURCE(filename,n)IMPLICIT NONEINTEGER nCHARACTER*80 filenameCHARACTER*20 aaOPEN(101,file=filenam

11、e,action=read)n=0DOREAD(101,*,END=100,ERR=100) aan=n+1ENDDO100 CLOSE(101)ENDSUBROUTINE!读取场源参数SUBROUTINE INPUT_SOURCE(sourcefile,SOURCE,m,n)IMPLICIT NONEINTEGER m,n,i,jREAL ax,ay,inclineCHARACTER*80 sourcefileREAL SOURCE(1:m,1:n)OPEN(12,FILE=sourcefile,ACTION=READ)READ(12,*)(SOURCE(i,j),j=1,n),i=1,m)

12、CLOSE(12)DO i=1,max=SOURCE(i,4);ay=SOURCE(i,5);incline=SOURCE(i,3)IF(ax+ay)=90.OR.(ax+ay)=270.OR.ABS(ax-ay)=90.)THENELSEPRINT*,错误：磁化方向参数不正确！请修改场源参数文件！PRINT*,出错行：,iSTOPENDIFIF(incline90)THENPRINT*,错误：磁倾角参数不正确！请修改场源参数文件！PRINT*,出错行：,iSTOPENDIFENDDOENDSUBROUTINE!读取测量点的个数SUBROUTINE GET_n_station(stationf

13、ile,plane,xmin,xmax,ymin,ymax,&mpoint,line,n_station)IMPLICIT NONEINTEGER plane,mpoint,line,n_stationREAL xmin,xmax,ymin,ymaxCHARACTER*80 stationfileIF(plane= = 1.OR.plane=3) THENCALL GET_GRID(stationfile,mpoint,line,xmin,xmax,ymin,ymax)n_station=mpoint*lineELSEmpoint=0;line=0;xmin=0;xmax=0;ymin=0;y

14、max=0CALL GET_NUMBER_SOURCE(stationfile,n_station)ENDIF!读取规则网测网网格参数SUBROUTINE GET_GRID(inputfile,mpoint,line,xmin,xmax,ymin,ymax)IMPLICIT NONECHARACTER*80 inputfileINTEGER mpoint,lineREAL xmin,xmax,ymin,ymaxOPEN(101,FILE=inputfile,ACTION=READ)READ(101,*)READ(101,*)mpoint,line,xmin,xmax,ymin,ymaxCLOS

15、E(101)ENDSUBROUTINE!读取测点坐标SUBROUTINE INPUT_STATION(stationfile,plane,xmin,xmax,ymin,ymax,&mpoint,line,XYZ,n_station)IMPLICIT NONEREAL xmin,xmax,ymin,ymax,dx,dy,tempINTEGER plane,mpoint,line,n_stationREAL XYZ(1:n_station,1:3)CHARACTER*80 stationfileREAL i,jOPEN(12,FILE=stationfile,ACTION=READ)IF(plan

16、e= = 1.OR.plane=3) THENDO i=1,5,1READ(12,*)ENDDOdx=(xmax-xmin)/(mpoint-1)dy=(ymax-ymin)/(line-1)DO j=1,lineDO i=1,mpointXYZ(j-1)*mpoint+i,1)=(i-1)*dx+xminXYZ(j-1)*mpoint+i,2)=(j-1)*dy+yminENDDOENDDOIF(plane= = 1) THENREAD(12,*) XYZ(1,3)XYZ(:,3)=XYZ(1,3)ELSEREAD(12,*) XYZ(:,3)ENDIFELSEIF(plane=2.OR.p

17、lane=4) THENIF(plane=2) THENREAD(12,*) XYZ(1,:)XYZ(:,3)=XYZ(1,3)DO i=2,n_stationREAD(12,*) XYZ(i,1),XYZ(i,2)ENDDOELSEDO i=1,n_stationREAD(12,*) XYZ(i,:)ENDDOENDIFELSEPRINT*,测网类型不合法！ENDIF!计算各测点的异常值SUBROUTINE CAL_ANOM(SOURCE,n_source,cs,XYZ,n_station,cxyz,&ANOM_GRA,ANOM_MAG,ANOM_POL,component)IMPLICIT

18、 NONEEXTERNAL Dg,ElementDg,DTREAL Dg,ElementDg,dg0,DT,dt0INTEGER n_source,cs,n_station,cxyz,i,j,componentREAL SOURCE(1:n_source,1:cs),XYZ(1:n_station,1:cxyz)REAL ANOM_GRA(1:n_station),ANOM_MAG(1:n_station)REAL ANOM_POL(1:n_station)ANOM_GRA=0.;ANOM_MAG=0.;ANOM_POL=0.IF(component= = 1) THENDO i=1,n_st

19、ation,1DO j=1,n_sourcedg0=Dg(XYZ(i,1)*1E6,XYZ(i,2)*1E6,XYZ(i,3)*1E6,SOURCE(j,6:7)*1E6&,SOURCE(j,8:9)*1E6,SOURCE(j,10:11)*1E6,SOURCE(j,1),&ElementDg)ANOM_GRA(i)=ANOM_GRA(i)+dg0dt0=DT(XYZ(i,1)*1E6,XYZ(i,2)*1E6,XYZ(i,3)*1E6,SOURCE(j,6:7)*1E6&,SOURCE(j,8:9)*1E6,SOURCE(j,10:11)*1E6,SOURCE(j,2)*1E-6,&SOUR

20、CE(j,3),SOURCE(j,4),SOURCE(j,5)ANOM_MAG(i)=ANOM_MAG(i)+dt0dt0=DT(XYZ(i,1)*1E6,XYZ(i,2)*1E6,XYZ(i,3)*1E6,SOURCE(j,6:7)*1E6&,SOURCE(j,8:9)*1E6,SOURCE(j,10:11)*1E6,SOURCE(j,2)*1E-6,&90.,SOURCE(j,4),SOURCE(j,5)ANOM_POL(i)=ANOM_POL(i)+dt0ENDDOENDDOELSEIF(component=2) THENPRINT*,重磁异常导数计算程序正在开发中.PRINT*,$目前

21、暂不支持此项功育性STOPENDIFENDSUBROUTINE!*计算某源点在计算点引起的重力异常程序集开始*FUNCTION Dg(X,Y,Z,X_Source,Y_Source,Z_Source,Density,ElementDg)IMPLICIT NONEREAL X,Y,Z,X_Source(1:2),Y_Source(1:2),Z_Source(1:2),&Density,ElementDgREAL Dg,ELMTINTEGER I,J,K,SIGNDg=0ELMT=0DO I=1,2DO J=1,2DO K=1,2SIGN=(-1)*(I+J+K)ELMT=ElementDg(X,

22、Y,Z,X_Source(I),&Y_Source(J),Z_Source(K),Density)Dg=Dg+ELMT*SIGNENDDOENDDOENDDOENDFUNCTION!计算积分表达式中的某一分量FUNCTION ElementDg(X,Y,Z,X_Source,Y_Source,Z_Source,Density)IMPLICIT NONEREAL GPARAMETER(G=6.672E-8)REAL ElementDg,X,Y,Z,X_Source,Y_Source,Z_Source,Density,R,T1,T2,T3R=SQRT(X_Source-X)*2+(Y_Source

23、-Y)*2+(Z_Source-Z)*2)T1=(X_Source-X)*LOG(Y_Source-Y)+R)T2=(Y_Source-Y)*LOG(X_Source-X)+R)T3=(Z_Source-Z)*ATAN(X_Source-X)*(Y_Source-Y)/(Z_Source-Z)*R)ElementDg=-G*Density*(T1+T2-T3)ENDFUNCTION!*计算某源点在计算点引起的重力异常程序集结束*!*计算某源点在计算点引起的重力异常程序集开始*FUNCTION DT(x,y,z,X_Source,Y_Source,Z_Source,moment,incline,

24、ax,ay)IMPLICIT NONEREAL x,y,z,moment,incline,ax,ay,DTREAL X_Source(1:2),Y_Source(1:2),Z_Source(1:2)INTEGER i,j,k,signREAL hax,hay,za,mx,my,mz,pi,r,xs,ys,zspi=3x=moment*COSD(incline)*COSD(ax)my=moment*COSD(incline)*COSD(ay)mz=moment*SIND(incline)hax=0.;hay=0.;za=0.DO i=1,2,1DO j=1,2,1DO

25、k=1,2,1sign=(-1)*(i+j+k)xs=X_Source(i);ys=Y_Source(j);zs=Z_Source(k)r=SQRT(x-xs)*2+(y-ys)*2+(z-zs)*2)hax=hax+sign*(-mx*ATAN(xs-x)*(ys-y)/(xs-x)*2+r*(zs-z)+& (zs-z)*2)+my*LOG(r+zs-z)+mz*LOG(r+ys-y)hay=hay+sign*(mx*LOG(r+zs-z)-my*ATAN(xs-x)*(ys-y)/(&(ys-y)*2+r*(zs-z)+(zs-z)*2)+mz*LOG(r+xs-x)za=za+sign

26、*(mx*LOG(r+ys-y)+my*LOG(r+xs-x)-mz*&ATAN(xs-x)*(ys-y)/r/(zs-z)ENDDOENDDOENDDODT=hax*COSD(incline)*COSD(ax)+hay*COSD(incline)*COSD(ay)+& za*SIND(incline)DT=DT/4/piENDFUNCTION!*计算某源点在计算点引起的重力异常程序集结束*!输出重磁异常到外部文件中SUBROUTINEOUTPUT(plane,XYZ,ANOM_GRA,ANOM_MAG,ANOM_POL,n_station&,xmin,xmax,ymin,ymax,mpoin

27、t,line,output_grafile,output_magfile&,output_polfile)INTEGER plane,n_station,mpoint,line,i,jREAL ANOM_GRA(1:n_station),ANOM_MAG(1:n_station)REAL ANOM_POL(1:n_station),XYZ(1:n_station,1:3)REAL xmin,xmax,ymin,ymax,gzmin,gzmax,mzmin,mzmax,pzmin,pzmaxREAL eigvCHARACTER*80 output_grafile,output_magfile,output_polfileeigv=1.7014100091878E+038OPEN(12,FILE=output_grafile,ACTION=WRITE)OPEN(13,FILE=output_magfile,ACTION=WRITE)OPEN(14,FILE=output_polfile,ACTION=