← Back to branch summary

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

« back to all changes in this revision

Viewing changes to prim/general/src/genxx2.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/general/src/genxx2.for
 
1
C===========================================================================
 
2
C Copyright (C) 1995-2008 European Southern Observatory (ESO)
 
3
C
 
4
C This program is free software; you can redistribute it and/or 
 
5
C modify it under the terms of the GNU General Public License as 
 
6
C published by the Free Software Foundation; either version 2 of 
 
7
C the License, or (at your option) any later version.
 
8
C
 
9
C This program is distributed in the hope that it will be useful,
 
10
C but WITHOUT ANY WARRANTY; without even the implied warranty of
 
11
C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 
12
C GNU General Public License for more details.
 
13
C
 
14
C You should have received a copy of the GNU General Public 
 
15
C License along with this program; if not, write to the Free 
 
16
C Software Foundation, Inc., 675 Massachusetss Ave, Cambridge, 
 
17
C MA 02139, USA.
 
18
C
 
19
C Correspondence concerning ESO-MIDAS should be addressed as follows:
 
20
C       Internet e-mail: midas@eso.org
 
21
C       Postal address: European Southern Observatory
 
22
C                       Data Management Division 
 
23
C                       Karl-Schwarzschild-Strasse 2
 
24
C                       D 85748 Garching bei Muenchen 
 
25
C                       GERMANY
 
26
C===========================================================================
 
27
C
 
28
C +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 
29
C 
 
30
C.VERSION
 
31
C 1.0           from 1986 orig code move to ST interfaces
 
32
C 
 
33
C 080122        last modif
 
34
C 
 
35
C ----------------------------------------------------------------------
 
36
C 
 
37
      SUBROUTINE GXMOVIT(A,B,AOFF,BOFF,NOPIX,NOLINS,AXSIZE,BXSIZE,RCUTS)
 
38
C 
 
39
C  copy stuff from array A to array B
 
40
C 
 
41
C     .................
 
42
C     .................          ..............
 
43
C     ....aaaaaaa......          .bbbbbbb......
 
44
C     ....aaaaaaa......   ==>    .bbbbbbb......
 
45
C     ....aaaaaaa......          .bbbbbbb......
 
46
C     .................          ..............
 
47
C     .................          ..............
 
48
C     .................
 
49
C 
 
50
      IMPLICIT NONE
 
51
C 
 
52
      REAL    A(*),B(*),RCUTS(2)
 
53
C 
 
54
      INTEGER AOFF,BOFF,NOPIX,NOLINS,AXSIZE,BXSIZE
 
55
      INTEGER NY,NX,AAOFF,BBOFF,KOFF,LOFF,LINLIM
 
56
C 
 
57
      AAOFF = AOFF
 
58
      BBOFF = BOFF                         !we don't modify original offsets...
 
59
C 
 
60
      IF ((NOPIX+BOFF-1).GT.BXSIZE) THEN
 
61
         LINLIM = BXSIZE - BOFF + 1
 
62
      ELSE
 
63
         LINLIM = NOPIX
 
64
      ENDIF
 
65
C 
 
66
      DO 2000 NY=1,NOLINS
 
67
         KOFF = AAOFF
 
68
         LOFF = BBOFF
 
69
         DO 1000 NX=1,LINLIM
 
70
            B(LOFF) = A(KOFF)
 
71
            IF (B(LOFF).LT.RCUTS(1)) RCUTS(1) = B(LOFF)
 
72
            IF (B(LOFF).GT.RCUTS(2)) RCUTS(2) = B(LOFF)
 
73
            LOFF = LOFF + 1
 
74
            KOFF = KOFF + 1
 
75
1000     CONTINUE
 
76
         AAOFF = AAOFF + AXSIZE
 
77
         BBOFF = BBOFF + BXSIZE
 
78
2000  CONTINUE
 
79
C      
 
80
      RETURN
 
81
      END
 
82
 
 
83
      SUBROUTINE GXDOIT(A,B,C,THR,INOUT,SIZE,TRUSIZ)
 
