C @(#)tdinter.for	17.1.1.1 (ESO-DMD) 01/25/02 17:47:14 
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  9:41 - 9 DEC 1987
C
C.LANGUAGE: F77+ESOext
C
C.AUTHOR: J.D.PONZ
C
C.KEYWORDS
C
C  TABLE, IMAGES, SPLINE, INTERPOLATION
C
C.PURPOSE
C
C  spline interpolation image to image in 1 dimension
C
C.ALGORITHM
C
C  Interpolate input data YI sampled at XI to
C  output data YO sampled at XO with binwidth.
C
C  Restriction:  Input and output independent variables are
C                assumed to be monotonic increasing or decreasing
C
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
C  WARNINGS : module to be tested. Replace allocated arrays by alloc memory
C
C-----------------------------------------------------------
C
      SUBROUTINE TDINII
      IMPLICIT   NONE
      INTEGER    NFPAR
      PARAMETER (NFPAR = 3)
C
      INTEGER MADRID
      INTEGER*8  IP1,IP2,IPQ,IPNTR,JPNTR
      INTEGER IACT,IDEGRE,IDERIV,IERR
      INTEGER IS,ISTAT,NAXISI,NAXISO
      INTEGER NDIMI,NDIMO,NFPACT,NI,NK,NN,NO
      INTEGER NPLUS, NBY, KUN, KNUL
      INTEGER NRDATA(1000),IMNO, IREF, IMOUT
      INTEGER NPIXI(3),NPIXO(3)
C
      REAL DELTA,FP,FP0,FPOLD,S,X1
      REAL RMIN,RMAX
      REAL CUTS(4)
      REAL FPARM(NFPAR),X0,WX
      REAL T(1000),C(1000)
C
      DOUBLE PRECISION STARTI(3), STEPI(3)
      DOUBLE PRECISION STARTO(3), STEPO(3)
C
      CHARACTER*60 INPIMA,OUTIMA,REFIMA
      CHARACTER*72 IDENTO,IDENTI
      CHARACTER*80 CUNITO,CUNITI
      CHARACTER    KBUF*40
C
      INCLUDE 'MID_INCLUDE:TABLES.INC'
      COMMON /VMR/MADRID(1)
      INCLUDE 'MID_INCLUDE:TABLED.INC'
      DATA INPIMA/'   '/
      DATA OUTIMA/'   '/
      DATA REFIMA/'   '/
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 STKRDR('INPUTR',1,NFPAR,NFPACT,FPARM,KUN,KNUL,ISTAT)
      IDEGRE = FPARM(2)
      IDERIV = FPARM(3)
      IF (IDEGRE.LT.1 .OR. IDEGRE.GT.5) THEN
          CALL STTPUT(' Degree parameter out of range',ISTAT)
          GO TO 10
      END IF
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,STARTO,KUN,KNUL,ISTAT)
      CALL STDRDD(IREF,'STEP',1,NAXISO,NN,STEPO,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 20
      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,IMNO,ISTAT)
      IF (NAXISI.GT.1 .AND. NPIXI(2).GT.1) THEN
          CALL STTPUT(' Only 1-D images are considered.',ISTAT)
          GO TO 10
      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,STARTO,STEPO,IDENTO,CUNITO,JPNTR,IMOUT,ISTAT)
C
C ... fill the arrays XI,WI
C
      NI     = NPIXI(1)  !   XI
      NBY    = 4*NI
      CALL TDMGET(NBY,IP1,IS)
      CALL TDMGET(NBY,IP2,IS)
      CALL TDMGET(NBY*6,IPQ,IS)
      X0     = STARTI(1)
      WX     = STEPI(1)
      X1     = X0 + (NI-1)*WX
      CALL IMVAL1(NI,X0,WX,MADRID(IP1),MADRID(IP2),MADRID(IPNTR),RMIN,
     +             RMAX)
C
C ... setup parameters
C
      DELTA  = RMAX - RMIN
      S      = DELTA* (DELTA*0.5+RMIN)*FPARM(1)*0.01
C S       = FPARM(1)
      S = ABS(S)
      IDEGRE = FPARM(2)
      IDERIV = FPARM(3)
