C===========================================================================
C Copyright (C) 1995-2005 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 STATIS
C
C++++++++++++++++++++++++++++++++++++++++++++++++++
C
C.LANGUAGE: F77+ESOext
C
C.AUTHOR: K.Banse
C
C.IDENTIFICATION
C  program STATIS	1.00	881028
C  K. Banse                    	ESO - Garching
C
C.KEYWORDS
C  maximum, minimum,mean value, std. dev., histogram
C
C.PURPOSE
C  calculate statistics of given bulk data frame
C
C.ALGORITHM
C  prepare input for options: COLUMN, ROW, TABLE, CURSOR 
C  do the actual calculations in subroutine ZMSTAT 
C  (former main program STATIST)
C
C.INPUT/OUTPUT
C  the following keys are used for I/O:
C
C  IN_A/C/1/120		input frame
C  P2/C/1/120           subframe area [...:...]
C           		+ for complete frame
C			ROW to work on rows of 2-dim frame
C			COLUMN to work on columns of 2-dim frame
C			PLANE to work on planes of 3-dim frame
C			table_name for table input
C			CURSOR for cursor defined subimages
C  DEFAULT/C/1/6	(1) = N, subarea is specified in P2
C                           = M, subareas are specified in P2 
C				 but merge them all together and do the
C				 statistics on that set
C			    = Y, use complete frame
C    			(2) = N, binsize is given in INPUTR(1)
C                           = Y, total no. of bins is given in INPUTR(1)
C			(3,4) = option as defined in command STATIST
C			(5) = Y/N, draw option for rectangles, only used 
C			      for CURSOR above
C    			(6) = P/N, for plotting histogram or not
C  OUT_A/C/1/120	for P2 = ROW, COLUMN, table_name or CURSOR
C			this holds the output_table, append_flag
C
C  OUTPUTR/R/1/12	will receive: Min,Max,Mean,Std,3rdMom,4thMom,Intensity
C                  	and Median,SmallestMode,Nobins,Binsize,Mode
C  OUTPUTI/I/1/7	(1) = no. valid of pixels in window
C               	(2,3,4) = window start pixels, set to 0, if not applic.
C               	(5,6,7) = window end pixels, set to 0, if not applic.
C
C.VERSIONS
C 
C 051111		last modif
C
C--------------------------------------------------
C
      IMPLICIT NONE
C
      INTEGER      IAV,STAT,N,NLOOP,CURMAX
      INTEGER      ST1,ST2,XY1(2),XY2(2),COOS(4)
      INTEGER      NPIX(3),NAXPIX(4),NMAL,NROW,NCOL1,SIZ(2)
      INTEGER*8    PNTRA,PNTRB,PNTRO
      INTEGER      SPIX(3),EPIX(3)
      INTEGER      TCOLNM(16),ICOLNM(6)
      INTEGER      IACTIO,COOFF,TBAPP
      INTEGER      TABNUL(6)
      INTEGER      UNI(1),NULO,MADRID(1)
      INTEGER      IMNO,IMNOB,IMNOC,INTID,OUTID,OUTNO
      INTEGER      INCOL,KBUF(25)
      INTEGER      NPXO(2),EXAMED
      INTEGER      EC,ED,EL
C
      DOUBLE PRECISION START(3),STEP(3),DVAL
      DOUBLE PRECISION STRTO(2),STPO(2)
      DOUBLE PRECISION DD1(4),DD2(4)
C 
      REAL         RVAL(16),XMEDIA,ZBINS(3)
C
      CHARACTER    ORFILE*120,INTABL*120,OUTTBL*120,OUTFRA*120
      CHARACTER    AREA*120,DEFAUL*6,OUTPUT*120,CCBUF*40
      CHARACTER    OPTION*5,FRMSTR*20,HISTORY*160
      CHARACTER    TABUNT*16,TABLAB(12)*16,CACT*8
      CHARACTER    INTLAB(6)*16
      CHARACTER    IDENT*72,CUNIT*64,MIDUNI*2,TMPFIL*14
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
      COMMON  /PLOTTO/  IDENT,CUNIT
C
      DATA      COOS     /-1,-1,-1,-1/
      DATA      RVAL     /16*0.0/
      DATA      TABUNT   /' '/
      DATA      TABLAB   /'MIN ','MAX ','MEAN ','STD ',
     +                    'MOM3 ','MOM4 ','TOT_INTENS',
     +                    'MEDIAN ','SMAL_MODE ',
     +                    'NO_BINS ','BIN_SIZE ', 'MODE '/
      DATA      INTLAB   /'XSTART ','YSTART ','XEND ','YEND ',
     +                    'ZSTART ','ZEND '/
      DATA      AREA     /' '/
      DATA      IDENT    /' '/, CUNIT /' '/
C
      INCLUDE 'MID_INCLUDE:ST_DAT.INC'
C 
C  init MIDAS environment + get current MIDAS unit
C
      CALL STSPRO('STATIST')
      CALL STPSET(F_FITS_PARM,1,STAT)              !ignore ESO.xyz FITS keywds
      CALL STKWRR('OUTPUTR',RVAL,1,14,UNI,STAT)    !clear keyword OUTPUTR
      CALL STKRDC('MID$SESS',1,11,2,IAV,MIDUNI,UNI,NULO,STAT)
      TMPFIL = 'stat'//MIDUNI//'dum.bdf '
C 
      DO 500, N=1,3
         NPIX(N) = 1
         SPIX(N) = 1
         EPIX(N) = 1
         START(N) = 0.D0
         STEP(N) = 1.D0
500   CONTINUE
      OUTID = -999                     !indicate that no output table/image
      OUTNO = -999                     !indicate that no output image/image
      INTID = -1
      FRMSTR(1:) = ' '
      EXAMED = 0
      KBUF(1) = -1                     !indicates an invalid row/column/plane
C
C  get info about histogram bins + the different options + format string
C 
      CALL STKRDC('DEFAULT',1,1,6,IAV,DEFAUL,UNI,NULO,STAT)
      CALL UPCAS(DEFAUL,DEFAUL)
      OPTION(1:4) = DEFAUL(1:4)
      OUTPUT(1:) = ' '
      CALL STKRDC('P3',1,1,20,IAV,OUTPUT,UNI,NULO,STAT)
      IF (OUTPUT(1:1) .EQ. '#') THEN
         OPTION(2:2) = 'Y'
         CALL GENCNV(OUTPUT(2:),2,1,ICOLNM,ZBINS,DVAL,IAV)
         IF (IAV .NE. 1) CALL STETER(34,'Invalid no. of bins...')
      ELSE
         OPTION(2:2) = 'N'
         CALL GENCNV(OUTPUT(1:),2,1,ICOLNM,ZBINS,DVAL,IAV)
         IF (IAV .NE. 1) CALL STETER(34,'Invalid binsize...')
      ENDIF
      OUTPUT(1:) = ' '
      CALL STKRDC('P4',1,1,120,IAV,OUTPUT,UNI,NULO,STAT)
      CALL GENCNV(OUTPUT,2,2,ICOLNM,ZBINS(2),DVAL,IAV)
      IF (IAV .NE. 2) CALL STETER(35,'Invalid excess bins...')
      CALL STKRDC('P8',1,1,20,IAV,FRMSTR,UNI,NULO,STAT)
C 
C  get name of frame to process + open it + get stand. descr.
C 
      CALL STKRDC('IN_A',1,1,120,IAV,ORFILE,UNI,NULO,STAT)
      CALL STFOPN(ORFILE,D_R4_FORMAT,0,F_IMA_TYPE,IMNO,STAT)
      CALL STDRDI(IMNO,'NAXIS',1,1,IAV,NAXPIX,UNI,NULO,STAT)
      IF (NAXPIX(1).GT.3) THEN
         CALL STTPUT('Warning: only 3 dims of frame processed ...',
     +               STAT)
         NAXPIX(1) = 3
      ENDIF
      CALL FPXWCO(0,IMNO,DD1,DD2,STAT)
      IF (STAT .GT. 0) CALL STETER(36,'Problems with WCS of frame...')
