C @(#)spesearch.for 17.1.1.1 (ESO-DMD) 01/25/02 17:56:15 C=========================================================================== C Copyright (C) 1995 European Southern Observatory (ESO) C C This program is free software; you can redistribute it and/or C modify it under the terms of the GNU General Public License as C published by the Free Software Foundation; either version 2 of C the License, or (at your option) any later version. C C This program is distributed in the hope that it will be useful, C but WITHOUT ANY WARRANTY; without even the implied warranty of C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the C GNU General Public License for more details. C C You should have received a copy of the GNU General Public C License along with this program; if not, write to the Free C Software Foundation, Inc., 675 Massachusetts Ave, Cambridge, C MA 02139, USA. C C Correspondence concerning ESO-MIDAS should be addressed as follows: C Internet e-mail: midas@eso.org C Postal address: European Southern Observatory C Data Management Division C Karl-Schwarzschild-Strasse 2 C D 85748 Garching bei Muenchen C GERMANY C=========================================================================== C C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C.COPYRIGHT: Copyright (c) 1987 European Southern Observatory, C all rights reserved C C.VERSION: 1.0 ESO-FORTRAN Conversion, AA 19:26 - 3 DEC 1987 C C.LANGUAGE: F77+ESOext C C.AUTHOR: D PONZ C C.IDENTIFICATION: C PROGRAM SPESEARCH C C.KEYWORD: C SEARCH, LINES C C.PURPOSE: C SEARCH FOR EMISSION/ABSORPTION LINES IN A 1D- OR 2D-IMAGE AND C CALCULATE PARAMETERS : C - LINE CENTER C - PEAK VALUE C C.ALGORITHM: C C LINES ARE DETECTED BY USING A MEDIAN + THRESHOLDING ALGORITHM C ABSORPTION LINES ARE TREATED AS NEGATIVE EMISSION C CENTER OF THE LINE IS COMPUTED AS C (1)- CENTER OF GRAVITY OF THE 2 HIGHEST PIXELS RELATIVE TO THE C THIRD ONE C (2)- CENTER OF THE GAUSSIAN FIT TO THE POINTS ABOVE THE BACKGROUND C (3)- PIXEL CORRESPONDING TO THE MAXIMUM C C.INPUT/OUTPUT: C THE FOLLOWING KEYWORDS SHOULD BE DEFINED : C P1 : IMAGE TO BE SEARCHED C INPUTR : SEARCH PARAMETERS C WINDOW, THRESHOLD, NAVER C P3 : TABLE FILE C P4 : SEARCH METHOD C GRAVITY - DEFAULT C GAUSSIAN C MAXIMUM C P5 : EMISSION/ABSORPTION C C MODIFICATIONS : C M Peron , DEc90, enables the comand SEARCH/RBR C------------------------------------------------------------------ C PROGRAM SPESEA IMPLICIT NONE INTEGER MD, NMAX PARAMETER (MD=2) PARAMETER (NMAX=1000) INTEGER MADRID INTEGER I,IAV,IC,IMETH,IMODE,IS,ISTAT INTEGER NAX,NOBJ,NAVER,ISIZE,KUN,KNUL,TID,IMNO INTEGER NPIX(MD) INTEGER*8 IPM,IPNTA C CHARACTER*2 METH,IEA,ACTION CHARACTER*80 INA CHARACTER*80 FOUTA CHARACTER*48 UNIT CHARACTER*72 INST CHARACTER*80 LINE,LINEM,LINEE C REAL RBUF(2,NMAX) REAL STRA(MD),STPA(MD),PARM(5) DOUBLE PRECISION DSTRA(2), DSTPA(2) C CHARACTER*16 LABEL(5),TUNIT(5) CHARACTER*6 FORM(5) C INCLUDE 'MID_INCLUDE:ST_DEF.INC' COMMON /VMR/MADRID(1) INCLUDE 'MID_INCLUDE:ST_DAT.INC' C DATA LABEL(1)/'X'/,TUNIT(1)/'World Coord'/,FORM(1)/'F10.2'/ DATA LABEL(2)/'Y'/,TUNIT(2)/'PIXEL'/,FORM(2)/'F10.2'/ DATA LABEL(3)/'PEAK'/,TUNIT(3)/'DN'/,FORM(3)/'E12.3'/ C C ... INITIATE MIDAS C CALL STSPRO('SPESEA') CALL STKRDC('P1',1,1,80,IAV,INA,KUN,KNUL,IS) CALL STIGET(INA,D_R4_FORMAT,F_I_MODE,F_IMA_TYPE, . MD,NAX,NPIX,DSTRA,DSTPA,INST,UNIT,IPNTA,IMNO,IS) NPIX(1) = MAX(1,NPIX(1)) NPIX(2) = MAX(1,NPIX(2)) NAX = 2 STRA(1) = DSTRA(1) STRA(2) = DSTRA(2) STPA(1) = DSTPA(1) STPA(2) = DSTPA(2) IF (NPIX(1).EQ.1) THEN NPIX(1) = NPIX(2) NPIX(2) = 1 END IF C C ... READ PARAMETERS C C NMAX = 800 CALL STKRDR('INPUTR',1,3,IAV,PARM,KUN,KNUL,IS) CALL STKRDC('P4',1,1,2,IAV,METH,KUN,KNUL,IS) CALL STKRDC('P5',1,1,1,IAV,IEA,KUN,KNUL,IS) CALL STKRDC('ACTION',1,1,2,IAV,ACTION,KUN,KNUL,IS) CALL FORUPC(METH,METH) CALL FORUPC(IEA,IEA) IF (METH(1:2).EQ.'GA') THEN IMETH = -1 LINEM = ' centering method : gaussian fit' ELSE IF (METH(1:2).EQ.'GR') THEN IMETH = 0 LINEM = ' centering method : center of gravity' ELSE IMETH = 1 LINEM = ' centering method : position of max(min)' END IF END IF IF (IEA(1:1).EQ.'E') THEN IMODE = 1 LINEE = ' search for emission lines' ELSE IMODE = 0 LINEE = ' search for absorption lines' END IF C C ... CREATE THE TABLE FILE PLUS COLUMNS C CALL STKRDC('P3',1,1,80,IAV,FOUTA,KUN,KNUL,IS) CALL TBTINI(FOUTA,F_TRANS,F_O_MODE,15,NMAX,TID,IS) CALL STDWRD(TID,'PIXEL',DSTPA,1,1,KUN,IS) DO 10 I = 1,3 CALL TBCINI(TID,D_R4_FORMAT,1,FORM(I), . TUNIT(I),LABEL(I),IC,IS) 10 CONTINUE CALL STTPUT(' search lines ',ISTAT) CALL STTPUT(' ------------',ISTAT) CALL STTPUT(' input image : '//INA,ISTAT) CALL STTPUT(' output table : '//FOUTA,ISTAT) CALL STTPUT(' input parameters ',ISTAT) WRITE (LINE,9000) NINT(PARM(1)) CALL STTPUT(LINE,ISTAT) WRITE (LINE,9010) PARM(2) CALL STTPUT(LINE,ISTAT) NAVER = PARM(3) IF (NAVER.GT.0) THEN WRITE (LINE,9030) NAVER CALL STTPUT(LINE,ISTAT) ENDIF CALL STTPUT(LINEM,ISTAT) CALL STTPUT(LINEE,ISTAT) C C ... DO THE SEARCH C IF (ACTION.EQ.'RB') THEN ISIZE = NPIX(1)*4 CALL TDMGET(ISIZE,IPM,ISTAT) CALL SRCOD3(MADRID(IPNTA),NAX,NPIX(1),NPIX(2),STRA,STPA, + NMAX,PARM,NOBJ,TID,RBUF,IMETH,IMODE,NAVER,MADRID(IPM)) CALL TDMFRE(ISIZE,IPM,ISTAT) ELSE IF (NAVER.LE.1) THEN CALL SRCOD1(MADRID(IPNTA),NAX,NPIX(1),NPIX(2),STRA,STPA, + NMAX,PARM,NOBJ,TID,RBUF,IMETH,IMODE) ELSE ISIZE = NPIX(1)*4 CALL TDMGET(ISIZE,IPM,ISTAT) CALL SRCOD2(MADRID(IPNTA),NAX,NPIX(1),NPIX(2),STRA,STPA, + NMAX,PARM,NOBJ,TID,RBUF,IMETH,IMODE,NAVER,MADRID(IPM)) CALL TDMFRE(ISIZE,IPM,ISTAT) ENDIF ENDIF WRITE (LINE,9020) NOBJ CALL STTPUT(LINE,ISTAT) C C ... FINISHED C CALL TBSINI(TID,IS) CALL TBTCLO(TID,IS) CALL STSEPI 9000 FORMAT (' search window : ',I3,' pixels') 9010 FORMAT (' detection threshold : ',F10.2,' DN') 9030 FORMAT (' average on : ',I3,' scan-lines') 9020 FORMAT (' no. of detections: ',I8) END SUBROUTINE SRCOD1(F,NAF,NPF1,NPF2,STRA,STPA,NMAX,PARM,NOBJ,TID, + RBUF,IMETH,IMODE) C C C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C.IDENTIFICATION C SUBROUTINE SRCOD1 VERSION 1 23 NOV 1982 C J.D.PONZ C C.KEYWORDS: C SEARCH C C.PURPOSE: C SEARCH SPECTRA FOR EMISSION LINES C C.ALGORITHM: C USE FINDM2 C C.INPUT/OUTPUT C F original image Input C NAF dimension of F Input C NPF number of pixels in F (x,y,z...) Input C STRA starting point in F Input C STPA step Input C PARM floating-point parameters of filter Input C NOBJ number of objects found Output C OBJ charateristics of objects Output C C------------------------------------------------------------- C C IMPLICIT NONE INTEGER NPF1, NPF2, NAF, NMAX, NOBJ, IMETH, IMODE REAL STRA(2),STPA(2),PARM(2),OBJ(6) REAL RBUF(2,NMAX),F(NPF1,NPF2) REAL XAUX(100) INTEGER IAUX(100), I, NSEQ, IS, NACT, TID INTEGER ICOL(5), NAUX, NT, NW, ISTAT, IAC, NW1, IORD REAL T CHARACTER LINE*80 C DATA ICOL/1,2,3,4,5/ C NAUX = 99 C C initialize variables C NT = NPF1 NW = PARM(1) T = PARM(2) NOBJ = 0 C C search features C IF (NPF2.GT.1) THEN CALL STTPUT(' seq.no. detected lines ',ISTAT) CALL STTPUT(' ------- -------------- ',ISTAT) END IF IAC = 3 NW1 = 2*NW + 1 NSEQ = 0 DO 20 IORD = 1, NPF2 CALL FNDMDN(F(1,IORD),NT,NW1,T,NMAX,RBUF,NACT,XAUX,IAUX, + NAUX, IMETH,IMODE,STRA(1),STPA(1)) IF (NACT.LE.0) CALL FNDMAX(F(1,IORD),NT,RBUF,NMAX,NACT, + STRA(1),STPA(1)) IF (NPF2.GT.1) THEN WRITE (LINE,9000) IORD,NACT CALL STTPUT(LINE,ISTAT) END IF NOBJ = NOBJ + NACT DO 10 I = 1,NACT OBJ(1) = RBUF(1,I) OBJ(2) = (IORD-1)*STPA(2) + STRA(2) OBJ(3) = RBUF(2,I) NSEQ = NSEQ + 1 CALL TBRWRR(TID,NSEQ,IAC,ICOL,OBJ,IS) 10 CONTINUE 20 CONTINUE IF (NPF2.GT.1) CALL STTPUT(' ----------------------- ',ISTAT) RETURN C 9000 FORMAT (1X,I8,2X,I8) END SUBROUTINE FNDMAX(X,N,R,NMAX,NACT,START,STEP) C C FIND MAXIMUM VALUE IN CASE OF HIGH THRESHOLD C IMPLICIT NONE INTEGER N, NMAX, NACT REAL X(N),R(2,NMAX) REAL START, STEP C INTEGER I, IP REAL XMAX C XMAX = 1.E-30 DO 10 I = 3,N - 2 IF (X(I).GT.XMAX) THEN XMAX = X(I) IP = I END IF 10 CONTINUE R(1,1) = START + (IP-1)*STEP R(2,1) = XMAX NACT = 1 RETURN END SUBROUTINE SRCOD3(F,NAF,NPF1,NPF2,STRA,STPA,NMAX,PARM,NOBJ, + TID,RBUF,IMETH,IMODE,NAVER,RVAL) C C C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C.IDENTIFICATION C SUBROUTINE SRCOD3 VERSION 1 23 NOV 1982 C J.D.PONZ C C.KEYWORDS: C SEARCH C C.PURPOSE: C SEARCH SPECTRA FOR EMISSION LINES C C.ALGORITHM: C USE FINDM2 C C.INPUT/OUTPUT C F original image Input C NAF dimension of F Input C NPF number of pixels in F (x,y,z...) Input C STRA starting point in F Input C STPA step Input C PARM floating-point parameters of filter Input C NOBJ number of objects found Output C OBJ charateristics of objects Output C C------------------------------------------------------------- C C IMPLICIT NONE INTEGER NPF1, NPF2, NAF, NMAX, NOBJ, IMETH, IMODE, NAVER REAL STRA(2),STPA(2),PARM(2),OBJ(6) REAL RBUF(2,NMAX), F(NPF1,NPF2), RVAL(NPF1) REAL XAUX(100) INTEGER IAUX(100), I, NSEQ, IS, NACT, TID INTEGER J, J1, J2, K INTEGER ICOL(5), NAUX, NT, NW, ISTAT, IAC, NW1, IORD REAL T CHARACTER LINE*80 C DATA ICOL/1,2,3,4,5/ C NAUX = 99 C NT = NPF1 NW = PARM(1) T = PARM(2) NOBJ = 0 C C search features C IF (NPF2.GT.1) THEN CALL STTPUT(' seq.no. detected lines ',ISTAT) CALL STTPUT(' ------- -------------- ',ISTAT) END IF IAC = 3 NW1 = 2*NW + 1 NSEQ = 0 DO 20 IORD = 1,NPF2 DO 1 J = 1, NPF1 RVAL(J) = 0. 1 CONTINUE J1 = MAX(IORD-NAVER,1) J2 = MIN(J1 + NAVER,NPF2) DO 3 J = J1, J2 DO 2 K = 1, NPF1 RVAL(K) = RVAL(K) + F(K,J)/(NAVER*2+1) 2 CONTINUE 3 CONTINUE CALL FNDMDN(RVAL,NT,NW1,T,NMAX,RBUF,NACT,XAUX,IAUX, + NAUX, IMETH,IMODE,STRA(1),STPA(1)) IF (NACT.LE.0) CALL FNDMAX(RVAL,NT,RBUF,NMAX,NACT, + STRA(1),STPA(1)) IF (NPF2.GT.1) THEN WRITE (LINE,9000) IORD,NACT CALL STTPUT(LINE,ISTAT) END IF NOBJ = NOBJ + NACT DO 10 I = 1,NACT OBJ(1) = RBUF(1,I) OBJ(2) = (IORD-1)*STPA(2) + STRA(2) OBJ(3) = RBUF(2,I) NSEQ = NSEQ + 1 CALL TBRWRR(TID,NSEQ,IAC,ICOL,OBJ,IS) 10 CONTINUE 20 CONTINUE IF (NPF2.GT.1) CALL STTPUT(' ----------------------- ',ISTAT) RETURN C 9000 FORMAT (1X,I8,2X,I8) END SUBROUTINE SRCOD2(F,NAF,NPF1,NPF2,STRA,STPA,NMAX,PARM,NOBJ, + TID,RBUF,IMETH,IMODE,NAVER,RVAL) C C C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C.IDENTIFICATION C SUBROUTINE SRCOD2 VERSION 1 23 NOV 1982 C J.D.PONZ C C.KEYWORDS: C SEARCH C C.PURPOSE: C SEARCH SPECTRA FOR EMISSION LINES C C.ALGORITHM: C USE FINDM2 C C.INPUT/OUTPUT C F original image Input C NAF dimension of F Input C NPF number of pixels in F (x,y,z...) Input C STRA starting point in F Input C STPA step Input C PARM floating-point parameters of filter Input C NOBJ number of objects found Output C OBJ charateristics of objects Output C C------------------------------------------------------------- C C IMPLICIT NONE INTEGER NPF1, NPF2, NAF, NMAX, NOBJ, IMETH, IMODE, NAVER REAL STRA(2),STPA(2),PARM(2),OBJ(6) REAL RBUF(2,NMAX), F(NPF1,NPF2), RVAL(NPF1) REAL XAUX(100) INTEGER IAUX(100), I, NSEQ, IS, NACT, TID INTEGER I1, J, J1, J2, K INTEGER ICOL(5), NAUX, NT, NW, ISTAT, IAC, NW1, IORD REAL T CHARACTER LINE*80 C DATA ICOL/1,2,3,4,5/ C NAUX = 99 C C initialize variables C NT = NPF1 NW = PARM(1) T = PARM(2) NOBJ = 0 C C search features C IF (NPF2.GT.1) THEN CALL STTPUT(' seq.no. detected lines ',ISTAT) CALL STTPUT(' ------- -------------- ',ISTAT) END IF IAC = 3 NW1 = 2*NW + 1 NSEQ = 0 I1 = NAVER/2+1 DO 20 IORD = I1, NPF2, NAVER DO 1 J = 1, NPF1 RVAL(J) = 0. 1 CONTINUE J1 = IORD - NAVER/2 J2 = J1 + NAVER - 1 DO 3 J = J1, J2 DO 2 K = 1, NPF1 RVAL(K) = RVAL(K) + F(K,J)/NAVER 2 CONTINUE 3 CONTINUE CALL FNDMDN(RVAL,NT,NW1,T,NMAX,RBUF,NACT,XAUX,IAUX, + NAUX, IMETH,IMODE,STRA(1),STPA(1)) IF (NACT.LE.0) CALL FNDMAX(RVAL,NT,RBUF,NMAX,NACT, + STRA(1),STPA(1)) IF (NPF2.GT.1) THEN WRITE (LINE,9000) IORD,NACT CALL STTPUT(LINE,ISTAT) END IF NOBJ = NOBJ + NACT DO 10 I = 1,NACT OBJ(1) = RBUF(1,I) OBJ(2) = (IORD-1)*STPA(2) + STRA(2) OBJ(3) = RBUF(2,I) NSEQ = NSEQ + 1 CALL TBRWRR(TID,NSEQ,IAC,ICOL,OBJ,IS) 10 CONTINUE 20 CONTINUE IF (NPF2.GT.1) CALL STTPUT(' ----------------------- ',ISTAT) RETURN C 9000 FORMAT (1X,I8,2X,I8) END SUBROUTINE FNDMDN(IA,LEN,WIDTH,THRES,NMAX,RBUF,NACT,FA,LINK,NAUX, + IMETH,IMODE,XSTR,XSTP) C+ C.COPYRIGHT: Copyright (c) 1987 European Southern Observatory, C all rights reserved C C.VERSION: 1.0 ESO-FORTRAN Conversion, AA 23:58 - 3 DEC 1987 C C.LANGUAGE: F77+ESOext C C.AUTHOR: D.PONZ C C.IDENTIFICATION C C SUBROUTINE FNDMDN C C.KEYWORDS C C SEARCH C C.PURPOSE C C SEARCH EMISSION/ABSORPTION FEATURES C C.ALGORITHM C C DATA IN A RUNNING WINDOW IS SORTED TO FORM A LOCAL HISTOGRAM C AN EMISSION(ABSORPTION) FEATURE IS FOUND WHEN THE CENTRAL C VALUE IN THE WINDOW IS GREATER THAN THE MEDIAN VALUE PLUS THE C THRESHOLD. C C.INPUT/OUTPUT C C IA REAL*4 INPUT ARRAY C LEN INTG*4 NUMBER OF INPUT PIXELS C WIDTH INTG*4 WINDOW WIDTH C THRES REAL*4 THRESHOLD C NMAX INTG*4 MAXIMUM NO. OF DETECTED LINES C RBUF REAL*4 RBUF(1,I) - CENTER C RBUF(2,I) - MAXIMUM VALUE C NACT INTG*4 ACTUAL NUMBER OF DETECTED LINES C FA(NAUX) REAL*4 AUXILIARY SPACE C LINK(NAUX) INTG*4 AUXILIARY SPACE C NAUX INTG*4 SPACE SIZE. NAUX .GT. WIDTH+1 C IMETH INTG*4 CENTERING METHOD AS C IMODE INTG*4 0 C XSTR REAL*4 STARTING COORD. C XSTP REAL*4 STEP C C- IMPLICIT NONE INTEGER NAUX,LINK(NAUX),MIN,MAX,WIDTH,LEN,MIDPT,FRONT,AVAIL INTEGER RANGE,NEWNOD,I,J,K,NODE,ITMP INTEGER NMAX,NACT,IMODE,IMETH INTEGER ISTART,IFAIL,NPIX,NOBJ,NLOBJ REAL FA(NAUX),IA(LEN),RBUF(2,NMAX),AM REAL THRES,T,XSTR,XSTP,CHICHK,XGO,XCNTR REAL ACOE(4),WORK1(200),WORK2(200),XPREV REAL DELTA LOGICAL IFIRST C C SET UP BEGINNING (MIN) AND END (MAX) OF FILTER ARRAY C ALSO SET UP POINTER TO AVAIL LIST AND LIST HEADER (FA(1),LINK(1)) C CHICHK = 0.005 IFIRST = .TRUE. T = ABS(THRES) MIN = 2 MAX = NAUX AVAIL = MIN FA(1) = -1.0 LINK(1) = 0 NOBJ = 0 C C SET UP FIRST FILTER WINDOW C RANGE = WIDTH - 1 DO 10 I = 1,RANGE,1 NODE = NEWNOD(IA(I),AVAIL,MIN,MAX,FA,LINK,NAUX) CALL INSNOD(NODE,FA,LINK,NAUX) 10 CONTINUE RANGE = WIDTH/2 C C FOR EACH POINT FROM WIDTH/2 TO LEN-WIDTH/2+1 C FIND ITS MEDIAN C FRONT = MIN C FRONT IS THE OLDEST POINT IN THE FILTER ARRAY K = RANGE C K IS THE POINT BEING PROCESSED MIDPT = RANGE + 1 C MIDPT IS THE MIDDLE OF THE WINDOW (IN FA) DO 70 J = WIDTH,LEN,1 NODE = NEWNOD(IA(J),AVAIL,MIN,MAX,FA,LINK,NAUX) CALL INSNOD(NODE,FA,LINK,NAUX) K = K + 1 C C FIND THE MIDPT - - IT IS THE MEDIAN VALUE C ITMP = LINK(1) DO 20 I = 2,MIDPT,1 ITMP = LINK(ITMP) 20 CONTINUE C C GET THE LINK WHICH POINTS TO THE VALUE OF THE MEDIAN C IF (ABS(FA(ITMP)-IA(K)).GT.T) THEN IF (IFIRST) THEN IFIRST = .FALSE. NLOBJ = 1 ISTART = K - 1 NPIX = 3 ELSE NPIX = NPIX + 1 END IF ELSE IF ( .NOT. IFIRST) THEN IFIRST = .TRUE. C C FIND LINE CENTER C XGO = XSTR + (ISTART-1)*XSTP IF (IMETH) 30,40,50 C C GAUSSIAN C 30 CALL SGAUS(IA(ISTART),WORK1,WORK2,IMODE,NPIX,IFAIL, + XGO,XSTP,XCNTR,CHICHK,1,AM,ACOE) GO TO 60 C C GRAVITY C 40 CALL GRAVT(IA(ISTART),NPIX,IMODE,IFAIL,XGO,XSTP,XCNTR, + AM) IF (IFAIL.EQ.0) GO TO 60 C C MAXIMUM C 50 CALL CNTRH(IA(ISTART),NPIX,IMODE,IFAIL,XGO,XSTP,XCNTR, + AM) 60 IF (IFAIL.EQ.0) THEN NOBJ = NOBJ + 1 RBUF(1,NOBJ) = XCNTR RBUF(2,NOBJ) = AM IF (NOBJ.EQ.NMAX) GO TO 80 END IF END IF END IF CALL DELNOD(FRONT,FA,LINK,NAUX) FRONT = FRONT + 1 IF (FRONT.GT.MAX) FRONT = MIN C C DELETE OLDEST AND ADVANCE PTR, FRONT ALSO WRAPS C 70 CONTINUE 80 CONTINUE C C CHECK DOUBLETS C XPREV = RBUF(1,1) I = 1 85 CONTINUE IF ( .NOT. (I.LT.NOBJ)) GO TO 110 DELTA = 2.*XSTP IF (ABS(XPREV-RBUF(1,I+1)).LE.DELTA) THEN DO 90 J = I,NOBJ - 2 RBUF(1,J) = RBUF(1,J+2) RBUF(2,J) = RBUF(2,J+2) 90 CONTINUE NOBJ = NOBJ - 2 ELSE I = I + 1 END IF XPREV = RBUF(1,I) GOTO 85 110 CONTINUE NACT = NOBJ RETURN END