C===========================================================================
C Copyright (C) 1995,2003 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
      PROGRAM AVESUB
C
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C.IDENTIFICATION:                
C  program AVESUB			version 3.00	881028
C  K. Banse				ESO - Garching
C
C.KEYWORDS:
C  statistics, subimage
C
C.PURPOSE:
C  calculate the average value (and rms) according to different methods 
C  of an area defined by the cursor rectangle or a table or by coordinates
C
C.ALGORITHM:
C  This info may also be stored either in a descriptor or a table...
C  
C.INPUT/OUTPUT:
C  
C  The following keywords are used:
C  IN_A/C/1/60				input frame
C  IN_B/C/1/60				input specs,
C					= +, for cursor input
C					= name or name/TABLE for table input
C  ACTION/C/1/1				averaging method,
C					K for kappa-sigma clipping,
C					M for median, A for simple average
C  INPUTC/C/1/50			specs for data storage
C					= :label for table or table,:label (a)
C					= descr for descriptor (b)
C					= + for just display of data
C  INPUTC/C/51/1			(b) append flag for descriptor,
C					= A, append data to descriptor
C					= +, start again at position 1 in descr
C					(a) center flag for table
C					= C, add also columns with center points
C					= +, no
C	
C  OUTPUTR/R/1/5			last displayed results
C  OUTPUTI/I/1/1			no. of loops in program
C	
C.VERSIONS
C  3.00		use new IDI interfaces
C	
C 031021	last modif
C 
C--------------------------------------------------------------------------
C
      IMPLICIT NONE
C      
      INTEGER      FELEM,IAV,INFLAG,ISTOP
      INTEGER      KAPITR,LEAP,NMAL,NROW,NN,N
      INTEGER      OUTFLG,STAT1,STAT2,STAT
      INTEGER      TPIX(2),MAXSIZ
      INTEGER      TBCOLN(7),OUTCOL
      INTEGER      SUBPIX(2,2),IPIXS(4),NPON(5)
      INTEGER      ICUR1(2),ICUR2(2)
      INTEGER*8    PNTRA,PNTRX
      INTEGER      NAXIS,NPIX(2),IMNOA
      INTEGER      TABNUL(7),COOS(4),CURSOF
      INTEGER      UNI(1),NULO,TID
      INTEGER      MADRID(1)
C
      CHARACTER*80 FRAME,TABLE
      CHARACTER    DESCR*48,APPFLG*1
      CHARACTER    CUNIT*48,IDENT*72,CBUF*80,ACTION*1
      CHARACTER    INFO*48
      CHARACTER    OUTUNI*16,OUTLAB(7)*16,DRAWY*8
C      
      REAL         PCUR1(6),PCUR2(6),PCEN(3),RDUM
      REAL         BGR,BGVL(5),RMSV(5),AVAREA
      REAL         REAL1(5),RBUF(10)
C 
      DOUBLE PRECISION START(2),STEP(2)
C
      LOGICAL      SELFLG
C 
      INCLUDE 'MID_INCLUDE:ST_DEF.INC'
C
      INCLUDE 'MID_INCLUDE:IDIDEV.INC'
      INCLUDE 'MID_INCLUDE:IDIMEM.INC'
C
      COMMON   /VMR/ MADRID
C 
      DATA  SUBPIX   /4*1/
      DATA  COOS     /-1,-1,-1,-1/,  CURSOF /0/
      DATA  TID      /-1/
      DATA  INFO     /'switch cursor(s) on - next time we exit... '/
C
      DATA      OUTLAB /'XSTART ','YSTART ','XEND ','YEND ',
     +                  ' ','XCEN ','YCEN '/
      DATA      OUTUNI /' '/
      DATA      REAL1  /5*0./
      DATA      CUNIT  /' '/,    IDENT  /' '/, DESCR /' '/
      DATA      CBUF  /' '/,  FRAME  /' '/, TABLE /' '/
C 
      INCLUDE 'MID_INCLUDE:ST_DAT.INC'
C
C  set up MIDAS environment + enable automatic error abort
      CALL STSPRO('AVESUB')
