C  @(#)tdcspl.for	19.1 (ESO-DMD) 02/25/03 14:11:16 
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  14:27 - 19 NOV 1987
C
C.LANGUAGE: F77+ESOext
C
C.AUTHOR: J.D.PONZ
C
C.IDENTIFICATION
C
C  program TDCSPL
C
C.PURPOSE
C
C  TABLE TO IMAGE CONVERSION
C  Execute the command
C  CONV/TAB image = table colx [z] refima SPLINE deg
C
C.KEYWORDS
C
C  table, image, conversion
C
C
C.ALGORITHM
C
C  spline interpolation
C
C  Ref.: P.Diercxx, 1982, Computer Graphics and Image Processing
C                         vol. 20, 171 - 184.
C
C.MODIFS
C     M Peron 071190 , in order to avoid a shift in the output frame
C  
C.VERSION
C 020807	last modif
C 
C-----------------------------------------------------------
C
      SUBROUTINE TDCSPL
      IMPLICIT NONE
C
      INTEGER     MADRID
      INTEGER PNVALS, IMNO
      INTEGER ICOL(3),NPTOT(100),NPIX(2)
      INTEGER NORDER(100)
      INTEGER NRDATA(1000)
      INTEGER STATUS, N1, IDEGR, NPAR, ISTAT, NN, NAXIS, NO, NB
      INTEGER NBY, TID, IS, NCOL, NROW, NSC, NAC, NAR, I, IREF
      INTEGER I1, I2, NDIM, KUN, KNUL
C 
      INTEGER*8  IP1, IP2, IP3, PNTR
C
      REAL S,SS,DELTA,YMAX,YMIN
      REAL STEP(2),START(2),RMIN,RMAX
      REAL CUTS(4),WSTART(100)
      REAL T(1000),C(1000)
C
      DOUBLE PRECISION DSTART(2), DSTEP(2)
C
      CHARACTER*16 MSG
      CHARACTER*80 TABLE,IMAGE,REFIMA
      CHARACTER*8 TYPE
      CHARACTER*17 COLREF(2)
      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 PNVALS/8/
      DATA MSG/'ERR:TCONVTBLxxxx'/
      DATA NPIX /1,1/
      DATA DSTEP /1.,1./
      DATA DSTART /0.,0./
C
C ... get input parameters + default
C
      CALL TDPGET(PNVALS,NPAR,ISTAT)
      IF (ISTAT.NE.0) GO TO 30
      IMAGE     = TPARBF(1)
      TABLE     = TPARBF(3)
      COLREF(1) = TPARBF(4)
      COLREF(2) = TPARBF(5)
      REFIMA = TPARBF(6)
      CALL STKRDR('INPUTR',1,1,NN,SS,KUN,KNUL,ISTAT)
      CALL STKRDR('INPUTR',2,1,NN,S,KUN,KNUL,ISTAT)
C      IDEGR = SS
C .................................... change parameters
      IDEGR = S
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)
      START(1) = DSTART(1)
      START(2) = DSTART(2)
      CALL STDRDD(IREF,'STEP',1,NAXIS,NN,DSTEP,KUN,KNUL,ISTAT)
      STEP(1) = DSTEP(1)
      STEP(2) = DSTEP(2)
      IF (NAXIS.GT.1 .AND. NPIX(2).GT.1) THEN
          TYPE   = ' '
          CALL STDFND(IREF,'WSTART',TYPE,NO,NB,ISTAT)
          IF (TYPE(1:1).EQ.' ') THEN
              CALL STTPUT(' Wrong reference image ',ISTAT)
              GO TO 30
          END IF
          CALL STDRDR
     .         (IREF,'WSTART',1,NPIX(2),NN,WSTART,KUN,KNUL,ISTAT)
          CALL STDRDI
     .         (IREF,'NPTOT',1,NPIX(2),NN,NPTOT,KUN,KNUL,ISTAT)
          CALL STDRDI
     .         (IREF,'NORDER',1,NPIX(2),NN,NORDER,KUN,KNUL,ISTAT)
          NDIM   = 2

      ELSE
          NPIX(2) = 1
          NDIM   = 1
      END IF
