C @(#)ccdscan.for	19.1 (ESO-IPG) 02/25/03 14:16:39
C $Id: ccdscan.for,v 1.1.1.1 2003/06/04 14:49:46 cguirao Exp $   
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 $Id: ccdscan.for,v 1.1.1.1 2003/06/04 14:49:46 cguirao Exp $
      PROGRAM OVSCAN
C --------------------------------------------------------------------
C.COPYRIGHT: Copyright (c) 1983 European Southern Observatory,
C                                         all rights reserved
C.IDENTIFICATION: CCDOVSCAN.FOR
C.LANGUAGE:  F77+ESOext
C.AUTHOR:    R.H. Warmels
C.KEYWORDS:  average, overscan
C.PURPOSE:   produce a 1d image from a 2d by averaging over rows or columns
C.ALGORITHM  add rows/columns
C.VERSION:   900202  KB   take care of options in lower case
C.VERSION:   ??????  MP   cosmectic changes
C.VERSION:   910417  xx   replace numbers by symbolic constants in ST calls
C.VERSION:   910417  xx   replace numbers by symbolic constants in ST calls
C.VERSION:   930226  RHW  modification for ccd package
C.VERSION:   930331  RHW  modification for average row wise
C.VERSION:   980609  SW   correct use of the selected AREA (OV_SEC)
C ------------------------------------------------------------------

      IMPLICIT NONE
C
      INTEGER MADRID(1)
      INTEGER NDIM,KUN,KNUL
      INTEGER I1,I2,L1,L2,NI,NL
      INTEGER NAXIS,IMNO
      INTEGER NPIX(3)
      INTEGER*8 PNTRA,PNTRB,PNTRC
      INTEGER TIDO
      INTEGER NCOL, NROW, COL(3)
      INTEGER NDUM
      INTEGER STATUS
      INTEGER IMAGE(4)
      INTEGER IAV, ILOG
      INTEGER I, ISMO
C variable in first part not used anymore 05/00 PN
C      INTEGER IR
      INTEGER IXMIN, IXMAX
      INTEGER NCUR, KA, NPAR
      INTEGER ACCESS, PLMODE
      INTEGER NBINS, NACT, NNN, IX
      INTEGER IFREQ(512)
C
      REAL    AMIN, AMAX, XMIN, XMAX, STP
      REAL    X0, X1, XT
      
      REAL    CINT(512), CR
      REAL    LHCUTS(4)
      REAL    FRAME(8)
      REAL    XPA(2), YPA(2)
      REAL    XMEAN, XSTD
      REAL    XCUR, YCUR
      REAL    SCALES(3),OFFSET(2),CLIP(4)

      DOUBLE PRECISION STEP(3)
      DOUBLE PRECISION START(3)
C
      CHARACTER*60     FRAM, AREA
      CHARACTER*60     TABLE
      CHARACTER*80     CINPUT
      CHARACTER*80     TEXT, SEL
      CHARACTER*17     COLUMN(3)
      CHARACTER*80     LABEL1, LABEL2, LABEL3
      CHARACTER*72     CUNIT
      CHARACTER*1      DIRECT
      CHARACTER*72     IDENT
      CHARACTER        IOP*1
      CHARACTER        TABUNI(3)*16,TABLAB(3)*16,TABFOR(3)*16
C 
      INCLUDE          'MID_INCLUDE:ST_DEF.INC'
      COMMON           /VMR/ MADRID
      INCLUDE          'MID_INCLUDE:ST_DAT.INC'
C
      DATA             SEL     /' '/
      DATA             ILOG    /0/
      DATA             COLUMN  /':X', ':FLUX',':FIT'/
      DATA             NDIM    /2/
      DATA             IOP     /'A'/
      DATA             TABUNI  /' ', ' ', ' '/
      DATA             TABLAB  /'X', 'FLUX', 'FIT'/
      DATA             TABFOR  /'G13.6', 'G13.6','G13.6'/
      DATA             SCALES  /0.0,0.0,0.0/
      DATA             CLIP    /0.15,0.95,0.1,0.95/
      DATA             OFFSET  /-999., -999./
C
 9001 FORMAT('*** WARNING: zero dynamics range in y; '//
     2       'value = ',G13.6)
 9002 FORMAT('[@', I5.5, ',@', I5.5, ':@', I5.5, ',@', I5.5, ']')
 9003 FORMAT(':X.GE.',I5.5,'.AND.:X.LE.',I5.5)       
C
C  *** initialize system
      CALL STSPRO('OVSCAN')
C
C  *** read parameter input
      CALL STKRDC('CCDIN',1,1,60,IAV,FRAM,KUN,KNUL,STATUS)
      CALL STKRDC('CCDOUT',1,1,60,IAV,TABLE,KUN,KNUL,STATUS)
      CALL STKRDC('OV_SEC',1,1,60,IAV,AREA,KUN,KNUL,STATUS)
      CALL STKRDC('DIRECT',1,1,1,IAV,DIRECT,KUN,KNUL,STATUS)     ! read direct.
      CALL STKRDC('OV_IMODE',1,1,3,IAV,CINPUT,KUN,KNUL,STATUS)
      CALL STKRDI('INPUTI',1,1,IAV,ISMO,KUN,KNUL,STATUS)
C
C *** read input frame
      CALL STIGET(FRAM,D_R4_FORMAT,F_I_MODE,F_IMA_TYPE,NDIM,NAXIS,
     +            NPIX,START,STEP,IDENT,CUNIT,PNTRA,IMNO,STATUS)
C 
      CALL EXTCOO(IMNO,AREA,NAXIS,NDUM,IMAGE(1),IMAGE(3),STATUS) !pixel coords.
      IF (STATUS.NE.0) THEN
         TEXT = '*** FATAL: invalid coordinate input ...'
         CALL STTPUT(TEXT,STATUS)
         CALL STSEPI
      ENDIF
      I1 = IMAGE(1)
      I2 = IMAGE(3)
      L1 = IMAGE(2)
      L2 = IMAGE(4)

C      type *,"IMAGE: ",IMAGE

      NI = I2-I1+1
      NL = L2-L1+1
C
C *** map output table and fill the column
      NCOL = 3
      CALL UPCAS(DIRECT,DIRECT)
      IF ((DIRECT.EQ.'R') .OR. (DIRECT.EQ.'X'))THEN
         DIRECT = 'X'
         FRAME(1) = FLOAT(L1)
         FRAME(2) = FLOAT(L2)
c         NROW     = NPIX(2)
         NROW     = NL

      ELSE IF ((DIRECT.EQ.'C') .OR. (DIRECT.EQ.'Y')) THEN
         DIRECT = 'Y'
         FRAME(1) = FLOAT(I1)
         FRAME(2) = FLOAT(I2)
c         NROW     = NPIX(1)
         NROW     = NI
      ENDIF
      CALL TBTINI(TABLE,F_TRANS,F_IO_MODE,NCOL,NROW,TIDO,STATUS)! create table
      DO 20 IAV = 1,NCOL                                       ! define columns
         CALL TBCINI(TIDO,D_R4_FORMAT,1,TABFOR(IAV),TABUNI(IAV),
     2               TABLAB(IAV),COL(IAV),STATUS)
   20 CONTINUE
C
      CALL TBCSER(TIDO,COLUMN(1),COL(1),STATUS)             ! find pixel column
      CALL TBCSER(TIDO,COLUMN(2),COL(2),STATUS)             ! find flux column
      CALL TBCMAP(TIDO,COL(1),PNTRB,STATUS)                 ! map x column
      CALL TBCMAP(TIDO,COL(2),PNTRC,STATUS)                 ! map flux column
      IF (DIRECT.EQ.'X') THEN                               ! column readout
         CALL AVERC(NPIX(1),NPIX(2),MADRID(PNTRA),MADRID(PNTRB),
     +              MADRID(PNTRC),I1,I2,IOP,ISMO,LHCUTS(3),LHCUTS(4))
         CALL TBIPUT(TIDO,NCOL,NROW,STATUS)                 ! put the columns
c 98.06.09-SW (swolf@eso.org):
c The following selection part is not necessary
c as the table consits only of valid rows!
c         DO 21 IR = 1,L1
c            CALL TBSPUT(TIDO,IR,.FALSE.,STATUS)
c   21    CONTINUE
c         DO 22 IR = L2,NROW
c            CALL TBSPUT(TIDO,IR,.FALSE.,STATUS)
c   22    CONTINUE
C 
C *** get the statistics
         CALL TDMXRS(TIDO,COL(2),NROW,0,XMIN,XMAX)
         IF (XMIN.EQ.XMAX) THEN
            XMAX   = XMIN + 1.
         END IF
         CALL TDSCAL(XMIN,XMAX,1.,X0,X1,IX,XT,NNN)
         AMIN    = X0*10.**IX
         AMAX    = X1*10.**IX
         STP     = XT
         NBINS   = (AMAX-AMIN)/STP + 2
         NBINS   = MIN(NBINS,512)
         CR      = AMIN
         CINT(1) = CR
         DO 23 I = 1,NBINS - 2
            CR        = CR + STP
            CINT(I+1) = CR
 23      CONTINUE
         CALL TDRSTA(TIDO,COL(2),NROW,NBINS,1,CINT,
     2               IFREQ,NACT,XMEAN,XSTD)

      ELSE                                                  ! row readout
         CALL AVERR(NPIX(1),NPIX(2),MADRID(PNTRA),MADRID(PNTRB),
     +              MADRID(PNTRC),L1,L2,IOP,ISMO,LHCUTS(3),LHCUTS(4))
         CALL TBIPUT(TIDO,NCOL,NROW,STATUS)                 ! put the columns
c 98.06.09-SW (swolf@eso.org):
c The following selection part is not necessary
c as the table consits only of valid rows!
c         DO 31 IR = 1,I1
c            CALL TBSPUT(TIDO,IR,.FALSE.,STATUS)
c   31    CONTINUE
c         DO 32 IR = I2,NROW
c            CALL TBSPUT(TIDO,IR,.FALSE.,STATUS)
c   32    CONTINUE
C 
C *** get the statistics
         CALL TDMXRS(TIDO,COL(2),NROW,0,XMIN,XMAX)
         IF (XMIN.EQ.XMAX) THEN
            XMAX   = XMIN + 1.
         END IF
         CALL TDSCAL(XMIN,XMAX,1.,X0,X1,IX,XT,NNN)
         AMIN   = X0*10.**IX
         AMAX   = X1*10.**IX
         STP    = XT
         NBINS  = (AMAX-AMIN)/STP + 2
         NBINS  = MIN(NBINS,512)
         CR     = AMIN
         CINT(1) = CR
         DO 33 I = 1,NBINS - 2
            CR        = CR + STP
            CINT(I+1) = CR
 33      CONTINUE
         CALL TDRSTA(TIDO,COL(2),NROW,NBINS,1,CINT,
     2               IFREQ,NACT,XMEAN,XSTD)
      ENDIF
      LHCUTS(3) = XMEAN-3.*XSTD
      LHCUTS(4) = XMEAN+3.*XSTD
C
C *** the interactive mode
      CALL UPCAS(CINPUT,CINPUT)
      IF (CINPUT(1:1).EQ.'Y') THEN
         IF (LHCUTS(3).EQ.LHCUTS(4)) THEN
            WRITE(TEXT,9001) LHCUTS(3)
            CALL STTPUT(TEXT,STATUS)
            FRAME(5) = 0.95*LHCUTS(3)
            FRAME(6) = 1.05*LHCUTS(4)
         ELSE
            FRAME(5) = LHCUTS(3)
            FRAME(6) = LHCUTS(4)
         ENDIF
C
C *** get the frame ranges in x and y
         CALL GETAXS('AUTO',FRAME(1))
         CALL GETAXS('AUTO',FRAME(5))
         CALL PTKWRR('XWNDL',4,FRAME(1))
         CALL PTKWRR('YWNDL',4,FRAME(5))
C
C *** statistics determined; start the display
         CALL PTKWRR('SCALE',2,SCALES)
         CALL PTKWRR('OFFS',2,OFFSET)
C
         ACCESS = 0
         PLMODE = 1
         CALL PTOPEN(' ','overscan',ACCESS,PLMODE)             ! open plot file
         CALL PTKWRR('CLPL',4,CLIP)
C
         IF (DIRECT.EQ.'X') THEN
            LABEL1 = 'Pixel Row Number '
            LABEL2 = 'Pixel Intensity'
         ELSE
            LABEL1 = 'Pixel Column Number '
            LABEL2 = 'Pixel Intensity'
         ENDIF
         LABEL3 = 'TITLE=Frame: '//TABLE
         CALL PTKWRR('CLPL',4,CLIP)
         CALL PTAXES(FRAME(1),FRAME(5),LABEL1,LABEL2,LABEL3)   ! plot the frame
         NPAR = 3
         CALL PLTBL(TIDO,NPAR,COL(1),COL(2),NROW)               ! plot the data
C
C ***    Now do the looping for 2 cursers
         NCUR  = 2
         DO 50 I = 1,NCUR 
            CALL PTGCUR(XCUR,YCUR,KA,STATUS)
            IF (XCUR.LT.1.0) THEN
               XCUR = 1.0
            ENDIF
            IF (XCUR.GT.FLOAT(NROW)) THEN
               XCUR = FLOAT(NROW)
            ENDIF 
            XPA(I) = XCUR
            YPA(I) = YCUR
   50    CONTINUE
         CALL SORT2(XPA,YPA,NCUR)
         CALL AGGPLM(XPA,YPA,2,4)
         IXMIN = INT(XPA(1))
         IXMAX = INT(XPA(2))  
         CALL PTCLOS()
         WRITE(SEL,9003) IXMIN,IXMAX
         CALL STKWRC('SCANSEL',1,SEL,1,80,KUN,STATUS)      ! set selection flag
      ENDIF
C
C *** free data
      CALL TBTCLO(TIDO,STATUS)
      CALL STSEPI
      END


      SUBROUTINE AVERR(NPIX1,NPIX2,Z,X,Y,L1,L2,IOP,ISM,YMIN,YMAX)
C
C row average from row I1 to I2
C
      IMPLICIT    NONE
      INTEGER     NPIX1, NPIX2
      REAL        Z(NPIX1,NPIX2)
      REAL        X(NPIX1),Y(NPIX1)
      INTEGER     L1, L2
      CHARACTER*1 IOP
      INTEGER     ISM
      REAL        YMIN, YMAX
C 
      INTEGER     NPMAX
      PARAMETER   (NPMAX=10000)
      INTEGER     I, L
      INTEGER     IX, ISX, ILIS, IOFF, IXG
      INTEGER     NAV
      REAL        RN
      REAL        DATA(NPMAX)
      REAL        AVE
C
      DO 10 I = 1,NPIX1
          X(I)   = I
          Y(I)   = 0.
   10 CONTINUE
      IF ((IOP.EQ.'S').OR.(IOP.EQ.'s')) THEN
          DO 30 I = 1,NPIX1
              DO 20 L = L1,L2
                  DATA(I)   = DATA(I) + Z(I,L)
   20         CONTINUE
   30     CONTINUE
      ELSE
          RN     = L2 - L1 + 1
          DO 50 I = 1,NPIX1
              DO 40 L = L1,L2
                  DATA(I)   = DATA(I) + Z(I,L)/RN
   40         CONTINUE
   50     CONTINUE
      END IF
C
C *** do the smooting
      IF (ISM.GT.1) THEN
         ILIS = 2*(ISM/2) + 1
         IOFF = ILIS/2
         DO 60 IX = 1, NPIX1
            AVE = 0.0
            NAV  = 0
            DO 61 ISX = 1,ILIS
               IXG = IX + ISX - IOFF - 1
               IF (IXG.LE.0 .OR. IXG.GT.NPIX1) THEN
                  GO TO 61
               ENDIF
               NAV = NAV + 1
               AVE = AVE + DATA(IXG)
   61       CONTINUE
            Y(IX) = AVE/FLOAT(NAV)
   60    CONTINUE

      ELSE
         DO 70 IX = 1, NPIX1
            Y(IX) = DATA(IX)
   70    CONTINUE
      ENDIF
C
      YMAX   = -1.E30
      YMIN   = -YMAX
      DO 80 I = 1,NPIX1
          YMIN   = MIN(YMIN,Y(I))
          YMAX   = MAX(YMAX,Y(I))
   80 CONTINUE
      RETURN
      END


      SUBROUTINE AVERC(NPIX1,NPIX2,Z,X,Y,I1,I2,IOP,ISM,YMIN,YMAX)
C
C Column average from column 1 to column 2
C
      IMPLICIT    NONE
      INTEGER     NPIX1, NPIX2
      REAL        Z(NPIX1,NPIX2)
      REAL        X(NPIX2), Y(NPIX2)
      INTEGER     I1, I2
      CHARACTER*1 IOP
      INTEGER     ISM
      REAL        YMIN, YMAX
C
      INTEGER     NPMAX
      PARAMETER   (NPMAX=10000)
      INTEGER     I, J
      INTEGER     IX, ISX, ILIS, IOFF, IXG
      INTEGER     NAV
      REAL        RN
      REAL        DATA(NPMAX)
      REAL        AVE
C
C *** do the integration
      DO 10 I = 1,NPIX2
          X(I)   = I
          Y(I)   = 0.
   10 CONTINUE
      IF ((IOP.EQ.'S').OR.(IOP.EQ.'s')) THEN
          DO 30 I = 1,NPIX2
              DO 20 J = I1,I2
                  DATA(I)   = DATA(I) + Z(J,I)
   20         CONTINUE
   30     CONTINUE
      ELSE
          RN     = I2 - I1 + 1
          DO 50 I = 1,NPIX2
              DO 40 J = I1,I2
                  DATA(I)   = DATA(I) + Z(J,I)/RN
   40         CONTINUE
   50     CONTINUE
      END IF
C
C *** do the smooting
      IF (ISM.GT.1) THEN
         ILIS = 2*(ISM/2) + 1
         IOFF = ILIS/2
         DO 60 IX = 1, NPIX2
            AVE = 0.0
            NAV  = 0
            DO 61 ISX = 1,ILIS
               IXG = IX + ISX - IOFF - 1
               IF (IXG.LE.0 .OR. IXG.GT.NPIX2) THEN
                  GO TO 61
               ENDIF
               NAV = NAV + 1
               AVE = AVE + DATA(IXG)
   61       CONTINUE
            Y(IX) = AVE/FLOAT(NAV)
   60    CONTINUE

      ELSE
         DO 70 IX = 1, NPIX2
            Y(IX) = DATA(IX)
   70    CONTINUE
      ENDIF
C
C *** get the minimum and maximum
      YMAX   = -1.E30
      YMIN   = -YMAX
      DO 80 I = 1,NPIX2
          YMIN   = MIN(YMIN,Y(I))
          YMAX   = MAX(YMAX,Y(I))
   80 CONTINUE
      RETURN
      END



      SUBROUTINE SORT2(RA,RB,N)                     
C+++
C.IDENTIFICATION 
C.AUTHOR:   M.Rosa, ESO-Garching, 831105
C.PURPOSE:  Sort string of values into ascending sequence of values 
C.ALGORITM: The art of scientific computing
C           Press et al., Cambridge university Press, 1986   
C.USE:      CALL SORT2(XC,YC,N)                              
C.VERSION:  910115 RHW IMPLICIT NONE added
C---
C     
      IMPLICIT   NONE
      INTEGER    N
      REAL       RA(N)
      REAL       RB(N)                                                         
C
      INTEGER    I, IR
      INTEGER    J
      INTEGER    L
      REAL       RRA, RRB

      L  = N/2 + 1
      IR = N

10    CONTINUE
      IF (L.GT.1) THEN 
         L      = L - 1
         RRA    = RA(L)
         RRB    = RB(L)
      ELSE 
         RRA    = RA(IR)
         RRB    = RB(IR)
         RA(IR) = RA(1) 
         RB(IR) = RB(1) 
         IR     = IR - 1

         IF (IR.EQ.1) THEN
            RA(1) = RRA   
            RB(1) = RRB   
            RETURN        
         ENDIF
      ENDIF

      I = L
      J = L+L

20    CONTINUE
      IF (J.LE.IR) THEN 
         IF (J.LT.IR) THEN
            IF (RA(J).LT.RA(J+1)) THEN
               J = J + 1
            ENDIF
         ENDIF

         IF (RRA.LT.RA(J)) THEN
            RA(I) = RA(J) 
            RB(I) = RB(J) 
            I     = J 
            J     = J + J
         ELSE            
            J     = IR + 1
         ENDIF
         GO TO 20 
      ENDIF

      RA(I) = RRA 
      RB(I) = RRB
      GO TO 10
      END


      SUBROUTINE PLTBL(TID,NCOLUM,I1,I2,NROW)
C +++
C.PURPOSE:      Plots one or two columns of a table                           
C.AUTHOR:       Rein H. Warmels                                               
C.COMMENTS:     none                                                          
C.VERSION:      931103 RHW Createsd for environment document
C ---
      IMPLICIT    NONE
      INTEGER     TID               ! Table identifier
      INTEGER     NCOLUM            ! # of columns to be plotted 
      INTEGER     I1                ! index to first column
      INTEGER     I2                ! index to second column
      INTEGER     NROW              ! number of row
C
      INTEGER     NPMAX
      PARAMETER   (NPMAX=100000)
      INTEGER     IFIRST, I, ISTAT, IR, IAC
      INTEGER     STYPE, LTYPE
      REAL        VX, VY
      REAL        XPS(NPMAX), YPS(NPMAX)
      LOGICAL     IPLOT,ISEL,INULL
      CHARACTER   TEXT*80
C
 9001 FORMAT('*** FATAL: Maximum number of table entries is ',I8)
C
      IFIRST = 0
      IR     = 0
C
      DO 10 I = 1,NROW
         IPLOT  = .TRUE.
         CALL TBSGET(TID,I,ISEL,ISTAT)
         IF (ISEL) THEN            
            CALL TBRRDR(TID,I,1,I1,VX,INULL,ISTAT)
            IF (INULL) THEN
               IPLOT  = .FALSE.
            ENDIF
C 
            CALL TBRRDR(TID,I,1,I2,VY,INULL,ISTAT)
            IF (INULL) THEN 
               IPLOT  = .FALSE. 
            ENDIF
            IF (IPLOT) THEN
               IR      = IR + 1
               IF (IR.GT.NPMAX) THEN
                  WRITE (TEXT,9001) NPMAX
                  CALL STTPUT(TEXT,ISTAT)
                  CALL STSEPI
               ENDIF
               XPS(IR) = VX
               YPS(IR) = VY
            ENDIF
         ENDIF 
   10 CONTINUE 
C
      CALL PTKRDI('STYPE',1,IAC,STYPE)
      CALL PTKRDI('LTYPE',1,IAC,LTYPE)
      IF (LTYPE.EQ.0 .AND. STYPE.EQ.0) THEN 
         CALL STTPUT('*** FATAL: LTYPE and STYPE '//                   
     2               'both equal 0: NO PLOT',ISTAT)
         CALL PTCLOS()
         CALL STSEPI
      ELSE 
         CALL PTDATA(STYPE,LTYPE,0,XPS,YPS,0.0,IR)
      ENDIF
C
C *** end of the plotting
      RETURN
      END
