C @(#)tdcspl.for 17.1.1.1 (ESO-DMD) 01/25/02 17:47:13 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 C 010706 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 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(I2+1:I2+16) = 'FLUX' IDENT = 'TABLE: '//TABLE//'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