C @(#)cmpsub.for 17.1.1.1 (ES0-DMD) 01/25/02 17:15:39 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+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C.IDENT: subroutine CMPSUB version 1.1 870315 C A. Kruszewski Obs. de Geneve C undersampled case version 2.0 880826 C A. Kruszewski Warsaw U. Obs. C.PURPOSE: Subtracts stellar components in a subframe copied from image frame. C.INPUT/OUTPUT C Input parameters C C A real*4 array input image array C JAPY integer*4 array pointers to imege lines C IBUF integer*4 limits of image buffer C I integer*4 x-coordinate of object C J integer*4 y-coordinate of object C L0 integer*4 offset to catalog objects C L1 integer*4 last object in catalog C LW integer*4 last object written to file C L integer*4 identification number of object C M integer*4 actuall number of objects C LSTP integer*4 array manage linked lists of pointers C NREG integer*4 number of regional linked lists C LHED integer*4 half edge size of subarray C TRSH real*4 limiting threshold C APSF real*4 array one-dimensional point spread function C FPSF real*4 array two-dimensional point spread function C IPSF integer*4 array pointers to two-dimensional p.s.f. C LPXL integer*4 number of 2-d pixels C LSBP integer*4 number of subpixels C NCAT integer*4 array integer parameters array C PMTR real*4 array real parameters array C PRCT real*4 array catalog of one-dimensional profiles C NP integer*4 dimension of catalog buffer C NPAS integer*4 actuall iteration number C HCUT real*4 saturation level C Output parameters C AS real*4 array data subarray after subtraction C JAPYS integer*4 array pointers to subarray lines C IBUFS integer*4 array limits of data in subarray C---------------------------------------------------------------------- SUBROUTINE CMPSUB (A, JAPY, IBUF, I, J, L0, L1, LW, & L, M, NREG, LSTP, LHED, TRSH, APSF, FPSF, & IPSF, LPXL, LSBP, NCAT, PMTR, PRCT, NPAS, & AS, JAPYS, IBUFS, HCUT, MAS ) C IMPLICIT NONE INCLUDE 'MID_REL_INCL:INVENT.INC/NOLIST' C INTEGER NP , NREG INTEGER JAPY(1) , IBUF(4), NCAT(NIPAR,MAXCNT) INTEGER JAPYS(2*MAXSUB+1), IBUFS(4) INTEGER LSTP(0:4, 0:NREG) INTEGER NCT(NIPAR) INTEGER MAS((-MAXSUB):MAXSUB,(-MAXSUB):MAXSUB) INTEGER I, IPSF(1), J, L, L0, L1, LHED INTEGER LPXL, LSBP, LW INTEGER M INTEGER NPAS INTEGER MDIS INTEGER MPXL INTEGER MSBP INTEGER IL1, IL2 INTEGER JL1, JL2 INTEGER MXS INTEGER ITMP INTEGER IS, IS1, IS2 INTEGER JOFF, JOFFS INTEGER JOF, JOFS INTEGER IL, IN INTEGER K, KK, KK2 INTEGER LK, LN, LZ, LD INTEGER JN, JL INTEGER IDIF INTEGER JDIF INTEGER IOPC INTEGER IINT C REAL AS((2*MAXSUB+1)**2) REAL A(1), APSF(0:MAXSUB) REAL CINT REAL CINT1 REAL FPSF(1) REAL HCUT REAL PRCT( 0:MAXSUB, MAXCNT), PMTR(NRPAR,MAXCNT) REAL PMT(NRPAR), PCT(0:MAXSUB) REAL SCALE REAL TRSH C C LOGICAL CENTER LOGICAL UNDER LOGICAL DONE C NP = MAXCNT MDIS = MAXSUB + LHED IF ( LPXL .GT. 0 .OR. LSBP .GT. 0 ) THEN UNDER = .TRUE. MPXL = 2 * LPXL + 1 MSBP = 2 * LSBP + 1 ELSE UNDER = .FALSE. ENDIF LHED = MIN( MAXSUB , LHED ) SCALE = 9.0 / ( 1.0 + 8.0*APSF(1) ) CINT = 0.0 C C ****** Prepare subarrays. C IL1 = MAX( I-LHED , IBUF(1) ) IL2 = MIN( I+LHED , IBUF(3) ) JL1 = MAX( J-LHED , IBUF(2) ) JL2 = MIN( J+LHED , IBUF(4) ) IBUFS(1) = IL1 - I IBUFS(2) = JL1 - J IBUFS(3) = IL2 - I IBUFS(4) = JL2 - J MXS = 2*MAXSUB + 1 IS1 = IBUFS(2) IS2 = IBUFS(4) DO 5 IS = IS1 , IS2 ITMP = IS - IS1 + 1 JAPYS(ITMP) = (ITMP-1) * MXS + MAXSUB + 1 5 CONTINUE JOFF = IBUF(2) - 1 JOFFS = IBUFS(2) - 1 DO 10 JL = JL1 , JL2 JOF = JAPY(JL-JOFF) JOFS = JAPYS(JL-J-JOFFS) DO 20 IL = IL1 , IL2 AS(JOFS+IL-I) = A(JOF+IL) 20 CONTINUE 10 CONTINUE C C *** Start subtracting contributions from close components. C K = 0 31 CONTINUE CALL GETLST( L , L0 , L1 , MDIS , NREG , & LSTP , NCAT , PMTR , PRCT , K , & NCT , PMT , PCT , DONE ) IF ( .NOT. DONE ) GOTO 30 C C *** Ignore the L-th object. C IF ( K .EQ. L ) GOTO 32 C C *** Look for close neighbour. C *** Check first if K-th object is in buffer. C c IF ( K .LE. L0 .OR. K .GT. L1 ) THEN c READ( ISF , REC=K ) NCT , PMT , PCT c ELSE LK = K - L0 DO 60 KK = 1 , NIPAR NCT(KK) = NCAT(KK,LK) 60 CONTINUE DO 70 KK = 1 , NRPAR PMT(KK) = PMTR(KK,LK) 70 CONTINUE KK2 = MIN( NCT(5) , MAXSUB ) DO 80 KK = 0 , KK2 PCT(KK) = PRCT(KK,LK) 80 CONTINUE c ENDIF IF ( PMT(2) .LE. TRSH ) THEN C C *** Too faint objects are not considered. C GOTO 32 ENDIF IN = NCT(1) JN = NCT(2) LN = NCT(5) IDIF = I-IN JDIF = J-JN LD = LN + LHED IF ( ABS(IDIF) .LE. LD .AND. ABS(JDIF) .LE. LD ) THEN C C *** Found one! C IOPC = -1 CALL STARSA( IOPC , AS , JAPYS , IBUFS , I , & J , LPXL , LSBP , NCT , PMT , & PCT , APSF , FPSF , NPAS , SCALE , & CINT1 ) IF ( CINT1 .GT. CINT ) THEN CINT = CINT1 IINT = K ENDIF ENDIF 32 CONTINUE IF ( DONE ) GOTO 31 30 CONTINUE C C *** Write down the number of the brightest stellar component C *** if it contribute more than 10 % to the central pixel. C LZ = L - L0 IF ( CINT .GT. 0.1*PMTR(2,LZ) ) THEN NCAT(10,LZ) = IINT ELSE NCAT(10,LZ) = 0 ENDIF C C *** Restore values of saturated points. C DO 40 JL = JL1 , JL2 JOFS = JAPYS(JL-J-JOFFS) DO 50 IL = IL1 ,IL2 IF ( MAS(IL-I,JL-J) .EQ. -1 ) THEN AS(JOFS+IL-I) = HCUT ENDIF 50 CONTINUE 40 CONTINUE C RETURN C END