C @(#)double.for	19.1 (ES0-DMD) 02/25/03 13:25:33
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.IDENTIFICATION
C  subroutine DOUBLE         version 1         820521
C  A. Kruszewski             ESO Garching
C  modified for FX           version 1.1       870219
C  A. Kruszewski             Obs. de Geneve
C  modified for MSDOS        version 2.0       880618
C  A. Kruszewski             Warsaw U. Obs.
C.KEYWORDS
C  search, double
C.PURPOSE
C  analyzes a region of multiple detections for presence of
C  two components
C.INPUT/OUTPUT
C  input arguments
C  A           real*4 array      image buffer
C  JAPY        integer*4 array   array of pointers
C  JOFF        integer*4         line offset to buffer start
C  JBE         integer*4         last line in buffer
C  IXYU        integer*4 array   limits of used area
C  IHED        integer*4         half-edge in pixels of used subarray
C  BGRD        real*4            local sky background
C  AVI         real*4            average x-coordinate in pixels
C  AVJ         real*4            average y-coordinate in pixels
C  THETA       real*4            position angle of major axis
C  ISIZ        integer*4         maximum size of an object
C  ADLM        real*4            minimum allowed separation of
C                                components
C  YFCT        real*4            ratio of sizes of rectangle which is
C                                aligned with major axis and used to
C                                scan object's region
C  DVAL        real*4            minimum depth of valley between
C                                components expressed in units
C                                of central value of a brighter
C                                component
C  output arguments
C  IDPM        integer*4 array   contains informations about number
C                                of detected components and theirs
C                                positions
C-----------------------------------------------------------------------
      SUBROUTINE DOUBLE(A, JAPY, JOFF, JBE, IXYU, IHED, BGRD,
     &           AVI, AVJ, THETA, ISIZ, ADLM, YFCT, DVAL,
     &           TRSH, IDPM)
C
      IMPLICIT NONE
      INCLUDE  'MID_REL_INCL:INVENT.INC'
C
      REAL      A(1)
      INTEGER   JAPY(1)
      INTEGER   JOFF
      INTEGER   JBE
      INTEGER   IXYU(4)
      INTEGER   IHED
      REAL      BGRD
      REAL      AVI
      REAL      AVJ 
      REAL      THETA
      INTEGER   ISIZ 
      REAL      ADLM 
      REAL      YFCT
      REAL      DVAL
      REAL      TRSH
      INTEGER   IDPM(5)
C
      INTEGER   I1, I2
      INTEGER   IDLM, IDX
      INTEGER   IHBX, IHXY
      INTEGER   IHBY 
      INTEGER   II 
      INTEGER   J1, J2 
      INTEGER   JDY, JJ 
      INTEGER   K , KK , KK2
      INTEGER   LCTN(3)
      INTEGER   NX(0:MAXSCN)
C
      REAL      CTHT
      REAL      DI, DJ, DX, DY
      REAL      EXTR(3)
      REAL      SCNX(0:MAXSCN), STHT
C
C ******      Define dimensions of scanning rectangle IHBX and IHBY.
C
      IHBX = MIN( (ISIZ/2+2) , MAXSCN/2 )
      IHBY = MIN( IHED , INT( FLOAT(IHBX) * YFCT ) )
      CTHT = COS(THETA)
      STHT = SIN(THETA)
C
C ******      Define region containing scanning rectangle.
C
      IHXY = INT( SQRT( FLOAT( IHBX*IHBX + IHBY*IHBY ) ) ) + 1
      I1 = MAX( IXYU(1) , NINT(AVI)-IHXY )
      I2 = MIN( IXYU(3) , NINT(AVI)+IHXY )
      J1 = MAX( JOFF+1 , NINT(AVJ)-IHXY )
      J2 = MIN( JBE , NINT(AVJ)+IHXY )
      KK = MIN( 2*IHBX , MAXSCN )
      KK2 = KK / 2
      DO 10 K = 0 , KK
          SCNX(K)=0.0
          NX(K)=0
   10      CONTINUE
C
C ******      Perform scanning.
C
      DO 12 JJ = J1 , J2
          DO 11 II = I1 , I2
            DI = FLOAT(II) - AVI
            DJ = FLOAT(JJ) - AVJ
            DX = DI*CTHT + DJ*STHT
            DY = DJ*CTHT - DI*STHT
            IDX = NINT( DX + 0.001 )
            JDY = NINT( DY + 0.001 )
            IF ( IDX. GE. -KK2 .AND. IDX .LE. KK2 .AND.
     &                JDY .GE. -IHBY .AND. JDY .LE. IHBY ) THEN
                     K = IDX + KK2
                SCNX(K) = SCNX(K) + A(JAPY(JJ-JOFF)+II)
                NX(K) = NX(K) + 1
            ENDIF
   11          CONTINUE
   12      CONTINUE 
        DO 14 K = 0 , KK
            IF ( NX(K) .GT. 0 ) THEN
                  SCNX(K) = SCNX(K) / NX(K)
            ELSE IF ( K .GT. 0 ) THEN
                  SCNX(K) = SCNX(K-1)
            END IF
  14      CONTINUE
C
C ******      Search for extrema in obtained scan.
C
      IDLM = INT(ADLM)
      CALL EXTREM( SCNX(0) , KK2 , IDLM , EXTR , LCTN )
C
C ******      Write an output array IDPM.
C
      IF ( ABS( LCTN(1) ) .LT. KK2 ) THEN
          IF ( ABS( LCTN(1)-LCTN(2) ) .GE. INT(ADLM) .AND.
     &                               ABS( LCTN(2) ) .LT. KK2 .AND.
     &                 EXTR(2)-EXTR(3) .GT. DVAL*(EXTR(1)-BGRD)
     &                 .AND. EXTR(2)-EXTR(3) .GT. TRSH ) THEN
              IDPM(1) = 2
              IDPM(4) = NINT( LCTN(2)*CTHT )
              IDPM(5) = NINT( LCTN(2)*STHT )
          ELSE
              IDPM(1)=1
          ENDIF
          IDPM(2) = NINT( LCTN(1)*CTHT )
          IDPM(3) = NINT( LCTN(1)*STHT )
      ELSE IF ( ABS( LCTN(1) ) .EQ. KK2 ) THEN
          IF ( ABS( LCTN(2) ) .LT. KK2 ) THEN 
              IDPM(1) = 1
              IDPM(2) = NINT( LCTN(2)*CTHT )
              IDPM(3) = NINT( LCTN(2)*STHT )
          ELSE
              IDPM(1) = 0
          ENDIF
      ENDIF
C
      RETURN
C
      END
