C @(#)spoext.for	17.1.1.1 (ES0-DMD) 01/25/02 17:54:01
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 @(#)spoext.for	17.1.1.1  (ESO)  01/25/02  17:54:01 
C +++++++++++++++++++++++++++++++++++++++++++++++++++++++
C .COPYRIGHT   (C) 1993 European Southern Observatory    
C .IDENT       .prg                                      
C .AUTHORS     Pascal Ballester (ESO/Garching)           
C              Cristian Levin   (ESO/La Silla)           
C .KEYWORDS    Spectroscopy, Long-Slit                   
C .PURPOSE                                               
C .VERSION     1.0  Package Creation  17-MAR-1993        
C -------------------------------------------------------

C @(#)spoext.for	4.2 (ESO-IPG) 11/20/92 14:47:36
          PROGRAM SPOEXT
C+
C
C.IDENTIFICATION
C
C  PROGRAM  SPOEXT                VERSION 1.0        880420
C                                        
C  Otmar Stahl                        Landessternwarte Heidelberg
C     M. Peron  , modif suggested by T. Oosterloo                                      900830
C.KEYWORDS
C
C  SPECTRA, CCD, DATA EXTRACTION
C
C.PURPOSE
C
C  EXTRACT STELLAR SPECTRA FROM A CCD IMAGE. THE INDIVIDUAL ROWS OF 
C  THE CCD IMAGE WHICH CONTAIN THE FLUX ARE ADDED. THE WEIGHTS ARE
C  CHOSEN SUCH THAT THE SUMMED SPECTRUM HAS THE OPTIMAL S/N-RATIO.
C  COSMIC RAY HITS ARE LOCATED BY ANALYSING THE PROFILE PERPENDICULAR 
C  TO THE DISPERSION AND MASKED OUT. 
C  
C.ALGORITHM
C
C  THE METHOD IS DESCRIBED IN DETAIL IN THE PAPER BY 
C  Keith Horne in PASP 98, pg. 609. 
C  THE MASK USED TO REMOVE THE COSMIC RAY HITS IS KEPT IN THE
C  MASK FILE cosmics.bdf.
C
C  EXTRACT/SPECTRUM  output input sky order,iter read-out-noise,g,thresh
C
C  OUT_A                output spectrum
C  IN_A                 input image (must be bias and dark subtracted 
C                       (if relevant) and divided by a normalized flat 
C                       field, not sky subtracted)
C  IN_B                 sky image (normally derived from IN_A in advance)
C  INPUTI               fit order for the polynomials, number of 
C                       iterations for cosmic ray removal
C  INPUTR                read-out-noise, conversion factor (e/ADU),
C                       threshold for cosmic ray removal
C
      IMPLICIT NONE
      CHARACTER*60    INFR,SKYFR,OUTFR
      CHARACTER       CUNIT*64,IDENT*72
      CHARACTER       CUNITS*64,IDENTS*72
      CHARACTER*80    OUTPUT
      INTEGER       NPIX(2),MPIX(2),IFIT(2),NUM
      INTEGER       IAV,ISTAT,NAXIS
      INTEGER*8     IPNTRA,IPNTRB,IPNTRC
      INTEGER       NBYTES
      INTEGER*8     IPNTR1,IPNTR2,IPNTR3,IPNTR4,IPNTR5,IPNTR6
      INTEGER       KUN,KNUL,IDI,IDS,IDO,MADRID,IDM,ISS
      REAL          CCD(3)
      DOUBLE PRECISION  START(2),STEP(2),STARTS(2),STEPS(2)
      INCLUDE       'MID_INCLUDE:ST_DEF.INC'
      COMMON/VMR/MADRID(1)
      INCLUDE       'MID_INCLUDE:ST_DAT.INC'
C
      CALL STSPRO('SPOEXT')
C
C get input and map images
C
      CALL STKRDC('IN_A',  1,1,60,IAV,INFR,KUN,KNUL,ISTAT)
      CALL STKRDC('IN_B',  1,1,60,IAV,SKYFR,KUN,KNUL,ISTAT)
      CALL STKRDC('OUT_A', 1,1,60,IAV,OUTFR,KUN,KNUL,ISTAT)
      CALL STKRDI('PUTI',1,2,IAV,IFIT,KUN,KNUL,ISTAT)
      CALL STKRDR('PUTR',1,3,IAV,CCD,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 STIGET(SKYFR,D_R4_FORMAT,F_I_MODE,F_IMA_TYPE,
     .2,NAXIS,MPIX,STARTS,STEPS,IDENTS,CUNITS,IPNTR6,IDS,ISTAT)
C
      IF ((MPIX(1).NE.NPIX(1)).OR.(MPIX(2).NE.NPIX(2))) THEN
          WRITE(OUTPUT,100) 
100       FORMAT('OBJECT AND SKY FRS TO NOT MATCH')
          CALL STTPUT(OUTPUT,ISTAT)
          GO TO 999
      ENDIF
C
      CALL STIPUT(OUTFR,D_R4_FORMAT,F_O_MODE,F_IMA_TYPE,
     .1,NPIX(1),START(1),STEP(1),IDENT,CUNIT,IPNTR2,IDO,ISTAT)
C
      CALL STIPUT('cosmics',D_R4_FORMAT,F_O_MODE,F_IMA_TYPE,
     .NAXIS,NPIX,START,STEP,IDENT,CUNIT,IPNTR3,IDM,ISTAT)
C
C reserve  space for arrays
C
      NUM = NPIX(1)*NPIX(2)*4
      NBYTES = NPIX(1)*4
      CALL TDMGET(NUM,IPNTR4,ISS)     
      CALL TDMGET(NUM,IPNTR5,ISS)     
      CALL TDMGET(NBYTES,IPNTRA,ISS)
      CALL TDMGET(NBYTES,IPNTRB,ISS)
      CALL TDMGET(NBYTES,IPNTRC,ISS)
C
C DO THE REAL WORK
C
      CALL XTRACT (MADRID(IPNTR1),MADRID(IPNTR2),MADRID(IPNTR3), 
     &MADRID(IPNTR4),MADRID(IPNTR5),MADRID(IPNTR6),MADRID(IPNTRA),
     &MADRID(IPNTRB),MADRID(IPNTRC),NPIX,IFIT,CCD)
C
C COPY HISTORY
C
C      CALL COPYDSC_ST(INFR,OUTFR,4,'HISTORY',ISTAT)
C
      CALL TDMFRE(NUM,IPNTR4,ISS)     
      CALL TDMFRE(NUM,IPNTR5,ISS)     
      CALL TDMFRE(NBYTES,IPNTRA,ISS)
      CALL TDMFRE(NBYTES,IPNTRB,ISS)
      CALL TDMFRE(NBYTES,IPNTRC,ISS)
999   CALL STSEPI
      END
C
C**************************************************************************
      SUBROUTINE XTRACT(IN,OUT,MASK,PROF,VARI,SKY,X,Y,Z,NPIX,IFIT,CCD)
C 
C THIS ROUTINE DOES THE EXTRACTION OF THE SPECTRUM
C
      IMPLICIT NONE
C      INTEGER        MASK(1) ! DECLARED AS INTEGER*2 !!!
      REAL            MASK(1) ! DECLARED AS INTEGER*2 !!!
      INTEGER        NPIX(2),IFIT(2)
      INTEGER        J,K,KK,L,II,M,N,IORD,ITER,IWIND
      REAL           IN(1),OUT(1),PROF(1),VARI(1),SKY(1)
      REAL           A(10),B(10),S(10),W(10)     
      REAL           RON,G,SIGMA,V0,SUM,CHISQ,YY,DIFF,DIFF2,QUOT
      REAL           SUM1,SUM2,RMAX
      REAL           X(1),Y(1),Z(1),CCD(3)
C                                                                       
      M = NPIX(1)
      N = NPIX(2)
C
      IORD = IFIT(1)
      ITER = IFIT(2)
      IWIND = 3
C
      RON   = CCD(1)
      G     = CCD(2)
      SIGMA = CCD(3)
C
      SIGMA = SIGMA*SIGMA
      V0 = RON**2/G**2
C
C  EXTRACT INITIAL SPECTRUM (JUST THE SUM) AND SET mask TO 1
C
      DO 20 J = 1,M
        SUM = 0.
        DO 10 K = 1,N
          II = (K-1)*M+J
          MASK(II) = 1
          SUM = SUM+IN(II)-SKY(II)
10      CONTINUE
        X(J) = J
        OUT(J) = SUM
20    CONTINUE
C
C COMPUTE SPATIAL PROFILE, I.E. THE 2-DIM SPECTRUM 
C IS DIVIDED (ROW BY ROW) BY THE SUM
C                                                 
      DO 40 J = 1,M
        DO 30 K = 1,N
          II = (K-1)*M+J
          IF (OUT(J).EQ.0) OUT(J)=1
          PROF(II) = (IN(II)-SKY(II))/OUT(J)
30      CONTINUE
40    CONTINUE
C
C  FILTER SPATIAL PROFILE AND FIT WITH POLYNOMIAL
C
      DO 70 K = 1,N
        DO 50 J = 1,M
          II = (K-1)*M+J
          Y(J) = PROF(II)
50      CONTINUE
C
C  MEDIAN FILTER OF ROWS
C
        CALL FILTER(Y,Z,M,IWIND)
C 
C  FIT WITH POLYNOMIAL 
C
        CALL LSORTH(X,Y,A,B,S,W,M,CHISQ,IORD)
        DO 60 J = 1,M
          II = (K-1)*M+J
          CALL POLY(X(J),YY,A,B,S,W,IORD)
          PROF(II) = YY
60      CONTINUE
70    CONTINUE
C
C NORMALIZE SPATIAL PROFILE AND ENFORCE POSITIVITY
C
      DO 100 J = 1,M
        SUM=0.
        DO 80 K = 1,N
          II = (K-1)*M+J
          IF(PROF(II).LT.0.) PROF(II)=0.
          SUM = SUM+PROF(II)
80      CONTINUE
        DO 90 K = 1,N
          II = (K-1)*M+J
          PROF(II) = PROF(II)/SUM
90      CONTINUE
100   CONTINUE
C
C ITERATE (in the present version, the spatial profile is not 
C iterated, but only the mask and the variance)
C                                                                       
      DO 160 L = 1,ITER                                                  
C                                                                       
C  compute variance estimate                 
C
        DO 120 J = 1,M                                               
          DO 110 K = 1,N                  
           II = (K-1)*M+J
           VARI(II) = V0 + ABS(OUT(J)*PROF(II)+SKY(II))/G
110       CONTINUE
120     CONTINUE
C                                                                       
C MASK COSMIC RAY HITS (ONLY ONE PER COLUMN AND ITERATION)
C
        DO 150 J = 1,M
          RMAX = 1.
          KK = 0
          DO 130 K = 1,N
            II = (K-1)*M+J                                                     
            DIFF = IN(II)-SKY(II)-OUT(J)*PROF(II)
            DIFF2 = DIFF*DIFF
            QUOT = (DIFF2/(SIGMA*VARI(II)))*MASK(II)
            IF(QUOT.GT.RMAX) THEN
              RMAX = QUOT
              KK = K
            ENDIF
130       CONTINUE
          IF(KK.NE.0) THEN
            II = (KK-1)*M+J
            MASK(II) = 0
          ENDIF
C
C COMPUTE WEIGHTED MEAN WITH COSMICS MASKED
C
          SUM1 = 0.
          SUM2 = 0.            
          DO 140 K = 1,N
            II = (K-1)*M+J
            SUM1 = SUM1 + (MASK(II)*PROF(II)*(IN(II)-SKY(II)))/VARI(II)
            SUM2 = SUM2 + (MASK(II)*PROF(II)*PROF(II))        /VARI(II)
140       CONTINUE
          OUT(J) = SUM1/SUM2
150     CONTINUE
160   CONTINUE
C                                                                       
      RETURN
      END                                                           
C
C****************************************************************************
      SUBROUTINE FILTER(ARR,Z,M,IWIND)
C
      IMPLICIT NONE
      INTEGER     II,IWIND,J,M,K
      REAL        ARR(1),Z(1),Y(50),XMED
C
      II = 2*IWIND+1
C
C  SORT DATA IN ARRAY X AND CALL MEDIAN FILTER
C
      DO 20 J = IWIND+1,M-IWIND
        DO 10 K = 1,II
          Y(K) = ARR(J+K-IWIND-1)
10      CONTINUE
C
        CALL MDIAN1(Y,II,XMED)
C
        Z(J) = XMED
20    CONTINUE
C
C  COPY DATA
C
      DO 30 J = IWIND+1,M-IWIND
         ARR(J) = Z(J) 
30    CONTINUE
C
      RETURN
      END
C
C**************************************************************************
      SUBROUTINE MDIAN1(ARR,N,XMED)
C
C COMPUTE MEDIAN VALUE OF AN ARRAY
C
      IMPLICIT NONE
      INTEGER        N,N2
      REAL           ARR(N),XMED
C
C  SORT DATA
C
      CALL SORT(N,ARR)
C
C  COMPUTE MEDIAN
C
      N2 = N/2
      IF(2*N2.EQ.N) THEN
        XMED = 0.5*(ARR(N2)+ARR(N2+1))
      ELSE
        XMED = ARR(N2+1)
      ENDIF
C
      RETURN
      END
C
C************************************************************************
      SUBROUTINE SORT(N,ARR)
C
C SORTS AN ARRAY
C 
      IMPLICIT NONE
      INTEGER     N,I,J
      REAL        A,ARR(N)
C
      DO 20 J = 2,N
        A = ARR(J)      
        DO 5 I = J-1,1,-1
          IF(ARR(I).LE.A) GOTO 10
          ARR(I+1) = ARR(I)
5       CONTINUE
      I = 0
10    ARR(I+1) = A
20    CONTINUE
C
      RETURN
      END
C
C*****************************************************************************
C        POLY : INPUT: XFIT AND POLYNOMIAL COEFFICIENTS COMPUTED
C              WITH LSORTH
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                 THE ARRAY Y IS RETURNED AS A VECTOR
C                 CONTAINING THE RESIDUALS.
C                 THE ROUTINE IS DESCRIBED IN PASP 91 (1979), PG. 546.    
C 
C                 THE POLYNOMIAL COEFFICIENTS 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. 
10      B(I)=0.
        P(1)=1.
        DO 20 N=1,M
        W(1)=W(1)+1.   
        S(1)=S(1)+Y(N)  
20      A(1)=A(1)+X(N)
        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