点元法重力异常正演.docx

点元法重力异常正演.docx

《点元法重力异常正演.docx》由会员分享，可在线阅读，更多相关《点元法重力异常正演.docx（2页珍藏版）》请在冰豆网上搜索。

 PROGRAM Point_Element !

点元法计算重力异常

 EXTERNAL Dg,ElementDg

 REAL X_Source（1:

2）,Y_Source（1:

2）,Z_Source（1:

2）,Density,Dg

 REAL :

:

 GraAnomaly=0.,X=0.,Y=0.,Z=0.

 INTEGER :

:

 K=0,I

 OPEN（1,FILE='SOURCE.DAT',ACTION='READ'）

 !

统计源点个数

 DO 

  READ（1,*,END=100）

  K=K+1

 ENDDO

100 REWIND

（1）

 PRINT*,K

 !

计算各个源点引起的异常并累加之  

 DO I=1,K

 READ（1,*） X_Source,Y_Source,Z_Source,Density

 GraAnomaly=GraAnomaly+

     &Dg（X,Y,Z,X_Source,Y_Source,Z_Source,Density,ElementDg）

 ENDDO

 CLOSE

（1）

 PRINT*,GraAnomaly

 ENDPROGRAM

 !

计算某源点在计算点引起的重力异常

 FUNCTION Dg（X,Y,Z,X_Source,Y_Source,Z_Source,Density,ElementDg）

 IMPLICIT NONE

 REAL X,Y,Z,X_Source（1:

2）,Y_Source（1:

2）,Z_Source（1:

2）,

     &Density,ElementDg

 REAL Dg,ELMT

 INTEGER I,J,K,SIGN

 Dg=0

 ELMT=0

 DO I=1,2

  DO J=1,2

   DO K=1,2

    SIGN=（‐1）**（I+J+K）

    ELMT=ElementDg（X,Y,Z,X_Source（I）,

     &Y_Source（J）,Z_Source（K）,Density）

    Dg=Dg+ELMT*SIGN

   ENDDO

  ENDDO

 ENDDO

 ENDFUNCTION

 !

计算积分表达式中的某一分量

 FUNCTION ElementDg（X,Y,Z,X_Source,Y_Source,Z_Source,Density）

 IMPLICIT NONE

 REAL G

 PARAMETER（G=6.672E‐11）

 REAL ElementDg,X,Y,Z,X_Source,Y_Source,Z_Source,Density,R,T1,T2,T3

 R=SQRT（（X_Source‐X）**2+（Y_Source‐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