C @(#)irskysub.for 17.1.1.1 (ES0-DMD) 01/25/02 17:53:06 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+++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C.COPYRIGHT (C) 1992 European Southern Observatory C.IDENT irskysub.for C.AUTHOR E. Oliva, Firenze-Arcetri C.KEYWORDS Spectroscopy, IRSPEC C C.PURPOSE Determine/apply "factor" and "shift" parameters C to optimize "obj-sky" operation. C C C.ALGORITHM C C C.INPUT/OUTPUT C C C.VERSION 1.0 Creation 02.09.1992 E. Oliva C.VERSION 1.1 25.10.1993 E. Oliva C START,STEP=1 tests eliminated (were useless and gave C problems with rounding errors). C C------------------------------------------------------- C PROGRAM SKYBUB IMPLICIT REAL(A-H,O-Z) IMPLICIT INTEGER(I-N) CHARACTER*64 CUNIT CHARACTER*72 IDENT CHARACTER*60 FOBJ,FSKY,FOUT,REFTAB C INTEGER NPIX(3) INTEGER*8 IPNTR,IAP,IBP,ICP C DOUBLE PRECISION START(3),STEP(3) C COMMON/PARAM/NX,NY,IAP,IBP,ICP,NWIND,JLIM(400) COMMON /VMR/ SS(1) COMMON/TALK/JYBAD(1000),NYBAD,STDX(1000),AVEX(1000),THRSX,IDEBUG C INCLUDE 'MID_INCLUDE:ST_DEF.INC' INCLUDE 'MID_INCLUDE:ST_DAT.INC' DATA MAXDIM/2/ IRET=1 CALL STSPRO('SKYSUB') C C GET INPUT (OBJ-DARK)/FLAT CLEANED FRAME AND MAP IT C CHECK ALSO WHETHER IT IS A 2D FRAME..... C CALL STKRDC('obj',1,1,60,IRET,FOBJ,KUNIT,KNUL,ISTAT) CALL CLNFRA(FOBJ,FOBJ,0) CALL STIGET(FOBJ,D_R4_FORMAT,F_I_MODE,F_IMA_TYPE,MAXDIM, , NAXIS,NPIX,START,STEP, , IDENT,CUNIT,IPNTR,NOBJ,ISTAT) IAP=IPNTR IF(NAXIS.NE.2) + CALL STETER(1,'Input frame must be 2-D') NX=NPIX(1) NY=NPIX(2) C C GET INPUT (SKY-DARK)/FLAT FRAME AND MAP IT C CHECK ALSO WHETHER IT HAS CORRECT NX,NY C CALL STKRDC('sky',1,1,60,IRET,FSKY,KUNIT,KNUL,ISTAT) CALL CLNFRA(FSKY,FSKY,0) CALL STIGET(FSKY,D_R4_FORMAT,F_I_MODE,F_IMA_TYPE,MAXDIM, , NAXIS,NPIX,START,STEP, , IDENT,CUNIT,IPNTR,NSKY,ISTAT) IBP=IPNTR IF(NPIX(1).NE.NX .OR. NPIX(2).NE.NY) + CALL STETER(1,'Sky frame has different size than object!') C C GET OUTPUT FRAME AND MAP IT C CALL STKRDC('out',1,1,60,IRET,FOUT,KUNIT,KNUL,ISTAT) CALL CLNFRA(FOUT,FOUT,0) CALL STIPUT(FOUT,D_R4_FORMAT,F_O_MODE,F_IMA_TYPE, , NAXIS,NPIX,START,STEP, , IDENT,CUNIT,IPNTR,NOUT,ISTAT) ICP=IPNTR C C Get "force_sky_to_zero" option C CALL STKRDI('sforce',1,1,IRET,ISFORCE,KUNIT,KNUL,ISTAT) C C Check whether a sky_table exists C CALL STKRDI('yestab',1,1,IRET,KKTAB,KUNIT,KNUL,ISTAT) IF(KKTAB.GT.0) THEN C C GET TABLE NAME (WHERE WINDOWS ARE DEFINED) C AND OPEN IT C CALL STKRDC('reft',1,1,60,IRET,REFTAB,KUNIT,KNUL,ISTAT) IOPEN=0 CALL TBTOPN(REFTAB,IOPEN,IDETAB,ISTAT) C C GET # CORRESPONDING TO XSTART,XEND,YSTART,YEND C CALL TBLSER(IDETAB,'XSTART',IX1,ISTAT) CALL TBLSER(IDETAB,'YSTART',IY1,ISTAT) CALL TBLSER(IDETAB,'XEND',IX2,ISTAT) CALL TBLSER(IDETAB,'YEND',IY2,ISTAT) C C GET INFORMATION UPON THE TABLE SIZE (N.B. NROW = # OF BAD PIXELS....) C CALL TBIGET(IDETAB,NCOL,NROW,NSORT,NCALL,NRALL,ISTAT) NWIND=NROW C C COPY TABLE ELEMENTS INTO VECTOR JLIM C (as xstart(1),xend(1),ystart(1),yend(1) etc..) C J=0 DO IROW=1,NWIND CALL TBERDR(IDETAB,IROW,IX1,VAL,NULL,ISTAT) J=J+1 JLIM(J)=VAL CALL TBERDR(IDETAB,IROW,IX2,VAL,NULL,ISTAT) J=J+1 JLIM(J)=VAL CALL TBERDR(IDETAB,IROW,IY1,VAL,NULL,ISTAT) J=J+1 JLIM(J)=VAL CALL TBERDR(IDETAB,IROW,IY2,VAL,NULL,ISTAT) J=J+1 JLIM(J)=VAL ENDDO CALL TBTCLO(IDETAB,ISTAT) ELSE C C In case no table exists C NWIND=0 ENDIF C C Get "MANUAL" values of fac,sh (if fac=0.0 then automatic) C CALL STKRDR('fac',1,1,IRET,FACT,KUNIT,KNUL,ISTAT) CALL STKRDR('sh',1,1,IRET,SHIFT,KUNIT,KNUL,ISTAT) C C Get debug flag C CALL STKRDI('debug',1,1,IRET,II,KUNIT,KNUL,ISTAT) IDEBUG=II C C Call routine which does the operations.... C CALL MIN_STD(FACT,SHIFT,ISFORCE) C C Store best values in keywords C CALL STKWRR('fac',FACT,1,1,KUNIT,ISTAT) CALL STKWRR('sh',SHIFT,1,1,KUNIT,ISTAT) CALL STDWRR(NOUT,'AB_FACTOR',FACT,1,1,KUNIT,ISTAT) CALL STDWRR(NOUT,'AB_SHIFT',SHIFT,1,1,KUNIT,ISTAT) C C RELEASE FILES, UPDATE KEYWORDS AND EXIT C CALL STSEPI END C C C SUBROUTINE MIN_STD(FAC,SH,ISFORCE) c c Used to find the values of shift (on image "A") and factor c (on image "B") which minimizes the c standard deviation of the pixels within the [x1,x2:y1,y2] c window(s). It also forces the average pixel c value in the windows to 0. c IMPLICIT REAL(A-H,O-Z) IMPLICIT INTEGER(I-N) C INTEGER*8 IAP,IBP,ICP C COMMON/PARAM/NX,NY,IAP,IBP,ICP,NWIND,JLIM(400) COMMON/TALK/JYBAD(1000),NYBAD,STDX(1000),AVEX(1000),THRSX,IDEBUG DIMENSION F(3) CHARACTER*80 OUTPUT c c c NB FACT>0 means that the user entered manual values for fac and sh c c c set first THRSX value c THRSX=2.0 CALL DIFF_AB(FAC,SH) IF(FAC.GT.0.0) THEN IF(NWIND.GT.0) CALL SKY_ZERO(ISFORCE) RETURN ENDIF CALL DEF_YBAD FAC=1.0 SH=0.0 DF=0.1 DS=0.1 FMIN=1E30 NITER=0 10 CONTINUE DO I=1,3 F(I)=FUNCT(FAC+(I-1)*DF,SH) ENDDO IF(F(1).GT.F(2) .AND. F(2).GT.F(3)) THEN FAC=FAC+2.*DF GO TO 10 ENDIF IF(F(1).LT.F(2) .AND. F(2).LT.F(3)) THEN FAC=FAC-DF GO TO 10 ENDIF IF(F(1).EQ.F(2) .AND. F(2).EQ.F(3)) THEN FAC=FAC+DF ELSE FAC=FAC+2*DF-DF*((F(3)-F(2))/(F(3)-2*F(2)+F(1))+0.50) ENDIF 11 CONTINUE DO I=1,3 F(I)=FUNCT(FAC,SH+(I-1)*DS) ENDDO IF(F(1).GT.F(2) .AND. F(2).GT.F(3)) THEN SH=SH+2.*DS GO TO 11 ENDIF IF(F(1).LT.F(2) .AND. F(2).LT.F(3)) THEN SH=SH-DS GO TO 11 ENDIF IF(F(1).EQ.F(2) .AND. F(2).EQ.F(3)) THEN SH=SH+DS ELSE SH=SH+2*DS-DS*((F(3)-F(2))/(F(3)-2*F(2)+F(1))+0.50) ENDIF FNEW=FUNCT(FAC,SH) ERR=ABS(FNEW/FMIN-1.) IF(IDEBUG.GT.0) THEN WRITE(OUTPUT,*) 'FAC,SH,FNEW',FAC,SH,FNEW CALL STTPUT(OUTPUT,ISTAT) ENDIF IF(ERR.GT.1) THEN FMIN=FNEW NITER=NITER+1 IF(NITER.GT.30) THEN CALL FAILED(1000) ELSE GO TO 10 ENDIF ENDIF IF(ERR.GT.1E-4) THEN DF=DF/2. DS=DS/2. DF=MAX(DF,1E-4) DS=MAX(DS,1E-4) FMIN=FNEW NITER=NITER+1 IF(NITER.GT.30) THEN CALL FAILED(1000) ELSE GO TO 10 ENDIF ENDIF CALL SKY_ZERO(ISFORCE) RETURN END FUNCTION FUNCT(FAC,SH) IMPLICIT REAL(A-H,O-Z) IMPLICIT INTEGER(I-N) CALL DIFF_AB(FAC,SH) CALL STD_WIND(AVEW,STDW) FUNCT=STDW RETURN END SUBROUTINE DEF_YBAD c c Subroutine which selects, whithin each window, the scan-lines c (Y=const) which are too noisy (e.g. which contain bad pixels). c The selection criterium is c sigma(Y_scanline) > THRSX*aver_sigma(all_Y_scanlines) c and go on iteratively. c c The positions of the bad scan-lines are stored in a rather weird c form: c c JYBAD = IY+(NW-1)*NY c c where c JYBAD is the stored value c IY is the Y coordinate of the bad scan line within the window NW c NW means that we are dealing with the NW-th window c NY is the number of Y pixels in the image c IMPLICIT REAL(A-H,O-Z) IMPLICIT INTEGER(I-N) C INTEGER*8 IAP,IBP,ICP C COMMON/PARAM/NX,NY,IAP,IBP,ICP,NWIND,JLIM(400) COMMON/TALK/JYBAD(1000),NYBAD,STDX(1000),AVEX(1000),THRSX,IDEBUG COMMON/VMR/SS(1) CHARACTER*80 OUTPUT CHARACTER*1 C1 IF(IDEBUG.GT.0) THEN CALL STTPUT('DEF_YBAD: defining bad-scan lines',ISTAT) ENDIF NYBAD=0 10 CONTINUE CALL STD_WIND(AVEW,STDW) IF(IDEBUG.GT.0) THEN WRITE(OUTPUT,*) + 'Entered with AVEW =',AVEW,' STDW =',STDW CALL STTPUT(OUTPUT,ISTAT) ENDIF NBAD=0 NSL=0 DO IW=1,NWIND IY1=JLIM(3+4*(IW-1)) IY2=JLIM(4+4*(IW-1)) DO IY=IY1,IY2 JYSHIFT=IY+(IW-1)*NY NSL=NSL+1 DIST=STDX(NSL)/STDW IF(DIST.GT.THRSX) THEN NBAD=NBAD+1 JYBAD(NBAD)=JYSHIFT IF(IDEBUG.GT.1) THEN WRITE(OUTPUT,*) 'Window #',IW,' Bad row at Y=',IY CALL STTPUT(OUTPUT,ISTAT) ENDIF ENDIF ENDDO ENDDO IF(IDEBUG.GT.0) THEN WRITE(OUTPUT,*) + 'Entered with NYBAD =',NYBAD,' found NBAD =',NBAD CALL STTPUT(OUTPUT,ISTAT) IF(IDEBUG.GT.1) THEN WRITE(*,*) 'Hit return to continue' READ(*,'(A1)') C1 ENDIF ENDIF IF(NBAD.GT.NYBAD) THEN NYBAD=NBAD GO TO 10 ENDIF RETURN END SUBROUTINE SKY_ZERO(ISFORCE) c c Subroutine which selects and excludes, whithin each window, the scan-lines c (Y=const) which contain signal well above (below) the typical sky noise c The selection criterium is c abs( aver(Y_scanline)-aver(all_Y_scan_lines) > c THRSY*aver_sigma(all_Y_scanlines) c and go on iteratively. c If, as sometimes happens, finds that all scan-lines are stars c retry with THRSY=THRSY*SQRT(2) and so forth. c c This routine assumes that image "C" (pointer ICP) was already created. c c This routine adds "bad scan lines" to JYBAD and must be called c only at the end of the program. c c This routine modify the image "C" and subtract the "good average" c of the sky windows from it. c c All the above operations are performed only if ISFORCE>1 c IMPLICIT REAL(A-H,O-Z) IMPLICIT INTEGER(I-N) C INTEGER*8 IAP,IBP,ICP C COMMON/PARAM/NX,NY,IAP,IBP,ICP,NWIND,JLIM(400) COMMON/TALK/JYBAD(1000),NYBAD,STDX(1000),AVEX(1000),THRSX,IDEBUG COMMON/VMR/SS(1) CHARACTER*80 OUTPUT CHARACTER*1 C1 IF(ISFORCE.LE.0) RETURN IF(IDEBUG.GT.0) THEN CALL STTPUT('SKY_ZERO: defining bad and star-Ylines',ISTAT) ENDIF THRSY=THRSX 100 CONTINUE NYBAD=0 10 CONTINUE CALL STD_WIND(AVEW,STDW) IF(IDEBUG.GT.0) THEN WRITE(OUTPUT,*) + 'Entered with AVEW =',AVEW,' STDW =',STDW CALL STTPUT(OUTPUT,ISTAT) ENDIF IF(STDW.LE.0.) THEN IF(IDEBUG.GT.0) CALL STTPUT('Expand Y-threshold',ISTAT) THRSY=THRSY*SQRT(2.) GO TO 100 ENDIF NBAD=0 NSL=0 DO IW=1,NWIND IY1=JLIM(3+4*(IW-1)) IY2=JLIM(4+4*(IW-1)) DO IY=IY1,IY2 JYSHIFT=IY+(IW-1)*NY NSL=NSL+1 DISTX=STDX(NSL)/STDW DISTY=ABS(AVEX(NSL)-AVEW)/STDW IF(DISTX.GT.THRSX .OR. DISTY.GT.THRSY) THEN NBAD=NBAD+1 JYBAD(NBAD)=JYSHIFT IF(IDEBUG.GT.1) THEN IF(DISTX.GT.THRSX) THEN WRITE(OUTPUT,*) 'Window #',IW,' Bad row at Y=',IY CALL STTPUT(OUTPUT,ISTAT) ENDIF IF(DISTY.GT.THRSY) THEN WRITE(OUTPUT,*) 'Window #',IW,' Star-row at Y=',IY CALL STTPUT(OUTPUT,ISTAT) ENDIF ENDIF ENDIF ENDDO ENDDO IF(IDEBUG.GT.0) THEN WRITE(OUTPUT,*) + 'Entered with NYBAD =',NYBAD,' found NBAD =',NBAD CALL STTPUT(OUTPUT,ISTAT) IF(IDEBUG.GT.1) THEN WRITE(*,*) 'Hit return to continue' READ(*,'(A1)') C1 ENDIF ENDIF IF(NBAD.GT.NYBAD) THEN NYBAD=NBAD GO TO 10 ENDIF c c Subtract the average value of the sky windows c DO IY=1,NY DO IX=1,NX IJC=IX+NX*(IY-1)+ICP-1 SS(IJC)=SS(IJC)-AVEW ENDDO ENDDO RETURN END SUBROUTINE STD_WIND(AVEW,STDW) c c c Define the average of the X-std (computed only along lines with Y=const) c of a a number of windows of the image "C" (pointer ICP). c Bad scan-lines (stored in JYBAD) are not considered. c The average and std of each scan-line are also stored inside STDX c and AVEX. c IMPLICIT REAL(A-H,O-Z) IMPLICIT INTEGER(I-N) C INTEGER*8 IAP,IBP,ICP C COMMON/PARAM/NX,NY,IAP,IBP,ICP,NWIND,JLIM(400) COMMON/TALK/JYBAD(1000),NYBAD,STDX(1000),AVEX(1000),THRSX,IDEBUG COMMON/VMR/SS(1) STDW=0. AVEW=0. NSL=0 DO IW=1,NWIND IX1=JLIM(1+4*(IW-1)) IX2=JLIM(2+4*(IW-1)) IY1=JLIM(3+4*(IW-1)) IY2=JLIM(4+4*(IW-1)) DO IY=IY1,IY2 c c N.B. JYSHIFT is used to codify the bad scan lines (see DEF_YBAD) c JYSHIFT=IY+(IW-1)*NY NSL=NSL+1 DO JY=1,NYBAD IF(JYBAD(JY).EQ.JYSHIFT) GO TO 11 ENDDO STDW=STDW+STDX(NSL)*STDX(NSL) AVEW=AVEW+AVEX(NSL) 11 CONTINUE ENDDO ENDDO STDW=SQRT(STDW/NSL) AVEW=AVEW/NSL RETURN END SUBROUTINE DIFF_AB(FAC,SH) c c Subroutine wich computed and stores into c the "C" image (i.e. SS(IJC)) the difference c between "A" (SS(IJA)) and "B" (SS(JJB)) with c "A" shifted by SH and "B" multiplied times FAC. c c N.B. IJA, IJB, IJB are the positions defined by the c pointers IAP,IBP and ICP c c It also computes and stores STDX and AVEX (standard deviation c and average, row by row, in the sky windows). c IMPLICIT REAL(A-H,O-Z) IMPLICIT INTEGER(I-N) C INTEGER*8 IAP,IBP,ICP INTEGER*8 IJ1,IJ2,IJ,IJB,IJC C COMMON/PARAM/NX,NY,IAP,IBP,ICP,NWIND,JLIM(400) COMMON/TALK/JYBAD(1000),NYBAD,STDX(1000),AVEX(1000),THRSX,IDEBUG COMMON/VMR/SS(1) !statt MADRID ... DO IY=1,NY DO IX=1,NX IX1=MAX(IX-1,1) IX2=MIN(IX+1,NX) IJ1=IX1+NX*(IY-1)+IAP-1 IJ2=IX2+NX*(IY-1)+IAP-1 IJ=IX+NX*(IY-1)+IAP-1 IF(SH.GT.0) ASH=(1-SH)*SS(IJ)+SH*SS(IJ2) IF(SH.LE.0) ASH=(1+SH)*SS(IJ)-SH*SS(IJ1) IJB=IX+NX*(IY-1)+IBP-1 IJC=IX+NX*(IY-1)+ICP-1 SS(IJC)=ASH-FAC*SS(IJB) ENDDO ENDDO NSL=0 DO IW=1,NWIND IX1=JLIM(1+4*(IW-1)) IX2=JLIM(2+4*(IW-1)) IY1=JLIM(3+4*(IW-1)) IY2=JLIM(4+4*(IW-1)) DO IY=IY1,IY2 c c N.B. JYSHIFT is used to codify the bad scan lines (see DEF_YBAD) c JYSHIFT=IY+(IW-1)*NY SUM=0.0 SUM2=0.0 NAVE=0 DO IX=IX1,IX2 IJC=IX+NX*(IY-1)+ICP-1 VAL=SS(IJC) SUM=SUM+VAL SUM2=SUM2+VAL*VAL NAVE=NAVE+1 ENDDO IF(NAVE.LE.1) CALL FAILED(-1) NSL=NSL+1 IF(NSL.GT.1000) CALL FAILED(7777) STDX(NSL)=SQRT((SUM2-SUM*SUM/NAVE)/(NAVE-1)) AVEX(NSL)=SUM/NAVE ENDDO ENDDO RETURN END SUBROUTINE FAILED(II) IMPLICIT REAL(A-H,O-Z) IMPLICIT INTEGER(I-N) c c Give error message depending on II c IF(II.EQ.7777) THEN CALL STTPUT + ('More than 1000 scan-lines in window(s)',ISTAT) CALL STTPUT + ('Need to modify skysub.for',ISTAT) GO TO 99 ENDIF IF(II.EQ.-1) THEN CALL STTPUT + ('Sky-windows must be at least 3 pixels wide in X',ISTAT) GO TO 99 ENDIF IF(II.EQ.-2) THEN CALL STTPUT + ('No good scan line in sky window(s)',ISTAT) GO TO 99 ENDIF IF(II.EQ.1000) THEN CALL STTPUT + ('WARNING! No convergency after 30 iterations',ISTAT) CALL STTPUT + ('May be better to try manually ',ISTAT) RETURN ENDIF 99 CONTINUE CALL STETER(1,' ') RETURN END