C
C ... init input table
C
      TID = -1
      CALL TBTOPN(TABLE,F_I_MODE,TID,STATUS)
      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
      DO 10, I = 1, 2
          CALL TBCSER(TID,COLREF(I),ICOL(I),ISTAT)
          IF (ICOL(I).EQ.-1) THEN
              ISTAT  = ERRCOL
              GO TO 30
          END IF
          CALL TBUGET(TID,ICOL(I),WORK,ISTAT)
          I1     = (I-1)*16 + 1
          I2     = I1 + 15
          CUNIT(I1:I2) = WORK(1:16)
   10 CONTINUE
C
C ... copy selected values of the table into work space
C
      NBY = 4*NROW
      CALL TDMGET(NBY,IP1,IS)
      CALL TDMGET(NBY,IP2,IS)
      CALL TDMGET(NBY,IP3,IS)
      CALL ACOPY(TID,ICOL(1),ICOL(2),NROW,
     +          MADRID(IP1),MADRID(IP2),MADRID(IP3),YMIN,YMAX,N1)
      IF (N1.LT.IDEGR+1) THEN
          CALL STTPUT(' Not enough points in table ',ISTAT)
          GO TO 20

      END IF
C
C ... map image
C
      CUNIT(1:16) = 'FLUX'
      I2 = INDEX(TABLE,' ')
      IF (I2.lt.1) I2 = LEN(TABLE)
      IDENT  = 'TABLE: '//TABLE(1:I2)//'COLS. :'//COLREF(1)//COLREF(2)
      CALL STIPUT(IMAGE,10,1,1,NDIM,NPIX,DSTART,DSTEP,IDENT,CUNIT,
     +            PNTR,IMNO,STATUS)
C
C ... interpolate from table values
C
      DELTA  = YMAX - YMIN
      SS     = DELTA* (DELTA*0.5+YMIN)*SS*0.01
      CALL AINTER(N1,MADRID(IP1),MADRID(IP2),NPIX(1),NPIX(2),
     +            MADRID(PNTR),START,STEP,WSTART,NPTOT,RMIN,RMAX,
     +            MADRID(IP3),IDEGR,SS,T,C,NRDATA)
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)
      IF (NPIX(2).GT.1) THEN
          CALL STDWRR(IMNO,'WSTART',WSTART,1,NPIX(2),KUN,ISTAT)
          CALL STDWRI(IMNO,'NPTOT',NPTOT,1,NPIX(2),KUN,ISTAT)
          CALL STDWRI(IMNO,'NORDER',NORDER,1,NPIX(2),KUN,ISTAT)
      END IF
C
C ... free memory
C
   20 CALL TDMFRE(NBY,IP1,IS)
      CALL TDMFRE(NBY,IP2,IS)
      CALL TDMFRE(NBY,IP3,IS)
C
C ... end
C
      CALL DSCUPT(IMNO,IMNO,' ',ISTAT)
      CALL TBTCLO(TID,ISTAT)
   30 IF (ISTAT.NE.0) THEN
          WRITE (MSG(13:16),9000) ISTAT
          CALL TDERRR(ISTAT,MSG,STATUS)
      END IF
      RETURN
 9000 FORMAT (I4)
      END

      SUBROUTINE ACOPY(TID,ICOL1,ICOL2,NROW,X1,Y1,W,YMIN,YMAX,N1)
C
C copy relevant part of the table
C
      IMPLICIT NONE
C
      INTEGER TID, ICOL1, ICOL2
      INTEGER NROW, N1, I
      INTEGER STATUS
C
      REAL YMIN, YMAX
      REAL X1(NROW),Y1(NROW), X, Y
      REAL W(NROW)
C 
      LOGICAL ISEL, INULL1, INULL2
