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 MODIF
C
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C.IDENTIFICATION              
C  program MODIF 			version 1.00	890630
C  K. Banse 				ESO - Garching
C
C.KEYWORDS    
C  modify pixels
C
C.PURPOSE
C  merge of MODPIX.FOR, MODAR.FOR, MODCUT.FOR
C
C.ALGORITHM
C  see the subroutines
C
C.INPUT/OUTPUT
C  the following keywords are used:
C  ACTION/C/1/2			action = MP for modify/pixel
C  				       = MA for modify/area
C  				       = MC for modify/cuts
C.VERSIONS
C  1.00		created 
C 
C 051014		last modif
C-------------------------------------------------------------------
C
      IMPLICIT NONE
C      
      INTEGER      IAV,STAT
      INTEGER      UNI(1),NULO
C 
      CHARACTER    ACTION*2
C
C  get into MIDAS environment
      CALL STSPRO('MODIF')
C 
      CALL STKRDC('ACTION',1,1,2,IAV,ACTION,UNI,NULO,STAT)
      CALL UPCAS(ACTION,ACTION)
C 
C  branch according to option

      IF (ACTION .EQ. 'MP') THEN
         CALL MODPIX(1)
      ELSE IF (ACTION .EQ. 'MA') THEN
         CALL MODPIX(2)
      ELSE IF (ACTION .EQ. 'MC') THEN
         CALL STKRDC('P2',1,1,1,IAV,ACTION,UNI,NULO,STAT)
         CALL MODCUT(ACTION)
      ENDIF
C   
      CALL STSEPI
      END

      SUBROUTINE MODPIX(MODMOD)
C
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C.IDENTIFICATION              
C  subroutine MODPIX 			version 1.00	881028
C  K. Banse 				ESO - Garching
C  1.50  980923 S.Wolf
C
C.KEYWORDS    
C  bulk data frame, bad pixels
C
C.PURPOSE
C  to replace a group of undefined or bad pixels values by surface 
C  interpolation of neighboring area. The bad pixel area is defined via the 
C  cursor.
C
C.ALGORITHM
C  a)   (original input from H. Waldthausen, had to be cleaned a lot...)
C  a surface is fitted using the pixels surrounding the cursor area
C  the surrounding area is chosen as 
C                       x-pixels = x-area factor * x-size of cursor rectangle
C                       y-pixels = y-area factor * y-size of cursor rectangle
C  the program forms an approximation to the weighted, least-squares Chebyshev 
C  series surface fit; the weights are derived from the Poisson error.
C
C  b)  (original input from A. Lauberts)
C  a LSQ surface is fitted using the pixels around a circle but inside
C  the cursor rectangle
C 
C.INPUT/OUTPUT
C             call as MODPIX(MODMOD)
C 
C  MODMOD:	Integer		1 = Modify/Pixel
C				2 = Modify/Area
C 
C  the following keywords are used:
C  a)
C  P1/C/1/60			input frame, table
C				or 
C				CURSOR
C  IN_A/C/1/60			input frame, if P1 = CURSOR
C  OUT_A/C/1/60			result frame
C  INPUTI/I/1/6			area factors in x,y
C				x,y degree of fitting polynomial, 
C				no. of iterations + overlay_draw_flag
C  NOISE/C/1/1                  Y or N - whether to add noise
C  IMASEC/C/1/100               allowed image section
C	
C  b)
C  P1/C/1/80			input frame, table
C				or 
C				CURSOR
C  IN_A/C/1/60			input frame, if P1 = CURSOR
C  OUT_A/C/1/60			result frame
C  INPUTI/I/1/2			degree of fitted SURFACs, draw flag
C  INPUTR/R/1/1			constant in case that INPUTI = 0
C 
C-------------------------------------------------------------------
C
      IMPLICIT NONE
C      
      INTEGER      MODMOD
      INTEGER      ICX,ICY,IPX,IPY
      INTEGER      ISX,ISY,IT,XDEG,YDEG
      INTEGER      IPIXS(4),ISUR(4)
      INTEGER      DEGREE,IAV,IDRAW
      INTEGER      N,N2,NAXIS,NCOL,NROW,NMAL
      INTEGER*8    PNTRA,PNTRB
      INTEGER      IMNOA,IMNOB,SIZE,LFLAGS(2)
      INTEGER      STAT1,STAT2,STAT,TABFLG
      INTEGER      NPIXA(3),KCOOR(4),WPIX(3)
      INTEGER      SPIX(3),EPIX(3)
      INTEGER      IXLO,IXUP,IYLO,IYUP,INPUTI(6),DAZOP
      INTEGER      ICUR1(2),ICUR2(2),WW(4)
      INTEGER      TID,TABCOLN(4)
      INTEGER      LOW1(3),HI(3),CUTIND,COOFF
      INTEGER      UNI(1),NUL,MADRID(1)
      INTEGER      XFIGU(6),YFIGU(8)
C      
      CHARACTER    FRAMEA*80,TABLE*80,FRAMEB*80
      CHARACTER*80  OUTPUT
      CHARACTER*100  CIMASEC
      CHARACTER      CUNITA*64,IDENTA*72,CACT*8,CNOISE*1
      CHARACTER      TABLAB(4)*16
      CHARACTER*40   ERMS12,ERMS13,ERMS15,ERMS16
      CHARACTER      INFO*44,TIP15*56,TIP12*28
      CHARACTER      CSYST
C      
      DOUBLE PRECISION STARTA(3),STEPA(3)
      DOUBLE PRECISION DD1(4),DD2(4)                     !for WCS conversions
C 
      REAL         RBUF(6),CUTS(4),COOR(4),STPINV(2)
      REAL         RCUTS(2),ZCUTS(2),PCUR1(6),PCUR2(6),CONST
      REAL         OFFSET,OMIN,CRMS,PCEN(6)
      REAL         INPUTR(8)
C      
      LOGICAL      NOISE,TABNUL(4),SELFLG
C 
      INCLUDE  'MID_INCLUDE:ST_DEF.INC'
C 
      INCLUDE  'MID_INCLUDE:IDIDEV.INC'
      INCLUDE  'MID_INCLUDE:IDIMEM.INC'
C 
      COMMON    /VMR/  MADRID
C
C  messages
      DATA
     +  INFO /' switch cursors on - next time we exit ... '/,
     +  ERMS12 /' ERR: surface dimension too large  '/,
     +  ERMS13 /' ERR: in mathematical subroutines '/,
     +  ERMS15 /' ERR: surface-column undersampled '/,
     +  ERMS16 /' ERR: area not fully in image section '/
      DATA
     +  TIP12 /' TIP: reduce cursor area '/,
     +  TIP15 /' TIP: enlarge cursor area, reduce degree and #iter.'/
C
      DATA    OUTPUT   /' '/, TABLE  /' '/,  FRAMEA  /' '/, FRAMEB /' '/
      DATA    TABLAB   /'XSTART ','YSTART ','XEND ','YEND '/
      DATA    KCOOR    /-1,-1,-1,-1/
      DATA    LOW1     /1,1,1/,   HI /1,1,1/
      DATA    NPIXA    /1,1,1/,   WPIX   /1,1,1/
      DATA    DAZOP  /1/,      COOFF /0/
      DATA    CUNITA   /' '/,     IDENTA /' '/
      DATA    SPIX     /1,1,0/,   EPIX   /1,1,0/
      DATA    PCUR1    /6*0.0/,   PCUR2  /6*0.0/
C 
      INCLUDE  'MID_INCLUDE:ST_DAT.INC'