C
C ... compute spline coeffs
C
      WRITE (KBUF,12345) S
      CALL STTPUT(KBUF,ISTAT)
      CALL SMOOT(MADRID(IP1),MADRID(IPNTR),MADRID(IP2),MADRID(IPQ),
     +           NI,X0,X1,IDEGRE,
     +           S,NK,T,C,FP,0,IERR,NRDATA,FP0,FPOLD,NPLUS)
      IF (IERR.EQ.-2) CALL STTPUT('Polynomial approximation',ISTAT)
      IF (IERR.EQ.1) THEN
          CALL STTPUT('Local storage exceeded (S too small)',ISTAT)
          GO TO 10
      END IF
      IF (IERR.EQ.2) THEN
          CALL STTPUT('Tolerance parameter too small',ISTAT)
          GO TO 10
      END IF
      IF (IERR.EQ.3) THEN
          CALL STTPUT('Maximum number of iter. exceeded',ISTAT)
          GO TO 10
      END IF
      IF (IERR.EQ.10) THEN
          CALL STTPUT('Invalid input arguments',ISTAT)
          GO TO 10
      END IF
      NO     = NPIXO(1)
      X0     = STARTO(1)
      WX     = STEPO(1)
C
C ... compute interpolated values
C
      CALL COMPIM(NO,X0,WX,MADRID(JPNTR),NK,IDEGRE,IDERIV,T,C,RMIN,RMAX)
C
C ... write cuts
C
      CUTS(1) = RMIN
      CUTS(2) = RMAX
      CUTS(3) = RMIN
      CUTS(4) = RMAX
      CALL STDWRR(IMOUT,'LHCUTS',CUTS,1,4,KUN,ISTAT)
      CALL DSCUPT(IMOUT,IMOUT,' ',ISTAT)
C
C ... free memory
C
      CALL TDMFRE(NBY,IP1,IS)
      CALL TDMFRE(NBY,IP2,IS)
      CALL TDMFRE(NBY*6,IPQ,IS)
C
C ... end
C
      CALL STFCLO(IMOUT,ISTAT)
   10 CALL STFCLO(IMNO, ISTAT)
   20 RETURN
C 
12345 FORMAT('S = ',F12.6)
      END

      SUBROUTINE TDINIT
      IMPLICIT NONE
      INTEGER     NFPAR
      PARAMETER  (NFPAR = 3)
C
      INTEGER MADRID
      INTEGER IACT,IC,IDEGRE,IDERIV,IERR,II,NBY, TID
      INTEGER*8   IP1,IP2,IPQ,IPNTR
      INTEGER IS,ISTAT,JSTAT,NAXIS
      INTEGER NCOL,NDIM,NI,NK,NPLUS,NROW,NSC,NAC,NAR
      INTEGER INDEX,IMNO
      INTEGER NFPACT,KUN,KNUL
      INTEGER ICOL(3),NPIX(3)
      INTEGER NRDATA(1000)
C
      REAL DELTA,FP,FP0,FPOLD,S,X1
      REAL RMIN,RMAX
      REAL FPARM(NFPAR),X0,WX
      REAL T(1000),C(1000)
C
      DOUBLE PRECISION START(3), STEP(3)
C
      CHARACTER*60 TABLE,IMAGE
      CHARACTER*16 MSG
      CHARACTER*17 COLREF(3),LABEL
      CHARACTER*53 OCOLS
      CHARACTER*72 IDENT
      CHARACTER*80 CUNIT
      CHARACTER   KBUF*40
C
      INCLUDE 'MID_INCLUDE:TABLES.INC'
      COMMON /VMR/MADRID(1)
      INCLUDE 'MID_INCLUDE:TABLED.INC'
C
      DATA MSG/'ERR:REBISTBxxx'/
      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 STKRDR('INPUTR',1,NFPAR,NFPACT,FPARM,KUN,KNUL,ISTAT)
      IDEGRE = FPARM(2)
      IDERIV = FPARM(3)
      IF (IDEGRE.LT.1 .OR. IDEGRE.GT.5) THEN
          CALL STTPUT('Parameter degree out of range',ISTAT)
          GO TO 20
      END IF
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'
          RETURN
      ELSE
          COLREF(1) = OCOLS(1:II-1)
          COLREF(2) = OCOLS(II+1:)
      END IF
