点元法重力异常正演.docx

1、 PROGRAMPoint_Element!点元法计算重力异常 EXTERNALDg,ElementDg REALX_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 ENDDO100REWIND(1) PRINT*,K !计算各个源点引起的异常并累加之 DOI=1,K READ(1,*)X_Source,Y

2、_Source,Z_Source,Density GraAnomaly=GraAnomaly+&Dg(X,Y,Z,X_Source,Y_Source,Z_Source,Density,ElementDg) ENDDO CLOSE(1) PRINT*,GraAnomaly ENDPROGRAM !计算某源点在计算点引起的重力异常 FUNCTIONDg(X,Y,Z,X_Source,Y_Source,Z_Source,Density,ElementDg) IMPLICITNONE REALX,Y,Z,X_Source(1:2),Y_Source(1:2),Z_Source(1:2),&Densit

3、y,ElementDg REALDg,ELMT INTEGERI,J,K,SIGN Dg=0 ELMT=0 DOI=1,2 DOJ=1,2 DOK=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 !计算积分表达式中的某一分量 FUNCTIONElementDg(X,Y,Z,X_Source,Y_Source,Z_Source,Density) IMPLICITNONE REAL

4、G PARAMETER(G=6.672E11) REALElementDg,X,Y,Z,X_Source,Y_Source,Z_Source,Density,R,T1,T2,T3 R=SQRT(X_SourceX)*2+(Y_SourceY)*2+(Z_SourceZ)*2) T1=(X_SourceX)*LOG(Y_SourceY)+R) T2=(Y_SourceY)*LOG(X_SourceX)+R) T3=(Z_SourceZ)*ATAN(X_SourceX)*(Y_SourceY)/(Z_SourceZ)*R) ElementDg=G*Density*(T1+T2T3) ENDFUNCTION