C @(#)tdnlrb.for	17.1.1.1 (ESO-DMD) 01/25/02 17:47:15 
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.COPYRIGHT: Copyright (c) 1987 European Southern Observatory,
C                                         all rights reserved
C
C.VERSION: 1.0  ESO-FORTRAN Conversion, AA  8:58 - 9 DEC 1987
C
C.LANGUAGE: F77+ESOext
C
C.IDENTIFICATION
C
C  SUBROUTINES TDRBII, TDRBIT, TDRBTI, TDRBTT
C
C.KEYWORDS
C
C  IMAGES, NONLINEAR REBIN, INTERPOLATION, INTEGRATION
C
C.PURPOSE
C
C  NONLINEAR REBIN IMAGE TO IMAGE IN 1 DIMENSION
C
C.ALGORITHM
C
C  Rebin input data YI sampled at XI with binwidth WI to
C  output data YO sampled at XO with binwidth WO.
C
C  Restriction:  Input and output independent variables are
C                assumed to be monotonic increasing or decreasing
C
C  Several FUNCTal relations between input and output independent
C  variables are foreseen, additional ones may be added in dummy call
C  U01 as follows:
C
C SUBROUTINE REB_U01 (DV,P,NP,NPC)
C   REAL*8 DV,P(NP)
C   DV = ...........
C   RETURN
C END
C
C  where DV is input-output value, P(NP) the parameter array and
C  NPC the nuber of active parameters in p(np).
C
C  According to the definition of IMAGES, pixels have constant
C  bin width (STEP) [WI;WO], are not overlapping and have filling
C  factors of exactly one. World coordinates [XI;XO] and flux values
C  [YI;YO] are pertinent
C  to the center of a pixel. In flux conservation mode (the default)
C  the flux is interpreted as the area under a smooth curve between
C  the boundaries of the pixel (i.e. a histogram representation of
C  a smooth curve).
C
C  Detailed flux conservation obtained by REAL*8 spline interpolation
C  and NAGLIB REAL*8 integration of input data over the projected
C  bin WI for fractions of input pixels and summation of flux in
C  entire input pixels covered within the valid domain of input data
C  i.e. from x(1)-0.5*step to x(n)+0.5*step. Several options available
C  to handle undefined input at the extremes and for NULLVALUE assigned
C  parts of input array.
C
C   Optimized for long arrays by search for valid range, optimum
C  direction or DO LOOPS and updating of loop boundaries - therefore
C  overdoing it for small arrays.
C
C.INPUT
C
C  Keys:  IN_A/C/1/8         input data array
C         IN_B/C/1/8         reference frame or dim of output
C         CFUNC/C/1/8        FUNCTal relation between indep.
C         INPUTD/D/1/12      parameters or FUNCT
C         COPT/C/1/8         options
C
C.OUTPUT
C
C  Key:   OUT_A/C/1/8        output data array
C 
C.MODIF 9-03-1992  Correct option COPT MP
C-----------------------------------------------------------
C
      SUBROUTINE TDRBII
      IMPLICIT NONE
      INTEGER   NFPAR
      PARAMETER (NFPAR = 12)
C
      INTEGER MADRID
      INTEGER I,IACT, IS,ISTAT
      INTEGER J,NAXISI,NAXISO,KUN,KNUL
      INTEGER NDIMI,NDIMO,NFPACT,NFUNC
      INTEGER NI,NINT,NN,NO,NBY,NBO
      INTEGER NPIXI(3),NPIXO(3),INOP,IMNO,IMNI,IREF
C 
      INTEGER*8   IP1,IP2,IP3,IPNTR, JP1,JP2,JPNTR
C 
      REAL    RMIN, RMAX, CUTS(4)
      REAL    X0,WX
C 
      DOUBLE PRECISION STEPI(3),STARTI(3)
      DOUBLE PRECISION DPARM(NFPAR),FPARM(NFPAR)
      DOUBLE PRECISION DSTART(3), DSTEP(3)
C 
      CHARACTER*3 FUN(9)
      CHARACTER*80 FUNCT,OPTION
      CHARACTER*60 INPIMA,OUTIMA,REFIMA
      CHARACTER*72 IDENTO,IDENTI
      CHARACTER*80 CUNITO,CUNITI
C
      INCLUDE 'MID_INCLUDE:TABLES.INC'
      COMMON /VMR/MADRID(1)
      INCLUDE 'MID_INCLUDE:TABLED.INC'
C
      DATA FUN/'LIN','POL','INV','EXP','DEX','LOG','DLG','IPO','U01'/
      DATA INPIMA/'  '/
      DATA OUTIMA/'  '/
      DATA REFIMA/'  '/
      DATA IDENTI/'  '/
      DATA CUNITI/'  '/
C
C ... get input parameters
C
      CALL STKRDC('OUT_A',1,1,60,IACT,OUTIMA,KUN,KNUL,ISTAT)
      CALL STKRDC('IN_A',1,1,60,IACT,INPIMA,KUN,KNUL,ISTAT)
      CALL STKRDC('IN_B',1,1,60,IACT,REFIMA,KUN,KNUL,ISTAT)
      CALL STKRDC('CFUNC',1,1,8,IACT,FUNCT,KUN,KNUL,ISTAT)
      CALL STKRDD('INPUTD',1,NFPAR,NFPACT,FPARM,KUN,KNUL,ISTAT)
      CALL STKRDC('COPT',1,1,8,IACT,OPTION,KUN,KNUL,ISTAT)
C
C ... det. interpolation/integration options
C  !   Not implemented
      INOP   = 4
C      IF (OPTION(1:3).EQ.'PIX') INOP   = 1
C      IF (OPTION(1:3).EQ.'LIN') INOP   = 2
C      IF (OPTION(1:3).EQ.'SPG') INOP   = 3
C .......................................... the actual values are:
      CALL FORUPC(OPTION,OPTION)
      IF (OPTION(1:1).EQ.'P') INOP   = 2
      IF (OPTION(1:1).EQ.'L') INOP   = 3
      IF (OPTION(1:1).EQ.'S') INOP   = 1
C .......................................... Bitte, bitte
C
C ... translate FUNCT into FUNCT number
C
      CALL FORUPC(FUNCT,FUNCT)
      NFUNC  = 0
      DO 10 I = 1,9
          IF (FUNCT(1:3).EQ.FUN(I)) NFUNC  = I
   10 CONTINUE
      IF (NFUNC.EQ.0) THEN
          CALL STTPUT(' Specified FUNCT non-existent...',ISTAT)
          GO TO 60
      END IF
C
C ... search for number of active parameters and DOUBLE
C
      DO 20 J = NFPAR,1,-1
          NFPACT = J
          IF (FPARM(J).NE.0.000) GO TO 30
   20 CONTINUE
   30 DO 40 J = 1,NFPACT
          DPARM(J) = FPARM(J)
   40 CONTINUE
C
C ... get parameters form reference image
C
      CALL STFOPN(REFIMA,D_R4_FORMAT,0,F_IMA_TYPE,IREF,ISTAT)
      CALL STDRDI(IREF,'NAXIS',1,1,NN,NAXISO,KUN,KNUL,ISTAT)
      CALL STDRDI(IREF,'NPIX',1,NAXISO,NN,NPIXO,KUN,KNUL,ISTAT)
      CALL STDRDD(IREF,'START',1,NAXISO,NN,DSTART,KUN,KNUL,ISTAT)
      CALL STDRDD(IREF,'STEP',1,NAXISO,NN,DSTEP,KUN,KNUL,ISTAT)
      CALL STDRDC(IREF,'IDENT',1,1,72,NN,IDENTO,KUN,KNUL,ISTAT)
      CALL STDRDC(IREF,'CUNIT',1,1,80,NN,CUNITO,KUN,KNUL,ISTAT)
      IF (NAXISO.GT.1 .AND. NPIXO(2).GT.1) THEN
          CALL STTPUT(' Only 1D images are considered ...',ISTAT)
          CALL STFCLO(IREF,ISTAT)
          GO TO 60
      ELSE
          NPIXO(2) = 1
          NDIMO  = 1
      END IF
      CALL STFCLO(IREF,ISTAT)
C
C ... get input image
C
      CALL STIGET(INPIMA,D_R4_FORMAT,F_I_MODE,F_IMA_TYPE,
     .2,NAXISI,NPIXI,STARTI,STEPI,IDENTI,CUNITI,IPNTR,IMNI,ISTAT)
      IF (NAXISI.GT.1 .AND. NPIXI(2).GT.1) THEN
          CALL STTPUT(' Only 1-D images are considered.',ISTAT)
          GO TO 50
      ELSE
          NPIXI(2) = 1
          NDIMI  = 1
      END IF
C
C ... map output image
C
      CALL STIPUT(OUTIMA,D_R4_FORMAT,F_O_MODE,F_IMA_TYPE,
     .NDIMO,NPIXO,DSTART,DSTEP,IDENTI,CUNITO,JPNTR,IMNO,ISTAT)
C
C ... fill the arrays XI,WI and XO,WO
C
      NI     = NPIXI(1)
      NBY    = 8*NI
      CALL TDMGET(NBY,IP1,IS)
      CALL TDMGET(NBY,IP2,IS)
      CALL TDMGET(NBY,IP3,IS)
      X0     = STARTI(1)
      WX     = STEPI(1)
      CALL IMVAL3(NI,X0,WX,MADRID(IPNTR),MADRID(IP1),MADRID(IP2),
     +             MADRID(IP3))
C
      NO     = NPIXO(1)
      NBO    = 8*NO
      CALL TDMGET(NBO,JP1,IS)
      CALL TDMGET(NBO,JP2,IS)
      X0     = DSTART(1)
      WX     = DSTEP(1)
      CALL IMVAL2(NO,X0,WX,MADRID(JP1),MADRID(JP2))
C
C  ...  Now, at last, we call the interpolating (and integrating) routines:
C
C  !!!! Watch out for VECI,VECD arrays in INTERP, dimension => NINT !!!!!!
C
      NINT   = 8
      CALL REBMET(NI,MADRID(IP1),MADRID(IP3),MADRID(IP2),NO,
     +             MADRID(JP1),MADRID(JP2),NFUNC,NFPAR,NFPACT,DPARM,
     +             INOP,NINT,MADRID(JPNTR),RMIN,RMAX)
C
C ... write cuts
C
      CUTS(1) = RMIN
      CUTS(2) = RMAX
      CUTS(3) = RMIN
      CUTS(4) = RMAX
      CALL STDWRR(IMNO,'LHCUTS',CUTS,1,4,KUN,ISTAT)
C
C ... free memory
C
      CALL TDMFRE(NBY,IP1,IS)
      CALL TDMFRE(NBY,IP2,IS)
      CALL TDMFRE(NBY,IP3,IS)
      CALL TDMFRE(NBO,JP1,IS)
      CALL TDMFRE(NBO,JP2,IS)
C
C ... end
C
      CALL DSCUPT(IMNO,IMNO,' ',ISTAT)
      CALL STFCLO(IMNO,ISTAT)
   50 CALL STFCLO(IMNI,ISTAT)
   60 RETURN
      END

C
      SUBROUTINE TDRBIT
      IMPLICIT NONE
      INTEGER    NFPAR
      PARAMETER (NFPAR = 12)
C
      LOGICAL BIN
      INTEGER MADRID
      INTEGER*8  IP1,IP2,IP3,IPNTR, JP1,JP2,JP3
      INTEGER I,IACT,IC,II
      INTEGER IS,ISTAT,J,JSTAT,TID,NAC,NAR
      INTEGER NAXIS,NCOL,NDIM,NFUNC,NI
      INTEGER NINT,NMAX,NO,NROW,NSC,NBY,NBY1
      INTEGER INDEX, KUN, KNUL
      INTEGER NFPACT,INOP,IMNO
      INTEGER ICOL(3),NPIX(3)
C
      REAL RMIN,RMAX
      REAL X0,WX
C
      DOUBLE PRECISION  STEP(3),START(3)
C
      DOUBLE PRECISION DPARM(NFPAR),FPARM(NFPAR)
      CHARACTER*3 FUN(9)
      CHARACTER*80 FUNCT1,OPTION
      CHARACTER*60 TABLE,IMAGE
      CHARACTER*16 MSG
      CHARACTER*17 COLREF(3),LABEL
      CHARACTER*53 OCOLS,CBUF
      CHARACTER*72 IDENT
      CHARACTER*80 CUNIT
C
      INCLUDE 'MID_INCLUDE:TABLES.INC'
      COMMON /VMR/MADRID(1)
      INCLUDE 'MID_INCLUDE:TABLED.INC'
C
      DATA MSG/'ERR:REBISTBxxx'/
      DATA FUN/'LIN','POL','INV','EXP','DEX','LOG','DLG','IPO','U01'/
      DATA TABLE/'  '/
      DATA IMAGE/'  '/
C
C ... get input parameters
C
      CALL STKRDC('OUT_A',1,1,60,IACT,TABLE,KUN,KNUL,ISTAT)
      CALL STKRDC('COLO',1,1,53,IACT,OCOLS,KUN,KNUL,ISTAT)
      CALL STKRDC('IN_A',1,1,60,IACT,IMAGE,KUN,KNUL,ISTAT)
      CALL STKRDC('CFUNC',1,1,8,IACT,FUNCT1,KUN,KNUL,ISTAT)
      CALL STKRDD('INPUTD',1,NFPAR,NFPACT,FPARM,KUN,KNUL,ISTAT)
      CALL STKRDC('COPT',1,1,12,IACT,OPTION,KUN,KNUL,ISTAT)
C
C ... det. integration/interpolation option
C  !   Not implemented
      INOP   = 4
C      IF (OPTION(1:3).EQ.'PIX') INOP   = 1
C      IF (OPTION(1:3).EQ.'LIN') INOP   = 2
C      IF (OPTION(1:3).EQ.'SPG') INOP   = 3
C .......................................... the actual values are:
      CALL FORUPC(OPTION,OPTION)
      IF (OPTION(1:1).EQ.'P') INOP   = 2
      IF (OPTION(1:1).EQ.'L') INOP   = 3
      IF (OPTION(1:1).EQ.'S') INOP   = 1
C .......................................... Bitte, bitte
C
C ... translate FUNCT1 into function number
C
      CALL FORUPC(FUNCT1,FUNCT1)
      NFUNC  = 0
      DO 10 I = 1,9
          IF (FUNCT1(1:3).EQ.FUN(I)) NFUNC  = I
   10 CONTINUE
      IF (NFUNC.EQ.0) THEN
          CALL STTPUT(' Specified function non-existent...',ISTAT)
          GO TO 70
      END IF
C
C ... search for number of active parameters and 'double-precision' them
C
      DO 20 J = NFPAR,1,-1
          NFPACT = J
          IF (FPARM(J).NE.0.000) GO TO 30
   20 CONTINUE
   30 DO 40 J = 1,NFPACT
          DPARM(J) = FPARM(J)
   40 CONTINUE
C
C ... sort out ind., dep. (and bin) variables
C
      II     = INDEX(OCOLS,',')
      IF (II.EQ.0) THEN
          WRITE (*,*) ' Specification of relevant input',
     +      ' table cols insufficient'
          GOTO 70
      ELSE
          COLREF(1) = OCOLS(1:II-1)
          CBUF   = OCOLS(II+1:)
          II     = INDEX(CBUF,',')
          IF (II.EQ.0) THEN
              COLREF(2) = CBUF(1:)
              BIN    = .FALSE.
          ELSE
              COLREF(2) = CBUF(1:II-1)
              COLREF(3) = CBUF(II+1:)
              BIN    = .TRUE.
          END IF
      END IF
C
C ... access input table
C
      CALL TBTOPN(TABLE,F_U_MODE,TID,ISTAT)
      IF (ISTAT.NE.0) GO TO 60
      CALL TBIGET(TID,NCOL,NROW,NSC,NAC,NAR,ISTAT)
      IF (ISTAT.NE.0) GO TO 60
C
C ... find input column
C
      CALL TBCSER(TID,COLREF(1),ICOL(1),ISTAT)
      IF (ICOL(1).EQ.-1) THEN
          CALL STTPUT('Error in column reference',IS)
          GO TO 60
      END IF
      IF (BIN) THEN
          CALL TBCSER(TID,COLREF(3),ICOL(3),ISTAT)
          IF (ICOL(3).EQ.-1) THEN
              CALL STTPUT('Error in column reference',IS)
              GO TO 60
          END IF
      END IF
C
C ... copy selected values of the table into work space
C
      NBY = 8*NROW
      CALL TDMGET(NBY,JP1,IS)
      CALL TDMGET(NBY,JP2,IS)
      CALL TDMGET(NBY,JP3,IS)
      IF (BIN) THEN
          CALL COPYIT(TID,ICOL(1),ICOL(3),NROW,
     +                MADRID(JP1),MADRID(JP3),NO)
          IF (NO.LE.3) THEN
              CALL STTPUT(' Not enough points in table ',ISTAT)
              GO TO 50
          END IF
      ELSE
          CALL BCPYIT(TID,ICOL(1),NROW,MADRID(JP1),
     +                 MADRID(JP3),NO)
          IF (NO.LE.3) THEN
              CALL STTPUT(' Not enough points in table ',ISTAT)
              GO TO 50
          END IF
      END IF
C
C ... get input image
C
      CALL STIGET(IMAGE,D_R4_FORMAT,F_I_MODE,F_IMA_TYPE,
     .     2,NAXIS,NPIX,START,STEP,IDENT,CUNIT,IPNTR,IMNO,ISTAT)
      IF (NAXIS.GT.1 .AND. NPIX(2).GT.1) THEN
          CALL STTPUT(' Only 1-dim images are considered.',ISTAT)
          GO TO 60
      ELSE
          NPIX(2) = 1
          NDIM   = 1
      END IF
C
C image input, then fill the XO and WO arrays OF INTERPOLn
C
      NI     = NPIX(1)
      NBY1 = 8*NI
      CALL TDMGET(NBY1,IP1,IS)
      CALL TDMGET(NBY1,IP2,IS)
      CALL TDMGET(NBY1,IP3,IS)
      X0     = START(1)
      WX     = STEP(1)
      CALL IMVAL4(NI,X0,WX,MADRID(IPNTR),MADRID(IP1),MADRID(IP2),
     +             MADRID(IP3))
C
C  ...  Now, at last, we call the interpolating (and integrating)
C	routines:
C  !!!! Watch out for VECI,VECD arrays in INTERP,
C	 dimension => NINT !!!!!!
C
      CALL REBMET(NI,MADRID(IP1),MADRID(IP2),MADRID(IP3),NO,
     +             MADRID(JP1),MADRID(JP3),NFUNC,NFPAR,NFPACT,DPARM,
     +             INOP,NINT,MADRID(JP2),RMIN,RMAX)
C
C ... output table
C
C ... find if column already exists
      CALL TBCSER(TID,COLREF(2),IC,ISTAT)
      IF (IC.EQ.-1) THEN
          LABEL  = COLREF(2) (2:)
          CALL TBCINI(TID,D_R4_FORMAT,1,'E16.8',CUNIT,
     .                LABEL,IC,ISTAT)
      END IF
      NMAX   = NO
      CALL OCOPY(TID,NMAX,IC,MADRID(JP2))
C
C ... free memory
C
      CALL TDMFRE(NBY,JP1,IS)
      CALL TDMFRE(NBY,JP2,IS)
      CALL TDMFRE(NBY,JP3,IS)
   50 CONTINUE
      CALL TDMFRE(NBY1,IP1,IS)
      CALL TDMFRE(NBY1,IP2,IS)
      CALL TDMFRE(NBY1,IP3,IS)
C
C ... end
C
      CALL TBIPUT(TID,0,NO,ISTAT)
      CALL TBSINI(TID,ISTAT)
      CALL DSCUPT(TID,TID,' ',ISTAT)
      CALL TBTCLO(TID,ISTAT)
   60 IF (ISTAT.NE.0) THEN
          WRITE (MSG(13:16),9000) ISTAT
          CALL TDERRR(ISTAT,MSG,JSTAT)
      END IF
   70 RETURN
 9000 FORMAT (I4)
      END


      SUBROUTINE TDRBTI
      IMPLICIT NONE
      INTEGER    NFPAR
      PARAMETER (NFPAR = 12)
C
      LOGICAL BIN
      INTEGER MADRID
      INTEGER*8  IP1,IP2,IP3,IPNTR, JP1,JP2
      INTEGER I,I1,I2,IACT,IB,II
      INTEGER IS,ISTAT,J,JSTAT,NAXIS
      INTEGER NCOL,NDIM,NFUNC,NI,NINT,NN,NO,NROW,NSC
      INTEGER INDEX,IMNO,IREF, KUN, KNUL
      INTEGER NFPACT,INOP,NAC,NAR,TID,NBY,NBY1
      INTEGER ICOL(3),NPIX(3)
      REAL    RMIN,RMAX
      REAL    CUTS(4)
      REAL    X0,WX
      DOUBLE PRECISION DPARM(NFPAR),FPARM(NFPAR)
      DOUBLE PRECISION DSTEP(3), DSTART(3)
      CHARACTER*3 FUN(9)
      CHARACTER*8 FUNCTN
      CHARACTER*80 OPTION
      CHARACTER*60 TABLE,IMAGE,REFIMA
      CHARACTER*16 MSG
      CHARACTER*17 COLREF(3)
      CHARACTER*53 INCOLS,CBUF
      CHARACTER*60 WORK
      CHARACTER*72 IDENT
      CHARACTER*80 CUNIT
C
      INCLUDE 'MID_INCLUDE:TABLES.INC'
      COMMON /VMR/MADRID(1)
      INCLUDE 'MID_INCLUDE:TABLED.INC'
C
      DATA MSG/'ERR:REBISTxxxx'/
      DATA FUN/'LIN','POL','INV','EXP','DEX','LOG','DLG','IPO','U01'/

      DATA TABLE/'  '/
      DATA IMAGE/'  '/
      DATA REFIMA/'  '/
C
C ... get input parameters
C
      CALL STKRDC('OUT_A',1,1,60,IACT,IMAGE,KUN,KNUL,ISTAT)
      CALL STKRDC('IN_A',1,1,60,IACT,TABLE,KUN,KNUL,ISTAT)
      CALL STKRDC('IN_B',1,1,60,IACT,REFIMA,KUN,KNUL,ISTAT)
      CALL STKRDC('COLI',1,1,53,IACT,INCOLS,KUN,KNUL,ISTAT)
      CALL STKRDC('CFUNC',1,1,8,IACT,FUNCTN,KUN,KNUL,ISTAT)
      CALL STKRDD('INPUTD',1,NFPAR,NFPACT,FPARM,KUN,KNUL,ISTAT)
      CALL STKRDC('COPT',1,1,12,IACT,OPTION,KUN,KNUL,ISTAT)
C
C ... determine integration/interpolation option
C  !   Not implemented
      INOP   = 4
C      IF (OPTION(1:3).EQ.'PIX') INOP   = 1
C      IF (OPTION(1:3).EQ.'LIN') INOP   = 2
C      IF (OPTION(1:3).EQ.'SPG') INOP   = 3
C .......................................... the actual values are:
      CALL FORUPC(OPTION,OPTION)
      IF (OPTION(1:1).EQ.'P') INOP   = 2
      IF (OPTION(1:1).EQ.'L') INOP   = 3
      IF (OPTION(1:1).EQ.'S') INOP   = 1
C .......................................... Bitte, bitte
C
C ... translate FUNCTN into function number
C
      CALL FORUPC(FUNCTN,FUNCTN)
      NFUNC  = 0
      DO 10 I = 1,9
          IF (FUNCTN(1:3).EQ.FUN(I)) NFUNC  = I
   10 CONTINUE
      IF (NFUNC.EQ.0) THEN
          CALL STTPUT(' Specified function non-existent...',ISTAT)
          GO TO 80
      END IF
C
C ... search for number of active parameters and 'double-precision' them
C
      DO 20 J = NFPAR,1,-1
          NFPACT = J
          IF (FPARM(J).NE.0.000) GO TO 30
   20 CONTINUE
   30 DO 40 J = 1,NFPACT
          DPARM(J) = FPARM(J)
   40 CONTINUE
C
C ... get parameters from reference image
C
      CALL STFOPN(REFIMA,D_R4_FORMAT,0,F_IMA_TYPE,IREF,ISTAT)
      CALL STDRDI(IREF,'NAXIS',1,1,NN,NAXIS,KUN,KNUL,ISTAT)
      CALL STDRDI(IREF,'NPIX',1,NAXIS,NN,NPIX,KUN,KNUL,ISTAT)
      CALL STDRDD(IREF,'START',1,NAXIS,NN,DSTART,KUN,KNUL,ISTAT)
      CALL STDRDD(IREF,'STEP',1,NAXIS,NN,DSTEP,KUN,KNUL,ISTAT)
      CALL STDRDC(IREF,'IDENT',1,1,72,NN,IDENT,KUN,KNUL,ISTAT)
      CALL STDRDC(IREF,'CUNIT',1,1,80,NN,CUNIT,KUN,KNUL,ISTAT)
      IF (NAXIS.GT.1 .AND. NPIX(2).GT.1) THEN
          CALL STTPUT(' Only 1D images are considered ...',ISTAT)
          CALL STFCLO(IREF,ISTAT)
          GO TO 70
      ELSE
          NPIX(2) = 1
          NDIM   = 1
      END IF
      CALL STFCLO(IREF,ISTAT)
C
C ... sort out ind., dep. (and bin) variables
C
      II     = INDEX(INCOLS,',')
      IF (II.EQ.0) THEN
          WRITE (*,*) ' Specification of relevant input',
     +      ' table cols insufficient'
          GOTO 70
      ELSE
          COLREF(1) = INCOLS(1:II-1)
          CBUF   = INCOLS(II+1:)
          II     = INDEX(CBUF,',')
          IF (II.EQ.0) THEN
              COLREF(2) = CBUF(1:)
              BIN    = .FALSE.
          ELSE
              COLREF(2) = CBUF(1:II-1)
              COLREF(3) = CBUF(II+1:)
              BIN    = .TRUE.
          END IF
      END IF
C
C ... init input table
C
      CALL TBTOPN(TABLE,F_I_MODE,TID,ISTAT)
      IF (ISTAT.NE.0) GO TO 70
      CALL TBIGET(TID,NCOL,NROW,NSC,NAC,NAR,ISTAT)
      IF (ISTAT.NE.0) GO TO 70
C
C ... find input column
C
      IB     = 2
      IF (BIN) IB     = 3
      DO 50 I = 1,IB
          CALL TBCSER(TID,COLREF(I),ICOL(I),ISTAT)
          IF (ICOL(I).EQ.-1) THEN
              CALL STTPUT('Error in column reference',IS)
              GO TO 70
          END IF
C
C ... copy units from table into output image
C
          IF (I.LE.2) THEN
              CALL TBUGET(TID,ICOL(I),WORK,ISTAT)
              I1     = (I-1)*16 + 1
              I2     = I1 + 15
              CUNIT(I1:I2) = WORK(1:16)
          END IF
   50 CONTINUE
C
C ... copy selected values of the table into work space
C
      NBY  = 8*NROW
      CALL TDMGET(NBY,IP1,IS)
      CALL TDMGET(NBY,IP2,IS)
      CALL TDMGET(NBY,IP3,IS)
      IF (BIN) THEN
          CALL BCOPYC(TID,ICOL(1),ICOL(2),ICOL(3),NROW,
     .         MADRID(IP1),MADRID(IP2),MADRID(IP3),NI)
          IF (NI.LE.3) THEN
              CALL STTPUT(' Not enough points in table ',ISTAT)
              GO TO 60
          END IF
      ELSE
          CALL BCOPYI(TID,ICOL(1),ICOL(2),NROW,
     .                MADRID(IP1),MADRID(IP2),MADRID(IP3),NI)
          IF (NI.LE.3) THEN
              CALL STTPUT(' Not enough points in table ',ISTAT)
              GO TO 60
          END IF
      END IF
C
C ... map image
C
      CUNIT(I2+1:I2+16) = 'FLUX'
      IDENT  = 'TABLE: '//TABLE//'COLS. :'//COLREF(1)//COLREF(2)
      CALL STIPUT(IMAGE,D_R4_FORMAT,F_O_MODE,F_IMA_TYPE,
     .NDIM,NPIX,DSTART,DSTEP,IDENT,CUNIT,IPNTR,IMNO,ISTAT)
C
C image output, then fill the XO and WO arrays OF INTERPOLn
C
      NO     = NPIX(1)
      NBY1  = 8*NO
      CALL TDMGET(NBY1,JP1,IS)
      CALL TDMGET(NBY1,JP2,IS)
      X0     = DSTART(1)
      WX     = DSTEP(1)
      CALL IMVALS(NO,X0,WX,MADRID(JP1),MADRID(JP2))
C
C ...  Now, at last, we call the interpolating (and integrating)
C ...  routines:
C !!!! Watch out for VECI,VECD arrays in INTERP,
C                                  dimension => NINT !!!!!!
C
      CALL REBMET(NI,MADRID(IP1),MADRID(IP2),MADRID(IP3),NO,
     +             MADRID(JP1),MADRID(JP2),NFUNC,NFPAR,NFPACT,DPARM,
     +             INOP,NINT,MADRID(IPNTR),RMIN,RMAX)
C
C ... write cuts
C
      CUTS(1) = RMIN
      CUTS(2) = RMAX
      CUTS(3) = RMIN
      CUTS(4) = RMAX
      CALL STDWRR(IMNO,'LHCUTS',CUTS,1,4,KUN,ISTAT)
      CALL DSCUPT(IMNO,IMNO,' ',ISTAT)
C
C ... free memory
C
      CALL TDMFRE(NBY1,JP1,IS)
      CALL TDMFRE(NBY1,JP2,IS)
   60 CONTINUE
      CALL TDMFRE(NBY,IP1,IS)
      CALL TDMFRE(NBY,IP2,IS)
      CALL TDMFRE(NBY,IP3,IS)
C
C ... end
C
      CALL TBTCLO(TID,ISTAT)
   70 IF (ISTAT.NE.0) THEN
          WRITE (MSG(13:16),9000) ISTAT
          CALL TDERRR(ISTAT,MSG,JSTAT)
      END IF
   80 RETURN
 9000 FORMAT (I4)
      END

C
      SUBROUTINE TDRBTT
      IMPLICIT NONE
      INTEGER MADRID
      INTEGER    NFPAR
      PARAMETER (NFPAR = 12)
C
      INTEGER*8  IP1,IP2,IP3,JP1, JP2,JP3
      INTEGER I,IACT,IB,IC,II,IS,ISTAT
      INTEGER J,JSTAT,NCO,NCOL
      INTEGER NFUNC,NI,NINT,NMAX,NO,NRO,NROW,NSC
      INTEGER TID, TID1, NAC,NAR,NBY,NBY1
      INTEGER INDEX, KUN, KNUL
      INTEGER NFPACT,INOP
      INTEGER ICOL(3),OCOL(3)
C
      REAL RMAX,RMIN
C
      LOGICAL BIN
C
      DOUBLE PRECISION DPARM(NFPAR)
      DOUBLE PRECISION FPARM(NFPAR)
C
      CHARACTER*3  FUN(9)
      CHARACTER*8  FUNCTN
      CHARACTER*80 OPTION
      CHARACTER*60 TABIN,TABOUT
      CHARACTER*16 MSG
      CHARACTER*17 COLRFI(3),COLRFO(3),LABEL
      CHARACTER*53 INCOLS,OUTCOL,CBUF
      CHARACTER*80 CUNIT
C
      INCLUDE 'MID_INCLUDE:TABLES.INC'
      COMMON /VMR/MADRID(1)
      INCLUDE 'MID_INCLUDE:TABLED.INC'
      DATA MSG/'ERR:REBISTBxxx'/
      DATA FUN/'LIN','POL','INV','EXP','DEX','LOG','DLG','IPO','U01'/

      DATA TABIN/'  '/
      DATA TABOUT/'  '/
C
C ... get input parameters
C
      CALL STKRDC('OUT_A',1,1,60,IACT,TABOUT,KUN,KNUL,ISTAT)
      CALL STKRDC('COLO',1,1,53,IACT,OUTCOL,KUN,KNUL,ISTAT)
      CALL STKRDC('IN_A',1,1,60,IACT,TABIN,KUN,KNUL,ISTAT)
      CALL STKRDC('COLI',1,1,53,IACT,INCOLS,KUN,KNUL,ISTAT)
      CALL STKRDC('CFUNC',1,1,8,IACT,FUNCTN,KUN,KNUL,ISTAT)
      CALL STKRDD('INPUTD',1,NFPAR,NFPACT,FPARM,KUN,KNUL,ISTAT)
      CALL STKRDC('COPT',1,1,12,IACT,OPTION,KUN,KNUL,ISTAT)
C
C ... determine integration/interpolation option
C  !   Not implemented
      INOP   = 4
C      IF (OPTION(1:3).EQ.'PIX') INOP   = 1
C      IF (OPTION(1:3).EQ.'LIN') INOP   = 2
C      IF (OPTION(1:3).EQ.'SPG') INOP   = 3
C .......................................... the actual values are:
      CALL FORUPC(OPTION,OPTION)
      IF (OPTION(1:1).EQ.'P') INOP   = 2
      IF (OPTION(1:1).EQ.'L') INOP   = 3
      IF (OPTION(1:1).EQ.'S') INOP   = 1
C .......................................... Bitte, bitte
C
C ... translate FUNCTN into function number
C
      CALL FORUPC(FUNCTN,FUNCTN)
      NFUNC  = 0
      DO 10 I = 1,9
          IF (FUNCTN(1:3).EQ.FUN(I)) NFUNC  = I
   10 CONTINUE
      IF (NFUNC.EQ.0) THEN
          CALL STTPUT(' Specified function non-existent...',ISTAT)
          GO TO 80
      END IF
C
C ... search for number of active parameters and 'double-precision' them
C
      DO 20 J = NFPAR,1,-1
          NFPACT = J
          IF (FPARM(J).NE.0.000) GO TO 30
   20 CONTINUE
   30 DO 40 J = 1,NFPACT
          DPARM(J) = FPARM(J)
   40 CONTINUE
C
C ... sort out ind., dep. (and bin) variables
C
      II     = INDEX(INCOLS,',')
      IF (II.EQ.0) THEN
          WRITE (*,*) ' Specification of relevant input',
     +      ' table cols insufficient'
          GOTO 70
      ELSE
          COLRFI(1) = INCOLS(1:II-1)
          CBUF   = INCOLS(II+1:)
          II     = INDEX(CBUF,',')
          IF (II.EQ.0) THEN
              COLRFI(2) = CBUF(1:)
              BIN    = .FALSE.
          ELSE
              COLRFI(2) = CBUF(1:II-1)
              COLRFI(3) = CBUF(II+1:)
              BIN    = .TRUE.
          END IF
      END IF
C
C ... init input table
C
      CALL TBTOPN(TABIN,F_I_MODE,TID,ISTAT)
      IF (ISTAT.NE.0) GO TO 70
      CALL TBIGET(TID,NCOL,NROW,NSC,NAC,NAR,ISTAT)
      IF (ISTAT.NE.0) GO TO 70
C
C ... find input column
C
      IB     = 2
      IF (BIN) IB     = 3
      DO 50 I = 1,IB
          CALL TBCSER(TID,COLRFI(I),ICOL(I),ISTAT)
          IF (ICOL(I).EQ.-1) THEN
              CALL STTPUT('Error in column reference',IS)
              GO TO 70
          END IF
   50 CONTINUE
C
C ... copy selected values of the table into work space
C
      NBY = 8*NROW
      CALL TDMGET(NBY,IP1,IS)
      CALL TDMGET(NBY,IP2,IS)
      CALL TDMGET(NBY,IP3,IS)
      IF (BIN) THEN
          CALL BCOPYC(TID,ICOL(1),ICOL(2),ICOL(3),NROW,
     .         MADRID(IP1),MADRID(IP2),MADRID(IP3),NI)
          IF (NI.LE.3) THEN
              CALL STTPUT(' Not enough points in table ',ISTAT)
              GO TO 60
          END IF
      ELSE
          CALL BCOPYI(TID,ICOL(1),ICOL(2),NROW,
     .         MADRID(IP1),MADRID(IP2),MADRID(IP3),NI)
          IF (NI.LE.3) THEN
              CALL STTPUT(' Not enough points in table ',ISTAT)
              GO TO 60
          END IF
      END IF
C
C ... sort out ind., dep. (and bin) variables for output table
C
      II     = INDEX(OUTCOL,',')
      IF (II.EQ.0) THEN
          WRITE (*,*) ' Specification of relevant output',
     +      ' table cols insufficient'
          GOTO 70
      ELSE
          COLRFO(1) = OUTCOL(1:II-1)
          CBUF   = OUTCOL(II+1:)
          II     = INDEX(CBUF,',')
          IF (II.EQ.0) THEN
              COLRFO(2) = CBUF(1:)
              BIN    = .FALSE.
          ELSE
              COLRFO(2) = CBUF(1:II-1)
              COLRFO(3) = CBUF(II+1:)
              BIN    = .TRUE.
          END IF
      END IF
C
C ... access output table
C
      CALL TBTOPN(TABOUT,F_U_MODE,TID1,ISTAT)
      IF (ISTAT.NE.0) GO TO 70
      CALL TBIGET(TID1,NCO,NRO,NSC,NAC,NAR,ISTAT)
      IF (ISTAT.NE.0) GO TO 70
C
C ... find output column
C
      CALL TBCSER(TID1,COLRFO(1),OCOL(1),ISTAT)
      IF (OCOL(1).EQ.-1) THEN
          CALL STTPUT('Error in column reference',IS)
          GO TO 70
      END IF
      IF (BIN) THEN
          CALL TBCSER(TID1,COLRFO(3),OCOL(3),ISTAT)
          IF (OCOL(3).EQ.-1) THEN
              CALL STTPUT('Error in column reference',IS)
              GO TO 70
          END IF
      END IF
C
C ... copy selected values of the table into work space
C
      NBY1 = 8*NRO
      CALL TDMGET(NBY1,JP1,IS)
      CALL TDMGET(NBY1,JP2,IS)
      CALL TDMGET(NBY1,JP3,IS)
      IF (BIN) THEN
          CALL COPYIT(TID1,OCOL(1),OCOL(3),NRO,
     +                MADRID(JP1),MADRID(JP3),NO)
          IF (NO.LE.3) THEN
              CALL STTPUT(' Not enough points in table ',ISTAT)
              GO TO 60
          END IF
      ELSE
          CALL BCPYIT(TID1,OCOL(1),NRO,MADRID(JP1),
     +                 MADRID(JP3),NO)
          IF (NO.LE.3) THEN
              CALL STTPUT(' Not enough points in table ',ISTAT)
              GO TO 60
          END IF
      END IF
C
C  ...  Now, at last, we call the interpolating (and integrating)
C       routines:
C  !!!! Watch out for VECI,VECD arrays in INTERP,
C       dimension => NINT !!!!!!
C
      CALL REBMET(NI,MADRID(IP1),MADRID(IP2),MADRID(IP3),NO,
     +             MADRID(JP1),MADRID(JP3),NFUNC,NFPAR,NFPACT,DPARM,
     +             INOP,NINT,MADRID(JP2),RMIN,RMAX)
C
C ... end
C
      CALL TBTCLO(TID,ISTAT)
C
C ... output table
C
C ... check if column already exists
C
      CALL TBCSER(TID1,COLRFO(2),IC,ISTAT)
      IF (IC.EQ.-1) THEN
          LABEL  = COLRFO(2) (2:)
          CUNIT  = '  '
          CALL TBCINI(TID1,D_R4_FORMAT,1,'E16.8',CUNIT,LABEL,IC,ISTAT)
      END IF
      NMAX   = NO
      CALL OCOPY(TID1,NMAX,IC,MADRID(JP2))
C
C ... free memory
C
      CALL TDMFRE(NBY1,JP1,IS)
      CALL TDMFRE(NBY1,JP2,IS)
      CALL TDMFRE(NBY1,JP3,IS)
   60 CONTINUE
      CALL TDMFRE(NBY,IP1,IS)
      CALL TDMFRE(NBY,IP2,IS)
      CALL TDMFRE(NBY,IP3,IS)
C
C ... end
C
      CALL TBIPUT(TID1,0,NO,ISTAT)
      CALL TBSINI(TID1,ISTAT)
      CALL DSCUPT(TID1,TID1,' ',ISTAT)
      CALL TBTCLO(TID1,ISTAT)
   70 IF (ISTAT.NE.0) THEN
          WRITE (MSG(13:16),9000) ISTAT
          CALL TDERRR(ISTAT,MSG,JSTAT)
      END IF
   80 RETURN
 9000 FORMAT (I4)
      END

      SUBROUTINE COPYIT(TID,IX,IW,NROW,X1,W1,N1)
C
      IMPLICIT NONE
      INTEGER  TID, NROW, N1, IX, IW
      DOUBLE PRECISION X1(NROW),W1(NROW)
C
      DOUBLE PRECISION XVAL, YVAL
      INTEGER I, ISTAT
      LOGICAL ISEL, INULL1, INULL2
C
      N1     = 0
      DO 10 I = 1,NROW
          CALL TBSGET(TID,I,ISEL,ISTAT)
          CALL TBERDD(TID,I,IX,XVAL,INULL1,ISTAT)
          CALL TBERDD(TID,I,IW,YVAL,INULL2,ISTAT)
          IF (ISEL .AND.(.NOT.INULL1).AND.(.NOT.INULL2)) THEN
              N1     = N1 + 1
              X1(N1) = XVAL
              W1(N1) = YVAL
          END IF
   10 CONTINUE
      RETURN

      END

      SUBROUTINE BCPYIT(TID,IX,NROW,X1,W1,N1)
C
      IMPLICIT NONE
      INTEGER  TID, IX, NROW, N1
      DOUBLE PRECISION X1(NROW),W1(NROW),XD,XM
C
      DOUBLE PRECISION XVAL
      INTEGER I, ISTAT
      LOGICAL ISEL, INULL
C
      N1     = 0
      DO 10 I = 1,NROW
          CALL TBSGET(TID,I,ISEL,ISTAT)
          CALL TBERDD(TID,I,IX,XVAL,INULL,ISTAT)
          IF (ISEL .AND.(.NOT.INULL)) THEN
              N1     = N1 + 1
              X1(N1) = XVAL
          END IF
   10 CONTINUE
      IF (X1(1).LT.X1(N1)) THEN
          XM     = 1.E+33
          DO 20 I = 2,N1
              XD     = X1(I) - X1(I-1)
              IF (XD.LT.XM) XM     = XD
   20     CONTINUE
      ELSE
          XM     = -1.E+33
          DO 30 I = 2,N1
              XD     = X1(I) - X1(I-1)
              IF (XD.GT.XM) XM     = XD
   30     CONTINUE
      END IF
      DO 40 I = 1,N1
          W1(I)  = XM
   40 CONTINUE
      RETURN
      END

      SUBROUTINE OCOPY(TID1,N,ICOL,ARRAY)
C
C ...  access output table column values from address
C
      IMPLICIT NONE
      INTEGER TID1, N, ICOL, L, ISTAT
      REAL ARRAY(N)
C
      DO 10 L = 1,N
          CALL TBEWRR(TID1,L,ICOL,ARRAY(L),ISTAT)
   10 CONTINUE
      RETURN
      END

      SUBROUTINE IMVAL2(NP,XS,WS,X,W)
C
C  Fill arrays X,W for image world coordinates and pixel size data
C
      IMPLICIT NONE
      INTEGER  I, NP
      REAL     XS, WS
      DOUBLE PRECISION X(NP),W(NP),WSS,XSS
C
      WSS    = WS
      XSS    = XS
      DO 10 I = 1,NP
          X(I)   = XSS + WSS* (I-1)
          W(I)   = WSS
   10 CONTINUE
      RETURN
      END

      SUBROUTINE IMVAL3(NP,XS,WS,YY,X,W,Y)
C
C  Fill arrays X,W for image world coordinates and pixel size data
C  Convert YY(real*4) into Y(real*8)
C
      IMPLICIT NONE
      INTEGER  NP, I
      DOUBLE PRECISION X(NP),W(NP),Y(NP),WSS,XSS
      REAL     XS, WS
      REAL YY(NP)
C
      WSS    = WS
      XSS    = XS
      DO 10 I = 1,NP
          X(I)   = XSS + WSS* (I-1)
          W(I)   = WSS
          Y(I)   = YY(I)
   10 CONTINUE
      RETURN
      END

      SUBROUTINE IMVAL4(NP,XS,WS,YS,X,Y,W)
C ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C copy IMAGE data into real*8
C       HERE: copy indep, dep and  bin
C -------------------------------------------------------------------
C
      IMPLICIT NONE
C
      INTEGER NP, I
      DOUBLE PRECISION X(NP),W(NP),Y(NP),WSS,XSS
      REAL YS(NP)
      REAL     XS, WS
C
      WSS    = WS
      XSS    = XS
      DO 10 I = 1,NP
          X(I)   = XSS + WSS* (I-1)
          W(I)   = WSS
          Y(I)   = YS(I)
   10 CONTINUE
      RETURN
      END

      SUBROUTINE BCOPYC(TID,IX,IY,IW,NROW,X1,Y1,W1,N1)
C
C  REAL*8, BIN CREATION
C
C copy relevant part of the table
C       HERE: copy indep, dep and  bin --> real*8
C             then create bin as the minimum of spacing
C
      IMPLICIT NONE
      INTEGER  TID, IX, IY, IW, NROW, N1
      DOUBLE PRECISION X1(NROW),Y1(NROW),W1(NROW)
C
      INTEGER I, ISTAT
      LOGICAL ISEL, INULL1, INULL2, INULL3
      DOUBLE PRECISION X, Y, W
C
C
      N1     = 0
      DO 10 I = 1,NROW
          CALL TBSGET(TID, I, ISEL, ISTAT)          
          CALL TBERDD(TID, I, IX, X, INULL1, ISTAT)          
          CALL TBERDD(TID, I, IY, Y, INULL2, ISTAT)          
          CALL TBERDD(TID, I, IW, W, INULL3, ISTAT)          
          IF (ISEL.AND.(.NOT.INULL1).AND.(.NOT.INULL2)
     .            .AND.(.NOT.INULL3)) THEN
              N1     = N1 + 1
              X1(N1) = X
              Y1(N1) = Y
              W1(N1) = W
          END IF
   10 CONTINUE
      RETURN
      END

      SUBROUTINE BCOPYI(TID,IX,IY,NROW,X1,Y1,W1,N1)
C
C  REAL*8, BIN CREATION
C
C copy relevant part of the table
C       HERE: copy indep, dep and  bin --> real*8
C             then create bin as the minimum of spacing
      IMPLICIT NONE
      INTEGER  TID, IX, IY, NROW, N1
      DOUBLE PRECISION X1(NROW),Y1(NROW),W1(NROW)
C
      LOGICAL  INULL1, INULL2, ISEL
      DOUBLE PRECISION X, Y, XM, XD
      INTEGER  I, ISTAT
C
C
      N1     = 0
      DO 10 I = 1,NROW
          CALL TBSGET(TID, I, ISEL, ISTAT)          
          CALL TBERDD(TID, I, IX, X, INULL1, ISTAT)          
          CALL TBERDD(TID, I, IY, Y, INULL2, ISTAT)          
          IF (ISEL.AND.(.NOT.INULL1).AND.(.NOT.INULL2)) THEN
              N1     = N1 + 1
              X1(N1) = X
              Y1(N1) = Y
          END IF
   10 CONTINUE
      IF (X1(1).LT.X1(N1)) THEN
          XM     = 1.D+33
          DO 20 I = 2,N1
              XD     = X1(I) - X1(I-1)
              IF (XD.LT.XM) XM     = XD
   20     CONTINUE
      ELSE
          XM     = -1.D+33
          DO 30 I = 2,N1
              XD     = X1(I) - X1(I-1)
              IF (XD.GT.XM) XM     = XD
   30     CONTINUE
      END IF
      DO 40 I = 1,N1
          W1(I)  = XM
   40 CONTINUE
      RETURN
      END

      SUBROUTINE IMVALS(NP,XS,WS,X,W)
C
C copy IMAGE data into real*8
C HERE: copy indep, dep and  bin
C
      IMPLICIT NONE
C
      INTEGER  NP, I
      REAL     XS, WS
      DOUBLE PRECISION X(NP),W(NP),WSS,XSS
C
      WSS    = DBLE(WS)
      XSS    = DBLE(XS)
      DO 10 I = 1,NP
          X(I)   = XSS + WSS* (I-1)
          W(I)   = WSS
   10 CONTINUE
      RETURN
      END