C
C
      N1     = 0
      DO 10, I = 1,NROW
          CALL TBSGET(TID,I,ISEL,STATUS)
          CALL TBERDR(TID,I,ICOL1,X,INULL1,STATUS)
          CALL TBERDR(TID,I,ICOL2,Y,INULL2,STATUS)
          IF (ISEL.AND. .NOT.INULL1 .AND. .NOT.INULL2) THEN
              N1     = N1 + 1
              X1(N1) = X
              Y1(N1) = Y
              W(N1)  = 1.
              IF (N1.EQ.1) THEN
                  YMIN   = Y
                  YMAX   = Y
              ELSE
                  IF (Y.GT.YMAX) YMAX   = Y
                  IF (Y.LT.YMIN) YMIN   = Y
              END IF
          END IF
   10 CONTINUE
      RETURN
      END

      SUBROUTINE AINTER(NROW,X,Y,NPIX1,NPIX2,F,START,STEP,WSTART,NPTOT,
     +                  RMIN,RMAX,W,IDEGR,S,T,C,NRDATA)
C
C call for each order to the interpolation routine
C
      IMPLICIT NONE
C 
      INTEGER MADRID
      INTEGER NROW,NPIX1, NPIX2,IDEGR
      INTEGER NPTOT(NPIX2)
      INTEGER ISTAT, NRDATA(1),NBYTES
      INTEGER IORD, NP, I, NK1, LD, NK, NPLUS,IERR
C 
      INTEGER*8   IPQ
C 
      REAL X(NROW),WSTART(NPIX2),F(NPIX1,NPIX2)
      REAL Y(NROW),START(2),STEP(2),RMIN,RMAX,S
      REAL W(1),T(1),C(1)
      REAL    DERIV, STARTX, STEPX, ENDX, VAL, XP, FP, FP0
      REAL    FPOLD, AMIN1, AMAX1
C 
      INCLUDE 'MID_INCLUDE:TABLES.INC'
      COMMON /VMR/MADRID(1)
      INCLUDE 'MID_INCLUDE:TABLED.INC'
C
      RMIN   = 0.
      RMAX   = 0.
      STARTX = START(1)
      STEPX  = STEP(1)
      NP     = NPIX1
      NBYTES = NROW*24  
C     NBYTES = NROW*6   M.Peron 23/07/93
      CALL TDMGET(NBYTES,IPQ,ISTAT)
C
C ... order loop
C
      DO 50, IORD = 1,NPIX2
          IF (NPIX2.GT.1) THEN
C              STARTX = WSTART(IORD)
              NP     = NPTOT(IORD)
          END IF
C
C ... loop for each point  W(1), IDEGR, S, T(1), C(1), NRDATA(1)
C
C          ENDX   = STARTX + (NP-1)*STEPX
           ENDX   = X(NROW)
           STARTX = X(1)
          CALL SMOOT(X,Y,W,MADRID(IPQ),NROW,STARTX,ENDX,IDEGR,S,NK,T,
     +               C,FP,0,IERR,NRDATA,FP0,FPOLD,NPLUS)
          IF (IERR.EQ.-2)
     .    CALL STTPUT('Polynomial approximation',ISTAT)
          IF (IERR.EQ.1) 
     .    CALL STTPUT('Local storage exceeded (S too small)',ISTAT)
          IF (IERR.EQ.2) 
     .    CALL STTPUT('Tolerance parameter too small',ISTAT)
          IF (IERR.EQ.3)
     .    CALL STTPUT('Maximum number of iter. exceeded',ISTAT)
          IF (IERR.EQ.10) 
     .    CALL STTPUT('Invalid input arguments',ISTAT)
          IF (IERR.GT.0) RETURN
          LD     = IDEGR + 1
          NK1    = NK - IDEGR - 1
          IF (NPIX2.EQ.1) STARTX = START(1) 
          DO 30, I = 1,NP
              XP     = STARTX + (I-1)*STEPX
   10         IF ((XP.LT.T(LD+1)) .OR. (LD.EQ.NK1)) GO TO 20
              LD     = LD + 1
              GO TO 10
   20         VAL    = DERIV(T,NK,C,NK1,0,XP,LD)
              RMIN   = AMIN1(RMIN,VAL)
              RMAX   = AMAX1(RMAX,VAL)
              F(I,IORD) = VAL
   30     CONTINUE
          DO 40, I = NP + 1,NPIX1
              F(I,IORD) = 0.
   40     CONTINUE
   50 CONTINUE
      CALL TDMFRE(NBYTES,IPQ,ISTAT)
      RETURN

      END