C
C ... access input/output table
C
      CALL TBTOPN(TABLE,F_U_MODE,TID,ISTAT)
      IF (ISTAT.NE.0) GO TO 20
      CALL TBIGET(TID,NCOL,NROW,NSC,NAC,NAR,ISTAT)
      IF (ISTAT.NE.0) GO TO 20
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 20
      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 20
      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)
      NBY    = 4*NI
      CALL TDMGET(NBY,IP1,IS)
      CALL TDMGET(NBY,IP2,IS)
      CALL TDMGET(NBY*6,IPQ,IS)
      X0     = START(1)
      WX     = STEP(1)
      X1     = X0 + (NI-1)*WX
      CALL IMVAL1(NI,X0,WX,MADRID(IP1),MADRID(IP2),MADRID(IPNTR),RMIN,
     +             RMAX)
C
C setup parameters
C
      DELTA  = RMAX - RMIN
      S      = DELTA*(DELTA*0.5+RMIN)*FPARM(1)*0.01
      S = ABS(S)
C
C ... compute spline coeffs
C
      WRITE (KBUF,12345) S
      CALL STTPUT(KBUF,ISTAT)
      CALL SMOOT(MADRID(IP1),MADRID(IPNTR),MADRID(IP2),MADRID(IPQ),
     +           NI,X0,X1,IDEGRE,
     +           S,NK,T,C,FP,0,IERR,NRDATA,FP0,FPOLD,NPLUS)
      IF (IERR.EQ.-2) CALL STTPUT('Polynomial approximation',ISTAT)
      IF (IERR.EQ.1) THEN
          CALL STTPUT('Local storage exceeded (S too small)',ISTAT)
          GO TO 20
      END IF
      IF (IERR.EQ.2) THEN
          CALL STTPUT('Tolerance parameter too small',ISTAT)
          GO TO 20
      END IF
      IF (IERR.EQ.3) THEN
          CALL STTPUT('Maximum number of iter. exceeded',ISTAT)
          GO TO 20
      END IF
      IF (IERR.EQ.10) THEN
          CALL STTPUT('Invalid input arguments',ISTAT)
          GO TO 20
      END IF
C
C ... find if column already exists
C
      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
C
C ... compute interpolated values
C
      CALL COMPTA(TID,ICOL(1),IC,NROW,NK,IDEGRE,IDERIV,T,C)
C
      CALL DSCUPT(TID,TID,' ',ISTAT)
      CALL TDMFRE(NBY,IP1,IS)
      CALL TDMFRE(NBY,IP2,IS)
      CALL TDMFRE(NBY*6,IPQ,IS)
C
C ... end
C
      CALL TBTCLO(TID,ISTAT)
   20 IF (ISTAT.NE.0) THEN
          WRITE (MSG(13:16),9000) ISTAT
          CALL TDERRR(ISTAT,MSG,JSTAT)
      END IF
      RETURN
C 
 9000 FORMAT (I4)
12345 FORMAT('S = ',F12.6)
      END

      SUBROUTINE TDINTI
      IMPLICIT NONE
      INTEGER   NFPAR
      PARAMETER (NFPAR = 3)
C
      INTEGER MADRID
      INTEGER I,I1,I2,IACT,IB,IDEGRE,IDERIV, NBY
      INTEGER*8    IP1,IP2,IP3,IPQ,IPNTR
      INTEGER IERR,II,IS,ISTAT
      INTEGER JSTAT,NAXIS,NCOL,NDIM,NI,NK
      INTEGER NN,NO,NPLUS,NROW,NSC,NAC,NAR,TID
      INTEGER INDEX,IMOUT
      INTEGER NFPACT,KUN,KNUL,IREF
      INTEGER ICOL(3),NPIX(3)
      INTEGER NRDATA(1000)
      REAL DELTA,FP,FPO,FPOLD,S,X1
      REAL RMIN,RMAX
      REAL CUTS(4)
      REAL FPARM(NFPAR),X0,WX
      REAL T(1000),C(1000)
      DOUBLE PRECISION START(3), STEP(3)
      CHARACTER*60 TABLE,IMAGE,REFIMA
      CHARACTER*16 MSG
      CHARACTER*17 COLREF(3)
      CHARACTER*53 INCOLS
      CHARACTER*60 WORK
      CHARACTER*72 IDENT
      CHARACTER*80 CUNIT
      CHARACTER    KBUF*40