C 
      CALL STDRDI(IMNO,'NPIX',1,NAXPIX(1),IAV,NPIX,UNI,NULO,STAT)
      CALL STDRDD(IMNO,'START',1,NAXPIX(1),IAV,START,UNI,NULO,STAT)
      CALL STDRDD(IMNO,'STEP',1,NAXPIX(1),IAV,STEP,UNI,NULO,STAT)
      IAV = 1
      DO 550, N=1,3
         NAXPIX(N+1) = NPIX(N)
         IAV = IAV * NPIX(N)
550   CONTINUE
      ICOLNM(1) = IMNO                          !use ICOLNM as buffer
      CALL STFINF(ORFILE,-4,ICOLNM,STAT)        !ICOLNM(4) = real no. of pixels
      IF ((IAV .GT. ICOLNM(4)).AND.(ICOLNM(3).NE.1))
     +   CALL STETER(37,'reference file with no data...')
C 
      IF (DEFAUL(6:6).EQ.'P') THEN
         IF (INDEX(DEFAUL(3:3),'RXSM') .GT. 0) THEN
            CALL STTPUT
     +      ('No histogram plotting possible for reduced statistics...',
     +       STAT)
            DEFAUL(6:6) = 'N'
         ELSE
            CALL STDRDC(IMNO,'IDENT',1,1,72,IAV,IDENT,UNI,NULO,STAT)
            N = (NAXPIX(1)+1) * 16
            CALL STDRDC(IMNO,'CUNIT',1,1,N,IAV,CUNIT,UNI,NULO,STAT)
         ENDIF
      ENDIF
C 
      IF ((OPTION(3:3).EQ.'G').OR.(OPTION(3:3).EQ.'X')
     +    .OR. (OPTION(3:3).EQ.'W')) THEN
         EXAMED = 1
      ENDIF
C 
C *********************
C
C  test the options: full frame, subframes, intervals,
C                    ROW, COLUMN, PLANE, CURSOR or table
C 
C  IACTIO specifies the different area options:
C    = 0, full frame; 
C    = 1, table;   = 2, rows;   = 3, columns;   = 4, cursor;   = 5, planes; 
C    = 6, subframe;   = 7, covering intervals
C 
C *********************
C 
      CALL STKRDC('P2',1,1,120,IAV,OUTPUT,UNI,NULO,STAT)
      OPTION(5:5) = 'N' 
      IF ( (DEFAUL(6:6).EQ.'P') .OR. (DEFAUL(1:1).EQ.'M') ) THEN
         OPTION(5:5) = 'H'             !we also have to fill descr. HISTOGRAM
      ENDIF
      INTABL(1:) = ' '
      CALL UPCAS(OUTPUT,AREA)
C
C  complete frame
C
      IF (DEFAUL(1:1).EQ.'Y') THEN
         IACTIO = 0
         GOTO 1111
      ENDIF