C
C  get input option (CURSOR or table) + read keys
      CSYST = 'W'               !coord.system default: world
      LFLAGS(1) = -99
      IXLO = 0
      IXUP = 0
      IYLO = 0
      IYUP = 0
      CALL STKRDC('P1',1,1,80,IAV,OUTPUT,UNI,NUL,STAT)
      N = INDEX(OUTPUT,',')
      IF (N.GT.0) THEN                               !table specified
         TABFLG = 2
         FRAMEA(1:80) = OUTPUT(1:N-1)//' '
         CALL CLNFRA(FRAMEA,FRAMEA,0)
         TABLE(1:80) = OUTPUT(N+1:)//' '
         N2 = INDEX(TABLE,',')
         IF (N2.GT.0) THEN                           !only for MODI/PIX implem.
            CSYST = TABLE(N2+1:N2+2)
            CALL UPCAS(CSYST,CSYST)
            TABLE(1:80) = TABLE(1:N2-1)//' '
            IF (CSYST.EQ.'P') THEN                   !pixel coord. selected
               DO 90, N=1,4
                  N2 = INDEX(TABLAB(N),' ')          !labels have to end
                  TABLAB(N)(1:) = TABLAB(N)(1:N2-1)//'PIX '
90             CONTINUE                              !with PIX
            ELSE IF (CSYST.NE.'W') THEN              !Wrong flag used
               OUTPUT(1:) = 'Wrong table flag used. Only PIXEL or'//
     +              ' WORLD is possible '
               CALL STETER(3,OUTPUT)
            ENDIF
         ELSE
            CSYST = 'W'                              !world coordinates
         ENDIF
         CALL CLNTAB(TABLE,TABLE,0)
         CALL TBTOPN(TABLE,F_I_MODE,TID,STAT)
         CALL TBIGET(TID,NCOL,NROW,N,N,N,STAT)       !get no. of columns + rows
C 
         DO 100, N=1,4
            CALL TBLSER(TID,TABLAB(N),TABCOLN(N),STAT)      !find columns...
            IF (TABCOLN(N).LE.0) THEN
               OUTPUT(1:) = 'column labelled '//TABLAB(N)//
     +                      'missing... '
               CALL STETER(1,OUTPUT)
            ENDIF
100      CONTINUE
      ELSE
C 
         TABFLG = 1
         CALL STKRDC('IN_A',1,1,80,IAV,FRAMEA,UNI,NUL,STAT)
      ENDIF
C      
C  get result frame + integer parameters
      CALL STKRDC('OUT_A',1,1,80,IAV,FRAMEB,UNI,NUL,STAT)
      CALL STKRDI('INPUTI',1,6,IAV,INPUTI,UNI,NUL,STAT)
      CALL STKRDR('INPUTR',1,3,IAV,INPUTR,UNI,NUL,STAT)
C
C get remaining parameters depending on method
      IF (MODMOD.EQ.1) THEN                          !get the noise flag
                                                     !(RNH)
         CALL STKRDC('NOISE',1,1,1,IAV,CNOISE,UNI,NUL,STAT)
         IF (CNOISE.EQ.'Y' .OR. CNOISE.EQ.'y') THEN
            NOISE=.TRUE.
         ELSE
            NOISE=.FALSE.
         ENDIF
                                                     !get allowed image section
         CALL STKRDC('IMASEC',1,1,100,IAV,CIMASEC,UNI,NUL,STAT)
C
C  verify area factors, degree of polynomial + no. of iterations
         INPUTR(1) = MAX(1.0,INPUTR(1))              !arfacts specified by
         INPUTR(2) = MAX(1.0,INPUTR(2))              !... MODI/PIX
         INPUTI(1) = MAX(1,INPUTI(1))                !only used by MODI/AREA
         INPUTI(2) = MAX(1,INPUTI(2))                !only used by MODI/AREA
         INPUTI(3) = MIN(10,MAX(0,INPUTI(3)))
         INPUTI(4) = MIN(10,MAX(0,INPUTI(4)))
         INPUTI(5) = MAX(0,INPUTI(5))
         IDRAW = INPUTI(6)
      ELSE
C
C  verify degree (only relevant for MODI/AREA)
         DEGREE = INPUTI(1)
         IF (DEGREE.LT.0) THEN
            DEGREE = 0
         ELSE IF (DEGREE.GT.2) THEN
            DEGREE = 2
         ENDIF
         IDRAW = INPUTI(2)
         IF (DEGREE.EQ.0)                   !for degree 0 get constant value
     +        CONST = INPUTR(1)             !changed by swolf@eso.org
c     +      CALL STKRDR('INPUTR',1,1,IAV,CONST,UNI,NUL,STAT)
      ENDIF
C
C  map input + result frame

      IF (FRAMEB(1:1).NE.'+') THEN
         CALL GENEQF(FRAMEA,FRAMEB,STAT) 
      ELSE
         FRAMEB(1:) = FRAMEA(1:)                  !result frame = input frame
         STAT = 1
      ENDIF
      IF (STAT.EQ.1) THEN
         CALL STIGET(FRAMEA,D_R4_FORMAT,F_IO_MODE,F_IMA_TYPE,
     +               2,NAXIS,NPIXA,STARTA,STEPA,IDENTA,
     +               CUNITA,PNTRA,IMNOA,STAT)
         IMNOB = IMNOA
         PNTRB = PNTRA
         SIZE = NPIXA(1)*NPIXA(2)
         IF (NAXIS.GT.2) THEN
            NAXIS = 2
            CALL STTPUT('Only 1. plane of frame used!',STAT)
         ENDIF
         CALL STTPUT
     +   ('Warning: Input frame will be updated...',STAT)
      ELSE
         CALL STIGET(FRAMEA,D_R4_FORMAT,F_I_MODE,F_IMA_TYPE,
     +               2,NAXIS,NPIXA,STARTA,STEPA,IDENTA,
     +               CUNITA,PNTRA,IMNOA,STAT)
         IF (NAXIS.GT.2) THEN
            NAXIS = 2
            CALL STTPUT('Only 1. plane of input frame used!',STAT)
         ENDIF
         CALL STIPUT(FRAMEB,D_R4_FORMAT,F_O_MODE,F_IMA_TYPE,
     +               NAXIS,NPIXA,STARTA,STEPA,IDENTA,
     +               CUNITA,PNTRB,IMNOB,STAT)
         SIZE = NPIXA(1)*NPIXA(2)
         OUTPUT(1:) = '... copying input frame to result frame ... '
         CALL STTDIS(OUTPUT,99,STAT)
         CALL COPYF(MADRID(PNTRA),MADRID(PNTRB),SIZE)
         IF (IDRAW .GE. 2) THEN                !copy to disk
            IPIXS(1) = 1                       !IPIXS will be reset later on...
            IPIXS(2) = 1
            CALL MEMDSK(IMNOB,0,MADRID(PNTRB),NPIXA,
     +                  IPIXS,NPIXA,STAT)
         ENDIF
      ENDIF
      LFLAGS(2) = D_R4_FORMAT                   !indicate real data type
C
C  get cuts of input image 
      CALL STDRDR(IMNOA,'LHCUTS',1,4,IAV,CUTS,UNI,NUL,STAT)
      IF (CUTS(1).GE.CUTS(2)) THEN
         CUTIND = 3
      ELSE
         CUTIND = 1
      ENDIF
      ZCUTS(1) = CUTS(3)
      ZCUTS(2) = CUTS(4)
      STPINV(1) = 1.0 / STEPA(1)
      STPINV(2) = 1.0 / STEPA(2)
C 
      IF (MODMOD.EQ.1) THEN
         IF (CIMASEC(1:1).NE.'+') THEN
            CALL EXTCOO(IMNOA,CIMASEC,3,IAV,SPIX,EPIX,STAT)
            IF (STAT.NE.0) 
     +         CALL STETER(2,'Invalid coords for image section...')
C 
         ELSE
            SPIX(1) = 1
            EPIX(1) = NPIXA(1)
            SPIX(2) = 1
            EPIX(2) = NPIXA(2)
         ENDIF
         IXLO = SPIX(1)         !lower x-coordinate
         IXUP = EPIX(1)         !upper x-coordinate
         IYLO = SPIX(2)         !lower y-coordinate
         IYUP = EPIX(2)         !upper y-coordinate of allowed image section
         OMIN = CUTS(3)
         IF (OMIN.LE.0.) THEN                    !set offset value...
            OFFSET = ABS(OMIN) + 1.
         ELSE
            OFFSET = 0.
         ENDIF
      ENDIF
