C @(#)prflcl.for	19.1 (ES0-DMD) 02/25/03 13:25:37
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
