← Back to branch summary

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

« back to all changes in this revision

Viewing changes to stdred/long/src/rebirbr.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:
stdred/long/src/rebirbr.for
 
1
C @(#)rebirbr.for       19.1 (ES0-DMD) 02/25/03 14:25:23
 
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 Massachusetss Ave, Cambridge, 
 
18
C MA 02139, USA.
 
19
C
 
20
C Corresponding 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 @(#)rebirbr.for       19.1  (ESO)  02/25/03  14:25:23 
 
30
C +++++++++++++++++++++++++++++++++++++++++++++++++++++++
 
31
C .COPYRIGHT   (C) 1993 European Southern Observatory    
 
32
C .IDENT       .prg                                      
 
33
C .AUTHORS     Pascal Ballester (ESO/Garching)           
 
34
C              Cristian Levin   (ESO/La Silla)           
 
35
C .KEYWORDS    Spectroscopy, Long-Slit                   
 
36
C .PURPOSE                                               
 
37
C .VERSION     1.0  Package Creation  17-MAR-1993        
 
38
C -------------------------------------------------------
 
39
 
 
40
C @(#)rebirbr.for       1.1 (ESO-IPG) 2/17/93 16:46:20
 
41
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 
42
C.COPYRIGHT: Copyright (c) 1987 European Southern Observatory,
 
43
C                                         all rights reserved
 
44
C
 
45
C.VERSION: 1.0  ESO-FORTRAN Conversion, AA  8:58 - 9 DEC 1987
 
46
C
 
47
C.LANGUAGE: F77+ESOext
 
48
C
 
49
C.IDENTIFICATION
 
50
C
 
51
C  SUBROUTINES TDRBII MODIFIED TO REBIN A 2D FRAME ROW BY ROW
 
52
C                                          (M.Pierre Nov. 1989)
 
53
C.KEYWORDS
 
54
C
 
55
C  IMAGES, NONLINEAR REBIN, INTERPOLATION, INTEGRATION
 
56
C
 
57
C.PURPOSE
 
58
C
 
59
C  NONLINEAR REBIN IMAGE TO IMAGE IN 2 DIMENSIONS
 
60
C
 
61
C.ALGORITHM
 
62
C
 
63
C  Rebin input data YI sampled at XI with binwidth WI to
 
64
C  output data YO sampled at XO with binwidth WO.
 
65
C
 
66
C  Restriction:  Input and output independent variables are
 
67
C                assumed to be monotonic increasing or decreasing
 
68
C
 
69
C  Several FUNCTal relations between input and output independent
 
70
C  variables are foreseen, additional ones may be added in dummy call
 
71
C  U01 as follows:
 
72
C
 
73
C SUBROUTINE REB_U01 (DV,P,NP,NPC)
 
74
C   REAL*8 DV,P(NP)
 
75
C   DV = ...........
 
76
C   RETURN
 
77
C END
 
78
C
 
79
C  where DV is input-output value, P(NP) the parameter array and
 
80
C  NPC the nuber of active parameters in p(np).
 
81
C
 
82
C  According to the definition of IMAGES, pixels have constant
 
83
C  bin width (STEP) [WI;WO], are not overlapping and have filling
 
84
C  factors of exactly one. World coordinates [XI;XO] and flux values
 
85
C  [YI;YO] are pertinent
 
86
C  to the center of a pixel. In flux conservation mode (the default)
 
87
C  the flux is interpreted as the area under a smooth curve between
 
88
C  the boundaries of the pixel (i.e. a histogram representation of
 
89
C  a smooth curve).
 
90
C
 
91
C  Detailed flux conservation obtained by REAL*8 spline interpolation
 
92
C  and  integration of input data over the projected
 
93
C  bin WI for fractions of input pixels and summation of flux in
 
94
C  entire input pixels covered within the valid domain of input data
 
95
C  i.e. from x(1)-0.5*step to x(n)+0.5*step. Several options available
 
96
C  to handle undefined input at the extremes and for NULLVALUE assigned
 
97
C  parts of input array.
 
98
C
 
99
C   Optimized for long arrays by search for valid range, optimum
 
100
C  direction or DO LOOPS and updating of loop boundaries - therefore
 
101
C  overdoing it for small arrays.
 
102
C
 
103
C.INPUT
 
104
C
 
105
C  Keys:  IN_A/C/1/8         input data array
 
106
C         IN_B/C/1/8         reference frame or dim of output
 
107
C         CFUNC/C/1/8        FUNCTal relation between indep.
 
108
C         INPUTD/D/1/12      parameters or FUNCT
 
109
C         COPT/C/1/8         options
 
110
C
 
111
C.OUTPUT
 
112
C
 
113
C  Key:   OUT_A/C/1/8        output data array
 
114
C
 
115
C.MODIFS
 
116
C M Peron sep 90 : remove stupid limit to the size of the output frame!!!!!
 
117
C--------------------------------------------------------------------------
 
118
C
 
119
       PROGRAMM REBIRBR
 
120
      IMPLICIT NONE
 
121
      INTEGER   NFPAR
 
122
      PARAMETER (NFPAR = 12)
 
123
C
 
124
      INTEGER    MADRID
 
125
      INTEGER    I,IACT
 
126
      INTEGER*8  IP1,IP2,IP3,IPNTR
 
127
      INTEGER    IS,ISTAT
 
128
      INTEGER*8  JP1,JP2,JP3,JPNTR
 
129
      INTEGER    NAXISI,NAXISO,KUN,KNUL
 
130
      INTEGER    NFPACT,NFUNC
 
131
      INTEGER    NI,NINT,NN,NO,NBY,NBO,NCOL,NTOT
 
132
      INTEGER NPIXI(3),NPIXO(3),INOP,IMNO,IMNI,IREF
 
133
       INTEGER JJ,TID,NROW,NSC,NAC,NAR,ICOL(12),STATUS
 
134
      REAL    RMIN, RMAX, CUTS(4)
 
135
      DOUBLE PRECISION STEPI(3),STARTI(3)
 
136
        REAL X0O,WXO,X0I,WXI
 
137
      DOUBLE PRECISION DPARM(NFPAR)
 
138
      DOUBLE PRECISION DSTART(3), DSTEP(3),YPOS
 
139
      CHARACTER*3 FUN(9)
 
140
      CHARACTER*60 INPIMA,OUTIMA,REFIMA,TBCOEFF
 
141
      CHARACTER*80 FUNCT,OPTION
 
142
      CHARACTER*72 IDENTO,IDENTI
 
143
      CHARACTER*80 CUNITO,CUNITI
 
144
       LOGICAL NULL(12)
 
145
C
 
146
      INCLUDE 'MID_INCLUDE:TABLES.INC'
 
147
      COMMON /VMR/MADRID(1)
 
148
      INCLUDE 'MID_INCLUDE:TABLED.INC'
 
149
C
 
150
      DATA FUN/'LIN','POL','INV','EXP','DEX','LOG','DLG','IPO','U01'/
 
151
 
 
152
C
 
153
C..... GET INTO MIDAS
 
154
       CALL STSPRO('REBLL')
 
155
C ... get input parameters
 
156
C
 
157
      CALL STKRDC('OUT_A',1,1,60,IACT,OUTIMA,KUN,KNUL,ISTAT)
 
158
      CALL STKRDC('IN_A', 1,1,60,IACT,INPIMA,KUN,KNUL,ISTAT)
 
159
      CALL STKRDC('IN_B', 1,1,60,IACT,REFIMA,KUN,KNUL,ISTAT)
 
160
      CALL STKRDC('CFUNC',1,1,8,IACT,FUNCT,KUN,KNUL,ISTAT)
 
161
      CALL STKRDC('IN_C', 1,1,60,IACT,TBCOEFF,KUN,KNUL,ISTAT)
 
162
      CALL STKRDC('COPT',1,1,8,IACT,OPTION,KUN,KNUL,ISTAT)
 
163
 
 
164
C
 
165
C ... det. interpolation/integration options
 
166
C  !   Not implemented
 
167
      INOP   = 4
 
168
C      IF (OPTION(1:3).EQ.'PIX') INOP   = 1
 
169
C      IF (OPTION(1:3).EQ.'LIN') INOP   = 2
 
170
C      IF (OPTION(1:3).EQ.'SPG') INOP   = 3
 
171
C .......................................... the actual values are:
 
172
      CALL FORUPC(OPTION,OPTION)
 
173
      IF (OPTION(1:3).EQ.'PIX') INOP   = 2
 
174
      IF (OPTION(1:3).EQ.'LIN') INOP   = 3
 
175
      IF (OPTION(1:3).EQ.'SPG') INOP   = 1
 
176
C .......................................... Bitte, bitte
 
177
C
 
178
C ... translate FUNCT into FUNCT number
 
179
C
 
180
      CALL FORUPC(FUNCT,FUNCT)
 
181
      NFUNC  = 0
 
182
      DO 10 I = 1,9
 
183
          IF (FUNCT(1:3).EQ.FUN(I)) NFUNC  = I
 
184
   10 CONTINUE
 
185
      IF (NFUNC.EQ.0) THEN
 
186
          CALL STTPUT(' Specified FUNCT non-existent...',ISTAT)
 
187
          GO TO 60
 
188
      END IF
 
189
C
 
190
C..... init the input table of dispersion coefficients
 
191
        CALL TBTOPN(TBCOEFF,F_U_MODE,TID,ISTAT)
 
192
       CALL TBIGET(TID,NCOL,NROW,NSC,NAC,NAR,ISTAT)
 
193
       NFPACT = NCOL-1
 
194
 
 
195
       DO 40 I = 1,NFPACT
 
196
       ICOL(I) = I+1
 
197
40       CONTINUE       
 
198
 
 
199
C
 
200
C ... get parameters form reference image
 
201
C
 
202
      CALL STFOPN(REFIMA,D_R4_FORMAT,0,F_IMA_TYPE,IREF,ISTAT)
 
203
      CALL STDRDI(IREF,'NAXIS',1,1,NN,NAXISO,KUN,KNUL,ISTAT)
 
204
      CALL STDRDI(IREF,'NPIX',1,NAXISO,NN,NPIXO,KUN,KNUL,ISTAT)
 
205
      CALL STDRDD(IREF,'START',1,NAXISO,NN,DSTART,KUN,KNUL,ISTAT)
 
206
      CALL STDRDD(IREF,'STEP',1,NAXISO,NN,DSTEP,KUN,KNUL,ISTAT)
 
207
      CALL STDRDC(IREF,'IDENT',1,1,72,NN,IDENTO,KUN,KNUL,ISTAT)
 
208
      CALL STDRDC(IREF,'CUNIT',1,1,80,NN,CUNITO,KUN,KNUL,ISTAT)
 
209
      CALL STFCLO(IREF,ISTAT)
 
210
       NTOT = NPIXO(1)*NPIXO(2)
 
211
C
 
212
C ... get input image
 
213
C
 
214
      CALL STIGET(INPIMA,D_R4_FORMAT,F_I_MODE,F_IMA_TYPE,
 
215
     .2,NAXISI,NPIXI,STARTI,STEPI,IDENTI,CUNITI,IPNTR,IMNI,ISTAT)
 
216
 
 
217
C..... check dimensions
 
218
       IF(NPIXI(2).NE.NROW) THEN
 
219
       CALL STTPUT('Input table and image have not 
 
220
     +       the same number of rows',ISTAT)
 
221
       GO TO 50
 
222
       ELSE 
 
223
       IF(NPIXI(2).NE.NPIXO(2)) THEN
 
224
       CALL STTPUT('Input and output images have not 
 
225
     +  the same number of rows',ISTAT)
 
226
       GO TO 50
 
227
             ENDIF
 
228
       ENDIF
 
229
       CALL TBRRDD(TID,1,1,1,YPOS,NULL,STATUS)
 
230
       IF (ABS(YPOS-DSTART(2)).GT.0.001) THEN
 
231
        CALL STTPUT(' WARNING:',ISTAT)
 
232
       CALL STTPUT('Y_position recorded in the 1st line
 
233
     +        of input table is not the same',ISTAT)
 
234
        CALL STTPUT('as Y_start of the  output image',ISTAT)
 
235
       ENDIF
 
236
C ... map output image
 
237
C
 
238
      CALL STIPUT(OUTIMA,D_R4_FORMAT,F_O_MODE,F_IMA_TYPE,
 
239
     .NAXISO,NPIXO,DSTART,DSTEP,IDENTO,CUNITO,JPNTR,IMNO,ISTAT)
 
240
C
 
241
 
 
242
C ... fill the arrays XI,WI and XO,WO
 
243
C
 
244
      NI     = NPIXI(1)
 
245
      NBY    = 8*NI
 
246
      CALL TDMGET(NBY,IP1,IS)
 
247
      CALL TDMGET(NBY,IP2,IS)
 
248
      CALL TDMGET(NBY,IP3,IS)
 
249
      X0O     = STARTI(1)
 
250
      WXO     = STEPI(1)
 
251
C
 
252
      NO     = NPIXO(1)
 
253
      NBO    = 8*NO
 
254
      CALL TDMGET(NBO,JP1,IS)
 
255
      CALL TDMGET(NBO,JP2,IS)
 
256
      CALL TDMGET(NBO,JP3,IS)
 
257
      X0I     = DSTART(1)
 
258
      WXI     = DSTEP(1)
 
259
 
 
260
 
 
261
C..... DO LOOP ON THE IMAGE ROWS
 
262
 
 
263
       DO 100       JJ = 1,NROW
 
264
 
 
265
      CALL IMVAL3(NI,X0O,WXO,MADRID(IPNTR),MADRID(IP1),MADRID(IP2),
 
266
     +             MADRID(IP3),JJ,NTOT)
 
267
      CALL IMVAL2(NO,X0I,WXI,MADRID(JP1),MADRID(JP2))
 
268
C
 
269
C...... Read dispersion coefficients for row JJ
 
270
       CALL TBRRDD(TID,JJ,NFPACT,ICOL,DPARM,NULL,STATUS)
 
271
       WRITE(6,*) 'ROW',JJ
 
272
C       TYPE *,'COEFF',DPARM(1),DPARM(2),DPARM(3),DPARM(4),DPARM(5)
 
273
 
 
274
C  ...  Now, at last, we call the interpolating (and integrating) routines:
 
275
C
 
276
C  !!!! Watch out for VECI,VECD arrays in INTERP, dimension => NINT !!!!!!
 
277
C
 
278
      NINT   = 8
 
279
      CALL REBMET(NI,MADRID(IP1),MADRID(IP3),MADRID(IP2),NO,
 
280
     +             MADRID(JP1),MADRID(JP2),NFUNC,NFPAR,NFPACT,DPARM,
 
281
     +             INOP,NINT,MADRID(JP3),RMIN,RMAX)
 
282
C
 
283
       CALL IMVAL5(NO,JJ,MADRID(JP3),MADRID(JPNTR),NTOT)
 
284
 
 
285
 
 
286
C
 
287
100       CONTINUE
 
288
 
 
289
      CALL TDMFRE(NBY,IP1,IS)
 
290
      CALL TDMFRE(NBY,IP2,IS)
 
291
      CALL TDMFRE(NBY,IP3,IS)
 
292
      CALL TDMFRE(NBO,JP1,IS)
 
293
      CALL TDMFRE(NBO,JP2,IS)
 
294
 
 
295
C ... write cuts
 
296
C
 
297
      CUTS(1) = RMIN
 
298
      CUTS(2) = RMAX
 
299
      CUTS(3) = RMIN
 
300
      CUTS(4) = RMAX
 
301
      CALL STDWRR(IMNO,'LHCUTS',CUTS,1,4,KUN,ISTAT)
 
302
C
 
303
C ... free memory
 
304
C
 
305
 
 
306
C ... end
 
307
C
 
308
      CALL STFCLO(IMNO,ISTAT)
 
309
   50 CALL STFCLO(IMNI,ISTAT)
 
310
       CALL TBTCLO(TID,ISTAT)
 
311
   60 CALL STSEPI
 
312
      END
 
313
 
 
314
C**************************************************************
 
315
      SUBROUTINE IMVAL2(NP,XS,WS,X,W)
 
316
C
 
317
C  Fill arrays X,W for image world coordinates and pixel size data
 
318
C
 
319
      IMPLICIT NONE
 
320
      INTEGER  I, NP
 
321
      REAL     XS, WS
 
322
      DOUBLE PRECISION X(NP),W(NP),WSS,XSS
 
323
C
 
324
      WSS    = WS
 
325
      XSS    = XS
 
326
      DO 10 I = 1,NP
 
327
          X(I)   = XSS + WSS* (I-1)
 
328
          W(I)   = WSS
 
329
   10 CONTINUE
 
330
      RETURN
 
331
      END
 
332
 
 
333
C------------------------------------------------------------------
 
334
      SUBROUTINE IMVAL3(NP,XS,WS,YY,X,W,Y,NL,NTOT)
 
335
C
 
336
C  Fill arrays X,W for image world coordinates and pixel size data
 
337
C  Convert YY(real*4) into Y(real*8)
 
338
C
 
339
      IMPLICIT NONE
 
340
      INTEGER  NP, I,NTOT
 
341
       INTEGER NL,IPT
 
342
      DOUBLE PRECISION X(NP),W(NP),Y(NP),WSS,XSS
 
343
      REAL     XS, WS
 
344
      REAL YY(NTOT)
 
345
C
 
346
      WSS    = WS
 
347
      XSS    = XS
 
348
      DO 10 I = 1,NP
 
349
       IPT = (NL-1)*NP+I
 
350
          X(I)   = XSS + WSS* (I-1)
 
351
          W(I)   = WSS
 
352
          Y(I)   = YY(IPT)
 
353
   10 CONTINUE
 
354
      RETURN
 
355
      END
 
356
 
 
357
C----------------------------------------------------------
 
358
       SUBROUTINE IMVAL5(N,NL,X,Y,NTOT)
 
359
C       FILLS THE BUFFER
 
360
       INTEGER I,N,NL,NTOT,NPT
 
361
 
 
362
       REAL*4 Y(NTOT),X(N)
 
363
       DO 10 I = 1,N
 
364
       NPT = (NL-1)*N+I
 
365
       Y(NPT) = X(I)
 
366
10        CONTINUE
 
367
       
 
368
       RETURN
 
369
       END
 
370
 
 
371