C  
C  fork according to input mode
      NMAL = 0
      IF (TABFLG.EQ.2) THEN
         IF (IDRAW.GE.2) THEN
            CALL DTOPEN(1,STAT)                  !check, if image is loaded
            CALL DTGICH(QDSPNO,QIMCH,OUTPUT,INPUTR,STAT)
            CALL GENEQF(FRAMEA,OUTPUT,STAT) 
            IF (STAT.NE.1) THEN
               CALL STTPUT
     +         ('Input image not lodaded => draw flag ignored...',
     +          STAT)
               IDRAW = 0
            ELSE
               CALL PIXXCV('INIT',IMNOA,RBUF,STAT)
            ENDIF
         ELSE
            DAZOP = 0                      !indicate, that we don't use display
         ENDIF
         IF (CSYST.EQ.'W')           !init wc --> frame pixel conversion
     +      CALL FPXWCO(0,IMNOA,DD1,DD2,STAT)
         GOTO 3000
      ENDIF
C
      IF ((IDRAW.EQ.1).OR.(IDRAW.EQ.3)) THEN
         CACT = 'NYYY?C0 '
      ELSE
         CACT = 'NNYY?C0 '
      ENDIF
C
C  for interactive (cursor) mode - attach ImageDisplay
      CALL DTOPEN(1,STAT)
      FRAMEA(1:) = ' '                         !to init GETCUR ...
      CALL SETCUR(QDSPNO,2,0,0,KCOOR,STAT)     !set cursor rectangle
      CALL CONCHA(QDSPNO,QOVCH,1,0)            !clear overlay channel
      CALL PIXXCV('INIT',IMNOA,RBUF,STAT)
C
C  input cursor positions and check status
1000  WRITE(OUTPUT,10000) NMAL+1
      CALL STTDIS(OUTPUT,99,STAT)
      CALL GETCUR(CACT,FRAMEA,
     +            ICUR1,PCUR1,COOR(1),RBUF(1),STAT1,
     +            ICUR2,PCUR2,COOR(3),RBUF(1),STAT2)
C
C  check cursor status  
      IF ((STAT1.EQ.0).AND.(STAT2.EQ.0)) THEN
         IF ((NMAL.EQ.0) .AND. (COOFF.EQ.0)) THEN
            COOFF = 1
            FRAMEA(1:) = ' '
            CALL STTPUT(INFO,STAT)
            GOTO 1000
         ELSE
            GOTO 5000
         ENDIF
      ENDIF
C      
C  now goto section, where we really do the work
      NMAL = NMAL + 1                        !update counter
      IF (MODMOD.EQ.1) THEN
         PCEN(1) = (ICUR1(1)+ICUR2(1)) * .5
         PCEN(2) = (ICUR1(2)+ICUR2(2)) * .5
         CALL PIXXCV('IRW',0,PCEN,STAT)
      ENDIF
      GOTO 4000
C      
C  here we get our input data from a table
C      
3000  NMAL = NMAL + 1
      IF (NMAL.GT.NROW) THEN
         CALL TBTCLO(TID,STAT)
         WRITE(OUTPUT,10001) (NMAL-1)
         CALL STTPUT(OUTPUT,STAT)
         GOTO 5000
      ENDIF
C      
C  read xstart, ystart, xend, yend from table
      CALL TBSGET(TID,NMAL,SELFLG,STAT)
      IF (.NOT.SELFLG) GOTO 3000
C 
      CALL TBRRDR(TID,NMAL,4,TABCOLN,COOR,TABNUL,STAT)
      IF (CSYST.EQ.'W') THEN          !world coordinates --> transfer to pixels
         DD1(1) = COOR(1)                !convert wc -> fp
         DD1(2) = COOR(2)
         CALL FPXWCO(-1,0,DD1,DD2,STAT)
         IF (STAT .GT. 0)
     +      CALL STETER(36,'Problems with WCS of frame...')
         PCUR1(1) = DD2(1)
         PCUR1(2) = DD2(2)
         DD1(1) = COOR(3)
         DD1(2) = COOR(4)
         CALL FPXWCO(-1,0,DD1,DD2,STAT)
         IF (STAT .GT. 0)
     +      CALL STETER(36,'Problems with WCS of frame...')
         PCUR2(1) = DD2(1)
         PCUR2(2) = DD2(2)
      ELSE                            !table entries area already in pixels
         PCUR1(1) = COOR(1)
         PCUR1(2) = COOR(2)
         PCUR2(1) = COOR(3)
         PCUR2(2) = COOR(4)
      END IF
      IF (IDRAW.GE.2) CALL PIXXCV('_RS',0,PCUR1,STAT)
      ICUR1(1) = NINT(PCUR1(5))
      ICUR1(2) = NINT(PCUR1(6))
C 
C  do the replacement
4000  IF (MODMOD.EQ.1) THEN           !define cursor area + surrounding surface 
         IPIXS(1) = NINT(PCUR1(1))    !real pixel no.s
         IPIXS(2) = NINT(PCUR2(1))
         IPIXS(3) = NINT(PCUR1(2))
         IPIXS(4) = NINT(PCUR2(2))
         IF ((IPIXS(2).LT.IXLO) .OR. (IPIXS(1).GT.IXUP) .OR.
     +       (IPIXS(4).LT.IYLO) .OR. (IPIXS(3).GT.IYUP)) THEN
            CALL STTPUT(ERMS16,STAT)
            GOTO (1000,3000),TABFLG
         ENDIF
         IPX = IPIXS(2) - IPIXS(1) + 1        !no. of pixels in cursor area...
         IPY = IPIXS(4) - IPIXS(3) + 1
         ICX = IPIXS(1) + IPX/2               !center of cursor area
         ICY = IPIXS(3) + IPY/2
C      
         ISX = NINT(INPUTR(1)*(IPX-1)) + 1     !no. of pixels in surface
         ISY = NINT(INPUTR(2)*(IPY-1)) + 1 
         ISUR(1) = MAX(IXLO,ICX - ISX/2)       !xlow
         ISUR(2) = MIN(IXUP,ISUR(1) + ISX - 1) !xhi
         ISUR(3) = MAX(IYLO,ICY - ISY/2)       !ylow
         ISUR(4) = MIN(IYUP,ISUR(3) + ISY - 1) !yhi
C
         XDEG = INPUTI(3)                     !erase spot
         YDEG = INPUTI(4)
         IT = INPUTI(5)
         CRMS = -1.
C 
         CALL ERASE(MADRID(PNTRB),NPIXA(1),NPIXA(2),IPIXS,ISUR,XDEG,
     +              YDEG,OFFSET,IT,CRMS,NOISE,STAT)
C
         IF (STAT.NE.0) THEN                        !check return status
            IF (STAT.GE.1 .AND. STAT.LE.5) THEN     !NAGLIB error
               WRITE(OUTPUT,111) STAT
               CALL STTPUT(OUTPUT,STAT)
               CALL STTPUT(ERMS13,STAT)
               GOTO 5500
            ELSE IF (STAT .EQ. 7) THEN             !surface dimension too large
               CALL STTPUT(ERMS12,STAT)
               CALL STTPUT(TIP12,STAT)
               GOTO (1000,3000),TABFLG
            ELSE IF (STAT.EQ.8) THEN                !subimage-row undersampled
               CALL STTPUT(ERMS15,STAT)
               CALL STTPUT(TIP15,STAT)
               GOTO (1000,3000),TABFLG
            ELSE IF (STAT.EQ.9) THEN                !subimage-not well defined
               CALL STTPUT(ERMS16,STAT)
               GOTO (1000,3000),TABFLG
            ENDIF
         ENDIF