C
C  get method 
      CALL STKRDC('ACTION',1,1,1,IAV,ACTION,UNI,NULO,STAT)
      IF (INDEX('AKM',ACTION).LE.0) ACTION = 'A'            !default to A(verage)
      CALL STKRDI('INPUTI',1,2,IAV,NPON,UNI,NULO,STAT) 
C      
      IF (ACTION.EQ.'M') THEN
         MAXSIZ = NPON(2)
         CALL STFXMP(MAXSIZ,D_R4_FORMAT,PNTRX,STAT)
      ELSE IF (ACTION.EQ.'K') THEN                    !for kappa-sigma clipping
         KAPITR = NPON(1)
      ENDIF
C      
C  get the subarea input options
      CALL STKRDC('IN_B',1,1,80,IAV,TABLE,UNI,NULO,STAT)
      IF (TABLE(1:1).EQ.'+') THEN
         INFLAG = 1                              !INFLAG = 1 for cursor
         DRAWY(1:7) = 'N YY?C0'
         CALL STKRDC('P4',1,1,1,IAV,DRAWY(2:2),UNI,NULO,STAT)
      ELSE
         CALL STKRDC('IN_A',1,1,80,IAV,FRAME,UNI,NULO,STAT)
         CALL CLNFRA(FRAME,FRAME,0)
         CALL STIGET(FRAME,D_R4_FORMAT,F_I_MODE,F_IMA_TYPE,
     +               2,NAXIS,NPIX,START,STEP,IDENT,
     +               CUNIT,PNTRA,IMNOA,STAT)
         IF (NAXIS.LT.2) 
     +      CALL STETER(77,'1-dim frames currently not supported...')
         CALL CLNTAB(TABLE,TABLE,0)
         INFLAG = 2                              !INFLAG = 2 for table
      ENDIF
C
C  get descriptor or table label for storage of data + output_options
      DESCR(1:1) = '+'
      OUTFLG = 1                          !OUTFLG = 1 for descriptor or nothing
      CALL STKRDC('INPUTC',1,1,51,IAV,CBUF,UNI,NULO,STAT)
      APPFLG = CBUF(51:51)
      IF (APPFLG .EQ. 'c') THEN                !take care of lower case options
         APPFLG = 'C'
      ELSE IF (APPFLG .EQ. 'a') THEN
         APPFLG = 'A'
      ENDIF
      CBUF(51:51) = ' '
      IF (CBUF(1:1).EQ.'+') GOTO 500           !neither table nor descriptor
C      
C  test for :label or table,:label or #no. or table,#no.
      IF (INDEX(CBUF,'#') .GT. 0) 
     +   CALL STETER(73,'We need a column label, no # sign, please.')
C 
      IF (CBUF(1:1).EQ.':') THEN
         N = 1
      ELSE
         N = INDEX(CBUF,',:')           
         IF (N.EQ.1) THEN
            CALL STETER(76,'Invalid syntax (missing table name)...')
         ELSE IF (N.GT.1) THEN
            CALL CLNTAB(CBUF(1:N-1),TABLE,0)      
            N = N + 1                !skip the ','
         ELSE
            GOTO 400                     !it's a descriptor
         ENDIF
      ENDIF
C      
C  handle table storage
      OUTFLG = 2                               !OUTFLG = 2 for table
      OUTLAB(5) = CBUF(N+1:50)
      IF (INFLAG.EQ.1) THEN                    !cursor input
         DO 80, N=1,7
            IF (N.NE.5) THEN
               IF (OUTLAB(5).EQ.OUTLAB(N))
     +            CALL STETER(78,'Invalid column label !')
            ENDIF
80       CONTINUE
C 
C  for cursor input create table with 100 rows
         IF (APPFLG.EQ.'C') THEN
            OUTCOL = 7
         ELSE
            OUTCOL = 5
         ENDIF
         CALL TBTINI(TABLE,0,F_O_MODE,OUTCOL,100,TID,STAT) 
         DO 100, N=1,OUTCOL
            CALL TBCINI(TID,D_R4_FORMAT,1,'E12.5',OUTUNI,OUTLAB(N),
     +                  TBCOLN(N),STAT)
100      CONTINUE
         GOTO 500
      ELSE
         IF (N.GT.1)           !for table input is table,:label not allowed...
     +      CALL STETER(79,'No table name possible in out_specs... ')
      ENDIF
