




版权说明:本文档由用户提供并上传,收益归属内容提供方,若内容存在侵权,请进行举报或认领
文档简介
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=
温馨提示
- 1. 本站所有资源如无特殊说明,都需要本地电脑安装OFFICE2007和PDF阅读器。图纸软件为CAD,CAXA,PROE,UG,SolidWorks等.压缩文件请下载最新的WinRAR软件解压。
- 2. 本站的文档不包含任何第三方提供的附件图纸等,如果需要附件,请联系上传者。文件的所有权益归上传用户所有。
- 3. 本站RAR压缩包中若带图纸,网页内容里面会有图纸预览,若没有图纸预览就没有图纸。
- 4. 未经权益所有人同意不得将文件中的内容挪作商业或盈利用途。
- 5. 人人文库网仅提供信息存储空间,仅对用户上传内容的表现方式做保护处理,对用户上传分享的文档内容本身不做任何修改或编辑,并不能对任何下载内容负责。
- 6. 下载文件中如有侵权或不适当内容,请与我们联系,我们立即纠正。
- 7. 本站不保证下载资源的准确性、安全性和完整性, 同时也不承担用户因使用这些下载资源对自己和他人造成任何形式的伤害或损失。
评论
0/150
提交评论