C @(#)clasfy.for	19.1 (ES0-DMD) 02/25/03 13:25:32
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.COPYRIGHT: Copyright (c) 1987 European Southern Observatory,
C                                         all rights reserved
C.IDENT:    subroutine CLASFY         version 1.2        870707
C           A. Kruszewski             ESO Garching
C.VERSION:  1.0  ESO-FORTRAN Conversion, AA  15:03 - 21 NOV 1987
C.LANGUAGE: F77+ESOext
C.AUTHOR:   A.KRUSZEWSKI
C.PURPOSE:
C.INPUT/OUTPUT:
C  input arguments
C  ICLS        integer*4 array      holds seed and current 
C  PCLA        real*4 array         array holding classifiers
C  BRGT        logical array        array indicating saturated objects
C  M           integer*4            number of objects
C  output arguments
C  PCLA        real*4 array      array holding corrected classifiers
C  ITER        integer*4         number of iteration
C  DONE        logical           convergence indicator
C-----------------------------------------------------------------------
C
      SUBROUTINE CLASFY(PCLA,BRGT,ICLS,M,ITER,DONE)
C
      IMPLICIT NONE
      INCLUDE 'MID_REL_INCL:INVENT.INC/NOLIST'
C
      INTEGER      I, IC, IN
      INTEGER      IW(0:MAXCNT)
      INTEGER      ICTL(7)
      INTEGER      ICLS(2,MAXCNT)
      INTEGER      ITLM(2)
      INTEGER      ITER
      INTEGER      IAC
      INTEGER      ICHG
      INTEGER      ISTAT
      INTEGER      IVAR
      INTEGER      M
      INTEGER      MADRID(1)
      INTEGER      J
      INTEGER      K
      INTEGER      KUN, KNUL
      INTEGER      L
      LOGICAL      BRGT(MAXCNT)
      LOGICAL      DONE
      REAL         DVGR(0:MAXCNT)
C     REAL         TINFO(13)
C     REAL         CLPR(3)
      REAL         CLDV(13,MAXCNT)
      REAL         X(0:MAXCNT)
      REAL         Y(0:MAXCNT)
      REAL         COVM(13,13)
      REAL         DIST(2,MAXCNT)
      REAL         PCLA(0:13,MAXCNT)
      REAL         SGLM
      REAL         SIGMA
      REAL         ZRLV
      REAL         GRDT
      REAL         SDIS
      REAL         GDIS
      REAL         TEMP
      CHARACTER*80 OUTPUT
C
      COMMON       /VMR/MADRID
C
      CALL STKRDR( 'INV_REAL' , 9 , 1 , IAC , SGLM , KUN , KNUL , ISTAT)
      CALL STKRDI( 'INV_INTG' , 17 , 2 , IAC , ITLM , KUN , KNUL ,ISTAT)
C
C         Stars and galaxies are treated seperately
C ****    (L=1 stars, L=2 galaxies) using current
C         classification resting in ICLS(2,*).
C
      DO 130 L = 1 , 2
C
C            Perform regression of parameters onto logarithm
C ****       of central brightness. At first the objects
C            of considered class are marked in array IW.
C
      IW(0)  = 0
      IN     = 0
      DO 10 I = 1,M
          IF ( .NOT. BRGT(I) .AND. ICLS(2,I).EQ.L) THEN
              IW(I)  = 1
              IN     = IN + 1

          ELSE
              IW(I)  = 0
          END IF

          X(I)   = PCLA(0,I) - 19.0
   10 CONTINUE
      IF (IN.GT.15) THEN
          DO 40 J = 1,13
              DO 20 I = 1,M
                  Y(I)   = PCLA(J,I)
   20         CONTINUE
              CALL GRADET(M,X,Y,IW,ZRLV,GRDT,DVGR,SIGMA)
              DO 30 I = 1,M
                  CLDV(J,I) = DVGR(I)
   30         CONTINUE
   40     CONTINUE
C
C       calculates covariance matrix COVM of parameter's deviations CLDV
C
          DO 60 J = 1,13
              DO 50 I = 1,13
                  COVM(I,J) = 0.0
   50         CONTINUE
   60     CONTINUE
          DO 90 J = 1,13
              DO 80 I = 1,13
                  IC     = 0
                  DO 70 K = 1,M
                      IF (IW(K).EQ.1) THEN
                          COVM(I,J) = COVM(I,J) + CLDV(I,K)*CLDV(J,K)
                          IC     = IC + 1
                      END IF

   70             CONTINUE
                  IF (L.EQ.1) THEN
                      ICTL(1) = IC

                  ELSE IF (L.EQ.2) THEN
                      ICTL(2) = IC
                  END IF

                  COVM(I,J) = COVM(I,J)/ (IC-13)
   80         CONTINUE
   90     CONTINUE
C            DO6912 I=1,13
C               TINFO(I)=COVM(I,I)
C            END DO
C
C       inverts covariance matrix
C
          CALL INVERT(COVM,13)
C
C       we can now calculate distances of particular objects from
C       the center of L-th class in the parameters space
C
          DO 120 K = 1,M
              IF ( .NOT. BRGT(K)) THEN
                  TEMP   = 0.0
                  DO 110 J = 1,13
                      DO 100 I = 1,13
                          TEMP   = TEMP + COVM(I,J)*CLDV(I,K)*CLDV(J,K)
  100                 CONTINUE
  110             CONTINUE
                  DIST(L,K) = TEMP
              END IF

  120     CONTINUE
      END IF

  130 CONTINUE
C
C classification is performed by comparing distances in parameter space
C    from centers of two classes of objects
C
      ICHG   = 0
      DO 140 K = 1,M
          IF ( .NOT. BRGT(K)) THEN
              SDIS   = DIST(1,K)
              GDIS   = DIST(2,K)
              IVAR   = ICLS(2,K)
              IF (SDIS.GT.SGLM .AND. GDIS.GT.SGLM) THEN
                  ICLS(2,K) = 0

              ELSE
                  IF (SDIS.LT.GDIS) THEN
                      ICLS(2,K) = 1

                  ELSE
                      ICLS(2,K) = 2
                  END IF

              END IF

              IF (ICLS(2,K).NE.IVAR) THEN
                  ICHG   = ICHG + 1
              END IF

          END IF

  140 CONTINUE
      ITER   = ITER + 1
      IF (ICHG.LE.ITLM(1) .OR. ITER.GT.ITLM(2)) THEN
          DONE   = .TRUE.
      END IF

      ICTL(3) = ICHG
      DO 150 I = 4,7
          ICTL(I) = 0
  150 CONTINUE
      WRITE (OUTPUT,'(3(A,I6))') ' Stars: ',ICTL(1),'   Galaxies: ',
     +  ICTL(2),'   Changes: ',ICTL(3)
      CALL STTPUT(OUTPUT,ISTAT)
C
      RETURN

      END