C 
         PCEN(5) = (COOR(1)+COOR(3)) * .5
         PCEN(6) = (COOR(2)+COOR(4)) * .5
         WRITE(OUTPUT,222) IPX,IPY,PCEN(5),PCEN(6)
         CALL STTPUT(OUTPUT,STAT)
C
      ELSE
         CALL STAREM(MADRID(PNTRB),NPIXA(1),NPIXA(2),
     +               STARTA,STEPA,COOR,CSYST,DEGREE,CONST,RCUTS)
         IF (ZCUTS(1).GT.RCUTS(1)) ZCUTS(1) = RCUTS(1)
         IF (ZCUTS(2).LT.RCUTS(2)) ZCUTS(2) = RCUTS(2)
      ENDIF
C
C  if required, reload modified subwindow
      IF (IDRAW.GE.2) THEN
         LOW1(1) = NINT(PCUR1(1))
         LOW1(2) = NINT(PCUR1(2))
         HI(1) = NINT(PCUR2(1))
         HI(2) = NINT(PCUR2(2))
         WPIX(1) = HI(1) - LOW1(1) + 1
         WPIX(2) = HI(2) - LOW1(2) + 1
         CALL MEMDSK(IMNOB,0,MADRID(PNTRB),NPIXA,LOW1,WPIX,STAT)
         WW(1) = ICUR1(1)
         WW(2) = ICUR1(2)
         WW(3) = SCALX
         WW(4) = SCALY
         CALL LOADWN(LFLAGS,IMNOB,NPIXA,LOW1,WPIX,WW,CUTS(CUTIND))
      ENDIF
C
      IF ((IDRAW.EQ.1) .OR. (IDRAW.EQ.3)) THEN
         WW(1) = ICUR1(1)
         WW(2) = ICUR1(2)
         WPIX(1) = HI(1) - LOW1(1) + 1
         WPIX(2) = HI(2) - LOW1(2) + 1
         IF (SCALX.GT.1) THEN
            WPIX(1) = WPIX(1)/SCALX
         ELSE IF (SCALX.LT.-1) THEN
            N2 = -SCALX
            WPIX(1) = WPIX(1)*N2
         ENDIF
         IF (SCALY.GT.1) THEN
            WPIX(2) = WPIX(2)/SCALY
         ELSE IF (SCALY.LT.-1) THEN
            N2 = -SCALY
            WPIX(2) = WPIX(2)*N2
         ENDIF
         WW(3) = WW(1) + WPIX(1) - 1
         WW(4) = WW(2)+ WPIX(2) - 1
         CALL BLDGRA('REC',WW,INPUTR,XFIGU,YFIGU,6,N2)
         CALL IIGPLY(QDSPNO,QOVCH,XFIGU,YFIGU,N2,255,1,STAT)
      ENDIF
C 
      IF (TABFLG.EQ.1) THEN
         GOTO 1000
      ELSE
         GOTO 3000
      ENDIF
C
C  get new dynamic range + write descriptor LHCUTS
5000  IF (CUTS(3).GT.ZCUTS(1)) CUTS(3) = ZCUTS(1)
      IF (CUTS(4).LT.ZCUTS(2)) CUTS(4) = ZCUTS(2)
C      
C  append HISTORY
      CALL DSCUPT(IMNOA,IMNOB,' ',STAT)
      CALL STDWRR(IMNOB,'LHCUTS',CUTS,1,4,UNI,STAT)      !but use old cuts...
C      
C  That's it folks...
5500  IF (DAZOP.EQ.1) THEN
         CALL DTCLOS(QDSPNO)
         IF (IDRAW.GT.0) CALL REFOVR(STAT)
      ENDIF
      OUTPUT(1:) = ' '
      CALL STTDIS(OUTPUT,99,STAT)
      RETURN
C      
C  format statements
10000 FORMAT(I3,', ready for cursor input...')
10001 FORMAT(I6,' rows processed.')
111   FORMAT('IFAIL = ',I3)
222   FORMAT('cleaned area:',I3,' x',I3,' pixels at xc = ',F10.3,
     +       ', yc = ',F10.3)
      END

      SUBROUTINE 
     + STAREM(A,NPIX1,NPIX2,STARTA,STEPA,COOR,CSYST,IDEG,CONST,RMAX)
C
      IMPLICIT NONE
C      
      INTEGER  NPIX1,NPIX2,IDEG
      INTEGER  I,I1,I2,J,J1,J2
      INTEGER  M,N,INC
      INTEGER  STAT
C 
      REAL     A(NPIX1,NPIX2),COOR(4),RMAX(2)
      REAL     X(200),Y(200),Z(200),B(7)
      REAL     AIJ,B1,B2,B3,B4,B5,B6,D2
      REAL     Q2,Q22,QX,QX2,QY,QY2
      REAL     X1,X2,XC,XX,Y1,Y2,Y22,YC,YY
      REAL     ZMAX,ZMIN
      REAL     STEP1,STEP2,START1,START2
      REAL     CONST,STEP11,STEP22
C 
      DOUBLE PRECISION STARTA(*),STEPA(*)
C
      CHARACTER CSYST*1, OUTPUT*80
C
222   FORMAT('clean area:',I3,' x',I3,' pixels at xc = ',F10.3,
     +       ', yc = ',F10.3)
C      
C  init
      STEP1=STEPA(1)
      STEP2=STEPA(2)
      START1=STARTA(1)
      START2=STARTA(2)
C      
C  compute
      X1=COOR(1)
      Y1=COOR(2)
      X2=COOR(3)
      Y2=COOR(4)
C      
      XC=(X1+X2)*0.5
      YC=(Y1+Y2)*0.5
      QX=X2-XC
      QY=Y2-YC
      QX2=QX**2
      QY2=QY**2
      Q2=QX*QY
      Q22=Q2**2
C      
C  define outer rectangular boundary
      IF (CSYST.EQ.'W') THEN
         J1=(Y1-START2)/STEP2+1.
         J2=(Y2-START2)/STEP2+1.
         IF (J1.LT.1) THEN
            J1=1
         ELSE IF (J2.GT.NPIX2) THEN
            J2=NPIX2
         ENDIF
C 
         I1=(X1-START1)/STEP1+1.
         I2=(X2-START1)/STEP1+1.
         IF (I1.LT.1) THEN
            I1=1
         ELSE IF (I2.GT.NPIX1) THEN
            I2=NPIX1
         ENDIF
         Y1 = FLOAT(J1-1)*STEP2+START2-YC
         X1 = FLOAT(I1-1)*STEP1+START1-XC
      ELSE
         J1=NINT(Y1)
         J2=NINT(Y2)
         I1=NINT(X1)
         I2=NINT(X2)
         IF (J1.LT.1) THEN
            J1=1
         ELSE IF (J2.GT.NPIX2) THEN
            J2=NPIX2
         ENDIF
         IF (I1.LT.1) THEN
            I1=1
         ELSE IF (I2.GT.NPIX1) THEN
            I2=NPIX1
         ENDIF
         Y1 = Y1 - YC
         X1 = X1 - XC
         STEP1 = 1.0
         STEP2 = 1.0
      ENDIF
C      
C  assemble points between outer ellipse and circumscribed rectangle
C  XX and YY are coordinates relative to the centre (XC,YC)
C 
      YY = Y1
C
      WRITE(OUTPUT,222) I2-I1+1,J2-J1+1,XC,YC
      CALL STTPUT(OUTPUT,STAT)
C      
C  if degree = 0, replace intensities inside ellipse by constant
      IF (IDEG.EQ.0) THEN
         DO 400, J=J1,J2
            Y22=YY*YY*QX2
            XX=X1
            DO 380, I=I1,I2
               D2 = (XX*XX*QY2) + Y22
               IF (D2.LE.Q22) THEN
                  A(I,J)=CONST
               ENDIF
               XX=XX+STEP1
380         CONTINUE
            YY=YY+STEP2
400      CONTINUE
C 
         RMAX(1) = CONST
         RMAX(2) = CONST
         RETURN                        !we're done...
      ENDIF
