← Back to branch summary

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

« back to all changes in this revision

Viewing changes to prim/display/src/avesub.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/display/src/avesub.for
 
1
C===========================================================================
 
2
C Copyright (C) 1995,2003 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 Massachusetts 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
      PROGRAM AVESUB
 
29
C
 
30
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 
31
C
 
32
C.IDENTIFICATION:                
 
33
C  program AVESUB                       version 3.00    881028
 
34
C  K. Banse                             ESO - Garching
 
35
C
 
36
C.KEYWORDS:
 
37
C  statistics, subimage
 
38
C
 
39
C.PURPOSE:
 
40
C  calculate the average value (and rms) according to different methods 
 
41
C  of an area defined by the cursor rectangle or a table or by coordinates
 
42
C
 
43
C.ALGORITHM:
 
44
C  This info may also be stored either in a descriptor or a table...
 
45
C  
 
46
C.INPUT/OUTPUT:
 
47
C  
 
48
C  The following keywords are used:
 
49
C  IN_A/C/1/60                          input frame
 
50
C  IN_B/C/1/60                          input specs,
 
51
C                                       = +, for cursor input
 
52
C                                       = name or name/TABLE for table input
 
53
C  ACTION/C/1/1                         averaging method,
 
54
C                                       K for kappa-sigma clipping,
 
55
C                                       M for median, A for simple average
 
56
C  INPUTC/C/1/50                        specs for data storage
 
57
C                                       = :label for table or table,:label (a)
 
58
C                                       = descr for descriptor (b)
 
59
C                                       = + for just display of data
 
60
C  INPUTC/C/51/1                        (b) append flag for descriptor,
 
61
C                                       = A, append data to descriptor
 
62
C                                       = +, start again at position 1 in descr
 
63
C                                       (a) center flag for table
 
64
C                                       = C, add also columns with center points
 
65
C                                       = +, no
 
66
C       
 
67
C  OUTPUTR/R/1/5                        last displayed results
 
68
C  OUTPUTI/I/1/1                        no. of loops in program
 
69
C       
 
70
C.VERSIONS
 
71
C  3.00         use new IDI interfaces
 
72
C       
 
73
C 031021        last modif
 
74
C 
 
75
C--------------------------------------------------------------------------
 
76
C
 
77
      IMPLICIT NONE
 
78
C      
 
79
      INTEGER      FELEM,IAV,INFLAG,ISTOP
 
80
      INTEGER      KAPITR,LEAP,NMAL,NROW,NN,N
 
81
      INTEGER      OUTFLG,STAT1,STAT2,STAT
 
82
      INTEGER      TPIX(2),MAXSIZ
 
83
      INTEGER      TBCOLN(7),OUTCOL
 
84
      INTEGER      SUBPIX(2,2),IPIXS(4),NPON(5)
 
85
      INTEGER      ICUR1(2),ICUR2(2)
 
86
      INTEGER*8    PNTRA,PNTRX
 
87
      INTEGER      NAXIS,NPIX(2),IMNOA
 
88
      INTEGER      TABNUL(7),COOS(4),CURSOF
 
89
      INTEGER      UNI(1),NULO,TID
 
90
      INTEGER      MADRID(1)
 
91
C
 
92
      CHARACTER*80 FRAME,TABLE
 
93
      CHARACTER    DESCR*48,APPFLG*1
 
94
      CHARACTER    CUNIT*48,IDENT*72,CBUF*80,ACTION*1
 
95
      CHARACTER    INFO*48
 
96
      CHARACTER    OUTUNI*16,OUTLAB(7)*16,DRAWY*8
 
97
C      
 
98
      REAL         PCUR1(6),PCUR2(6),PCEN(3),RDUM
 
99
      REAL         BGR,BGVL(5),RMSV(5),AVAREA
 
100
      REAL         REAL1(5),RBUF(10)
 
101
C 
 
102
      DOUBLE PRECISION START(2),STEP(2)
 
103
C
 
104
      LOGICAL      SELFLG
 
105
C 
 
106
      INCLUDE 'MID_INCLUDE:ST_DEF.INC'
 
107
C
 
108
      INCLUDE 'MID_INCLUDE:IDIDEV.INC'
 