C
      INCLUDE  'MID_INCLUDE:TABLES.INC'
      COMMON /VMR/MADRID(1)
      INCLUDE  'MID_INCLUDE:TABLED.INC'
C
      DATA MSG/'ERR:REBISTBxxx'/
      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('COLI',1,1,53,IACT,INCOLS,KUN,KNUL,ISTAT)
      CALL STKRDC('IN_B',1,1,60,IACT,REFIMA,KUN,KNUL,ISTAT)
      CALL STKRDR('INPUTR',1,NFPAR,NFPACT,FPARM,KUN,KNUL,ISTAT)
      IDEGRE = FPARM(2)
      IDERIV = FPARM(3)
      IF (IDEGRE.LT.1 .OR. IDEGRE.GT.5) THEN
          CALL STTPUT('Parameter degree out of range',ISTAT)
          GO TO 30
      END IF
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,START,KUN,KNUL,ISTAT)
      CALL STDRDD(IREF,'STEP',1,NAXIS,NN,STEP,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 30
      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'
          RETURN
      ELSE
          COLREF(1) = INCOLS(1:II-1)
          COLREF(2) = INCOLS(II+1:)
      END IF
C
C ... init input table
C
      CALL TBTOPN(TABLE,F_I_MODE,TID,ISTAT)
      IF (ISTAT.NE.0) GO TO 30
      CALL TBIGET(TID,NCOL,NROW,NSC,NAC,NAR,ISTAT)
      IF (ISTAT.NE.0) GO TO 30
C
C ... find input column
C
      IB     = 2
      DO 10 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 30
          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
   10 CONTINUE
C
C ... copy selected values of the table into work space
C  !   indep
      NBY    = 4*NROW
      CALL TDMGET(NBY,IP1,IS)
      CALL TDMGET(NBY,IP2,IS)
      CALL TDMGET(NBY,IP3,IS)
      CALL TDMGET(NBY*6,IPQ,IS)
      CALL BCOPYT(TID,ICOL(1),ICOL(2),NROW,
     +            MADRID(IP1),MADRID(IP2),MADRID(IP3),
     +            NI,RMIN,RMAX,X0,X1)
      IF (NI.LE.3) THEN
          CALL STTPUT(' Not enough points in table ',ISTAT)
          GO TO 20
      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,START,STEP,IDENT,CUNIT,IPNTR,IMOUT,ISTAT)
C
C setup parameters
C
      DELTA  = RMAX - RMIN
      S      = DELTA* (DELTA*0.5+RMIN)*FPARM(1)*0.01
      S      = ABS(S)
C
C compute spline coeffs.
C
      WRITE (KBUF,12345) S
      CALL STTPUT(KBUF,ISTAT)
      CALL SMOOT(MADRID(IP1),MADRID(IP2),MADRID(IP3),MADRID(IPQ),
     +           NI,X0,X1,IDEGRE,S,
     +           NK,T,C,FP,0,IERR,NRDATA,FPO,FPOLD,NPLUS)
      IF (IERR.EQ.-2) CALL STTPUT('Polynomial approximation',ISTAT)
      IF (IERR.EQ.1) THEN
          CALL STTPUT('Local storage exceeded (S too small)',ISTAT)
          GO TO 30
      END IF
      IF (IERR.EQ.2) THEN
          CALL STTPUT('Tolerance parameter too small',ISTAT)
          GO TO 30
      END IF
      IF (IERR.EQ.3) THEN
          CALL STTPUT('Maximum number of iter. exceeded',ISTAT)
          GO TO 30
      END IF
      IF (IERR.EQ.10) THEN
          CALL STTPUT('Invalid input arguments',ISTAT)
          GO TO 30
      END IF
C
C image output, then fill the XO and WO arrays OF INTERPOLn
C
      NO     = NPIX(1)
      X0     = START(1)
      WX     = STEP(1)
      X1     = X0 + (NO-1)*WX
      CALL COMPIM(NO,X0,WX,MADRID(IPNTR),NK,IDEGRE,IDERIV,T,C,RMIN,RMAX)
