← Back to branch summary

~ubuntu-branches/debian/jessie/eso-midas/jessie

« back to all changes in this revision

Viewing changes to prim/table/libsrc/tdconv.for

  • Committer: Package Import Robot
  • Author(s): Ole Streicher
  • Date: 2014-04-22 14:44:58 UTC
  • Revision ID: package-import@ubuntu.com-20140422144458-okiwi1assxkkiz39
Tags: upstream-13.09pl1.2+dfsg
Import upstream version 13.09pl1.2+dfsg

Show diffs side-by-side

added added

removed removed

Lines of Context:
prim/table/libsrc/tdconv.for
 
1
C @(#)tdconv.for        19.1 (ESO-IPG) 02/25/03 14:11:11
 
2
C===========================================================================
 
3
C Copyright (C) 1995 European Southern Observatory (ESO)
 
4
C
 
5
C This program is free software; you can redistribute it and/or 
 
6
C modify it under the terms of the GNU General Public License as 
 
7
C published by the Free Software Foundation; either version 2 of 
 
8
C the License, or (at your option) any later version.
 
9
C
 
10
C This program is distributed in the hope that it will be useful,
 
11
C but WITHOUT ANY WARRANTY; without even the implied warranty of
 
12
C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 
13
C GNU General Public License for more details.
 
14
C
 
15
C You should have received a copy of the GNU General Public 
 
16
C License along with this program; if not, write to the Free 
 
17
C Software Foundation, Inc., 675 Massachusetts Ave, Cambridge, 
 
18
C MA 02139, USA.
 
19
C
 
20
C Correspondence concerning ESO-MIDAS should be addressed as follows:
 
21
C       Internet e-mail: midas@eso.org
 
22
C       Postal address: European Southern Observatory
 
23
C                       Data Management Division 
 
24
C                       Karl-Schwarzschild-Strasse 2
 
25
C                       D 85748 Garching bei Muenchen 
 
26
C                       GERMANY
 
27
C===========================================================================
 
28
C
 
29
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 
30
C.COPYRIGHT: Copyright (c) 1987 European Southern Observatory,
 
31
C                                         all rights reserved
 
32
C
 
33
C.VERSION: 1.0  ESO-FORTRAN Conversion, AA  14:22 - 19 NOV 1987
 
34
C
 
35
C.LANGUAGE: F77+ESOext
 
36
C
 
37
C.AUTHOR: J.D.PONZ
 
38
C  900219  KB, throw out SX calls...
 
39
C
 
40
C.IDENTIFICATION
 
41
C  program TDCONV
 
42
C
 
43
C.PURPOSE
 
44
C  Table to Image conversion
 
45
C  Execute the command
 
46
C  CONV/TAB image = table colx [z] refima method
 
47
C                 for method = SPLINE                 
 
48
C
 
49
C.KEYWORDS
 
50
C  table, image, conversion
 
51
C
 
52
C
 
53
C.ALGORITHM
 
54
C  Use table interface routines for I/O
 
55
C   The conversion from table format into image is done
 
56
C  by spline interpolation.
 
57
C  The interpolation scheme is based on Hermite polynomials
 
58
C  Ref .: Publ. of the Dominion Astrophys. Obs. 204, 1982.
 
59
C
 
60
C 021129                last modif
 
61
C 
 
62
C-----------------------------------------------------------
 
63
C
 
64
      SUBROUTINE TDCONV
 
65
      IMPLICIT NONE
 
66
C
 
67
      INTEGER     MADRID
 
68
      INTEGER PNVALS, STATUS, N1, NBY, IMNO,REF
 
69
      INTEGER ICOL(3),NPTOT(100),NPIX(2),KUN,KNUL
 
70
      INTEGER NORDER(100)
 
71
      INTEGER NPAR, ISTAT, NN, NAXIS, NO, NB, NDIM, TID, IS
 
72
      INTEGER NCOL, NROW, NSC, NAC, NAR, I, I1, I2
 
73
      INTEGER*8   PNTR, IP1, IP2
 
74
C
 
75
      DOUBLE PRECISION STEP(2),START(2)
 
76
      REAL   RMIN,RMAX, CUTS(4)
 
77
      DOUBLE PRECISION WSTART(100)
 
78
C
 
79
      CHARACTER*80 TABLE,IMAGE,REFIMA
 
80
      CHARACTER*17 COLREF(2)
 
81
      CHARACTER*60 WORK
 
82
      CHARACTER*16 MSG
 
83
      CHARACTER*8  TYPE
 
84
      CHARACTER*72 IDENT
 
85
      CHARACTER*80 CUNIT
 
86
C
 
87
      INCLUDE 'MID_INCLUDE:TABLES.INC'
 
88
      COMMON /VMR/MADRID(1)
 
89
      INCLUDE 'MID_INCLUDE:TABLED.INC'
 
90
C
 
91
      DATA PNVALS/8/
 
92
      DATA MSG/'ERR:TCONVTBLxxxx'/
 
93
      DATA NPIX /1,1/, START /0.,0./, STEP /1.,1./
 
94
C
 
95
C ... get input parameters + default
 
96
C
 
97
      CALL TDPGET(PNVALS,NPAR,ISTAT)
 
98
      IF (ISTAT.NE.0) GO TO 30
 
99
      IMAGE     = TPARBF(1)
 
100
      TABLE     = TPARBF(3)
 
101
      COLREF(1) = TPARBF(4)
 
102
      COLREF(2) = TPARBF(5)
 
103
      REFIMA = TPARBF(6)
 
104
      CALL STFOPN(REFIMA,D_R4_FORMAT,0,F_IMA_TYPE,REF,ISTAT)
 
105
      CALL STDRDI(REF,'NAXIS',1,1,NN,NAXIS,KUN,KNUL,ISTAT)
 
106
      CALL STDRDI(REF,'NPIX',1,NAXIS,NN,NPIX,KUN,KNUL,ISTAT)
 
107
      CALL STDRDD(REF,'START',1,NAXIS,NN,START,KUN,KNUL,ISTAT)
 
108
      CALL STDRDD(REF,'STEP',1,NAXIS,NN,STEP,KUN,KNUL,ISTAT)
 
109
      IF (NAXIS.GT.1 .AND. NPIX(2).GT.1) THEN
 
110
          TYPE   = ' '
 
111
          CALL STDFND(REF,'WSTART',TYPE,NO,NB,ISTAT)
 
112
          IF (TYPE(1:1).EQ.' ') THEN
 
113
              CALL STTPUT(' Wrong reference image ',ISTAT)
 
114
              GO TO 30
 
115
          END IF
 
116
          CALL STDRDD(REF,'WSTART',1,NPIX(2),NN,WSTART,
 
117
     .                KUN,KNUL,ISTAT)
 
118
          CALL STDRDI(REF,'NPTOT',1,NPIX(2),NN,NPTOT,
 
119
     .                KUN,KNUL,ISTAT)
 
120
          CALL STDRDI(REF,'NORDER',1,NPIX(2),NN,NORDER,
 
121
     .                KUN,KNUL,ISTAT)
 
122
          NDIM   = 2
 
123
 
 
124
      ELSE
 
125
          NPIX(2) = 1
 
126
          NDIM   = 1
 
127
      END IF
 
128
C CALL EXIMAG(REFIMA,ISTAT)
 
129
C
 
130
C ... init input table
 
131
C
 
132
      TID = -1
 
133
      CALL TBTOPN(TABLE,F_I_MODE,TID,STATUS)
 
134
      IF (ISTAT.NE.0) GO TO 30
 
135
      CALL TBIGET(TID,NCOL,NROW,NSC,NAC,NAR,ISTAT)
 
136
      IF (ISTAT.NE.0) GO TO 30
 
137
C
 
138
C ... find input column
 
139
C
 
140
      DO 10 I = 1,2
 
141
          CALL TBCSER(TID,COLREF(I),ICOL(I),ISTAT)
 
142
          IF (ICOL(I).EQ.-1) THEN
 
143
              ISTAT  = ERRCOL
 
144
              GO TO 30
 
145
          END IF
 
146
          CALL TBUGET(TID,ICOL(I),WORK,ISTAT)
 
147
CC old: wrong, because CUNIT(1:16) is for unit of pixel values
 
148
CC          I1     = (I-1)*16 + 1
 
149
CC PN 8/99: add offset of 16: (I-1)*16 + 1 + 16 = 
 
150
          I1     = I*16 + 1
 
151
          I2     = I1 + 15
 
152
          CUNIT(I1:I2) = WORK(1:16)
 
153
   10 CONTINUE
 
154
C
 
155
C ... copy selected values of the table into work space
 
156
C
 
157
      NBY = 4*NROW
 
158
      CALL TDMGET(NBY, IP1, IS)
 
159
      CALL TDMGET(NBY, IP2, IS)
 
160
      CALL COPY(TID,ICOL(1),ICOL(2),NROW,MADRID(IP1),MADRID(IP2),N1)
 
161
      IF (N1.LE.3) THEN
 
162
          CALL STTPUT(' Not enough points in table ',ISTAT)
 
163
          GO TO 20
 
164
      END IF
 
165
C
 
166
C ... map image
 
167
C
 
168
CC old: wrong, because FLUX is the unit of the values in the image
 
169
CC      CUNIT(I2+1:I2+16) = 'FLUX'
 
170
CC PN 8/99: put this at the beginning
 
171
      CUNIT(1:16) = 'FLUX'
 
172
      I2 = INDEX(TABLE,' ')
 
173
      IF (I2.lt.1) I2 = LEN(TABLE)
 
174
      IDENT  = 'TABLE: '//TABLE(1:I2)//'COLS. :'//COLREF(1)//COLREF(2)
 
175
      CALL STIPUT(IMAGE,D_R4_FORMAT,F_O_MODE,F_IMA_TYPE,
 
176
     .            NDIM,NPIX,START,STEP,IDENT,CUNIT,
 
177
     .            PNTR,IMNO,STATUS)
 
178
C
 
179
C ... interpolate from table values
 
180
C
 
181
      CALL INTERP(N1,MADRID(IP1),MADRID(IP2),NPIX(1),NPIX(2),
 
182
     +            MADRID(PNTR),START,STEP,WSTART,NPTOT,RMIN,RMAX)
 
183
C
 
184
C ... write cuts
 
185
C
 
186
      CUTS(1) = RMIN
 
187
      CUTS(2) = RMAX
 
188
      CUTS(3) = RMIN
 
189
      CUTS(4) = RMAX
 
190
      CALL STDWRR(IMNO,'LHCUTS',CUTS,1,4,KUN,ISTAT)
 
191
      IF (NPIX(2).GT.1) THEN
 
192
          CALL STDWRD(IMNO,'WSTART',WSTART,1,NPIX(2),KUN,ISTAT)
 
193
          CALL STDWRI(IMNO,'NPTOT',NPTOT,1,NPIX(2),KUN,ISTAT)
 
194
          CALL STDWRI(IMNO,'NORDER',NORDER,1,NPIX(2),KUN,ISTAT)
 
195
      END IF
 
196
C
 
197
C ... free memory
 
198
C
 
199
   20 CALL TDMFRE(NBY, IP1, IS)
 
200
      CALL TDMFRE(NBY, IP2, IS)
 
201
C
 
202
C ... end
 
203
C
 
204
      CALL DSCUPT(IMNO,IMNO,' ',ISTAT)
 
205
      CALL TBTCLO(TID,ISTAT)
 
206
   30 IF (ISTAT.NE.0) THEN
 
207
          WRITE (MSG(13:16),9000) ISTAT
 
208
          CALL TDERRR(ISTAT,MSG,STATUS)
 
209
      END IF
 
210
      RETURN
 
211
 9000 FORMAT (I4)
 
212
      END
 
213
 
 
214
      SUBROUTINE COPY(TID,ICOL1,ICOL2,NROW,X1,Y1,N1)
 
215
C
 
216
C copy relevant part of the table
 
217
C
 
218
      IMPLICIT NONE
 
219
C
 
220
      INTEGER TID, ICOL1, ICOL2
 
221
      INTEGER NROW, N1, I
 
222
C
 
223
      INTEGER STATUS
 
224
      LOGICAL ISEL, INULL1, INULL2
 
225
      REAL X,Y,X1(NROW),Y1(NROW)
 
226
C
 
227
      N1     = 0
 
228
      DO 10 I = 1,NROW
 
229
          CALL TBSGET(TID,I,ISEL,STATUS)
 
230
          CALL TBERDR(TID,I,ICOL1,X,INULL1,STATUS)
 
231
          CALL TBERDR(TID,I,ICOL2,Y,INULL2,STATUS)
 
232
          IF (ISEL .AND. .NOT.INULL1 .AND. .NOT.INULL2) THEN
 
233
              N1     = N1 + 1
 
234
              X1(N1) = X
 
235
              Y1(N1) = Y
 
236
          END IF
 
237
   10 CONTINUE
 
238
      RETURN
 
239
      END
 
240
 
 
241
      SUBROUTINE INTERP(NROW,X,Y,NPIX1,NPIX2,F,START,STEP,WSTART,NPTOT,
 
242
     +                  RMIN,RMAX)
 
243
C
 
244
C      IMPLICIT NONE
 
245
C call for each order to the interpolation routine
 
246
C
 
247
      INTEGER NROW,NPIX1,NPIX2,IORD,I,NP,FLAG,I0
 
248
      REAL X(NROW),F(NPIX1,NPIX2),RMIN,RMAX
 
249
      REAL Y(NROW),STARTX,STEPX,VAL,XP
 
250
      INTEGER NPTOT(100)
 
251
      DOUBLE PRECISION WSTART(100)
 
252
      DOUBLE PRECISION START(2),STEP(2)
 
253
C
 
254
      RMIN   = +99999.9
 
255
      RMAX   = -99999.9
 
256
      STARTX = START(1)
 
257
      STEPX  = STEP(1)
 
258
      NP     = NPIX1
 
259
C
 
260
C ... order loop
 
261
C
 
262
      DO 30 IORD = 1,NPIX2
 
263
          IF (NPIX2.GT.1) THEN
 
264
              STARTX = WSTART(IORD)
 
265
              NP     = NPTOT(IORD)
 
266
          END IF
 
267
C
 
268
C ... loop for each point
 
269
C
 
270
          FLAG = 1
 
271
          CALL INTEP(FLAG,STARTX,VAL,X,Y,NROW,I0)
 
272
          RMIN   = AMIN1(RMIN,VAL)
 
273
          RMAX   = AMAX1(RMAX,VAL)
 
274
          F(1,IORD) = VAL
 
275
          FLAG = 0
 
276
          DO 10 I = 2,NP
 
277
              XP     = STARTX + (I-1)*STEPX
 
278
              CALL INTEP(FLAG,XP,VAL,X,Y,NROW,I0)
 
279
              RMIN   = AMIN1(RMIN,VAL)
 
280
              RMAX   = AMAX1(RMAX,VAL)
 
281
              F(I,IORD) = VAL
 
282
   10     CONTINUE
 
283
          DO 20 I = NP + 1,NPIX1
 
284
              F(I,IORD) = 0.
 
285
   20     CONTINUE
 
286
   30 CONTINUE
 
287
      RETURN
 
288
 
 
289
      END
 
290
 
 
291
      SUBROUTINE INTEP(FLAG,XP,P,X,Y,N,I0)
 
292
C
 
293
C Interpolate a function P at a given value XP using
 
294
C      IMPLICIT NONE
 
295
C the table of values (X,F).
 
296
C Spline interpolation scheme based on Hermite polynomials.
 
297
C Ref.: U.S. Airforce Surveys in Geophysics, no.272
 
298
C       Publ of the Dominion Astrophys. Obs. 204, 1982
 
299
C For random values of XP use INTEP, or after the first
 
300
C call to INTEP for incresing/decreasing values of XP consistent
 
301
C with the vector X, use EINTP
 
302
C
 
303
C Arguments :
 
304
C XP Input argument value
 
305
C P Output interpolated value
 
306
C X Input vector of independent values
 
307
C F Input vector of dependent values
 
308
C N Input number of points in (X,F) vector
 
309
C I0 will be saved in the calling routine
 
310
C
 
311
C global variables
 
312
      INTEGER FLAG,N,I0
 
313
C local variables
 
314
      INTEGER I,N1,IUP, K
 
315
C 
 
316
C global variables
 
317
      REAL Y(N),X(N),XP,P
 
318
C local variables      
 
319
      REAL FP1,FP2,XPI,XPI1
 
320
      REAL LP1,LP2,L1,L2
 
321
C
 
322
      SAVE
 
323
C
 
324
      IUP = 0
 
325
      IF (X(2).LT.X(1)) IUP = 1 
 
326
      N1 = N - 1 
 
327
 
 
328
C FLAG is only 1 when INTEP is called the first time
 
329
      IF (FLAG .EQ. 1) THEN 
 
330
         I0 = 1
 
331
      ENDIF
 
332
 
 
333
C  This part has to be outside of IF FLAG.EQ.1 loop 
 
334
C  to avoid problems on DEC machines:
 
335
      IF ((XP.GT.X(N).AND.IUP.EQ.0) .OR. 
 
336
     +       (XP.LT.X(N).AND.IUP.EQ.1) .OR.
 
337
     +       (XP.LT.X(1).AND.IUP.EQ.0) .OR. 
 
338
     +       (XP.GT.X(1).AND.IUP.EQ.1)) THEN
 
339
             P = 0.0
 
340
             RETURN
 
341
      ENDIF    
 
342
 
 
343
      DO 40 I = 1,N
 
344
         IF ((XP.LT.X(I) .AND. IUP.EQ.0) .OR. 
 
345
     +        (XP.GT.X(I) .AND. IUP.EQ.1)) THEN
 
346
            K = I
 
347
            GOTO 50
 
348
         ENDIF
 
349
 40   CONTINUE 
 
350
      P=0.0
 
351
      RETURN     
 
352
      
 
353
 50   I  = K - 1 
 
354
      IF (I.NE.I0-1) THEN
 
355
         I0     = I + 1
 
356
         LP1    = 1./ (X(I)-X(I+1))
 
357
         LP2    = 1./ (X(I+1)-X(I))
 
358
      ENDIF
 
359
      
 
360
      IF (I.NE.1) THEN
 
361
         FP1    = (Y(I+1)-Y(I-1))/ (X(I+1)-X(I-1))
 
362
         FP2 = 0.0 
 
363
      ELSE
 
364
         FP1    = (Y(2)-Y(1))/ (X(2)-X(1))
 
365
      ENDIF
 
366
 
 
367
C leave the IF order like in original code      
 
368
      IF (I.GE.N1) THEN
 
369
         FP2    = (Y(N)-Y(N-1))/ (X(N)-X(N-1))
 
370
      ELSE
 
371
         FP2    = (Y(I+2)-Y(I))/ (X(I+2)-X(I))
 
372
      ENDIF
 
373
      
 
374
      XPI1   = XP - X(I+1) 
 
375
      XPI    = XP - X(I)
 
376
      L1     = XPI1*LP1
 
377
      L2     = XPI*LP2
 
378
      P      = Y(I)* (1.-2.*LP1*XPI)*L1*L1 +
 
379
     +     Y(I+1)* (1.-2.*LP2*XPI1)*L2*L2 + FP2*XPI1*L2*L2 +
 
380
     +     FP1*XPI*L1*L1
 
381
 
 
382
      RETURN
 
383
      END