109
      INCLUDE 'MID_INCLUDE:IDIMEM.INC'
 
110
C
 
111
      COMMON   /VMR/ MADRID
 
112
C 
 
113
      DATA  SUBPIX   /4*1/
 
114
      DATA  COOS     /-1,-1,-1,-1/,  CURSOF /0/
 
115
      DATA  TID      /-1/
 
116
      DATA  INFO     /'switch cursor(s) on - next time we exit... '/
 
117
C
 
118
      DATA      OUTLAB /'XSTART ','YSTART ','XEND ','YEND ',
 
119
     +                  ' ','XCEN ','YCEN '/
 
120
      DATA      OUTUNI /' '/
 
121
      DATA      REAL1  /5*0./
 
122
      DATA      CUNIT  /' '/,    IDENT  /' '/, DESCR /' '/
 
123
      DATA      CBUF  /' '/,  FRAME  /' '/, TABLE /' '/
 
124
C 
 
125
      INCLUDE 'MID_INCLUDE:ST_DAT.INC'
 
126
C
 
127
C  set up MIDAS environment + enable automatic error abort
 
128
      CALL STSPRO('AVESUB')
 
129
C
 
130
C  get method 
 
131
      CALL STKRDC('ACTION',1,1,1,IAV,ACTION,UNI,NULO,STAT)
 
132
      IF (INDEX('AKM',ACTION).LE.0) ACTION = 'A'            !default to A(verage)
 
133
      CALL STKRDI('INPUTI',1,2,IAV,NPON,UNI,NULO,STAT) 
 
134
C      
 
135
      IF (ACTION.EQ.'M') THEN
 
136
         MAXSIZ = NPON(2)
 
137
         CALL STFXMP(MAXSIZ,D_R4_FORMAT,PNTRX,STAT)
 
138
      ELSE IF (ACTION.EQ.'K') THEN                    !for kappa-sigma clipping
 
139
         KAPITR = NPON(1)
 
140
      ENDIF
 
141
C      
 
142
C  get the subarea input options
 
143
      CALL STKRDC('IN_B',1,1,80,IAV,TABLE,UNI,NULO,STAT)
 
144
      IF (TABLE(1:1).EQ.'+') THEN
 
145
         INFLAG = 1                              !INFLAG = 1 for cursor
 
146
         DRAWY(1:7) = 'N YY?C0'
 
147
         CALL STKRDC('P4',1,1,1,IAV,DRAWY(2:2),UNI,NULO,STAT)
 
148
      ELSE
 
149
         CALL STKRDC('IN_A',1,1,80,IAV,FRAME,UNI,NULO,STAT)
 
150
         CALL CLNFRA(FRAME,FRAME,0)
 
151
         CALL STIGET(FRAME,D_R4_FORMAT,F_I_MODE,F_IMA_TYPE,
 
152
     +               2,NAXIS,NPIX,START,STEP,IDENT,
 
153
     +               CUNIT,PNTRA,IMNOA,STAT)
 
154
         IF (NAXIS.LT.2) 
 
155
     +      CALL STETER(77,'1-dim frames currently not supported...')
 
156
         CALL CLNTAB(TABLE,TABLE,0)
 
157
         INFLAG = 2                              !INFLAG = 2 for table
 
158
      ENDIF
 
159
C
 
160
C  get descriptor or table label for storage of data + output_options
 
161
      DESCR(1:1) = '+'
 
162
      OUTFLG = 1                          !OUTFLG = 1 for descriptor or nothing
 
163
      CALL STKRDC('INPUTC',1,1,51,IAV,CBUF,UNI,NULO,STAT)
 
164
      APPFLG = CBUF(51:51)
 
165
      IF (APPFLG .EQ. 'c') THEN                !take care of lower case options
 
166
         APPFLG = 'C'
 
167
      ELSE IF (APPFLG .EQ. 'a') THEN
 
168
         APPFLG = 'A'
 
169
      ENDIF
 
170
      CBUF(51:51) = ' '
 
171
      IF (CBUF(1:1).EQ.'+') GOTO 500           !neither table nor descriptor
 
172
C      
 
173
C  test for :label or table,:label or #no. or table,#no.
 
