C @(#)ffnorm.for 17.1.1.1 (ES0-DMD) 01/25/02 17:53:59 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 @(#)ffnorm.for 17.1.1.1 (ESO) 01/25/02 17:53:59 C +++++++++++++++++++++++++++++++++++++++++++++++++++++++ C .COPYRIGHT (C) 1993 European Southern Observatory C .IDENT ffnorm.for C .AUTHORS Otmar STAHL (Landessternwarte Heidelberg) C .KEYWORDS Spectroscopy, Long-Slit, Flat-Field, C Normalization C .PURPOSE C C Normalize a flat field image which is extended along columns. C The image is divided by a polynomial along the rows, in order C to achieve flux values close to unity. The bias value is subtracted C before the division. The normalized flat field and the polynomial C (in the form of a spectrum) are returned. C C .ALGORITHM C C THE PROGRAM FIRST LOCATES THE EDGES OF THE SLIT IN THE DIRECTION C PERPENDICULAR TO THE DISPERSION BY LOOKING FOR A FLUX DECREASE. C THEN THE MEAN OF THE SPECTRUM IS TAKEN AND FITTED WITH A POLYNOMIAL. C C .INPUT/OUTPUT C C NORMALIZE/FLAT norm_image fit_func input bias fit_order C C OUT_A normalized output frame C OUT_B function used for normalization C IN_A input frame C INPUTR bias value to subtract C INPUTI order of the polynomial fit C C .VERSION 1.0 From package Long 22-APR-1988 C 2.0 Package Creation 17-MAR-1993 C ------------------------------------------------------- IMPLICIT NONE INTEGER NPIX(2),IAV,ISTAT,NAXIS,IORD,NUM INTEGER*8 IPNTR1,IPNTR2,IPNTR3,IPNTR4,IPNTR5 INTEGER KUN, KNUL, MADRID, IDI, IDO, IDF, ISS CHARACTER*60 INFR,OUTFR,FITTED CHARACTER CUNIT*64,IDENT*72 REAL BIAS DOUBLE PRECISION START(2),STEP(2) INCLUDE 'MID_INCLUDE:ST_DEF.INC' COMMON/VMR/MADRID(1) INCLUDE 'MID_INCLUDE:ST_DAT.INC' C CALL STSPRO('FFNORM') C C GET INPUT AND MAP IMAGES C CALL STKRDC('IN_A', 1,1,60,IAV,INFR,KUN,KNUL,ISTAT) CALL STKRDC('OUT_A',1,1,60,IAV,OUTFR,KUN,KNUL,ISTAT) CALL STKRDC('OUT_B',1,1,60,IAV,FITTED,KUN,KNUL,ISTAT) CALL STKRDR('INPUTR',1,1,IAV,BIAS,KUN,KNUL,ISTAT) CALL STKRDI('INPUTI',1,1,IAV,IORD,KUN,KNUL,ISTAT) C CALL STIGET(INFR, D_R4_FORMAT, F_I_MODE, F_IMA_TYPE, .2,NAXIS,NPIX,START,STEP,IDENT,CUNIT,IPNTR1,IDI,ISTAT) C CALL STIPUT(OUTFR,D_R4_FORMAT, F_O_MODE, F_IMA_TYPE, .NAXIS,NPIX,START,STEP,IDENT,CUNIT,IPNTR2,IDO,ISTAT) C CALL STIPUT(FITTED,D_R4_FORMAT, F_O_MODE, F_IMA_TYPE, .1,NPIX(1),START(1),STEP(1),IDENT,CUNIT,IPNTR3,IDF,ISTAT) C C GET VIRTUAL MEMORY FOR ARRAYS C NUM = NPIX(1) IF(NUM.LT.NPIX(2)) NUM = NPIX(2) NUM = NUM*4 CALL TDMGET(NUM,IPNTR4,ISS) CALL TDMGET(NUM,IPNTR5,ISS) C C NORMALIZE FRAME C CALL NORM(MADRID(IPNTR1),MADRID(IPNTR2),MADRID(IPNTR3), . MADRID(IPNTR4),MADRID(IPNTR5),NPIX,BIAS,IORD) C C COPY HISTOR C C CALL COPYDSC_ST(INFR,OUTFR,4,'HISTORY',ISTAT) C CALL COPYDSC_ST(INFR,FITTED,4,'HISTORY',ISTAT) C CALL TDMFRE(NUM,IPNTR4,ISS) CALL TDMFRE(NUM,IPNTR5,ISS) CALL STSEPI END C C************************************************************************* SUBROUTINE NORM(IN,OUT,BLAZE,X,Y,NPIX,BIAS,IORD) C C normalize the image in dispersion correction C IMPLICIT NONE INTEGER NPIX(2),M,N,J,II,K,IMAX,J1,J2,IORD REAL IN(1),OUT(1),BLAZE(1),FRAC,SUM,BIAS,RMAX,RMIN REAL X(1),Y(1),CHISQ,YY REAL A(10),B(10),S(10),W(10) C C fixed parameter: C frac: used to locate the edges of the useful flat field (frac*max flux) C FRAC = 0.5 M = NPIX(1) N = NPIX(2) C C DETERMINE PROFILE PERPENDICULAR TO DISPERSION C DO 20 J = 1,N SUM = 0. II = (J-1)*M DO 10 K = 1,M SUM = SUM+IN(II+K) 10 CONTINUE Y(J) = SUM-N*BIAS 20 CONTINUE C C LOCATE BEST EXPOSED PART OF SPECTRUM C RMAX = 0. DO 30 J = 1,N IF(Y(J).GT.RMAX) THEN IMAX = J RMAX = Y(J) ENDIF 30 CONTINUE C RMIN = RMAX*FRAC J1 = 1 J2 = N DO 40 J = IMAX,N J2 = J IF(Y(J).LE.RMIN) GOTO 100 40 CONTINUE 100 J2 = J2-1 DO 50 J = IMAX,1,-1 J1 = J IF(Y(J).LE.RMIN) GOTO 200 50 CONTINUE 200 J1 = J1+1 II = J2-J1+1 C C find mean spectrum (average perpendicular to dispersion) C DO 70 J = 1,M SUM = 0. DO 60 K = J1,J2 SUM = SUM+IN((K-1)*M+J) 60 CONTINUE X(J) = J Y(J) = SUM/II-BIAS 70 CONTINUE C C FIT POLYNOMIAL C CALL LSORTH(X,Y,A,B,S,W,M,CHISQ,IORD) C DO 80 J = 1,M CALL POLY(X(J),YY,A,B,S,W,IORD) BLAZE(J) = YY 80 CONTINUE C C NORMALIZE C DO 90 J = 1,M DO 95 K = 1,N II = (K-1)*M+J OUT(II) = (IN(II)-BIAS)/BLAZE(J) 95 CONTINUE 90 CONTINUE C C RETURN END C C***************************************************************************** C POLY : EVALUATES A POLYNOMIAL AT A POSITION XFIT C INPUT: XFIT AND POLYNOMIAL COEFFICIENTS COMPUTED WITH LSORTH C C OUTPUT: YFIT C SUBROUTINE POLY(XFIT,YFIT,A,B,S,W,IORD) C IMPLICIT NONE INTEGER K,IORD REAL A(10),B(10),S(10),W(10),P(10) REAL XFIT,YFIT C P(1)=1. YFIT=S(1) C P(2)=XFIT-A(1) YFIT=YFIT+S(2)*P(2) C DO 10 K=3,IORD+1 P(K)=(XFIT-A(K-1))*(P(K-1))-B(K-1)*P(K-2) YFIT=YFIT+S(K)*P(K) 10 CONTINUE RETURN END C C*************************************************************************** C LSORTH : FITS ORTHOGONAL POLYNOMIAL OF MAXIMUM DEGREE IORD C THROUGH X-Y TABLE. AN F-TEST DECIDES WHETHER THE FIT C STOPS ALREADY AT A LOWER ORDER THAN IORD. C INPUT: THE ARRAYS X AND Y OF SIZE M CONTAIN THE VARIABLES. C THE ARRAY Y IS RETURNED AS A VECTOR WITH THE RESIDUALS. C THE ROUTINE IS DESCRIBED IN PASP 91 (1979), PG. 546. C C OUTPUT: THE POLYNOMIAL COEFFICIENTS AND CHISQ ARE RETURNED C SUBROUTINE LSORTH (X,Y,A,B,S,W,M,CHISQ,IORD) C IMPLICIT NONE INTEGER IMAX,IMAXST,IFAUTO,IORD,I,N,M,I1,IFF,K,K1 INTEGER IFTEST REAL X(1),Y(1),A(10),B(10),S(10),W(10),P(10) REAL CHISQ,DEL,XNU,F,F95,DCHI,P2 C DATA IMAXST/10/,IFTEST/2/ IMAX=IMAXST IFAUTO=0 IF(IORD .LT. 10) GOTO 110 IORD=9 IFAUTO=1 110 CONTINUE IF(IORD .NE. 0) IMAX=MAX0(IABS(IORD)+1,2) DO 10 I=1,IMAXST W(I)=0. S(I)=0. A(I)=0. B(I)=0. 10 CONTINUE P(1)=1. DO 20 N=1,M W(1)=W(1)+1. S(1)=S(1)+Y(N) A(1)=A(1)+X(N) 20 CONTINUE S(1)=S(1)/W(1) A(1)=A(1)/W(1) I=1 XNU=M-1 30 IFF=1 40 I1=I IF(IMAX .GT. I) I = I+1 CHISQ=0. DO 70 N=1,M P(2)=X(N)-A(1) IF(I .EQ. 2) GOTO 60 DO 50 K=3,I K1=K-1 50 P(K)=(X(N)-A(K1))*P(K1)-B(K1)*P(K-2) 60 DEL=Y(N)-S(I1)*P(I1) Y(N)=DEL CHISQ=CHISQ+DEL*DEL IF(IMAX .LE. I1) GOTO 70 P2=P(I)**2 S(I)=S(I)+DEL*P(I) A(I)=A(I)+X(N)*P2 W(I)=W(I)+P2 70 CONTINUE IF(IMAX .LE. I1) GOTO 80 A(I)=A(I)/W(I) B(I)=W(I)/W(I1) S(I)=S(I)/W(I) XNU=XNU-1. DCHI=S(I)**2*W(I) IF(DCHI .GE. CHISQ) GOTO 30 F=XNU*DCHI/(CHISQ-DCHI) F95=3.84+(10.+(12.+(30.+105./XNU/XNU)/XNU)/XNU)/XNU IF(F .GT. F95) GOTO 30 IF(IFAUTO.EQ.0) GOTO 30 ! F-TEST OFF XNU=XNU+1. IFF=IFF+1 S(I)=0. IF(IFF .LE. IFTEST) GOTO 40 80 IORD=MIN0(IMAX-1,I1)-IFF+1 RETURN END