C @(#)prflcl.for 17.1.1.1 (ES0-DMD) 01/25/02 17:15:44 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 Massachusetss Ave, Cambridge, C MA 02139, USA. C C Corresponding 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.IDENTIFICATION C subroutine PRFLCL version 1.3 830915 C A. Kruszewski ESO Garching C modified version 1.4 870303 C A. Kruszewski Obs. de Geneve C modified for MSDOS version 2.0 880830 C A. Kruszewski Warsaw U. Obs. C.PURPOSE C cleans two-dimensional profile "PRFL" from influence of blended C components C puts informations about newly found components into array "MCM0" C.INPUT/OUTPUT C input arguments C PRFL real*4 array two-dimensional profile C ICNT integer*4 array numbers of pixels used to calculate PRFL C AVER real*4 array one-dimensional profile C LIM integer*4 radial size of profiles in pixels C IDIF integer*4 distance in pixels to the nearest C catalogued object C LOUT integer*4 output controlling flag C RARR real*4 array values of real keywords C output arguments C PRFL real*4 array cleaned two-dimensional profile C MCM0 integer*4 array array holding informations about C newly found objects C KCLN integer*4 size in pixels of good non-cleaned C profile C KSAT integer*4 size in pixels of saturated part C of profile C----------------------------------------------------------------------- SUBROUTINE PRFLCL(PRFL, ICNT, AVPR, IAVR, LIM, & BGRD, IDIF, LDBG, RARR, MCM0, & KCLN, KSAT) C IMPLICIT NONE INCLUDE 'MID_REL_INCL:INVENT.INC/NOLIST' C REAL PRFL(8,0:MAXSUB) INTEGER ICNT(8,0:MAXSUB) REAL AVPR(0:MAXSUB) INTEGER IAVR(0:MAXSUB) INTEGER LIM REAL BGRD INTEGER IDIF INTEGER LDBG REAL RARR(64) INTEGER MCM0(21) INTEGER KCLN INTEGER KSAT C INTEGER ICPY(0:9,0:MAXSUB) INTEGER ISTAT INTEGER INDX(8,0:MAXSUB) C INTEGER IIND(0:9) INTEGER IOUT(8) INTEGER IB, IR, IC, ID, IK, II INTEGER IDLM INTEGER LDIF INTEGER K, KM, KN, KS, KK, K1 INTEGER L, L4 INTEGER LIM2 INTEGER JJ, JJ1 C REAL CPY(8,0:MAXSUB) REAL CLPR1, CLPR2, CLPR3, CLPR4, CLPR5, CLPR6 REAL T1(2) REAL DFPR(MAXSUB) REAL FAMP(3,MAXSUB) REAL TEMP(MAXSUB) REAL COEF(7,MAXSUB) REAL ATRSH, AVER REAL BROBJ, BTRSH REAL PZERO, PMIN, PMAX REAL FMIN1, FMIN2 REAL FCTR REAL CORR REAL CTRSH REAL SCALE REAL PR0 REAL HCUT REAL TMP REAL SUM REAL TRSH REAL TRLIM C CHARACTER*80 TEXT LOGICAL BLEND,NEW,BRIGHT LOGICAL LOG01,LOG02,LOG03,LOG04,LOG05,LOG06,LOG07 LOGICAL LOG08,LOG09,LOG10,LOG11,LOG12,LOG13,LOG14 LOGICAL LOG15,LOG16 C SCALE = 1000.0 / BGRD HCUT = RARR(2) TRSH = RARR(3) IF (PRFL(1,0) .GT. 0.0) THEN TRLIM = TRSH**0.8 * PRFL(1,0)**0.2 ELSE TRLIM = TRSH ENDIF C CLPR1 = RARR(49) CLPR2 = RARR(50) IF (CLPR1 .GT. 0.1) THEN CLPR3 = 1.0 / CLPR1 ELSE CLPR3 = 10.0 ENDIF C IF (CLPR1 .GT. 0.1) THEN CLPR4 = 1.0 / CLPR1 ELSE CLPR4 = 10.0 ENDIF C CLPR6 = RARR(51) IF (CLPR6 .LT. 1.0) THEN CLPR6 = 1.0 ENDIF C CLPR5 = 1.0 / CLPR6 ATRSH = CLPR2 * TRSH BTRSH = -ATRSH ATRSH = ATRSH CTRSH = TRSH / CLPR2 BROBJ = RARR(47) IDLM = INT(RARR(42)) LDIF = IDIF + IDLM C C ****** Mark empty places in the profile. C DO 10 K = 0,MAXSUB DO 20 L = 1,8 IF (ICNT(L,K) .EQ. 0) THEN INDX(L,K) = -1 ELSE INDX(L,K) = 0 ENDIF 20 CONTINUE 10 CONTINUE KM = 0 DO 30 K = 1,21 MCM0(K) = 0 30 CONTINUE C C ****** Find largest saturated ring number KSAT. C c HHCUT = 0.7 * (HCUT-BGRD) c K = -1 C c 40 CONTINUE c K = K + 1 c IF (IAVR(K) .EQ. 0 .OR. AVPR(K) .GT. HHCUT) GOTO 40 c KSAT = K - 1 C C ****** Calculate Fourier amplitudes over octants FAMP. C CALL OCTFRR( PRFL , ICNT , AVPR , LIM , COEF , FAMP ) C C ****** Outputs average profiles AVPR and C ****** Fourier amplitudes FAMP into screen. C IF (LDBG .GE. 3) THEN WRITE (TEXT,'(A)') 'Average profile and first three'// & ' Fourier amplitudes' CALL STTPUT(TEXT,ISTAT) LIM2 = MIN(10,LIM) DO 45 K = 1,LIM2 IOUT(1) = NINT(AVPR(K) * SCALE) DO 46 L = 1,3 IOUT(L+1) = NINT (FAMP(L,K) * SCALE) 46 CONTINUE WRITE(TEXT,'(4I12)') (IOUT(L),L=1,4) CALL STTPUT(TEXT,ISTAT) 45 CONTINUE ENDIF C C ****** Look for faint component in each octant separately. C KS = MIN(KSAT+2,LIM) DO 50 L = 1,8 C C *** Fix starting values of various quantities. C BLEND = .FALSE. BRIGHT = .FALSE. IB = 0 K = KS PZERO = MAX(PRFL(L,0),PRFL(L,K-1),PRFL(L,K)) PMIN = PRFL(L,K+1) PMAX = PRFL(L,K+1) IF (K .EQ. 1) THEN FMIN1 = MIN(FAMP(1,1),FAMP(1,2)) FMIN2 = MIN(FAMP(2,1),FAMP(2,2)) ELSE FMIN1 = FAMP(1,K) FMIN2 = FAMP(2,K) ENDIF DFPR(K) = CLPR1 * (PRFL(L,K) - AVPR(K)) C C *** Start DO WHILE loop. C 61 CONTINUE IF (K .GE. LIM) GOTO 60 K = K + 1 IF (FAMP(1,K) .LT. FMIN1) THEN FMIN1 = FAMP(1,K) ENDIF IF (FAMP(2,K) .LT. FMIN2) THEN FMIN2 = FAMP(2,K) ENDIF IF (PRFL(L,K) .GT. PMAX) THEN PMAX = PRFL(L,K) IF (PMIN .LT. PZERO .AND. PMAX .GT. PZERO) THEN BRIGHT = .TRUE. ELSE BRIGHT = .FALSE. ENDIF ELSE BRIGHT = .FALSE. ENDIF IF (PRFL(L,K) .LT. PMIN) THEN PMIN = PRFL(L,K) ENDIF DFPR(K) = CLPR1 * (PRFL(L,K) - AVPR(K)) C C *** Detect inverted profile. C LOG01 = PRFL(L,K).GT.2.0*PMIN LOG02 = PRFL(L,K).GT.ATRSH LOG03 = DFPR(K).GE.DFPR(K-1) LOG04 = DFPR(K).GE.CLPR3*AVPR(K) LOG05 = PRFL(L,K) .GE. PRFL(L,K-1) LOG06 = DFPR(K).GT.CLPR4*(AVPR(K-1)-AVPR(K)) LOG07 = PRFL(L,K).GT.ATRSH LOG08 = FAMP(1,K).GT.FMIN1 LOG09 = FAMP(1,K).EQ.FMIN1 LOG10 = FAMP(1,K).LT.CTRSH LOG11 = FAMP(2,K).GT.FMIN2 LOG12 = CLPR6.GE.5.0 LOG13 = FAMP(2,K).LT.CTRSH LOG14 = MIN(FAMP(1,K),FAMP(2,K)).LT.CTRSH LOG15 = FAMP(1,K).GT.CLPR5*FAMP(2,K) LOG16 = FAMP(1,K).LT.CLPR6*FAMP(2,K) IF ( (LOG01 .AND. LOG02) .OR. & (((LOG03 .AND. LOG04) .OR. LOG05) .AND. LOG06 .AND. & LOG07 .AND. ((LOG08 .OR. (LOG09 .AND. LOG10)) .AND. & (LOG11 .OR. LOG12 .OR. LOG13) .AND. & (LOG14 .OR. (LOG15 .AND. LOG16))))) THEN C C *** Check if an inverted profile C *** has not been detected before. C IF (.NOT.BLEND .AND. IB.EQ.0) THEN C C *** Component is detected for the first time. C BLEND=.TRUE. C C **** Marks profile points affected C **** by component by means C **** of setting INDX(L,K) to 1. C IF ( PRFL(L,K) .GE. PRFL(L,K-2) ) THEN IF ( K .GT. 2 .AND. INDX(L,K-1) .EQ. 0 ) THEN INDX(L,K-1) = 1 ENDIF IB = 3 ELSE IB = 2 ENDIF IF ( INDX(L,K) .EQ. 0 ) THEN INDX(L,K) = 1 ENDIF ELSE IB = IB + 1 IF ( INDX(L,K) .EQ. 0 ) THEN INDX(L,K) = 1 ENDIF ENDIF ELSE IF (BLEND) THEN BLEND = .FALSE. C C *** Looks if found component can be C *** added to the list of objects. C IF (K .GT. MAX(IDLM,KSAT+2) .AND. & K .LE. LDIF .AND. K .LT. LIM .AND. c & PRFL(L,K) .GT. TRLIM .AND. & INDX(L,K-2) .NE. -1 .AND. & INDX(L,K-1) .NE. -1 .AND. INDX(L,K) & .NE. -1 ) THEN NEW = .TRUE. DO 70 KN = 1,KM IF (ABS(MCM0(2*KN)-L) .LT. 1 .AND. & ABS(MCM0(2*KN+1)-K).LT.3) THEN NEW = .FALSE. ENDIF 70 CONTINUE IF (NEW .AND. KM .LT. 10 .AND. K .GT. 2) THEN C C ***** An information about C ***** new object is written C ***** into array MCM0. C KM = KM + 1 MCM0(2*KM) = L MCM0(2*KM+1) = K - 1 ENDIF ENDIF ENDIF C IF (IB .GT. 0 .AND. INDX(L,K) .EQ. 0 ) THEN INDX(L,K) = 1 IB = IB - 1 ENDIF ENDIF GOTO 61 60 CONTINUE 50 CONTINUE C MCM0(1) = KM IF (LDBG .GT. 3) THEN WRITE(TEXT,'(A,I2,3X,20I3)') ' New objects: ',MCM0 CALL STTPUT(TEXT,ISTAT) ENDIF C C *** Expand a set of points which are affected by C *** components by including also neighbouring points. C DO 340 K = 0 , LIM DO 350 L = 1 , 8 ICPY(L,K) = INDX(L,K) 350 CONTINUE ICPY(0,K) = 0 ICPY(9,K) = 0 340 CONTINUE DO 80 K = 1 , LIM DO 100 L = 1 , 8 IF ( INDX(L,K-1) .EQ. 1 ) THEN ICPY(L-1,K) = 1 ICPY(L+1,K) = 1 IF ( K .GT. 4 ) THEN ICPY(L,K) = 1 ENDIF ENDIF IF ( K .LT. LIM .AND. K .GT. 4 .AND. INDX(L,K+1) & .EQ. 1 ) THEN ICPY(L,K) = 1 ENDIF 100 CONTINUE IF ( INDX(1,K) .NE. -1 ) THEN INDX(1,K) = MAX( ICPY(1,K) , ICPY(9,K) ) ENDIF IF ( INDX(8,K) .NE. -1 ) THEN INDX(8,K) = MAX( ICPY(0,K) , ICPY(8,K) ) ENDIF 80 CONTINUE DO 90 K = 1 , LIM DO 110 L = 2 , 7 IF ( INDX(L,K) .NE. -1 ) THEN INDX(L,K) = ICPY(L,K) ENDIF 110 CONTINUE 90 CONTINUE C C *** Remove low points from a set of affected points. C DO 120 K = 1 , LIM DO 130 L = 1 , 8 TMP = PRFL(L,K) IF ( TMP .LT. ATRSH ) THEN IF ( TMP .GT. BTRSH .AND. ICNT(L,K) .GT. 0 ) THEN INDX(L,K) = 0 ELSE IF ( TMP .LE. BTRSH ) THEN INDX(L,K) = 1 ENDIF ENDIF 130 CONTINUE 120 CONTINUE C C *** Copy the profile into scratch array CPY. C DO 140 L = 0 , LIM DO 150 K = 1 , 8 CPY(K,L) = PRFL(K,L) 150 CONTINUE 140 CONTINUE C C *** Find extend of non cleaned profile KCLN. C K = 0 KCLN = 0 IR = 0 161 CONTINUE IF (IR .NE. 0 .OR. K .GE. LIM) GOTO 160 KCLN = K K = K + 1 IR=MAX(INDX(1,K),INDX(2,K),INDX(3,K),INDX(4,K) , & INDX(5,K),INDX(6,K),INDX(7,K),INDX(8,K)) GOTO 161 160 CONTINUE C C *** Clean the profile. C DO 170 K = 1,LIM IC = 0 SUM = 0 DO 180 L = 1,8 IF (INDX(L,K) .EQ. 0) THEN SUM = SUM + PRFL(L,K) IC = IC + 1 ENDIF 180 CONTINUE IF (IC .GT. 0) THEN AVER = SUM / IC ELSE AVER = MIN(PRFL(1,K),PRFL(2,K) , & PRFL(3,K),PRFL(4,K),PRFL(5,K) , & PRFL(6,K),PRFL(7,K),PRFL(8,K)) END IF DO 190 L = 1,8 IF (INDX(L,K) .EQ. 1) THEN IF (L .LE. 4) THEN L4 = L + 4 ELSE L4 = L - 4 ENDIF IF (INDX(L4,K) .EQ. 0) THEN PRFL(L,K) = PRFL(L4,K) ELSE PRFL(L,K) = AVER ENDIF ENDIF 190 CONTINUE 170 CONTINUE C C ****** Check if cleaning is done properly. C DO 200 L = 1,8 DO 210 K = 1,LIM TEMP(K) = CPY(L,K) - PRFL(L,K) 210 CONTINUE IC = 0 ID = 0 IK = 0 C C Finds if there is a bright neighbour in this octant. C C PMAX=MAX(CPY(L,0),CPY(L,1)) C K=2 C BRIGHT=.FALSE. C DO WHILE (K.LE.LIM.AND.(.NOT.BRIGHT)) C IF (CPY(L,K).GT.PMAX) THEN C BRIGHT=.TRUE. C END IF C K=K+1 C END DO C DO 230 K = 1 , LIM IF ( TEMP(K) .LT. 0.0 .AND. ICNT(L,K) .GT. 0 & .AND. CPY(L,K) .GT. BTRSH ) THEN PRFL(L,K) = CPY(L,K) INDX(L,K) = 0 IC = 0 ID = 0 ELSE IF (TEMP(K).GT.0.0 .AND. IK.EQ.0) THEN IC = IC + 1 IF ( IC .GT. 1 .AND. TEMP(K) .GT. TEMP(K-1) ) THEN IF (ID .EQ. 0) THEN ID = 1 T1(1)=TEMP(K-1)/TEMP(K) ELSE ID = 0 T1(2)=TEMP(K-1)/TEMP(K) IK = K - 3 FCTR = SQRT(T1(1)*T1(2)) FCTR = MIN(BROBJ,FCTR) ENDIF ELSE ID = 0 ENDIF ELSE IC = 0 ID = 0 ENDIF 230 CONTINUE C C ****** Correct central part if bright neighbour. C IF (IK .GT. 0) THEN CORR = TEMP(IK+1) DO 240 K = IK,-1,-1 CORR = CORR * FCTR IF (CORR .GT. 0.3 * TRSH) THEN IF (K .GE. 0) THEN PRFL(L,K)=PRFL(L,K)-CORR ELSE IF (L .LE. 4) THEN L4 = L + 4 ELSE L4 = L - 4 ENDIF KK = ABS(K) PRFL(L4,KK) = PRFL(L4,KK)-CORR ENDIF IF (K .EQ. 0) THEN IF (L .EQ. 1) THEN PRFL(8,1) = PRFL(8,1)-CORR PRFL(2,1) = PRFL(2,1)-CORR ELSE IF (L.EQ.8) THEN PRFL(7,1) = PRFL(7,1)-CORR PRFL(1,1) = PRFL(1,1)-CORR ELSE PRFL(L-1,1) = PRFL(L-1,1)-CORR PRFL(L+1,1) = PRFL(L+1,1)-CORR ENDIF ENDIF ENDIF 240 CONTINUE ENDIF 200 CONTINUE C PR0 = MIN(PRFL(1,0),PRFL(2,0),PRFL(3,0),PRFL(4,0) , & PRFL(5,0),PRFL(6,0),PRFL(7,0),PRFL(8,0)) DO 250 L = 1,8 PRFL(L,0) = PR0 250 CONTINUE C C ****** Write informations about points affected by components. C IF (LDBG .GE. 3) THEN WRITE(TEXT,'(A)') 'Bad points are marked by 1`s' CALL STTPUT(TEXT,ISTAT) JJ1 = MIN(10,LIM-1) DO 260 JJ = 0,JJ1 WRITE(TEXT,'(8I9)') (INDX(II,JJ),II = 1,8) CALL STTPUT(TEXT,ISTAT) 260 CONTINUE ENDIF C C ****** Write cleaned profile to the screen. C IF (LDBG .GE. 3) THEN WRITE(TEXT,'(A)') 'Corrected 8 octant profile' CALL STTPUT(TEXT,ISTAT) K1 = MIN(9,LIM-1) DO 270 K = 0,K1 DO 280 L = 1,8 IOUT(L) = NINT(PRFL(L,K) * SCALE) 280 CONTINUE WRITE(TEXT, '(8I9)') IOUT CALL STTPUT(TEXT,ISTAT) 270 CONTINUE ENDIF C RETURN END