C
C ... write cuts
C
      CUTS(1) = RMIN
      CUTS(2) = RMAX
      CUTS(3) = RMIN
      CUTS(4) = RMAX
      CALL STDWRR(IMOUT,'LHCUTS',CUTS,1,4,KUN,ISTAT)
      CALL DSCUPT(IMOUT,IMOUT,' ',ISTAT)
C
C ... free memory
C
C
   20 CONTINUE
      CALL TDMFRE(NBY,IP1,IS)
      CALL TDMFRE(NBY,IP2,IS)
      CALL TDMFRE(NBY,IP3,IS)
      CALL TDMFRE(NBY*6,IPQ,IS) 
C
C ... end
C
      CALL TBTCLO(TID,ISTAT)
   30 IF (ISTAT.NE.0) THEN
          WRITE (MSG(13:16),9000) ISTAT
          CALL TDERRR(ISTAT,MSG,JSTAT)
      END IF
      RETURN
C 
 9000 FORMAT (I4)
12345 FORMAT('S = ',F12.6)
      END

      SUBROUTINE TDINTT
      IMPLICIT NONE
      INTEGER   NFPAR
      PARAMETER (NFPAR = 3)
C
      INTEGER MADRID
      INTEGER I,IACT,IB,IC,IDEGRE,IDERIV,IERR,II
      INTEGER ISTAT,JSTAT,NBY
      INTEGER NCO,NCOL,NI,NK,NPLUS,NRO,NROW,NSC
      INTEGER TID1,TID2,NAC,NAR,INDEX
      INTEGER NFPACT,KUN,KNUL
      INTEGER ICOL(3),OCOL(3)
      INTEGER NRDATA(1000)
C 
      INTEGER*8  IP1,IP2,IP3,IPQ
C 
      REAL DELTA,FP,FP0,FPOLD
      REAL RMAX,RMIN,S,X0,X1
      REAL FPARM(NFPAR)
      REAL T(1000),C(1000)
C 
      CHARACTER*60 TABIN,TABOUT
      CHARACTER*16 MSG
      CHARACTER*17 COLRFI(3),COLRFO(3),LABEL
      CHARACTER*53 INCOLS,OUTCLS
      CHARACTER*80 CUNIT
      CHARACTER    KBUF*40
C
      INCLUDE 'MID_INCLUDE:TABLES.INC'
      COMMON/VMR/MADRID(1)
      INCLUDE 'MID_INCLUDE:TABLED.INC'
C
      DATA MSG/'ERR:REBISTBxxx'/
      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,OUTCLS,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 STKRDR('INPUTR',1,NFPAR,NFPACT,FPARM,KUN,KNUL,ISTAT)
      IDEGRE = FPARM(2)
      IDERIV = FPARM(3)
      IF (IDEGRE.LT.1 .OR. IDEGRE.GT.5) THEN
          CALL STTPUT('Parameter degree out of range',ISTAT)
          GO TO 30
      END IF
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'
          RETURN
      ELSE
          COLRFI(1) = INCOLS(1:II-1)
          COLRFI(2) = INCOLS(II+1:)
      END IF
C
C ... init input table
C
      CALL TBTOPN(TABIN,F_I_MODE,TID1,ISTAT)
      IF (ISTAT.NE.0) GO TO 30
      CALL TBIGET(TID1,NCOL,NROW,NSC,NAC,NAR,ISTAT)
      IF (ISTAT.NE.0) GO TO 30
C
C ... find input column
C
      IB     = 2
      DO 10 I = 1,IB
          CALL TBCSER(TID1,COLRFI(I),ICOL(I),ISTAT)
          IF (ICOL(I).EQ.-1) THEN
              CALL STTPUT('Error in column reference',ISTAT)
              GO TO 30
          END IF
   10 CONTINUE
C
C ... copy selected values of the table into work space
C
      NBY    = 4*NROW
      CALL TDMGET(NBY,IP1,ISTAT)
      CALL TDMGET(NBY,IP2,ISTAT)
      CALL TDMGET(NBY,IP3,ISTAT)
      CALL TDMGET(NBY*6,IPQ,ISTAT)
      CALL BCOPYT(TID1,ICOL(1),ICOL(2),NROW,
     +            MADRID(IP1),MADRID(IP2),MADRID(IP3),
     +            NI,RMIN,RMAX,X0,X1)
      IF (NI.LE.3) THEN
          CALL STTPUT(' Not enough points in table ',ISTAT)
          GO TO 20
      END IF
