C @(#)feet.for 17.1.1.1 (ES0-DMD) 01/25/02 17:18: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 SUBROUTINE FEET(MASK,ATEMP,IDEG,FLMSK,REJ,NMX,NMY,RXSTR,RYSTR, + RXSTP,RYSTP,NSX,NSY,SXSTR,SYSTR,SXSTP,SYSTP,ICONT, + RMEAN,SIGM1,SKEW1,RKURT,C,JER) C ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C.PURPOSE: Calculate the coefficients of a 2D polinomial fitting the C sky background,then give the first moments of the distribution C of O-C residuals C.ALGORITHM: Lleast squares 2-dim. polinomial fit C.INPUT: MASK : frame masking bright objects in the field C ATEMP : temporary buffer containing the rebinned input image C FLMSK : flag = .TRUE. if MASK exists C REJ : value assigned to rejected pixels in frame ATEMP C NMX,... usual descriptors of the frame MASK C NSX,... descriptors of the frame-buffer ATEMP C IDEG : polinomial degree C.OUTPUT: ICONT : number of processed 'good' points C RMEAN : mean value of residuals image-pol.fit for 'good' points C SIGM1 : standard deviation of residuals C SKEW1 : 3rd moment of the distribution of residuals (skewness) C RKURT : 4th moment of the distribution of residuals (kurtosis) C C : up to 21 coefficients of the polinomial surface C JER : error return ; 0=OK 1=interpolation failure C -------------------------------------------------------------------- IMPLICIT NONE C INTEGER NMX INTEGER NMY REAL MASK(NMX,NMY) REAL ATEMP(128,128) INTEGER IDEG LOGICAL FLMSK REAL REJ DOUBLE PRECISION RXSTR DOUBLE PRECISION RYSTR DOUBLE PRECISION RXSTP DOUBLE PRECISION RYSTP INTEGER NSX INTEGER NSY DOUBLE PRECISION SXSTR DOUBLE PRECISION SYSTR DOUBLE PRECISION SXSTP DOUBLE PRECISION SYSTP INTEGER ICONT REAL RMEAN REAL SIGM1 REAL SKEW1 REAL RKURT DOUBLE PRECISION C(21) INTEGER JER C INTEGER I, I1, I2 INTEGER I1MSK, I2MSK INTEGER ICX, ICY INTEGER IF1, IF2 INTEGER IPERC INTEGER IW1, IW2, IW3, IW4 INTEGER IX, IY INTEGER J, J1, J2 INTEGER J1MSK, J2MSK INTEGER NC, NC1 INTEGER NCV(5) REAL ANP REAL BUP(128) REAL BLW(128) REAL CONT REAL DCX, DCY REAL D1, D2, D3, D4 REAL DSQRT REAL POLIN REAL RES, RES2 REAL X1, X2 REAL Y1, Y2 REAL V REAL VUP1, VUP2 REAL VLW1, VLW2 DOUBLE PRECISION Z(64,10),A(21,22),B(21,22),CC(6) DOUBLE PRECISION W1,W2,W3,W4,DMEAN DOUBLE PRECISION XD,YD,SIGMA,SKEW,DKURT C EXTERNAL POLIN EXTERNAL SHORT,LONG C DATA NCV/3,6,10,15,21/ C C *** initialize the powers and arrays C DO 20 I = 1,64 Z(I,1) = I - 0.5 DO 10 J = 2,10 Z(I,J) = Z(I,J-1)*Z(I,1) 10 CONTINUE 20 CONTINUE C DO 40 I = 1,21 DO 30 J = 1,22 A(I,J) = 0. 30 CONTINUE 40 CONTINUE C DO 50 I = 1,21 C(I) = 0.D0 50 CONTINUE C C *** initialize variables C ICX = NSX/2 ICY = NSY/2 C C *** loop: load image values (testing mask and rejected pixels) C IF1 = 1 IF2 = 128* (NSY-1) + 1 C DO 90 J1 = 1,ICY J2 = NSY - J1 + 1 IY = ICY - J1 + 1 Y1 = SYSTR + SYSTP* (J1-1.) Y2 = SYSTR + SYSTP* (J2-1.) IF (FLMSK) THEN ! nearest pixel in mask frame J1MSK = (Y1-RYSTR)/RYSTP + 1.5 J2MSK = (Y2-RYSTR)/RYSTP + 1.5 C C *** mask is extrapolated outside C IF (J1MSK.LT.1) THEN J1MSK = 1 ELSE IF (J1MSK.GT.NMY) THEN J1MSK = NMY END IF IF (J2MSK.LT.1) THEN J2MSK = 1 ELSE IF (J2MSK.GT.NMY) THEN J2MSK = NMY END IF END IF C C *** read needed rows from temporary array C CALL ROW(ATEMP,BUP,IF1,NSX) CALL ROW(ATEMP,BLW,IF2,NSX) IF1 = IF1 + 128 IF2 = IF2 - 128 C DO 70 IX = 1,ICX I1 = ICX - IX + 1 I2 = ICX + IX VUP1 = 1. VUP2 = 1. VLW1 = 1. VLW2 = 1. C IF (FLMSK) THEN X1 = SXSTR + SXSTP* (I1-1.) X2 = SXSTR + SXSTP* (I2-1.) I1MSK = (X1-RXSTR)/RXSTP + 1.5 I2MSK = (X2-RXSTR)/RXSTP + 1.5 IF (I1MSK.LT.1) THEN I1MSK = 1 ELSE IF (I1MSK.GT.NMX) THEN I1MSK = NMX END IF IF (I2MSK.LT.1) THEN I2MSK = 1 ELSE IF (I2MSK.GT.NMX) THEN I2MSK = NMX END IF C C *** find out mask values - assumed floating p. C VUP1 = MASK(I1MSK,J1MSK) VUP2 = MASK(I2MSK,J1MSK) VLW1 = MASK(I1MSK,J2MSK) VLW2 = MASK(I2MSK,J2MSK) END IF C C *** a pixel is 'bad' if rejected via mask frame OR temp. buffer C IF (BUP(I1).EQ.REJ) VUP1 = -1. IF (BUP(I2).EQ.REJ) VUP2 = -1. IF (BLW(I1).EQ.REJ) VLW1 = -1. IF (BLW(I2).EQ.REJ) VLW2 = -1. C IF ((VUP1.GT.0.) .AND. (VUP2.GT.0.) .AND. + (VLW1.GT.0.) .AND. (VLW2.GT.0)) THEN D1 = BUP(I1) D2 = BUP(I2) D3 = BLW(I1) D4 = BLW(I2) CALL SHORT(D1,D2,D3,D4,IX,IY,Z,A) ELSE IW2 = 1 IF (VUP2.LE.0.) THEN BUP(I2) = 0. IW2 = 0 END IF IW1 = 1 IF (VUP1.LE.0.) THEN BUP(I1) = 0. IW1 = 0 END IF IW4 = 1 IF (VLW2.LE.0.) THEN BLW(I2) = 0. IW4 = 0 END IF IW3 = 1 IF (VLW1.LE.0.) THEN BLW(I1) = 0. IW3 = 0 END IF W1 = IW1 + IW2 + IW3 + IW4 IF (W1.EQ.0.0D0) GO TO 60 W2 = IW2 + IW4 - IW1 - IW3 W3 = IW1 + IW2 - IW3 - IW4 W4 = IW2 + IW3 - IW1 - IW4 D1 = BUP(I1) D2 = BUP(I2) D3 = BLW(I1) D4 = BLW(I2) CALL LONG(W1,W2,W3,W4,D1,D2,D3,D4,IX,IY,Z,A) END IF 60 CONTINUE 70 CONTINUE 90 CONTINUE C C *** fill matrix "a" using symmetries C CALL FLMTX(A) ! processed points ANP = 4.0D0*A(1,1) DO 110 I = 1,20 I1 = I + 1 DO 100 J = I1,21 A(J,I) = A(I,J) 100 CONTINUE 110 CONTINUE C C *** calculate coefficients --->array c C NC = NCV(IDEG) NC1 = NC + 1 DO 130 I = 1,NC DO 120 J = 1,NC B(I,J) = A(I,J) 120 CONTINUE B(I,NC1) = A(I,22) 130 CONTINUE C CALL SYST(NC,B,C,JER) !interpolation failure IF (JER.NE.0) RETURN C C *** calculate parameters of distribution of residuals C initialize variables C DCX = (NSX-1.)/2. DCY = (NSY-1.)/2. DMEAN = 0.D0 SIGMA = 0.D0 SKEW = 0.D0 DKURT = 0.D0 ICONT = 0 C C *** loop : calculate parameters C DO 170 J = 1,NSY IY = J - 1 YD = DCY - IY CALL COEFY(YD,C,CC,IDEG) Y1 = SYSTR + SYSTP*FLOAT(IY) IF (FLMSK) THEN J1MSK = (Y1-RYSTR)/RYSTP + 1.5 IF (J1MSK.LT.1) THEN J1MSK = 1 ELSE IF (J1MSK.GT.NMY) THEN J1MSK = NMY END IF END IF C DO 150 I = 1,NSX IX = I - 1 XD = IX - DCX V = 1. IF (FLMSK) THEN X1 = SXSTR + SXSTP*FLOAT(IX) I1MSK = (X1-RXSTR)/RXSTP + 1.5 IF (I1MSK.LT.1) THEN I1MSK = 1 ELSE IF (I1MSK.GT.NMX) THEN I1MSK = NMX END IF V = MASK(I1MSK,J1MSK) END IF IF (ATEMP(I,J).EQ.REJ) V = -1 IF (V.LE.0) GO TO 140 ICONT = ICONT + 1 RES = ATEMP(I,J) - POLIN(XD,CC,IDEG) RES2 = RES*RES DMEAN = DMEAN + RES SIGMA = SIGMA + RES2 SKEW = SKEW + RES2*RES DKURT = DKURT + RES2*RES2 140 CONTINUE 150 CONTINUE 170 CONTINUE C CONT = FLOAT(ICONT) DMEAN = DMEAN/DBLE(CONT) SIGMA = DSQRT(SIGMA/DBLE(CONT)) SKEW = SKEW/ (DBLE(CONT)* (SIGMA**3)) DKURT = DKURT/ (DBLE(CONT)* (SIGMA**4)) RMEAN = SNGL(DMEAN) SIGM1 = SNGL(SIGMA) SKEW1 = SNGL(SKEW) RKURT = SNGL(DKURT) IPERC = (CONT/ (NSX*NSY))*100. C RETURN END