← Back to branch summary

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

« back to all changes in this revision

Viewing changes to applic/display/src/perspec.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:
applic/display/src/perspec.for
 
1
C @(#)perspec.for       13.1.1.1 (ES0-DMD) 06/02/98 18:20:55
 
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
      PROGRAM PERSP
 
30
C
 
31
C++++++++++++++++++++++++++++++++++++++++++++++++++
 
32
C
 
33
C.IDENTIFICATION
 
34
C  program PERSP                        version 1.00    920103
 
35
C  K. Banse                             ESO - Garching
 
36
C  1.10  930512         2.00    940104
 
37
C
 
38
C.KEYWORDS
 
39
C  perspective, rebinning
 
40
C
 
41
C.INPUT/OUTPUT
 
42
C  the following keywords are used:
 
43
C
 
44
C  IN_A/C/1/60                  input frame
 
45
C  OUT_A/C/1/60                 output frame
 
46
C  INPUTR/1/2                   rotation angle in degrees, in [0,360]
 
47
C                               tilt angle in degrees in [0,90]
 
48
C
 
49
C.VERSIONS
 
50
C  1.00   creation
 
51
C  1.10   STKWRR had 1 param. too many !
 
52
C  2.00   add z-perspective plot
 
53
C       
 
54
C--------------------------------------------------
 
55
C
 
56
      IMPLICIT NONE
 
57
C
 
58
      INTEGER      IAV,N,NAXISA,NAXISC
 
59
      INTEGER*8    PNTRA,PNTRB,PNTRC
 
60
      INTEGER      IMNOA,IMNOB,IMNOC,STAT
 
61
      INTEGER      NPIXA(3),NPIXB(2),NPIXC(2),KZ,KNO
 
62
      INTEGER      N1,N2,NSIZ,NPERSP,DIM3C,KDIFF,PLANES(20)
 
63
      INTEGER      UNI(1),NULO,MADRID(1)
 
64
      INTEGER      ISCAL
 
65
C
 
66
      CHARACTER*60 FRAMEA,FRAMEB
 
67
      CHARACTER    CUNITA*64,IDENTA*72,STR*80
 
68
C
 
69
      DOUBLE PRECISION STEPA(3),STARTA(3),STARTB(2),STEPB(2),STEPC(2)
 
70
      DOUBLE PRECISION ST(2)
 
71
C
 
72
      REAL         RNULL
 
73
      REAL         ANGDEG,ANGLE,AC,AS,X1,X2,X3,Y1,Y2,Y3,X4,Y4
 
74
      REAL         PC, PCDEG
 
75
      REAL         PIXROT(2),ROTOFF(2),EN(2),RST(2)
 
76
      REAL         QUOT,RQUOT,R1,R2,RROT(4),CUTS(4)
 
77
      REAL         FACTO,RBUF(20),ZCUTS(2)
 
78
C
 
79
      COMMON   /VMR/  MADRID
 
80
C
 
81
      INCLUDE  'MID_INCLUDE:ST_DEF.INC'
 
82
C 
 
83
      EQUIVALENCE (X2,RBUF(1)),(Y2,RBUF(2)),(X4,RBUF(3)),(Y4,RBUF(4))
 
84
      EQUIVALENCE (X3,RBUF(5)),(Y3,RBUF(6)),(X1,RBUF(7)),(Y1,RBUF(8))
 
85
C
 
86
      DATA     FACTO  /0.0174533/               !  Pi / 180.
 
87
      DATA     NPIXA  /1,1,1/
 
88
      DATA     IDENTA /' '/, CUNITA /' '/       !necessary ... !!!
 
89
C
 
90
      INCLUDE  'MID_INCLUDE:ST_DAT.INC'
 
91
C 
 
92
C  get into MIDAS
 
93
      CALL STSPRO('PERSP')
 
94
C
 
95
C  get input frame + map it
 
96
      CALL STKRDC('IN_A',1,1,60,IAV,FRAMEA,UNI,NULO,STAT)
 
97
      CALL STIGET(FRAMEA,D_R4_FORMAT,F_I_MODE,F_IMA_TYPE,3,
 
98
     +            NAXISA,NPIXA,STARTA,STEPA,
 
99
     +            IDENTA,CUNITA,PNTRA,IMNOA,STAT)
 
100
      IF (NAXISA.EQ.1) 
 
101
     +   CALL STETER(1,'We need a 2-dim image ...')
 
102
      IF ((NPIXA(1).EQ.1) .OR. (NPIXA(2).EQ.1))
 
103
     +   CALL STETER(1,'We need a 2-dim image ...')
 
104
      CALL STDRDR(IMNOA,'LHCUTS',1,4,IAV,CUTS,UNI,NULO,STAT)
 
105
      ZCUTS(1) = CUTS(3)
 
106
      ZCUTS(2) = CUTS(4)
 
107
      RNULL = CUTS(3) - 0.01*(CUTS(4)-CUTS(3))                !null val < min
 
108
C 
 
109
      R1 = NPIXA(1) - 1.
 
110
      R2 = NPIXA(2) - 1.
 
111
      KDIFF = 12
 
112
      DIM3C = 1
 
113
C
 
114
C  get name of output frame + plane no.s (or action code)
 
115
      CALL STKRDC('OUT_A',1,1,60,IAV,FRAMEB,UNI,NULO,STAT)
 
116
      CALL STKRDC('INPUTC',1,1,80,IAV,STR,UNI,NULO,STAT)
 
117
C
 
118
C  get rotation angle + perspective angle + (scaling size)
 
119
      CALL STKRDR('INPUTR',1,4,IAV,RROT,UNI,NULO,STAT)
 
120
      IF (STR(1:4).EQ.'zper') THEN
 
121
         PLANES(1) = 1
 
122
         PLANES(2) = 0
 
123
         KNO = 0
 
124
         ISCAL = NINT(RROT(3))
 
125
         IF (ISCAL.LT.1) 
 
126
     +      CALL STETER(1,'Invalid scaling size (should be > 0)...')
 
127
         ZCUTS(2) = RROT(4)
 
128
         IF (ZCUTS(2).LE.ZCUTS(1)) CALL STETER(1,
 
129
     +      'Invalid scale limit (should be > LHCUTS(3))...')
 
130
C 
 
131
C  get planes to work on
 
132
      ELSE IF ((STR(1:1).EQ.'A') .OR. (STR(1:1).EQ.'a')) THEN
 
133
         PLANES(1) = -1                    !indicate all planes
 
134
         KNO = NPIXA(3) - 1                !no. of planes - 1
 
135
      ELSE
 
136
         CALL GENCNV(STR,1,20,PLANES,QUOT,STEPC,IAV)
 
137
         IF (IAV.LE.0)
 
138
     +      CALL STETER(5,'Invalid planes specified ...')
 
139
         IF ((PLANES(1).LE.0).OR.(PLANES(1).GT.NPIXA(3)))
 
140
     +      CALL STETER(6,'specified plane outside limits... ')
 
141
         N = PLANES(1) - 1                 !update pointer
 
142
         KNO = IAV - 1                     !no. of planes - 1
 
143
         PNTRA = PNTRA + N*NPIXA(1)*NPIXA(2)
 
144
         IAV = IAV + 1
 
145
         PLANES(IAV) = 0                   !mark the end
 
146
      ENDIF
 
147
C 
 
148
      ANGDEG = RROT(1)
 
149
      PCDEG = RROT(2)
 
150
C 
 
151
      ANGLE = ANGDEG * FACTO
 
152
      N = (NPIXA(1) + 1) / 2                  !use center pixel
 
153
      PIXROT(1) = N - 1
 
154
      N = (NPIXA(2) + 1) / 2                  !use center pixel
 
155
      PIXROT(2) = N - 1
 
156
C 
 
157
      AC = COS(ANGLE)
 
158
      AS = SIN(ANGLE)
 
159
      IF ((ANGDEG.LT.0.1).OR.(ANGDEG.GT.359.9)) THEN
 
160
         NPIXB(1) = NPIXA(1)                   !no rotation
 
161
         NPIXB(2) = NPIXA(2)
 
162
         STARTB(1) = STARTA(1)
 
163
         STARTB(2) = STARTA(2)
 
164
         STEPB(1) = STEPA(1)
 
165
         STEPB(2) = STEPA(2)
 
166
         N1 = 0
 
167
C 
 
168
         RBUF(1) = 0.
 
169
         RBUF(2) = R2
 
170
         RBUF(3) = 0.
 
171
         RBUF(4) = 0.
 
172
         RBUF(5) = R1
 
173
         RBUF(6) = 0.
 
174
         RBUF(7) = R1
 
175
         RBUF(8) = R2
 
176
C 
 
177
         GOTO 3000
 
178
      ELSE
 
179
         N1 = 1                                !yes. we do
 
180
      ENDIF
 
181
C
 
182
C  to get size of output frame, rotate corners around lower left corner
 
183
C  (assumed = 0,0)
 
184
C 
 
185
      CALL ROTA1(R1,R2,AC,AS,X1,Y1)         !(NX-1,NY-1)
 
186
      CALL ROTA1(0.,R2,AC,AS,X2,Y2)         !(0,NY-1)
 
187
      CALL ROTA1(R1,0.,AC,AS,X3,Y3)         !(NX-1,0)
 
188
      X4 = 0.
 
189
      Y4 = 0.                               !rotation point stays fixed
 
190
C
 
191
C  get result frame "corners"
 
192
      ST(1) = MIN(0.,X1,X2,X3)             !rotated x-start
 
193
      EN(1) = MAX(0.,X1,X2,X3)             !rotated x-end
 
194
      ST(2) = MIN(0.,Y1,Y2,Y3)             !rotated y-start
 
195
      EN(2) = MAX(0.,Y1,Y2,Y3)             !rotated y-end
 
196
C
 
197
C  from coords. of rotation point determine shift after rotation
 
198
      QUOT = 1. - AC
 
199
      ROTOFF(1) = PIXROT(1)*QUOT + PIXROT(2)*AS
 
200
      ROTOFF(2) = PIXROT(2)*QUOT - PIXROT(1)*AS
 
201
      CALL ROTA1(PIXROT(1),PIXROT(2),AC,AS,RROT(1),RROT(2))
 
202
C
 
203
C  choose corners such, that rotation point is hit ...
 
204
      DO 2000, N=1,2
 
205
         STEPB(N) = STEPA(N)
 
206
         QUOT =  RROT(N) - ST(N)           !distance of fixed p. to rotation p.
 
207
         RQUOT = FLOAT(INT(QUOT)) + 1.     !is forced to multiples of steps
 
208
         ST(N) = RROT(N) - RQUOT
 
209
         QUOT = EN(N) - ST(N)              !distance of start to end point
 
210
         RQUOT = FLOAT(INT(QUOT)) + 1.     !in multiples of steps
 
211
         EN(N) = ST(N) + RQUOT
 
212
         NPIXB(N) = NINT(RQUOT) + 1        !no. of points in new image
 
213
2000  CONTINUE
 
214
C
 
215
C  now update start of rotated image, so that rotation point stays fixed
 
216
      RST(1) = ST(1) + ROTOFF(1)
 
217
      RST(2) = ST(2) + ROTOFF(2)
 
218
      STARTB(1) = STARTA(1) + ( RST(1)*STEPA(1) )      !we work in 0,0 space
 
219
      STARTB(2) = STARTA(2) + ( RST(2)*STEPA(2) )
 
220
      DO 2200, N=1,7,2                       
 
221
         RBUF(N) = RBUF(N) - ST(1)              !relate to new frame
 
222
         RBUF(N+1) = RBUF(N+1) - ST(2)
 
223
2200  CONTINUE
 
224
C
 
225
C  create and map intermediate frame
 
226
C 
 
227
3000  NSIZ = NPIXB(1) * NPIXB(2)
 
228
C 
 
229
C  do the rotation
 
230
      IF (N1.EQ.1) THEN
 
231
         IF (DIM3C.EQ.1) THEN
 
232
            CALL STFCRE('midpersp',D_R4_FORMAT,F_X_MODE,1,NSIZ,
 
233
     +                  IMNOB,STAT)
 
234
            CALL STFMAP(IMNOB,F_X_MODE,1,NSIZ,IAV,PNTRB,STAT)
 
235
         ENDIF
 
236
         CALL ROTA2(MADRID(PNTRA),MADRID(PNTRB),NPIXA,NPIXB,
 
237
     +              ST,AC,AS,RNULL)
 
238
      ENDIF
 
239
C 
 
240
C  create result frame
 
241
      IF (DIM3C.EQ.1) THEN
 
242
         NPIXC(1) = NPIXB(1)
 
243
         STEPC(1) = STEPB(1)
 
244
         IF (PCDEG.LE.0.1) THEN                 !no tilt at all ...
 
245
            N2 = 0
 
246
            PC = 1
 
247
            NPIXC(2) = NPIXB(2)
 
248
         ELSE IF (PCDEG.GT.89.5) THEN           !total side view => single line
 
249
            N2 = 2
 
250
            NPIXC(2) = 1
 
251
            PC = 0
 
252
         ELSE
 
253
            N2 = 1
 
254
            PC = COS(PCDEG * FACTO)
 
255
            NPIXC(2) = INT (NPIXB(2) * PC)
 
256
            NPERSP = NPIXC(2)
 
257
            STEPC(2) = STEPB(2) * PC
 
258
         ENDIF
 
259
C
 
260
C  move corners to world coords.
 
261
         DO 4100, N=1,7,2
 
262
            RBUF(N) = STARTB(1) + (RBUF(N)*STEPC(1))
 
263
            IAV = INT (RBUF(N+1)*PC)
 
264
            IF (IAV.GE.NPIXC(2)) IAV = NPIXC(2) - 1
 
265
            RBUF(N+1) = STARTB(2) + (IAV*STEPC(2))
 
266
4100     CONTINUE
 
267
         IF (KNO.GT.0) THEN                      !if more than one plane
 
268
            QUOT = STEPC(2) * ((NPERSP+KDIFF)*KNO)
 
269
            RBUF(9) = RBUF(2) + QUOT
 
270
            RBUF(10) = RBUF(4) + QUOT
 
271
            RBUF(11) = RBUF(6) + QUOT
 
272
            RBUF(12) = RBUF(8) + QUOT
 
273
            RBUF(13) = KNO + 1                   !no. of planes
 
274
            RBUF(14) = KDIFF
 
275
            RBUF(15) = STEPC(2) * KDIFF
 
276
            RBUF(16) = STEPC(2) * (NPERSP+KDIFF)
 
277
C 
 
278
            NAXISC = 2
 
279
            NPIXC(2) = (KNO+1)*NPIXC(2) + KNO*KDIFF
 
280
            RQUOT = STARTB(2) + (NPIXC(2)-1)*STEPC(2)         !y-end
 
281
            IF (STEPC(2).GT.0.) THEN
 
282
               IF (RBUF(9).GT.RQUOT) RBUF(9) = RQUOT
 
283
               IF (RBUF(10).GT.RQUOT) RBUF(10) = RQUOT
 
284
               IF (RBUF(11).GT.RQUOT) RBUF(11) = RQUOT
 
285
               IF (RBUF(12).GT.RQUOT) RBUF(12) = RQUOT
 
286
            ELSE
 
287
               IF (RBUF(9).LT.RQUOT) RBUF(9) = RQUOT
 
288
               IF (RBUF(10).LT.RQUOT) RBUF(10) = RQUOT
 
289
               IF (RBUF(11).LT.RQUOT) RBUF(11) = RQUOT
 
290
               IF (RBUF(12).LT.RQUOT) RBUF(12) = RQUOT
 
291
            ENDIF
 
292
         ELSE
 
293
            NAXISC = NAXISA
 
294
         ENDIF
 
295
C 
 
296
         CALL STKWRR('OUTPUTR',RBUF,1,20,UNI,STAT)
 
297
         IF (STR(1:4).EQ.'zper') THEN
 
298
            IDENTA(1:) = 'perspective z-plot of input frame '
 
299
            NAXISC = 2
 
300
            NPIXC(2) = NPIXC(2)+ISCAL
 
301
            NPERSP = NPIXC(2)                     !in case we reset it later on
 
302
         ELSE
 
303
            IDENTA(1:) = 
 
304
     +      'perspective plot of selected planes of input frame '
 
305
         ENDIF
 
306
         CALL STIPUT(FRAMEB,D_R4_FORMAT,F_O_MODE,F_IMA_TYPE,
 
307
     +               NAXISC,NPIXC,STARTB,STEPC,
 
308
     +               IDENTA,CUNITA,PNTRC,IMNOC,STAT)
 
309
      ENDIF
 
310
C
 
311
C  and do the perspective rebinning
 
312
      IF (N2.EQ.1) THEN                         !real job
 
313
         NPIXC(2) = NPERSP                      !reset to original value (3-dim)
 
314
         IF (N1.EQ.1) THEN
 
315
            IF (STR(1:4).EQ.'zper') THEN
 
316
               CALL JPERSP(MADRID(PNTRB),MADRID(PNTRC),
 
317
     +                     NPIXB,NPIXC,PC,RNULL,ZCUTS,ISCAL)
 
318
               CALL STDWRI(IMNOC,'NPIX',NPIXC(2),2,1,UNI,STAT)
 
319
            ELSE
 
320
               CALL KPERSP(MADRID(PNTRB),MADRID(PNTRC),
 
321
     +                     NPIXB,NPIXC,PC)
 
322
            ENDIF
 
323
         ELSE
 
324
            IF (STR(1:4).EQ.'zper') THEN
 
325
               CALL JPERSP(MADRID(PNTRA),MADRID(PNTRC),
 
326
     +                     NPIXB,NPIXC,PC,RNULL,ZCUTS,ISCAL)
 
327
               CALL STDWRI(IMNOC,'NPIX',NPIXC(2),2,1,UNI,STAT)
 
328
            ELSE
 
329
               CALL KPERSP(MADRID(PNTRA),MADRID(PNTRC),
 
330
     +                     NPIXB,NPIXC,PC)
 
331
            ENDIF
 
332
         ENDIF
 
333
      ELSE IF (N2.EQ.0) THEN                    !just rotated image
 
334
         IF (N1.EQ.1) THEN
 
335
            CALL COPYF(MADRID(PNTRB),MADRID(PNTRC),NSIZ)
 
336
         ELSE
 
337
            CALL COPYF(MADRID(PNTRA),MADRID(PNTRC),NSIZ)
 
338
         ENDIF
 
339
      ELSE                                      !only sideline
 
340
         IF (N1.EQ.1) THEN
 
341
            CALL COPYF(MADRID(PNTRB),MADRID(PNTRC),NPIXB(1))
 
342
         ELSE
 
343
            CALL COPYF(MADRID(PNTRA),MADRID(PNTRC),NPIXB(1))
 
344
         ENDIF
 
345
      ENDIF
 
346
C 
 
347
C  check for 3-dim case
 
348
      DIM3C = DIM3C + 1
 
349
      IF (PLANES(1).EQ.-1) THEN
 
350
         IF (DIM3C.GT.NPIXA(3)) THEN
 
351
            GOTO 9000
 
352
         ELSE
 
353
            KZ = 1
 
354
         ENDIF
 
355
      ELSE
 
356
         IF (PLANES(DIM3C).EQ.0) THEN
 
357
            GOTO 9000
 
358
         ELSE
 
359
            IF ((PLANES(DIM3C).LE.0).OR.(PLANES(DIM3C).GT.NPIXA(3)))
 
360
     +         CALL STETER(6,'specified plane outside limits... ')
 
361
            KZ = PLANES(DIM3C) - PLANES(DIM3C-1)
 
362
         ENDIF
 
363
      ENDIF
 
364
      PNTRA = PNTRA + (NPIXA(1)*NPIXA(2))*KZ
 
365
      PNTRC = PNTRC + (NPIXC(1) * (NPIXC(2) + KDIFF))
 
366
      GOTO 3000
 
367
C
 
368
9000  CALL STDWRR(IMNOC,'LHCUTS',CUTS,1,4,UNI,STAT)
 
369
      CALL STSEPI
 
370
      END
 
371
 
 
372
      SUBROUTINE ROTA1(X,Y,AC,AS,XOUT,YOUT)
 
373
C
 
374
      IMPLICIT NONE
 
375
C
 
376
      REAL     X,Y,AC,AS,XOUT,YOUT
 
377
C
 
378
      XOUT = X*AC - Y*AS
 
379
      YOUT = X*AS + Y*AC
 
380
C
 
381
      RETURN
 
382
      END
 
383
 
 
384
      SUBROUTINE KPERSP(A,B,NPIXA,NPIXB,AC)
 
385
C
 
386
C  K. Banse            920103
 
387
C
 
388
C  do "backward" perspective
 
389
C 
 
390
C  XPERS = X
 
391
C  YPERS = Y/AC
 
392
C
 
393
      IMPLICIT NONE
 
394
C
 
395
      INTEGER      INDXB,KSTX,KY
 
396
      INTEGER      NO1,NO3
 
397
      INTEGER      NX,NY,SIZEA,XDIMA,YDIMA
 
398
      INTEGER      NPIXA(2),NPIXB(2)
 
399
C
 
400
      REAL         A(*),B(*)
 
401
      REAL         AC,YFRAC
 
402
C 
 
403
      DOUBLE PRECISION  RYAC,ACINV
 
404
C
 
405
C  init
 
406
      XDIMA = NPIXA(1)
 
407
      YDIMA = NPIXA(2)
 
408
      SIZEA = XDIMA * YDIMA
 
409
C
 
410
      ACINV = 1.D0 / AC
 
411
      INDXB = 0
 
412
      RYAC = ACINV
 
413
C
 
414
C  loop through result lines
 
415
      DO 3000, NY=1,NPIXB(2)
 
416
         DO 1000, NX=1,NPIXB(1)
 
417
            INDXB = INDXB + 1              !loop along result line
 
418
            KY = INT(RYAC)
 
419
            KSTX = XDIMA*KY                !get 1. pixel in this line - 1
 
420
            NO1 = NX + KSTX                !get pixel close to integer part of
 
421
            NO3 = NO1 + XDIMA               !get low pixel in next line up
 
422
            IF (NO3.GT.SIZEA) THEN          !are we out of grid in y?
 
423
               B(INDXB) = A(NO1)
 
424
            ELSE
 
425
               YFRAC = RYAC - KY
 
426
               B(INDXB) = A(NO1) + (YFRAC * (A(NO3) - A(NO1)))
 
427
            ENDIF
 
428
1000     CONTINUE
 
429
         RYAC = RYAC + ACINV
 
430
3000  CONTINUE
 
431
C 
 
432
      RETURN
 
433
      END
 
434
 
 
435
      SUBROUTINE JPERSP(A,B,NPIXA,NPIXB,AC,ZNULL,FMINMA,ISCAL)
 
436
C
 
437
C
 
438
C  do "backward" z-perspective
 
439
C
 
440
C  XPERS = X
 
441
C  YPERS = Y/AC 
 
442
C 
 
443
C  Z = YPERS * ISCAL  =>      XPERS,YPERS+n   = Z
 
444
C                             XPERS,YPERS+n-1 = ...
 
445
C                                  ...
 
446
C                             XPERS,YPERS     = Min
 
447
C
 
448
      IMPLICIT NONE
 
449
C
 
450
      INTEGER      INDXB,KSTX,KY,ISCAL
 
451
      INTEGER      NO1,NO3,N
 
452
      INTEGER      NX,NY,SIZEA,XDIMA,YDIMA,XDIMB
 
453
      INTEGER      NPIXA(2),NPIXB(2)
 
454
      INTEGER      KSCAL,IZ,KNDXB,KK
 
455
C
 
456
      REAL         A(*),B(*)
 
457
      REAL         AC,FMINMA(2),ZNULL
 
458
      REAL         YFRAC,RMIN,RMAX
 
459
      REAL         Z,DZ,SCALF,CMIN,FF
 
460
      REAL         VALU
 
461
C
 
462
      DOUBLE PRECISION  RYAC,ACINV
 
463
C
 
464
C  init
 
465
      XDIMA = NPIXA(1)
 
466
      YDIMA = NPIXA(2)
 
467
      SIZEA = XDIMA * YDIMA
 
468
      XDIMB = NPIXB(1)
 
469
C
 
470
      RMIN = FMINMA(1)
 
471
      RMAX = FMINMA(2)
 
472
      ACINV = 1.D0 / AC
 
473
      INDXB = NPIXB(1)*NPIXB(2)
 
474
      KSCAL = ISCAL - 1
 
475
      SCALF = KSCAL / (RMAX-RMIN)
 
476
      CMIN = RMIN*ISCAL
 
477
      KK = NPIXB(2)-ISCAL               !reserved space for spikes of last line
 
478
C
 
479
C  loop through result lines
 
480
      DO 300, NY=NPIXB(2),KK+1,-1
 
481
         DO 100, NX=1,XDIMB
 
482
            B(INDXB) = ZNULL
 
483
            INDXB = INDXB - 1              !loop along result line
 
484
100      CONTINUE
 
485
300   CONTINUE
 
486
C 
 
487
      RYAC = KK * ACINV
 
488
      DO 3000, NY=KK,1,-1
 
489
C
 
490
         DO 1000, NX=1,XDIMB
 
491
            KY = INT(RYAC)
 
492
            KSTX = XDIMA*KY                !get 1. pixel in this line - 1
 
493
            NO1 = NX + KSTX                !get pixel close to integer part of
 
494
            NO3 = NO1 + XDIMA              !get low pixel in next line up
 
495
            IF (NO3.GT.SIZEA) THEN         !are we out of grid in y?
 
496
               Z = A(NO1)
 
497
            ELSE
 
498
               YFRAC = RYAC - KY
 
499
               Z = A(NO1) + (YFRAC * (A(NO3) - A(NO1)))
 
500
            ENDIF
 
501
            IF (Z .GT. RMAX) Z = RMAX
 
502
C 
 
503
C  get height from intensity
 
504
            DZ = Z - RMIN
 
505
            IZ = NINT(DZ*SCALF)
 
506
            IF (IZ .GE. 1) THEN
 
507
               B(INDXB) = RMIN
 
508
               FF = DZ/IZ
 
509
               KNDXB = INDXB + XDIMB
 
510
               VALU = FF + RMIN
 
511
               DO 800, N=1,IZ
 
512
                  B(KNDXB) = VALU
 
513
                  VALU = VALU + FF
 
514
                  KNDXB = KNDXB + XDIMB
 
515
800            CONTINUE
 
516
            ELSE
 
517
               B(INDXB) = Z
 
518
            ENDIF
 
519
C
 
520
            INDXB = INDXB - 1              !loop along result line
 
521
1000     CONTINUE
 
522
C
 
523
         RYAC = RYAC - ACINV
 
524
3000  CONTINUE
 
525
C
 
526
C  get rid of superfluous lines in the end
 
527
      INDXB = NPIXB(1)*NPIXB(2)
 
528
      DO 4000, NY=NPIXB(2),KK+1,-1
 
529
         DO 3300, NX=1,XDIMB
 
530
            IF (B(INDXB) .GT. ZNULL) THEN
 
531
               NPIXB(2) = NY
 
532
               RETURN
 
533
            ENDIF
 
534
            INDXB = INDXB - 1
 
535
3300     CONTINUE
 
536
4000  CONTINUE
 
537
C 
 
538
      RETURN
 
539
      END
 
540
 
 
541
      SUBROUTINE ROTA2(A,B,NPIXA,NPIXB,START,AC,AS,ZNULL)
 
542
C
 
543
C  K. Banse            861029, 880902, 910204
 
544
C
 
545
C  do "backward" rotation, i.e. solve the system XROT = X*AC - Y*AS
 
546
C                                       YROT = X*AS + Y*AC for X,Y
 
547
C  to yield X = XROT*AC + YROT*AS
 
548
C          Y = YROT*AC - XROT*AS  which takes AC**2 + AS**2 = 1 into account....
 
549
C
 
550
      IMPLICIT NONE
 
551
C
 
552
      INTEGER      INDXB,KX,KY,KSTX
 
553
      INTEGER      NO1,NO2,NO3,NO4
 
554
      INTEGER      NX,NY,SIZEA,XDIMA,YDIMA
 
555
      INTEGER      NPIXA(2),NPIXB(2)
 
556
C
 
557
      REAL         A(*),B(*)
 
558
      REAL         AC,AS,ZNULL
 
559
      REAL         XT,YT,XFRAC,YFRAC,XLAST,YLAST
 
560
      REAL         XF,RX
 
561
C
 
562
      DOUBLE PRECISION  START(2)
 
563
      DOUBLE PRECISION  RYAC,RYAS
 
564
C
 
565
C  init
 
566
      XF = START(1)
 
567
      XDIMA = NPIXA(1)
 
568
      YDIMA = NPIXA(2)
 
569
      XLAST = XDIMA - 1.
 
570
      YLAST = YDIMA - 1.
 
571
      SIZEA = XDIMA * YDIMA
 
572
      INDXB = 0
 
573
C
 
574
      RYAS = START(2) * AS                 !init RYAS,RYAC according to start-y
 
575
      RYAC = START(2) * AC
 
576
C 
 
577
C  loop through result lines
 
578
      DO 3000, NY=1,NPIXB(2)
 
579
         RX = XF                           !init RX to start-x
 
580
C
 
581
         DO 1000, NX=1,NPIXB(1)
 
582
            INDXB = INDXB + 1              !loop along result line
 
583
            XT = RX*AC + RYAS
 
584
            IF ((XT.LT.0.).OR.(XT.GT.XLAST)) THEN
 
585
               B(INDXB) = ZNULL
 
586
               GOTO 1000
 
587
            ENDIF
 
588
            YT = RYAC - RX*AS
 
589
            IF ((YT.LT.0.).OR.(YT.GT.YLAST)) THEN
 
590
               B(INDXB) = ZNULL
 
591
               GOTO 1000
 
592
            ENDIF
 
593
C
 
594
C  interpolate...
 
595
            XT = XT + 1.                  !adjust 0,0 to 1,1 ...
 
596
            YT = YT + 1.
 
597
            KX = INT(XT)
 
598
            KY = INT(YT)
 
599
            KSTX = XDIMA*(KY-1)          !get 1. pixel in this line - 1
 
600
            NO1 = KX + KSTX          !get pixel close to int. part of real pixel
 
601
C
 
602
            NO2 = NO1 + 1                  !get next pixel in x-direction
 
603
CC            IF ((NO2-KSTX).GT.XDIMA) THEN  !are we out of grid in x?
 
604
            IF (KX.GE.XDIMA) THEN          !are we out of grid in x?
 
605
               IF (NO2.GT.SIZEA) THEN      !also out in y?
 
606
                  B(INDXB) = A(NO1)        !yes.
 
607
               ELSE
 
608
                  NO3 = NO1 + XDIMA        !no.
 
609
                  YFRAC = YT - KY
 
610
                  B(INDXB) = A(NO1) + ( YFRAC * (A(NO3) - A(NO1)) )
 
611
               ENDIF
 
612
            ELSE
 
613
               XFRAC = XT - KX
 
614
               NO3 = NO1 + XDIMA               !get low pixel in next line up
 
615
               IF (NO3.GT.SIZEA) THEN          !are we out of grid in y?
 
616
                  B(INDXB) = A(NO1) + ( XFRAC * (A(NO2) - A(NO1)) )
 
617
               ELSE
 
618
                  NO4 = NO3 + 1
 
619
                  YFRAC = YT - KY
 
620
                  B(INDXB) = A(NO1) +
 
621
     +              XFRAC * (A(NO2) - A(NO1)) +
 
622
     +              YFRAC * (A(NO3) - A(NO1)) +
 
623
     +              XFRAC * YFRAC * (A(NO1) - A(NO2) - A(NO3) + A(NO4))
 
624
               ENDIF
 
625
            ENDIF
 
626
C
 
627
1000     RX = RX + 1.0                        !move to next x-coord
 
628
C
 
629
         RYAS = RYAS + AS                     !move to next y-coord.
 
630
         RYAC = RYAC + AC       
 
631
3000  CONTINUE
 
632
C 
 
633
      RETURN
 
634
      END