C
C  area option = [...
C 
      IF (OUTPUT(1:1).EQ.'[') THEN
         N = INDEX(ORFILE,']')
         IF (N .GT. 0) THEN
            IAV = LEN(ORFILE)
            IF ((N .EQ. IAV) .OR. (ORFILE(N+1:N+1) .EQ. ' '))
     +         CALL STETER(33,'double subframe option not supported...')
         ENDIF
C
C  test, if subframe or interval option
         N = INDEX(OUTPUT,':')
         IF (N.GT.1) THEN    
            IACTIO = 6                          !for subframe(s)
         ELSE
C 
C  covering intervals 
C 
            IACTIO = 7 
            N = INDEX(AREA,'X')                    !look for 123x456
            IF (N.LT.3) THEN
               N = INDEX(AREA,'*')                 !or for 123*456
               IF (N.LT.3) 
     +            CALL STETER(2,'invalid interval option...')
            ENDIF
            AREA(N:N) = ','
            NMAL = INDEX(AREA,']')
            AREA(NMAL:NMAL) = ','
            CALL GENCNV(AREA(2:),1,2,SIZ,RVAL,DVAL,IAV)
            IF (IAV .NE. 2) 
     +         CALL STETER(34,'Invalid x/ysize of intervals...')
            IF (NPIX(3).GT.1) CALL STTPUT
     +      ('Only 1. plane of 3-dim frame will be processed...',STAT)
            NMAL = (NPIX(1)-1)/SIZ(1) + 1
            NROW = (NPIX(2)-1)/SIZ(2) + 1
         ENDIF
         GOTO 1111
      ENDIF
C 
C  rows
C 
      IF (AREA(1:3).EQ.'ROW') THEN
         IF (AREA(4:4).EQ.' ') THEN
            NROW = NPIX(2)
         ELSE IF (AREA(4:4).EQ.',') THEN
            CALL GENCNV(AREA(5:),1,25,KBUF,RVAL,DVAL,NROW)
            IF (NROW .LT. 1) CALL STETER(34,'Invalid row no.')
            DO 900,N=1,NROW
               IF ((KBUF(N).LT.1).OR.(KBUF(N).GT.NPIX(2)))
     +             CALL STETER(34,'Invalid row no...')
900         CONTINUE
         ELSE
            GOTO 999                         !table name starting with ROW..
         ENDIF
         IACTIO = 2
         IF (NPIX(3).GT.1) CALL STTPUT
     +      ('Only 1. plane of 3-dim frame will be processed...',STAT)
         GOTO 1111
      ENDIF
C 
C  planes
C 
      IF (AREA(1:5).EQ.'PLANE') THEN
         IF (AREA(6:6).EQ.' ') THEN
            NROW = NPIX(3)
         ELSE IF (AREA(6:6).EQ.',') THEN
            CALL GENCNV(AREA(7:),1,25,KBUF,RVAL,DVAL,NROW)
            IF (NROW .LT. 1) CALL STETER(34,'Invalid plane no.')
            DO 930,N=1,NROW
               IF ((KBUF(N).LT.1).OR.(KBUF(N).GT.NPIX(3)))
     +             CALL STETER(34,'Invalid plane no...')
930         CONTINUE
         ELSE
            GOTO 999                         !table name starting with PLANE..
         ENDIF
         IACTIO = 5
         IF (NPIX(3).EQ.1) THEN
            CALL STTPUT
     +      ('Single plane of 2-dim frame will be processed...',STAT)
            NROW = 1
         ENDIF
         GOTO 1111
      ENDIF
C 
C  columns
C 
      IF (AREA(1:6).EQ.'COLUMN') THEN
         IF (AREA(7:7).EQ.' ') THEN
            NROW = NPIX(1)
         ELSE IF (AREA(7:7).EQ.',') THEN
            CALL GENCNV(AREA(8:),1,25,KBUF,RVAL,DVAL,NROW)
            IF (NROW .LT. 1) CALL STETER(34,'Invalid column no.')
            DO 960,N=1,NROW
               IF ((KBUF(N).LT.1).OR.(KBUF(N).GT.NPIX(1)))
     +             CALL STETER(34,'Invalid column no...')
960         CONTINUE
         ELSE
            GOTO 999                         !table name starting with COLUMN..
         ENDIF
         IACTIO = 3
         IF (NPIX(1).EQ.1) 
     +      CALL STETER(34,'Invalid option: 1-dim input frame...')
         IF (NPIX(3).GT.1) CALL STTPUT
     +      ('Only 1. plane of 3-dim frame will be processed...',STAT)
         GOTO 1111
      ENDIF
C 
C  cursor input
C 
      IF ((AREA(1:7).EQ.'CURSOR ') .OR.
     +         (AREA(1:7).EQ.'CURSOR,')) THEN
         IACTIO = 4
         IF (AREA(7:7) .EQ. ',') THEN
            CALL GENCNV(AREA(8:),1,1,CURMAX,RVAL,DVAL,IAV)
            IF ((IAV.LT.1) .OR. (CURMAX.LT.1))
     +         CALL STETER(2,'Invalid no. for max. cursor input...')
            NROW = CURMAX
         ELSE
            CURMAX = 99999
            NROW = 100
         ENDIF
         GOTO 1111
      ENDIF
C 
C  everything else is interpreted as a table name
C 
999   IACTIO = 1                                !indicate table input
      CALL CLNTAB(OUTPUT,INTABL,0)
      CALL TBTOPN(INTABL,F_I_MODE,INTID,STAT)
C
C  get number of rows we'll work on + look for the required columns
      CALL TBIGET(INTID,N,NROW,N,N,N,STAT)
      DO 1000,N=1,4
         CALL TBLSER(INTID,INTLAB(N),ICOLNM(N),STAT)
         IF (ICOLNM(N).LE.0) THEN
            ST1 = INDEX(INTABL,' ')
            IF (ST1.LE.0) ST1 = LEN(INTABL)
            ST2 = INDEX(INTLAB(N),' ')
            IF (ST2.LE.0) ST2 = LEN(INTLAB(N))
            WRITE(OUTPUT,10030) INTLAB(N)(1:ST2),INTABL(1:ST1)
            CALL STETER(3,OUTPUT)
         ENDIF
1000  CONTINUE
C 
C  if NPIX(3) > 1, look for 3-dim table entries
      INCOL = 4
      IF (NPIX(3).GT.1) THEN
         CALL TBLSER(INTID,INTLAB(5),ICOLNM(5),STAT)
         IF (ICOLNM(5).GT.0) THEN
            CALL TBLSER(INTID,INTLAB(6),ICOLNM(6),STAT)
            IF (ICOLNM(6).GT.0) INCOL = 6
         ENDIF
      ENDIF
C 
C ***********
C
C  after the area option let's look for an output table or image
C 
C ***********
C
1111  IF ((DEFAUL(4:4).NE.'N').AND.(IACTIO.NE.4)) 
     +   CALL FRAMOU(ORFILE)
      CALL STKRDC('OUT_A',1,1,120,IAV,OUTPUT,UNI,NULO,STAT)
      TBAPP = 0                    !point to start of table
      ST2 = 0                      !default append flag to 0
      N = INDEX(OUTPUT,',')
      IF (N.GT.1) THEN
         IF ((OUTPUT(N:N+1).EQ.',i').OR.(OUTTBL(N:N+1).EQ.',I')) THEN
            CALL CLNFRA(OUTPUT(1:N-1),OUTFRA,0)   !save the file name
            OUTTBL(1:) = 'middummXyZ.tbl '
            IF (INTABL.EQ.OUTTBL) OUTTBL(1:) = 'middummXYZ.tbl '
C 
            IAV = INDEX(ORFILE,' ')
            WRITE(HISTORY,20010) ORFILE(1:IAV-1),OUTPUT
            OUTNO = 999
            OUTPUT(1:) = OUTPUT(N+1:)
C 
C  then we also need descr. IDENT + CUNIT
            CALL STDRDC(IMNO,'IDENT',1,1,72,IAV,IDENT,UNI,NULO,STAT)
            N = (NAXPIX(1)+1) * 16
            CALL STDRDC(IMNO,'CUNIT',1,1,N,IAV,CUNIT,UNI,NULO,STAT)
C 
            N = INDEX(OUTPUT,',')        !look for qualifier: name,I,qualif
            IF (N.GT.1) THEN
               CALL UPCAS(OUTPUT(N+1:),CCBUF)
               IF (CCBUF(1:1).EQ.'M') THEN
                  IF (CCBUF(2:2).EQ.'I') THEN
                     OUTNO = 1                           !minimum
                  ELSE IF (CCBUF(2:2).EQ.'A') THEN
                     OUTNO = 2                           !maximum
                  ELSE IF (CCBUF(2:2).EQ.'E') THEN
                     IF (CCBUF(3:3).EQ.'A') THEN
                        OUTNO = 3                        !mean
                     ELSE
                        OUTNO = 8                        !median
                     ENDIF
                  ELSE IF (CCBUF(2:4).EQ.'ODE') THEN
                     IF (CCBUF(5:5).EQ.'1') THEN
                        OUTNO = 9                        !1st mode
                     ELSE
                        OUTNO = 12                       !mode
                     ENDIF
                  ELSE IF (CCBUF(2:3).EQ.'OM') THEN
                     IF (CCBUF(4:4).EQ.'3') THEN
                        OUTNO = 5                        !3rd moment
                     ELSE
                        OUTNO = 6                        !4th moment
                     ENDIF
                  ENDIF
               ELSE IF (CCBUF(1:2).EQ.'ST') THEN
                  OUTNO = 4                             !stddev
               ELSE IF (CCBUF(1:2).EQ.'TO') THEN
                  OUTNO = 7                             !total intensity
               ENDIF
            ENDIF
C 
            GOTO 1490                             !skip the rest
         ENDIF
C 
         IF ((OUTPUT(N:N+1).EQ.',a').OR.(OUTTBL(N:N+1).EQ.',A')) THEN
            ST2 = 1
         ENDIF
         CALL CLNTAB(OUTPUT(1:N-1),OUTTBL,0)
C 
         CALL STECNT('GET',EC,EL,ED)
         CALL STECNT('PUT',1,0,0)            !disable errors ...
         CALL TBTOPN(OUTTBL,F_IO_MODE,OUTID,STAT)
         IF (STAT.NE.0) ST2 = 0              !clear append flag again
         CALL STECNT('PUT',EC,EL,ED)
      ELSE
         CALL CLNTAB(OUTPUT,OUTTBL,0)
      ENDIF
C
C  create NROW * NCOL1 table to have already columns allocated for later use 
1490  IF (OUTTBL(1:1).NE.'+') THEN
         IF (OPTION(3:3).EQ.'M') THEN
            NCOL1 = 2                     !only 2 columns for Minmax
         ELSE IF (OPTION(3:3).EQ.'H') THEN
            NCOL1 = 2                     !only 2 columns for Histogram
         ELSE IF (OPTION(3:3).EQ.'S') THEN
            NCOL1 = 4                     !only 4 columns for Short
         ELSE IF (OPTION(3:3).EQ.'R') THEN
            NCOL1 = 7                     !only 7 columns for Reduced
         ELSE IF (OPTION(3:3).EQ.'X') THEN
            NCOL1 = 8                     !only 8 columns for XReduced
         ELSE
            NCOL1 = 12
         ENDIF
C 
         IF (INTABL.NE.OUTTBL) THEN
            IF (ST2.EQ.0) THEN			!build new table
               CALL TBTINI(OUTTBL,0,F_O_MODE,NCOL1+4,NROW,OUTID,STAT)
               DO 1500 N=1,NCOL1                        !and define the columns
                  CALL TBCINI(OUTID,D_R4_FORMAT,1,'G12.6',TABUNT,
     +                        TABLAB(N),TCOLNM(N),STAT)
1500           CONTINUE
C 
               IF (IACTIO.EQ.4) THEN          !4 more columns for CURSOR option
                  DO 1600 N=1,4
                     IAV = N + NCOL1
                     CALL TBCINI(OUTID,D_R4_FORMAT,1,'G12.6',TABUNT,
     +                           INTLAB(N),TCOLNM(IAV),STAT)
1600              CONTINUE
               ENDIF
            ELSE                                !append to output table
C 
               CALL TBIGET(OUTID,N,TBAPP,N,N,N,STAT)      !get no. of rows
               DO 1700 N=1,NCOL1
                  CALL TBLSER(OUTID,TABLAB(N),TCOLNM(N),STAT)
                  IF (TCOLNM(N).EQ.-1) THEN
                     ST1 = INDEX(OUTTBL,' ')
                     IF (ST1.LE.0) ST1 = LEN(OUTTBL)
                     ST2 = INDEX(TABLAB(N),' ')
                     IF (ST2.LE.0) ST2 = LEN(TABLAB(N))
                     WRITE(OUTPUT,10020) TABLAB(N)(1:ST2),OUTTBL(1:ST1)
                     CALL STETER(3,OUTPUT)
                  ENDIF
1700           CONTINUE
C
               IF (IACTIO.EQ.4) THEN          !4 more columns for CURSOR option
                  DO 1770 N=1,4
                     IAV = N + NCOL1
                     CALL TBLSER(OUTID,INTLAB(N),TCOLNM(IAV),STAT)
                     IF (TCOLNM(IAV).EQ.-1) THEN
                        ST1 = INDEX(OUTTBL,' ')
                        IF (ST1.LE.0) ST1 = LEN(OUTTBL)
                        ST2 = INDEX(INTLAB(N),' ')
                        IF (ST2.LE.0) ST2 = LEN(INTLAB(N))
                        WRITE(OUTPUT,10020) 
     +                       INTLAB(N)(1:ST2),OUTTBL(1:ST1)
                        CALL STETER(3,OUTPUT)
                     ENDIF
1770              CONTINUE
               ENDIF
            ENDIF
C
C  we use same table for input + output
         ELSE
            OUTID = INTID
            DO 1800, N=1,NCOL1
               CALL TBLSER(OUTID,TABLAB(N),TCOLNM(N),STAT)
               IF (TCOLNM(N).EQ.-1)
     +            CALL TBCINI(OUTID,D_R4_FORMAT,1,'G12.6',TABUNT,
     +                        TABLAB(N),TCOLNM(N),STAT)
1800        CONTINUE
         ENDIF
      ENDIF
C 
C ******************
C 
C  now do the statistics for every area option
C 
C ******************
C
C ---------
C
C  1) process complete frame 
C
C ---------
C
      IF (IACTIO.EQ.0) THEN
         OPTION(5:5) = 'Y'
         IF (EXAMED.EQ.1) THEN
            CALL EXMED(IACTIO,IMNO,ZBINS(2),NPIX,
     +                 SPIX,EPIX,XMEDIA,STAT)
            IF (STAT.NE.0) CALL STETER(35,
     +         'Problems with exact median calculation...')
         ENDIF
         CALL ZMSTAT(IMNO,AREA,NAXPIX,START,STEP,
     +               ZBINS,FRMSTR,OPTION)
         IF (OUTID.NE.-999) THEN
            N = 1
            CALL STKRDR('OUTPUTR',1,12,IAV,RVAL,UNI,NULO,STAT)
            CALL TBRWRR(OUTID,N+TBAPP,NCOL1,TCOLNM,RVAL,STAT)
         ENDIF
         IF (DEFAUL(6:6).EQ.'P') CALL PLOTI(IMNO)