174
      IF (INDEX(CBUF,'#') .GT. 0) 
 
175
     +   CALL STETER(73,'We need a column label, no # sign, please.')
 
176
C 
 
177
      IF (CBUF(1:1).EQ.':') THEN
 
178
         N = 1
 
179
      ELSE
 
180
         N = INDEX(CBUF,',:')           
 
181
         IF (N.EQ.1) THEN
 
182
            CALL STETER(76,'Invalid syntax (missing table name)...')
 
183
         ELSE IF (N.GT.1) THEN
 
184
            CALL CLNTAB(CBUF(1:N-1),TABLE,0)      
 
185
            N = N + 1                !skip the ','
 
186
         ELSE
 
187
            GOTO 400                     !it's a descriptor
 
188
         ENDIF
 
189
      ENDIF
 
190
C      
 
191
C  handle table storage
 
192
      OUTFLG = 2                               !OUTFLG = 2 for table
 
193
      OUTLAB(5) = CBUF(N+1:50)
 
194
      IF (INFLAG.EQ.1) THEN                    !cursor input
 
195
         DO 80, N=1,7
 
196
            IF (N.NE.5) THEN
 
197
               IF (OUTLAB(5).EQ.OUTLAB(N))
 
198
     +            CALL STETER(78,'Invalid column label !')
 
199
            ENDIF
 
200
80       CONTINUE
 
201
C 
 
202
C  for cursor input create table with 100 rows
 
203
         IF (APPFLG.EQ.'C') THEN
 
204
            OUTCOL = 7
 
205
         ELSE
 
206
            OUTCOL = 5
 
207
         ENDIF
 
208
         CALL TBTINI(TABLE,0,F_O_MODE,OUTCOL,100,TID,STAT) 
 
209
         DO 100, N=1,OUTCOL
 
210
            CALL TBCINI(TID,D_R4_FORMAT,1,'E12.5',OUTUNI,OUTLAB(N),
 
211
     +                  TBCOLN(N),STAT)
 
212
100      CONTINUE
 
213
         GOTO 500
 
214
      ELSE
 
215
         IF (N.GT.1)           !for table input is table,:label not allowed...
 
216
     +      CALL STETER(79,'No table name possible in out_specs... ')
 
217
      ENDIF
 
218
C      
 
219
C  for existing table, test, if label(s) already there
 
220
      CALL TBTOPN(TABLE,F_IO_MODE,TID,STAT)
 
221
      IF (APPFLG.EQ.'C') THEN
 
222
         NN = 7
 
223
      ELSE
 
224
         NN = 5
 
225
      ENDIF
 
226
      DO 150, N=5,NN
 
227
         CALL TBLSER(TID,OUTLAB(N),TBCOLN(N),STAT)
 
228
         IF (TBCOLN(N).LE.0)                          !No. We have to define it...
 
229
     +      CALL TBCINI(TID,D_R4_FORMAT,1,'E12.5',OUTUNI,OUTLAB(N),
 
230
     +                  TBCOLN(N),STAT)
 
231
150   CONTINUE
 
232
      GOTO 500
 
233
C      
 
234
C  handle descriptor storage
 
235
400   DESCR(1:48) = CBUF(1:48)
 
236
      IF (APPFLG.EQ.'A') THEN                   !see, if data is to be appended
 
237
         FELEM = -1
 
238
      ELSE
 
239
         FELEM = 1
 
240
      ENDIF
 
241
C      
 
242
C  branch according to input option
 
243
500   NMAL = 0                                  !init counter
 
244
      IF (INFLAG.EQ.2) GOTO 2000
 
245
C      
 
246
C  a) cursor defined subimage(s) - attach Image Display + get all info
 
247
C 
 
248
      CALL DTOPEN(1,STAT)
 
249
      CALL SETCUR(QDSPNO,2,1,2,COOS,STAT)       !set up cursor rectangle
 
250
      FRAME(1:) = ' '
 
251
C
 
252
C  input cursor positions and check status
 
253
1000  CALL GETCUR(DRAWY,FRAME,                  !draw box in overlay channel
 
254
     +            ICUR1,PCUR1(3),PCUR1(5),RDUM,STAT1,
 
255
     +            ICUR2,PCUR2(3),PCUR2(5),RDUM,STAT2)
 
