C @(#)spskyfit.for	19.1 (ES0-DMD) 02/25/03 14:25:24
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 @(#)spskyfit.for	19.1  (ESO)  02/25/03  14:25:24 
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 @(#)skyfit.for	2.1.1.1 (ESO-IPG) 9/27/91 19:06:46
         PROGRAM SKYFIT
C+
C
C.IDENTIFICATION
C
C  PROGRAM  SKYFIT                            VERSION 1.0        880420 OS
C
C  Otmar Stahl                          Landessternwarte Heidelberg
C
C.KEYWORDS
C
C  SPECTRA, SKY DETERMINATION
C
C.PURPOSE
C
C  DETERMINE A SKY FRAME  BY FITTING A POLYNOMIAL TO EVERY COLUMN
C  OF A FRAME. THE SKY IS MEASURED IN TWO SELECTED SKY WINDOWS. 
C
C.ALGORITHM
C
C  THE SKY IS TAKEN IN TWO WINDOWS (NORMALLY ONE BELOW AND ONE ABOVE 
C  THE OBJECT. THESE DATA ARE FITTED WITH A POLYNOMIAL. THE COEFFI-
C  CIENTS OF THE POLYNOMIAL ARE USED TO CREATE A SKY FRAME OF THE SIZE
C  OF THE INPUT FRAME. TWO MODES ARE POSSIBLE: 
C  1) FIT THE AVERAGE SPATIAL PROFILE AND SCALE TO THE COLUMNS AND
C  2) FIT EVERY COLUMN INDEPENDENTLY
C
C.INPUT/OUTPUT
C
C  FIT/COULMN OUTPUT INPUT A1,A2,B1,B2 ORDER MODE GAIN,RON,THRESHOLD RADIUS
C
C  OUT_A            output frame
C  IN_A             input frame
C  INPUTI/1/7       sky windows, two window can be selected, typically
C                   one below and one above the object spectrum, 
C                   the order of the fit and the mode.
C
      IMPLICIT NONE
      INTEGER          STAT,NPIX(2),IPAR(10),I,IAV,NAXIS
      INTEGER*8        IPN,IPNR,IPNTRA,IPNTRB,IPNTRC
      INTEGER          ICHECK
      INTEGER          MADRID, KUN, KNUL, IDI, IDO, ISS
      INTEGER          NBYT1, NBYT2
      REAL             RPAR(3)
      CHARACTER*80        INFILE,OUTFIL
      CHARACTER*48       CUNIT
      CHARACTER*72       IDENT
      CHARACTER*80       OUTPUT
      DOUBLE PRECISION START(2),STEP(2)
      INCLUDE          'MID_INCLUDE:ST_DEF.INC'
      COMMON/VMR/MADRID(1)
      INCLUDE          'MID_INCLUDE:ST_DAT.INC'
C
C  INITIALIZE MIDAS ENVIRONMENT
C
      CALL STSPRO('SKYFIT')
C
C  GET PARAMETERS 
C
      CALL STKRDC('IN_A',1,1,80,IAV,INFILE,KUN,KNUL,STAT)
      CALL STKRDC('OUT_A',1,1,80,IAV,OUTFIL,KUN,KNUL,STAT)
      CALL STKRDI('INPUTI',1,7,IAV,IPAR,KUN,KNUL,STAT)
      CALL STKRDR('INPUTR',1,3,IAV,RPAR,KUN,KNUL,STAT)

C
C  MAP INPUT FILE
C
      CALL STIGET(INFILE,D_R4_FORMAT,F_I_MODE,F_IMA_TYPE,
     .2,NAXIS,NPIX,START,STEP,IDENT,CUNIT,IPN,IDI,STAT)
C
C CHECK LIMITS
C
      ICHECK = 0
      IF(IPAR(1) .GT. IPAR(2)) ICHECK = 1
      IF(IPAR(3) .GT. IPAR(4)) ICHECK = 1
      IF(IPAR(1) .LT. 1)       ICHECK = 1
      IF(IPAR(4) .GT. NPIX(2)) ICHECK = 1 
      IF (ICHECK.EQ.1) THEN
          WRITE(OUTPUT,100)
          CALL STTPUT(OUTPUT,STAT)
          WRITE(OUTPUT,200) (IPAR(I),I=1,4)
          CALL STTPUT(OUTPUT,STAT)
          WRITE(OUTPUT,300) NPIX(2)
          CALL STTPUT(OUTPUT,STAT)
          GOTO 999
      ENDIF
C
C  MAP OUTPUT FILE
C
      CALL STIPUT(OUTFIL,D_R4_FORMAT,F_O_MODE,F_IMA_TYPE,
     .2,NPIX,START,STEP,IDENT,CUNIT,IPNR,IDO,STAT)
C
C  get virtual memory for arrays
C
      NBYT1 = NPIX(2)*4
      NBYT2 = NPIX(1)*4
      CALL TDMGET(NBYT1,IPNTRA,ISS)
      CALL TDMGET(NBYT1,IPNTRB,ISS)
      CALL TDMGET(NBYT2,IPNTRC,ISS)
C
C      IPAR(7)=NPIX(1)
C      IPAR(8)=NPIX(2)
C
C DO THE FITTING IN THE SELECTED MODE
C
      IF (IPAR(6).EQ.0) THEN      
           CALL FITPL1(NPIX(1),NPIX(2),MADRID(IPN),MADRID(IPNR),
     .     MADRID(IPNTRA),MADRID(IPNTRB),MADRID(IPNTRC),IPAR)
      ELSE
           CALL FITPL2(NPIX(1),NPIX(2),MADRID(IPN),MADRID(IPNR),
     .     MADRID(IPNTRA),MADRID(IPNTRB),MADRID(IPNTRC),IPAR,RPAR)
      ENDIF
C
      CALL TDMFRE(NBYT1,IPNTRA,ISS)
      CALL TDMFRE(NBYT1,IPNTRB,ISS)
      CALL TDMFRE(NBYT2,IPNTRC,ISS)
999   CALL STSEPI
100     FORMAT('ERROR IN LIMITS')
200     FORMAT('USED WINDOWS',4(1X,I4))
300     FORMAT('ALLOWED LIMITS ARE 1',1X,I4)
      END

C******************************************************************************
C FIRST VERSION: FIT ONE POLYNOMIAL FOR ALL COLUMNS. ASSUMES THAT 
C THE NORMALIZED STATIAL PROFILE IS CONSTANT WITH X-POSITION. 
C 
      SUBROUTINE FITPL1(IX,IY,IN,SKY,Y,FLUX,X,IPAR)
C
      IMPLICIT NONE
      INTEGER         IPAR(10),I,N,J,NUM,ORDER,IX,IY
      REAL            IN(IX,IY),SKY(IX,IY)
      REAL            A(10),B(10),S(10),W(10)
      REAL            X(1),Y(1),FLUX(1),YY,CHISQ,SUM
C
      ORDER = IPAR(5)
      NUM = IPAR(2)-IPAR(1)+IPAR(4)-IPAR(3)+2
C
C------mean profile perpendicular to dispersion
C
      N = 0
      DO 20 J=IPAR(1),IPAR(2)
        N = N+1
        Y(N) = J
        FLUX(N) = 0.
        DO 10 I = 1,IX
          FLUX(N) = FLUX(N)+IN(I,J)
10      CONTINUE 
20    CONTINUE
      DO 40 J=IPAR(3),IPAR(4)
        N = N+1
        Y(N) = J
        FLUX(N) = 0.
        DO 30 I = 1,IX
          FLUX(N) = FLUX(N)+IN(I,J)
30      CONTINUE
40    CONTINUE
C
C normalize and enforce positivity
C
      SUM = 0.
      DO 50 I = 1,NUM
        IF(FLUX(I).LT.0.) FLUX(I) = 0.
        SUM = SUM+FLUX(I)
50    CONTINUE
      DO 60 I = 1,NUM
        FLUX(I) = FLUX(I)/SUM
60    CONTINUE
C
C------mean sky spectrum as scaling factor
C
      DO 90 I = 1,IX
        X(I) = 0.
        DO  70 J = IPAR(1),IPAR(2)
          X(I) = X(I)+IN(I,J)
70      CONTINUE
        DO 80 J = IPAR(3),IPAR(4)
          X(I) = X(I)+IN(I,J)
80      CONTINUE
90    CONTINUE
C
C the case of ORDER = 0 is treated seperately
C
      IF (ORDER.LE.0) THEN
        SUM = 0
        DO 100 J = 1,NUM
          SUM = SUM+FLUX(J)
100     CONTINUE
        SUM = SUM/NUM
        DO 120 J =  1,IY
          DO 110 I = 1,IX
            SKY(I,J) = SUM*X(I)
110       CONTINUE
120     CONTINUE
      ELSE
C
C fit polynomial and scale it for all rows in the column
C
        CALL LSORTH(Y,FLUX,A,B,S,W,NUM,CHISQ,ORDER)
        DO 140 J =  1,IY
          YY = J
          CALL POLY(YY,FLUX(J),A,B,S,W,ORDER)
          DO 130 I = 1,IX
            SKY(I,J) = FLUX(J)*X(I)
130       CONTINUE
140     CONTINUE
      ENDIF
C
      RETURN
      END
C
C******************************************************************************
C SECOND VERSION: FIT AN INDEPENDENT POLYNOMIAL FOR EVERY COLUMN
C
      SUBROUTINE FITPL2(IX,IY,IN,SKY,Y,FLUX,X,IPAR,RPAR)
C
      IMPLICIT NONE
      INTEGER         IPAR(10),I,N,J,NUM,ORDER,IX,IY
      INTEGER         NS,COSMIC(4000),RADIUS
      REAL            IN(IX,IY),SKY(IX,IY),RPAR(3)
      REAL            A(10),B(10),S(10),W(10)
      REAL            X(1),Y(1),FLUX(1),YY,CHISQ,SUM
      REAL            WDT(4000),WDTY(4000)
      REAL            GAIN,RON,THRESHOLD
C

      ORDER = IPAR(5)
      RADIUS = IPAR(7)
      NUM = IPAR(2)-IPAR(1)+IPAR(4)-IPAR(3)+2
      GAIN = RPAR(1)
      RON = RPAR(2)
      THRESHOLD = RPAR(3)

C
C-------COLLECT SKY VALUES
C
      DO 60 I=1,IX
        N = 0
        NS = 0

          DO 10 J=IPAR(1),IPAR(2)
            N = N+1
            WDT(N) = IN(I,J)
            WDTY(N) = J
10        CONTINUE
          CALL FCOSMIC(WDT,WDTY,COSMIC,N,GAIN,RON,THRESHOLD,RADIUS)
          DO 12 J=RADIUS+1,N-RADIUS
            IF (COSMIC(J) .EQ. 0) THEN
               NS = NS+1
               FLUX(NS) = WDT(J)
               Y(NS) = WDTY(J)
            ENDIF
 12      CONTINUE

         N=0
         DO 20  J=IPAR(3),IPAR(4)
            N = N+1
            WDT(N)=IN(I,J)
            WDTY(N) = J
 20      CONTINUE
         CALL FCOSMIC(WDT,WDTY,COSMIC,N,GAIN,RON,THRESHOLD,RADIUS)
         DO 22 J=RADIUS+1,N-RADIUS
            IF (COSMIC(J) .EQ. 0) THEN
               NS = NS+1
               FLUX(NS) = WDT(J)
               Y(NS) = WDTY(J)
            ENDIF
 22      CONTINUE
         NUM = NS
C
C the case of ORDER = 0 is treated seperately
C
          IF (ORDER.LE.0) THEN
            SUM = 0
            DO 30 J = 1,NUM
              SUM = SUM+FLUX(J)
30          CONTINUE
            SUM = SUM/NUM
            DO 40 J = 1,IY
              SKY(I,J) = SUM
40          CONTINUE
          ELSE
C
C fit polynomial and evaluate it for all rows in the column
C
            CALL LSORTH(Y,FLUX,A,B,S,W,NUM,CHISQ,ORDER)
            DO 50  J = 1,IY
              YY = J
              CALL POLY(YY,SKY(I,J),A,B,S,W,ORDER)
50          CONTINUE
          ENDIF
60      CONTINUE
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. 
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




C
C**********************************************************************
C      FCOSMIC: SEARCH FOR COSMICS DEFINED AS PIXELS WHOSE VALUES ARE 
C               THRESHOLD TIMES THE STANDARD DEVIATION ABOVE THE MEDIAN 
C               AVERAGE ON PIXELS WITHIN RADIUS (TAKING INTO ACCOUNT 
C               THE RON AND THE CONVERSION FACTOR) 
C
C       INPUT:  ARRAYS WDT (VALUES), WDTY (POSITIONS), N (NUMBER OF VALUES)
C               RON, CONVERSION FACTOR, THRESHOLD, RADIUS
C
C       OUTPUT: COSMIC = INTEGER ARRAY WITH SUSPECTED COSMIC = 1

        SUBROUTINE FCOSMIC(WDT,WDTY,COSMIC,N,GAIN,RON,THRESHOLD,RADIUS)


        INTEGER N, RADIUS, COSMIC(1),I,IS,IE,K,NC,J
C
        REAL WDT(1),WDTY(1),GAIN,RON,THRESHOLD,RONSQ
        REAL AVR,DIFF,ARR(4000),AA,VAL,SGA
C

        RONSQ = RON*RON
C
        DO 10 I = 1,N 
           COSMIC(I) = 0
 10     CONTINUE

        DO 20 I=RADIUS+1,N-RADIUS

           AVR = 0
           DIFF = 0

           IS = I-RADIUS
           IE = I+RADIUS

           NC = 0

           DO 30 K = IS,IE
              IF (K .NE. I ) THEN 
                 NC = NC+1
                 ARR(NC) = WDT(K)
              ENDIF
 30        CONTINUE

           DO 40 K = 2,NC
              AA = ARR(K)
              DO 45 J=K-1,1,-1
                 IF (ARR(J) .LE. AA ) GO TO 90
                 ARR(J+1) = ARR(J)
 45           CONTINUE
              J = 0
 90           ARR(J+1) = AA
 40        CONTINUE

           AVR = ARR(RADIUS)
           SGA = AVR/GAIN+RONSQ

           DIFF = WDT(I)-AVR
           IF (SGA.GT.0.) THEN
               VAL = THRESHOLD*SQRT(SGA)
           ELSE
               VAL = 0.
           ENDIF

           IF (ABS(DIFF) .GT. VAL) COSMIC(I) = 1

 20     CONTINUE


        RETURN

        END
