C @(#)appmtr.for	19.1 (ES0-DMD) 02/25/03 13:25:31
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
C.IDENTIFICATION
C
C  subroutine APPMTR         version 1.2       831004
C  A. Kruszewski             ESO Garching
C
C  modified                  version 1.3       870127
C  A. Kruszewski             Obs. de Geneve
C
C.KEYWORDS
C
C  classifiers, magnitudes
C
C.PURPOSE
C
C  calculates various kinds of magnitudes and sizes for an object
C  with given average profile "AVPR"
C
C.INPUT/OUTPUT
C
C  input arguments
C
C  AVPR        real*4 array      one-dimensional profile
C  LIM         integer*4         radial extend of determined profile
C  ISOP        integer*4         isophotal size flag
C  RARR        real*4 array      real keywords array
C  PCTL        real*4 array      array holding classifiers
C
C  output arguments
C
C  AVPR        real*4 array      modified input array
C  PCTL        real*4 array      array holding classifiers
C
C-----------------------------------------------------------------------
      SUBROUTINE APPMTR(AVPR, LIM, ISOP, RARR, APSF, PCTL)
C
      IMPLICIT  NONE
      INCLUDE   'MID_REL_INCL:INVENT.INC/NOLIST'
C
      REAL      AVPR(0:MAXSUB)
      INTEGER   LIM
      INTEGER   ISOP
      REAL      RARR(64)
      REAL      APSF(0:MAXSUB)
      REAL      PCTL(NRPAR)
C
      INTEGER   IC
      INTEGER   IW( 0:3 )
      INTEGER   K, KM
      INTEGER   L, LL
      INTEGER   NSUM
      INTEGER   NPXL( 0:MAXSUB )
      REAL      ALP , AKML, AKM
      REAL      AVPR0, AVNPK
      REAL      BSUM, BSUM1
      REAL      DETA
      REAL      DFAV
      REAL      DVGR ( 0:3 )
      REAL      ETA, ETA1
      REAL      FRCT
      REAL      PETA
      REAL      RISP
      REAL      SBRH( 0:MAXSUB )
      REAL      SIGMA
      REAL      SMTR( 2 )
      REAL      STMT
      REAL      SUMM, SUM0, SUM1, SUM2, SUM3
      REAL      TRSH
      REAL      VAR
      REAL      XRAD( 0:3 ) , YMAG( 0:3 )
C     REAL      X(5)
      REAL      ZRLV 
      REAL      ZRMG
C
      DATA NPXL/1,8,12,16,32,28,40,40,48,68,56,72,68,88,88,84,112,
     &          112,112,116,112,144,140,144,144,168,164,160,184,172,
     &          200,192,188,208,224,224,228,224,248,236,264,248,264,
     &          276,264,288,276,304,304,312,316/
C
C ***
      TRSH = RARR(3)
      ZRMG = RARR(11)
      SMTR(1) = RARR(12)
      SMTR(2) = RARR(13)
      PETA = RARR(53)
      RISP = PCTL(15)
      ETA1 = 0.0
      BSUM1 = 0.0
C
C *** Calculate quantities integrated down to isophote TRSH.
C
      AVPR0 = AVPR(0)
      SUM0 = AVPR0
      SUM1 = 0.333 * AVPR0
      SUM2 = 0.1 * AVPR0
      SUM3 = 10.0 * AVPR0
      K = 1
   11 CONTINUE
          IF ( AVPR(K) .LT. TRSH .OR. K .GE. LIM ) GOTO 10
          AVNPK = AVPR(K) * NPXL(K)
          SUM0 = SUM0 + AVNPK
          SUM1 = SUM1 + AVNPK * K
          SUM2 = SUM2 + AVNPK * K * K
          SUM3 = SUM3 + AVNPK / ( K * K )
          K = K + 1
          GOTO 11
   10 CONTINUE
      DFAV = AVPR(K-1) - AVPR(K)
      IF ( DFAV .GT. 0.01*TRSH .AND. K .LT. LIM ) THEN
          FRCT = ( AVPR(K-1) - TRSH ) / DFAV
      ELSE
          FRCT = 1.0
      ENDIF
      AVNPK = AVPR(K) * NPXL(K)
      SUM0 = SUM0 + FRCT * AVNPK
      SUM1 = SUM1 + FRCT * AVNPK * K
      SUM2 = SUM2 + FRCT * AVNPK * K * K
      SUM3 = SUM3 + FRCT * AVNPK / ( K * K )
      IF ( SUM0 .GT. 0.1 * TRSH ) THEN
          PCTL(5) = ZRMG - 2.5 * LOG10(SUM0)
          PCTL(6) = SUM1 / SUM0
          PCTL(7) = SUM2 / SUM0
          PCTL(18) = SUM3 / SUM0
          PCTL(15) = FLOAT(K) - 1.0 + FRCT
      ELSE
          PCTL(5) = ZRMG - 2.5 * ( LOG10(TRSH) - 1.0 )
      ENDIF
      IF ( PCTL(7) .GT. 0.001 ) THEN
          PCTL(7) = SQRT( PCTL(7) )
      ELSE
          PCTL(7) = 0.0
      ENDIF
      IF ( PCTL(18) .GT. 0.001 ) THEN
          PCTL(18) = ( PCTL(18) ) ** (-0.5)
      ELSE
          PCTL(18) = 0.0
      ENDIF
C
C *** If control parameter ISOP is equal 1 then isophotal
C *** magnitude is calculated inside the isophotal radius RISP.
C
      IF ( ISOP .EQ. 1 ) THEN
          AVPR0 = AVPR(0)
          SUM0 = AVPR0
          K = 1
   21     CONTINUE
              IF ( FLOAT(K) .GE. RISP .OR. K .GE. LIM ) GOTO 20 
              AVNPK = AVPR(K) * NPXL(K)
              SUM0 = SUM0 + AVNPK
              K = K + 1
              GOTO 21
   20     CONTINUE
          FRCT = RISP - K + 1
          AVNPK = AVPR(K) * NPXL(K)
          SUM0 = SUM0 + FRCT * AVNPK
          IF ( SUM0 .GT. 0.1 * TRSH ) THEN
              PCTL(5) = ZRMG - 2.5 * LOG10(SUM0)
          ELSE
              PCTL(5) = ZRMG - 2.5 * ( LOG10(TRSH) - 1.0 )
          ENDIF
          PCTL(15) = RISP
      ENDIF
C
C *** Calculate metric magnitudes and logarithmic
C *** derivative ALPHA at the first metric size. 
C
      DO 30 LL = 1 , 2
          STMT = SMTR(LL)
C
C ***     Calculate smallest radius KM out of 4 rings used
C ***     to determine metric magnitude and gradient alpha.
C
          VAR = STMT - 1.49
          KM = MAX( MIN( INT(VAR) , LIM-3 ) , 0 )
          IF ( KM .GE. 0 ) THEN
C
C ***         Calculate how much is standard
C ***         metric STMT shifted from KM.
C
              AKM = STMT - 0.5 - KM
C
C ***         Calculates cumulative magnitudes YMAG
C ***         as function of logarithmic ring radius
C ***         XRAD for four radii starting with KM.
C
              SUMM = 0.0
              DO 40 L = 0 , KM
                  SUMM = SUMM + NPXL(L) * AVPR(L)
   40         CONTINUE
              IC = 0
              IF ( SUMM .GT. 0.1*TRSH ) THEN
                  YMAG(0) = ZRMG - 2.5 * LOG10(SUMM)
                  IF ( KM .GT. 0 ) THEN
                      XRAD(0) = LOG10( FLOAT(KM) )
                  ELSE
                      XRAD(0) = -0.5
                  ENDIF
              ELSE
                  IC = 1
              ENDIF
              DO 50 L = 1 , 3
                  SUMM = SUMM + NPXL(KM+L) * AVPR(KM+L)
                  IF ( SUMM .GT. 0.1*TRSH ) THEN
                      YMAG(L) = ZRMG - 2.5 * LOG10(SUMM)
                      XRAD(L) = LOG10( FLOAT(KM+L) )
                  ELSE
                      IC=1
                  ENDIF
   50         CONTINUE
              DO 60 L = 0 , 3
                  IW(L) = 1
   60         CONTINUE
              IF ( IC .EQ. 0 ) THEN
                  CALL GRADET( 3 , XRAD , YMAG , IW , ALP ,
     &                         ZRLV , DVGR , SIGMA )
                  AKML = LOG10( FLOAT(KM) + AKM )
                  IF ( LL .EQ. 1 ) THEN
                      PCTL(13) = ZRLV + AKML * ALP
                  ELSE
                      PCTL(20) = ZRLV + AKML * ALP
                  ENDIF
              ELSE
                  IF ( LL .EQ. 1 ) THEN
                      PCTL(13) = ZRMG + 50.0
                      ALP = 0.0
                  ELSE
                      PCTL(20) = ZRMG + 50.0
                  ENDIF
              ENDIF
              IF ( LL .EQ. 1 ) THEN
                  PCTL(19) = -0.4 * ALP
              ENDIF
          ELSE
              IF ( LL .EQ. 1 ) THEN
                  PCTL(13) = ZRMG + 50.0
                  PCTL(19) = 0.0
              ELSE
                  PCTL(20) = ZRMG + 50.0
              ENDIF
          ENDIF
   30 CONTINUE
C
C *** Calculate Petrosian magnitude and radius.
C
      NSUM = 1
      BSUM = AVPR(0)
      SBRH(0) = AVPR(0)
      ETA = 0.0
      K = 1
   71 CONTINUE
          IF ( ETA .GT. PETA .OR. K .GE. LIM ) GOTO 70
          NSUM = NSUM + NPXL(K)
          BSUM = BSUM + NPXL(K) * AVPR(K)
          SBRH(K) = BSUM / NSUM
          IF ( MIN( SBRH(K)+SBRH(K-1) , AVPR(K) ) .GT. 0.1*TRSH ) THEN
              ETA = 2.5 * LOG10( 0.5 * (SBRH(K)+SBRH(K-1)) / AVPR(K) )
          ELSE
              ETA = PETA + 1.0
          ENDIF
          IF ( ETA .LT. PETA .AND. K .LT. LIM-1 ) THEN
              BSUM1 = BSUM
              ETA1 = ETA
          ENDIF
          K = K + 1
          GOTO 71
   70 CONTINUE
      DETA = ETA - ETA1
      IF ( DETA .GT. 0.001 ) THEN
          FRCT = (PETA-ETA1) / DETA
      ELSE
          FRCT = 0.5
      ENDIF
      IF ( K .GE. LIM .AND. FRCT .GT. 1.0 ) THEN
          FRCT = 1.0
      ENDIF
      BSUM = BSUM1 + FRCT * NPXL(K-1) * AVPR(K-1)
      PCTL(17) = FLOAT(K) - 2.0 + FRCT
      IF ( BSUM .GT. 0.1*TRSH ) THEN
          PCTL(16) = ZRMG - 2.5 * LOG10(BSUM)
      ELSE
          PCTL(16) = ZRMG + 50.0
      ENDIF
C
      RETURN
C
      END
C