256
C
 
257
C  check cursor status  
 
258
      IF ( (STAT1.EQ.0) .AND. (STAT2.EQ.0) ) THEN
 
259
         IF (NMAL.EQ.0) THEN
 
260
            IF (CURSOF.EQ.1) GOTO 9000
 
261
            CALL STTPUT(INFO,STAT)
 
262
            FRAME(1:) = ' '
 
263
            CURSOF = 1
 
264
            GOTO 1000
 
265
         ELSE
 
266
            GOTO 9000                             !cursor loop terminated ...
 
267
         ENDIF
 
268
      ENDIF
 
269
C      
 
270
      NMAL = NMAL + 1                             !update counter
 
271
C      
 
272
C  fill output buffer
 
273
      REAL1(1) = PCUR1(5)
 
274
      REAL1(2) = PCUR1(6)
 
275
      REAL1(3) = PCUR2(5)
 
276
      REAL1(4) = PCUR2(6)
 
277
C      
 
278
C  finally map displayed frame
 
279
      CALL STIGET(FRAME,D_R4_FORMAT,F_I_MODE,F_IMA_TYPE,2,
 
280
     +            NAXIS,NPIX,START,STEP,
 
281
     +            IDENT,CUNIT,PNTRA,IMNOA,STAT)
 
282
      GOTO 3000
 
283
C      
 
284
C  b) table defined subimage(s) 
 
285
C 
 
286
2000  IF (OUTFLG.NE.2) CALL TBTOPN(TABLE,F_IO_MODE,TID,STAT)
 
287
      CALL TBIGET(TID,N,NROW,N,N,N,STAT)              !get total no. of rows
 
288
      DO 2020, N=1,4
 
289
         CALL TBLSER(TID,OUTLAB(N),TBCOLN(N),STAT)
 
290
2020  CONTINUE
 
291
C      
 
292
C  here we loop
 
293
2100  IF (NMAL.GE.NROW) GOTO 9000
 
294
      NMAL = NMAL + 1
 
295
      CALL TBSGET(TID,NMAL,SELFLG,STAT)             !only use selected rows
 
296
      IF (.NOT.SELFLG) GOTO 2100
 
297
C      
 
298
C  get next row of values
 
299
      CALL TBRRDR(TID,NMAL,4,TBCOLN,REAL1(1),TABNUL,STAT)
 
300
      PCUR1(3) = (REAL1(1)-START(1))/STEP(1) + 1.
 
301
      PCUR1(4) = (REAL1(2)-START(2))/STEP(2) + 1.
 
302
      PCUR2(3) = (REAL1(3)-START(1))/STEP(1) + 1.
 
303
      PCUR2(4) = (REAL1(4)-START(2))/STEP(2) + 1.
 
304
C 
 
305
C  -----------------------------------
 
306
C  common section for option a) and b)
 
307
C  -----------------------------------
 
308
C 
 
309
3000  SUBPIX(1,1) = PCUR1(3)
 
310
      SUBPIX(2,1) = PCUR1(4)
 
311
      SUBPIX(1,2) = PCUR2(3)
 
312
      SUBPIX(2,2) = PCUR2(4)
 
313
C      
 
314
      IF (APPFLG.EQ.'C') THEN
 
315
         N = (SUBPIX(1,1)+SUBPIX(1,2)+1) / 2
 
316
         PCEN(2) = (N-1)*STEP(1) + START(1)
 
317
         N = (SUBPIX(2,1)+SUBPIX(2,2)+1) / 2
 
318
         PCEN(3) = (N-1)*STEP(2) + START(2)
 
319
      ENDIF
 
320
C      
 
321
C  and finally let's do the averaging 
 
322
C      
 
323
      IF (ACTION.EQ.'K') THEN
 
324
         IPIXS(1) = SUBPIX(1,1)
 
325
         IPIXS(2) = SUBPIX(1,2)
 
326
         IPIXS(3) = SUBPIX(2,1)
 
327
         IPIXS(4) = SUBPIX(2,2)
 
328
         LEAP = 0
 
329
         ISTOP = 0
 