84
C
 
85
      IMPLICIT NONE
 
86
C
 
87
      INTEGER    INOUT,SIZE,TRUSIZ
 
88
      INTEGER    N,NOUT
 
89
C
 
90
      REAL       A(SIZE),B(SIZE),C(*),THR(2)
 
91
C
 
92
      NOUT = 0
 
93
      IF(INOUT .EQ. 0) THEN             !use .GE. for backw. compatibility
 
94
         DO 1000, N=1,SIZE
 
95
            IF (B(N).GE.THR(1)) THEN
 
96
               NOUT = NOUT + 1
 
97
               C(NOUT) = A(N)
 
98
            ENDIF
 
99
1000     CONTINUE
 
100
C
 
101
      ELSE IF(INOUT .EQ. 1) THEN
 
102
         DO 2000, N=1,SIZE
 
103
            IF ((B(N).GE.THR(1)) .AND.
 
104
     +          (B(N).LE.THR(2))) THEN
 
105
               NOUT = NOUT + 1
 
106
               C(NOUT) = A(N)
 
107
            ENDIF
 
108
2000     CONTINUE
 
109
C 
 
110
      ELSE
 
111
         DO 5000, N=1,SIZE
 
112
            IF ((B(N).LT.THR(1)) .OR.
 
113
     +          (B(N).GT.THR(2))) THEN
 
114
               NOUT = NOUT + 1
 
115
               C(NOUT) = A(N)
 
116
            ENDIF
 
117
5000     CONTINUE
 
118
      ENDIF
 
119
C
 
120
      TRUSIZ = NOUT
 
121
      RETURN
 
122
      END
 
123
 
 
124
      SUBROUTINE GXIZAMAP(NPIXA,NPIXB,STARTB,STEPB,A,B,C,FLAG)
 
125
C
 
126
      IMPLICIT NONE
 
127
C
 
128
      REAL         A(*),B(*),C(*)
 
129
      REAL         RVAL,ENDB,DIFF,DIFF1,DIFF2
 
130
C
 
131
      DOUBLE PRECISION STARTB,STEPB
 
132
C
 
133
      INTEGER      NPIXA(*),NPIXB(*),FLAG(*)
 
134
      INTEGER      N,KIDX,TOTAL,IOFF,NN,MCLOSE
 
135
C
 
136
      TOTAL = NPIXA(1) * NPIXA(2) * NPIXA(3)
 
137
      IF (FLAG(2).EQ.1) GOTO 3000
 
138
      IF (FLAG(2).EQ.2) GOTO 5000
 
139
C
 
140
C -------------------------------------------------------C
 
141
C
 
142
C  here for equidistant map space
 
143
C
 
144
C -------------------------------------------------------C
 
145
C
 
146
      ENDB = STARTB + ((NPIXB(1) - 1) * STEPB)
 
147
      IF (FLAG(1).EQ.1) GOTO 1000
 
148
C
 
149
      DO 400 N=1,TOTAL
 
150
         RVAL = A(N)
 
151
C
 
152
         IF (RVAL.LE.STARTB) THEN
 
153
            KIDX = 1
 
154
         ELSE IF (RVAL.GE.ENDB) THEN
 
155
            KIDX = NPIXB(1)
 
156
         ELSE
 
157
            KIDX = NINT((RVAL-STARTB)/STEPB) + 1
 
158
         ENDIF
 
159
C
 
160
         C(N) = B(KIDX)
 
161
400   CONTINUE
 
162
C
 
163
      RETURN
 
164
C
 
165
C  here we leave all pixels with intensities outside the map coordinate space
 
166
C  as they are
 
167
C
 
168
1000  DO 1400 N=1,TOTAL
 
169
         RVAL = A(N)
 
170
C
 
171
         IF ( (RVAL.LT.STARTB) .OR. (RVAL.GT.ENDB) ) THEN
 
172
            C(N) = A(N)
 