C      
C  for existing table, test, if label(s) already there
      CALL TBTOPN(TABLE,F_IO_MODE,TID,STAT)
      IF (APPFLG.EQ.'C') THEN
         NN = 7
      ELSE
         NN = 5
      ENDIF
      DO 150, N=5,NN
         CALL TBLSER(TID,OUTLAB(N),TBCOLN(N),STAT)
         IF (TBCOLN(N).LE.0)                          !No. We have to define it...
     +      CALL TBCINI(TID,D_R4_FORMAT,1,'E12.5',OUTUNI,OUTLAB(N),
     +                  TBCOLN(N),STAT)
150   CONTINUE
      GOTO 500
C      
C  handle descriptor storage
400   DESCR(1:48) = CBUF(1:48)
      IF (APPFLG.EQ.'A') THEN                   !see, if data is to be appended
         FELEM = -1
      ELSE
         FELEM = 1
      ENDIF
C      
C  branch according to input option
500   NMAL = 0                                  !init counter
      IF (INFLAG.EQ.2) GOTO 2000
C      
C  a) cursor defined subimage(s) - attach Image Display + get all info
C 
      CALL DTOPEN(1,STAT)
      CALL SETCUR(QDSPNO,2,1,2,COOS,STAT)       !set up cursor rectangle
      FRAME(1:) = ' '
C
C  input cursor positions and check status
1000  CALL GETCUR(DRAWY,FRAME,                  !draw box in overlay channel
     +            ICUR1,PCUR1(3),PCUR1(5),RDUM,STAT1,
     +            ICUR2,PCUR2(3),PCUR2(5),RDUM,STAT2)
C
C  check cursor status  
      IF ( (STAT1.EQ.0) .AND. (STAT2.EQ.0) ) THEN
         IF (NMAL.EQ.0) THEN
            IF (CURSOF.EQ.1) GOTO 9000
            CALL STTPUT(INFO,STAT)
            FRAME(1:) = ' '
            CURSOF = 1
            GOTO 1000
         ELSE
            GOTO 9000                             !cursor loop terminated ...
         ENDIF
      ENDIF
C      
      NMAL = NMAL + 1                             !update counter
C      
C  fill output buffer
      REAL1(1) = PCUR1(5)
      REAL1(2) = PCUR1(6)
      REAL1(3) = PCUR2(5)
      REAL1(4) = PCUR2(6)
C      
C  finally map displayed frame
      CALL STIGET(FRAME,D_R4_FORMAT,F_I_MODE,F_IMA_TYPE,2,
     +            NAXIS,NPIX,START,STEP,
     +            IDENT,CUNIT,PNTRA,IMNOA,STAT)
      GOTO 3000
C      
C  b) table defined subimage(s) 
C 
2000  IF (OUTFLG.NE.2) CALL TBTOPN(TABLE,F_IO_MODE,TID,STAT)
      CALL TBIGET(TID,N,NROW,N,N,N,STAT)              !get total no. of rows
      DO 2020, N=1,4
         CALL TBLSER(TID,OUTLAB(N),TBCOLN(N),STAT)
2020  CONTINUE
C      
C  here we loop
2100  IF (NMAL.GE.NROW) GOTO 9000
      NMAL = NMAL + 1
      CALL TBSGET(TID,NMAL,SELFLG,STAT)             !only use selected rows
      IF (.NOT.SELFLG) GOTO 2100
C      
C  get next row of values
      CALL TBRRDR(TID,NMAL,4,TBCOLN,REAL1(1),TABNUL,STAT)
      PCUR1(3) = (REAL1(1)-START(1))/STEP(1) + 1.
      PCUR1(4) = (REAL1(2)-START(2))/STEP(2) + 1.
      PCUR2(3) = (REAL1(3)-START(1))/STEP(1) + 1.
      PCUR2(4) = (REAL1(4)-START(2))/STEP(2) + 1.
C 
C  -----------------------------------
C  common section for option a) and b)
C  -----------------------------------
C 
3000  SUBPIX(1,1) = PCUR1(3)
      SUBPIX(2,1) = PCUR1(4)
      SUBPIX(1,2) = PCUR2(3)
      SUBPIX(2,2) = PCUR2(4)
