C @(#)necback.for	19.1 (ES0-DMD) 02/25/03 14:20:22
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
C
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C. IDENTIFICATION
C
C   Program    ECHBACK            version 1.00 : 27.06.90
C
C   Author : P. BALLESTER       ESO - Garching
C
C. KEYWORDS
C   Echelle spectroscopy background rebinning
C
C. PURPOSE
C   Fit the background and generates a frame
C
C. ALGORITHM
C Adapted from W. Verschueren, H. Hensberge, Astr. Astroph. 240,216(1990)
C "Order extraction and background subtraction on CCD Caspec Echelle Spectra"
C
C. INPUT/OUTPUT
C   The following keywords are used :
C     IN_A/C/1/60           IN: input table (type ORDER)
C     IN_B/C/1/60           IN: reference image
C     INPUTI/I/1/1          IN: half-size of the measure window
C     INPUTI/I/2/1          IN: degree of the spline polynomials
C     INPUTR/R/1/1          IN: smoothing factor of the spline interpolation
C     OUT_A/C/1/60          OUT: output image (background estimate)
C
C. VERSIONS
C   0.00   27.06.90   from version 1.20 of 19.02.90 of REBIN.FOR
C   1.00   17.07.90   Reference version
C   1.10   19.07.90   Speed improvements : from 9.26 Mmbbe to 2.06 Mmbbe
C          22.06.99   Measurement mode (MMODE) added (swolf@eso.org)
C. BUGS
C   Since the STDWRC standard interface does not work properly for felem=-1,
C     the first line of the HISTORY descriptor of the output frame is
C     overwritten (17.07.90)
C
C-----------------------------------------------------------------------------
C
      IMPLICIT NONE
C
      LOGICAL FLAG,NULL
C     LOGICAL DBG
C
      INTEGER   NAXIS,ROW,ABSROW,NPOINT,DEGREE
      INTEGER   WIDTH,TYP,DUN,MMODE
C     INTEGER   COL   !  DBG
      INTEGER   IMNOA,IMNOB
      INTEGER*8 PNTRA,PNTRB
      INTEGER   NPIX(2),ERROR
      INTEGER   IAV,UNI(1),NULO,STAT
      INTEGER   MADRID(1)
      INTEGER   TID,NCINT,NRINT,NACINT,NARINT,NSINT
      INTEGER   POSX,POSY,POSBK
      INTEGER   PSXSTA,PSXEND,PSYSTA,PSYEND
      REAL    XSTA,XEND,YSTA,YEND
C
      CHARACTER*60  ORDER,REFIMA,BAKIMA
      CHARACTER*80  STRING,HISTO
      CHARACTER     CUNIT*48,IDENT*72
      CHARACTER*16  LABX,LABY,LABBK,LBXSTA,LBXEND
      CHARACTER*16  LBYSTA,LBYEND,FORM,UNIT
C
      PARAMETER (NPOINT=20000)
      REAL      BACK(3,NPOINT),SMOOTH,MEASUR,CUTS(4)
C
      DOUBLE PRECISION STEP(2),START(2)
C
      INCLUDE 'MID_INCLUDE:ST_DEF.INC'
C
      COMMON /VMR/ MADRID
C
      INCLUDE 'MID_INCLUDE:ST_DAT.INC'
C
      DATA LABX,LABY,LABBK /'X' , 'YBKG' , 'BKG'/
      DATA LBXSTA,LBXEND,LBYSTA,LBYEND /'XSTA','XEND','YSTA','YEND'/
C
C --- 0. MIDAS environment set up and enables automatic error abort
C
      CALL STSPRO ('ECHBAK')
C
C --- Set DEBUG flag or not
C
C      DBG = (.FALSE.)
C
C --- 1. Read input/output descriptors
C
      CALL STKRDC ('IN_A',1,1,60,IAV,ORDER,UNI,NULO,STAT)
      CALL STKRDC ('IN_B',1,1,60,IAV,REFIMA,UNI,NULO,STAT)
      CALL STKRDC ('OUT_A',1,1,60,IAV,BAKIMA,UNI,NULO,STAT)
C
C --- Translate names
C
      CALL CLNTAB (ORDER,ORDER,0)
      CALL CLNFRA (REFIMA,REFIMA,0)
      CALL CLNFRA (BAKIMA,BAKIMA,0)
C
      CALL STKRDI ('INPUTI',1,1,IAV,WIDTH,UNI,NULO,STAT)
      CALL STKRDR ('INPUTR',1,1,IAV,SMOOTH,UNI,NULO,STAT)
      CALL STKRDI ('INPUTI',2,1,IAV,DEGREE,UNI,NULO,STAT)
      CALL STKRDI ('INPUTI',3,1,IAV,MMODE,UNI,NULO,STAT)
C
C --- Display read descriptors
C
C       IF (DBG) THEN
C        STRING = '*** Read descriptors : '
C        CALL STTPUT (STRING,STAT)
C        STRING = ORDER(1:25)//REFIMA(1:25)//BAKIMA(1:25)
C        CALL STTPUT (STRING,STAT)
C       ENDIF
C
C --- 2. Get reference frame and map it.
C
      CALL STIGET (REFIMA,D_R4_FORMAT,F_I_MODE,F_IMA_TYPE,
     +             2,NAXIS,NPIX,START,STEP,IDENT,CUNIT,
     +             PNTRA,IMNOA,STAT)
C
C --- A call to STETER will abort the program...
C
      IF (NAXIS.NE.2)
     +     CALL STETER (1,'Input frame must have two dimensions...')
C
C --- 3.1. Get input table and open it
C
      CALL TBTOPN (ORDER,F_IO_MODE,TID,STAT)
C
C --- 3.2. Read table descriptors
C
      CALL TBIGET (TID,NCINT,NRINT,NSINT,NACINT,NARINT,STAT)
      IF (NRINT.EQ.0) THEN
          STRING = '*** FATAL : There are no data in the table'
          CALL STETER (9,STRING)
      ENDIF
C
C --- 4. Table columns processing
C
C --- 4.1. Search for columns defining background positions
C
      CALL TBLSER (TID,LABX ,POSX, STAT)
      IF (POSX.EQ.-1) THEN
         STRING = 'Column label '//LABX//'not found in table '//ORDER
         CALL STTPUT (STRING,STAT)
      ENDIF
C
      CALL TBLSER (TID,LABY ,POSY, STAT)
      IF (POSY.EQ.-1) THEN
         STRING = 'Column label '//LABY//'not found in table '//ORDER
         CALL STTPUT (STRING,STAT)
      ENDIF

      CALL TBLSER (TID,LBXSTA ,PSXSTA, STAT)
      CALL TBLSER (TID,LBXEND ,PSXEND, STAT)
      CALL TBLSER (TID,LBYSTA ,PSYSTA, STAT)
      CALL TBLSER (TID,LBYEND ,PSYEND, STAT)
C
C --- 4.2. Create column for measured background
C          Check that the column do not already exist and if not, creates
C          a new one.
C
      CALL TBLSER (TID,LABBK,POSBK,STAT)

      IF (POSBK.GT.0) THEN
         STRING = 'Warning : Column label '//LABBK//'already exists'
         CALL STTPUT (STRING,STAT)
         STRING = '          Not created'
         CALL STTPUT (STRING,STAT)
      ELSE
         TYP  = D_R4_FORMAT
         FORM = 'E12.3'
         UNIT = 'PIXELS'
         CALL TBCINI (TID,TYP,1,FORM,UNIT,LABBK,POSBK,STAT)
      ENDIF
C
C --- Displays columns positions
C
C      IF (DBG) THEN
C       STRING = LABX(1:16)//LABY(1:16)//LABBK(1:16)
C       CALL STTPUT (STRING,STAT)
C       WRITE (STRING,200) POSX,POSY,POSBK
C 200   FORMAT (3(I3,13X))
C       CALL STTPUT (STRING,STAT)
C      ENDIF
C
C --- 5. Prepares array BACK by reading selected position values in the table
C         and measures local background level on these positions.
C
      ABSROW=0
      DO 100 ROW = 1,NRINT
C
C --- 5.1. Read selection flag
C
         CALL TBSGET (TID,ROW,FLAG,STAT)

         IF (STAT.NE.0) THEN
          WRITE (STRING,350) ROW,ORDER
350       FORMAT ('*** Fatal : Cannot access row ',I5,' of table ',A40)
          CALL STETER (9,STRING)
         ENDIF
C
         IF (FLAG) THEN

             ABSROW = ABSROW + 1

             IF (ABSROW.GT.NPOINT) THEN
          STRING = 'Do you really want to process such a big table ???'
          CALL STTPUT (STRING,STAT)
          STRING= '*** Warning : Buffer overflow in module FITBACK.FOR'
          CALL STETER (9,STRING)
             ENDIF
C
C --- 5.2. Read positions
C
             CALL TBERDR (TID,ROW,POSX, BACK(1,ABSROW),NULL,STAT)
             CALL TBERDR (TID,ROW,POSY, BACK(2,ABSROW),NULL,STAT)
             CALL TBERDR (TID,ROW,PSXSTA, XSTA,NULL,STAT)
             CALL TBERDR (TID,ROW,PSXEND, XEND,NULL,STAT)
             CALL TBERDR (TID,ROW,PSYSTA, YSTA,NULL,STAT)
             CALL TBERDR (TID,ROW,PSYEND, YEND,NULL,STAT)
C
C --- 5.3. Measures background
C
             BACK (3,ABSROW) = MEASUR (MADRID(PNTRA),NPIX(1),NPIX(2),
     +             XSTA,XEND,YSTA,YEND,WIDTH,MMODE)
C
C --- 5.4. Write measured background on the table
C
             CALL TBEWRR (TID,ROW,POSBK,BACK(3,ABSROW),STAT)
C
         ENDIF

C         IF (DBG) THEN
C          WRITE (STRING,300) ROW,FLAG,ABSROW,(BACK(COL,ABSROW),COL=1,3)
C300       FORMAT(3I5,3F12.3)
C          CALL STTPUT (STRING,STAT)
C         ENDIF

100   CONTINUE
C
C --- 5.5 If DEGREE<0 EXIT
C      
      IF (DEGREE.LT.0) GOTO 1000

C
C --- 6. Get output frame
C
      IF (REFIMA.EQ.BAKIMA) THEN
         CALL STETER (2,'Input and Output frames must be different...')
      ENDIF
C
      CALL STIPUT (BAKIMA,D_R4_FORMAT,F_O_MODE,F_IMA_TYPE,
     +             NAXIS,NPIX,START,STEP,IDENT,CUNIT,
     +             PNTRB,IMNOB,STAT)
C
C --- 7. Update descriptors of the output frame
C
C --- 7.1 Copy all descriptors from the input frame to the output frame
C
      CALL STDCOP (IMNOA,IMNOB,1,' ',STAT)
C
C --- 7.2. Write specific IDENT on the output frame
C
      WRITE (HISTO,800) WIDTH,SMOOTH,DEGREE
800   FORMAT ('Background estimate -  Half-width:',I3,' Smooth:',E13.4,
     +        ' Degree:',I3)
      CALL STTPUT (HISTO,STAT)
C
C
      CALL STDWRC (IMNOB,'ORIGIN',1,REFIMA,1,60,DUN,STAT)
      CALL STDWRC (IMNOB,'HISTORY',1,HISTO,1,80,DUN,STAT)
C
C      CALL STDWRC (IMNOB,'HISTORY',1,HISTO,-1,80,DUN,STAT)
C. ---> STDWRC does not work properly for felem = -1. (17.07.90)
C
C
C --- 8. Call the routine which actually process the data
C
C     - Check initial conditions
C     - Do the smoothing spline interpolation along the columns.
C     - Interpolates background row by row for residual undefined positions
C     - Writes interpolated values on the frame
C     - Returns an error code (0 = O.K.)
C  !!! The subroutine modifies the value of SMOOTH !!!
C
      CALL FITBAK (NPOINT,BACK,ABSROW,MADRID(PNTRB),NPIX(1),NPIX(2),
     +             DEGREE,SMOOTH,ERROR,CUTS)
C
      IF (ERROR.NE.0) THEN
         STRING = '*** An error occured during execution'
         CALL STETER (9,STRING)
      ENDIF
C
C --- 8.2 Finish descriptor updating
C
      CALL STDWRR (IMNOB,'LHCUTS',CUTS,1,4,DUN,STAT)
C
C --- 9. Release files, update keywords and exit
C
 1000 CONTINUE
      CALL STSEPI
C
      STOP
      END




      SUBROUTINE FITBAK (DIMTAB,TABLE,NBVAL,IMAGE,NX,NY,DEGREE,SMOOTH,
     +                   ERR,CUTS)
C
C     
C
C     DESCRIPTION :
C
C     This subroutine do the actual work of inter- or extrapolation
C     using a smoothing spline algorithm to fit the data in mainly
C     two steps :
C
C     1. First, assuming that background positions are defined column
C     by column, one fits the data for all defined columns. This operation
C     involves the following steps.
C
C       - Search for a defined column in the table and load X,Y and W arrays
C            which represents respectively the row position, the background
C            measured value and the relative weight of the value (here W=1.0).
C       - Updates the array of defined columns (TABFLG)
C       - If necessary, forces the values on the edges of the column
C       - Computes knots and b-spline coefficients
C       - Computes interpolated values along the whole column and writes
C            them on the frame.
C
C     2. Residual undefined pixels are interpolated using a row by row
C     spline interpolation. This time, one forces the spline interpolation
C     to pass exactly by the already defined values (smoothing factor s=0.).
C     This is done in a similar way of the previous step :
C
C       - In order to accelerate computation :
C          One loads an array of residual undefined columns positions (UNDCOL)
C          Since this array is constant, so are the knots positions of the
C          spline interpolation, which are computed only once, for the
C          first row.
C       - Computes the ratio of the job already done and displays it.
C       - Computes b-spline coefficients for each row
C       - Computes interpolated values along the whole row and writes
C          them on the frame.
C
C     If something is wrong during the execution, an error message is
C     displayed, trying to fix the problem. An error code is returned to
C     the driver program and the execution immediately stopped. 
C
C
C     INPUT/OUTPUT :
C
C     ERR   is the only output variable
C     IMAGE is an input/output variable
C     All others variables are input data
C
C     DIMTAB     INTEGER      Second dimension of the array TABLE
C                             (First dimension is supposed to be 3)
C     TABLE      REAL array   Col. 1 : X posiitons of background
C                             Col. 2 : Y posiitons of background
C                             Col. 3 : measured values of background
C     NBVAL      INTEGER      Number of significant data in array TABLE
C     IMAGE      REAL array   Adress of the output image
C     NX         INTEGER      X-dimension of the output image
C     NY         INTEGER      Y-dimension of the output image
C     DEGREE     INTEGER      Degree of the spline polynomials
C     SMOOTH     REAL         Smoothing factor of the spline interpolation
C     ERR        INTEGER      Error return code (0 = O.K.)
C
C     VERSIONS :
C
C     V. 0.00   27.06.90   Creation
C     V. 1.00   17.06.90   Reference version
C     V. 1.10   19.07.90   Remove MASK, introduce CUTS computation
C          Speed improv. : (use PREPOS, remove IMAGE init. and W init)
C
C --- 1. Variables declarations
C
      IMPLICIT NONE
C
C --- 1.1. Input/Output variables
C
      INTEGER   DIMTAB
      REAL      TABLE(3,DIMTAB),CUTS(4)
      INTEGER   NBVAL,NX,NY,DEGREE,ERR
      REAL      IMAGE(1),SMOOTH
C
C --- 1.2. Internal variables
C
      INTEGER   M,NEST,LEVEL
      PARAMETER (M=10000,NEST=4096)
      REAL      X(M),Y(M),W(M),Q(M,6),XB,XE,KNOTS(NEST),COEF(NEST),FP
      REAL      NRDATA(NEST),FPO,FPOLD,NPLUS,THRES
      INTEGER   IOPT,NBKNOT,IER,NBUND,PREPOS
      REAL      ARG,DERIV,UNDCOL(M),MEAN,RATIO
      INTEGER   PTIMA,PTTAB,STAT,FLAG,NK1,NU,PASS1
      INTEGER   ROW,COL,POSARG,POSCOL,POSTAB,NPTS,TABFLG(M,2)
      INTEGER   INTCOL,INTROW
      REAL      BIGSMO
      CHARACTER STRING*80
C
C --- 2. Initialisations
C
       IF (NX.GT.M.OR.NY.GT.M) THEN
        STRING = 'Do you really want to process such a big image ???'
        CALL STTPUT (STRING,STAT)
        STRING = '*** Fatal : Parameter M in module FITBAK is too small'
        CALL STTPUT (STRING,STAT)
        IF (NY.GT.NX) NX = NY
        WRITE (STRING,5100) NX,M
5100    FORMAT ('Minimum value is : ',I8,' instead of :',I8)
        CALL STTPUT (STRING,STAT)
        ERR = 10		! No "science" error
        RETURN
       ENDIF
C
C
      DO 50 COL = 1 , M
         TABFLG (COL,1) = 0	! Number of data for interpolation
         TABFLG (COL,2) = 0	! Interpolation done (=1)
         W (COL) = 1.0		! Weights of values
50    CONTINUE

      ERR    = 0	! Normal return
      LEVEL  = 2	! Level of messages when calling IERSPL
C      SMOOTH = 0.	! Smoothing factor for this part of the algorithm
      INTCOL  = 0
      INTROW  = 0
      BIGSMO  = 1.0E+20
C
C --- 3. Building the table of values
C
100   CONTINUE		! Beginning of the loop on table values
C
C --- 3.1. Identification of a new column position not yet interpolated
C
      FLAG = 0
      DO 110 PTTAB = 1 , NBVAL
         COL = NINT (TABLE(1,PTTAB))
         IF (TABFLG(COL,1).EQ.0) THEN
            FLAG   = 1
            POSTAB = PTTAB
            GOTO 150
         ENDIF
110   CONTINUE
150   CONTINUE
C
C --- 3.2. If FLAG#0, loading of table data for the given column number
C
      IF (FLAG.EQ.0) THEN
cd         STRING = '*** No new columns'
cd         CALL STTPUT (STRING,STAT)
         GOTO 1000
      ELSE

         NPTS = 0	! Number of ORDER table values for column COL
         DO 200 PTTAB = POSTAB , NBVAL
C
            POSCOL = NINT (TABLE(1,PTTAB))
            IF (COL.EQ.POSCOL) THEN
C
                NPTS = NPTS + 1

       IF (NPTS.GT.M) THEN
        STRING = '*** Fatal : Parameter M in module FITBAK is too small'
        CALL STTPUT (STRING,STAT)
        WRITE (STRING,5100) NPTS,M
        ERR = 10		! No "science" error
        RETURN
       ENDIF

                X(NPTS) = TABLE (2,PTTAB)	! Y position, in pixels
                Y(NPTS) = TABLE (3,PTTAB)	! Background value
cd                W(NPTS) = 1.0			! Weight
            ENDIF
200      CONTINUE
       ENDIF
C
C --- 3.3. Update the array of interpolated columns and check number of data
C
         TABFLG(COL,1) = NPTS
cd         WRITE (STRING,5000) COL,TABFLG(COL,1)
cd         CALL STTPUT (STRING,STAT)
cd5000     FORMAT ('Process column : ',I5,' (',I4,' pixels)')
C
        IF (TABFLG(COL,1).LE.DEGREE) THEN
         STRING = '*** Warning : Number insufficient for interpolation'
         CALL STTPUT (STRING,STAT)
         GOTO 100	! Search for another column number
        ENDIF
C
C --- 4. Smoothing spline interpolation along columns
C
             IOPT = 0		! Restart all computations
C
C --- 4.1. Evaluating minimum abscisse XB and maximum XE and eventually,
C          forces the values on the edges of the column
C
C     CALL SORSPL (X,Y,W,NPTS) ! If it is necessary to sort the data
C
             IF (NINT(X(1)).NE.1) THEN
                DO 350 PTTAB = NPTS , 1 , -1
                  X(PTTAB+1) = X(PTTAB)
                  Y(PTTAB+1) = Y(PTTAB)
cd                  W(PTTAB+1) = W(PTTAB)	! Speed optim.
350             CONTINUE
                X(1) = 1.0
                Y(1) = Y(2)
cd                W(1) = W(2)			! Speed optim.
                NPTS = NPTS + 1
             ENDIF

             IF (NINT(X(NPTS)).NE.NY) THEN
                X(NPTS+1) = FLOAT(NY)
                Y(NPTS+1) = Y(NPTS)
cd                W(NPTS+1) = W(NPTS)		! Speed optim.
                NPTS = NPTS + 1
             ENDIF

             XB = 1.0
             XE = FLOAT(NY)

C
C --- 4.2. Computation of knots sequence and b-spline coefficients
C
             CALL SMOOT (X,Y,W,Q,NPTS,XB,XE,DEGREE,SMOOTH,NBKNOT,
     +                    KNOTS,COEF,FP,IOPT,IER,NRDATA,FPO,FPOLD,NPLUS)
C
C --- 4.3. Error management
C
             IF (IER.EQ.3) THEN  !   No convergence achieved

             CALL SMOOT (X,Y,W,Q,NPTS,XB,XE,DEGREE,BIGSMO,NBKNOT,
     +                    KNOTS,COEF,FP,0,IER,NRDATA,FPO,FPOLD,NPLUS)
             PASS1  = 0
             INTCOL = INTCOL + 1

             ENDIF

             CALL IERSPL (IER,LEVEL)	! Displays error messages

             IF (IER.GT.0) THEN

              ERR = 1		! Error during spline interpolation
5200          FORMAT ('Error while processing column ',I6)
              WRITE(STRING,5200) POSCOL
              CALL STTPUT(STRING,STAT)
              RETURN

             ELSE

              TABFLG (COL,2) = 1 ! Interpolation correctly done

             ENDIF
C
C --- 5. Computation of interpolated values
C
      NU  = 0				! Order of the derivative
      PREPOS = DEGREE + 1

      DO 400 ARG = XB , XE , 1.0
C
C --- 5.1. Find the position of the argument within the knots sequence
C
         CALL ARGSPL (ARG,KNOTS,NBKNOT,DEGREE,PREPOS,POSARG,NK1)
C
C --- 5.2. Computes interpolated values
C
C        COL is already defined within the main loop
         ROW = NINT(ARG)
         PTIMA = (ROW-1)*NX+COL
         IMAGE(PTIMA)  = DERIV (KNOTS,NBKNOT,COEF,NK1,NU,ARG,POSARG)

cd         WRITE (STRING,390) POSARG,ARG,COL,ROW,IMAGE(PTIMA)
cd390      FORMAT (I3,F5.1,2I5,F15.2)
cd         CALL STTPUT (STRING,STAT)

400   CONTINUE
C
C --- End of loop on table values
C
      GOTO 100		! Search for another column number
1000  CONTINUE		! Continue the process
      IF (INTCOL.GT.0) THEN
         WRITE (STRING,395) INTCOL
395      FORMAT(1X,I5,' columns interpolated by a polynomial')
         CALL STTPUT(STRING,STAT)
      ENDIF         
C
C --- 6. Interpolating along the rows
C
C --- 6.0. Initialisations
C
      THRES = 20.
      PASS1 = 0		! Indicates the first pass in the loop on ROW

      XB = 1.0		! Boundaries of the interpolation
      XE = FLOAT(NX)
C      SMOOTH = 0.	! Smoothing factor for this part of the algorithm

      NBUND = 0		! Initialises number of undefined columns
      DO 1050 COL = 1 , NX
         IF (TABFLG(COL,2).EQ.0) THEN
            NBUND = NBUND + 1
            UNDCOL (NBUND) = FLOAT(COL)
         ENDIF
1050  CONTINUE
      IF (NBUND.EQ.0) RETURN ! All columns are defined.
C                              That's no more work to do...

      DO 1100 ROW = 1 , NY

C
C --- 
C
         RATIO = 100.*FLOAT(ROW)/FLOAT(NY)
         IF (RATIO.GE.THRES) THEN
             WRITE (STRING,1150) THRES
1150         FORMAT ('Background interpolation:',F5.1,' % done')
             CALL STTPUT (STRING,STAT)
             THRES = THRES + 20.0
         ENDIF

         NPTS = 0
         DO 1200 COL = 1 , NX
C
C --- 6.1. Load values for the defined columns in arrays X,Y,W
C
           IF (TABFLG(COL,2).NE.0) THEN

              NPTS = NPTS + 1
              X(NPTS) = FLOAT (COL)
              PTIMA = (ROW-1)*NX+COL
              Y(NPTS) = IMAGE(PTIMA)
cd              W(NPTS) = 1.0		! Speed optim.
cd              TYPE*,NPTS,X(NPTS),Y(NPTS)
           ENDIF
1200     CONTINUE
C
C --- 6.2. Computes b-splines coefficients.
C

         IF (PASS1.EQ.0) THEN 	! On the first pass, ...
            IOPT = 0		! Restart all computations, ...
            PASS1 = 1
         ELSE			! then ...
            IOPT = 1		! starts with the knots found at the
         ENDIF			! last call.
         
         CALL SMOOT (X,Y,W,Q,NPTS,XB,XE,DEGREE,SMOOTH,NBKNOT,
     +                KNOTS,COEF,FP,IOPT,IER,NRDATA,FPO,FPOLD,NPLUS)
C
C --- 6.3. Error management
C
             IF (IER.EQ.3) THEN  !   No convergence achieved

             CALL SMOOT (X,Y,W,Q,NPTS,XB,XE,DEGREE,BIGSMO,NBKNOT,
     +                    KNOTS,COEF,FP,0,IER,NRDATA,FPO,FPOLD,NPLUS)
             PASS1 = 0
             INTROW = INTROW + 1

             ENDIF

         CALL IERSPL (IER,LEVEL)	! Displays error messages

         IF (IER.GT.0) THEN

              ERR = 1		! Error during spline interpolation
C5300          FORMAT ('Error while processing row ',I6)
              WRITE(STRING,5200) ROW
              CALL STTPUT(STRING,STAT)
              RETURN

         ENDIF
C
C --- 6.4. Computes interpolated values 
C          (Only where it is not already defined)
C
      NU  = 0				! Order of the derivative
      PREPOS = DEGREE + 1

C      DO 1400 NPTS = 1 , NBUND
C         ARG = UNDCOL (NPTS) ! The argument is an undefined column position
      DO 1400 ARG = 1.0 , NX
C
C --- 6.4.1. Find the position of the argument within the knots sequence
C
         CALL ARGSPL (ARG,KNOTS,NBKNOT,DEGREE,PREPOS,POSARG,NK1)
C
C --- 6.4.2. Computes interpolated values
C
C        ROW is already defined within the main loop
         COL = NINT(ARG)
         PTIMA = (ROW-1)*NX+COL
         IMAGE(PTIMA)  = DERIV (KNOTS,NBKNOT,COEF,NK1,NU,ARG,POSARG)

1400   CONTINUE		! End of the loop on ARG

1100   CONTINUE		! End of the loop on ROW

      IF (INTROW.GT.0) THEN
         WRITE (STRING,396) INTROW
396      FORMAT(1X,I5,' rows interpolated by a polynomial')
         CALL STTPUT(STRING,STAT)
      ENDIF         
C
C --- 7. CUTS computation
C
      CUTS(3) = IMAGE(1)	! Minimum value
      CUTS(4) = IMAGE(1)	! Maximum value

      DO 2000 PTIMA = 2 , NX*NY

       IF (IMAGE(PTIMA).LT.CUTS(3)) CUTS(3)=IMAGE(PTIMA)
       IF (IMAGE(PTIMA).GT.CUTS(4)) CUTS(4)=IMAGE(PTIMA)

2000  CONTINUE

      MEAN = (CUTS(3)+CUTS(4))/2.
      CUTS(3) = MEAN - (MEAN-CUTS(3))*1.0
      CUTS(4) = MEAN + (CUTS(4)-MEAN)*1.0

      CUTS(1) = CUTS(3)
      CUTS(2) = CUTS(4)
C
C --- 8. That's all folks          
C
      RETURN
      END










       FUNCTION MEASUR (IMAGE,NX,NY,XSTA,XEND,YSTA,YEND,WIDTH,MMODE)
C
C ----------------------------------------------------------------------------
C
C.  IDENTIFICATION
C
C       Author : P. Ballester    ESO - Garching
C
C       Versions :           1.00    30.06.90
C
C.  PURPOSE
C
C      This function measures on the frame the median value of the signal
C           on the neighbourhood of the position (POSX,POY), inside
C           a sub-frame of dimension (2*WIDX+1,2*WIDY+1)
C
C.  RESTRICTIONS : 
C
C       (2*WIDX+1)*(2*WIDY+1) < 225 (Parameter MAXPIX)
C       Executable only within MIDAS environment (STTPUT,STETER)
C
C.  INPUTS :
C
C       IMAGE       Frames on which the signal is measured
C       NX,NY       Dimensions of the frame
C       POSX,POSY   Central position of the sub-image
C       WIDX        X Half-size of the window
C       WIDY        Y Half-size of the window
C
C.  OUTPUTS :
C
C       MEASUR      Mean value of the signal in the sub-image
C
C.  ALGORITHM :
C
C      One checks that the central position is defined in the frame. Then
C          the pixels of the neighbourhood of the central position are
C          summed up (excepted the ones which are not defined in the
C          frame) and the median value is returned.
C
C ----------------------------------------------------------------------------
C
C --- 1. Variables definition
C
       IMPLICIT NONE
C
       LOGICAL FLAG
       INTEGER MAXPIX,NPIX,MEDMIN,MEDMAX,KC,MMODE
       PARAMETER (MAXPIX=1024)
       REAL    MEASUR,IMAGE(1),HISTO(MAXPIX),BUFFER
       INTEGER NX,NY,WIDTH
       INTEGER PTIMA,COL,ROW,POSY,KR, STP, NXSTA, NXEND
       REAL    XSTA,XEND,YSTA,YEND,SLOPE,INTER
       CHARACTER STRING*80
C
C --- 3. Computes mean value
C

       SLOPE  = (YEND-YSTA)/(XEND-XSTA)
       INTER  = ((YEND+YSTA)-SLOPE*(XEND+XSTA))/2.

       STP = 1
       IF (XSTA.GT.XEND) THEN 
           NXSTA = NINT(XEND)
           NXEND = NINT(XSTA)
       ELSE
           NXSTA = NINT(XSTA)
           NXEND = NINT(XEND)
       ENDIF

       MEASUR = 0.
       NPIX   = 0

       DO  100  COL = NXSTA, NXEND

       POSY = NINT(COL*SLOPE+INTER)

       DO  50  KR = -WIDTH,WIDTH

       ROW  = POSY + KR

       IF (COL.LT.1.OR.COL.GT.NX.OR.ROW.LT.1.OR.ROW.GT.NY) GOTO 50

       NPIX  = NPIX + 1

       IF (NPIX.GT.MAXPIX) THEN
           STRING = '*** Fatal : Buffer HISTO is too small in MEASUR'
           CALL STETER (9,STRING)
       ENDIF

       PTIMA = (ROW-1)*NX + COL

       HISTO(NPIX) = IMAGE (PTIMA)

50     CONTINUE
100    CONTINUE
C
C --- 4. Sort the values
C
      IF (NPIX.EQ.1) THEN
          MEASUR = HISTO(1)
          RETURN
      ENDIF

150   CONTINUE
      FLAG = (.FALSE.)
      DO 200 KC = 1, NPIX-1
         IF (HISTO(KC).GT.HISTO(KC+1)) THEN
             BUFFER = HISTO (KC+1)
             HISTO(KC+1) = HISTO(KC)
             HISTO(KC)   = BUFFER
             FLAG = (.TRUE.)
         ENDIF
200    CONTINUE
       IF (FLAG) GOTO 150
C
C --- 5.a) Return the MINIMUM if desired
C
       IF (MMODE.EQ.1) THEN
          MEASUR = HISTO(1)
          RETURN
       ENDIF
C
C --- 5.b) Define the median value
C
       IF (2*NINT(NPIX/2.).EQ.NPIX) THEN
            MEDMIN = NPIX/2
            MEDMAX = NPIX/2 + 1
            MEASUR = (HISTO(MEDMIN)+HISTO(MEDMAX))/2.
       ELSE
            MEDMIN = (NPIX-1)/2+1
            MEASUR = HISTO(MEDMIN)
       ENDIF

C       write(6,*) 'made it till here and measured ', MEASUR

       RETURN
       END













