C @(#)spskyfit.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 @(#)spskyfit.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 @(#)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