C @(#)functions.for	16.1 (ESO-DMD) 06/19/01 15:22:45
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 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
      SUBROUTINE IFILL(A,NAXIS,NPIX,START,STEP,RC,NO,RMIN,RMAX)
C
C++++++++++++++++++++++++++++++++++++++++++++++++++
C
C.COPYRIGHT: Copyright (c) 1988 European Southern Observatory,
C                                         all rights reserved
C
C.LANGUAGE: F77+ESOext
C
C.AUTHOR: K.Banse
C
C.IDENTIFICATION
C  subroutine IFILL		version 2.30	881028
C             IMFUNC			2.00	881028
C             GAUSS			2.05	881028
C  K. Banse                        ESO - Garching
C  source code on FUNCTIONS.FOR
C
C.KEYWORDS
C  images, masks
C
C.PURPOSE
C  create a MIDAS image or mask frame
C
C.ALGORITHM
C  use MIDAS interfaces
C
C.INPUT/OUTPUT
C  call as IFILL(A,NAXIS,NPIX,START,STEP,RC,NO,RMIN,RMAX)
C          IMFUNC(METHOD,A,NAXIS,NPIX,START,STEP,RC,RMIN,RMAX)
C          GAUSS(A,NAXIS,NPIX,START,STEP,RC,NO,RMIN,RMAX)
C          RANDOM(METHOD,IROOT,A,NAXIS,NPIX,RC,RMIN,RMAX)
C      
C.VERSIONS
C  2.50 (IFILL)    split for 1-dim and 2-dim option
C  2.10 (IMFUNCT)  add method=3 for circle
C  3.00		work on precision for pixels (IMFILL)
C 
C  010802		last modif
C      
C--------------------------------------------------
C      
      IMPLICIT NONE
C      
      REAL         A(*),RC(*),RMIN,RMAX
C 
      DOUBLE PRECISION START(*),STEP(*)
      DOUBLE PRECISION XSTART,YSTART,XSTEP,YSTEP
      DOUBLE PRECISION WX,WY,CY,TY,FYY,PIX
C      
      INTEGER      NAXIS,NPIX(*),NO
      INTEGER      IOFF,N,NDIM,NX,NY
      INTEGER      XDIM,YDIM
C      
C -----
C      
C  branch according to dimension of image
C      
C -----
C      
C  for NAXIS > 2 only constant option supported ...
C      
      IF (NAXIS.GT.2) THEN
         NDIM = 1
         DO 100 N=1,NAXIS
            NDIM = NPIX(N) * NDIM
100      CONTINUE
C      
         DO 150 N=1,NDIM
            A(N) = RC(1)
150      CONTINUE
         RMIN = RC(1)
         RMAX = RMIN
         RETURN
      ENDIF
C      
C  for NAXIS = 2 do it elsewhere...
      IF (NAXIS.EQ.2) GOTO 10000
C      
C  here we handle NAXIS = 1: that is the polynomial a + bx + cxx
C      
      XDIM = NPIX(1)
      XSTART = START(1)
      XSTEP = STEP(1)
C      
C  branch according to no. of coefficients (= order of polynomial + 1)
      IF (NO.GT.3) NO = 3
      GOTO (1000,2000,3000),NO
C      
C  just constant
1000  DO 1100 N=1,XDIM
         A(N) = RC(1)
1100  CONTINUE
      RMIN = RC(1)
      RMAX = RMIN
      RETURN
C      
C  p(x) = a + bx
2000  WX = XSTART
      DO 2100 NX=1,XDIM
         PIX = RC(1) + RC(2)*WX
         IF (RMIN.GT.PIX) RMIN = PIX
         IF (RMAX.LT.PIX) RMAX = PIX
         A(NX) = PIX
         WX = WX + XSTEP
2100  CONTINUE
      RETURN
C      
C  p(x) = a + bx + cxx
3000  WX = XSTART 
      DO 3100 NX=1,XDIM
         PIX = RC(1) + RC(2)*WX + RC(3)*WX*WX
         IF (RMIN.GT.PIX) RMIN = PIX
         IF (RMAX.LT.PIX) RMAX = PIX
         A(NX) = PIX
         WX = WX + XSTEP
3100  CONTINUE
      RETURN
C      
C  here we handle NAXIS = 2, i.e. the polynomial a + bx + cy + dxy + exx + fyy
C      
10000 XDIM = NPIX(1)
      YDIM = NPIX(2)
      XSTART = START(1)
      YSTART = START(2)
      XSTEP = STEP(1)
      YSTEP = STEP(2)
      IOFF = 0
C      
C  branch according to no. of coefficients (= order of polynomial + 1)
      GOTO (11000,12000,13000,14000,15000,16000),NO
C      
C  just constant
11000 NDIM = 1
      DO 11010 N=1,NAXIS
         NDIM = NPIX(N) * NDIM
11010 CONTINUE
C      
      DO 11020 N=1,NDIM
         A(N) = RC(1)
11020 CONTINUE
      RMIN = RC(1)
      RMAX = RMIN
      RETURN
C      
C  p(x,y) = a + bx
12000 DO 12500 NY=1,YDIM
         WX = XSTART
         DO 12100 NX=1,XDIM
            PIX = RC(1) + RC(2)*WX
            IF (RMIN.GT.PIX) RMIN = PIX
            IF (RMAX.LT.PIX) RMAX = PIX
            A(NX+IOFF) = PIX
            WX = WX + XSTEP
12100    CONTINUE
         IOFF = IOFF + XDIM
12500 CONTINUE
      RETURN
C      
C  p(x,y) = a + bx + cy
13000 WY = YSTART
      DO 13500 NY=1,YDIM
         CY = WY * RC(3)
         WX = XSTART 
         DO 13100 NX=1,XDIM
            PIX = RC(1) + RC(2)*WX + CY
            IF (RMIN.GT.PIX) RMIN = PIX
            IF (RMAX.LT.PIX) RMAX = PIX
            A(NX+IOFF) = PIX
            WX = WX + XSTEP
13100    CONTINUE
         WY = WY + YSTEP
         IOFF = IOFF + XDIM
13500 CONTINUE
      RETURN
C      
C  p(x,y) = a + bx + cy + dxy
14000 WY = YSTART
      DO 14500 NY=1,YDIM
         CY = WY * RC(3)
         TY = WY * RC(4)
         WX = XSTART 
         DO 14100 NX=1,XDIM
            PIX = RC(1) + RC(2)*WX + CY + TY*WX
            IF (RMIN.GT.PIX) RMIN = PIX
            IF (RMAX.LT.PIX) RMAX = PIX
            A(NX+IOFF) = PIX
            WX = WX + XSTEP
14100    CONTINUE
         WY = WY + YSTEP
         IOFF = IOFF + XDIM
14500 CONTINUE
      RETURN
C      
C  p(x,y) = a + bx + cy + dxy + exx
15000 WY = YSTART
      DO 15500 NY=1,YDIM
         CY = WY * RC(3)
         TY = WY * RC(4)
         WX = XSTART 
         DO 15100 NX=1,XDIM
            PIX = RC(1) + RC(2)*WX + CY + WX*(TY + RC(5)*WX)
            IF (RMIN.GT.PIX) RMIN = PIX
            IF (RMAX.LT.PIX) RMAX = PIX
            A(NX+IOFF) = PIX
            WX = WX + XSTEP
15100    CONTINUE
         IOFF = IOFF + XDIM
         WY = WY + YSTEP
15500 CONTINUE
      RETURN
C      
C  p(x,y) = a + bx + cy + dxy + exx + fyy
16000 WY = YSTART
      DO 16500 NY=1,YDIM
         CY = WY * RC(3)
         TY = WY * RC(4)
         FYY = WY * WY * RC(6)
         WX = XSTART
         DO 16100 NX=1,XDIM
            PIX = RC(1) + RC(2)*WX + CY + WX*(TY + RC(5)*WX) + FYY
            IF (RMIN.GT.PIX) RMIN = PIX
            IF (RMAX.LT.PIX) RMAX = PIX
            A(NX+IOFF) = PIX
            WX = WX + XSTEP
16100    CONTINUE
         IOFF = IOFF + XDIM
         WY = WY + YSTEP
16500 CONTINUE
      RETURN
C      
      END

      SUBROUTINE IMFUNC(METHOD,A,NAXIS,NPIX,START,STEP,RC,RMIN,RMAX)
C      
      IMPLICIT NONE
C 
      REAL      A(*),RC(*),RMIN,RMAX
C 
      DOUBLE PRECISION XSTEP,XSTART,YSTEP,YSTART
      DOUBLE PRECISION CX,CY,CPA,SPA,CDI,SC,AMP
      DOUBLE PRECISION X1,Y1,R,PIX,DX,DY,FACT
      DOUBLE PRECISION START(*),STEP(*)
C 
      INTEGER   METHOD,NAXIS,NPIX(*)
      INTEGER   XDIM,YDIM,IPAIR,NH
      INTEGER   I,NX,NY,IRC
C 
      DATA  FACT /0.0174532925/                   !FACT = PI / 180.
C      
C  initialize + branch according to method
      XDIM = NPIX(1)
      XSTART = START(1)
      XSTEP = STEP(1)
C 
      IF (NAXIS.GE.2) THEN                  !but we only fill xy-plane...
         YDIM = NPIX(2)
         YSTART = START(2)
         YSTEP = STEP(2)
      ELSE
         YDIM = 1
         YSTART = 0.
         YSTEP = 1.
      ENDIF
C      
      IPAIR = 1 - MOD(XDIM,2)                  !get center coords.
      NH = (XDIM-1)/2      
      CX = XSTART + NH*XSTEP + IPAIR*XSTEP*0.5
      IPAIR = 1 - MOD(YDIM,2)      
      NH = (YDIM-1)/2      
      CY = YSTART + NH*YSTEP + IPAIR*YSTEP*0.5
C      
      GOTO (1000,2000,3000,4000,5000,6000,7000,8000),METHOD
C      
C +++++
C  build exponential disk
C +++++
C      
1000  CPA = COS(RC(3)*FACT)
      SPA = SIN(RC(3)*FACT)
      CDI = COS(RC(4)*FACT)
      SC = -1./RC(2)
      AMP = RC(1)
C      
C  now loop
      I = 1
      DY = YSTART - CY
      DO 1500 NY=1,YDIM
         DX = XSTART - CX
         DO 1100 NX=1,XDIM
            X1 = SPA*DX - CPA*DY
            Y1 = CDI * (CPA*DX + SPA*DY)
            R = DSQRT(X1*X1 + Y1*Y1)
            PIX = AMP * EXP(R*SC)
            IF (RMIN.GT.PIX) RMIN = PIX
            IF (RMAX.LT.PIX) RMAX = PIX
            A(I) = PIX
            I = I + 1
            DX = DX + XSTEP
1100     CONTINUE
         DY = DY + YSTEP
1500  CONTINUE
      RETURN
C      
C +++++
C  build frame according to radius law
C +++++
C      
2000  CPA = COS(RC(3)*FACT)
      SPA = SIN(RC(3)*FACT)
      CDI = COS(RC(4)*FACT)
      SC = 1./RC(2)
      AMP = RC(1)
C      
C  now loop
      I = 1
      DY = YSTART - CY
      DO 2500 NY=1,YDIM
         DX = XSTART - CX
         DO 2100 NX=1,XDIM
            X1 = SPA*DX - CPA*DY
            Y1 = CDI * (CPA*DX + SPA*DY)
            R = DSQRT(X1*X1 + Y1*Y1)
            PIX = AMP * 10.**(-3.33*(R*SC)**0.25)
            IF (RMIN.GT.PIX) RMIN = PIX
            IF (RMAX.LT.PIX) RMAX = PIX
            A(I) = PIX
            I = I + 1
            DX = DX + XSTEP
2100     CONTINUE
         DY = DY + YSTEP
2500  CONTINUE
      RETURN
C      
C +++++
C  build frame with circle in the center
C +++++
C      
3000  AMP = RC(1) * RC(1)                  !radius squared...
C      
C  now loop
      I = 1
      DY = YSTART
      DO 3500 NY=1,YDIM
         DX = XSTART 
         DO 3100 NX=1,XDIM
            SC =  ((DX-CX)*(DX-CX)) + ((DY-CY)*(DY-CY))
            IF (SC.LE.AMP) THEN
             PIX = RC(2)
            ELSE
             PIX = RC(3)
            ENDIF
            IF (RMIN.GT.PIX) RMIN = PIX
            IF (RMAX.LT.PIX) RMAX = PIX
            A(I) = PIX
            I = I + 1
            DX = DX + XSTEP
3100     CONTINUE
         DY = DY + YSTEP
3500  CONTINUE
      RETURN
C      
C +++++
C  build frame with ellips in the center
C +++++
C      
4000  CPA = 1./(RC(1)*RC(1))                  !get 1/a**2, 1/b**2
      SPA = 1./(RC(2)*RC(2))
C      
C  now loop
      I = 1
      DY = YSTART
      DO 4500 NY=1,YDIM
         DX = XSTART 
         DO 4100 NX=1,XDIM
            SC =  ( (DX-CX)*(DX-CX)*CPA ) + ( (DY-CY)*(DY-CY)*SPA )
            IF (SC.LE.1.) THEN
               PIX = RC(3)
            ELSE
               PIX = RC(4)
            ENDIF
            IF (RMIN.GT.PIX) RMIN = PIX
            IF (RMAX.LT.PIX) RMAX = PIX
            A(I) = PIX
            I = I + 1
            DX = DX + XSTEP
4100     CONTINUE
         DY = DY + YSTEP
4500  CONTINUE
      RETURN
C      
C +++++
C  build Butterworth lowpass filter
C +++++
C      
5000  IRC = NINT(RC(2)) * 2                  !get 2n
      RC(1) = 1./RC(1)                       !get 1/Do
      IF ( (IRC.EQ.2) .AND.
     +     ((ABS(RC(3))-1.).LE.10.E-12) ) GOTO 5500            !simple filter
C      
C  now loop
      I = 1
      DY = YSTART
      DO 5200 NY=1,YDIM
         DX = XSTART 
         DO 5100 NX=1,XDIM
            CPA = (DX-CX)
            SPA = (DY-CY)
            SC =  DSQRT( CPA*CPA + SPA*SPA )           !get distance from center
            PIX = 1. + ( RC(3) * (SC*RC(1))**IRC )
            PIX = 1./PIX
            IF (RMIN.GT.PIX) RMIN = PIX
            IF (RMAX.LT.PIX) RMAX = PIX
            A(I) = PIX
            I = I + 1
            DX = DX + XSTEP
5100     CONTINUE
         DY = DY + YSTEP
5200  CONTINUE
      RETURN
C      
C  here we loop through the simple version
5500  I = 1
      DY = YSTART
      DO 5700 NY=1,YDIM
         DX = XSTART 
         DO 5600 NX=1,XDIM
            CPA = (DX-CX)
            SPA = (DY-CY)
            SC =  DSQRT( CPA*CPA + SPA*SPA )           !get distance from center
            SC = SC * RC(1)
            PIX = 1. / ( 1. + (SC * SC) )
            IF (RMIN.GT.PIX) RMIN = PIX
            IF (RMAX.LT.PIX) RMAX = PIX
            A(I) = PIX
            I = I + 1
            DX = DX + XSTEP
5600     CONTINUE
         DY = DY + YSTEP
5700  CONTINUE
      RETURN
C      
C +++++
C  build exponential lowpass filter
C +++++
C      
6000  IRC = NINT(RC(2))                   !get n
      RC(1) = 1./RC(1)                    !get 1/Do
      IF ( (IRC.EQ.1) .AND.
     +     ((ABS(RC(3))-1.).LE.10.E-12) ) GOTO 6500            !simple filter
C      
C  now loop
      I = 1
      DY = YSTART
      DO 6200 NY=1,YDIM
         DX = XSTART 
         DO 6100 NX=1,XDIM
            CPA = (DX-CX)
            SPA = (DY-CY)
            SC =  DSQRT( CPA*CPA + SPA*SPA )           !get distance from center
            PIX = DEXP ( -(RC(3) * (SC*RC(1))**IRC) )
            IF (RMIN.GT.PIX) RMIN = PIX
            IF (RMAX.LT.PIX) RMAX = PIX
            A(I) = PIX
            I = I + 1
            DX = DX + XSTEP
6100     CONTINUE
         DY = DY + YSTEP
6200  CONTINUE
      RETURN
C      
C  here we loop through the simple version
6500  I = 1
      DY = YSTART
      DO 6700 NY=1,YDIM
         DX = XSTART 
         DO 6600 NX=1,XDIM
            CPA = (DX-CX)
            SPA = (DY-CY)
            SC =  DSQRT( CPA*CPA + SPA*SPA )           !get distance from center
            PIX = DEXP ( -(SC*RC(1)) )
            IF (RMIN.GT.PIX) RMIN = PIX
            IF (RMAX.LT.PIX) RMAX = PIX
            A(I) = PIX
            I = I + 1
            DX = DX + XSTEP
6600     CONTINUE
         DY = DY + YSTEP
6700  CONTINUE
      RETURN
C      
C +++++
C  build Butterworth highpass filter
C +++++
C      
7000  IRC = NINT(RC(2)) * 2                  !get 2n
      AMP = MIN(ABS(XSTEP),ABS(YSTEP))       !take minimum of steps
      IF ( (IRC.EQ.2) .AND.
     +     ((ABS(RC(3))-1.).LE.10.E-12) ) GOTO 7500            !simple filter
C      
C  now loop
      I = 1
      DY = YSTART
      DO 7200 NY=1,YDIM
         DX = XSTART 
         DO 7100 NX=1,XDIM
            CPA = (DX-CX)
            SPA = (DY-CY)
            SC =  DSQRT( CPA*CPA + SPA*SPA )           !get distance from center
            IF (SC.GE.AMP) THEN
               PIX = 1. + ( RC(3) * (RC(1)/SC)**IRC )
               PIX = 1.0D0/PIX
            ELSE
               PIX = 0.0D0
            ENDIF
            IF (RMIN.GT.PIX) RMIN = PIX
            IF (RMAX.LT.PIX) RMAX = PIX
            A(I) = PIX
            I = I + 1
            DX = DX + XSTEP
7100     CONTINUE
         DY = DY + YSTEP
7200  CONTINUE
      RETURN
C      
C  here we loop through the simple version
7500  I = 1
      DY = YSTART
      DO 7700 NY=1,YDIM
         DX = XSTART 
         DO 7600 NX=1,XDIM
            CPA = (DX-CX)
            SPA = (DY-CY)
            SC =  DSQRT( CPA*CPA + SPA*SPA )           !get distance from center
            IF (SC.GE.AMP) THEN
               PIX = 1. + ( RC(1) / SC )
               PIX = 1.0D0/PIX
            ELSE
               PIX = 0.0D0
            ENDIF
            IF (RMIN.GT.PIX) RMIN = PIX
            IF (RMAX.LT.PIX) RMAX = PIX
            A(I) = PIX
            I = I + 1
            DX = DX + XSTEP
7600     CONTINUE
         DY = DY + YSTEP
7700  CONTINUE
      RETURN
C      
C +++++
C  build exponential highpass filter
C +++++
C      
8000  IRC = NINT(RC(2))                     !get n
      AMP = MIN(ABS(XSTEP),ABS(YSTEP))      !take minimum of steps
      IF ( (IRC.EQ.1) .AND.
     +     ((ABS(RC(3))-1.).LE.10.E-12) ) GOTO 8500            !simple filter
C      
C  now loop
      I = 1
      DY = YSTART
      DO 8200 NY=1,YDIM
         DX = XSTART 
         DO 8100 NX=1,XDIM
            CPA = (DX-CX)
            SPA = (DY-CY)
            SC =  DSQRT( CPA*CPA + SPA*SPA )           !get distance from center
            IF (SC.GE.AMP) THEN
               PIX = DEXP ( -(RC(3) * (RC(1)/SC)**IRC) )
            ELSE
               PIX = 0.
            ENDIF
            IF (RMIN.GT.PIX) RMIN = PIX
            IF (RMAX.LT.PIX) RMAX = PIX
            A(I) = PIX
            I = I + 1
            DX = DX + XSTEP
8100     CONTINUE
         DY = DY + YSTEP
8200  CONTINUE
      RETURN
C      
C  here we loop through the simple version
8500  I = 1
      DY = YSTART
      DO 8700 NY=1,YDIM
         DX = XSTART 
         DO 8600 NX=1,XDIM
            CPA = (DX-CX)
            SPA = (DY-CY)
            SC =  DSQRT( CPA*CPA + SPA*SPA )           !get distance from center
            IF (SC.GE.AMP) THEN
               PIX = DEXP ( -(RC(1)/SC) )
            ELSE
               PIX = 0.
            ENDIF
            IF (RMIN.GT.PIX) RMIN = PIX
            IF (RMAX.LT.PIX) RMAX = PIX
            A(I) = PIX
            I = I + 1
            DX = DX + XSTEP
8600     CONTINUE
         DY = DY + YSTEP
8700  CONTINUE
      RETURN
C      
      END

      SUBROUTINE GAUSS(A,NAXIS,NPIX,START,STEP,RC,NO,RMIN,RMAX)
C      
      IMPLICIT NONE
C 
      REAL         A(*),RC(*),RMIN,RMAX
C 
      DOUBLE PRECISION PI2,WX,WY,WXX,WYY,W1,W2,F1,PIX
      DOUBLE PRECISION X,Y
      DOUBLE PRECISION START(*),STEP(*)
C 
      INTEGER      NAXIS,NPIX(*),NO
      INTEGER      NOFF,N,NX,NY
C      
      DATA   PI2 /6.2831853071/                   ! PI2 = 2 * PI
C      
C  branch according to NAXIS
      GOTO (1000,2000),NAXIS
C      
C  1-dim gaussian function:
C  f(x) = 1/(sigma*sqrt(2pi)) * exp(-(x-mean)**2/2sigma**2
C 
1000  W1 = DSQRT(PI2)
      W2 = 1./(2.*RC(2)*RC(2))
      F1 = 1./(RC(2)*W1)
      WX = START(1) - RC(1)                !do 1st pixel twice for MIN, MAX
      WXX = WX * WX
      PIX = F1 * DEXP(-WXX*W2)
      RMIN = PIX
      RMAX = RMIN
C 
      X = START(1)
      DO 1100, N=1,NPIX(1)
         WX = X - RC(1)
         WXX = WX * WX
         PIX = F1 * DEXP(-WXX*W2)
         IF (RMIN.GT.PIX) RMIN = PIX
         IF (RMAX.LT.PIX) RMAX = PIX
         A(N) = PIX
         X = X + STEP(1)
1100  CONTINUE
      RETURN
C      
C  2-dim gaussian function: 
C  f(x,y) = 1/(xsigma*ysigma*2pi) * 
C           exp(-(x-xmean)**2/2xsigma**2 + -(y-ymean)**2/2ysigma**2)
C 
2000  F1 = 1./ (RC(2)*RC(4)*PI2)
      W1 = 1./ (2. * RC(2) * RC(2))
      W2 = 1./ (2. * RC(4) * RC(4))
C 
      WX = START(1) - RC(1)                    !do 1st pixel twice for MIN, MAX
      WXX = WX * WX
      WY = START(2) - RC(3)
      WYY = WY * WY
      PIX = F1 * DEXP(-(WXX*W1 + WYY*W2))
      RMIN = PIX
      RMAX = RMIN
C 
      NOFF = 0
      Y = START(2)
C 
      DO 2200, NY=1,NPIX(2)
         WY = Y - RC(3)
         WYY = WY * WY
         X = START(1)
C 
         DO 2100 NX=1,NPIX(1)
            WX = X - RC(1)
            WXX = WX * WX
            PIX = F1 * DEXP(-(WXX*W1 + WYY*W2))
            IF (RMIN.GT.PIX) RMIN = PIX
            IF (RMAX.LT.PIX) RMAX = PIX
            A(NX+NOFF) = PIX
            X = X + STEP(1)
2100     CONTINUE
C 
         Y = Y + STEP(2)
         NOFF = NOFF + NPIX(1)
2200  CONTINUE
C 
      RETURN
      END