C @(#)statis.for 17.1.1.1 (ESO-DMD) 01/25/02 17:40:02 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 STATIS C C++++++++++++++++++++++++++++++++++++++++++++++++++ C C.COPYRIGHT: Copyright (c) 1988 European Southern Observatory, C all rights reserved C C.LANGUAGE: F77+ESOext C C.AUTHOR: K.Banse C C.IDENTIFICATION C program STATIS 1.00 881028 C K. Banse ESO - Garching C C.KEYWORDS C maximum, minimum,mean value, std. dev., histogram C C.PURPOSE C calculate statistics of given bulk data frame C C.ALGORITHM C prepare input for options: COLUMN, ROW, TABLE, CURSOR C do the actual calculations in subroutine ZMSTAT C (former main program STATIST) C C.INPUT/OUTPUT C the following keys are used for I/O: C C IN_A/C/1/80 input frame C P2/C/1/80 subframe area [...:...] C + for complete frame C ROW to work on rows of 2-dim frame C COLUMN to work on columns of 2-dim frame C PLANE to work on planes of 3-dim frame C table_name for table input C CURSOR for cursor defined subimages C DEFAULT/C/1/6 (1) = N, subarea is specified in P2 C = M, subareas are specified in P2 C but merge them all together and do the C statistics on that set C = Y, use complete frame C (2) = N, binsize is given in INPUTR(1) C = Y, total no. of bins is given in INPUTR(1) C (3,4) = option as defined in command STATIST C (5) = Y/N, draw option for rectangles, only used C for CURSOR above C (6) = P/N, for plotting histogram or not C OUT_A/C/1/80 for P2 = ROW, COLUMN, table_name or CURSOR C this holds the output_table, append_flag C C OUTPUTR/R/1/12 will receive: Min,Max,Mean,Std,3rdMom,4thMom,Intensity C and Median,SmallestMode,Nobins,Binsize,Mode C OUTPUTI/I/1/7 (1) = no. valid of pixels in window C (2,3,4) = window start pixels, set to 0, if not applic. C (5,6,7) = window end pixels, set to 0, if not applic. C C.VERSIONS C C 010410 last modif C C-------------------------------------------------- C IMPLICIT NONE C INTEGER IAV,STAT,N,NLOOP,CURMAX INTEGER ST1,ST2,XY1(2),XY2(2),COOS(4) INTEGER NPIX(3),NAXPIX(4),NMAL,NROW,NCOL1,SIZ(2) INTEGER*8 PNTRA,PNTRB,PNTRO INTEGER SPIX(3),EPIX(3) INTEGER TCOLNM(16),ICOLNM(6) INTEGER IACTIO,COOFF,TBAPP INTEGER TABNUL(6) INTEGER UNI(1),NULO,MADRID(1) INTEGER IMNO,IMNOB,IMNOC,INTID,OUTID,OUTNO INTEGER INCOL,KBUF(25) INTEGER NPXO(2),EXAMED INTEGER EC,ED,EL C DOUBLE PRECISION START(3),STEP(3),DVAL DOUBLE PRECISION STRTO(2),STPO(2) DOUBLE PRECISION DD1(4),DD2(4) C REAL RVAL(16),XMEDIA,ZBINS(3) C CHARACTER ORFILE*80,INTABL*80,OUTTBL*80,OUTFRA*80 CHARACTER AREA*100,DEFAUL*6,OUTPUT*100,CCBUF*40 CHARACTER OPTION*5,FRMSTR*20,HISTORY*160 CHARACTER TABUNT*16,TABLAB(12)*16,CACT*8 CHARACTER INTLAB(6)*16 CHARACTER IDENT*72,CUNIT*64,MIDUNI*2,TMPFIL*14 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 COMMON /PLOTTO/ IDENT,CUNIT C DATA COOS /-1,-1,-1,-1/ DATA RVAL /16*0.0/ DATA TABUNT /' '/ DATA TABLAB /'MIN ','MAX ','MEAN ','STD ', + 'MOM3 ','MOM4 ','TOT_INTENS', + 'MEDIAN ','SMAL_MODE ', + 'NO_BINS ','BIN_SIZE ', 'MODE '/ DATA INTLAB /'XSTART ','YSTART ','XEND ','YEND ', + 'ZSTART ','ZEND '/ DATA AREA /' '/ DATA IDENT /' '/, CUNIT /' '/ C INCLUDE 'MID_INCLUDE:ST_DAT.INC' C C init MIDAS environment + get current MIDAS unit C CALL STSPRO('STATIST') CALL STPSET(F_FITS_PARM,1,STAT) !ignore ESO.xyz FITS keywds CALL STKWRR('OUTPUTR',RVAL,1,14,UNI,STAT) !clear keyword OUTPUTR CALL STKRDC('MID$SESS',1,11,2,IAV,MIDUNI,UNI,NULO,STAT) TMPFIL = 'stat'//MIDUNI//'dum.bdf ' C DO 500, N=1,3 NPIX(N) = 1 SPIX(N) = 1 EPIX(N) = 1 START(N) = 0.D0 STEP(N) = 1.D0 500 CONTINUE OUTID = -999 !indicate that no output table/image OUTNO = -999 !indicate that no output image/image INTID = -1 FRMSTR(1:) = ' ' EXAMED = 0 KBUF(1) = -1 !indicates an invalid row/column/plane C C get info about histogram bins + the different options + format string C CALL STKRDC('DEFAULT',1,1,6,IAV,DEFAUL,UNI,NULO,STAT) CALL UPCAS(DEFAUL,DEFAUL) OPTION(1:4) = DEFAUL(1:4) OUTPUT(1:) = ' ' CALL STKRDC('P3',1,1,20,IAV,OUTPUT,UNI,NULO,STAT) IF (OUTPUT(1:1) .EQ. '#') THEN OPTION(2:2) = 'Y' CALL GENCNV(OUTPUT(2:),2,1,ICOLNM,ZBINS,DVAL,IAV) IF (IAV .NE. 1) CALL STETER(34,'Invalid no. of bins...') ELSE OPTION(2:2) = 'N' CALL GENCNV(OUTPUT(1:),2,1,ICOLNM,ZBINS,DVAL,IAV) IF (IAV .NE. 1) CALL STETER(34,'Invalid binsize...') ENDIF OUTPUT(1:) = ' ' CALL STKRDC('P4',1,1,80,IAV,OUTPUT,UNI,NULO,STAT) CALL GENCNV(OUTPUT,2,2,ICOLNM,ZBINS(2),DVAL,IAV) IF (IAV .NE. 2) CALL STETER(35,'Invalid excess bins...') CALL STKRDC('P8',1,1,20,IAV,FRMSTR,UNI,NULO,STAT) C C get name of frame to process + open it + get stand. descr. C CALL STKRDC('IN_A',1,1,80,IAV,ORFILE,UNI,NULO,STAT) CALL STFOPN(ORFILE,D_R4_FORMAT,0,F_IMA_TYPE,IMNO,STAT) CALL STDRDI(IMNO,'NAXIS',1,1,IAV,NAXPIX,UNI,NULO,STAT) IF (NAXPIX(1).GT.3) THEN CALL STTPUT('Warning: only 3 dims of frame processed ...', + STAT) NAXPIX(1) = 3 ENDIF CALL FPXWCO(0,IMNO,DD1,DD2,STAT) IF (STAT .GT. 0) CALL STETER(36,'Problems with WCS of frame...') C CALL STDRDI(IMNO,'NPIX',1,NAXPIX(1),IAV,NPIX,UNI,NULO,STAT) CALL STDRDD(IMNO,'START',1,NAXPIX(1),IAV,START,UNI,NULO,STAT) CALL STDRDD(IMNO,'STEP',1,NAXPIX(1),IAV,STEP,UNI,NULO,STAT) IAV = 1 DO 550, N=1,3 NAXPIX(N+1) = NPIX(N) IAV = IAV * NPIX(N) 550 CONTINUE ICOLNM(1) = IMNO !use ICOLNM as buffer CALL STFINF(ORFILE,-4,ICOLNM,STAT) !ICOLNM(4) = real no. of pixels IF ((IAV .GT. ICOLNM(4)).AND.(ICOLNM(3).NE.1)) + CALL STETER(37,'reference file with no data...') C IF (DEFAUL(6:6).EQ.'P') THEN IF (INDEX(DEFAUL(3:3),'RXSM') .GT. 0) THEN CALL STTPUT + ('No histogram plotting possible for reduced statistics...', + STAT) DEFAUL(6:6) = 'N' ELSE CALL STDRDC(IMNO,'IDENT',1,1,72,IAV,IDENT,UNI,NULO,STAT) N = (NAXPIX(1)+1) * 16 CALL STDRDC(IMNO,'CUNIT',1,1,N,IAV,CUNIT,UNI,NULO,STAT) ENDIF ENDIF C IF ((OPTION(3:3).EQ.'G').OR.(OPTION(3:3).EQ.'X') + .OR. (OPTION(3:3).EQ.'W')) THEN EXAMED = 1 ENDIF C C ********************* C C test the options: full frame, subframes, intervals, C ROW, COLUMN, PLANE, CURSOR or table C C IACTIO specifies the different area options: C = 0, full frame; C = 1, table; = 2, rows; = 3, columns; = 4, cursor; = 5, planes; C = 6, subframe; = 7, covering intervals C C ********************* C CALL STKRDC('P2',1,1,100,IAV,OUTPUT,UNI,NULO,STAT) OPTION(5:5) = 'N' IF ( (DEFAUL(6:6).EQ.'P') .OR. (DEFAUL(1:1).EQ.'M') ) THEN OPTION(5:5) = 'H' !we also have to fill descr. HISTOGRAM ENDIF INTABL(1:) = ' ' CALL UPCAS(OUTPUT,AREA) C C complete frame C IF (DEFAUL(1:1).EQ.'Y') THEN IACTIO = 0 GOTO 1111 ENDIF C C area option = [... C IF (OUTPUT(1:1).EQ.'[') THEN N = INDEX(ORFILE,']') IF (N .GT. 0) THEN IAV = LEN(ORFILE) IF ((N .EQ. IAV) .OR. (ORFILE(N+1:N+1) .EQ. ' ')) + CALL STETER(33,'double subframe option not supported...') ENDIF C C test, if subframe or interval option N = INDEX(OUTPUT,':') IF (N.GT.1) THEN IACTIO = 6 !for subframe(s) ELSE C C covering intervals C IACTIO = 7 N = INDEX(AREA,'X') !look for 123x456 IF (N.LT.3) THEN N = INDEX(AREA,'*') !or for 123*456 IF (N.LT.3) + CALL STETER(2,'invalid interval option...') ENDIF AREA(N:N) = ',' NMAL = INDEX(AREA,']') AREA(NMAL:NMAL) = ',' CALL GENCNV(AREA(2:),1,2,SIZ,RVAL,DVAL,IAV) IF (IAV .NE. 2) + CALL STETER(34,'Invalid x/ysize of intervals...') IF (NPIX(3).GT.1) CALL STTPUT + ('Only 1. plane of 3-dim frame will be processed...',STAT) NMAL = (NPIX(1)-1)/SIZ(1) + 1 NROW = (NPIX(2)-1)/SIZ(2) + 1 ENDIF GOTO 1111 ENDIF C C rows C IF (AREA(1:3).EQ.'ROW') THEN IF (AREA(4:4).EQ.' ') THEN NROW = NPIX(2) ELSE IF (AREA(4:4).EQ.',') THEN CALL GENCNV(AREA(5:),1,25,KBUF,RVAL,DVAL,NROW) IF (NROW .LT. 1) CALL STETER(34,'Invalid row no.') DO 900,N=1,NROW IF ((KBUF(N).LT.1).OR.(KBUF(N).GT.NPIX(2))) + CALL STETER(34,'Invalid row no...') 900 CONTINUE ELSE GOTO 999 !table name starting with ROW.. ENDIF IACTIO = 2 IF (NPIX(3).GT.1) CALL STTPUT + ('Only 1. plane of 3-dim frame will be processed...',STAT) GOTO 1111 ENDIF C C planes C IF (AREA(1:4).EQ.'PLANE') THEN IF (AREA(6:6).EQ.' ') THEN NROW = NPIX(3) ELSE IF (AREA(6:6).EQ.',') THEN CALL GENCNV(AREA(7:),1,25,KBUF,RVAL,DVAL,NROW) IF (NROW .LT. 1) CALL STETER(34,'Invalid plane no.') DO 930,N=1,NROW IF ((KBUF(N).LT.1).OR.(KBUF(N).GT.NPIX(3))) + CALL STETER(34,'Invalid plane no...') 930 CONTINUE ELSE GOTO 999 !table name starting with PLANE.. ENDIF IACTIO = 5 IF (NPIX(3).EQ.1) THEN CALL STTPUT + ('Single plane of 2-dim frame will be processed...',STAT) NROW = 1 ENDIF GOTO 1111 ENDIF C C columns C IF (AREA(1:6).EQ.'COLUMN') THEN IF (AREA(7:7).EQ.' ') THEN NROW = NPIX(1) ELSE IF (AREA(7:7).EQ.',') THEN CALL GENCNV(AREA(8:),1,25,KBUF,RVAL,DVAL,NROW) IF (NROW .LT. 1) CALL STETER(34,'Invalid column no.') DO 960,N=1,NROW IF ((KBUF(N).LT.1).OR.(KBUF(N).GT.NPIX(1))) + CALL STETER(34,'Invalid column no...') 960 CONTINUE ELSE GOTO 999 !table name starting with COLUMN.. ENDIF IACTIO = 3 IF (NPIX(1).EQ.1) + CALL STETER(34,'Invalid option: 1-dim input frame...') IF (NPIX(3).GT.1) CALL STTPUT + ('Only 1. plane of 3-dim frame will be processed...',STAT) GOTO 1111 ENDIF C C cursor input C IF ((AREA(1:7).EQ.'CURSOR ') .OR. + (AREA(1:7).EQ.'CURSOR,')) THEN IACTIO = 4 IF (AREA(7:7) .EQ. ',') THEN CALL GENCNV(AREA(8:),1,1,CURMAX,RVAL,DVAL,IAV) IF ((IAV.LT.1) .OR. (CURMAX.LT.1)) + CALL STETER(2,'Invalid no. for max. cursor input...') NROW = CURMAX ELSE CURMAX = 99999 NROW = 100 ENDIF GOTO 1111 ENDIF C C everything else is interpreted as a table name C 999 IACTIO = 1 !indicate table input CALL CLNTAB(OUTPUT,INTABL,0) CALL TBTOPN(INTABL,F_I_MODE,INTID,STAT) C C get number of rows we'll work on + look for the required columns CALL TBIGET(INTID,N,NROW,N,N,N,STAT) DO 1000,N=1,4 CALL TBLSER(INTID,INTLAB(N),ICOLNM(N),STAT) IF (ICOLNM(N).LE.0) THEN ST1 = INDEX(INTABL,' ') IF (ST1.LE.0) ST1 = LEN(INTABL) ST2 = INDEX(INTLAB(N),' ') IF (ST2.LE.0) ST2 = LEN(INTLAB(N)) WRITE(OUTPUT,10030) INTLAB(N)(1:ST2),INTABL(1:ST1) CALL STETER(3,OUTPUT) ENDIF 1000 CONTINUE C C if NPIX(3) > 1, look for 3-dim table entries INCOL = 4 IF (NPIX(3).GT.1) THEN CALL TBLSER(INTID,INTLAB(5),ICOLNM(5),STAT) IF (ICOLNM(5).GT.0) THEN CALL TBLSER(INTID,INTLAB(6),ICOLNM(6),STAT) IF (ICOLNM(6).GT.0) INCOL = 6 ENDIF ENDIF C C *********** C C after the area option let's look for an output table or image C C *********** C 1111 IF ((DEFAUL(4:4).NE.'N').AND.(IACTIO.NE.4)) + CALL FRAMOU(ORFILE) CALL STKRDC('OUT_A',1,1,80,IAV,OUTPUT,UNI,NULO,STAT) TBAPP = 0 !point to start of table ST2 = 0 !default append flag to 0 N = INDEX(OUTPUT,',') IF (N.GT.1) THEN IF ((OUTPUT(N:N+1).EQ.',i').OR.(OUTTBL(N:N+1).EQ.',I')) THEN CALL CLNFRA(OUTPUT(1:N-1),OUTFRA,0) !save the file name OUTTBL(1:) = 'middummXyZ.tbl ' IF (INTABL.EQ.OUTTBL) OUTTBL(1:) = 'middummXYZ.tbl ' C IAV = INDEX(ORFILE,' ') WRITE(HISTORY,20010) ORFILE(1:IAV-1),OUTPUT OUTNO = 999 OUTPUT(1:) = OUTPUT(N+1:) C C then we also need descr. IDENT + CUNIT CALL STDRDC(IMNO,'IDENT',1,1,72,IAV,IDENT,UNI,NULO,STAT) N = (NAXPIX(1)+1) * 16 CALL STDRDC(IMNO,'CUNIT',1,1,N,IAV,CUNIT,UNI,NULO,STAT) C N = INDEX(OUTPUT,',') !look for qualifier: name,I,qualif IF (N.GT.1) THEN CALL UPCAS(OUTPUT(N+1:),CCBUF) IF (CCBUF(1:1).EQ.'M') THEN IF (CCBUF(2:2).EQ.'I') THEN OUTNO = 1 !minimum ELSE IF (CCBUF(2:2).EQ.'A') THEN OUTNO = 2 !maximum ELSE IF (CCBUF(2:2).EQ.'E') THEN IF (CCBUF(3:3).EQ.'A') THEN OUTNO = 3 !mean ELSE OUTNO = 8 !median ENDIF ELSE IF (CCBUF(2:4).EQ.'ODE') THEN IF (CCBUF(5:5).EQ.'1') THEN OUTNO = 9 !1st mode ELSE OUTNO = 12 !mode ENDIF ELSE IF (CCBUF(2:3).EQ.'OM') THEN IF (CCBUF(4:4).EQ.'3') THEN OUTNO = 5 !3rd moment ELSE OUTNO = 6 !4th moment ENDIF ENDIF ELSE IF (CCBUF(1:2).EQ.'ST') THEN OUTNO = 4 !stddev ELSE IF (CCBUF(1:2).EQ.'TO') THEN OUTNO = 7 !total intensity ENDIF ENDIF C GOTO 1490 !skip the rest ENDIF C IF ((OUTPUT(N:N+1).EQ.',a').OR.(OUTTBL(N:N+1).EQ.',A')) THEN ST2 = 1 ENDIF CALL CLNTAB(OUTPUT(1:N-1),OUTTBL,0) C CALL STECNT('GET',EC,EL,ED) CALL STECNT('PUT',1,0,0) !disable errors ... CALL TBTOPN(OUTTBL,F_IO_MODE,OUTID,STAT) IF (STAT.NE.0) ST2 = 0 !clear append flag again CALL STECNT('PUT',EC,EL,ED) ELSE CALL CLNTAB(OUTPUT,OUTTBL,0) ENDIF C C create NROW * NCOL1 table to have already columns allocated for later use 1490 IF (OUTTBL(1:1).NE.'+') THEN IF (OPTION(3:3).EQ.'M') THEN NCOL1 = 2 !only 2 columns for Minmax ELSE IF (OPTION(3:3).EQ.'H') THEN NCOL1 = 2 !only 2 columns for Histogram ELSE IF (OPTION(3:3).EQ.'S') THEN NCOL1 = 4 !only 4 columns for Short ELSE IF (OPTION(3:3).EQ.'R') THEN NCOL1 = 7 !only 7 columns for Reduced ELSE IF (OPTION(3:3).EQ.'X') THEN NCOL1 = 8 !only 8 columns for XReduced ELSE NCOL1 = 12 ENDIF C IF (INTABL.NE.OUTTBL) THEN IF (ST2.EQ.0) THEN !build new table CALL TBTINI(OUTTBL,0,F_O_MODE,NCOL1+4,NROW,OUTID,STAT) DO 1500 N=1,NCOL1 !and define the columns CALL TBCINI(OUTID,D_R4_FORMAT,1,'G12.6',TABUNT, + TABLAB(N),TCOLNM(N),STAT) 1500 CONTINUE C IF (IACTIO.EQ.4) THEN !4 more columns for CURSOR option DO 1600 N=1,4 IAV = N + NCOL1 CALL TBCINI(OUTID,D_R4_FORMAT,1,'G12.6',TABUNT, + INTLAB(N),TCOLNM(IAV),STAT) 1600 CONTINUE ENDIF ELSE !append to output table C CALL TBIGET(OUTID,N,TBAPP,N,N,N,STAT) !get no. of rows DO 1700 N=1,NCOL1 CALL TBLSER(OUTID,TABLAB(N),TCOLNM(N),STAT) IF (TCOLNM(N).EQ.-1) THEN ST1 = INDEX(OUTTBL,' ') IF (ST1.LE.0) ST1 = LEN(OUTTBL) ST2 = INDEX(TABLAB(N),' ') IF (ST2.LE.0) ST2 = LEN(TABLAB(N)) WRITE(OUTPUT,10020) TABLAB(N)(1:ST2),OUTTBL(1:ST1) CALL STETER(3,OUTPUT) ENDIF 1700 CONTINUE C IF (IACTIO.EQ.4) THEN !4 more columns for CURSOR option DO 1770 N=1,4 IAV = N + NCOL1 CALL TBLSER(OUTID,INTLAB(N),TCOLNM(IAV),STAT) IF (TCOLNM(IAV).EQ.-1) THEN ST1 = INDEX(OUTTBL,' ') IF (ST1.LE.0) ST1 = LEN(OUTTBL) ST2 = INDEX(INTLAB(N),' ') IF (ST2.LE.0) ST2 = LEN(INTLAB(N)) WRITE(OUTPUT,10020) + INTLAB(N)(1:ST2),OUTTBL(1:ST1) CALL STETER(3,OUTPUT) ENDIF 1770 CONTINUE ENDIF ENDIF C C we use same table for input + output ELSE OUTID = INTID DO 1800, N=1,NCOL1 CALL TBLSER(OUTID,TABLAB(N),TCOLNM(N),STAT) IF (TCOLNM(N).EQ.-1) + CALL TBCINI(OUTID,D_R4_FORMAT,1,'G12.6',TABUNT, + TABLAB(N),TCOLNM(N),STAT) 1800 CONTINUE ENDIF ENDIF C C ****************** C C now do the statistics for every area option C C ****************** C C --------- C C 1) process complete frame C C --------- C IF (IACTIO.EQ.0) THEN OPTION(5:5) = 'Y' IF (EXAMED.EQ.1) THEN CALL EXMED(IACTIO,IMNO,ZBINS(2),NPIX, + SPIX,EPIX,XMEDIA,STAT) IF (STAT.NE.0) CALL STETER(35, + 'Problems with exact median calculation...') ENDIF CALL ZMSTAT(IMNO,AREA,NAXPIX,START,STEP, + ZBINS,FRMSTR,OPTION) IF (OUTID.NE.-999) THEN N = 1 CALL STKRDR('OUTPUTR',1,12,IAV,RVAL,UNI,NULO,STAT) CALL TBRWRR(OUTID,N+TBAPP,NCOL1,TCOLNM,RVAL,STAT) ENDIF IF (DEFAUL(6:6).EQ.'P') CALL PLOTI(IMNO) C C --------- C C 2) handle subframe(s) mode C C --------- C ELSE IF (IACTIO.EQ.6) THEN NROW = 1 OUTPUT(1:) = AREA(1:) 700 N = INDEX(OUTPUT,'],[') IF (N.GT.1) THEN NROW = NROW + 1 OUTPUT(1:) = OUTPUT(N+2:) GOTO 700 ENDIF IF (NROW.EQ.1) OPTION(5:5) = 'Y' !that's how it was before... C C if Merge option we first get the no. of pixels involved + then create C an artificial image with all those pixels IF (DEFAUL(1:1).EQ.'M') THEN HISTORY(1:100) = AREA(1:) !save AREA string SIZ(1) = 0 DO 770, NLOOP=1,NROW !loop through subframes N = INDEX(AREA(1:),']') OUTPUT(1:) = AREA(1:N)//' ' CALL EXTCOO(IMNO,OUTPUT,3,IAV,SPIX,EPIX,STAT) SIZ(2) = (EPIX(3) - SPIX(3) + 1) * + (EPIX(2) - SPIX(2) + 1) * + (EPIX(1) - SPIX(1) + 1) SIZ(1) = SIZ(1) + SIZ(2) AREA(1:) = AREA(N+2:)//' ' 770 CONTINUE CALL ARTIMA(0,IMNO,NPIX,AREA,SIZ(1),IMNOC,STAT) AREA(1:) = HISTORY(1:100) DO 800, NLOOP=1,NROW N = INDEX(AREA(1:),']') OUTPUT(1:) = AREA(1:N)//' ' CALL ARTIMA(1,IMNO,NPIX,OUTPUT,0,IMNOC,STAT) CALL STTPUT(OUTPUT,STAT) AREA(1:) = AREA(N+2:)//' ' 800 CONTINUE GOTO 8000 !now proceed as if full frame ENDIF C C loop through all subframes DO 850, NLOOP=1,NROW N = INDEX(AREA(1:),']') OUTPUT(1:) = AREA(1:N)//' ' IF (EXAMED.EQ.1) THEN CALL EXTCOO(IMNO,OUTPUT,3,IAV,SPIX,EPIX,STAT) IF (STAT.NE.0) CALL STETER(34, + 'Invalid syntax in subframe ...') IF (IAV.LT.3) THEN CALL EXMED(IACTIO,IMNO,ZBINS(2),NPIX,SPIX, + EPIX,XMEDIA,STAT) IF (STAT.NE.0) CALL STETER(35, + 'Problems with exact median calculation...') ELSE CALL STETER(27, + 'G- or X-option not supported for 3-d frames, yet') ENDIF ENDIF CALL ZMSTAT(IMNO,OUTPUT,NAXPIX,START,STEP, + ZBINS,FRMSTR,OPTION) IF (OUTID.NE.-999) THEN CALL STKRDR('OUTPUTR',1,12,IAV,RVAL,UNI,NULO,STAT) CALL TBRWRR(OUTID,NLOOP+TBAPP,NCOL1,TCOLNM,RVAL,STAT) ENDIF IF (DEFAUL(6:6).EQ.'P') CALL PLOTI(IMNO) AREA(1:) = AREA(N+2:)//' ' 850 CONTINUE C C --------- C C 3) handle interval mode C C --------- C ELSE IF (IACTIO.EQ.7) THEN N = 1 SPIX(2) = 1 C IF (NAXPIX(1).GE.3) THEN NPIX(3) = 1 !work on 1st plane only C DO 2200, ST2=1,NROW !loop over all intervals SPIX(1) = 1 EPIX(2) = SPIX(2)+SIZ(2)-1 IF (EPIX(2).GT.NPIX(2)) EPIX(2) = NPIX(2) C DO 2100, ST1=1,NMAL EPIX(1) = SPIX(1)+SIZ(1)-1 IF (EPIX(1).GT.NPIX(1)) EPIX(1) = NPIX(1) IF (EXAMED.EQ.1) THEN CALL EXMED(IACTIO,IMNO,ZBINS(2),NPIX,SPIX,EPIX, + XMEDIA,STAT) IF (STAT.NE.0) CALL STETER(35, + 'Problems with exact median calculation...') ENDIF WRITE(AREA,10009) SPIX(1),SPIX(2),EPIX(1),EPIX(2) CALL ZMSTAT(IMNO,AREA,NAXPIX,START, + STEP,ZBINS,FRMSTR,OPTION) IF (OUTID.NE.-999) THEN CALL STKRDR('OUTPUTR',1,12,IAV,RVAL,UNI,NULO,STAT) CALL TBRWRR(OUTID,N+TBAPP,NCOL1,TCOLNM,RVAL,STAT) N = N + 1 ENDIF IF (DEFAUL(6:6).EQ.'P') CALL PLOTI(IMNO) SPIX(1) = SPIX(1) + SIZ(1) 2100 CONTINUE SPIX(2) = SPIX(2) + SIZ(2) 2200 CONTINUE C ELSE IF (NAXPIX(1).EQ.2) THEN DO 2400, ST2=1,NROW !loop over all intervals SPIX(1) = 1 EPIX(2) = SPIX(2)+SIZ(2)-1 IF (EPIX(2).GT.NPIX(2)) EPIX(2) = NPIX(2) C DO 2300, ST1=1,NMAL EPIX(1) = SPIX(1)+SIZ(1)-1 IF (EPIX(1).GT.NPIX(1)) EPIX(1) = NPIX(1) IF (EXAMED.EQ.1) THEN CALL EXMED(IACTIO,IMNO,ZBINS(2),NPIX,SPIX,EPIX, + XMEDIA,STAT) IF (STAT.NE.0) CALL STETER(35, + 'Problems with exact median calculation...') ENDIF WRITE(AREA,10008) SPIX(1),SPIX(2),EPIX(1),EPIX(2) CALL ZMSTAT(IMNO,AREA,NAXPIX,START, + STEP,ZBINS,FRMSTR,OPTION) IF (OUTID.NE.-999) THEN CALL STKRDR('OUTPUTR',1,12,IAV,RVAL,UNI,NULO,STAT) CALL TBRWRR(OUTID,N+TBAPP,NCOL1,TCOLNM,RVAL,STAT) N = N + 1 ENDIF IF (DEFAUL(6:6).EQ.'P') CALL PLOTI(IMNO) SPIX(1) = SPIX(1) + SIZ(1) 2300 CONTINUE SPIX(2) = SPIX(2) + SIZ(2) 2400 CONTINUE C ELSE SPIX(1) = 1 SPIX(2) = 1 EPIX(2) = 1 DO 2500, ST2=1,NMAL !loop over all intervals EPIX(1) = SPIX(1)+SIZ(1)-1 IF (EPIX(1).GT.NPIX(1)) EPIX(1) = NPIX(1) IF (EXAMED.EQ.1) THEN CALL EXMED(IACTIO,IMNO,ZBINS(2),NPIX,SPIX,EPIX, + XMEDIA,STAT) IF (STAT.NE.0) CALL STETER(35, + 'Problems with exact median calculation...') ENDIF WRITE(AREA,10007) SPIX(1),EPIX(1) CALL ZMSTAT(IMNO,AREA,NAXPIX,START, + STEP,ZBINS,FRMSTR,OPTION) IF (OUTID.NE.-999) THEN CALL STKRDR('OUTPUTR',1,12,IAV,RVAL,UNI,NULO,STAT) CALL TBRWRR(OUTID,N+TBAPP,NCOL1,TCOLNM,RVAL,STAT) N = N + 1 ENDIF IF (DEFAUL(6:6).EQ.'P') CALL PLOTI(IMNO) SPIX(1) = SPIX(1) + SIZ(1) 2500 CONTINUE ENDIF N = NMAL * NROW WRITE(OUTPUT,10011) N !display no. of intervals processed CALL STTPUT(OUTPUT,STAT) C C --------- C C 4) handle row modes C C --------- C ELSE IF (IACTIO.EQ.2) THEN CACT(1:) = 'RRRR ' C C if Merge option we first get the no. of pixels involved + then create C an artificial image with all those pixels IF (DEFAUL(1:1).EQ.'M') THEN SIZ(1) = NROW * NPIX(1) CALL ARTIMA(0,IMNO,NPIX,AREA,SIZ(1),IMNOC,STAT) DO 2700, NLOOP=1,NROW IF (KBUF(1).GT.0) THEN N = KBUF(NLOOP) ELSE N = NLOOP ENDIF IF (NPIX(3).GT.1) THEN WRITE(AREA,10004) N,N ELSE WRITE(AREA,10003) N,N ENDIF CALL ARTIMA(1,IMNO,NPIX,AREA,0,IMNOC,STAT) WRITE(OUTPUT,20101) N CALL STTPUT(OUTPUT,STAT) 2700 CONTINUE GOTO 8000 !now proceed as if full frame ENDIF C C loop through all rows DO 2800, NLOOP=1,NROW IF (KBUF(1).GT.0) THEN N = KBUF(NLOOP) ELSE N = NLOOP ENDIF WRITE(AREA,10000) CACT(1:1),N,N IF (EXAMED.EQ.1) THEN SPIX(1) = 1 SPIX(2) = N EPIX(1) = NPIX(1) EPIX(2) = N CALL EXMED(IACTIO,IMNO,ZBINS(2),NPIX,SPIX,EPIX, + XMEDIA,STAT) IF (STAT.NE.0) CALL STETER(35, + 'Problems with exact median calculation...') ENDIF CALL ZMSTAT(IMNO,AREA,NAXPIX,START,STEP, + ZBINS,FRMSTR,OPTION) IF (OUTID.NE.-999) THEN CALL STKRDR('OUTPUTR',1,12,IAV,RVAL,UNI,NULO,STAT) CALL TBRWRR(OUTID,NLOOP+TBAPP,NCOL1,TCOLNM,RVAL,STAT) ENDIF IF (DEFAUL(6:6).EQ.'P') CALL PLOTI(IMNO) 2800 CONTINUE C C C --------- C C 5) handle column modes C C --------- C ELSE IF (IACTIO.EQ.3) THEN CACT(1:) = 'CCCC ' C C transpose frame + switch descriptors SIZ(1) = 1 DO 2900, N=1,NAXPIX(1) SIZ(1) = SIZ(1) * NPIX(N) 2900 CONTINUE CALL STFMAP(IMNO,F_I_MODE,1,SIZ(1),IAV,PNTRA,STAT) CALL STFCRE(TMPFIL,D_R4_FORMAT,F_O_MODE,F_IMA_TYPE, + SIZ(1),IMNOB,STAT) CALL STFMAP(IMNOB,F_O_MODE,1,SIZ(1),IAV,PNTRB,STAT) SIZ(1) = 128 SIZ(2) = 256 CALL LINCOL(MADRID(PNTRA),NPIX,SIZ,MADRID(PNTRB)) CALL STFUNM(IMNO,STAT) CALL STFUNM(IMNOB,STAT) C C we have to write these descriptors to the transposed frame, C because of function Convcoo inside Zstats (ZMSTAT) C CALL STDWRI(IMNOB,'NAXIS',NAXPIX(1),1,1,UNI,STAT) N = NAXPIX(2) NAXPIX(2) = NAXPIX(3) NAXPIX(3) = N CALL STDWRI(IMNOB,'NPIX',NAXPIX(2),1,2,UNI,STAT) DVAL = START(1) START(1) = START(2) START(2) = DVAL CALL STDWRD(IMNOB,'STEP',STEP,1,2,UNI,STAT) DVAL = STEP(1) STEP(1) = STEP(2) STEP(2) = DVAL CALL STDWRD(IMNOB,'START',START,1,2,UNI,STAT) NPIX(1) = NAXPIX(2) !NAXPIX already swapped NPIX(2) = NAXPIX(3) OUTPUT(1:) = CUNIT(16:) !switch CUNIT (needed for Wcoords) CUNIT(16:31) = OUTPUT(32:47) CUNIT(32:47) = OUTPUT(16:31) N = (NAXPIX(1)+1) * 16 CALL STDWRC(IMNOB,'CUNIT',1,CUNIT,1,N,UNI,STAT) C C if Merge option we first get the no. of pixels involved + then create C an artificial image with all those pixels IF (DEFAUL(1:1).EQ.'M') THEN SIZ(1) = NROW * NPIX(1) CALL ARTIMA(0,IMNOB,NPIX,AREA,SIZ(1),IMNOC,STAT) DO 2950, NLOOP=1,NROW IF (KBUF(1).GT.0) THEN N = KBUF(NLOOP) ELSE N = NLOOP ENDIF IF (NPIX(3).GT.1) THEN WRITE(AREA,10004) N,N ELSE WRITE(AREA,10003) N,N ENDIF CALL ARTIMA(1,IMNOB,NPIX,AREA,0,IMNOC,STAT) WRITE(OUTPUT,20102) N CALL STTPUT(OUTPUT,STAT) 2950 CONTINUE GOTO 8000 !now proceed as if full frame ENDIF C C now loop through rows (columns of original frame) DO 3000, NLOOP=1,NROW IF (KBUF(1).GT.0) THEN N = KBUF(NLOOP) ELSE N = NLOOP ENDIF WRITE(AREA,10000) CACT(1:1),N,N IF (EXAMED.EQ.1) THEN SPIX(1) = 1 SPIX(2) = N EPIX(1) = NAXPIX(2) !x-dim EPIX(2) = N CALL EXMED(IACTIO,IMNOB,ZBINS(2),NPIX,SPIX,EPIX, + XMEDIA,STAT) IF (STAT.NE.0) CALL STETER(35, + 'Problems with exact median calculation...') ENDIF CALL ZMSTAT(IMNOB,AREA,NAXPIX,START,STEP, + ZBINS,FRMSTR,OPTION) IF (OUTID.NE.-999) THEN CALL STKRDR('OUTPUTR',1,12,IAV,RVAL,UNI,NULO,STAT) CALL TBRWRR(OUTID,NLOOP+TBAPP,NCOL1,TCOLNM,RVAL,STAT) ENDIF IF (DEFAUL(6:6).EQ.'P') CALL PLOTI(IMNO) 3000 CONTINUE CALL STFDEL(TMPFIL,STAT) C C --------- C C 6) in table mode, we need some more table massaging C C --------- C ELSE IF (IACTIO.EQ.1) THEN C C if Merge option we first have to get the no. of pixels involved IF (DEFAUL(1:1).EQ.'M') THEN SIZ(1) = 0 DO 4100, N=1,NROW CALL TBSGET(INTID,N,SELFLG,STAT) !only use selected rows IF (SELFLG) THEN CALL TBRRDR(INTID,N,INCOL,ICOLNM(1),RVAL,TABNUL,STAT) IF (INCOL.EQ.4) THEN WRITE(AREA,10001) RVAL(1),RVAL(2),RVAL(3),RVAL(4) ELSE WRITE(AREA,10002) + RVAL(1),RVAL(2),RVAL(5),RVAL(3),RVAL(4),RVAL(6) ENDIF CALL EXTCOO(IMNO,AREA,3,IAV,SPIX,EPIX,STAT) SIZ(2) = (EPIX(3) - SPIX(3) + 1) * + (EPIX(2) - SPIX(2) + 1) * + (EPIX(1) - SPIX(1) + 1) SIZ(1) = SIZ(1) + SIZ(2) ENDIF 4100 CONTINUE CALL ARTIMA(0,IMNO,NPIX,AREA,SIZ(1),IMNOC,STAT) C C then create artificial image with all those pixels DO 4200, N=1,NROW CALL TBSGET(INTID,N,SELFLG,STAT) !only use selected rows IF (SELFLG) THEN CALL TBRRDR(INTID,N,INCOL,ICOLNM(1),RVAL,TABNUL,STAT) IF (INCOL.EQ.4) THEN WRITE(AREA,10001) RVAL(1),RVAL(2),RVAL(3),RVAL(4) ELSE WRITE(AREA,10002) + RVAL(1),RVAL(2),RVAL(5),RVAL(3),RVAL(4),RVAL(6) ENDIF CALL BLANKO(AREA) CALL STTPUT(AREA,STAT) CALL ARTIMA(1,IMNO,NPIX,AREA,0,IMNOC,STAT) ENDIF 4200 CONTINUE GOTO 8000 !now proceed as if full frame ENDIF C C now loop through input table rows NMAL = 0 DO 4400, N=1,NROW C CALL TBSGET(INTID,N,SELFLG,STAT) !only use selected rows IF (.NOT.SELFLG) GOTO 4400 C NMAL = NMAL + 1 AREA(1:) = ' ' CALL TBRRDR(INTID,N,INCOL,ICOLNM(1),RVAL,TABNUL,STAT) C IF (INCOL.EQ.4) THEN WRITE(AREA,10001) RVAL(1),RVAL(2),RVAL(3),RVAL(4) ELSE WRITE(AREA,10002) + RVAL(1),RVAL(2),RVAL(5),RVAL(3),RVAL(4),RVAL(6) ENDIF CALL BLANKO(AREA) IF (OPTION(4:4) .NE. 'N') THEN OUTPUT(1:) = ' ' CALL STTPUT(OUTPUT,STAT) ENDIF IF (EXAMED.EQ.1) THEN IF (INCOL.EQ.4) THEN DD1(1) = RVAL(1) !convert wc -> fp DD1(2) = RVAL(2) CALL FPXWCO(-1,0,DD1,DD2,STAT) IF (STAT .GT. 0) + CALL STETER(36,'Problems with WCS of frame...') SPIX(1) = DD2(1) SPIX(2) = DD2(2) DD1(1) = RVAL(3) DD1(2) = RVAL(4) CALL FPXWCO(-1,0,DD1,DD2,STAT) IF (STAT .GT. 0) + CALL STETER(36,'Problems with WCS of frame...') EPIX(1) = DD2(1) EPIX(2) = DD2(2) CALL EXMED(IACTIO,IMNO,ZBINS(2),NPIX,SPIX,EPIX, + XMEDIA,STAT) IF (STAT.NE.0) CALL STETER(35, + 'Problems with exact median calculation...') ELSE CALL STETER(27, + 'G- or X-option not supported for 3-d frames, yet') ENDIF ENDIF CALL ZMSTAT(IMNO,AREA,NAXPIX,START,STEP, + ZBINS,FRMSTR,OPTION) IF (OUTID.NE.-999) THEN CALL STKRDR('OUTPUTR',1,12,IAV,RVAL,UNI,NULO,STAT) CALL TBRWRR(OUTID,N+TBAPP,NCOL1,TCOLNM,RVAL,STAT) ENDIF IF (DEFAUL(6:6).EQ.'P') CALL PLOTI(IMNO) 4400 CONTINUE C CALL TBTCLO(INTID,STAT) WRITE(OUTPUT,10010) NMAL !display no. of rows processed CALL STTPUT(OUTPUT,STAT) C C --------- C C 7) handle plane modes C C --------- C ELSE IF (IACTIO.EQ.5) THEN CACT(1:) = 'PPPP ' C C if Merge option we first get the no. of pixels involved + then create C an artificial image with all those pixels IF (DEFAUL(1:1).EQ.'M') THEN SIZ(1) = NROW * NPIX(1)*NPIX(2) CALL ARTIMA(0,IMNO,NPIX,AREA,SIZ(1),IMNOC,STAT) DO 4700, NLOOP=1,NROW IF (KBUF(1).GT.0) THEN N = KBUF(NLOOP) ELSE N = NLOOP ENDIF WRITE(AREA,10005) N,N CALL ARTIMA(1,IMNO,NPIX,AREA,0,IMNOC,STAT) WRITE(OUTPUT,20100) N CALL STTPUT(OUTPUT,STAT) 4700 CONTINUE GOTO 8000 !now proceed as if full frame ENDIF C C loop through planes DO 5000, NLOOP=1,NROW IF (KBUF(1).GT.0) THEN N = KBUF(NLOOP) ELSE N = NLOOP ENDIF AREA(1:) = ' ' WRITE(AREA,20000) CACT(1:1),N,N IF (EXAMED.EQ.1) THEN SPIX(1) = 1 SPIX(2) = 1 SPIX(3) = N EPIX(1) = NPIX(1) EPIX(2) = NPIX(2) EPIX(3) = N CALL EXMED(IACTIO,IMNO,ZBINS(2),NPIX,SPIX,EPIX, + XMEDIA,STAT) IF (STAT.NE.0) CALL STETER(35, + 'Problems with exact median calculation...') ENDIF CALL ZMSTAT(IMNO,AREA,NAXPIX,START,STEP, + ZBINS,FRMSTR,OPTION) IF (OUTID.NE.-999) THEN CALL STKRDR('OUTPUTR',1,12,IAV,RVAL,UNI,NULO,STAT) CALL TBRWRR(OUTID,NLOOP+TBAPP,NCOL1,TCOLNM,RVAL,STAT) ENDIF IF (DEFAUL(6:6).EQ.'P') CALL PLOTI(IMNO) 5000 CONTINUE C C --------- C C 8) finally process cursor mode C C --------- C ELSE CALL DTOPEN(1,STAT) CALL SETCUR(QDSPNO,2,0,0,COOS,STAT) !set up cursor window CACT = 'NNYY?C0 ' IF (DEFAUL(5:5).NE.'N') CACT(2:2) = 'Y' ORFILE(1:) = ' ' COOFF = 0 C C if Merge option we create the art. image with 100 000 pixels first C we may have to extend the image later on... IF (DEFAUL(1:1).EQ.'M') THEN SIZ(1) = 100000 CALL ARTIMA(0,IMNO,NPIX,AREA,SIZ(1),IMNOC,STAT) ENDIF C C now loop through cursor inputs DO 6000, N=1,CURMAX 5800 CALL GETCUR(CACT,ORFILE, + XY1,RVAL(1),RVAL(13),RVAL(10),ST1, + XY2,RVAL(3),RVAL(15),RVAL(11),ST2) C IF ((ST1.EQ.0).AND.(ST2.EQ.0)) THEN IF ((N.EQ.1).AND.(COOFF.EQ.0)) THEN COOFF = 1 ORFILE(1:) = ' ' CALL STTDIS + ('switch cursors on - next time we exit...',99,STAT) GOTO 5800 ELSE GOTO 6600 !we're finished ENDIF ENDIF C AREA(1:) = ' ' WRITE(AREA,10001) RVAL(13),RVAL(14),RVAL(15),RVAL(16) CALL BLANKO(AREA) IF (DEFAUL(1:1).EQ.'M') THEN CALL STTPUT(AREA,STAT) 5880 CALL ARTIMA(1,IMNO,NPIX,AREA,0,IMNOC,STAT) IF (STAT .EQ. -4) THEN !expand art. image SIZ(1) = SIZ(1) + SIZ(1) CALL ARTIMA(2,IMNO,NPIX,OUTPUT,SIZ(1),IMNOC,STAT) GOTO 5880 !redo last filling ENDIF ELSE IF (EXAMED.EQ.1) THEN SPIX(1) = NINT(RVAL(1)) SPIX(2) = NINT(RVAL(2)) EPIX(1) = NINT(RVAL(3)) EPIX(2) = NINT(RVAL(4)) CALL EXMED(IACTIO,IMNO,ZBINS(2),NPIX, + SPIX,EPIX,XMEDIA,STAT) IF (STAT.NE.0) CALL STETER(35, + 'Problems with exact median calculation...') ENDIF CALL ZMSTAT(IMNO,AREA,NAXPIX,START,STEP, + ZBINS,FRMSTR,OPTION) IF (OUTID.NE.-999) THEN CALL STKRDR('OUTPUTR',1,12,IAV,RVAL,UNI,NULO,STAT) NROW = N+TBAPP CALL TBRWRR(OUTID,NROW,NCOL1,TCOLNM,RVAL,STAT) IAV = NCOL1 + 1 !first column after statistic stuff CALL TBRWRR(OUTID,NROW,4,TCOLNM(IAV),RVAL(13),STAT) ENDIF IF (DEFAUL(6:6).EQ.'P') CALL PLOTI(IMNO) ENDIF 6000 CONTINUE C 6600 CALL DTCLOS(QDSPNO) IF (DEFAUL(1:1).EQ.'M') THEN !get exact size CALL ARTIMA(3,IMNO,NPIX,AREA,SIZ(1),IMNOC,STAT) GOTO 8000 !do the statistics now ENDIF C ENDIF GOTO 9000 C C --------- C C here from Merge option - proceed as if full frame C C --------- C 8000 NPIX(1) = SIZ(1) NPIX(2) = 1 NPIX(3) = 1 IF (EXAMED.EQ.1) THEN CALL EXMED(0,IMNOC,ZBINS(2),NPIX,SPIX,EPIX,XMEDIA,STAT) IF (STAT.NE.0) + CALL STETER(35,'Problems with exact median calculation...') ENDIF NAXPIX(1) = 1 NAXPIX(2) = NPIX(1) NAXPIX(3) = 1 NAXPIX(4) = 1 OPTION(5:5) = 'N' !don't try to write descr. START(1) = 0.D0 STEP(1) = 1.D0 CALL ZMSTAT(IMNOC,AREA,NAXPIX,START,STEP,ZBINS,FRMSTR,OPTION) IF (OUTID.NE.-999) THEN N = 1 CALL STKRDR('OUTPUTR',1,12,IAV,RVAL,UNI,NULO,STAT) CALL TBRWRR(OUTID,N+TBAPP,NCOL1,TCOLNM,RVAL,STAT) ENDIF IF (DEFAUL(6:6).EQ.'P') CALL PLOTI(IMNO) C C C *********************** C C everything joins here C C *********************** C C 9000 IF ((OUTID.NE.-999) .AND. (OUTID.NE.INTID)) THEN IF (OUTNO .GT. 0) THEN !copy temp table to outimage CALL TBIGET(OUTID,NPXO(1),NPXO(2),N,N,N,STAT) STRTO(1) = 0.0 STPO(1) = 1.0 IF (OUTNO.EQ.999) THEN STRTO(2) = 0.0 STPO(2) = 1.0 NAXPIX(1) = 2 !2dim result image N = 1 ELSE NAXPIX(1) = 1 !1dim result image N = 2 ENDIF CALL STIPUT(OUTFRA,D_R4_FORMAT,F_O_MODE,F_IMA_TYPE, + NAXPIX(1),NPXO(N),STRTO,STPO, + IDENT,CUNIT,PNTRO,IMNOB,STAT) CALL TB2IMA(OUTNO,OUTID,NPXO,MADRID(PNTRO)) !table -> image CALL STDWRC(IMNOB,'HISTORY',1,HISTORY,1,80,UNI,STAT) IF (IACTIO.EQ.7) THEN !save xsize, ysize of intv. ICOLNM(1) = SIZ(1) ICOLNM(2) = SIZ(2) ICOLNM(3) = NMAL ICOLNM(4) = NROW CALL STDWRI(IMNOB,'INTVALS',ICOLNM,1,4,UNI,STAT) CALL STECNT('GET',EC,EL,ED) CALL STECNT('PUT',1,0,0) !disable errors ... CALL STKRDC('MYCOMND',1,1,80,IAV,HISTORY,UNI,NULO,STAT) IF (STAT.EQ.0) + CALL STDWRC(IMNOB,'HISTORY',1,HISTORY,1,80,UNI,STAT) CALL STECNT('PUT',EC,EL,ED) ENDIF CALL STFCLO(IMNOB,STAT) ELSE CALL TBSINI(OUTID,STAT) !init selection flags CALL TBTCLO(OUTID,STAT) !close table correctly ENDIF ENDIF CALL STSEPI C 10000 FORMAT(A,'[<,@',I5.5,':>,@',I5.5,'] ') 10001 FORMAT('[',G12.6,',',G12.6,':',G12.6,',',G12.6,'] ') 10002 FORMAT +('[',G12.6,',',G12.6,',',G12.6,':',G12.6,',',G12.6,',',G12.6,']') 10003 FORMAT('[<,@',I5.5,':>,@',I5.5,'] ') 10004 FORMAT('[<,@',I5.5,',<:>,@',I5.5,',<] ') 10005 FORMAT('[<,<,@',I5.5,':>,>,@',I5.5,'] ') 10007 FORMAT('[@',I7.7,':@',I7.7,'] ') 10008 FORMAT('[@',I7.7,',@',I7.7,':@',I7.7,',@',I7.7,'] ') 10009 FORMAT('[@',I7.7,',@',I7.7,',<:@',I7.7,'@',I7.7,',<] ') 10010 FORMAT(I7,' table entries processed.') 10011 FORMAT(I7,' intervals processed.') 10020 FORMAT('column :',A,'is missing in output table ',A) 10030 FORMAT('column :',A,'is missing in input table ',A) 20000 FORMAT(A,'[<,<,@',I5.5,':>,>,@',I5.5,'] ') 20010 FORMAT('from STATIST/IMAGE ',A,' ... ',A) 20100 FORMAT('plane ',I4) 20101 FORMAT('row ',I6) 20102 FORMAT('column ',I6) C END C C copy all or just on column of the result table C to a 2/1-dim image C SUBROUTINE TB2IMA(OUTNO,TID,NCRO,A) C IMPLICIT NONE C INTEGER OUTNO,TID,NCRO(*) INTEGER INDX,N,NCOL,NROW,STAT INTEGER ICOLNM(12) !max no. of columns INTEGER TABNUL(12) C REAL A(*),RBUF(12) C NCOL = NCRO(1) NROW = NCRO(2) DO 500, N=1,NCOL ICOLNM(N) = N 500 CONTINUE C IF (OUTNO.EQ.999) THEN INDX = 1 DO 1000, N=1,NROW CALL TBRRDR(TID,N,NCOL,ICOLNM(1),A(INDX),TABNUL,STAT) INDX = INDX + NCOL 1000 CONTINUE C ELSE DO 1500, N=1,12 RBUF(N) = 0.0 1500 CONTINUE INDX = 1 DO 2000, N=1,NROW CALL TBRRDR(TID,N,NCOL,ICOLNM(1),RBUF,TABNUL,STAT) A(INDX) = RBUF(OUTNO) INDX = INDX + 1 2000 CONTINUE ENDIF C RETURN END SUBROUTINE EXMED(ACTIO,INO,ZBINS,NPIX,SPIX,EPIX,VALMED,STAT) C INTEGER ACTIO,INO,NPIX(*),SPIX(*),EPIX(*),STAT INTEGER IAV,N,CHUNKA,CHUNKB,NN INTEGER NA,NB,INSIZ,OUTSIZ,UNI(1),IBUF(7) INTEGER*8 PNTRX,PNTRY INTEGER IMNOX,IMNOY,MADRID(1) C REAL ZBINS(*),VALMED REAL RBUF(7) C INCLUDE 'MID_INCLUDE:ST_DEF.INC' COMMON /VMR/ MADRID INCLUDE 'MID_INCLUDE:ST_DAT.INC' C C plane option C IF (ACTIO.EQ.5) THEN INSIZ = NPIX(1) * NPIX(2) OUTSIZ = INSIZ NA = ((SPIX(3)-1)*NPIX(1)*NPIX(2)) + 1 CALL STFCRE('middstat2',D_R4_FORMAT,F_X_MODE,F_IMA_TYPE, + OUTSIZ,IMNOY,STAT) CALL STFMAP(IMNOY,F_X_MODE,1,OUTSIZ,IAV,PNTRY,STAT) CALL STFGET(INO,NA,INSIZ,IAV,MADRID(PNTRY),STAT) C C subframe input C ELSE IF (ACTIO.NE.0) THEN NN = EPIX(2)-SPIX(2)+1 INSIZ = NPIX(1) * NN !size of full Y-strip CHUNKA = NPIX(1) CHUNKB = EPIX(1) - SPIX(1) + 1 OUTSIZ = CHUNKB * NN NA = (SPIX(2)-1)*NPIX(1) + 1 C CALL STFCRE('middstat1',D_R4_FORMAT,F_X_MODE,F_IMA_TYPE, + CHUNKA,IMNOX,STAT) CALL STFMAP(IMNOX,F_X_MODE,1,CHUNKA,IAV,PNTRX,STAT) CALL STFCRE('middstat2',D_R4_FORMAT,F_X_MODE,F_IMA_TYPE, + OUTSIZ,IMNOY,STAT) CALL STFMAP(IMNOY,F_X_MODE,1,OUTSIZ,IAV,PNTRY,STAT) NB = 1 DO 500,N=1,NN CALL STFGET(INO,NA,CHUNKA,IAV,MADRID(PNTRX),STAT) CALL DATMOV(MADRID(PNTRX),NB,CHUNKB,SPIX(1),MADRID(PNTRY)) NA = NA + CHUNKA NB = NB + CHUNKB 500 CONTINUE CALL STFCLO(IMNOX,STAT) C C full frame C ELSE INSIZ = NPIX(1) * NPIX(2) * NPIX(3) OUTSIZ = INSIZ CALL STFCRE('middstat2',D_R4_FORMAT,F_X_MODE,F_IMA_TYPE, + OUTSIZ,IMNOY,STAT) CALL STFMAP(IMNOY,F_X_MODE,1,OUTSIZ,IAV,PNTRY,STAT) CALL STFGET(INO,1,INSIZ,IAV,MADRID(PNTRY),STAT) ENDIF C C sort array + return median C NA = (OUTSIZ+1)/2 !index of median CALL SORTME(MADRID(PNTRY),ZBINS,OUTSIZ,NA,VALMED,STAT) C IF (STAT.EQ.-1) THEN CALL STTPUT('No pixels found with data in given interval...',N) DO 900, N=1,7 RBUF(N) = -1.0 IBUF(N) = -1 900 CONTINUE CALL STKWRR('OUTPUTR',RBUF,1,7,UNI,N) CALL STKWRI('OUTPUTI',IBUF,1,7,UNI,N) CALL STFCLO(IMNOY,N) !keep value of STAT ELSE CALL STKWRR('OUTPUTR',VALMED,8,1,UNI,STAT) CALL STFCLO(IMNOY,STAT) ENDIF C RETURN END SUBROUTINE PLOTI(IMNO) C IMPLICIT NONE C INTEGER IMNO INTEGER UNI(1),NULO,STAT,IAV C REAL RV C CALL STKRDR('OUTPUTR',10,1,IAV,RV,UNI,NULO,STAT) IF (RV.GT.0.99) THEN CALL PLOHI(IMNO) !only if we got a histogram ELSE CALL STTPUT('No histogram calculated => no plotting...',STAT) ENDIF C RETURN END SUBROUTINE DATMOV(A,OFF,CHUNK,SPIX,Y) C IMPLICIT NONE C INTEGER OFF,CHUNK,SPIX C REAL A(*),Y(*) C INTEGER N,NN C NN = OFF DO 500,N=SPIX,SPIX+CHUNK-1 Y(NN) = A(N) NN = NN + 1 500 CONTINUE C RETURN END