173
         ELSE
 
174
            KIDX = NINT((RVAL-STARTB)/STEPB) + 1
 
175
            C(N) = B(KIDX)
 
176
         ENDIF
 
177
1400  CONTINUE
 
178
C
 
179
      RETURN
 
180
C
 
181
C -------------------------------------------------------C
 
182
C
 
183
C  here we have the coordinates stored in the 1. line of the map frame,
 
184
C  intensities follow in the 2. line
 
185
C
 
186
C -------------------------------------------------------C
 
187
C
 
188
3000  IOFF = NPIXB(1)
 
189
      STARTB = B(1)
 
190
      ENDB = B(IOFF)
 
191
      IF (FLAG(1).EQ.1) GOTO 4000
 
192
C
 
193
      DO 3500 N=1,TOTAL
 
194
         RVAL = A(N)
 
195
C
 
196
         IF (RVAL.LE.STARTB) THEN
 
197
            KIDX = 1 + IOFF
 
198
         ELSE IF (RVAL.GE.ENDB) THEN
 
199
            KIDX = IOFF + IOFF
 
200
         ELSE
 
201
            MCLOSE = 1
 
202
            DIFF1 = ABS(RVAL-STARTB)
 
203
C
 
204
            DO 3200 NN=2,IOFF                  !find closest coordinate
 
205
               DIFF = ABS(RVAL-B(NN))
 
206
               IF (DIFF.LT.DIFF1) THEN
 
207
                  DIFF1 = DIFF
 
208
                  MCLOSE = NN
 
209
               ELSE
 
210
                  KIDX = MCLOSE + IOFF
 
211
                  GOTO 3333
 
212
               ENDIF
 
213
3200        CONTINUE
 
214
         ENDIF
 
215
C
 
216
3333     C(N) = B(KIDX)
 
217
3500  CONTINUE
 
218
C
 
219
      RETURN
 
220
C
 
221
C  here we leave all pixels with intensities outside the map coordinate space
 
222
C  as they are
 
223
C
 
224
4000  DO 4500 N=1,TOTAL
 
225
         RVAL = A(N)
 
226
C
 
227
         IF ( (RVAL.LT.STARTB) .OR. (RVAL.GT.ENDB) ) THEN
 
228
            C(N) = A(N)
 
229
C
 
230
         ELSE
 
231
            MCLOSE = 1
 
232
            DIFF1 = ABS(RVAL-STARTB)
 
233
C
 
234
            DO 4200 NN=2,IOFF                  !find closest coordinate
 
235
               DIFF = ABS(RVAL-B(NN))
 
236
               IF (DIFF.LT.DIFF1) THEN
 
237
                  DIFF1 = DIFF
 
238
                  MCLOSE = NN
 
239
               ELSE
 
240
                  KIDX = MCLOSE + IOFF
 
241
                  GOTO 4333
 
242
               ENDIF
 
243
4200        CONTINUE
 
244
4333        C(N) = B(KIDX)
 
245
         ENDIF
 
246
C
 
247
4500  CONTINUE
 
248
C
 
249
      RETURN
 
250
C
 
251
C -------------------------------------------------------C
 
252
C
 
253
C  here we have the intervals stored in the 1. line of the map frame,
 
254
C  intensities follow in the 2. line (for start,end - in between interpolation ...)
 
255
C
 
256
C -------------------------------------------------------C
 
257
C
 
258
5000  IOFF = NPIXB(1)
 
259
      STARTB = B(1)
 
260
      ENDB = B(IOFF)
 
261
      IF (FLAG(1).EQ.1) GOTO 6000
 
262
C
 
263
      DO 5555 N=1,TOTAL
 
264
         RVAL = A(N)
 
265
C
 
266
         IF (RVAL.LT.STARTB) THEN
 
267
            RVAL = STARTB
 
268
         ELSE IF (RVAL.GT.ENDB) THEN
 
269
            RVAL = ENDB
 
270
         ENDIF
 
271
C
 