C
C ... sort out ind., dep. (and bin) variables for output table
C
      II     = INDEX(OUTCLS,',')
      IF (II.EQ.0) THEN
          WRITE (*,*) ' Specification of relevant output',
     +      ' table cols insufficient'
          RETURN
      ELSE
          COLRFO(1) = OUTCLS(1:II-1)
          COLRFO(2) = OUTCLS(II+1:)
      END IF
C
C ... access output table
C
      CALL TBTOPN(TABOUT,F_U_MODE,TID2,ISTAT)
      IF (ISTAT.NE.0) GO TO 30
      CALL TBIGET(TID2,NCO,NRO,NSC,NAC,NAR,ISTAT)
      IF (ISTAT.NE.0) GO TO 30
C
C ... find output column
C
      CALL TBCSER(TID2,COLRFO(1),OCOL(1),ISTAT)
      IF (OCOL(1).EQ.-1) THEN
          CALL STTPUT('Error in column reference',ISTAT)
          GO TO 30
      END IF
C
      CALL TBCSER(TID2,COLRFO(2),IC,ISTAT)
      IF (IC.EQ.-1) THEN
          LABEL  = COLRFO(2) (2:)
          CUNIT  = '  '
          CALL TBCINI(TID2,D_R4_FORMAT,1,'E16.8',
     .                CUNIT,LABEL,IC,ISTAT)
      END IF
C
C ... setup parameters
C
      DELTA  = RMAX - RMIN
      S      = DELTA* (DELTA*0.5+RMIN)*FPARM(1)*0.01
      S      = ABS(S)
C
C ... compute spline coeffs
C
      WRITE (KBUF,12345) S
      CALL STTPUT(KBUF,ISTAT)
      CALL SMOOT(MADRID(IP1),MADRID(IP2),MADRID(IP3),MADRID(IPQ),
     +            NI,X0,X1,IDEGRE,S,
     +           NK,T,C,FP,0,IERR,NRDATA,FP0,FPOLD,NPLUS)
      IF (IERR.EQ.-2) CALL STTPUT('Polynonial approximation',ISTAT)
      IF (IERR.EQ.1) THEN
          CALL STTPUT('Local storage exceeded (S too small)',ISTAT)
          GO TO 30
      END IF
      IF (IERR.EQ.2) THEN
          CALL STTPUT('Tolerance parameter too small',ISTAT)
          GO TO 30
      END IF
      IF (IERR.EQ.3) THEN
          CALL STTPUT('Maximum number of iter. exceeded',ISTAT)
          GO TO 30
      END IF
      IF (IERR.EQ.10) THEN
          CALL STTPUT('Invalid input arguments',ISTAT)
          GO TO 30
      END IF
C
C ... compute interpolated values
C
      CALL COMPTA(TID2,OCOL(1),IC,NRO,NK,IDEGRE,IDERIV,T,C)
C
C ... end
C
      CALL TBTCLO(TID1,ISTAT)
C
C ... free memory
C
C
   20 CONTINUE
      CALL TDMFRE(NBY,IP1,ISTAT)
      CALL TDMFRE(NBY,IP2,ISTAT)
      CALL TDMFRE(NBY,IP3,ISTAT)
      CALL TDMFRE(NBY*6,IPQ,ISTAT)
C
C ... end
C
      CALL DSCUPT(TID2,TID2,' ',ISTAT)
      CALL TBTCLO(TID2,ISTAT)
   30 IF (ISTAT.NE.0) THEN
          WRITE (MSG(13:16),9000) ISTAT
          CALL TDERRR(ISTAT,MSG,JSTAT)
      END IF
      RETURN
C 
 9000 FORMAT (I4)
12345 FORMAT('S = ',F12.6)
      END

      SUBROUTINE IMVAL1(NP,XS,WS,X,W,Y,YMIN,YMAX)
C
C  Subroutine IMVAL1
C
C
C  Fill arrays X,W for image world coordinates and pixel size data
C
      IMPLICIT NONE
      INTEGER  NP, I
      REAL     WS, XS, YMIN, YMAX
      REAL X(NP),W(NP),Y(NP)
