C LAST UPDATE: 
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C******************************
C**     LFN_SPOT_RD          **
C******************************
C
CD-Routine: Get Radially Masked Spot Position - LFN_SPOT_RD
C This routine finds a spot position by determining the centre of gravity 
C of intensities lying within a radially masked ellipse allowing for a 
C constant intensity threshold. It is based on code from the program GENLAUE 
C including the radial masking code of T.J. Greenhough and A.K. Shrive 
C (University of Keele)
CD-end
CD-Fortran:
      SUBROUTINE LFN_SPOT_RD (IPACK, IPLATE, IXC, IYC, 
     +            IXSIZ, IYSIZ, IDATA, 
     +            IWORK, ITHRESH, IMG, ITYPE, IORD, NFOFF, 
     +            NXRAST, NYRAST,
     +            X_CEN, Y_CEN, NOVER, 
     +            CGX, CGY, IERR)
CD-end
      INCLUDE 'implicit_none.finc'
C
C Purpose:    Find the centre of gravity of intensities lying within a
C             radially masked ellipse allowing for a constant intensity 
C             threshold.
C
C Author:     Based on code from GENLAUE including radial masking code
C             of T.J. Greengough and A.K. Shrive (Univ. of Keele)
C
C Arguments:
C
      INTEGER IPACK
      INTEGER IPLATE
      INTEGER IXC
      INTEGER IYC
      INTEGER IXSIZ
      INTEGER IYSIZ
      INTEGER IDATA(-IYSIZ:IYSIZ,-IXSIZ:IXSIZ)
      INTEGER IWORK(*)
      INTEGER ITHRESH
      INTEGER IMG(*)
      INTEGER ITYPE
      INTEGER IORD
      INTEGER NFOFF
      INTEGER NXRAST
      INTEGER NYRAST
      REAL X_CEN
      REAL Y_CEN
      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 at 
C                  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  IORD    i (R)   Image data order flag (1-8) 
C  NFOFF   i (R)   Offset in pixels between lines in stored image
C  NXRAST  i (R)   No. of x-rasters
C  NYRAST  i (R)   No. y-rasters
C  X_CEN   r (R)   Corrected 'x' centre position (rasters)
C  Y_CEN   r (R)   Corrected 'y' centre position (rasters)
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
C parameter statements
C

C
C GLOBALS:
C

C
C LOCALS:
C
      INTEGER IXM
      INTEGER IYM
      INTEGER XPN
      INTEGER YPN
      INTEGER IMARK
      INTEGER IX
      INTEGER IY
      INTEGER INTENS
      INTEGER IFLAG
      INTEGER KERR
      INTEGER ND
      REAL SPOTFAC
      REAL SPOTL
      REAL SPOTW
      REAL CTOF
      REAL MM_RAST_X
      REAL MM_RAST_Y
      REAL SUMX
      REAL SUMY
      REAL SUM_INTENS   
      REAL PI
      REAL CMASK
      REAL RADIUS
      REAL AMAJOR
      REAL BMINOR
      REAL THETA
      REAL XPRIME
      REAL YPRIME
      REAL ELLIP
      CHARACTER*2 ERRSTR

