C @(#)ffnorm.for	19.1 (ES0-DMD) 02/25/03 14:25:21
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	19.1  (ESO)  02/25/03  14:25:21 
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