C 
C ---------
C
C  2) handle subframe(s) mode
C 
C ---------
C
      ELSE IF (IACTIO.EQ.6) THEN
         NROW = 1
         OUTPUT(1:) = AREA(1:) 
700      N = INDEX(OUTPUT,'],[')
         IF (N.GT.1) THEN
            NROW = NROW + 1
            OUTPUT(1:) = OUTPUT(N+2:)
            GOTO 700
         ENDIF
         IF (NROW.EQ.1) OPTION(5:5) = 'Y'        !that's how it was before...
C
C  if Merge option we first get the no. of pixels involved + then create
C  an artificial image with all those pixels
         IF (DEFAUL(1:1).EQ.'M') THEN
            HISTORY(1:100) = AREA(1:)            !save AREA string
            SIZ(1) = 0
            DO 770, NLOOP=1,NROW                 !loop through subframes
               N = INDEX(AREA(1:),']')
               OUTPUT(1:) = AREA(1:N)//' '
               CALL EXTCOO(IMNO,OUTPUT,3,IAV,SPIX,EPIX,STAT)
               SIZ(2) = (EPIX(3) - SPIX(3) + 1) *
     +                  (EPIX(2) - SPIX(2) + 1) *
     +                  (EPIX(1) - SPIX(1) + 1)
               SIZ(1) = SIZ(1) + SIZ(2)
               AREA(1:) = AREA(N+2:)//' '
770         CONTINUE
            CALL ARTIMA(0,IMNO,NPIX,AREA,SIZ(1),IMNOC,STAT)
            AREA(1:) = HISTORY(1:100)
            DO 800, NLOOP=1,NROW
               N = INDEX(AREA(1:),']')
               OUTPUT(1:) = AREA(1:N)//' '
               CALL ARTIMA(1,IMNO,NPIX,OUTPUT,0,IMNOC,STAT)
               CALL STTPUT(OUTPUT,STAT)
               AREA(1:) = AREA(N+2:)//' '
800         CONTINUE
            GOTO 8000                          !now proceed as if full frame
         ENDIF
C
C  loop through all subframes
         DO 850, NLOOP=1,NROW
            N = INDEX(AREA(1:),']')
            OUTPUT(1:) = AREA(1:N)//' '
            IF (EXAMED.EQ.1) THEN
               CALL EXTCOO(IMNO,OUTPUT,3,IAV,SPIX,EPIX,STAT)
               IF (STAT.NE.0) CALL STETER(34,
     +              'Invalid syntax in subframe ...')
               IF (IAV.LT.3) THEN
                  CALL EXMED(IACTIO,IMNO,ZBINS(2),NPIX,SPIX,
     +                       EPIX,XMEDIA,STAT)
                  IF (STAT.NE.0) CALL STETER(35,
     +               'Problems with exact median calculation...')
               ELSE
                  CALL STETER(27,
     +            'G- or X-option not supported for 3-d frames, yet')
               ENDIF
            ENDIF
            CALL ZMSTAT(IMNO,OUTPUT,NAXPIX,START,STEP,
     +                  ZBINS,FRMSTR,OPTION)
           IF (OUTID.NE.-999) THEN
              CALL STKRDR('OUTPUTR',1,12,IAV,RVAL,UNI,NULO,STAT)
              CALL TBRWRR(OUTID,NLOOP+TBAPP,NCOL1,TCOLNM,RVAL,STAT)
            ENDIF
            IF (DEFAUL(6:6).EQ.'P') CALL PLOTI(IMNO)
            AREA(1:) = AREA(N+2:)//' '
850      CONTINUE
C 
C ---------
C
C  3) handle interval mode
C 
C ---------
C
      ELSE IF (IACTIO.EQ.7) THEN
         N = 1
         SPIX(2) = 1
C 
         IF (NAXPIX(1).GE.3) THEN
            NPIX(3) = 1                         !work on 1st plane only
C 
            DO 2200, ST2=1,NROW                   !loop over all intervals
               SPIX(1) = 1
               EPIX(2) = SPIX(2)+SIZ(2)-1
               IF (EPIX(2).GT.NPIX(2)) EPIX(2) = NPIX(2)
C 
               DO 2100, ST1=1,NMAL
                  EPIX(1) = SPIX(1)+SIZ(1)-1
                  IF (EPIX(1).GT.NPIX(1)) EPIX(1) = NPIX(1)
                  IF (EXAMED.EQ.1) THEN
                     CALL EXMED(IACTIO,IMNO,ZBINS(2),NPIX,SPIX,EPIX,
     +                          XMEDIA,STAT)
                     IF (STAT.NE.0) CALL STETER(35,
     +                  'Problems with exact median calculation...')
                  ENDIF
                  WRITE(AREA,10009) SPIX(1),SPIX(2),EPIX(1),EPIX(2)
                  CALL ZMSTAT(IMNO,AREA,NAXPIX,START,
     +                        STEP,ZBINS,FRMSTR,OPTION)
                  IF (OUTID.NE.-999) THEN
                     CALL STKRDR('OUTPUTR',1,12,IAV,RVAL,UNI,NULO,STAT)
                     CALL TBRWRR(OUTID,N+TBAPP,NCOL1,TCOLNM,RVAL,STAT)
                     N = N + 1
                     ENDIF
                  IF (DEFAUL(6:6).EQ.'P') CALL PLOTI(IMNO)
                  SPIX(1) = SPIX(1) + SIZ(1)
2100           CONTINUE
               SPIX(2) = SPIX(2) + SIZ(2)
2200        CONTINUE
C 
         ELSE IF (NAXPIX(1).EQ.2) THEN
            DO 2400, ST2=1,NROW                   !loop over all intervals
               SPIX(1) = 1
               EPIX(2) = SPIX(2)+SIZ(2)-1
               IF (EPIX(2).GT.NPIX(2)) EPIX(2) = NPIX(2)
C 
               DO 2300, ST1=1,NMAL
                  EPIX(1) = SPIX(1)+SIZ(1)-1
                  IF (EPIX(1).GT.NPIX(1)) EPIX(1) = NPIX(1)
                  IF (EXAMED.EQ.1) THEN
                     CALL EXMED(IACTIO,IMNO,ZBINS(2),NPIX,SPIX,EPIX,
     +                          XMEDIA,STAT)
                     IF (STAT.NE.0) CALL STETER(35,
     +                  'Problems with exact median calculation...')
                  ENDIF
                  WRITE(AREA,10008) SPIX(1),SPIX(2),EPIX(1),EPIX(2)
                  CALL ZMSTAT(IMNO,AREA,NAXPIX,START,
     +                        STEP,ZBINS,FRMSTR,OPTION)
                  IF (OUTID.NE.-999) THEN
                     CALL STKRDR('OUTPUTR',1,12,IAV,RVAL,UNI,NULO,STAT)
                     CALL TBRWRR(OUTID,N+TBAPP,NCOL1,TCOLNM,RVAL,STAT)
                     N = N + 1
                  ENDIF
                  IF (DEFAUL(6:6).EQ.'P') CALL PLOTI(IMNO)
                  SPIX(1) = SPIX(1) + SIZ(1)