C 
C  compute approx max no of points (the exact no depends on ellipse centre) 
C  if necessary take coarser sample
C 
      M = (I2-I1+1)*(J2-J1+1) - NINT(3.14159*Q2/STEP1**2)      !no. of points
      INC = 1 + (M/200)
      STEP11 = STEP1*FLOAT(INC)
      STEP22 = STEP2*FLOAT(INC)
      M = 0
      ZMIN = 1.E10
      ZMAX = -ZMIN
      DO 2000, J=J1,J2,INC
         Y22=YY*YY*QX2
         XX=X1
         DO 1000, I=I1,I2,INC
            D2 = (XX*XX*QY2) + Y22
            IF(D2.LT.Q22) GOTO 990
            M=M+1
            X(M)=XX
            Y(M)=YY
            AIJ=A(I,J)
            Z(M)=AIJ
            ZMIN=AMIN1(ZMIN,AIJ)
            ZMAX=AMAX1(ZMAX,AIJ)
            IF(M.EQ.200) GOTO 3000
990         XX=XX+STEP11
1000     CONTINUE
         YY=YY+STEP22
2000  CONTINUE
C      
C  test M
      RMAX(1) = ZMIN
      RMAX(2) = ZMAX
      IF (M.LT.4) RETURN
C      
C  compute N
      IF ( (IDEG.EQ.1) .OR.
     +     ((ZMAX-ZMIN).LT.(ZMAX+ZMIN)*0.2) ) THEN
         N=3
      ELSE
         N=6
      ENDIF
C
C  fit max 2nd order SURFAC to M points 
C  if M < 12, then use a plane only (degree = 1)
C  the 4 corner points should always be included
C  otherwise skip replacement
C
3000  CALL SURFAC(X,Y,Z,B,M,N)
      B1 = B(1)
      IF (N.NE.3.AND.(B1.LT.ZMIN.OR.B1.GT.ZMAX)) THEN
         N=3
         GOTO 3000
      ENDIF
C      
      B2=B(2)
      B3=B(3)
      B4=B(4)
      B5=B(5)
      B6=B(6)
C      
C  replace intensities inside ellipse
      YY=Y1
      DO 4000, J=J1,J2
         Y22=YY*YY*QX2
         XX=X1
         DO 3800, I=I1,I2
            D2 = (XX*XX*QY2) + Y22
            IF (D2.LE.Q22) THEN
               A(I,J)=B1 + XX*(B2+XX*B4+YY*B5) + YY*(B3+YY*B6)
               IF (RMAX(1) .GT. A(I,J)) THEN
                  RMAX(1) = A(I,J)
               ELSE IF (RMAX(2) .LT. A(I,J)) THEN
                  RMAX(2) = A(I,J)
               ENDIF
            ENDIF
            XX=XX+STEP1
3800     CONTINUE
         YY=YY+STEP2
4000  CONTINUE
C      
      RETURN
      END

      SUBROUTINE SURFAC(X,Y,Z,B,M,N)
C 
      IMPLICIT NONE
C 
      INTEGER   M,N
      INTEGER   N1,K,I,J
C 
      DOUBLE PRECISION  A(7,7),C(7)
C 
      REAL   X(M),Y(M),Z(M),B(7)
C  
C      LSQ fitting of surface
C      Z = B1 + B2*X + B3*Y + B4*X*X + B5*X*Y + B6*Y*Y
C      to M points
C
      IF (M.LT.12) N=3
      N1 = N+1
      DO 2, I=1,7
         DO 1, J=1,7
            A(I,J)=0.D0
1        CONTINUE
         B(I)=0.
2     CONTINUE
C 
      DO 5, K=1,M
         C(1)=1.D0
         C(2)=DBLE(X(K))                
         C(3)=DBLE(Y(K))
         C(N1)=DBLE(Z(K))
         IF (N.EQ.6) THEN
            C(4)=C(2)**2
            C(5)=C(2)*C(3)
            C(6)=C(3)**2
         ENDIF
C 
         DO 4, I=1,N1
            DO 3, J=1,N1
               A(I,J)=A(I,J)+C(I)*C(J)
3           CONTINUE
4        CONTINUE     
5     CONTINUE
C 
      CALL SIMUL(N,A,B)
C 
      RETURN
      END   

      SUBROUTINE SIMUL(N,A,X)
C 
C   the GAUSS-JORDAN complete elimination method is employed
C   with the maximum pivot strategy
C 
      IMPLICIT NONE
C 
      INTEGER    N,I,J,K,IROWI,IROWK,JCOLI,JCOLK
      INTEGER    ISCAN,JSCAN,KM1,MAX 
      INTEGER    IROW(6),JCOL(6),STAT
C 
      DOUBLE PRECISION  A(7,7),EPS,AIJCK,PIVOT,FPIVOT
C 
      REAL       X(N)
C 
      DATA   EPS  /1.D-15/
      DATA   IROW /6*0/, JCOL /6*0/
C 
      MAX=N+1   
C 
C  begin elimination procedure   
C 
      DO 19,K=1,N                 
         KM1=K-1                     
         PIVOT=0.D0                    !search for the pivot element  
         DO 10,I=1,N                 
            DO 11,J=1,N                 
C 
C search irow and jcol arrays for invalid pivot subscripts
               IF (K.EQ.1) GOTO 9            
               DO 8,ISCAN=1,KM1            
                  DO 7,JSCAN=1,KM1            
                     IF (I.EQ.IROW(ISCAN)) GOTO 11 
                     IF (J.EQ.JCOL(JSCAN)) GOTO 11 
7                 CONTINUE                    
8              CONTINUE                    
9              CONTINUE                    
               IF (DABS(A(I,J)).LE.DABS(PIVOT)) GOTO 11
               PIVOT=A(I,J)
               IROW(K)=I   
               JCOL(K)=J   
11          CONTINUE    
10       CONTINUE    
C 
C  insure that selected pivot is larger than eps
         IF (DABS(PIVOT).LE.EPS) THEN
            CALL STTPUT(' Warning! Singular Matrix',STAT)
            RETURN                      
         ELSE
            FPIVOT = 1.D0 / PIVOT
         ENDIF
C 
         IROWK=IROW(K)               
         JCOLK=JCOL(K)               
C 
C   normalize pivot row elements  
         DO 14,J=1,MAX               
            A(IROWK,J)=A(IROWK,J)*FPIVOT 
14       CONTINUE
C 
C carry out elimination   
C 
         A(IROWK,JCOLK)=FPIVOT     
         DO 18,I=1,N                 
            AIJCK=A(I,JCOLK)            
            IF (I.EQ.IROWK) GOTO 18       
            A(I,JCOLK) = -(AIJCK*FPIVOT)
            DO 17,J=1,MAX               
               IF (J.NE.JCOLK) A(I,J)=A(I,J)-AIJCK*A(IROWK,J) 
17          CONTINUE              
18       CONTINUE
19    CONTINUE   
C 
C order solution values   
C 
      DO 20,I=1,N           
         IROWI=IROW(I)         
         JCOLI=JCOL(I)         
         X(JCOLI)= A(IROWI,MAX)
20    CONTINUE
C 
      RETURN
      END

      SUBROUTINE MODCUT(ACTIO)
C
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C.IDENTIFICATION              
C  subroutine MODCUT		version 1.00	890119
C  K. Banse			ESO - Garching
C
C.KEYWORDS    
C  bulk data frame, scaling
C
C.PURPOSE
C  to rescale subimages on the ImageDisplay
C
C.ALGORITM
C  extract subimage via cursor window + apply new cut values
C
C.INPUT/OUTPUT
C  ACTIO	char. array	1 = C - cursor option, else full image
C 
C  the following keywords are used:
C  INPUTI/I/1/4                 no_of_cursors,draw_flag 
C                               and xsiz,ysiz if single cursor
C 
C  last used cut values are written to
C  OUTPUTR/R/1/2
C
C-------------------------------------------------------------------
C
      IMPLICIT NONE
