C @(#)tdconv.for	19.1 (ESO-IPG) 02/25/03 14:11:11
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:22 - 19 NOV 1987
C
C.LANGUAGE: F77+ESOext
C
C.AUTHOR: J.D.PONZ
C  900219  KB, throw out SX calls...
C
C.IDENTIFICATION
C  program TDCONV
C
C.PURPOSE
C  Table to Image conversion
C  Execute the command
C  CONV/TAB image = table colx [z] refima method
C                 for method = SPLINE                 
C
C.KEYWORDS
C  table, image, conversion
C
C
C.ALGORITHM
C  Use table interface routines for I/O
C   The conversion from table format into image is done
C  by spline interpolation.
C  The interpolation scheme is based on Hermite polynomials
C  Ref .: Publ. of the Dominion Astrophys. Obs. 204, 1982.
C
C 021129		last modif
C 
C-----------------------------------------------------------
C
      SUBROUTINE TDCONV
      IMPLICIT NONE
C
      INTEGER     MADRID
      INTEGER PNVALS, STATUS, N1, NBY, IMNO,REF
      INTEGER ICOL(3),NPTOT(100),NPIX(2),KUN,KNUL
      INTEGER NORDER(100)
      INTEGER NPAR, ISTAT, NN, NAXIS, NO, NB, NDIM, TID, IS
      INTEGER NCOL, NROW, NSC, NAC, NAR, I, I1, I2
      INTEGER*8   PNTR, IP1, IP2
C
      DOUBLE PRECISION STEP(2),START(2)
      REAL   RMIN,RMAX, CUTS(4)
      DOUBLE PRECISION WSTART(100)
C
      CHARACTER*80 TABLE,IMAGE,REFIMA
      CHARACTER*17 COLREF(2)
      CHARACTER*60 WORK
      CHARACTER*16 MSG
      CHARACTER*8  TYPE
      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/, START /0.,0./, STEP /1.,1./
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 STFOPN(REFIMA,D_R4_FORMAT,0,F_IMA_TYPE,REF,ISTAT)
      CALL STDRDI(REF,'NAXIS',1,1,NN,NAXIS,KUN,KNUL,ISTAT)
      CALL STDRDI(REF,'NPIX',1,NAXIS,NN,NPIX,KUN,KNUL,ISTAT)
      CALL STDRDD(REF,'START',1,NAXIS,NN,START,KUN,KNUL,ISTAT)
      CALL STDRDD(REF,'STEP',1,NAXIS,NN,STEP,KUN,KNUL,ISTAT)
      IF (NAXIS.GT.1 .AND. NPIX(2).GT.1) THEN
          TYPE   = ' '
          CALL STDFND(REF,'WSTART',TYPE,NO,NB,ISTAT)
          IF (TYPE(1:1).EQ.' ') THEN
              CALL STTPUT(' Wrong reference image ',ISTAT)
              GO TO 30
          END IF
          CALL STDRDD(REF,'WSTART',1,NPIX(2),NN,WSTART,
     .                KUN,KNUL,ISTAT)
          CALL STDRDI(REF,'NPTOT',1,NPIX(2),NN,NPTOT,
     .                KUN,KNUL,ISTAT)
          CALL STDRDI(REF,'NORDER',1,NPIX(2),NN,NORDER,
     .                KUN,KNUL,ISTAT)
          NDIM   = 2

      ELSE
          NPIX(2) = 1
          NDIM   = 1
      END IF
C CALL EXIMAG(REFIMA,ISTAT)
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)
CC old: wrong, because CUNIT(1:16) is for unit of pixel values
CC          I1     = (I-1)*16 + 1
CC PN 8/99: add offset of 16: (I-1)*16 + 1 + 16 = 
          I1     = I*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 COPY(TID,ICOL(1),ICOL(2),NROW,MADRID(IP1),MADRID(IP2),N1)
      IF (N1.LE.3) THEN
          CALL STTPUT(' Not enough points in table ',ISTAT)
          GO TO 20
      END IF
C
C ... map image
C
CC old: wrong, because FLUX is the unit of the values in the image
CC      CUNIT(I2+1:I2+16) = 'FLUX'
CC PN 8/99: put this at the beginning
      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,D_R4_FORMAT,F_O_MODE,F_IMA_TYPE,
     .            NDIM,NPIX,START,STEP,IDENT,CUNIT,
     .            PNTR,IMNO,STATUS)
C
C ... interpolate from table values
C
      CALL INTERP(N1,MADRID(IP1),MADRID(IP2),NPIX(1),NPIX(2),
     +            MADRID(PNTR),START,STEP,WSTART,NPTOT,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)
      IF (NPIX(2).GT.1) THEN
          CALL STDWRD(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)
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 COPY(TID,ICOL1,ICOL2,NROW,X1,Y1,N1)
C
C copy relevant part of the table
C
      IMPLICIT NONE
C
      INTEGER TID, ICOL1, ICOL2
      INTEGER NROW, N1, I
C
      INTEGER STATUS
      LOGICAL ISEL, INULL1, INULL2
      REAL X,Y,X1(NROW),Y1(NROW)
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
          END IF
   10 CONTINUE
      RETURN
      END

      SUBROUTINE INTERP(NROW,X,Y,NPIX1,NPIX2,F,START,STEP,WSTART,NPTOT,
     +                  RMIN,RMAX)