2300           CONTINUE
               SPIX(2) = SPIX(2) + SIZ(2)
2400        CONTINUE
C 
         ELSE
            SPIX(1) = 1
            SPIX(2) = 1
            EPIX(2) = 1
            DO 2500, ST2=1,NMAL                   !loop over all intervals
               EPIX(1) = SPIX(1)+SIZ(1)-1
               IF (EPIX(1).GT.NPIX(1)) EPIX(1) = NPIX(1)
               IF (EXAMED.EQ.1) THEN
                  CALL EXMED(IACTIO,IMNO,ZBINS(2),NPIX,SPIX,EPIX,
     +                       XMEDIA,STAT)
                  IF (STAT.NE.0) CALL STETER(35,
     +               'Problems with exact median calculation...')
               ENDIF
               WRITE(AREA,10007) SPIX(1),EPIX(1)
               CALL ZMSTAT(IMNO,AREA,NAXPIX,START,
     +                     STEP,ZBINS,FRMSTR,OPTION)
               IF (OUTID.NE.-999) THEN
                  CALL STKRDR('OUTPUTR',1,12,IAV,RVAL,UNI,NULO,STAT)
                  CALL TBRWRR(OUTID,N+TBAPP,NCOL1,TCOLNM,RVAL,STAT)
                  N = N + 1
               ENDIF
               IF (DEFAUL(6:6).EQ.'P') CALL PLOTI(IMNO)
               SPIX(1) = SPIX(1) + SIZ(1)
2500        CONTINUE
         ENDIF
         N = NMAL * NROW
         WRITE(OUTPUT,10011) N             !display no. of intervals processed
         CALL STTPUT(OUTPUT,STAT)
C 
C ---------
C
C  4) handle row modes
C 
C ---------
C
      ELSE IF (IACTIO.EQ.2) THEN
         CACT(1:) = 'RRRR    '
C 
C  if Merge option we first get the no. of pixels involved + then create 
C  an artificial image with all those pixels
         IF (DEFAUL(1:1).EQ.'M') THEN
            SIZ(1) = NROW * NPIX(1)
            CALL ARTIMA(0,IMNO,NPIX,AREA,SIZ(1),IMNOC,STAT)
            DO 2700, NLOOP=1,NROW
               IF (KBUF(1).GT.0) THEN
                  N = KBUF(NLOOP)
               ELSE
                  N = NLOOP 
               ENDIF
               IF (NPIX(3).GT.1) THEN
                  WRITE(AREA,10004) N,N
               ELSE
                  WRITE(AREA,10003) N,N
               ENDIF
               CALL ARTIMA(1,IMNO,NPIX,AREA,0,IMNOC,STAT)
               WRITE(OUTPUT,20101) N
               CALL STTPUT(OUTPUT,STAT)
2700        CONTINUE
            GOTO 8000                           !now proceed as if full frame
         ENDIF
C
C  loop through all rows
         DO 2800, NLOOP=1,NROW
            IF (KBUF(1).GT.0) THEN
               N = KBUF(NLOOP)
            ELSE
               N = NLOOP 
            ENDIF
            WRITE(AREA,10000) CACT(1:1),N,N
            IF (EXAMED.EQ.1) THEN
               SPIX(1) = 1
               SPIX(2) = N
               EPIX(1) = NPIX(1)
               EPIX(2) = N
               CALL EXMED(IACTIO,IMNO,ZBINS(2),NPIX,SPIX,EPIX,
     +                    XMEDIA,STAT)
               IF (STAT.NE.0) CALL STETER(35,
     +            'Problems with exact median calculation...')
            ENDIF
            CALL ZMSTAT(IMNO,AREA,NAXPIX,START,STEP,
     +                  ZBINS,FRMSTR,OPTION)
            IF (OUTID.NE.-999) THEN
               CALL STKRDR('OUTPUTR',1,12,IAV,RVAL,UNI,NULO,STAT) 
               CALL TBRWRR(OUTID,NLOOP+TBAPP,NCOL1,TCOLNM,RVAL,STAT)
            ENDIF
            IF (DEFAUL(6:6).EQ.'P') CALL PLOTI(IMNO)
2800     CONTINUE
C 
C 
C ---------
C
C  5) handle column modes
C 
C ---------
C
      ELSE IF (IACTIO.EQ.3) THEN
         CACT(1:) = 'CCCC    '
C
C  transpose frame + switch descriptors
         SIZ(1) = 1
         DO 2900, N=1,NAXPIX(1)
            SIZ(1) = SIZ(1) * NPIX(N)
2900     CONTINUE
         CALL STFMAP(IMNO,F_I_MODE,1,SIZ(1),IAV,PNTRA,STAT)
         CALL STFCRE(TMPFIL,D_R4_FORMAT,F_O_MODE,F_IMA_TYPE,
     +               SIZ(1),IMNOB,STAT)
         CALL STFMAP(IMNOB,F_O_MODE,1,SIZ(1),IAV,PNTRB,STAT)
         SIZ(1) = 128
         SIZ(2) = 256
         CALL LINCOL(MADRID(PNTRA),NPIX,SIZ,MADRID(PNTRB))
         CALL STFUNM(IMNO,STAT)
         CALL STFUNM(IMNOB,STAT)
C
C  we have to write these descriptors to the transposed frame,
C  because of function Convcoo inside Zstats (ZMSTAT)
C 
         CALL STDWRI(IMNOB,'NAXIS',NAXPIX(1),1,1,UNI,STAT)
         N = NAXPIX(2)
         NAXPIX(2) = NAXPIX(3)
         NAXPIX(3) = N
         CALL STDWRI(IMNOB,'NPIX',NAXPIX(2),1,2,UNI,STAT)
         DVAL = START(1)
         START(1) = START(2)
         START(2) = DVAL
         CALL STDWRD(IMNOB,'STEP',STEP,1,2,UNI,STAT)
         DVAL = STEP(1)
         STEP(1) = STEP(2)
         STEP(2) = DVAL
         CALL STDWRD(IMNOB,'START',START,1,2,UNI,STAT)
         NPIX(1) = NAXPIX(2)                 !NAXPIX already swapped
         NPIX(2) = NAXPIX(3)
         OUTPUT(1:) = CUNIT(16:)             !switch CUNIT (needed for Wcoords)
         CUNIT(16:31) = OUTPUT(32:47)
         CUNIT(32:47) = OUTPUT(16:31)
         N = (NAXPIX(1)+1) * 16
         CALL STDWRC(IMNOB,'CUNIT',1,CUNIT,1,N,UNI,STAT)
C
C  if Merge option we first get the no. of pixels involved + then create
C  an artificial image with all those pixels
         IF (DEFAUL(1:1).EQ.'M') THEN
            SIZ(1) = NROW * NPIX(1)
            CALL ARTIMA(0,IMNOB,NPIX,AREA,SIZ(1),IMNOC,STAT)
            DO 2950, NLOOP=1,NROW
               IF (KBUF(1).GT.0) THEN
                  N = KBUF(NLOOP)
               ELSE
                  N = NLOOP
               ENDIF
               IF (NPIX(3).GT.1) THEN
                  WRITE(AREA,10004) N,N
               ELSE
                  WRITE(AREA,10003) N,N
               ENDIF
               CALL ARTIMA(1,IMNOB,NPIX,AREA,0,IMNOC,STAT)
               WRITE(OUTPUT,20102) N
               CALL STTPUT(OUTPUT,STAT)
2950        CONTINUE
            GOTO 8000                           !now proceed as if full frame
         ENDIF