C      
      IF (APPFLG.EQ.'C') THEN
         N = (SUBPIX(1,1)+SUBPIX(1,2)+1) / 2
         PCEN(2) = (N-1)*STEP(1) + START(1)
         N = (SUBPIX(2,1)+SUBPIX(2,2)+1) / 2
         PCEN(3) = (N-1)*STEP(2) + START(2)
      ENDIF
C      
C  and finally let's do the averaging 
C      
      IF (ACTION.EQ.'K') THEN
         IPIXS(1) = SUBPIX(1,1)
         IPIXS(2) = SUBPIX(1,2)
         IPIXS(3) = SUBPIX(2,1)
         IPIXS(4) = SUBPIX(2,2)
         LEAP = 0
         ISTOP = 0
         CALL BGVAL(1,MADRID(PNTRA),NPIX(1),NPIX(2),IPIXS,KAPITR,
     +              2,LEAP,BGVL,RMSV,NPON,ISTOP,N)
         BGR = BGVL(1)
CC       SIG = RMSV(1)
CC       NPN = NPON(1)
      ELSE
         TPIX(1) = SUBPIX(1,2) - SUBPIX(1,1) + 1            !get size of subimage
         TPIX(2) = SUBPIX(2,2) - SUBPIX(2,1) + 1
         IF ((ACTION.EQ.'M').AND.(TPIX(1)*TPIX(2).GT.MAXSIZ)) THEN
            CALL STTPUT('subimage too large for internal buffers -'//
     +                  ' skipped...',STAT)
            IF (INFLAG.EQ.1) NMAL = NMAL - 1                !modify counter
            GOTO 7000
         ENDIF
         BGR = AVAREA(ACTION,MADRID(PNTRA),NPIX(1),SUBPIX,TPIX,
     +                MADRID(PNTRX))
      ENDIF
C
C  display results
      IF (NMAL.EQ.1) THEN                               !here for the 1. time...?
         CBUF(1:) = '   xstart        ystart          xend           '//
     +              'yend        average_val'
         CALL STTPUT(CBUF,STAT)
      ENDIF
      REAL1(5) = BGR
      WRITE(CBUF,20000) (REAL1(N),N=1,5)
      CALL STTPUT(CBUF,STAT)
C
C  fill table, if applicable
      IF (OUTFLG.EQ.2) THEN
         IF (INFLAG.EQ.2) THEN                  !table input
            PCEN(1) = BGR
            IF (APPFLG.EQ.'C') THEN
               N = 3
            ELSE
               N = 1
            ENDIF                               !row already exists
            CALL TBRWRR(TID,NMAL,N,TBCOLN(5),PCEN,STAT)
         ELSE                                    !cursor input
            DO 3500, N=1,5
               RBUF(N) = REAL1(N)
3500        CONTINUE
            IF (APPFLG.EQ.'C') THEN
               RBUF(6) = PCEN(2)
               RBUF(7) = PCEN(3)
            ENDIF
            CALL TBRWRR(TID,NMAL,OUTCOL,TBCOLN,RBUF,STAT)
         ENDIF
      ELSE
C      
C  fill descriptor, if applicable
         IF (DESCR(1:1).NE.'+') THEN
            CALL STDWRR(IMNOA,DESCR,REAL1,FELEM,5,UNI,STAT) 
            FELEM = -1
         ENDIF
      ENDIF
C      
C  store info also in keyword OUTPUTR
      CALL STKWRR('OUTPUTR',REAL1(1),1,5,UNI,STAT)

C      
C  and loop again
7000  GOTO (1000,2100),INFLAG
C  
C  That's it folks...
9000  IF ((TID.NE.-1) .AND. (OUTFLG.EQ.2)) THEN
         CALL TBSINI(TID,STAT)
         CALL TBTCLO(TID,STAT)
      ENDIF
C      
C  save no. of coordinates obtained for subsequent applications
      CALL STKWRI('OUTPUTI',NMAL,1,1,UNI,STAT)
C      
      IF (INFLAG.EQ.1) CALL DTCLOS(QDSPNO)
      CALL STSEPI
C      
C  format statements
20000 FORMAT(1X,5G15.7)
C 
      END         
