C LAST UPDATE: 
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C**************************
C**     LFN_SPOT_CG      **
C**************************
C
CD-Routine: Get Spot Position - LFN_SPOT_CG
C This routine finds a spot position by determining the centre of gravity 
C of intensities lying within a defined box allowing for a constant 
C intensity threshold. It is based on a routine from the Daresbury Laboratory 
C Laue suite program GENLAUE as of September 1990)
CD-end
CD-Fortran:
      SUBROUTINE LFN_SPOT_CG (IXC, IYC, IXSIZ, IYSIZ, IDATA, IWORK,
     +                        ITHRESH, IMG, ITYPE, IORDER, NFOFF, 
     +                        NXRAST, NYRAST, NOVER, CGX, CGY, IERR)
CD-end
      INCLUDE 'implicit_none.finc'
C
C Purpose:    Find the centre of gravity of intensities lying within a
C             defined box allowing for a constant intensity threshold.
C
C Author:     John W. Campbell, 
C             (Based on a routine from the Daresbury Laboratory Laue
C              suite program GENLAUE as of September 1990)
C
C Arguments:
C
      INTEGER IXC
      INTEGER IYC
      INTEGER IXSIZ
      INTEGER IYSIZ
      INTEGER IDATA(-IYSIZ:IYSIZ,-IXSIZ:IXSIZ)
      INTEGER IWORK(*)
      INTEGER ITHRESH
      INTEGER IMG(*)
      INTEGER ITYPE
      INTEGER IORDER
      INTEGER NFOFF
      INTEGER NXRAST
      INTEGER NYRAST
      INTEGER NOVER
      REAL CGX
      REAL CGY
      INTEGER IERR
CD-Parameters:
C  IXC     i (R)   X pixel position for box centre (numbered from 1 up)
C  IYC     i (R)   Y pixel position for box centre
C  IXSIZ   i (R)   X box parameter (Box size will be (2*IXSIZ+1)*(2*IYSIZ+1)
C  IYSIZ   i (R)   Y box parameter (See previous parameter)
C  IDATA   i (W)   Array to hold read in data. Must be dimensioned to 
C                  at least (2*IXSIZ+1) * (2*IYSIZ+1) elements
C  IWORK   i (R)   Work array (dimensioned to at least larger of
C                  (2*IXSIZ+1) * (2*IYSIZ+1) elements
C  ITHRESH i (R)   Intensity threshold
C  IMG     i (R)   In memory image data
C  ITYPE   i (R)   Data type in stored image 1=unsigned byte,
C                  2=unsigned two-byte
C  IORDER  i (R)   Image data order flag (1-8) 
C  NFOFF   i (R)   Offset in pixels between lines in stored image
C  NXRAST  i (R)   The number of x rasters in the image
C  NYRAST  i (R)   The number of y-rasters in the image
C  NOVER   i (W)   Returns the number of pixels over the threshold value
C  CGX     r (W)   Returns the x coordinate (in raster units) of the 
C                  c of g
C  CGY     r (W)   Returns the y coordinate (in raster units) of the 
C                  c of g
C  IERR    i (W)   Error flag =0 OK, 
C                        =-1 flag error, =-2 byte handling not available
C                         =1 box partly out of range
C                         =2 box completely out of range of image.
CD-end
C Declarations
C ============
C
C GLOBALS:
C

C
C LOCALS:
C
      INTEGER IX
      INTEGER IY
      INTEGER INTENS
      REAL SUMX
      REAL SUMY
      REAL SUM_INTENS   

C     INTEGER IX      !Loop counter over x points
C     INTEGER IY      !Loop counter over y points
C     INTEGER INTENS  !Pixel intensity
C     REAL SUMX       !Sum for cgx calcn.
C     REAL SUMY       !Sum for cgy calcn.
C     REAL SUM_INTENS !Sum of intensities over threshold
C
C-------------------------------------------------------------------------------
C
C====== Initialisations
C
      NOVER = 0
      CGX = IXC
      CGY = IYC
C
C====== Read spot
C
      CALL LFN_IMG_RECT (IMG, ITYPE, IORDER, NFOFF, 
     +                   IXC-IXSIZ, IYC-IYSIZ, 2*IXSIZ+1, 2*IYSIZ+1,
     +                   NXRAST, NYRAST, 0, IDATA, IWORK, IERR)
      IF (IERR.LT.0.OR.IERR.GT.1) RETURN
C
C====== Find Centre of Gravity
C
      SUMX = 0.0
      SUMY = 0.0
      SUM_INTENS = 0.0
      DO 200 IX = -IXSIZ, IXSIZ
         DO 100 IY = -IYSIZ, IYSIZ
            INTENS = MAX (0, IDATA(IY,IX) - ITHRESH)
            IF (INTENS.GT.0) NOVER = NOVER + 1
            SUMX = SUMX + IX * INTENS
            SUMY = SUMY + IY * INTENS
            SUM_INTENS = SUM_INTENS + INTENS
100      CONTINUE
200   CONTINUE
      IF (SUM_INTENS.EQ.0.0) RETURN
      CGX = IXC + SUMX / SUM_INTENS
      CGY = IYC + SUMY / SUM_INTENS
      RETURN
      END