C      
      INTEGER      IMNO,IAV,N,NMAL,NAXPIX(4),KEY
      INTEGER      STAT1,STAT2,STAT,LFLAGS(2)
      INTEGER      NPIXA(3),KCOOR(4),WPIX(3)
      INTEGER      LOW1(3),HI(3),IFIRST
      INTEGER      BOX(4),COOFF,UNI(1),NULO
C      
      CHARACTER    ACTIO*(*)
      CHARACTER    FRAMEA*100
      CHARACTER    INFO*44,AREA*80,FRMSTR*8,OPTION*5,CACT*8
C      
      DOUBLE PRECISION STARTA(3),STEPA(3)
C 
      REAL         RR,RVAL1,RVAL2,COOR(4),RCUTS(2)
      REAL         PCUR(4),ZBINS(3),RBUF(8),ROLD(4)
C 
      INCLUDE  'MID_INCLUDE:ST_DEF.INC'
      INCLUDE  'MID_INCLUDE:IDIDEV.INC'
      INCLUDE  'MID_INCLUDE:IDIMEM.INC'
C 
      DATA     INFO
     +    /' switch cursors on - next time we exit ... '/
C
      DATA    FRAMEA   /'         '/
      DATA    KCOOR    /-1,-1,-1,-1/
      DATA    LOW1     /1,1,1/,     HI   /1,1,1/
      DATA    NPIXA    /1,1,1/,     WPIX /1,1,1/
      DATA    COOFF    /0/,         IFIRST /0/
      DATA    FRMSTR   /'E15.6'/
      DATA    ZBINS    /256.0,0.0,0.0/
C 
      INCLUDE  'MID_INCLUDE:ST_DAT.INC'
C      
C  attach Image Display + get loaded image
      CALL DTOPEN(1,STAT)
      CALL DTGICH(QDSPNO,QIMCH,FRAMEA,RBUF,STAT)
C 
      CALL STFOPN(FRAMEA,D_R4_FORMAT,0,F_IMA_TYPE,IMNO,STAT)
      CALL STDRDI(IMNO,'NAXIS',1,1,IAV,NAXPIX,UNI,NULO,STAT)
      CALL STDRDI(IMNO,'NPIX',1,NAXPIX(1),IAV,NPIXA,UNI,NULO,STAT)
      CALL STDRDD(IMNO,'START',1,NAXPIX(1),IAV,STARTA,UNI,NULO,STAT)
      CALL STDRDD(IMNO,'STEP',1,NAXPIX(1),IAV,STEPA,UNI,NULO,STAT)
      DO 100, N=1,3
         NAXPIX(N+1) = NPIXA(N)
100   CONTINUE
C
      LFLAGS(2) = D_R4_FORMAT     !tell LOADWN that we opened as R4 data
      IF (ACTIO(1:1).EQ.'C') GOTO 490
C 
C ---------------------------------------
C 
C  here we handle the full image
C 
C ---------------------------------------
C 
      OPTION = 'YYHNH'            !full frame, histogram only
      BOX(1) = SSPX
      BOX(2) = SSPY
      LOW1(1) = SFPX
      LOW1(2) = SFPY
      WPIX(1) = NSX
      WPIX(2) = NSY
      LFLAGS(1) = -100            !whole frame in display window
C 
400   AREA(1:) = ' '
      WRITE(AREA,10001) COOR(1),COOR(2),COOR(3),COOR(4)
      CALL BLANKO(AREA)           !make statistics in displayed area
      CALL ZMSTAT(IMNO,AREA,NAXPIX,STARTA,STEPA,ZBINS,FRMSTR,OPTION)
C
440   CALL PLOHI(IMNO)
      IF (IFIRST.EQ.0) THEN
         RVAL1 = 50
         RVAL2 = 50
         IFIRST = 1
      ELSE
         CALL PTDATA(20,0,0,ROLD(1),ROLD(2),0.0,1)    !show old cuts
         CALL PTDATA(20,0,0,ROLD(3),ROLD(4),0.0,1)
      ENDIF
      CALL PTGCUR(RVAL1,RVAL2,KEY,STAT)               !graphCursor input
      IF (KEY.EQ.32) GOTO 8900                        !gCursor exit
C 
      ROLD(1) = RVAL1
      ROLD(2) = RVAL2
      CALL PTDATA(4,0,0,RVAL1,RVAL2,0.0,1)
      RCUTS(1) = RVAL1
      CALL PTGCUR(RVAL1,RVAL2,KEY,STAT)
      IF (KEY.EQ.32) GOTO 8900                        !gCursor exit
C 
      ROLD(3) = RVAL1
      ROLD(4) = RVAL2
      CALL PTDATA(4,0,0,RVAL1,RVAL2,0.0,1)
      RCUTS(2) = RVAL1
      IF (RCUTS(1).GT.RCUTS(2)) THEN
         RR = RCUTS(1)
         RCUTS(1) = RCUTS(2)
         RCUTS(2) = RR
      ENDIF
      WRITE(AREA,10010) RCUTS(1),RCUTS(2)
      CALL STTDIS(AREA,99,STAT)
      CALL STDWRR(IMNO,'LHCUTS',RCUTS,1,2,UNI,STAT)
C
C  load image with new cuts
C
      CALL LOADWN(LFLAGS,IMNO,NPIXA,LOW1,WPIX,BOX,RCUTS)
      CALL DAZVIS(QDSPNO,QIMCH,2,1)
      GOTO 440
C 
C ---------------------------------------
C 
C  here we handle cursor defined sub-frames
C 
C ---------------------------------------
C 
490   LFLAGS(1) = -99             !subframes in display window
      CACT = 'YNYY?C0 '
      CALL SETCUR(QDSPNO,2,0,2,KCOOR,STAT)     !set cursor shape
      FRAMEA(1:) = ' '                         !to init GETCUR ...
      NMAL = -1
      OPTION = 'NYHNH'            !subframe, histogram only
C
C  input cursor positions and check status
500   NMAL = NMAL + 1
550   CALL GETCUR(CACT,FRAMEA,
     +            BOX(1),PCUR(1),COOR(1),RVAL1,STAT1,
     +            BOX(3),PCUR(3),COOR(3),RVAL2,STAT2)
C
C  check cursor status  
      IF ((STAT1.EQ.0).AND.(STAT2.EQ.0)) THEN
         IF ((NMAL.EQ.0) .AND. (COOFF.EQ.0)) THEN
            COOFF = 1
            FRAMEA(1:) = ' '
            CALL STTPUT(INFO,STAT)
            GOTO 550
         ELSE
            GOTO 9000                          !get out
         ENDIF
      ENDIF
C
      AREA(1:) = ' '
      WRITE(AREA,10001) COOR(1),COOR(2),COOR(3),COOR(4)
      CALL BLANKO(AREA)
      CALL ZMSTAT(IMNO,AREA,NAXPIX,STARTA,STEPA,ZBINS,FRMSTR,OPTION)
C 
700   CALL PLOHI(IMNO)
      IF (IFIRST.EQ.0) THEN
         RVAL1 = 50
         RVAL2 = 50
         IFIRST = 1
      ELSE
         CALL PTDATA(20,0,0,ROLD(1),ROLD(2),0.0,1)    !show old cuts
         CALL PTDATA(20,0,0,ROLD(3),ROLD(4),0.0,1)
      ENDIF
      CALL PTGCUR(RVAL1,RVAL2,KEY,STAT)
      IF (KEY.EQ.32) GOTO 1000                        !move to next subframe
C 
      ROLD(1) = RVAL1
      ROLD(2) = RVAL2
      CALL PTDATA(4,0,0,RVAL1,RVAL2,0.0,1)
      RCUTS(1) = RVAL1
      CALL PTGCUR(RVAL1,RVAL2,KEY,STAT)
      IF (KEY.EQ.32) GOTO 1000                        !move to next subframe
