C @(#)ccdscan.for 17.1.1.1 (ESO-IPG) 01/25/02 17:49:37 C $Id: ccdscan.for,v 1.2 1998/08/04 11:26:00 swolf Exp $ 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 Massachusetss Ave, Cambridge, C MA 02139, USA. C C Corresponding 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 C $Id: ccdscan.for,v 1.2 1998/08/04 11:26:00 swolf Exp $ PROGRAM OVSCAN C -------------------------------------------------------------------- C.COPYRIGHT: Copyright (c) 1983 European Southern Observatory, C all rights reserved C.IDENTIFICATION: CCDOVSCAN.FOR C.LANGUAGE: F77+ESOext C.AUTHOR: R.H. Warmels C.KEYWORDS: average, overscan C.PURPOSE: produce a 1d image from a 2d by averaging over rows or columns C.ALGORITHM add rows/columns C.VERSION: 900202 KB take care of options in lower case C.VERSION: ?????? MP cosmectic changes C.VERSION: 910417 xx replace numbers by symbolic constants in ST calls C.VERSION: 910417 xx replace numbers by symbolic constants in ST calls C.VERSION: 930226 RHW modification for ccd package C.VERSION: 930331 RHW modification for average row wise C.VERSION: 980609 SW correct use of the selected AREA (OV_SEC) C ------------------------------------------------------------------ IMPLICIT NONE C INTEGER MADRID(1) INTEGER NDIM,KUN,KNUL INTEGER I1,I2,L1,L2,NI,NL INTEGER NAXIS,IMNO INTEGER NPIX(3) INTEGER*8 PNTRA,PNTRB,PNTRC INTEGER TIDO INTEGER NCOL, NROW, COL(3) INTEGER NDUM INTEGER STATUS INTEGER IMAGE(4) INTEGER IAV, ILOG INTEGER I, ISMO C variable in first part not used anymore 05/00 PN C INTEGER IR INTEGER IXMIN, IXMAX INTEGER NCUR, KA, NPAR INTEGER ACCESS, PLMODE INTEGER NBINS, NACT, NNN, IX INTEGER IFREQ(512) C REAL AMIN, AMAX, XMIN, XMAX, STP REAL X0, X1, XT REAL CINT(512), CR REAL LHCUTS(4) REAL FRAME(8) REAL XPA(2), YPA(2) REAL XMEAN, XSTD REAL XCUR, YCUR REAL SCALES(3),OFFSET(2),CLIP(4) DOUBLE PRECISION STEP(3) DOUBLE PRECISION START(3) C CHARACTER*60 FRAM, AREA CHARACTER*60 TABLE CHARACTER*80 CINPUT CHARACTER*80 TEXT, SEL CHARACTER*17 COLUMN(3) CHARACTER*80 LABEL1, LABEL2, LABEL3 CHARACTER*72 CUNIT CHARACTER*1 DIRECT CHARACTER*72 IDENT CHARACTER IOP*1 CHARACTER TABUNI(3)*16,TABLAB(3)*16,TABFOR(3)*16 C INCLUDE 'MID_INCLUDE:ST_DEF.INC' COMMON /VMR/ MADRID INCLUDE 'MID_INCLUDE:ST_DAT.INC' C DATA SEL /' '/ DATA ILOG /0/ DATA COLUMN /':X', ':FLUX',':FIT'/ DATA NDIM /2/ DATA IOP /'A'/ DATA TABUNI /' ', ' ', ' '/ DATA TABLAB /'X', 'FLUX', 'FIT'/ DATA TABFOR /'G13.6', 'G13.6','G13.6'/ DATA SCALES /0.0,0.0,0.0/ DATA CLIP /0.15,0.95,0.1,0.95/ DATA OFFSET /-999., -999./ C 9001 FORMAT('*** WARNING: zero dynamics range in y; '// 2 'value = ',G13.6) 9002 FORMAT('[@', I5.5, ',@', I5.5, ':@', I5.5, ',@', I5.5, ']') 9003 FORMAT(':X.GE.',I5.5,'.AND.:X.LE.',I5.5) C C *** initialize system CALL STSPRO('OVSCAN') C C *** read parameter input CALL STKRDC('CCDIN',1,1,60,IAV,FRAM,KUN,KNUL,STATUS) CALL STKRDC('CCDOUT',1,1,60,IAV,TABLE,KUN,KNUL,STATUS) CALL STKRDC('OV_SEC',1,1,60,IAV,AREA,KUN,KNUL,STATUS) CALL STKRDC('DIRECT',1,1,1,IAV,DIRECT,KUN,KNUL,STATUS) ! read direct. CALL STKRDC('OV_IMODE',1,1,3,IAV,CINPUT,KUN,KNUL,STATUS) CALL STKRDI('INPUTI',1,1,IAV,ISMO,KUN,KNUL,STATUS) C C *** read input frame CALL STIGET(FRAM,D_R4_FORMAT,F_I_MODE,F_IMA_TYPE,NDIM,NAXIS, + NPIX,START,STEP,IDENT,CUNIT,PNTRA,IMNO,STATUS) C CALL EXTCOO(IMNO,AREA,NAXIS,NDUM,IMAGE(1),IMAGE(3),STATUS) !pixel coords. IF (STATUS.NE.0) THEN TEXT = '*** FATAL: invalid coordinate input ...' CALL STTPUT(TEXT,STATUS) CALL STSEPI ENDIF I1 = IMAGE(1) I2 = IMAGE(3) L1 = IMAGE(2) L2 = IMAGE(4) C type *,"IMAGE: ",IMAGE NI = I2-I1+1 NL = L2-L1+1 C C *** map output table and fill the column NCOL = 3 CALL UPCAS(DIRECT,DIRECT) IF ((DIRECT.EQ.'R') .OR. (DIRECT.EQ.'X'))THEN DIRECT = 'X' FRAME(1) = FLOAT(L1) FRAME(2) = FLOAT(L2) c NROW = NPIX(2) NROW = NL ELSE IF ((DIRECT.EQ.'C') .OR. (DIRECT.EQ.'Y')) THEN DIRECT = 'Y' FRAME(1) = FLOAT(I1) FRAME(2) = FLOAT(I2) c NROW = NPIX(1) NROW = NI ENDIF CALL TBTINI(TABLE,F_TRANS,F_IO_MODE,NCOL,NROW,TIDO,STATUS)! create table DO 20 IAV = 1,NCOL ! define columns CALL TBCINI(TIDO,D_R4_FORMAT,1,TABFOR(IAV),TABUNI(IAV), 2 TABLAB(IAV),COL(IAV),STATUS) 20 CONTINUE C CALL TBCSER(TIDO,COLUMN(1),COL(1),STATUS) ! find pixel column CALL TBCSER(TIDO,COLUMN(2),COL(2),STATUS) ! find flux column CALL TBCMAP(TIDO,COL(1),PNTRB,STATUS) ! map x column CALL TBCMAP(TIDO,COL(2),PNTRC,STATUS) ! map flux column IF (DIRECT.EQ.'X') THEN ! column readout CALL AVERC(NPIX(1),NPIX(2),MADRID(PNTRA),MADRID(PNTRB), + MADRID(PNTRC),I1,I2,IOP,ISMO,LHCUTS(3),LHCUTS(4)) CALL TBIPUT(TIDO,NCOL,NROW,STATUS) ! put the columns c 98.06.09-SW (swolf@eso.org): c The following selection part is not necessary c as the table consits only of valid rows! c DO 21 IR = 1,L1 c CALL TBSPUT(TIDO,IR,.FALSE.,STATUS) c 21 CONTINUE c DO 22 IR = L2,NROW c CALL TBSPUT(TIDO,IR,.FALSE.,STATUS) c 22 CONTINUE C C *** get the statistics CALL TDMXRS(TIDO,COL(2),NROW,0,XMIN,XMAX) IF (XMIN.EQ.XMAX) THEN XMAX = XMIN + 1. END IF CALL TDSCAL(XMIN,XMAX,1.,X0,X1,IX,XT,NNN) AMIN = X0*10.**IX AMAX = X1*10.**IX STP = XT NBINS = (AMAX-AMIN)/STP + 2 NBINS = MIN(NBINS,512) CR = AMIN CINT(1) = CR DO 23 I = 1,NBINS - 2 CR = CR + STP CINT(I+1) = CR 23 CONTINUE CALL TDRSTA(TIDO,COL(2),NROW,NBINS,1,CINT, 2 IFREQ,NACT,XMEAN,XSTD) ELSE ! row readout CALL AVERR(NPIX(1),NPIX(2),MADRID(PNTRA),MADRID(PNTRB), + MADRID(PNTRC),L1,L2,IOP,ISMO,LHCUTS(3),LHCUTS(4)) CALL TBIPUT(TIDO,NCOL,NROW,STATUS) ! put the columns c 98.06.09-SW (swolf@eso.org): c The following selection part is not necessary c as the table consits only of valid rows! c DO 31 IR = 1,I1 c CALL TBSPUT(TIDO,IR,.FALSE.,STATUS) c 31 CONTINUE c DO 32 IR = I2,NROW c CALL TBSPUT(TIDO,IR,.FALSE.,STATUS) c 32 CONTINUE C C *** get the statistics CALL TDMXRS(TIDO,COL(2),NROW,0,XMIN,XMAX) IF (XMIN.EQ.XMAX) THEN XMAX = XMIN + 1. END IF CALL TDSCAL(XMIN,XMAX,1.,X0,X1,IX,XT,NNN) AMIN = X0*10.**IX AMAX = X1*10.**IX STP = XT NBINS = (AMAX-AMIN)/STP + 2 NBINS = MIN(NBINS,512) CR = AMIN CINT(1) = CR DO 33 I = 1,NBINS - 2 CR = CR + STP CINT(I+1) = CR 33 CONTINUE CALL TDRSTA(TIDO,COL(2),NROW,NBINS,1,CINT, 2 IFREQ,NACT,XMEAN,XSTD) ENDIF LHCUTS(3) = XMEAN-3.*XSTD LHCUTS(4) = XMEAN+3.*XSTD C C *** the interactive mode CALL UPCAS(CINPUT,CINPUT) IF (CINPUT(1:1).EQ.'Y') THEN IF (LHCUTS(3).EQ.LHCUTS(4)) THEN WRITE(TEXT,9001) LHCUTS(3) CALL STTPUT(TEXT,STATUS) FRAME(5) = 0.95*LHCUTS(3) FRAME(6) = 1.05*LHCUTS(4) ELSE FRAME(5) = LHCUTS(3) FRAME(6) = LHCUTS(4) ENDIF C C *** get the frame ranges in x and y CALL GETAXS('AUTO',FRAME(1)) CALL GETAXS('AUTO',FRAME(5)) CALL PTKWRR('XWNDL',4,FRAME(1)) CALL PTKWRR('YWNDL',4,FRAME(5)) C C *** statistics determined; start the display CALL PTKWRR('SCALE',2,SCALES) CALL PTKWRR('OFFS',2,OFFSET) C ACCESS = 0 PLMODE = 1 CALL PTOPEN(' ','overscan',ACCESS,PLMODE) ! open plot file CALL PTKWRR('CLPL',4,CLIP) C IF (DIRECT.EQ.'X') THEN LABEL1 = 'Pixel Row Number ' LABEL2 = 'Pixel Intensity' ELSE LABEL1 = 'Pixel Column Number ' LABEL2 = 'Pixel Intensity' ENDIF LABEL3 = 'TITLE=Frame: '//TABLE CALL PTKWRR('CLPL',4,CLIP) CALL PTAXES(FRAME(1),FRAME(5),LABEL1,LABEL2,LABEL3) ! plot the frame NPAR = 3 CALL PLTBL(TIDO,NPAR,COL(1),COL(2),NROW) ! plot the data C C *** Now do the looping for 2 cursers NCUR = 2 DO 50 I = 1,NCUR CALL PTGCUR(XCUR,YCUR,KA,STATUS) IF (XCUR.LT.1.0) THEN XCUR = 1.0 ENDIF IF (XCUR.GT.FLOAT(NROW)) THEN XCUR = FLOAT(NROW) ENDIF XPA(I) = XCUR YPA(I) = YCUR 50 CONTINUE CALL SORT2(XPA,YPA,NCUR) CALL AGGPLM(XPA,YPA,2,4) IXMIN = INT(XPA(1)) IXMAX = INT(XPA(2)) CALL PTCLOS() WRITE(SEL,9003) IXMIN,IXMAX CALL STKWRC('SCANSEL',1,SEL,1,80,KUN,STATUS) ! set selection flag ENDIF C C *** free data CALL TBTCLO(TIDO,STATUS) CALL STSEPI END SUBROUTINE AVERR(NPIX1,NPIX2,Z,X,Y,L1,L2,IOP,ISM,YMIN,YMAX) C C row average from row I1 to I2 C IMPLICIT NONE INTEGER NPIX1, NPIX2 REAL Z(NPIX1,NPIX2) REAL X(NPIX1),Y(NPIX1) INTEGER L1, L2 CHARACTER*1 IOP INTEGER ISM REAL YMIN, YMAX C INTEGER NPMAX PARAMETER (NPMAX=10000) INTEGER I, L INTEGER IX, ISX, ILIS, IOFF, IXG INTEGER NAV REAL RN REAL DATA(NPMAX) REAL AVE C DO 10 I = 1,NPIX1 X(I) = I Y(I) = 0. 10 CONTINUE IF ((IOP.EQ.'S').OR.(IOP.EQ.'s')) THEN DO 30 I = 1,NPIX1 DO 20 L = L1,L2 DATA(I) = DATA(I) + Z(I,L) 20 CONTINUE 30 CONTINUE ELSE RN = L2 - L1 + 1 DO 50 I = 1,NPIX1 DO 40 L = L1,L2 DATA(I) = DATA(I) + Z(I,L)/RN 40 CONTINUE 50 CONTINUE END IF C C *** do the smooting IF (ISM.GT.1) THEN ILIS = 2*(ISM/2) + 1 IOFF = ILIS/2 DO 60 IX = 1, NPIX1 AVE = 0.0 NAV = 0 DO 61 ISX = 1,ILIS IXG = IX + ISX - IOFF - 1 IF (IXG.LE.0 .OR. IXG.GT.NPIX1) THEN GO TO 61 ENDIF NAV = NAV + 1 AVE = AVE + DATA(IXG) 61 CONTINUE Y(IX) = AVE/FLOAT(NAV) 60 CONTINUE ELSE DO 70 IX = 1, NPIX1 Y(IX) = DATA(IX) 70 CONTINUE ENDIF C YMAX = -1.E30 YMIN = -YMAX DO 80 I = 1,NPIX1 YMIN = MIN(YMIN,Y(I)) YMAX = MAX(YMAX,Y(I)) 80 CONTINUE RETURN END SUBROUTINE AVERC(NPIX1,NPIX2,Z,X,Y,I1,I2,IOP,ISM,YMIN,YMAX) C C Column average from column 1 to column 2 C IMPLICIT NONE INTEGER NPIX1, NPIX2 REAL Z(NPIX1,NPIX2) REAL X(NPIX2), Y(NPIX2) INTEGER I1, I2 CHARACTER*1 IOP INTEGER ISM REAL YMIN, YMAX C INTEGER NPMAX PARAMETER (NPMAX=10000) INTEGER I, J INTEGER IX, ISX, ILIS, IOFF, IXG INTEGER NAV REAL RN REAL DATA(NPMAX) REAL AVE C C *** do the integration DO 10 I = 1,NPIX2 X(I) = I Y(I) = 0. 10 CONTINUE IF ((IOP.EQ.'S').OR.(IOP.EQ.'s')) THEN DO 30 I = 1,NPIX2 DO 20 J = I1,I2 DATA(I) = DATA(I) + Z(J,I) 20 CONTINUE 30 CONTINUE ELSE RN = I2 - I1 + 1 DO 50 I = 1,NPIX2 DO 40 J = I1,I2 DATA(I) = DATA(I) + Z(J,I)/RN 40 CONTINUE 50 CONTINUE END IF C C *** do the smooting IF (ISM.GT.1) THEN ILIS = 2*(ISM/2) + 1 IOFF = ILIS/2 DO 60 IX = 1, NPIX2 AVE = 0.0 NAV = 0 DO 61 ISX = 1,ILIS IXG = IX + ISX - IOFF - 1 IF (IXG.LE.0 .OR. IXG.GT.NPIX2) THEN GO TO 61 ENDIF NAV = NAV + 1 AVE = AVE + DATA(IXG) 61 CONTINUE Y(IX) = AVE/FLOAT(NAV) 60 CONTINUE ELSE DO 70 IX = 1, NPIX2 Y(IX) = DATA(IX) 70 CONTINUE ENDIF C C *** get the minimum and maximum YMAX = -1.E30 YMIN = -YMAX DO 80 I = 1,NPIX2 YMIN = MIN(YMIN,Y(I)) YMAX = MAX(YMAX,Y(I)) 80 CONTINUE RETURN END SUBROUTINE SORT2(RA,RB,N) C+++ C.IDENTIFICATION C.AUTHOR: M.Rosa, ESO-Garching, 831105 C.PURPOSE: Sort string of values into ascending sequence of values C.ALGORITM: The art of scientific computing C Press et al., Cambridge university Press, 1986 C.USE: CALL SORT2(XC,YC,N) C.VERSION: 910115 RHW IMPLICIT NONE added C--- C IMPLICIT NONE INTEGER N REAL RA(N) REAL RB(N) C INTEGER I, IR INTEGER J INTEGER L REAL RRA, RRB L = N/2 + 1 IR = N 10 CONTINUE IF (L.GT.1) THEN L = L - 1 RRA = RA(L) RRB = RB(L) ELSE RRA = RA(IR) RRB = RB(IR) RA(IR) = RA(1) RB(IR) = RB(1) IR = IR - 1 IF (IR.EQ.1) THEN RA(1) = RRA RB(1) = RRB RETURN ENDIF ENDIF I = L J = L+L 20 CONTINUE IF (J.LE.IR) THEN IF (J.LT.IR) THEN IF (RA(J).LT.RA(J+1)) THEN J = J + 1 ENDIF ENDIF IF (RRA.LT.RA(J)) THEN RA(I) = RA(J) RB(I) = RB(J) I = J J = J + J ELSE J = IR + 1 ENDIF GO TO 20 ENDIF RA(I) = RRA RB(I) = RRB GO TO 10 END SUBROUTINE PLTBL(TID,NCOLUM,I1,I2,NROW) C +++ C.PURPOSE: Plots one or two columns of a table C.AUTHOR: Rein H. Warmels C.COMMENTS: none C.VERSION: 931103 RHW Createsd for environment document C --- IMPLICIT NONE INTEGER TID ! Table identifier INTEGER NCOLUM ! # of columns to be plotted INTEGER I1 ! index to first column INTEGER I2 ! index to second column INTEGER NROW ! number of row C INTEGER NPMAX PARAMETER (NPMAX=100000) INTEGER IFIRST, I, ISTAT, IR, IAC INTEGER STYPE, LTYPE REAL VX, VY REAL XPS(NPMAX), YPS(NPMAX) LOGICAL IPLOT,ISEL,INULL CHARACTER TEXT*80 C 9001 FORMAT('*** FATAL: Maximum number of table entries is ',I8) C IFIRST = 0 IR = 0 C DO 10 I = 1,NROW IPLOT = .TRUE. CALL TBSGET(TID,I,ISEL,ISTAT) IF (ISEL) THEN CALL TBRRDR(TID,I,1,I1,VX,INULL,ISTAT) IF (INULL) THEN IPLOT = .FALSE. ENDIF C CALL TBRRDR(TID,I,1,I2,VY,INULL,ISTAT) IF (INULL) THEN IPLOT = .FALSE. ENDIF IF (IPLOT) THEN IR = IR + 1 IF (IR.GT.NPMAX) THEN WRITE (TEXT,9001) NPMAX CALL STTPUT(TEXT,ISTAT) CALL STSEPI ENDIF XPS(IR) = VX YPS(IR) = VY ENDIF ENDIF 10 CONTINUE C CALL PTKRDI('STYPE',1,IAC,STYPE) CALL PTKRDI('LTYPE',1,IAC,LTYPE) IF (LTYPE.EQ.0 .AND. STYPE.EQ.0) THEN CALL STTPUT('*** FATAL: LTYPE and STYPE '// 2 'both equal 0: NO PLOT',ISTAT) CALL PTCLOS() CALL STSEPI ELSE CALL PTDATA(STYPE,LTYPE,0,XPS,YPS,0.0,IR) ENDIF C C *** end of the plotting RETURN END