330
         CALL BGVAL(1,MADRID(PNTRA),NPIX(1),NPIX(2),IPIXS,KAPITR,
 
331
     +              2,LEAP,BGVL,RMSV,NPON,ISTOP,N)
 
332
         BGR = BGVL(1)
 
333
CC       SIG = RMSV(1)
 
334
CC       NPN = NPON(1)
 
335
      ELSE
 
336
         TPIX(1) = SUBPIX(1,2) - SUBPIX(1,1) + 1            !get size of subimage
 
337
         TPIX(2) = SUBPIX(2,2) - SUBPIX(2,1) + 1
 
338
         IF ((ACTION.EQ.'M').AND.(TPIX(1)*TPIX(2).GT.MAXSIZ)) THEN
 
339
            CALL STTPUT('subimage too large for internal buffers -'//
 
340
     +                  ' skipped...',STAT)
 
341
            IF (INFLAG.EQ.1) NMAL = NMAL - 1                !modify counter
 
342
            GOTO 7000
 
343
         ENDIF
 
344
         BGR = AVAREA(ACTION,MADRID(PNTRA),NPIX(1),SUBPIX,TPIX,
 
345
     +                MADRID(PNTRX))
 
346
      ENDIF
 
347
C
 
348
C  display results
 
349
      IF (NMAL.EQ.1) THEN                               !here for the 1. time...?
 
350
         CBUF(1:) = '   xstart        ystart          xend           '//
 
351
     +              'yend        average_val'
 
352
         CALL STTPUT(CBUF,STAT)
 
353
      ENDIF
 
354
      REAL1(5) = BGR
 
355
      WRITE(CBUF,20000) (REAL1(N),N=1,5)
 
356
      CALL STTPUT(CBUF,STAT)
 
357
C
 
358
C  fill table, if applicable
 
359
      IF (OUTFLG.EQ.2) THEN
 
360
         IF (INFLAG.EQ.2) THEN                  !table input
 
361
            PCEN(1) = BGR
 
362
            IF (APPFLG.EQ.'C') THEN
 
363
               N = 3
 
364
            ELSE
 
365
               N = 1
 
366
            ENDIF                               !row already exists
 
367
            CALL TBRWRR(TID,NMAL,N,TBCOLN(5),PCEN,STAT)
 
368
         ELSE                                    !cursor input
 
369
            DO 3500, N=1,5
 
370
               RBUF(N) = REAL1(N)
 
371
3500        CONTINUE
 
372
            IF (APPFLG.EQ.'C') THEN
 
373
               RBUF(6) = PCEN(2)
 
374
               RBUF(7) = PCEN(3)
 
375
            ENDIF
 
376
            CALL TBRWRR(TID,NMAL,OUTCOL,TBCOLN,RBUF,STAT)
 
377
         ENDIF
 
378
      ELSE
 
379
C      
 
380
C  fill descriptor, if applicable
 
381
         IF (DESCR(1:1).NE.'+') THEN
 
382
            CALL STDWRR(IMNOA,DESCR,REAL1,FELEM,5,UNI,STAT) 
 
383
            FELEM = -1
 
384
         ENDIF
 
385
      ENDIF
 
386
C      
 
387
C  store info also in keyword OUTPUTR
 
388
      CALL STKWRR('OUTPUTR',REAL1(1),1,5,UNI,STAT)
 
389
 
 
390
C      
 
391
C  and loop again
 
392
7000  GOTO (1000,2100),INFLAG
 
393
C  
 
394
C  That's it folks...
 
395
9000  IF ((TID.NE.-1) .AND. (OUTFLG.EQ.2)) THEN
 
396
         CALL TBSINI(TID,STAT)
 
397
         CALL TBTCLO(TID,STAT)
 
398
      ENDIF
 
399
C      
 
400
C  save no. of coordinates obtained for subsequent applications
 
401
      CALL STKWRI('OUTPUTI',NMAL,1,1,UNI,STAT)
 
402
C      
 
403
      IF (INFLAG.EQ.1) CALL DTCLOS(QDSPNO)
 
404
      CALL STSEPI
 
405
C      
 
406
C  format statements
 
407
20000 FORMAT(1X,5G15.7)
 
408
C 
 
409
      END