C @(#)tdconv.for 17.1.1.1 (ESO-IPG) 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: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 C program TDCONV C C.PURPOSE C 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 C table, image, conversion C C C.ALGORITHM C 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 010110 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 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' 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, . 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