滑动平均法.docx
《滑动平均法.docx》由会员分享，可在线阅读，更多相关《滑动平均法.docx（17页珍藏版）》请在冰豆网上搜索。
滑动平均法
PROGRAMHDPJ
CHARACTER*（100）cmdfile,infile,regionfile,localfile,filetype
realeigval
INTEGERiteration
WRITE（*,*）'INPUTCMDFILE-NAME:
'
READ（*,*）cmdfile
CALLREADCMD（cmdfile,filetype,infile,regionfile,localfile,mline,
*nline,radx,rady,iteration,eps,eigval）
IF（（filetype.eq.'grd'）.or.（filetype.eq.'GRD'））THEN
CALLPROCESSGRD（INFILE,REGIONFILE,LOCALFILE,MLINE,NLINE,RADX,
*RADY,ITERATION,EPS,EIGVAL）
ELSEIF（（filetype.eq.'bln'）.or.（filetype.eq.'BLN'））THEN
CALLPROCESSBLN（INFILE,REGIONFILE,LOCALFILE,MLINE,RADX,
*ITERATION,EPS,EIGVAL）
ELSEIF（（filetype.eq.'xyz'）.or.（filetype.eq.'XYZ'））THEN
CALLPROCESSXYZ（INFILE,REGIONFILE,LOCALFILE,MLINE,RADX,
*RADY,ITERATION,EPS,EIGVAL）
ELSE
WRITE（*,*）'无此类型'
ENDIF
ENDPROGRAM
!
参数输入子程序
SUBROUTINEREADCMD（cmdfile,filetype,infile,regionfile,localfile,
*mline,nline,radx,rady,iteration,eps,eigval）
CHARACTER*（*）cmdfile,filetype,infile,regionfile,localfile
OPEN（10,file=cmdfile,status='old',
*access='sequential',form='formatted'）
READ（10,*）filetype
READ（10,*）infile
READ（10,*）regionfile
READ（10,*）localfile
READ（10,*）radx,rady
READ（10,*）iteration
READ（10,*）eps
READ（10,*）eigval
CLOSE（10）
IF（（filetype.eq.'grd'）.or.（filetype.eq.'GRD'））THEN
CALLMNGRD（INFILE,MLINE,NLINE）
ELSEIF（（filetype.eq.'bln'）.or.（filetype.eq.'BLN'））THEN
CALLMNBLN（INFILE,MLINE）
ELSEIF（（filetype.eq.'xyz'）.or.（filetype.eq.'XYZ'））THEN
CALLMNXYZ（INFILE,MLINE）
ELSE
ENDIF
ENDSUBROUTINE
!
*******************************************GRD**************************************************
!
个数输入子程序
SUBROUTINEMNGRD（INFILE,mline,nline）
CHARACTER*（*）INFILE
OPEN（10,file=infile,status='old',access='sequential',
*form='formatted'）
read（10,*）
read（10,*）mline,nline
close（10）
ENDSUBROUTINE
!
GRD类型数据处理子程序
SUBROUTINEPROCESSGRD（INFILE,REGIONFILE,LOCALFILE,MLINE,NLINE,
*RADX,RADY,ITERATION,EPS,EIGVAL）
character*（*）infile,regionfile,localfile
REAL,ALLOCATABLE:
:
first（:
:
）,region（:
:
）,local（:
:
）,
*refirst（:
:
）
ALLOCATE（first（mline,nline）,region（mline,nline）,
*local（mline,nline）,refirst（mline,nline））
callreadgrd（infile,mline,nline,xmin,xmax,ymin,
*ymax,first）
dx=（xmax-xmin）/（mline-1）
dy=（ymax-ymin）/（nline-1）
mr=int（radx/dx+0.5）
nr=int（rady/dy+0.5）
callITERATION_GRD（mline,nline,MR,NR,ITERATION,EPS,EIGVAL,
*first,refirst,region,local）
callwritegrd（regionfile,mline,nline,xmin,xmax,
*ymin,ymax,eigval,region）
callwritegrd（localfile,mline,nline,xmin,xmax,
*ymin,ymax,eigval,local）
DEALLOCATE（first,region,local,refirst）
ENDSUBROUTINE
!
读入grd文件
subroutinereadgrd（infile,mline,nline,xmin,xmax,ymin,
*ymax,first）
character*（*）infile
REALfirst（mline,nline）
open（10,file=infile,status='old',access='sequential',
*form='formatted'）
read（10,*）
read（10,*）
read（10,*）xmin,xmax
READ（10,*）ymin,ymax
READ（10,*）
read（10,*）（（（first（i,j））,I=1,Mline）,J=1,Nline）
close（10）
endsubroutine
!
输出grd文件
subroutinewritegrd（outfile,mline,nline,xmin,xmax,
*ymin,ymax,eigval,outgrd）
character*（*）outfile
REALoutgrd（mline,nline）
zmin=huge（zmin）
zmax=-huge（zmax）
DOj=1,nline,1
DOi=1,mline,1
IF（outgrd（i,j）.lt.0.5*eigval）THEN
zmin=min（zmin,outgrd（i,j））
zmax=max（zmax,outgrd（I,j））
ENDIF
ENDDO
ENDDO
open（11,file=outfile,status='unknown',
*access='sequential',form='formatted'）
write（11,'（a）'）'DSAA'
WRITE（11,*）MLINE,nline
write（11,*）xmin,xmax
write（11,*）ymin,ymax
write（11,*）Zmin,zmax
DOj=1,nline
write（11,*）（outgrd（i,j）,i=1,mline）
ENDDO
close（11）
endsubroutine
!
滑动平均法（grd）
SUBROUTINEITERATION_GRD（mline,nline,MR,NR,ITERATION,EPS,EIGVAL,
*first,refirst,region,local）
realfirst（mline,nline）,region（mline,nline）,
*local（mline,nline）,refirst（mline,nline）
EPSG=0.;ITE=1;REFIRST=FIRST
DOWHILE（EPSG>EPS.OR.ITE<=ITERATION）
EPSG=0.
ITE=ITE+1
DOI=1,MLINE,1
MRMIN=MAX（I-MR,1）
MRMAX=MIN（I+MR,MLINE）
DOJ=1,NLINE,1
IF（REFIRST（I,J）<0.5*EIGVAL）THEN
NRMIN=MAX（J-NR,1）
NRMAX=MIN（J+NR,NLINE）
SUM=0.
NUM=0
DOK=MRMIN,MRMAX,1
DOL=NRMIN,NRMAX,1
IF（REFIRST（K,L）<0.5*EIGVAL）THEN
SUM=SUM+REFIRST（K,L）
NUM=NUM+1
ENDIF
ENDDO
ENDDO
REGION（I,J）=SUM/NUM
EPSG=MAX（EPSG,ABS（REFIRST（I,J）-REGION（I,J）））
ELSE
REGION（I,J）=REFIRST（I,J）
ENDIF
ENDDO
ENDDO
REFIRST=REGION
ENDDO
DOJ=1,NLINE,1
DOI=1,MLINE,1
IF（FIRST（I,J）<0.5*EIGVAL）THEN
LOCAL（I,J）=FIRST（I,J）-REGION（I,J）
ELSE
LOCAL（I,J）=FIRST（I,J）
ENDIF
ENDDO
ENDDO
ENDSUBROUTINE
!
*******************************************BLN**************************************************
!
BLN个数输入程序
SUBROUTINEMNBLN（INFILE,MLINE）
CHARACTER*（*）INFILE
OPEN（10,file=infile,status='old',access='sequential',
*form='formatted'）
read（10,*）mline
close（10）
ENDSUBROUTINE
!
bln文件处理子程序
SUBROUTINEPROCESSBLN（INFILE,REGIONFILE,LOCALFILE,MLINE,RADX,
*ITERATION,EPS,EIGVAL）
character*（*）infile,regionfile,localfile
REAL,ALLOCATABLE:
:
first（:
）,region（:
）,local（:
）,
*refirst（:
）,X（:
）
ALLOCATE（first（mline）,region（mline）,
*local（mline）,refirst（mline）,X（mline））
callreadbln（infile,mline,X,first）
CALLORDER_BLN（X,FIRST,MLINE）
callITERATION_BLN（mline,RADX,ITERATION,EPS,EIGVAL,X,
*first,refirst,region,local）
callwriteBLN（infile,mline,X,first）
callwriteBLN（regionfile,mline,X,region）
callwriteBLN（localfile,mline,x,local）
DEALLOCATE（X,first,region,local,refirst）
ENDSUBROUTINE
!
读取BLN文件
SUBROUTINEreadbln（infile,mline,X,first）
character*（*）infile
REALX（mline）,first（mline）
open（10,file=infile,status='old',access='sequential',
*form='formatted'）
read（10,*）
read（10,*）（（X（I）,first（i））,I=1,Mline）
close（10）
endsubroutine
!
对BLN文件排序
SUBROUTINEORDER_BLN（X,FIRST,MLINE）
REALX（mline）,first（mline）
INTEGERFLAG
FLAG=1
DOWHILE（FLAG==1）
FLAG=0
DOI=1,MLINE-1,1
IF（X（I+1）FLAG=1
TEMP=X（I）
X（I）=X（I+1）
X（I+1）=TEMP
TEMP=FIRST（I）
FIRST（I）=FIRST（I+1）
FIRST（I+1）=TEMP
ENDIF
ENDDO
ENDDO
endsubroutine
!
滑动平均法（BLN）
SUBROUTINEITERATION_BLN（mline,RADX,ITERATION,EPS,EIGVAL,X,
*first,refirst,region,local）
realX（MLINE）,first（mline）,region（mline）,
*local（mline）,refirst（mline）
EPSG=0.;ITE=1;REFIRST=FIRST
DOWHILE（EPSG>EPS.OR.ITE<=ITERATION）
EPSG=0.
ITE=ITE+1
DOI=1,MLINE,1
IF（REFIRST（I）<0.5*EIGVAL）THEN
J=I;SUM=0;NUM=0
DOWHILE（（X（J）+RADX）>=X（I）.AND.J>=1）
IF（REFIRST（J）<0.5*EIGVAL）THEN
SUM=SUM+REFIRST（J）
NUM=NUM+1
ENDIF
J=J-1
ENDDO
J=I+1
DOWHILE（（X（I）+RADX）>=X（J）.AND.J<=MLINE）
IF（REFIRST（J）<0.5*EIGVAL）THEN
SUM=SUM+REFIRST（J）
NUM=NUM+1
ENDIF
J=J+1
ENDDO
REGION（I）=SUM/NUM
EPSG=MAX（EPSG,ABS（REFIRST（I）-REGION（I）））
ELSE
REGION（I）=REFIRST（I）
ENDIF
ENDDO
REFIRST=REGION
ENDDO
DOI=1,MLINE,1
IF（FIRST（I）<0.5*EIGVAL）THEN
LOCAL（I）=FIRST（I）-REGION（I）
ELSE
LOCAL（I）=FIRST（I）
ENDIF
ENDDO
ENDSUBROUTINE
!
输出BLN文件
subroutinewriteBLN（outfile,mline,x,outBLN）
character*（*）outfile
REALoutBLN（mline）,X（MLINE）
open（11,file=outfile,status='unknown',
*access='sequential',form='formatted'）
WRITE（11,*）MLINE
DOI=1,MLINE,1
WRITE（11,*）X（I）,OUTBLN（I）
ENDDO
close（11）
endsubroutine
!
*******************************************xyz**************************************************
!
xyz个数检测程序
SUBROUTINEMNXYZ（filename,number）
character*（*）filename
logicalunit_open,file_exist
integernunit_in
nunit_in=11
unit_open=.TRUE.
dowhile（unit_open）
nunit_in=nunit_in+1
INQUIRE（unit=nunit_in,opened=unit_open）
enddo
INQUIRE（file=filename,exist=file_exist）
if（file_exist）then
open（unit=nunit_in,file=filename,status='old'）
number=0
DOwhile（.not.EOF（nunit_in））
read（nunit_in,*,end=200,err=200）x,f
number=number+1
200ENDDO
close（unit=nunit_in）
else
write（*,'（A,A,A）'）'文件:
',TRIM（filename）,'不存在!
'
endif
ENDsubroutine
!
xyz文件处理子程序
SUBROUTINEPROCESSXYZ（INFILE,REGIONFILE,LOCALFILE,MLINE,
*RADX,RADY,ITERATION,EPS,EIGVAL）
character*（*）infile,regionfile,localfile
REAL,ALLOCATABLE:
:
first（:
）,region（:
）,local（:
）,
*refirst（:
）,XY（:
:
）
ALLOCATE（first（mline）,region（mline）,
*local（mline）,refirst（mline）,XY（2,mline））
callreadxyz（infile,mline,XY,first）
CALLORDER_XYZ（XY,FIRST,MLINE）
callITERATION_XYZ（mline,RADX,RADY,ITERATION,EPS,EIGVAL,XY,
*first,refirst,region,local）
callwriteXYZ（regionfile,mline,XY,region）
callwriteXYZ（localfile,mline,XY,local）
DEALLOCATE（XY,first,region,local,refirst）
ENDSUBROUTINE
!
读取XYZ文件
SUBROUTINEreadxyz（infile,mline,XY,first）
character*（*）infile
REALfirst（mline）,XY（2,MLINE）
open（10,file=infile,status='old',access='sequential',
*form='formatted'）
DOI=1,MLINE,1
READ（10,*）XY（1,I）,XY（2,I）,FIRST（I）
ENDDO
close（10）
ENDSUBROUTINE
!
对XYZ文件排序
SUBROUTINEORDER_XYZ（XY,FIRST,MLINE）
REALXY（2,mline）,first（mline）
INTEGERFLAG
FLAG=1
DOWHILE（FLAG==1）
FLAG=0
DOI=1,MLINE-1,1
IF（XY（1,I+1）FLAG=1
TEMP=XY（1,I）
XY（1,I）=XY（1,I+1）
XY（1,I+1）=TEMP
TEMP=XY（2,I）
XY（2,I）=XY（2,I+1）
XY（2,I+1）=TEMP
TEMP=FIRST（I）
FIRST（I）=FIRST（I+1）
FIRST（I+1）=TEMP
ELSEIF（XY（1,I+1）==XY（1,I）.AND.XY（2,I+1）FLAG=1
TEMP=XY（2,I）
XY（2,I）=XY（2,I+1）
XY（2,I+1）=TEMP
TEMP=FIRST（I）
FIRST（I）=FIRST（I+1）
FIRST（I+1）=TEMP
ENDIF
ENDDO
ENDDO
endsubroutine
!
滑动平均法（XYZ）
SUBROUTINEITERATION_XYZ（mline,RADX,RADY,ITERATION,EPS,EIGVAL,XY,
*first,refirst,region,local）
realXY（2,MLINE）,first（mline）,region（mline）,
*local（mline）,refirst（mline）
EPSG=0.;ITE=1;REFIRST=FIRST
DOWHILE（EPSG>EPS.OR.ITE<=ITERATION）
EPSG=0.
升级会员 