C 
      ROLD(3) = RVAL1
      ROLD(4) = RVAL2
      CALL PTDATA(4,0,0,RVAL1,RVAL2,0.0,1)
      RCUTS(2) = RVAL1
      IF (RCUTS(1).GT.RCUTS(2)) THEN
         RR = RCUTS(1)
         RCUTS(1) = RCUTS(2)
         RCUTS(2) = RR
      ENDIF
      WRITE(AREA,10010) RCUTS(1),RCUTS(2)
      CALL STTDIS(AREA,99,STAT)
C 
C  load sub-image with new cuts
C 
      LOW1(1) = NINT(PCUR(1))
      LOW1(2) = NINT(PCUR(2))
      HI(1) = NINT(PCUR(3))
      HI(2) = NINT(PCUR(4))
      WPIX(1) = HI(1) - LOW1(1) + 1
      WPIX(2) = HI(2) - LOW1(2) + 1
      CALL LOADWN(LFLAGS,IMNO,NPIXA,LOW1,WPIX,BOX,RCUTS)
      GOTO 700
C 
1000  CALL STTPUT(AREA,STAT)
      CALL STTDIS
     + ('Move to next subframe                           ',99,STAT)
      IFIRST = 0
      GOTO 500
C      
C  That's it folks...
8900  CALL STTPUT(AREA,STAT)                  !write permanently
      CALL STDWRR(IMNO,'LHCUTS',RCUTS,1,2,UNI,STAT)    !update descr. LHCUTS
C 
9000  CALL STTDIS(' ',-9,STAT)
      CALL DTCLOS(QDSPNO)
      CALL STKWRR('OUTPUTR',RCUTS,1,2,UNI,STAT)         !save last cuts
      RETURN
C      
10001 FORMAT('[',G12.6,',',G12.6,':',G12.6,',',G12.6,'] ')
10010 FORMAT('Cut values used: ',2F12.5)
      END

      SUBROUTINE ERASE(CIMG,NPIX,NLIN,ICUR,IMAP,
     +                 K,L,OFS,IT,CRMS,NOISE,ISTAT)
C
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C.IDENTIFICATION 
C subroutine ERASE     		
C (H.Waldthausen)
C K. Banse			version 1.00	880922
C R. Hook                       version 1.10    890621
C        - Less restrictive area of fit
C        - Optional noise according to new logical flag (NOISE)
C
C.TYPE        subroutine
C
C.KEYWORDS    interpolation
C
C.PURPOSE
C
C	This routine removes an unwanted detail in an image. 
C	The surface is 
C       an approximation to the weighted,least- squares Chebyshev
C	series surface fit.The weights are derived from the Poisson
C	error for the two first modes. 
C
C.ALGORITHM   1) read cursor coordinates
C             2) construct surface area around cursor
C	      3) replace bad pixels inside cursor
C	      4) optionally add noise
C
C.MODULES CALLS   E02CAF,E02CBF  from NAGLIB
C
C.INPUT/OUPUT   parameters
C
C	->CIMG  input image to be updated
C	->NPIX  x-pixels and 
C	->NLIN  y-lines of CIMG
C	->IMAP  surface area indices
C	->ICUR  cursor area indices
C	->K     surface fit order in x
C	->L     surface fit order in y
C	->OFS   add this pixel offset to CIMG 
C	->IT    number of iterations 
C	->CRMS  rms of CIMG  (if pos. take this value)  
C	->NOISE logical flag for whether patch should be filled with noise
C	<-ISTAT output status 
C               : 0 O.K.
C               : POS 1-5 failure of E02CAF or E02CBF routines
C		: 6 fit order exeeded
C		: 7 surface-size larger than 102400 pixels
C		: 8 subimage-row undersampled
C		: 9 subimage-bounds are wrong
C
C CHANGED:
C     28.4.98 swolf@eso.org - blow area up: 60x60 --> 320x320
C     24.9.98 swolf@eso.org - check areas (xlo<xup, ylo<yup)
C------------------------------------------------------------------
C					 	
      IMPLICIT NONE
C 
      INTEGER   KMAX,MAXLIN,MAXPIX,MTOTAL,NA,NW
      PARAMETER (KMAX = 10)
      PARAMETER (MAXLIN = 4096)           !max 4096 lines
      PARAMETER (MAXPIX = 1024)           !max 1024 pixels per line
      PARAMETER (MTOTAL = 204800)         !but only 1024*200 in total
      PARAMETER (NA = (KMAX+1)*(KMAX+1))
      PARAMETER (NW = MTOTAL+2*MAXLIN+2*MAXLIN*(KMAX+2)+5*(1+KMAX))
C 
      INTEGER   NPIX,NLIN,ICUR(4),IMAP(4),M(MAXLIN)
      INTEGER   K,L,IT,ISTAT
      INTEGER   JOXA,JOXE,JOYA,JOYE
      INTEGER   JDXA,JDXE,JDYA,JDYE
      INTEGER   KP,KI,KX,KY,KPO,LPO,NPC,NLC,NPP,NPY
      INTEGER   MPC,MLC,JY,JX,KA,KE
      INTEGER   IN,IS,I,ITER
C 
      LOGICAL	NOISE
C
      CHARACTER CBUF*40
C
      REAL      ZZ,CRMS,OFS,CIMG(NPIX,NLIN),RAN1
C
      DOUBLE PRECISION  X(MTOTAL),Y(MAXLIN),F(MTOTAL),W(MTOTAL),A(NA)
      DOUBLE PRECISION  XMIN(MAXLIN),XMAX(MAXLIN),UX(1),UY(1),WORK(NW)
      DOUBLE PRECISION  FF(MTOTAL),YMI,YMA,XMI,XMA,YBS,SQ,DEL,DELC
      DOUBLE PRECISION  ARG,RMS,AKAP,AP
C
C     dimensions of arrays used by E02CAF and E02CBF are a function
C     of MAXLIN, MAXPIX, MTOTAL + KMAX 
C     currently MAXLIN = 1024, MAXPIX = 4096, MTOTAL = 1024*200, KMAX = 10
C     to change, use the following recipe:
C     dimension MAXLIN for array M,Y,XMIN,XMAX
C     dimension MTOTAL for array X,F,W
C     dimension (KMAX+1)**2  for array A 
C     dimension NW           for array WORK
C     NW = MTOTAL + 2*max(MAXLIN,MAXPIX) + 2*MAXLIN*(KMAX+2) + 5*(1+KMAX)
C      
      ISTAT = 0
C
C  initialize 
      JOXA = IMAP(1)
      JOXE = IMAP(2)
      JOYA = IMAP(3)
      JOYE = IMAP(4)
      JDXA = ICUR(1)
      JDXE = ICUR(2)
      JDYA = ICUR(3)
      JDYE = ICUR(4) 
      KPO = K + 1
      LPO = L + 1
C
C
C**      automatic kap*sig surface fit
C 
C
C  check size of surrounding surface...
      NPC = JOXE - JOXA + 1
      NLC = JOYE - JOYA + 1
      IF ( (NPC.GT.MAXPIX) .OR.
     +     (NLC.GT.MAXLIN) .OR.
     +     ((NPC*NLC).GT.MTOTAL) ) THEN
         ISTAT = 7
         RETURN
      ENDIF
      MPC = JDXE - JDXA + 1
      MLC = JDYE - JDYA + 1
C
C  make sure, that cursor area is within surrounding area + not undersampled...
C  This has been modified to be less restrictive in V1.1 (RNH)
      IF ( (JDXA.LT.JOXA) .OR.
     +     (JDXE.GT.JOXE) .OR.
     +     (JDYA.LT.JOYA) .OR.
     +     (JDYE.GT.JOYE) .OR.
     +     (NPC*NLC-MPC*MLC) .LT. LPO+KPO) THEN
         ISTAT = 8
         RETURN
      ENDIF
C
C  check if areas are well defined (SW)
      IF ( (JDXA.GE.JDXE) .OR.
     +     (JDYA.GE.JDYE) .OR.
     +     (JOXA.GE.JOXE) .OR.
     +     (JOYA.GE.JOYE) ) THEN
         ISTAT = 9
         RETURN
      ENDIF