C
C  now loop through rows (columns of original frame)
         DO 3000, NLOOP=1,NROW
            IF (KBUF(1).GT.0) THEN
               N = KBUF(NLOOP)
            ELSE
               N = NLOOP
            ENDIF
            WRITE(AREA,10000) CACT(1:1),N,N
            IF (EXAMED.EQ.1) THEN
               SPIX(1) = 1
               SPIX(2) = N
               EPIX(1) = NAXPIX(2)                  !x-dim
               EPIX(2) = N
               CALL EXMED(IACTIO,IMNOB,ZBINS(2),NPIX,SPIX,EPIX,
     +                    XMEDIA,STAT)
               IF (STAT.NE.0) CALL STETER(35,
     +            'Problems with exact median calculation...')
            ENDIF
            CALL ZMSTAT(IMNOB,AREA,NAXPIX,START,STEP,
     +                  ZBINS,FRMSTR,OPTION)
            IF (OUTID.NE.-999) THEN
               CALL STKRDR('OUTPUTR',1,12,IAV,RVAL,UNI,NULO,STAT) 
               CALL TBRWRR(OUTID,NLOOP+TBAPP,NCOL1,TCOLNM,RVAL,STAT)
            ENDIF
            IF (DEFAUL(6:6).EQ.'P') CALL PLOTI(IMNO)
3000     CONTINUE
         CALL STFDEL(TMPFIL,STAT)
C 
C ---------
C
C  6) in table mode, we need some more table massaging
C 
C ---------
C
      ELSE IF (IACTIO.EQ.1) THEN
C 
C  if Merge option we first have to get the no. of pixels involved
         IF (DEFAUL(1:1).EQ.'M') THEN
            SIZ(1) = 0
            DO 4100, N=1,NROW
               CALL TBSGET(INTID,N,SELFLG,STAT)        !only use selected rows
               IF (SELFLG) THEN
                  CALL TBRRDR(INTID,N,INCOL,ICOLNM(1),RVAL,TABNUL,STAT)
                  IF (INCOL.EQ.4) THEN
                     WRITE(AREA,10001) RVAL(1),RVAL(2),RVAL(3),RVAL(4)
                  ELSE
                     WRITE(AREA,10002) 
     +               RVAL(1),RVAL(2),RVAL(5),RVAL(3),RVAL(4),RVAL(6)
                  ENDIF
                  CALL EXTCOO(IMNO,AREA,3,IAV,SPIX,EPIX,STAT)
                  SIZ(2) = (EPIX(3) - SPIX(3) + 1) * 
     +                     (EPIX(2) - SPIX(2) + 1) *
     +                     (EPIX(1) - SPIX(1) + 1)
                  SIZ(1) = SIZ(1) + SIZ(2)
               ENDIF
4100        CONTINUE
            CALL ARTIMA(0,IMNO,NPIX,AREA,SIZ(1),IMNOC,STAT)
C 
C  then create artificial image with all those pixels
            DO 4200, N=1,NROW
               CALL TBSGET(INTID,N,SELFLG,STAT)        !only use selected rows
               IF (SELFLG) THEN
                  CALL TBRRDR(INTID,N,INCOL,ICOLNM(1),RVAL,TABNUL,STAT)
                  IF (INCOL.EQ.4) THEN
                     WRITE(AREA,10001) RVAL(1),RVAL(2),RVAL(3),RVAL(4)
                  ELSE
                     WRITE(AREA,10002) 
     +               RVAL(1),RVAL(2),RVAL(5),RVAL(3),RVAL(4),RVAL(6)
                  ENDIF
                  CALL BLANKO(AREA)
                  CALL STTPUT(AREA,STAT)
                  CALL ARTIMA(1,IMNO,NPIX,AREA,0,IMNOC,STAT)
               ENDIF
4200        CONTINUE
            GOTO 8000                           !now proceed as if full frame
         ENDIF
C 
C  now loop through input table rows
         NMAL = 0
         DO 4400, N=1,NROW
C 
            CALL TBSGET(INTID,N,SELFLG,STAT)        !only use selected rows
            IF (.NOT.SELFLG) GOTO 4400
C 
            NMAL = NMAL + 1
            AREA(1:) = ' '
            CALL TBRRDR(INTID,N,INCOL,ICOLNM(1),RVAL,TABNUL,STAT)
C
            IF (INCOL.EQ.4) THEN
               WRITE(AREA,10001) RVAL(1),RVAL(2),RVAL(3),RVAL(4)
            ELSE
               WRITE(AREA,10002) 
     +              RVAL(1),RVAL(2),RVAL(5),RVAL(3),RVAL(4),RVAL(6)
            ENDIF
            CALL BLANKO(AREA)
            IF (OPTION(4:4) .NE. 'N') THEN
               OUTPUT(1:) = ' '
               CALL STTPUT(OUTPUT,STAT)
            ENDIF
            IF (EXAMED.EQ.1) THEN
               IF (INCOL.EQ.4) THEN
                  DD1(1) = RVAL(1)                !convert wc -> fp
                  DD1(2) = RVAL(2)
                  CALL FPXWCO(-1,0,DD1,DD2,STAT)
                  IF (STAT .GT. 0) 
     +               CALL STETER(36,'Problems with WCS of frame...')
                  SPIX(1) = DD2(1)
                  SPIX(2) = DD2(2)
                  DD1(1) = RVAL(3)
                  DD1(2) = RVAL(4)
                  CALL FPXWCO(-1,0,DD1,DD2,STAT)
                  IF (STAT .GT. 0) 
     +               CALL STETER(36,'Problems with WCS of frame...')
                  EPIX(1) = DD2(1)
                  EPIX(2) = DD2(2)
                  CALL EXMED(IACTIO,IMNO,ZBINS(2),NPIX,SPIX,EPIX,
     +                       XMEDIA,STAT)
                  IF (STAT.NE.0) CALL STETER(35,
     +               'Problems with exact median calculation...')
               ELSE
                  CALL STETER(27,
     +            'G- or X-option not supported for 3-d frames, yet')
               ENDIF
            ENDIF
            CALL ZMSTAT(IMNO,AREA,NAXPIX,START,STEP,
     +                  ZBINS,FRMSTR,OPTION)
            IF (OUTID.NE.-999) THEN
               CALL STKRDR('OUTPUTR',1,12,IAV,RVAL,UNI,NULO,STAT) 
               CALL TBRWRR(OUTID,N+TBAPP,NCOL1,TCOLNM,RVAL,STAT)
            ENDIF
            IF (DEFAUL(6:6).EQ.'P') CALL PLOTI(IMNO)
4400     CONTINUE
C
         CALL TBTCLO(INTID,STAT)
         WRITE(OUTPUT,10010) NMAL             !display no. of rows processed
         CALL STTPUT(OUTPUT,STAT)
C
C ---------
C
C  7) handle plane modes
C
C ---------
C
      ELSE IF (IACTIO.EQ.5) THEN
         CACT(1:) = 'PPPP    '
C
C  if Merge option we first get the no. of pixels involved + then create
C  an artificial image with all those pixels
         IF (DEFAUL(1:1).EQ.'M') THEN
            SIZ(1) = NROW * NPIX(1)*NPIX(2)
            CALL ARTIMA(0,IMNO,NPIX,AREA,SIZ(1),IMNOC,STAT)
            DO 4700, NLOOP=1,NROW
               IF (KBUF(1).GT.0) THEN
                  N = KBUF(NLOOP)
               ELSE
                  N = NLOOP
               ENDIF
               WRITE(AREA,10005) N,N
               CALL ARTIMA(1,IMNO,NPIX,AREA,0,IMNOC,STAT)
               WRITE(OUTPUT,20100) N
               CALL STTPUT(OUTPUT,STAT)
4700        CONTINUE
            GOTO 8000                           !now proceed as if full frame
         ENDIF
C
C  loop through planes
         DO 5000, NLOOP=1,NROW
            IF (KBUF(1).GT.0) THEN
               N = KBUF(NLOOP)
            ELSE
               N = NLOOP
            ENDIF
            AREA(1:) = ' '
            WRITE(AREA,20000) CACT(1:1),N,N
            IF (EXAMED.EQ.1) THEN
               SPIX(1) = 1
               SPIX(2) = 1
               SPIX(3) = N
               EPIX(1) = NPIX(1)
               EPIX(2) = NPIX(2)
               EPIX(3) = N
               CALL EXMED(IACTIO,IMNO,ZBINS(2),NPIX,SPIX,EPIX,
     +                    XMEDIA,STAT)
               IF (STAT.NE.0) CALL STETER(35,
     +            'Problems with exact median calculation...')
            ENDIF
            CALL ZMSTAT(IMNO,AREA,NAXPIX,START,STEP,
     +                  ZBINS,FRMSTR,OPTION)
            IF (OUTID.NE.-999) THEN
               CALL STKRDR('OUTPUTR',1,12,IAV,RVAL,UNI,NULO,STAT)
               CALL TBRWRR(OUTID,NLOOP+TBAPP,NCOL1,TCOLNM,RVAL,STAT)
            ENDIF
            IF (DEFAUL(6:6).EQ.'P') CALL PLOTI(IMNO)