C     INTEGER IX      !Loop counter over x points
C     INTEGER IY      !Loop counter over y points
C     INTEGER INTENS  !Pixel intensity
C     INTEGER IXM     !x offset of ellipse centre from film centre in rasters
C     INTEGER IYM     !y offset of ellipse centre from film centre in rasters
C     INTEGER XPN     !x offset of current pixel from film centre in rasters
C     INTEGER YPN     !y offset of current pixel from film centre in rasters
C     INTEGER IMARK   !Flag =1 include pixel as part of spot, =0 do not
C     INTEGER IFLAG   !Parameter status flag
C     INTEGER KERR    !Error flag
C     INTEGER ND      !No. of decimal places flag
C     REAL SPOTFAC    !Spot expansion factor
C     REAL SPOTL      !Spot length
C     REAL SPOTW      !Spot width
C     REAL CTOF       !Crystral-film distance
C     REAL MM_RAST_X  !Conversion factor mm to x-rasters
C     REAL MM_RAST_Y  !Conversion factor mm to y-rasters
C     REAL SUMX       !Sum for cgx calcn.
C     REAL SUMY       !Sum for cgy calcn.
C     REAL SUM_INTENS !Sum of intensities over threshold
C     REAL PI         !Pi
C     REAL CMASK      !1/2 streak length (uncorrected) in rasters
C     REAL RADIUS     !Distance of ellipse centre from film centre (rasters)
C     REAL AMAJOR     !Ellipse major axis (corrected) in rasters
C     REAL BMINOR     !Ellipse minor axis in rasters
C     REAL THETA      !Angle of ellipse major axis wrt. xf (radians)
C     REAL XPRIME     !x' coordinate of current point wrt ellipse axes (rasters)
C     REAL YPRIME     !y' coordinate of current point wrt ellipse axes (rasters)
C     REAL ELLIP      !Value of ellipse function x**2/a**+y**2/b**2
C     CHARACTER*2 ERRSTR !Error string (dummy)
C
C-------------------------------------------------------------------------------
C
C====== Initialisations
C
      NOVER = 0
      CGX = IXC
      CGY = IYC
      CALL LDM_GET_SPOTF (IPACK, IPLATE, SPOTFAC, IFLAG, ND, 
     +                    KERR, ERRSTR)
      CALL LDM_GET_CTOF (IPACK, IPLATE, CTOF, IFLAG, ND,
     +                   KERR, ERRSTR)
      CALL LDM_GET_SPOTL (IPACK, IPLATE, SPOTL, IFLAG, ND,
     +                    KERR, ERRSTR)
      CALL LDM_GET_SPOTW (IPACK, IPLATE, SPOTW, IFLAG, ND,
     +                    KERR, ERRSTR)
      CALL LDM_MMRAST (IPACK, IPLATE, MM_RAST_X, MM_RAST_Y)
C
C====== Read spot
C
      CALL LFN_IMG_RECT (IMG, ITYPE, IORD, 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
C
C====== Get offsets in pixels of box centre from pattern centre
C
      IXM = NINT(FLOAT(IXC)-X_CEN)
      IYM = NINT(FLOAT(IYC)-Y_CEN)
C
C====== Set up parameters for radial masking calculation
C
      PI = 4.0*ATAN(1.0)
      CMASK = SPOTL*MM_RAST_X/2.0
      RADIUS = SQRT(FLOAT(IXM**2)+FLOAT(IYM**2))
      AMAJOR = CMASK*(1.0+SPOTFAC*RADIUS/
     +         (CTOF*MM_RAST_X))
      BMINOR = SPOTW*MM_RAST_X/2.0
      IF (IXM.EQ.0) THEN
         THETA = PI/2.0
         IF (IYM.LT.0) THEN
            THETA = 1.5*PI
         END IF
         GO TO 50
      END IF

      THETA = ATAN(FLOAT(IYM)/FLOAT(IXM))
      IF (IXM.LT.0) THEN
         THETA = THETA + PI
      ELSE IF (IYM.LT.0) THEN
         THETA = THETA + 2.0*PI
      END IF

50    CONTINUE
C
C====== Loop through pixels within box and see if each is within ellipse
C
      DO 200 IX = -IXSIZ, IXSIZ
         DO 100 IY = -IYSIZ, IYSIZ
C
C====== Get offsets in pixels of current point from pattern centre
C
            XPN = NINT(FLOAT(IXC+IX)-X_CEN)
            YPN = NINT(FLOAT(IYC+IY)-Y_CEN)

            IMARK = 0
            XPRIME = (FLOAT(XPN)-IXM)*COS(THETA)
            XPRIME = XPRIME + (FLOAT(YPN)-IYM)*SIN(THETA)
            YPRIME = (FLOAT(YPN)-IYM)*COS(THETA)
            YPRIME = YPRIME - (FLOAT(XPN)-IXM)*SIN(THETA)
            ELLIP = (XPRIME**2/AMAJOR**2) + (YPRIME**2/BMINOR**2)
            IF (ELLIP.LE.1.0) IMARK = 1
C
C====== Include pixel intensity if within ellipse
C
            IF (IMARK.EQ.1) THEN
               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
            END IF
100      CONTINUE
200   CONTINUE
      IF (SUM_INTENS.EQ.0.0) RETURN
      CGX = IXC + SUMX / SUM_INTENS
      CGY = IYC + SUMY / SUM_INTENS
      RETURN
      END