C
C  define param. and subr-matrices
      YMI = JOYA
      YMA = JOYE
      XMI = JOXA
      XMA = JOXE
      IS = 0
      DO 200,JY=JOYA,JOYE
         IS = IS + 1
         Y(IS) = JY
         M(IS) = NPC
         XMIN(IS) = XMI
         XMAX(IS) = XMA
200   CONTINUE
C
C  data and weight matrices
      I = 1
      DO 600,JY = JOYA,JOYE
         DO 400,JX = JOXA,JOXE
            X(I) = JX
            F(I) = DBLE(CIMG(JX,JY))
            W(I) = 1.D0
            IF ( (JY.GE.JDYA) .AND. (JY.LE.JDYE) .AND.
     +           (JX.GE.JDXA) .AND. (JX.LE.JDXE) )     W(I) = 0.D0
            I = I + 1
400      CONTINUE
600   CONTINUE
C
C  fit a first surface to all points
C 
      CALL E02CAF(M,NLC,K,L,X,Y,F,W,MTOTAL,A,NA,XMIN,XMAX,UX,1,UY,1,
     +            WORK,NW,ISTAT)
      IF (ISTAT.NE.0) RETURN
C
C  compute initial rms 
      KA = 1
      NPP = 0
      DEL = 0.0D0
      DO 1000,JY = JOYA,JOYE
         YBS = JY
         KE = KA + NPC - 1
         CALL E02CBF(KA,KE,K,L,X,XMI,XMA,YBS,YMI,YMA,
     +               FF,A,NA,WORK,NW,ISTAT)
         DO 800,JX = KA,KE
            NPP = NPP + 1
            SQ = F(JX) - FF(JX) 
            DELC = SQ * SQ                  !DELC = (F-FF)**2
            DEL = DEL + DELC
800      CONTINUE
         KA = KE + 1
1000  CONTINUE
C
      AP = NPP - 1
      DEL = DEL/AP
      ARG = DMAX1(1.D-37,DEL)
      RMS = DSQRT(ARG)
      AKAP = 2.0D0*RMS
C
C  now iterate
C 
      IF (IT.LE.0) GO TO 8000
C 
      DO 4000,ITER=1,IT
         KA = 1
         NPP = 0
         DEL = 0.0
         DO 2000,JY=JOYA,JOYE
            YBS = JY
            KE = KA + NPC - 1
            CALL E02CBF(KA,KE,K,L,X,XMI,XMA,YBS,YMI,YMA,FF,A,NA,WORK,
     +                  NW,ISTAT)
            NPY = 0
            DO 1500,JX=KA,KE
               W(JX) = 0.0D0
               NPY = NPY + 1
               IF ( DABS(F(JX) - FF(JX)) .LE. AKAP) THEN
                  W(JX) = 1.0D0
                  NPP = NPP + 1
                  NPY = NPY - 1
                  SQ = F(JX) - FF(JX) 
                  DELC = SQ * SQ                  !DELC = (F-FF)**2
                  DEL = DEL + DELC
               ENDIF
1500        CONTINUE
C 
            KP = NPC - NPY
            IF (KP.LT.KPO)  THEN
               CBUF(1:) = 'Warning: No convergence ... '
               CALL STTPUT(CBUF,KI)
               GOTO 8000
            ELSE
               KA = KE + 1
            ENDIF
2000     CONTINUE
C
C  find new rms
         AP = NPP - 1
         DEL = DEL/AP
         ARG = DMAX1(1.D-37,DEL)
         RMS = DSQRT(ARG)
         AKAP = 2.0D0*RMS
C
C  fit new surface
         CALL E02CAF(M,NLC,K,L,X,Y,F,W,MTOTAL,A,NA,XMIN,XMAX,UX,1,UY,1,
     +               WORK,NW,ISTAT)
         IF (ISTAT.NE.0) RETURN
4000  CONTINUE
C
C  flag "bad" pixels...
8000  KI = 1
      DO 8800,KY=JDYA,JDYE
         DO 8600,KX=JDXA,JDXE
            W(KI) = KX               !W now contains indices of bad pixels
            KI = KI + 1
8600     CONTINUE
8800  CONTINUE
C      
C  and calculate new pixels
      KA = 1
      DO 9000,KY=JDYA,JDYE
         KE = KA + MPC - 1
         YBS = KY 
         CALL E02CBF(KA,KE,K,L,W,XMI,XMA,YBS,YMI,YMA,FF,A,NA,WORK,
     +               NW,ISTAT)
         IF (ISTAT.NE.0) RETURN
         KA = KE + 1
9000  CONTINUE
C
C  add noise if optional flag is set (RNH)
C 
      IF (NOISE) THEN
         OFS = RAN1(-JDXA)                       !init random generator
         IS = 1
      ENDIF
      IF (CRMS.GT.0.) RMS = DBLE(CRMS)
C 
      DO 9500,KY=JDYA,JDYE
         IN = (KY-JDYA) * MPC - JDXA + 1
         IF (NOISE) THEN
            DO 9100,KX=JDXA,JDXE
               ZZ = SNGL(FF(IN+KX))
               OFS = RAN1(0)                          !IN is in ]0,1]
               IF (OFS.LT.0.5) THEN
                  IF (IS.GT.0) THEN
                     CIMG(KX,KY) = (1.0 + 0.1*OFS) * ZZ
                     IS = 0
                  ELSE
                     CIMG(KX,KY) = (1.0 - 0.1*OFS) * ZZ
                     IS = 1
                  ENDIF
               ELSE
                  CIMG(KX,KY) = ZZ
               ENDIF
9100        CONTINUE
            L = 1000.*OFS                              !new root for RAN1
            OFS = RAN1(-L)
         ELSE
            DO 9200,KX=JDXA,JDXE
               CIMG(KX,KY) = SNGL(FF(IN+KX))
9200        CONTINUE
         ENDIF
9500  CONTINUE
C
      RETURN
      END

      REAL FUNCTION RAN1(IDUM)
C 
C    returns an uniform deviate between 0.0 and 1.0. 
C    input par. IDUM, set IDUM to a negative value to initialize the sequence
C    K. Banse    941024			(from Numerical Recipes)
C
      IMPLICIT NONE
C 
      INTEGER IDUM
      INTEGER IFF,IX1,IX2,IX3,IC1,IC2,IC3
      INTEGER J,IA1,IA2,IA3,M1,M2,M3
C 
      REAL R(97),RM1,RM2
C 
      SAVE R,IX1,IX2,IX3
C 
      DATA M1/259200/,IA1/7141/,IC1/54773/
      DATA M2/134456/,IA2/8121/,IC2/28411/
      DATA M3/243000/,IA3/4561/,IC3/51349/
      DATA IFF/0/
C 
      RM1=1./M1
      RM2=1./M2
      IF (IDUM.LT.0.OR.IFF.EQ.0) THEN
         IFF = 1
         IX1 = MOD(IC1-IDUM,M1)
         IX1 = MOD(IA1*IX1+IC1,M1)
         IX2 = MOD(IX1,M2)
         IX1 = MOD(IA1*IX1+IC1,M1)
         IX3 = MOD(IX1,M3)
         DO 10, J=1,97
            IX1 = MOD(IA1*IX1+IC1,M1)
            IX2 = MOD(IA2*IX2+IC2,M2)
            R(J) = (FLOAT(IX1)+FLOAT(IX2)*RM2)*RM1
10       CONTINUE
      ENDIF
C 
      IX1 = MOD(IA1*IX1+IC1,M1)
      IX2 = MOD(IA2*IX2+IC2,M2)
      IX3 = MOD(IA3*IX3+IC3,M3)
      J = 1+(97*IX3)/M3
      RAN1 = R(J)
      R(J) = (FLOAT(IX1)+FLOAT(IX2)*RM2)*RM1
      RETURN
C 
      END       