5000     CONTINUE
C
C ---------
C
C  8) finally process cursor mode
C 
C ---------
C
      ELSE
         CALL DTOPEN(1,STAT)
         CALL SETCUR(QDSPNO,2,0,0,COOS,STAT)             !set up cursor window
         CACT = 'NNYY?C0 '
         IF (DEFAUL(5:5).NE.'N') CACT(2:2) = 'Y'
         ORFILE(1:) = ' '
         COOFF = 0
C
C  if Merge option we create the art. image with 100 000 pixels first
C  we may have to extend the image later on...
         IF (DEFAUL(1:1).EQ.'M') THEN
            SIZ(1) = 100000
            CALL ARTIMA(0,IMNO,NPIX,AREA,SIZ(1),IMNOC,STAT)
         ENDIF
C 
C  now loop through cursor inputs
         DO 6000, N=1,CURMAX
5800        CALL GETCUR(CACT,ORFILE,
     +                  XY1,RVAL(1),RVAL(13),RVAL(10),ST1,
     +                  XY2,RVAL(3),RVAL(15),RVAL(11),ST2)
C 
            IF ((ST1.EQ.0).AND.(ST2.EQ.0)) THEN
               IF ((N.EQ.1).AND.(COOFF.EQ.0)) THEN
                  COOFF = 1
                  ORFILE(1:) = ' '
                  CALL STTDIS
     +            ('switch cursors on - next time we exit...',99,STAT)
                  GOTO 5800
               ELSE
                  GOTO 6600         !we're finished
               ENDIF
            ENDIF
C 
            AREA(1:) = ' '
            IF (NAXPIX(1) .GT. 2) THEN
               WRITE(AREA,11001) 
     +         RVAL(13),RVAL(14),ZPLANE,RVAL(15),RVAL(16),ZPLANE
            ELSE
               WRITE(AREA,10001) RVAL(13),RVAL(14),RVAL(15),RVAL(16)
            ENDIF
            CALL BLANKO(AREA)                  !get rid of blanks in string
            IF (DEFAUL(1:1).EQ.'M') THEN
               CALL STTPUT(AREA,STAT)
5880           CALL ARTIMA(1,IMNO,NPIX,AREA,0,IMNOC,STAT)
               IF (STAT .EQ. -4) THEN                 !expand art. image
                  SIZ(1) = SIZ(1) + SIZ(1)
                  CALL ARTIMA(2,IMNO,NPIX,OUTPUT,SIZ(1),IMNOC,STAT)
                  GOTO 5880                           !redo last filling
               ENDIF
            ELSE
               IF (EXAMED.EQ.1) THEN
                  SPIX(1) = NINT(RVAL(1))
                  SPIX(2) = NINT(RVAL(2))
                  EPIX(1) = NINT(RVAL(3))
                  EPIX(2) = NINT(RVAL(4))
                  CALL EXMED(IACTIO,IMNO,ZBINS(2),NPIX,
     +                       SPIX,EPIX,XMEDIA,STAT)
                  IF (STAT.NE.0) CALL STETER(35,
     +               'Problems with exact median calculation...')
               ENDIF
               CALL ZMSTAT(IMNO,AREA,NAXPIX,START,STEP,
     +                     ZBINS,FRMSTR,OPTION)
               IF (OUTID.NE.-999) THEN
                  CALL STKRDR('OUTPUTR',1,12,IAV,RVAL,UNI,NULO,STAT) 
                  NROW = N+TBAPP
                  CALL TBRWRR(OUTID,NROW,NCOL1,TCOLNM,RVAL,STAT)
                  IAV = NCOL1 + 1           !first column after statistic stuff
                  CALL TBRWRR(OUTID,NROW,4,TCOLNM(IAV),RVAL(13),STAT)
               ENDIF
               IF (DEFAUL(6:6).EQ.'P') CALL PLOTI(IMNO)
            ENDIF
6000     CONTINUE
C 
6600     CALL DTCLOS(QDSPNO)
         IF (DEFAUL(1:1).EQ.'M') THEN                 !get exact size
            CALL ARTIMA(3,IMNO,NPIX,AREA,SIZ(1),IMNOC,STAT)
            GOTO 8000                                 !do the statistics now 
         ENDIF
C 
      ENDIF
      GOTO 9000
C
C ---------
C
C  here from Merge option - proceed as if full frame
C
C ---------
C
8000  NPIX(1) = SIZ(1)
      NPIX(2) = 1
      NPIX(3) = 1
      IF (EXAMED.EQ.1) THEN
         CALL EXMED(0,IMNOC,ZBINS(2),NPIX,SPIX,EPIX,XMEDIA,STAT)
         IF (STAT.NE.0) 
     +      CALL STETER(35,'Problems with exact median calculation...')
      ENDIF
      NAXPIX(1) = 1
      NAXPIX(2) = NPIX(1)
      NAXPIX(3) = 1
      NAXPIX(4) = 1
      OPTION(5:5) = 'N'                   !don't try to write descr.
      START(1) = 0.D0
      STEP(1) = 1.D0
      CALL ZMSTAT(IMNOC,AREA,NAXPIX,START,STEP,ZBINS,FRMSTR,OPTION)
      IF (OUTID.NE.-999) THEN
         N = 1
         CALL STKRDR('OUTPUTR',1,12,IAV,RVAL,UNI,NULO,STAT)
         CALL TBRWRR(OUTID,N+TBAPP,NCOL1,TCOLNM,RVAL,STAT)
      ENDIF
      IF (DEFAUL(6:6).EQ.'P') CALL PLOTI(IMNO)
C 
C 
C ***********************
C
C  everything joins here
C
C ***********************
C 
C
9000  IF ((OUTID.NE.-999) .AND. (OUTID.NE.INTID)) THEN
         IF (OUTNO .GT. 0) THEN               !copy temp table to outimage
            CALL TBIGET(OUTID,NPXO(1),NPXO(2),N,N,N,STAT)
            STRTO(1) = 0.0
            STPO(1) = 1.0
            IF (OUTNO.EQ.999) THEN
               STRTO(2) = 0.0
               STPO(2) = 1.0
               NAXPIX(1) = 2                       !2dim result image
               N = 1
            ELSE
               NAXPIX(1) = 1                       !1dim result image
               N = 2
            ENDIF
            CALL STIPUT(OUTFRA,D_R4_FORMAT,F_O_MODE,F_IMA_TYPE,
     +            NAXPIX(1),NPXO(N),STRTO,STPO,
     +            IDENT,CUNIT,PNTRO,IMNOB,STAT)
            CALL TB2IMA(OUTNO,OUTID,NPXO,MADRID(PNTRO))   !table -> image
            CALL STDWRC(IMNOB,'HISTORY',1,HISTORY,1,80,UNI,STAT)
            IF (IACTIO.EQ.7) THEN                  !save xsize, ysize of intv.
               ICOLNM(1) = SIZ(1)
               ICOLNM(2) = SIZ(2)
               ICOLNM(3) = NMAL
               ICOLNM(4) = NROW
               CALL STDWRI(IMNOB,'INTVALS',ICOLNM,1,4,UNI,STAT)
               CALL STECNT('GET',EC,EL,ED)
               CALL STECNT('PUT',1,0,0)            !disable errors ...
               CALL STKRDC('MYCOMND',1,1,80,IAV,HISTORY,UNI,NULO,STAT)
               IF (STAT.EQ.0) 
     +            CALL STDWRC(IMNOB,'HISTORY',1,HISTORY,1,80,UNI,STAT)
               CALL STECNT('PUT',EC,EL,ED)
            ENDIF
            CALL STFCLO(IMNOB,STAT)
         ELSE
            CALL TBSINI(OUTID,STAT)               !init selection flags
            CALL TBTCLO(OUTID,STAT)               !close table correctly
         ENDIF
      ENDIF
      CALL STSEPI