C
C      IMPLICIT NONE
C call for each order to the interpolation routine
C
      INTEGER NROW,NPIX1,NPIX2,IORD,I,NP,FLAG,I0
      REAL X(NROW),F(NPIX1,NPIX2),RMIN,RMAX
      REAL Y(NROW),STARTX,STEPX,VAL,XP
      INTEGER NPTOT(100)
      DOUBLE PRECISION WSTART(100)
      DOUBLE PRECISION START(2),STEP(2)
C
      RMIN   = +99999.9
      RMAX   = -99999.9
      STARTX = START(1)
      STEPX  = STEP(1)
      NP     = NPIX1
C
C ... order loop
C
      DO 30 IORD = 1,NPIX2
          IF (NPIX2.GT.1) THEN
              STARTX = WSTART(IORD)
              NP     = NPTOT(IORD)
          END IF
C
C ... loop for each point
C
          FLAG = 1
          CALL INTEP(FLAG,STARTX,VAL,X,Y,NROW,I0)
          RMIN   = AMIN1(RMIN,VAL)
          RMAX   = AMAX1(RMAX,VAL)
          F(1,IORD) = VAL
          FLAG = 0
          DO 10 I = 2,NP
              XP     = STARTX + (I-1)*STEPX
              CALL INTEP(FLAG,XP,VAL,X,Y,NROW,I0)
              RMIN   = AMIN1(RMIN,VAL)
              RMAX   = AMAX1(RMAX,VAL)
              F(I,IORD) = VAL
   10     CONTINUE
          DO 20 I = NP + 1,NPIX1
              F(I,IORD) = 0.
   20     CONTINUE
   30 CONTINUE
      RETURN

      END

      SUBROUTINE INTEP(FLAG,XP,P,X,Y,N,I0)
C
C Interpolate a function P at a given value XP using
C      IMPLICIT NONE
C the table of values (X,F).
C Spline interpolation scheme based on Hermite polynomials.
C Ref.: U.S. Airforce Surveys in Geophysics, no.272
C       Publ of the Dominion Astrophys. Obs. 204, 1982
C For random values of XP use INTEP, or after the first
C call to INTEP for incresing/decreasing values of XP consistent
C with the vector X, use EINTP
C
C Arguments :
C XP Input argument value
C P Output interpolated value
C X Input vector of independent values
C F Input vector of dependent values
C N Input number of points in (X,F) vector
C I0 will be saved in the calling routine
C
C global variables
      INTEGER FLAG,N,I0
C local variables
      INTEGER I,N1,IUP, K
C 
C global variables
      REAL Y(N),X(N),XP,P
C local variables      
      REAL FP1,FP2,XPI,XPI1
      REAL LP1,LP2,L1,L2
C
      SAVE
C
      IUP = 0
      IF (X(2).LT.X(1)) IUP = 1 
      N1 = N - 1 

C FLAG is only 1 when INTEP is called the first time
      IF (FLAG .EQ. 1) THEN 
         I0 = 1
      ENDIF

C  This part has to be outside of IF FLAG.EQ.1 loop 
C  to avoid problems on DEC machines:
      IF ((XP.GT.X(N).AND.IUP.EQ.0) .OR. 
     +       (XP.LT.X(N).AND.IUP.EQ.1) .OR.
     +       (XP.LT.X(1).AND.IUP.EQ.0) .OR. 
     +       (XP.GT.X(1).AND.IUP.EQ.1)) THEN
             P = 0.0
             RETURN
      ENDIF    

      DO 40 I = 1,N
         IF ((XP.LT.X(I) .AND. IUP.EQ.0) .OR. 
     +        (XP.GT.X(I) .AND. IUP.EQ.1)) THEN
            K = I
            GOTO 50
         ENDIF
 40   CONTINUE 
      P=0.0
      RETURN     
      
 50   I  = K - 1 
      IF (I.NE.I0-1) THEN
         I0     = I + 1
         LP1    = 1./ (X(I)-X(I+1))
         LP2    = 1./ (X(I+1)-X(I))
      ENDIF
      
      IF (I.NE.1) THEN
         FP1    = (Y(I+1)-Y(I-1))/ (X(I+1)-X(I-1))
         FP2 = 0.0 
      ELSE
         FP1    = (Y(2)-Y(1))/ (X(2)-X(1))
      ENDIF
 
C leave the IF order like in original code      
      IF (I.GE.N1) THEN
         FP2    = (Y(N)-Y(N-1))/ (X(N)-X(N-1))
      ELSE
         FP2    = (Y(I+2)-Y(I))/ (X(I+2)-X(I))
      ENDIF
      
      XPI1   = XP - X(I+1) 
      XPI    = XP - X(I)
      L1     = XPI1*LP1
      L2     = XPI*LP2
      P      = Y(I)* (1.-2.*LP1*XPI)*L1*L1 +
     +     Y(I+1)* (1.-2.*LP2*XPI1)*L2*L2 + FP2*XPI1*L2*L2 +
     +     FP1*XPI*L1*L1

      RETURN
      END
