C @(#)ccdillum.for	19.1 (ES0-DMD) 02/25/03 14:16:38
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
      PROGRAM MKILLU
C++++++++++++++++++++++++++++++++++++++++++++++++++
C.IDENTIFICATION: ccdillum.FOR
C.KEYWORDS:       box smoothing
C.PURPOSE:        box smoothing with growing boxcar size
C.VERSION:        930810 RHW Creation
C--------------------------------------------------
      IMPLICIT     NONE

      INTEGER      IAV
      INTEGER      STATUS,MADRID
      INTEGER      NAXIS,NPIX(2)
      INTEGER      IMNOC,IMNOA
      INTEGER*8    PNTRA,PNTRC
      INTEGER      KUN,KNUL
      INTEGER      XBMIN, XBMAX
      INTEGER      YBMIN, YBMAX
C
      CHARACTER*60 FRAMEA,FRAMEC
      CHARACTER*72 IDENT
      CHARACTER*48 CUNIT
      CHARACTER*3  ACTION
C
      DOUBLE PRECISION START(2),STEP(2)
C
      REAL         AVER
      REAL         XBOX(2), YBOX(2)
      REAL         XBMINR, XBMAXR
      REAL         YBMINR, YBMAXR
      REAL         CLIP(2)
      REAL         LOWSIG,HIGSIG
C
      INCLUDE 'MID_INCLUDE:ST_DEF.INC'
      COMMON/VMR/MADRID(1)
      INCLUDE 'MID_INCLUDE:ST_DAT.INC'
C
C *** do the work
      CALL STSPRO('ILLUM')

      CALL STKRDC('IN_A',1,1,60,IAV,FRAMEA,KUN,KNUL,STATUS)
      CALL STIGET(FRAMEA,D_R4_FORMAT,F_I_MODE,F_IMA_TYPE,
     1                2,NAXIS,NPIX,START,STEP,IDENT,CUNIT,
     2                PNTRA,IMNOA,STATUS)
      CALL STKRDC('OUT_A',1,1,60,IAV,FRAMEC,KUN,KNUL,STATUS)
      CALL STIPUT(FRAMEC,D_R4_FORMAT,F_O_MODE,F_IMA_TYPE,
     1            NAXIS,NPIX,START,STEP,IDENT,CUNIT,PNTRC,
     2            IMNOC,STATUS)
C
C *** smoothing boxes
      CALL STKRDR('XBOX',1,2,IAV,XBOX,KUN,KNUL,STATUS)
      CALL STKRDR('YBOX',1,2,IAV,YBOX,KUN,KNUL,STATUS)
      XBMINR = XBOX(1)
      XBMAXR = XBOX(2)
      YBMINR = YBOX(1)
      YBMAXR = YBOX(2)
C
C *** check the smoothing boxes
      IF (XBMINR.LT.1.) THEN
         XBMINR = XBMINR*NPIX(1)
      ENDIF 
      IF (XBMAXR.LT.1.) THEN
         XBMAXR = XBMAXR*NPIX(1)
      ENDIF      
      IF (YBMINR.LT.1.) THEN
         YBMINR = YBMINR*NPIX(2)
      ENDIF      
      IF (YBMAXR.LT.1.) THEN
         YBMAXR = YBMAXR*NPIX(2)
      ENDIF
C
      XBMIN = INT(MAX(2,MIN(NPIX(1),NINT(MIN(XBMINR,XBMAXR)))))
      XBMAX = INT(MAX(2,MIN(NPIX(1),NINT(MAX(XBMINR,XBMAXR)))))
      YBMIN = INT(MAX(2,MIN(NPIX(2),NINT(MIN(YBMINR,YBMAXR)))))
      YBMAX = INT(MAX(2,MIN(NPIX(2),NINT(MAX(YBMINR,YBMAXR)))))
C
C *** to clip or not to clip
      CALL STKRDC('CLIP',1,1,3,IAV,ACTION,KUN,KNUL,STATUS)
      CALL UPCAS(ACTION,ACTION)
      IF (ACTION(1:1).EQ.'Y') THEN
         CALL STKRDR('SIGMA',1,2,IAV,CLIP,KUN,KNUL,STATUS)
         LOWSIG = CLIP(1)
         HIGSIG = CLIP(2)
      ENDIF
C
C *** test on clipping
      IF (ACTION(1:1).EQ.'Y') THEN
         CALL ILLUM(PNTRA,PNTRC,NPIX,XBMIN,XBMAX,YBMIN,YBMAX,
     2              LOWSIG,HIGSIG,AVER)
      ELSE
         CALL QILLUM(PNTRA,PNTRC,NPIX,XBMIN,XBMAX,YBMIN,YBMAX,AVER)
      ENDIF
C
      CALL STDCOP(IMNOA,IMNOC,1,' ',STATUS) 
      CALL STDWRR(IMNOC,'CCDMEAN',AVER,1,1,KUN,STATUS)
      CALL STSEPI
      END


      SUBROUTINE ILLUM(IN,OUT,NPIX,XBMIN,XBMAX,YBMIN,YBMAX,
     2                 LOW,HIGH,MEAN)
C
      IMPLICIT NONE

      INTEGER       IN, OUT
      INTEGER       NPIX(2)
      INTEGER       XBMIN,XBMAX,YBMIN,YBMAX
      REAL          LOW, HIGH
      REAL          MEAN
C
      INTEGER       LINEIN, LINEOUT
      INTEGER       IMNOS, IMNOG, IMNOP
      INTEGER*8     IPOINT, PTRS, SUM, AVG, OUTPUT
      INTEGER*8     PTR
      INTEGER       I, NITER, NREJ
      INTEGER       MADRID(1)
      INTEGER       NCOLS,NLINES
      INTEGER       STAT, IAV
      INTEGER       YBOX2

      REAL          SCALE
      REAL          ILLSUM
      EXTERNAL      ILLSUM
      INTEGER       BXCLN
      EXTERNAL      BXCLN
      INCLUDE       'MID_INCLUDE:ST_DEF.INC'
      COMMON        /VMR/MADRID
      INCLUDE       'MID_INCLUDE:ST_DAT.INC'
      DATA          NITER/10/
C
C *** here create vitual memory
      NCOLS  = NPIX(1)
      NLINES = NPIX(2)
C      
      CALL STFCRE('dummsumm',D_R4_FORMAT,F_X_MODE,F_IMA_TYPE,
     2            NCOLS,IMNOS,STAT)                          ! create sum array
      CALL STFMAP(IMNOS,F_X_MODE,1,NCOLS,IAV,SUM,STAT)
C
      CALL STFCRE('dummavg',D_R4_FORMAT,F_X_MODE,F_IMA_TYPE,
     2            NCOLS,IMNOG,STAT)                        !create output array
      CALL STFMAP(IMNOG,F_X_MODE,1,NCOLS,IAV,AVG,STAT)
      OUTPUT = AVG
C
      CALL STFCRE('dummpntr',D_I4_FORMAT,F_X_MODE,F_IMA_TYPE,
     2            YBMAX,IMNOP,STAT)                      ! create pointer array
      CALL STFMAP(IMNOP,F_X_MODE,1,YBMAX,IAV,PTRS,STAT)
C
C *** now start the work
C      IF (YBMAX.LT.NLINES) THEN
C         YBMAX = NLINES
C      ENDIF
C
C *** accumulate the minimum y box
      CALL ILLZER(MADRID(SUM),NCOLS)                   ! initialize sum array
C
C *** acumulate the minimum y box
      LINEIN = 0
   10 CONTINUE
      IF (LINEIN.LT.YBMIN) THEN
         LINEIN = LINEIN + 1
         IPOINT = IN + (LINEIN-1)*NCOLS
         CALL ILLADD(MADRID(IPOINT),MADRID(SUM),MADRID(SUM),NCOLS)
         PTR = PTRS + MOD(LINEIN,YBMAX) 
         MADRID(PTR) = IPOINT
         GOTO 10
      ENDIF
C
C *** output the minimum y box
      YBOX2 = YBMIN
      SCALE = FLOAT(YBMIN)
      CALL AGBCAR(MADRID(SUM),MADRID(AVG),NCOLS,XBMIN,XBMAX,SCALE)
C
C *** Iteratively clean the initial lines
      PTR = PTRS
      IF (YBOX2 .NE. YBMAX) THEN
         PTR = PTR+1
      ENDIF
C
      DO I = 1, NITER
         NREJ = 0
         DO LINEOUT = 1,LINEIN
            IPOINT = MADRID(PTR+LINEOUT-1)
            NREJ   = NREJ + BXCLN(MADRID(IPOINT),MADRID(AVG),
     2                            MADRID(SUM),NCOLS,LOW,HIGH)
         ENDDO
         IF (NREJ.GT.0) THEN
            CALL AGBCAR(MADRID(SUM),MADRID(AVG),NCOLS,XBMIN,XBMAX,SCALE)
         ELSE
            GOTO 20
         ENDIF
      ENDDO
C
   20 CONTINUE
      YBOX2 = (YBMIN+1)/2
      LINEOUT = 0
          
   30 CONTINUE
      IF (LINEOUT.LT.YBOX2) THEN
         LINEOUT = LINEOUT + 1
         IPOINT  = OUT + (LINEOUT-1)*NCOLS
         CALL COPYF(MADRID(OUTPUT),MADRID(IPOINT),NCOLS)
         GO TO 30
      ENDIF
      MEAN = YBOX2*ILLSUM(MADRID(OUTPUT),NCOLS)            ! calculate the MEAN
C
C *** increase the y box size by steps of 2 until the maximum szie.
   40 CONTINUE
      IF (LINEIN.LT.YBMAX) THEN
         LINEIN = LINEIN + 1
         IPOINT = IN + (LINEIN-1)*NCOLS
         CALL ILLADD(MADRID(SUM),MADRID(IPOINT),MADRID(SUM),NCOLS) 
         PTR = PTRS + MOD(LINEIN,YBMAX)
         MADRID(PTR) = IPOINT
         SCALE = SCALE + 1
C
         NREJ  = BXCLN(MADRID(IPOINT),MADRID(AVG),
     2                 MADRID(SUM),NCOLS,LOW,HIGH)
         CALL AGBCAR(MADRID(SUM),MADRID(AVG),NCOLS,XBMIN,XBMAX,SCALE)
C
         LINEIN = LINEIN + 1
         IPOINT = IN + (LINEIN-1)*NCOLS
         CALL ILLADD(MADRID(SUM),MADRID(IPOINT),MADRID(SUM),NCOLS) 
         PTR    = PTRS + MOD(LINEIN,YBMAX)
         MADRID(PTR) = IPOINT
C
         NREJ  = BXCLN(MADRID(IPOINT),MADRID(AVG),MADRID(SUM),
     2           NCOLS,LOW,HIGH)
         SCALE = SCALE + 1
         CALL AGBCAR(MADRID(SUM),MADRID(AVG),NCOLS,XBMIN,XBMAX,SCALE)
C
         LINEOUT = LINEOUT + 1
         IPOINT = OUT + (LINEOUT-1)*NCOLS
         CALL COPYF(MADRID(AVG),MADRID(IPOINT),NCOLS)
         MEAN = MEAN + ILLSUM(MADRID(IPOINT),NCOLS)            ! calculate MEAN
         GOTO 40
      ENDIF
C
C *** for each line subtract the last line from the sum, add the 
C     next line to the sum, and output a line
   50 CONTINUE
      IF (LINEIN.LT.NLINES) THEN
         LINEIN = LINEIN + 1
         PTR = PTRS + MOD(LINEIN,YBMAX)
         IPOINT = MADRID(PTR)
         CALL ILLSUB(MADRID(SUM),MADRID(IPOINT),MADRID(SUM),NCOLS) 
         IPOINT = IN + (LINEIN-1)*NCOLS
         CALL ILLADD(MADRID(SUM),MADRID(IPOINT),MADRID(SUM),NCOLS) 
         MADRID(PTR) = IPOINT
C
         NREJ  = BXCLN(MADRID(IPOINT),MADRID(AVG),MADRID(SUM),
     2           NCOLS,LOW,HIGH)
C
         LINEOUT = LINEOUT + 1
         IPOINT = OUT + (LINEOUT-1)*NCOLS
         CALL AGBCAR(MADRID(SUM),MADRID(AVG),NCOLS,XBMIN,XBMAX,SCALE)
         CALL COPYF(MADRID(AVG),MADRID(IPOINT),NCOLS)
         MEAN = MEAN + ILLSUM(MADRID(IPOINT),NCOLS)        ! calculate the MEAN
         GOTO 50
      ENDIF
C
C *** decrease the y box in steps of 2 until mimimum y box
   60 CONTINUE
      IF (LINEOUT.LT.(NLINES-YBOX2)) THEN
         LINEIN = LINEIN + 1
         PTR = PTRS + MOD(LINEIN,YBMAX)
         IPOINT = MADRID(PTR)
         CALL ILLSUB(MADRID(SUM),MADRID(IPOINT),MADRID(SUM),NCOLS) 
         LINEIN = LINEIN + 1
         PTR = PTRS + MOD(LINEIN,YBMAX)
         IPOINT = MADRID(PTR)
         CALL ILLSUB(MADRID(SUM),MADRID(IPOINT),MADRID(SUM),NCOLS) 
         SCALE = SCALE-2
C
         LINEOUT  = LINEOUT + 1
         IPOINT = OUT + (LINEOUT-1)*NCOLS
         CALL AGBCAR(MADRID(SUM),MADRID(IPOINT),NCOLS,XBMIN,XBMAX,SCALE)
         MEAN = MEAN + ILLSUM(MADRID(IPOINT),NCOLS)        ! calculate the MEAN
         GOTO 60
      ENDIF
C
C *** output the last liens of the minimum y box size
      CALL AGBCAR(MADRID(SUM),MADRID(AVG),NCOLS,XBMIN,XBMAX,SCALE)
   70 CONTINUE
      IF (LINEOUT.LT.NLINES) THEN
         LINEOUT = LINEOUT + 1
         IPOINT = OUT + (LINEOUT-1)*NCOLS
         CALL COPYF(MADRID(OUTPUT),MADRID(IPOINT),NCOLS)
         GOTO 70
      ENDIF
C
C *** write the scaling vector
      MEAN = MEAN + (NLINES-LINEOUT)*ILLSUM(MADRID(AVG),NCOLS)
      MEAN = MEAN/(NCOLS*NLINES)
C
      RETURN
      END 



      SUBROUTINE QILLUM(IN,OUT,NPIX,XBMIN,XBMAX,YBMIN,YBMAX,MEAN)
C
      IMPLICIT NONE

      INTEGER       IN, OUT
      INTEGER       NPIX(2)
      INTEGER       XBMIN,XBMAX,YBMIN,YBMAX
      REAL          MEAN
C
      INTEGER       LINEIN, LINEOUT
      INTEGER       IMNOS, IMNOO, IMNOP
      INTEGER       IPOINT, PTRS, SUM, OUTPUT
      INTEGER       PTR
      INTEGER       MADRID(1)
      INTEGER       NCOLS,NLINES
      INTEGER       STAT, IAV
      INTEGER       YBOX1

      REAL          SCALE
      REAL          ILLSUM
      EXTERNAL      ILLSUM
      INCLUDE       'MID_INCLUDE:ST_DEF.INC'
      COMMON        /VMR/MADRID
      INCLUDE       'MID_INCLUDE:ST_DAT.INC'
C
C *** here create vitual memory
      NCOLS  = NPIX(1)
      NLINES = NPIX(2)
C      
      CALL STFCRE('dummsumm',D_R4_FORMAT,F_X_MODE,F_IMA_TYPE,
     2            NCOLS,IMNOS,STAT)                          ! create sum array
      CALL STFMAP(IMNOS,F_X_MODE,1,NCOLS,IAV,SUM,STAT)
C
      CALL STFCRE('dummoutp',D_R4_FORMAT,F_X_MODE,F_IMA_TYPE,
     2            NCOLS,IMNOO,STAT)                        !create output array
      CALL STFMAP(IMNOO,F_X_MODE,1,NCOLS,IAV,OUTPUT,STAT)
C
      CALL STFCRE('dummpntr',D_I4_FORMAT,F_X_MODE,F_IMA_TYPE,
     2            YBMAX,IMNOP,STAT)                      ! create pointer array
      CALL STFMAP(IMNOP,F_X_MODE,1,YBMAX,IAV,PTRS,STAT)
C
C *** now start the work
C      IF (YBMAX.LT.NLINES) THEN
C         YBMAX = NLINES
C      ENDIF
C
C *** accumulate the minimum y box
      CALL ILLZER(MADRID(SUM),NCOLS)                   ! initialize sum array
C
C *** acumulate the minimum y box
      LINEIN = 0
   10 CONTINUE
      IF (LINEIN.LT.YBMIN) THEN
         LINEIN = LINEIN + 1
         IPOINT = IN + (LINEIN-1)*NCOLS
         CALL ILLADD(MADRID(IPOINT),MADRID(SUM),MADRID(SUM),NCOLS)
         PTR = PTRS + MOD(LINEIN,YBMAX) 
         MADRID(PTR) = IPOINT
         GOTO 10
      ENDIF
C
C *** output the minimum y box
      YBOX1 = (YBMIN+1)/2
      SCALE = FLOAT(YBMIN)
      CALL AGBCAR(MADRID(SUM),MADRID(OUTPUT),NCOLS,XBMIN,XBMAX,SCALE)
      LINEOUT = 0
   20 CONTINUE
      IF (LINEOUT.LT.YBOX1) THEN
         LINEOUT = LINEOUT + 1
         IPOINT  = OUT + (LINEOUT-1)*NCOLS
         CALL COPYF(MADRID(OUTPUT),MADRID(IPOINT),NCOLS)
         GO TO 20
      ENDIF
      MEAN = YBOX1*ILLSUM(MADRID(OUTPUT),NCOLS)            ! calculate the MEAN
C
C *** increase the y box size by steps of 2 until the maximum szie.
   30 CONTINUE
      IF (LINEIN.LT.YBMAX) THEN
         LINEIN = LINEIN + 1
         IPOINT = IN + (LINEIN-1)*NCOLS
         CALL ILLADD(MADRID(SUM),MADRID(IPOINT),MADRID(SUM),NCOLS) 
         PTR    = PTRS + MOD(LINEIN,YBMAX)
         MADRID(PTR) = IPOINT
         LINEIN = LINEIN + 1
         IPOINT = IN + (LINEIN-1)*NCOLS
         CALL ILLADD(MADRID(SUM),MADRID(IPOINT),MADRID(SUM),NCOLS) 
         PTR    = PTRS + MOD(LINEIN,YBMAX)
         MADRID(PTR) = IPOINT
C
         SCALE = SCALE + 2
         LINEOUT = LINEOUT + 1
         IPOINT = OUT + (LINEOUT-1)*NCOLS
         CALL AGBCAR(MADRID(SUM),MADRID(IPOINT),NCOLS,XBMIN,XBMAX,SCALE)
         MEAN = MEAN + ILLSUM(MADRID(IPOINT),NCOLS)            ! calculate MEAN
         GOTO 30
      ENDIF
C
C *** for each line subtract the last line from the sum, add the 
C     next line to the sum, and output a line
   40 CONTINUE
      IF (LINEIN.LT.NLINES) THEN
         LINEIN = LINEIN + 1
         PTR = PTRS + MOD(LINEIN,YBMAX)
         IPOINT = MADRID(PTR)
         CALL ILLSUB(MADRID(SUM),MADRID(IPOINT),MADRID(SUM),NCOLS) 
         IPOINT = IN + (LINEIN-1)*NCOLS
         CALL ILLADD(MADRID(SUM),MADRID(IPOINT),MADRID(SUM),NCOLS) 
         MADRID(PTR) = IPOINT
C
         LINEOUT = LINEOUT + 1
         IPOINT = OUT + (LINEOUT-1)*NCOLS
         CALL AGBCAR(MADRID(SUM),MADRID(IPOINT),NCOLS,XBMIN,XBMAX,SCALE)
         MEAN = MEAN + ILLSUM(MADRID(IPOINT),NCOLS)        ! calculate the MEAN
         GOTO 40
      ENDIF
C
C *** decrease the y box in steps of 2 until mimimum y box
   50 CONTINUE
      IF (LINEOUT.LT.(NLINES-YBOX1)) THEN
         LINEIN = LINEIN + 1
         PTR = PTRS + MOD(LINEIN,YBMAX)
         IPOINT = MADRID(PTR)
         CALL ILLSUB(MADRID(SUM),MADRID(IPOINT),MADRID(SUM),NCOLS) 
         LINEIN = LINEIN + 1
         PTR = PTRS + MOD(LINEIN,YBMAX)
         IPOINT = MADRID(PTR)
         CALL ILLSUB(MADRID(SUM),MADRID(IPOINT),MADRID(SUM),NCOLS) 
C
         LINEOUT  = LINEOUT + 1
         SCALE = SCALE -2
         IPOINT = OUT + (LINEOUT-1)*NCOLS
         CALL AGBCAR(MADRID(SUM),MADRID(IPOINT),NCOLS,XBMIN,XBMAX,SCALE)
         MEAN = MEAN + ILLSUM(MADRID(IPOINT),NCOLS)        ! calculate the MEAN
         GOTO 50
      ENDIF
C
C *** output the last liens of the minimum y box size
      CALL AGBCAR(MADRID(SUM),MADRID(OUTPUT),NCOLS,XBMIN,XBMAX,SCALE)
   60 CONTINUE
      IF (LINEOUT.LT.NLINES) THEN
         LINEOUT = LINEOUT + 1
         IPOINT = OUT + (LINEOUT-1)*NCOLS
         CALL COPYF(MADRID(OUTPUT),MADRID(IPOINT),NCOLS)
         GOTO 60
      ENDIF
C
C *** write the scaling vector
      MEAN = MEAN + (NLINES-LINEOUT)*ILLSUM(MADRID(OUTPUT),NCOLS)
      MEAN = MEAN/(NCOLS*NLINES)
C
      RETURN
      END 
C
C
C
      INTEGER FUNCTION BXCLN(DATA,BXAVG,SUM,NCOLS,LOW,HIGH)
C
C ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C   Reject data values from thr sum for the next boxcar average which
C   exceed the minimum and maximum residual values from the current
C   boxcar average. This excludes data from the moving avarge before it 
C   enters the average.
C-----------------------------------------------------------------------
C
      IMPLICIT NONE
   
      REAL     DATA(1)
      REAL     BXAVG(1)
      REAL     SUM(1)
      INTEGER  NCOLS
      REAL     LOW
      REAL     HIGH
C
      INTEGER  I,NREJ
      REAL     RMS,RESID,MINRES,MAXRES
C
      RMS    = 0.0
      DO I = 1,NCOLS
         RMS = RMS + (DATA(I)-BXAVG(1))**2
      ENDDO
      RMS    = SQRT(RMS/NCOLS)
      MINRES = -LOW*RMS
      MAXRES = HIGH*RMS
C
      NREJ = 0
      DO I = 1,NCOLS
         RESID = DATA(I)-BXAVG(I)
         IF ((RESID.LT.MINRES) .OR. (RESID.GT.MAXRES)) THEN
            DATA(I) = BXAVG(I)
            SUM(I)  = SUM(I)-RESID
            NREJ    = NREJ+1
         ENDIF
      ENDDO
C
      BXCLN = NREJ
      RETURN
      END
C
C
C
      SUBROUTINE ILLZER(A,NC)
C
      IMPLICIT NONE
      REAL     A(1)
      INTEGER  NC
      INTEGER  IC
C
      DO IC = 1,NC
         A(IC) = 0.0
      ENDDO
      RETURN
      END
C
C
C
      SUBROUTINE ILLADD(A,B,C,NC)
C
      IMPLICIT NONE
      REAL     A(1),B(1),C(1)
      INTEGER  NC
      INTEGER  IC
C
      DO IC = 1,NC
         C(IC) = A(IC)+B(IC)
      ENDDO
      RETURN
      END
C
C
C
      SUBROUTINE ILLSUB(A,B,C,NC)
C
      IMPLICIT NONE
      REAL     A(1),B(1),C(1)
      INTEGER  NC
      INTEGER  IC
C
      DO IC = 1,NC
         C(IC) = A(IC)-B(IC)
      ENDDO
      RETURN
      END
C
C
C
      REAL FUNCTION ILLSUM(A,NC)
C
      IMPLICIT NONE
      REAL     A(1)
      INTEGER  NC
      INTEGER  IC
C
      ILLSUM = 0.0       
      DO IC = 1,NC
         ILLSUM = ILLSUM + A(IC)
      ENDDO

      RETURN
      END
C
C
C

      SUBROUTINE AGBCAR(INA,OUTA,NCOLS,XBMIN,XBMAX,YBOX)
C 
C +++ Vector growing boxcar smooth. Taken from the IRAF CCD package
C
      REAL     INA(1)
      REAL     OUTA(1)
      INTEGER  NCOLS
      INTEGER  XBMIN,XBMAX
      REAL     YBOX
C
      INTEGER  COLIN, COLOUT, LASTCOL
      INTEGER  NPIX, XBMIN2
      REAL     SUM,OUTPUT
C
      XBMIN2 = (XBMIN+1)/2
      COLIN  = 0 
      SUM    = 0.0
   10 CONTINUE
      IF (COLIN.LT.XBMIN) THEN
         COLIN = COLIN + 1
         SUM   =  SUM + INA(COLIN)
         GOTO 10
      ENDIF
C
      NPIX   = XBMIN * YBOX
      OUTPUT = SUM/NPIX
      COLOUT = 0
   20 CONTINUE
      IF (COLOUT.LT.XBMIN2) THEN
         COLOUT       = COLOUT + 1
         OUTA(COLOUT) = OUTPUT
         GO TO 20
      ENDIF 
C
   30 CONTINUE
      IF (COLIN.LT.XBMAX) THEN
         COLIN  = COLIN + 1
         SUM    = SUM + INA(COLIN)
         COLIN  = COLIN + 1
         SUM    = SUM + INA(COLIN)
         NPIX   = NPIX + 2*YBOX
         COLOUT = COLOUT + 1
         OUTA(COLOUT) = SUM/NPIX
         GO TO 30
      ENDIF
C
      LASTCOL = 0
   40 CONTINUE
      IF (COLIN.LT.NCOLS) THEN
         COLIN   = COLIN + 1
         LASTCOL = LASTCOL +1
         SUM     = SUM + INA(COLIN) - INA(LASTCOL)
         COLOUT  = COLOUT + 1
         OUTA(COLOUT) = SUM/NPIX
         GO TO 40
      ENDIF
C
   50 CONTINUE
      IF (COLOUT.LT. NCOLS - XBMIN2) THEN
         LASTCOL     = LASTCOL + 1
         SUM         = SUM - INA(LASTCOL)
         LASTCOL     = LASTCOL + 1
         SUM         = SUM - INA(LASTCOL)
         NPIX        = NPIX - 2*YBOX
         COLOUT      = COLOUT + 1
         OUTA(COLOUT)= SUM/NPIX
         GO TO 50
      ENDIF
C
      OUTPUT = SUM/NPIX
   70 CONTINUE
      IF (COLOUT.LT.NCOLS) THEN
         COLOUT = COLOUT + 1
         OUTA(COLOUT) = OUTPUT
         GO TO 70
      ENDIF

      RETURN
      END