C
10000 FORMAT(A,'[<,@',I5.5,':>,@',I5.5,'] ')
10001 FORMAT('[',G12.6,',',G12.6,':',G12.6,',',G12.6,'] ')
10002 FORMAT
     +('[',G12.6,',',G12.6,',',G12.6,':',G12.6,',',G12.6,',',G12.6,']')
10003 FORMAT('[<,@',I5.5,':>,@',I5.5,'] ')
10004 FORMAT('[<,@',I5.5,',<:>,@',I5.5,',<] ')
10005 FORMAT('[<,<,@',I5.5,':>,>,@',I5.5,'] ')
10007 FORMAT('[@',I7.7,':@',I7.7,'] ')
10008 FORMAT('[@',I7.7,',@',I7.7,':@',I7.7,',@',I7.7,'] ')
10009 FORMAT('[@',I7.7,',@',I7.7,',<:@',I7.7,'@',I7.7,',<] ')
10010 FORMAT(I7,' table entries processed.')
10011 FORMAT(I7,' intervals processed.')
10020 FORMAT('column :',A,'is missing in output table ',A)
10030 FORMAT('column :',A,'is missing in input table ',A)
11001 FORMAT
     +('[',G12.6,',',G12.6,',@',I4,':',G12.6,',',G12.6,',@',I4,']')
20000 FORMAT(A,'[<,<,@',I5.5,':>,>,@',I5.5,'] ')
20010 FORMAT('from STATIST/IMAGE ',A,' ... ',A)
20100 FORMAT('plane ',I4)
20101 FORMAT('row ',I6)
20102 FORMAT('column ',I6)
C
      END

C 
C  copy all or just on column of the result table
C  to a 2/1-dim image
C 
      SUBROUTINE TB2IMA(OUTNO,TID,NCRO,A)
C 
      IMPLICIT NONE
C 
      INTEGER      OUTNO,TID,NCRO(*)
      INTEGER      INDX,N,NCOL,NROW,STAT
      INTEGER      ICOLNM(12)                !max no. of columns
      INTEGER      TABNUL(12)
C 
      REAL         A(*),RBUF(12)
C 
      NCOL = NCRO(1)
      NROW = NCRO(2)
      DO 500, N=1,NCOL
         ICOLNM(N) = N
500   CONTINUE
C 
      IF (OUTNO.EQ.999) THEN
         INDX = 1
         DO 1000, N=1,NROW
            CALL TBRRDR(TID,N,NCOL,ICOLNM(1),A(INDX),TABNUL,STAT)
            INDX = INDX + NCOL
1000     CONTINUE
C 
      ELSE
         DO 1500, N=1,12
            RBUF(N) = 0.0
1500     CONTINUE
         INDX = 1
         DO 2000, N=1,NROW
            CALL TBRRDR(TID,N,NCOL,ICOLNM(1),RBUF,TABNUL,STAT)
            A(INDX) = RBUF(OUTNO)
            INDX = INDX + 1
2000     CONTINUE
      ENDIF
C 
      RETURN
      END

      SUBROUTINE EXMED(ACTIO,INO,ZBINS,NPIX,SPIX,EPIX,VALMED,STAT)
C 
      INTEGER      ACTIO,INO,NPIX(*),SPIX(*),EPIX(*),STAT
      INTEGER      IAV,N,CHUNKA,CHUNKB,NN
      INTEGER      NA,NB,INSIZ,OUTSIZ,UNI(1),IBUF(7)
      INTEGER*8    PNTRX,PNTRY
      INTEGER      IMNOX,IMNOY,MADRID(1)
C 
      REAL         ZBINS(*),VALMED
      REAL         RBUF(7)
C
      INCLUDE 'MID_INCLUDE:ST_DEF.INC'
      COMMON  /VMR/ MADRID
      INCLUDE 'MID_INCLUDE:ST_DAT.INC'
C 
C  plane option
C 
      IF (ACTIO.EQ.5) THEN       
         INSIZ = NPIX(1) * NPIX(2)
         OUTSIZ = INSIZ
         NA = ((SPIX(3)-1)*NPIX(1)*NPIX(2)) + 1
         CALL STFCRE('middstat2',D_R4_FORMAT,F_X_MODE,F_IMA_TYPE,
     +               OUTSIZ,IMNOY,STAT)
         CALL STFMAP(IMNOY,F_X_MODE,1,OUTSIZ,IAV,PNTRY,STAT)
         CALL STFGET(INO,NA,INSIZ,IAV,MADRID(PNTRY),STAT)
C 
C  subframe input
C 
      ELSE IF (ACTIO.NE.0) THEN 
         NN = EPIX(2)-SPIX(2)+1
         INSIZ = NPIX(1) * NN                 !size of full Y-strip
         CHUNKA = NPIX(1)
         CHUNKB = EPIX(1) - SPIX(1) + 1
         OUTSIZ = CHUNKB * NN
         NA = (SPIX(2)-1)*NPIX(1) + 1
C 
         CALL STFCRE('middstat1',D_R4_FORMAT,F_X_MODE,F_IMA_TYPE,
     +               CHUNKA,IMNOX,STAT)
         CALL STFMAP(IMNOX,F_X_MODE,1,CHUNKA,IAV,PNTRX,STAT)
         CALL STFCRE('middstat2',D_R4_FORMAT,F_X_MODE,F_IMA_TYPE,
     +               OUTSIZ,IMNOY,STAT)
         CALL STFMAP(IMNOY,F_X_MODE,1,OUTSIZ,IAV,PNTRY,STAT)
         NB = 1
         DO 500,N=1,NN
            CALL STFGET(INO,NA,CHUNKA,IAV,MADRID(PNTRX),STAT)
            CALL DATMOV(MADRID(PNTRX),NB,CHUNKB,SPIX(1),MADRID(PNTRY))
            NA = NA + CHUNKA
            NB = NB + CHUNKB
500      CONTINUE
         CALL STFCLO(IMNOX,STAT)
C 
C  full frame
C 
      ELSE      
         INSIZ = NPIX(1) * NPIX(2) * NPIX(3)
         OUTSIZ = INSIZ
         CALL STFCRE('middstat2',D_R4_FORMAT,F_X_MODE,F_IMA_TYPE,
     +               OUTSIZ,IMNOY,STAT)
         CALL STFMAP(IMNOY,F_X_MODE,1,OUTSIZ,IAV,PNTRY,STAT)
         CALL STFGET(INO,1,INSIZ,IAV,MADRID(PNTRY),STAT)
      ENDIF
C 
C  sort array + return median
C 
      NA = (OUTSIZ+1)/2                             !index of median
      CALL SORTME(MADRID(PNTRY),ZBINS,OUTSIZ,NA,VALMED,STAT)
C 
      IF (STAT.EQ.-1) THEN
         CALL STTPUT('No pixels found with data in given interval...',N)
         DO 900, N=1,7
            RBUF(N) = -1.0
            IBUF(N) = -1
900      CONTINUE
         CALL STKWRR('OUTPUTR',RBUF,1,7,UNI,N)
         CALL STKWRI('OUTPUTI',IBUF,1,7,UNI,N)
         CALL STFCLO(IMNOY,N)                       !keep value of STAT
      ELSE
         CALL STKWRR('OUTPUTR',VALMED,8,1,UNI,STAT)
         CALL STFCLO(IMNOY,STAT)
      ENDIF
C 
      RETURN
      END

      SUBROUTINE PLOTI(IMNO)
C 
      IMPLICIT NONE
C
      INTEGER      IMNO
      INTEGER      UNI(1),NULO,STAT,IAV
C 
      REAL         RV
C 
      CALL STKRDR('OUTPUTR',10,1,IAV,RV,UNI,NULO,STAT) 
      IF (RV.GT.0.99) THEN
         CALL PLOHI(IMNO)          !only if we got a histogram
      ELSE
         CALL STTPUT('No histogram calculated => no plotting...',STAT)
      ENDIF
C 
      RETURN
      END

      SUBROUTINE DATMOV(A,OFF,CHUNK,SPIX,Y)
C 
      IMPLICIT NONE
C 
      INTEGER      OFF,CHUNK,SPIX
C 
      REAL         A(*),Y(*)
C 
      INTEGER   N,NN
C 
      NN = OFF
      DO 500,N=SPIX,SPIX+CHUNK-1
         Y(NN) = A(N)
         NN = NN + 1 
500   CONTINUE
C 
      RETURN
      END