272
         DO 5200 NN=1,IOFF,2                 !find enclosing interval
 
273
            IF (RVAL.LE.B(NN+1)) THEN
 
274
                KIDX = NN
 
275
                GOTO 5333
 
276
            ELSE IF (RVAL.LT.B(NN+2)) THEN
 
277
                C(N) = A(N)
 
278
                GOTO 5555                    !no interval there ...!
 
279
            ENDIF
 
280
5200     CONTINUE
 
281
C
 
282
5333     DIFF = B(KIDX+1) - B(KIDX)
 
283
         DIFF1 = (RVAL - B(KIDX)) / DIFF
 
284
         DIFF2 = (B(KIDX+1) - RVAL) / DIFF
 
285
         KIDX = KIDX + IOFF                  !point to data line
 
286
         C(N) = ( B(KIDX+1) * DIFF1 ) + ( B(KIDX) * DIFF2 )
 
287
5555  CONTINUE
 
288
C
 
289
      RETURN
 
290
C
 
291
C  here we leave all pixels with intensities outside the map coordinate space
 
292
C  as they are
 
293
C
 
294
6000  DO 6555 N=1,TOTAL
 
295
         RVAL = A(N)
 
296
C
 
297
         IF (RVAL.LT.STARTB) THEN
 
298
            C(N) = A(N)
 
299
            GOTO 6555
 
300
         ELSE IF (RVAL.GT.ENDB) THEN
 
301
            C(N) = A(N)
 
302
            GOTO 6555
 
303
         ENDIF
 
304
C
 
305
         DO 6200 NN=1,IOFF,2                  !find enclosing interval
 
306
            IF (RVAL.LE.B(NN+1)) THEN
 
307
               KIDX = NN
 
308
               GOTO 6333
 
309
            ELSE IF (RVAL.LT.B(NN+2)) THEN
 
310
               C(N) = A(N)
 
311
               GOTO 6555                      !no interval there ...!
 
312
            ENDIF
 
313
6200     CONTINUE
 
314
C
 
315
6333     DIFF = B(KIDX+1) - B(KIDX)
 
316
         DIFF1 = (RVAL - B(KIDX)) / DIFF
 
317
         DIFF2 = (B(KIDX+1) - RVAL) / DIFF
 
318
         KIDX = KIDX + IOFF                   !point to data line
 
319
         C(N) = ( B(KIDX+1) * DIFF1 ) + ( B(KIDX) * DIFF2 )
 
320
6555  CONTINUE
 
321
C
 
322
      RETURN
 
323
      END
 
324
 
 
325
      SUBROUTINE GXMATRX(A,NPIX,C)
 
326
C
 
327
C +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 
328
C
 
329
C  straight forward transposition of an image (seen as a matrix)
 
330
C  this however implies, that the first element is in the upper left corner
 
331
C  since matrices are usually written like that ...
 
332
C
 
333
C -------------------------------------------------------------------
 
334
C
 
335
      IMPLICIT NONE
 
336
C
 
337
      INTEGER      NPIX(*)
 
338
      INTEGER      NX,NY
 
339
      INTEGER      XDIM,YDIM,OFFA,OFFC,WORK
 
340
C
 
341
      REAL         A(*),C(*)
 
342
C
 
343
C  init
 
344
      XDIM = NPIX(1)
 
345
      YDIM = NPIX(2)
 
346
      OFFC = 0
 
347
      WORK = XDIM * (YDIM-1)
 
348
C
 
349
C  loop through rows in output matrix
 
350
      DO 1000 NX=XDIM,1,-1
 
351
         OFFA = WORK
 
352
C
 
353
         DO 400 NY=1,YDIM                          !loop till middle of row
 
354
            C(NY+OFFC) = A(NX+OFFA)
 
355
            OFFA = OFFA - XDIM
 
356
400      CONTINUE
 
357
C
 
358
         OFFC = OFFC + YDIM
 
359
1000  CONTINUE
 
360
      RETURN
 
361
C
 
362
      END