C
      YMIN   = Y(1)
      YMAX   = Y(1)
      DO 10 I = 1,NP
          X(I)   = XS + WS* (I-1)
          W(I)   = 1.
          YMIN   = AMIN1(YMIN,Y(I))
          YMAX   = AMAX1(YMAX,Y(I))
   10 CONTINUE
      RETURN
      END

      SUBROUTINE COMPTA(TID,IX,IY,NROW,NK,IDEGRE,IDERIV,T,C)
C
      IMPLICIT NONE
C
      INTEGER   NROW, IDEGRE, IDERIV, IX, IY, NK, TID
      REAL X, Y
      REAL T(1),C(1)
C
      INTEGER I, LD, NK1, STATUS
      REAL    DERIV
      LOGICAL ISEL, INULL
C
      LD     = IDEGRE + 1
      NK1    = NK - IDEGRE - 1
C
      DO 30 I = 1,NROW
         CALL TBSGET(TID, I, ISEL, STATUS)
          IF (ISEL) THEN
              CALL TBERDR(TID, I, IX, X, INULL, STATUS)
              IF (.NOT.INULL) THEN
   10             IF ((X.LT.T(LD+1)) .OR. (LD.EQ.NK1)) GO TO 20
                  LD     = LD + 1
                  GO TO 10
   20             Y   = DERIV(T,NK,C,NK1,IDERIV,X,LD)
                  CALL TBEWRR(TID, I, IY, Y, STATUS)
              ELSE
                  CALL TBEDEL(TID, I, IY, STATUS)
              END IF
          END IF
   30 CONTINUE
      RETURN
      END

      SUBROUTINE COMPIM(NO,X0,WX,Y,NK,IDEGRE,IDERIV,T,C,RMIN,RMAX)
C
C COMPUTE INTERPOLATED VALUES
      IMPLICIT NONE
C
      INTEGER NO, IDEGRE, IDERIV, NK
      REAL Y(NO),T(1),C(1)
C
      INTEGER    LD, I, NK1
      REAL       XP, X0, WX, RMIN, RMAX
      REAL       DERIV
C
      LD     = IDEGRE + 1
      NK1    = NK - IDEGRE - 1
      DO 30 I = 1,NO
          XP     = X0 + (I-1)*WX
   10     IF ((XP.LT.T(LD+1)) .OR. (LD.EQ.NK1)) GO TO 20
          LD     = LD + 1
          GO TO 10
   20     Y(I)   = DERIV(T,NK,C,NK1,IDERIV,XP,LD)
          IF (I.EQ.1) THEN
              RMIN   = Y(1)
              RMAX   = Y(1)
          ELSE
              RMIN   = AMIN1(RMIN,Y(I))
              RMAX   = AMAX1(RMAX,Y(I))
          END IF
   30 CONTINUE
      RETURN
      END
 
      SUBROUTINE BCOPYT(TID,IX,IY,NROW,X1,Y1,W1,N1,RMIN,RMAX,XS,XE)
C
C COPY TABLE VALUES
C
      IMPLICIT NONE
      INTEGER  NROW, IX, IY
      REAL X,Y
      REAL X1(NROW),Y1(NROW),W1(NROW)
      REAL RMIN, RMAX, XS, XE
C
      INTEGER N1, I, STATUS, TID
      LOGICAL ISEL, INULL1, INULL2
C
      N1     = 0
      DO 10 I = 1,NROW
          CALL TBSGET(TID, I, ISEL, STATUS)
          CALL TBERDR(TID, I, IX, X, INULL1, STATUS)
          CALL TBERDR(TID, I, IY, Y, INULL2, STATUS)
          IF (ISEL.AND.(.NOT.INULL1).AND.(.NOT.INULL2)) THEN
              N1     = N1 + 1
              X1(N1) = X
              Y1(N1) = Y
              W1(N1) = 1.
              IF (N1.EQ.1) THEN
                  RMIN   = Y1(1)
                  RMAX   = Y1(1)
                  XS     = X1(1)
              ELSE
                  RMIN   = AMIN1(RMIN,Y1(N1))
                  RMAX   = AMAX1(RMAX,Y1(N1))
                  XE     = X1(N1)
              END IF
          END IF
   10 CONTINUE
      RETURN
      END