C @(#)avesub.for 17.1.1.1 (ES0-DMD) 01/25/02 17:39:57 C=========================================================================== C Copyright (C) 1995 European Southern Observatory (ESO) C C This program is free software; you can redistribute it and/or C modify it under the terms of the GNU General Public License as C published by the Free Software Foundation; either version 2 of C the License, or (at your option) any later version. C C This program is distributed in the hope that it will be useful, C but WITHOUT ANY WARRANTY; without even the implied warranty of C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the C GNU General Public License for more details. C C You should have received a copy of the GNU General Public C License along with this program; if not, write to the Free C Software Foundation, Inc., 675 Massachusetts Ave, Cambridge, C MA 02139, USA. C C Correspondence concerning ESO-MIDAS should be addressed as follows: C Internet e-mail: midas@eso.org C Postal address: European Southern Observatory C Data Management Division C Karl-Schwarzschild-Strasse 2 C D 85748 Garching bei Muenchen C GERMANY C=========================================================================== C PROGRAM AVESUB C C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C.IDENTIFICATION: C program AVESUB version 3.00 881028 C K. Banse ESO - Garching C 3.02 900531 3.20 900719 C C.KEYWORDS: C statistics, subimage C C.PURPOSE: C calculate the average value (and rms) according to different methods C of an area defined by the cursor rectangle or a table or by coordinates C C.ALGORITHM: C This info may also be stored either in a descriptor or a table... C C.INPUT/OUTPUT: C C The following keywords are used: C IN_A/C/1/60 input frame C IN_B/C/1/60 input specs, C = +, for cursor input C = name or name/TABLE for table input C ACTION/C/1/1 averaging method, C K for kappa-sigma clipping, C M for median, A for simple average C INPUTC/C/1/50 specs for data storage C = :label for table or table,:label (a) C = descr/DESCR for descriptor (b) C = + for just display of data C INPUTC/C/51/1 (b) append flag for descriptor, C = A, append data to descriptor C = +, start again at position 1 in descr C (a) center flag for table C = C, add also columns with center points C = +, no C C OUTPUTR/R/1/5 last displayed results C OUTPUTI/I/1/1 no. of loops in program C C.VERSIONS C 3.00 use new IDI interfaces C 3.02 use DAZCLO only if DAZOPN was done... C 3.20 fix all the bugs reported in the Midas Verification Form C C 990805 C-------------------------------------------------------------------------- C IMPLICIT NONE C INTEGER FELEM,IAV,INFLAG,ISTOP INTEGER KAPITR,LEAP,NMAL,NROW,NN,N INTEGER OUTFLG,STAT1,STAT2,STAT INTEGER TPIX(2),MAXSIZ INTEGER TBCOLN(7),OUTCOL INTEGER SUBPIX(2,2),IPIXS(4),NPON(5) INTEGER ICUR1(2),ICUR2(2) INTEGER*8 PNTRA,PNTRX INTEGER NAXIS,NPIX(2),IMNOA INTEGER TABNUL(7),COOS(4),CURSOF INTEGER UNI(1),NULO,TID INTEGER MADRID(1) C CHARACTER*80 FRAME,TABLE CHARACTER DESCR*15,APPFLG*1 CHARACTER CUNIT*48,IDENT*72,CBUF*80,ACTION*1 CHARACTER INFO*48 CHARACTER OUTUNI*16,OUTLAB(7)*16,DRAWY*8 C REAL PCUR1(6),PCUR2(6),PCEN(3),RDUM REAL BGR,BGVL(5),RMSV(5),AVAREA REAL REAL1(5),RBUF(10) C DOUBLE PRECISION START(2),STEP(2) C LOGICAL SELFLG C INCLUDE 'MID_INCLUDE:ST_DEF.INC' C INCLUDE 'MID_INCLUDE:IDIDEV.INC' INCLUDE 'MID_INCLUDE:IDIMEM.INC' C COMMON /VMR/ MADRID C DATA SUBPIX /4*1/ DATA COOS /-1,-1,-1,-1/, CURSOF /0/ DATA TID /-1/ DATA INFO /'switch cursor(s) on - next time we exit... '/ C DATA OUTLAB /'XSTART ','YSTART ','XEND ','YEND ', + ' ','XCEN ','YCEN '/ DATA OUTUNI /' '/ DATA REAL1 /5*0./ DATA CUNIT /' '/, IDENT /' '/ DATA CBUF /' '/, FRAME /' '/, TABLE /' '/ C INCLUDE 'MID_INCLUDE:ST_DAT.INC' C C set up MIDAS environment + enable automatic error abort CALL STSPRO('AVESUB') C C get method CALL STKRDC('ACTION',1,1,1,IAV,ACTION,UNI,NULO,STAT) IF (INDEX('AKM',ACTION).LE.0) ACTION = 'A' !default to A(verage) CALL STKRDI('INPUTI',1,2,IAV,NPON,UNI,NULO,STAT) C IF (ACTION.EQ.'M') THEN MAXSIZ = NPON(2) CALL STFXMP(MAXSIZ,D_R4_FORMAT,PNTRX,STAT) ELSE IF (ACTION.EQ.'K') THEN !for kappa-sigma clipping KAPITR = NPON(1) ENDIF C C get the subarea input options CALL STKRDC('IN_B',1,1,80,IAV,TABLE,UNI,NULO,STAT) IF (TABLE(1:1).EQ.'+') THEN INFLAG = 1 !INFLAG = 1 for cursor DRAWY(1:7) = 'N YY?C0' CALL STKRDC('P4',1,1,1,IAV,DRAWY(2:2),UNI,NULO,STAT) ELSE CALL STKRDC('IN_A',1,1,80,IAV,FRAME,UNI,NULO,STAT) CALL CLNFRA(FRAME,FRAME,0) CALL STIGET(FRAME,D_R4_FORMAT,F_I_MODE,F_IMA_TYPE, + 2,NAXIS,NPIX,START,STEP,IDENT, + CUNIT,PNTRA,IMNOA,STAT) IF (NAXIS.LT.2) + CALL STETER(77,'1-dim frames currently not supported...') CALL CLNTAB(TABLE,TABLE,0) INFLAG = 2 !INFLAG = 2 for table ENDIF C C get descriptor or table label for storage of data + output_options OUTFLG = 1 !OUTFLG = 1 for descriptor or nothing CALL STKRDC('INPUTC',1,1,51,IAV,CBUF,UNI,NULO,STAT) APPFLG = CBUF(51:51) IF (APPFLG .EQ. 'c') THEN !take care of lower case options APPFLG = 'C' ELSE IF (APPFLG .EQ. 'a') THEN APPFLG = 'A' ENDIF CBUF(51:51) = ' ' IF (CBUF(1:1).EQ.'+') GOTO 500 !neither table nor descriptor C C test for :label or table,:label or #no. or table,#no. IF (INDEX(CBUF,'#') .GT. 0) + CALL STETER(73,'We need a column label, no # sign, please.') C IF (CBUF(1:1).EQ.':') THEN N = 1 ELSE N = INDEX(CBUF,',:') IF (N.EQ.1) THEN CALL STETER(76,'Invalid syntax (missing table name)...') ELSE IF (N.GT.1) THEN CALL CLNTAB(CBUF(1:N-1),TABLE,0) N = N + 1 !skip the ',' ELSE GOTO 400 !it's a descriptor ENDIF ENDIF C C handle table storage OUTFLG = 2 !OUTFLG = 2 for table OUTLAB(5) = CBUF(N+1:50) IF (INFLAG.EQ.1) THEN !cursor input DO 80, N=1,7 IF (N.NE.5) THEN IF (OUTLAB(5).EQ.OUTLAB(N)) + CALL STETER(78,'Invalid column label !') ENDIF 80 CONTINUE C C for cursor input create table with 100 rows IF (APPFLG.EQ.'C') THEN OUTCOL = 7 ELSE OUTCOL = 5 ENDIF CALL TBTINI(TABLE,0,F_O_MODE,OUTCOL,100,TID,STAT) DO 100, N=1,OUTCOL CALL TBCINI(TID,D_R4_FORMAT,1,'E12.5',OUTUNI,OUTLAB(N), + TBCOLN(N),STAT) 100 CONTINUE GOTO 500 ELSE IF (N.GT.1) !for table input is table,:label not allowed... + CALL STETER(79,'No table name possible in out_specs... ') ENDIF C C for existing table, test, if label(s) already there CALL TBTOPN(TABLE,F_IO_MODE,TID,STAT) IF (APPFLG.EQ.'C') THEN NN = 7 ELSE NN = 5 ENDIF DO 150, N=5,NN CALL TBLSER(TID,OUTLAB(N),TBCOLN(N),STAT) IF (TBCOLN(N).LE.0) !No. We have to define it... + CALL TBCINI(TID,D_R4_FORMAT,1,'E12.5',OUTUNI,OUTLAB(N), + TBCOLN(N),STAT) 150 CONTINUE GOTO 500 C C handle descriptor storage 400 DESCR(1:15) = CBUF(1:15) IF (APPFLG.EQ.'A') THEN !see, if data is to be appended FELEM = -1 ELSE FELEM = 1 ENDIF C C branch according to input option 500 NMAL = 0 !init counter IF (INFLAG.EQ.2) GOTO 2000 C C a) cursor defined subimage(s) - attach Image Display + get all info C CALL DTOPEN(1,STAT) CALL SETCUR(QDSPNO,2,1,2,COOS,STAT) !set up cursor rectangle FRAME(1:) = ' ' C C input cursor positions and check status 1000 CALL GETCUR(DRAWY,FRAME, !draw box in overlay channel + ICUR1,PCUR1(3),PCUR1(5),RDUM,STAT1, + ICUR2,PCUR2(3),PCUR2(5),RDUM,STAT2) C C check cursor status IF ( (STAT1.EQ.0) .AND. (STAT2.EQ.0) ) THEN IF (NMAL.EQ.0) THEN IF (CURSOF.EQ.1) GOTO 9000 CALL STTPUT(INFO,STAT) FRAME(1:) = ' ' CURSOF = 1 GOTO 1000 ELSE GOTO 9000 !cursor loop terminated ... ENDIF ENDIF C NMAL = NMAL + 1 !update counter C C fill output buffer REAL1(1) = PCUR1(5) REAL1(2) = PCUR1(6) REAL1(3) = PCUR2(5) REAL1(4) = PCUR2(6) C C finally map displayed frame CALL STIGET(FRAME,D_R4_FORMAT,F_I_MODE,F_IMA_TYPE,2, + NAXIS,NPIX,START,STEP, + IDENT,CUNIT,PNTRA,IMNOA,STAT) GOTO 3000 C C b) table defined subimage(s) C 2000 IF (OUTFLG.NE.2) CALL TBTOPN(TABLE,F_IO_MODE,TID,STAT) CALL TBIGET(TID,N,NROW,N,N,N,STAT) !get total no. of rows DO 2020, N=1,4 CALL TBLSER(TID,OUTLAB(N),TBCOLN(N),STAT) 2020 CONTINUE C C here we loop 2100 IF (NMAL.GE.NROW) GOTO 9000 NMAL = NMAL + 1 CALL TBSGET(TID,NMAL,SELFLG,STAT) !only use selected rows IF (.NOT.SELFLG) GOTO 2100 C C get next row of values CALL TBRRDR(TID,NMAL,4,TBCOLN,REAL1(1),TABNUL,STAT) PCUR1(3) = (REAL1(1)-START(1))/STEP(1) + 1. PCUR1(4) = (REAL1(2)-START(2))/STEP(2) + 1. PCUR2(3) = (REAL1(3)-START(1))/STEP(1) + 1. PCUR2(4) = (REAL1(4)-START(2))/STEP(2) + 1. C C ----------------------------------- C common section for option a) and b) C ----------------------------------- C 3000 SUBPIX(1,1) = PCUR1(3) SUBPIX(2,1) = PCUR1(4) SUBPIX(1,2) = PCUR2(3) SUBPIX(2,2) = PCUR2(4) C IF (APPFLG.EQ.'C') THEN N = (SUBPIX(1,1)+SUBPIX(1,2)+1) / 2 PCEN(2) = (N-1)*STEP(1) + START(1) N = (SUBPIX(2,1)+SUBPIX(2,2)+1) / 2 PCEN(3) = (N-1)*STEP(2) + START(2) ENDIF C C and finally let's do the averaging C IF (ACTION.EQ.'K') THEN IPIXS(1) = SUBPIX(1,1) IPIXS(2) = SUBPIX(1,2) IPIXS(3) = SUBPIX(2,1) IPIXS(4) = SUBPIX(2,2) LEAP = 0 ISTOP = 0 CALL BGVAL(1,MADRID(PNTRA),NPIX(1),NPIX(2),IPIXS,KAPITR, + 2,LEAP,BGVL,RMSV,NPON,ISTOP,N) BGR = BGVL(1) CC SIG = RMSV(1) CC NPN = NPON(1) ELSE TPIX(1) = SUBPIX(1,2) - SUBPIX(1,1) + 1 !get size of subimage TPIX(2) = SUBPIX(2,2) - SUBPIX(2,1) + 1 IF ((ACTION.EQ.'M').AND.(TPIX(1)*TPIX(2).GT.MAXSIZ)) THEN CALL STTPUT('subimage too large for internal buffers -'// + ' skipped...',STAT) IF (INFLAG.EQ.1) NMAL = NMAL - 1 !modify counter GOTO 7000 ENDIF BGR = AVAREA(ACTION,MADRID(PNTRA),NPIX(1),SUBPIX,TPIX, + MADRID(PNTRX)) ENDIF C C display results IF (NMAL.EQ.1) THEN !here for the 1. time...? CBUF(1:) = ' xstart ystart xend '// + 'yend average_val' CALL STTPUT(CBUF,STAT) ENDIF REAL1(5) = BGR WRITE(CBUF,20000) (REAL1(N),N=1,5) CALL STTPUT(CBUF,STAT) C C fill table, if applicable IF (OUTFLG.EQ.2) THEN IF (INFLAG.EQ.2) THEN !table input PCEN(1) = BGR IF (APPFLG.EQ.'C') THEN N = 3 ELSE N = 1 ENDIF !row already exists CALL TBRWRR(TID,NMAL,N,TBCOLN(5),PCEN,STAT) ELSE !cursor input DO 3500, N=1,5 RBUF(N) = REAL1(N) 3500 CONTINUE IF (APPFLG.EQ.'C') THEN RBUF(6) = PCEN(2) RBUF(7) = PCEN(3) ENDIF CALL TBRWRR(TID,NMAL,OUTCOL,TBCOLN,RBUF,STAT) ENDIF ELSE C C fill descriptor, if applicable IF (DESCR(1:1).NE.'+') THEN CALL STDWRR(IMNOA,DESCR,REAL1,FELEM,5,UNI,STAT) FELEM = -1 ENDIF ENDIF C C store info also in keyword OUTPUTR CALL STKWRR('OUTPUTR',REAL1(1),1,5,UNI,STAT) C C and loop again 7000 GOTO (1000,2100),INFLAG C C That's it folks... 9000 IF ((TID.NE.-1) .AND. (OUTFLG.EQ.2)) THEN CALL TBSINI(TID,STAT) CALL TBTCLO(TID,STAT) ENDIF C C save no. of coordinates obtained for subsequent applications CALL STKWRI('OUTPUTI',NMAL,1,1,UNI,STAT) C IF (INFLAG.EQ.1) CALL DTCLOS(QDSPNO) CALL STSEPI C C format statements 20000 FORMAT(1X,5G15.7) C END