← Back to branch summary

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

« back to all changes in this revision

Viewing changes to prim/table/libsrc/tdcspl.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/tdcspl.for
 
1
C  @(#)tdcspl.for       19.1 (ESO-DMD) 02/25/03 14:11:16 
 
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
 
 
30
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 
31
C.COPYRIGHT: Copyright (c) 1987 European Southern Observatory,
 
32
C                                         all rights reserved
 
33
C
 
34
C.VERSION: 1.0  ESO-FORTRAN Conversion, AA  14:27 - 19 NOV 1987
 
35
C
 
36
C.LANGUAGE: F77+ESOext
 
37
C
 
38
C.AUTHOR: J.D.PONZ
 
39
C
 
40
C.IDENTIFICATION
 
41
C
 
42
C  program TDCSPL
 
43
C
 
44
C.PURPOSE
 
45
C
 
46
C  TABLE TO IMAGE CONVERSION
 
47
C  Execute the command
 
48
C  CONV/TAB image = table colx [z] refima SPLINE deg
 
49
C
 
50
C.KEYWORDS
 
51
C
 
52
C  table, image, conversion
 
53
C
 
54
C
 
55
C.ALGORITHM
 
56
C
 
57
C  spline interpolation
 
58
C
 
59
C  Ref.: P.Diercxx, 1982, Computer Graphics and Image Processing
 
60
C                         vol. 20, 171 - 184.
 
61
C
 
62
C.MODIFS
 
63
C     M Peron 071190 , in order to avoid a shift in the output frame
 
64
C  
 
65
C.VERSION
 
66
C 020807        last modif
 
67
C 
 
68
C-----------------------------------------------------------
 
69
C
 
70
      SUBROUTINE TDCSPL
 
71
      IMPLICIT NONE
 
72
C
 
73
      INTEGER     MADRID
 
74
      INTEGER PNVALS, IMNO
 
75
      INTEGER ICOL(3),NPTOT(100),NPIX(2)
 
76
      INTEGER NORDER(100)
 
77
      INTEGER NRDATA(1000)
 
78
      INTEGER STATUS, N1, IDEGR, NPAR, ISTAT, NN, NAXIS, NO, NB
 
79
      INTEGER NBY, TID, IS, NCOL, NROW, NSC, NAC, NAR, I, IREF
 
80
      INTEGER I1, I2, NDIM, KUN, KNUL
 
81
C 
 
82
      INTEGER*8  IP1, IP2, IP3, PNTR
 
83
C
 
84
      REAL S,SS,DELTA,YMAX,YMIN
 
85
      REAL STEP(2),START(2),RMIN,RMAX
 
86
      REAL CUTS(4),WSTART(100)
 
87
      REAL T(1000),C(1000)
 
88
C
 
89
      DOUBLE PRECISION DSTART(2), DSTEP(2)
 
90
C
 
91
      CHARACTER*16 MSG
 
92
      CHARACTER*80 TABLE,IMAGE,REFIMA
 
93
      CHARACTER*8 TYPE
 
94
      CHARACTER*17 COLREF(2)
 
95
      CHARACTER*60 WORK
 
96
      CHARACTER*72 IDENT
 
97
      CHARACTER*80 CUNIT
 
98
C
 
99
      INCLUDE 'MID_INCLUDE:TABLES.INC'
 
100
      COMMON /VMR/MADRID(1)
 
101
      INCLUDE 'MID_INCLUDE:TABLED.INC'
 
102
C
 
103
      DATA PNVALS/8/
 
104
      DATA MSG/'ERR:TCONVTBLxxxx'/
 
105
      DATA NPIX /1,1/
 
106
      DATA DSTEP /1.,1./
 
107
      DATA DSTART /0.,0./
 
108
C
 
109
C ... get input parameters + default
 
110
C
 
111
      CALL TDPGET(PNVALS,NPAR,ISTAT)
 
112
      IF (ISTAT.NE.0) GO TO 30
 
113
      IMAGE     = TPARBF(1)
 
114
      TABLE     = TPARBF(3)
 
115
      COLREF(1) = TPARBF(4)
 
116
      COLREF(2) = TPARBF(5)
 
117
      REFIMA = TPARBF(6)
 
118
      CALL STKRDR('INPUTR',1,1,NN,SS,KUN,KNUL,ISTAT)
 
119
      CALL STKRDR('INPUTR',2,1,NN,S,KUN,KNUL,ISTAT)
 
120
C      IDEGR = SS
 
121
C .................................... change parameters
 
122
      IDEGR = S
 
123
C 
 
124
      CALL STFOPN(REFIMA,D_R4_FORMAT,0,F_IMA_TYPE,IREF,ISTAT)
 
125
      CALL STDRDI(IREF,'NAXIS',1,1,NN,NAXIS,KUN,KNUL,ISTAT)
 
126
      CALL STDRDI(IREF,'NPIX',1,NAXIS,NN,NPIX,KUN,KNUL,ISTAT)
 
127
      CALL STDRDD(IREF,'START',1,NAXIS,NN,DSTART,KUN,KNUL,ISTAT)
 
128
      START(1) = DSTART(1)
 
129
      START(2) = DSTART(2)
 
130
      CALL STDRDD(IREF,'STEP',1,NAXIS,NN,DSTEP,KUN,KNUL,ISTAT)
 
131
      STEP(1) = DSTEP(1)
 
132
      STEP(2) = DSTEP(2)
 
133
      IF (NAXIS.GT.1 .AND. NPIX(2).GT.1) THEN
 
134
          TYPE   = ' '
 
135
          CALL STDFND(IREF,'WSTART',TYPE,NO,NB,ISTAT)
 
136
          IF (TYPE(1:1).EQ.' ') THEN
 
137
              CALL STTPUT(' Wrong reference image ',ISTAT)
 
138
              GO TO 30
 
139
          END IF
 
140
          CALL STDRDR
 
141
     .         (IREF,'WSTART',1,NPIX(2),NN,WSTART,KUN,KNUL,ISTAT)
 
142
          CALL STDRDI
 
143
     .         (IREF,'NPTOT',1,NPIX(2),NN,NPTOT,KUN,KNUL,ISTAT)
 
144
          CALL STDRDI
 
145
     .         (IREF,'NORDER',1,NPIX(2),NN,NORDER,KUN,KNUL,ISTAT)
 
146
          NDIM   = 2
 
147
 
 
148
      ELSE
 
149
          NPIX(2) = 1
 
150
          NDIM   = 1
 
151
      END IF
 
152
C
 
153
C ... init input table
 
154
C
 
155
      TID = -1
 
156
      CALL TBTOPN(TABLE,F_I_MODE,TID,STATUS)
 
157
      IF (ISTAT.NE.0) GO TO 30
 
158
      CALL TBIGET(TID,NCOL,NROW,NSC,NAC,NAR,ISTAT)
 
159
      IF (ISTAT.NE.0) GO TO 30
 
160
C
 
161
C ... find input column
 
162
C
 
163
      DO 10, I = 1, 2
 
164
          CALL TBCSER(TID,COLREF(I),ICOL(I),ISTAT)
 
165
          IF (ICOL(I).EQ.-1) THEN
 
166
              ISTAT  = ERRCOL
 
167
              GO TO 30
 
168
          END IF
 
169
          CALL TBUGET(TID,ICOL(I),WORK,ISTAT)
 
170
          I1     = (I-1)*16 + 1
 
171
          I2     = I1 + 15
 
172
          CUNIT(I1:I2) = WORK(1:16)
 
173
   10 CONTINUE
 
174
C
 
175
C ... copy selected values of the table into work space
 
176
C
 
177
      NBY = 4*NROW
 
178
      CALL TDMGET(NBY,IP1,IS)
 
179
      CALL TDMGET(NBY,IP2,IS)
 
180
      CALL TDMGET(NBY,IP3,IS)
 
181
      CALL ACOPY(TID,ICOL(1),ICOL(2),NROW,
 
182
     +          MADRID(IP1),MADRID(IP2),MADRID(IP3),YMIN,YMAX,N1)
 
183
      IF (N1.LT.IDEGR+1) THEN
 
184
          CALL STTPUT(' Not enough points in table ',ISTAT)
 
185
          GO TO 20
 
186
 
 
187
      END IF
 
188
C
 
189
C ... map image
 
190
C
 
191
      CUNIT(1:16) = 'FLUX'
 
192
      I2 = INDEX(TABLE,' ')
 
193
      IF (I2.lt.1) I2 = LEN(TABLE)
 
194
      IDENT  = 'TABLE: '//TABLE(1:I2)//'COLS. :'//COLREF(1)//COLREF(2)
 
195
      CALL STIPUT(IMAGE,10,1,1,NDIM,NPIX,DSTART,DSTEP,IDENT,CUNIT,
 
196
     +            PNTR,IMNO,STATUS)
 
197
C
 
198
C ... interpolate from table values
 
199
C
 
200
      DELTA  = YMAX - YMIN
 
201
      SS     = DELTA* (DELTA*0.5+YMIN)*SS*0.01
 
202
      CALL AINTER(N1,MADRID(IP1),MADRID(IP2),NPIX(1),NPIX(2),
 
203
     +            MADRID(PNTR),START,STEP,WSTART,NPTOT,RMIN,RMAX,
 
204
     +            MADRID(IP3),IDEGR,SS,T,C,NRDATA)
 
205
C
 
206
C ... write cuts
 
207
C
 
208
      CUTS(1) = RMIN
 
209
      CUTS(2) = RMAX
 
210
      CUTS(3) = RMIN
 
211
      CUTS(4) = RMAX
 
212
      CALL STDWRR(IMNO,'LHCUTS',CUTS,1,4,KUN,ISTAT)
 
213
      IF (NPIX(2).GT.1) THEN
 
214
          CALL STDWRR(IMNO,'WSTART',WSTART,1,NPIX(2),KUN,ISTAT)
 
215
          CALL STDWRI(IMNO,'NPTOT',NPTOT,1,NPIX(2),KUN,ISTAT)
 
216
          CALL STDWRI(IMNO,'NORDER',NORDER,1,NPIX(2),KUN,ISTAT)
 
217
      END IF
 
218
C
 
219
C ... free memory
 
220
C
 
221
   20 CALL TDMFRE(NBY,IP1,IS)
 
222
      CALL TDMFRE(NBY,IP2,IS)
 
223
      CALL TDMFRE(NBY,IP3,IS)
 
224
C
 
225
C ... end
 
226
C
 
227
      CALL DSCUPT(IMNO,IMNO,' ',ISTAT)
 
228
      CALL TBTCLO(TID,ISTAT)
 
229
   30 IF (ISTAT.NE.0) THEN
 
230
          WRITE (MSG(13:16),9000) ISTAT
 
231
          CALL TDERRR(ISTAT,MSG,STATUS)
 
232
      END IF
 
233
      RETURN
 
234
 9000 FORMAT (I4)
 
235
      END
 
236
 
 
237
      SUBROUTINE ACOPY(TID,ICOL1,ICOL2,NROW,X1,Y1,W,YMIN,YMAX,N1)
 
238
C
 
239
C copy relevant part of the table
 
240
C
 
241
      IMPLICIT NONE
 
242
C
 
243
      INTEGER TID, ICOL1, ICOL2
 
244
      INTEGER NROW, N1, I
 
245
      INTEGER STATUS
 
246
C
 
247
      REAL YMIN, YMAX
 
248
      REAL X1(NROW),Y1(NROW), X, Y
 
249
      REAL W(NROW)
 
250
C 
 
251
      LOGICAL ISEL, INULL1, INULL2
 
252
C
 
253
C
 
254
      N1     = 0
 
255
      DO 10, I = 1,NROW
 
256
          CALL TBSGET(TID,I,ISEL,STATUS)
 
257
          CALL TBERDR(TID,I,ICOL1,X,INULL1,STATUS)
 
258
          CALL TBERDR(TID,I,ICOL2,Y,INULL2,STATUS)
 
259
          IF (ISEL.AND. .NOT.INULL1 .AND. .NOT.INULL2) THEN
 
260
              N1     = N1 + 1
 
261
              X1(N1) = X
 
262
              Y1(N1) = Y
 
263
              W(N1)  = 1.
 
264
              IF (N1.EQ.1) THEN
 
265
                  YMIN   = Y
 
266
                  YMAX   = Y
 
267
              ELSE
 
268
                  IF (Y.GT.YMAX) YMAX   = Y
 
269
                  IF (Y.LT.YMIN) YMIN   = Y
 
270
              END IF
 
271
          END IF
 
272
   10 CONTINUE
 
273
      RETURN
 
274
      END
 
275
 
 
276
      SUBROUTINE AINTER(NROW,X,Y,NPIX1,NPIX2,F,START,STEP,WSTART,NPTOT,
 
277
     +                  RMIN,RMAX,W,IDEGR,S,T,C,NRDATA)
 
278
C
 
279
C call for each order to the interpolation routine
 
280
C
 
281
      IMPLICIT NONE
 
282
C 
 
283
      INTEGER MADRID
 
284
      INTEGER NROW,NPIX1, NPIX2,IDEGR
 
285
      INTEGER NPTOT(NPIX2)
 
286
      INTEGER ISTAT, NRDATA(1),NBYTES
 
287
      INTEGER IORD, NP, I, NK1, LD, NK, NPLUS,IERR
 
288
C 
 
289
      INTEGER*8   IPQ
 
290
C 
 
291
      REAL X(NROW),WSTART(NPIX2),F(NPIX1,NPIX2)
 
292
      REAL Y(NROW),START(2),STEP(2),RMIN,RMAX,S
 
293
      REAL W(1),T(1),C(1)
 
294
      REAL    DERIV, STARTX, STEPX, ENDX, VAL, XP, FP, FP0
 
295
      REAL    FPOLD, AMIN1, AMAX1
 
296
C 
 
297
      INCLUDE 'MID_INCLUDE:TABLES.INC'
 
298
      COMMON /VMR/MADRID(1)
 
299
      INCLUDE 'MID_INCLUDE:TABLED.INC'
 
300
C
 
301
      RMIN   = 0.
 
302
      RMAX   = 0.
 
303
      STARTX = START(1)
 
304
      STEPX  = STEP(1)
 
305
      NP     = NPIX1
 
306
      NBYTES = NROW*24  
 
307
C     NBYTES = NROW*6   M.Peron 23/07/93
 
308
      CALL TDMGET(NBYTES,IPQ,ISTAT)
 
309
C
 
310
C ... order loop
 
311
C
 
312
      DO 50, IORD = 1,NPIX2
 
313
          IF (NPIX2.GT.1) THEN
 
314
C              STARTX = WSTART(IORD)
 
315
              NP     = NPTOT(IORD)
 
316
          END IF
 
317
C
 
318
C ... loop for each point  W(1), IDEGR, S, T(1), C(1), NRDATA(1)
 
319
C
 
320
C          ENDX   = STARTX + (NP-1)*STEPX
 
321
           ENDX   = X(NROW)
 
322
           STARTX = X(1)
 
323
          CALL SMOOT(X,Y,W,MADRID(IPQ),NROW,STARTX,ENDX,IDEGR,S,NK,T,
 
324
     +               C,FP,0,IERR,NRDATA,FP0,FPOLD,NPLUS)
 
325
          IF (IERR.EQ.-2)
 
326
     .    CALL STTPUT('Polynomial approximation',ISTAT)
 
327
          IF (IERR.EQ.1) 
 
328
     .    CALL STTPUT('Local storage exceeded (S too small)',ISTAT)
 
329
          IF (IERR.EQ.2) 
 
330
     .    CALL STTPUT('Tolerance parameter too small',ISTAT)
 
331
          IF (IERR.EQ.3)
 
332
     .    CALL STTPUT('Maximum number of iter. exceeded',ISTAT)
 
333
          IF (IERR.EQ.10) 
 
334
     .    CALL STTPUT('Invalid input arguments',ISTAT)
 
335
          IF (IERR.GT.0) RETURN
 
336
          LD     = IDEGR + 1
 
337
          NK1    = NK - IDEGR - 1
 
338
          IF (NPIX2.EQ.1) STARTX = START(1) 
 
339
          DO 30, I = 1,NP
 
340
              XP     = STARTX + (I-1)*STEPX
 
341
   10         IF ((XP.LT.T(LD+1)) .OR. (LD.EQ.NK1)) GO TO 20
 
342
              LD     = LD + 1
 
343
              GO TO 10
 
344
   20         VAL    = DERIV(T,NK,C,NK1,0,XP,LD)
 
345
              RMIN   = AMIN1(RMIN,VAL)
 
346
              RMAX   = AMAX1(RMAX,VAL)
 
347
              F(I,IORD) = VAL
 
348
   30     CONTINUE
 
349
          DO 40, I = NP + 1,NPIX1
 
350
              F(I,IORD) = 0.
 
351
   40     CONTINUE
 
352
   50 CONTINUE
 
353
      CALL TDMFRE(NBYTES,IPQ,ISTAT)
 
354
      RETURN
 
355
 
 
356
      END