C @(#)tdaver.for	19.1 (ESO-DMD) 02/25/03 14:11: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.COPYRIGHT: Copyright (c) 1987 European Southern Observatory,
C                                         all rights reserved
C
C.VERSION: 1.1  ESO-FORTRAN Conversion, AA  14:09 - 19 NOV 1987
C
C.LANGUAGE: F77+ESOext
C
C.AUTHOR: J.D.PONZ
C  900219  KB, throw out SX calls...
C  900320  MP modify call to STIGET
C
C
C
C.IDENTIFICATION:
C  PROGRAM TDAVER
C
C.KEYWORDS:
C  LOCAL  STATISTICS
C
C.PURPOSE:
C  calculate the mean value and rms of an area defined by
C  the table positions.
C
C  Implements the command
C
C  AVERAGE/TAB image table x[,y] z [size]
C
C 020807		last modif
C 
C -------------------------------------------------------------------
C
C
      SUBROUTINE TDAVER
      IMPLICIT NONE
C
      INTEGER     MADRID
      INTEGER     ICOL(2),OCOL,NPIX(2)
      INTEGER     ISTAT, IAV, NP, TID, NCOL, NROW, NSC, NAC, NAR
      INTEGER     IND1, IND2, NINCOL, LEN, NAXIS, I
      INTEGER*8   IPNTR
      INTEGER     IX, IX1, IX2, IY, IY1, IY2, NT
      INTEGER     INDEX,KNUL,KUN, TYPE, TYPE1
      INTEGER     IMNO
C
      REAL VALUE(2), BGVL,STARTA(2),STEPA(2)
      DOUBLE PRECISION DSTARTA(2),DSTEPA(2) 
      LOGICAL NULL(2),SELECT
C
      CHARACTER*80 TABLE, LINE
      CHARACTER*16 LABEL, UNIT
      CHARACTER CUNIT*64,IDENT*72
      CHARACTER FRAME*64,ICOLR*17,OCOLR*17
      CHARACTER FORM*8
      CHARACTER FORM1*8
C
      INCLUDE 'MID_INCLUDE:TABLES.INC'
      COMMON /VMR/MADRID(1)
      INCLUDE 'MID_INCLUDE:TABLED.INC'
C
      DATA      UNIT/'DN '/, FORM/'G15.7 '/
      DATA      NPIX /1,1/, DSTARTA /0.,0./, DSTEPA /1.,1./
C
      TYPE = D_R4_FORMAT
C
C ... INITIALIZE MIDAS
C
      CALL STKRDI('INPUTI',1,1,IAV,NP,KUN,KNUL,ISTAT)
C
C ... GET PARAMS
C
      FRAME(1:)  = ' '
      LINE(1:)   = ' '
      TABLE(1:)  = ' '
      CALL STKRDC('P1',1,1,64,IAV,FRAME,KUN,KNUL,ISTAT)
      CALL STKRDC('P2',1,1,64,IAV,TABLE,KUN,KNUL,ISTAT)
      CALL STKRDC('P3',1,1,36,IAV,LINE,KUN,KNUL,ISTAT)
      CALL STKRDC('P4',1,1,17,IAV,OCOLR,KUN,KNUL,ISTAT)
C
C ... READ TABLE
C
      TID = -1
      CALL TBTOPN(TABLE,F_U_MODE,TID,ISTAT)
      CALL TBIGET(TID,NCOL,NROW,NSC,NAC,NAR,ISTAT)
C
C ... DECODE PARAMETERS AND FIND TABLE COLUMNS
C
      IND1   = INDEX(LINE,',') - 1
      IND2   = INDEX(LINE,' ') - 1
      NINCOL = 2
      IF (IND1.EQ.-1) THEN
          IND1   = INDEX(LINE,' ') - 1
          NINCOL = 1
          NPIX(2) = 1
      END IF
      ICOLR = LINE(1:IND1)
      CALL TBCSER(TID,ICOLR,ICOL(1),ISTAT)
      IF (ICOL(1).LE.0) THEN
          CALL STTPUT(' input column not found ...',ISTAT)
          GO TO 20
      END IF
      CALL TBFGET(TID,ICOL(1),FORM1,LEN,TYPE1,ISTAT)
      IF (TYPE1.EQ.D_C_FORMAT) THEN
          CALL STTPUT(' wrong input column format ...',ISTAT)
          GO TO 20
      END IF
      IF (NINCOL.EQ.2) THEN
          ICOLR = LINE(IND1+2:IND2)
          CALL TBCSER(TID,ICOLR,ICOL(2),ISTAT)
          IF (ICOL(2).LE.0) THEN
              CALL STTPUT(' input column not found ...',ISTAT)
              GO TO 20
          END IF
          CALL TBFGET(TID,ICOL(2),FORM1,LEN,TYPE1,ISTAT)
          IF (TYPE1.EQ.D_C_FORMAT) THEN
              CALL STTPUT(' wrong input column format ...',ISTAT)
              GO TO 20
          END IF
      END IF
      CALL TBCSER(TID,OCOLR,OCOL,ISTAT)
      IF (OCOL.LE.0) THEN
          LABEL  = OCOLR(2:)
          CALL TBCINI(TID,TYPE,1,FORM,UNIT,LABEL,OCOL,ISTAT)
      ELSE
          CALL TBFGET(TID,OCOL,FORM1,LEN,TYPE1,ISTAT)
          IF (TYPE1.EQ.D_C_FORMAT) THEN
              CALL STTPUT(' wrong output column format ...',ISTAT)
              GO TO 20
          END IF
      END IF
C
C ... READ IMAGE
C
      CALL STIGET(FRAME,D_R4_FORMAT,F_I_MODE,F_IMA_TYPE,
     +            2,NAXIS,NPIX,DSTARTA,DSTEPA,IDENT,CUNIT,
     +            IPNTR,IMNO,ISTAT)
C
C ... ITERATION ON POSITIONS
C
      DO 30 I =1,2
        STARTA(I) = DSTARTA(I)
        STEPA(I) = DSTEPA(I)
30    CONTINUE
      NULL(2) = .FALSE.
      DO 10 I = 1,NROW
C
C ...     READ VALUES AND SELECT FLAG
C
          CALL TBRRDR(TID,I,NINCOL,ICOL,VALUE,NULL,ISTAT)
          CALL TBSGET(TID,I,SELECT,ISTAT)
          IF (SELECT .AND. (.NOT.NULL(1)) .AND. (.NOT.NULL(2))) THEN
              IX     = (VALUE(1)-STARTA(1))/STEPA(1) + 1
              IX1    = MAX(IX-NP,1)
              IX2    = MIN(IX+NP,NPIX(1))
              IF (NINCOL.EQ.2) THEN
                  IY     = (VALUE(2)-STARTA(2))/STEPA(2) + 1
                  IY1    = MAX(IY-NP,1)
                  IY2    = MIN(IY+NP,NPIX(2))
              ELSE
                  IY     = 1
                  IY1    = 1
                  IY2    = 1
              END IF
C
C ...       COMPUTE MEAN VALUE AND WRITE RESULT
C
              IF (IX.GE.1 .AND. IX.LE.NPIX(1) .AND. IY.GE.1 .AND. IY.LE.
     +            NPIX(2)) THEN
                  NT     = (IX2-IX1+1)* (IY2-IY1+1)
                  CALL MEAN(MADRID(IPNTR),NPIX(1),NPIX(2),IX1,IX2,IY1,
     +                      IY2,NT,BGVL)
                  CALL TBEWRR(TID,I,OCOL,BGVL,ISTAT)
              ELSE
                  CALL TBEDEL(TID,I,OCOL,ISTAT)
              END IF

          ELSE
              CALL TBEDEL(TID,I,OCOL,ISTAT)
          END IF
   10 CONTINUE
C
C ... HISTORY AND BYE
C
      CALL TDHSTR(TID,ISTAT)
   20 CONTINUE
      CALL TBTCLO(TID,ISTAT)
      RETURN
      END

      SUBROUTINE MEAN(A,NPIX1,NPIX2,IXA,IXE,IYA,IYE,NP,VAL)
C
C COMPUTE THE MEAN VALUE OF THE SQUARE CENTERED AT (IX, IY),
C WITH DIMENSION 2*NP+1 PIXELS
C
      IMPLICIT NONE
C
      INTEGER   NPIX1, NPIX2, IXA, IXE, IYA, IYE, NP
      INTEGER   I, J
C
      REAL      VAL, A(NPIX1,NPIX2)
C
      DOUBLE PRECISION SMEAN
C
      SMEAN  = 0.D0
      DO 20 I = IYA,IYE
          DO 10 J = IXA,IXE
              SMEAN  = SMEAN + A(J,I)
   10     CONTINUE
   20 CONTINUE
      VAL    = SMEAN/NP
      RETURN
      END
