C LAST UPDATE:
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C************************
C**      LRL_OVLP      **
C************************
C
CD-Routine: Find and Mark Spatial Overlaps - LRL_OVLP
C This routine is used to determine an mark the spatial overlaps in a
C Laue Reflection List. A number of ways of specifying the overlap
C criterion are available.
CD-end
CD-Fortran:
      SUBROUTINE  LRL_OVLP (IOVTYP, IERR, ERRSTR)
CD-end
      INCLUDE 'implicit_none.finc'
C
C Purpose: Sort reflexions in LRL to find and mark spatial overlaps
C
C          (Used as a service routine for LRL_GEN but may be
C           called directly if desired)
C
C Author:  Based on version from GENLAUE by Clifton et al.
C
C Arguments:
C
      INTEGER IOVTYP
      INTEGER IERR
      CHARACTER*(*) ERRSTR
CD-Parameters:
C IOVTYP    i (R)   Flag for calculation of spatial overlaps
C                        = 1, analytical calculation (spot + unused
C                             pixel + border) v spot
C                        = 2, analytical calculation (spot + border) v spot
C                        = 3, analytical calculation spot v spot
C                        = 4, mask method (spot + unused pixel + border) 
C                             v spot
C                        = 5, mask method (spot + border) v spot
C                        = 6, mask method spot v spot
C
C IERR      i (W)  Error flag = 0, OK
C                             = 1, Invalid overlaps type flag
C                             = 2, Too many spots in a strip for overlaps 
C                                  check
C ERRSTR    c (W)  Error message string (max 80 characters)
CD-end
C====== Parameter statements:
C
      INCLUDE 'lrl_maxss.finc'
C
C====== GLOBALS:
C
      INCLUDE 'lrl_nodmax.finc'
      INCLUDE 'lrl_data.finc'
C
C====== FUNCTION TYPES:
C 
      LOGICAL LRL_OVQ
      LOGICAL LRL_OVM
      INTEGER NBITST
      EXTERNAL LRL_BDAT
C
C====== LOCALS:
C
      LOGICAL    RAD_MASK
      LOGICAL    OVLP
      INTEGER    IFLAG
      INTEGER    ND
      INTEGER    I
      INTEGER    J
      INTEGER    II
      INTEGER    ISPOT
      INTEGER    ISPOT1
      INTEGER    ISPOT2
      INTEGER    NOD_IDX
      INTEGER    N
      INTEGER    NUM_STRIPS
      INTEGER    STRIP
      INTEGER    NH
      INTEGER    NK
      INTEGER    NL
      REAL       XMIN
      REAL       XMAX
      REAL       X_STEP
      REAL       XFD
      REAL       YFD
      REAL       X1
      REAL       Y1
      REAL       X2
      REAL       Y2
      REAL       SPOT_L
      REAL       SPOT_W
      REAL       SPOT_B
      REAL       SPOTFAC
      REAL       MM_RAST_X
      REAL       MM_RAST_Y
      REAL       RASTOMM
      REAL       CTOF
      REAL       RMAX
      REAL       EPSILON
      REAL       DELTAM
      REAL       DELTA2
      REAL       EPSILON2
      REAL       D2
      REAL       X_DIFF
      REAL       Y_DIFF
C
C------------------------------------------------------------------------------
C
      IERR = 0
      ERRSTR = ' '
      IF (NUM_SPOTS.LT.2) RETURN
      DO 10 I = 1, NUM_SPOTS
         CALL STBITS(FLAGS(I),0,4,0)
10    CONTINUE
      IF (IOVTYP.LT.0.OR.IOVTYP.GT.6) THEN
         IERR = 1
         ERRSTR = 'Invalid overlaps calculation flag'
         RETURN
      END IF
      OVLAP_STATUS = IOVTYP
C
C====== Get delta value or ellipse max. dimension if radial masking
C
      CALL LDM_MMRAST (IPK, IPL, MM_RAST_X, MM_RAST_Y)
      RASTOMM = 1.0/MM_RAST_X

      RAD_MASK = .FALSE.
      CALL LDM_GET_SPOTW (IPK, IPL, SPOT_W, IFLAG, ND, 
     +                    IERR, ERRSTR)
      IF (SPOT_W.GT.0.0) RAD_MASK = .TRUE.
      CALL LDM_GET_SPOTL (IPK, IPL, SPOT_L, IFLAG, ND, 
     +                    IERR, ERRSTR)
      IF (RAD_MASK) THEN
         DELTAM = SPOT_L
         CALL LDM_GET_SPOTF (IPK, IPL, SPOTFAC, IFLAG, ND, 
     +                       IERR, ERRSTR)
         CALL LDM_GET_RMAX (RMAX, IFLAG, ND)
         CALL LDM_GET_CTOF (IPK, IPL, CTOF, IFLAG, ND,
     +                         IERR, ERRSTR)
         IF (SPOTFAC.GT.0.0) THEN
            DELTAM = DELTAM*(1.0+SPOTFAC*RMAX/CTOF)       
         END IF
         IF (SPOT_W.GT.DELTAM) DELTAM = SPOT_W
         CALL LDM_GET_SPOTB (IPK, IPL, SPOT_B, IFLAG, ND,
     +                       IERR, ERRSTR)
         DELTAM = DELTAM + (1.0+SPOT_B)*RASTOMM
      ELSE
         CALL LDM_GET_SPOTD (IPK, IPL, DELTAM, IFLAG, ND,
     +                       IERR, ERRSTR)
      END IF
C
C====== Adjust delta to a minimum of 10 micron
C
      DELTAM = MAX(DELTAM,0.01)
      DELTA2 = DELTAM**2
      CALL LDM_GET_EPS (EPSILON, IFLAG, ND)
      EPSILON2 = EPSILON**2

      NUM_SPAT = 0
      NUM_CLOS = 0

C
C====== Sort spot list on 'yfd' 
C
      CALL LRL_SORT (4)
C
C====== To allow for fact that some strips will contain more spots
C====== than others, fill strips half full on average.
C====== Besides this, it seems to be necessary to take a few extra strips.
C====== (Currently 2 extra)
C====== Changed again for Janos  28 Feb 87 (/3 rather than /2)
C
      NUM_STRIPS = NUM_SPOTS/(MAX_SPOTS_IN_STRIP/3.0) + 1 + 1
C
C====== Loop through the strips
C
      X_STEP = (XEND-XSTART)/NUM_STRIPS
      XMIN = XSTART - DELTAM - X_STEP
      XMAX = XMIN + X_STEP + 2.0*DELTAM
      DO 70 STRIP = 1 , NUM_STRIPS
         N = 0
         XMIN = XMIN + X_STEP
         XMAX = XMAX + X_STEP
         ISPOT = 0
C
C====== Read through & sort reflexions
C
 30      ISPOT = ISPOT + 1
         IF (ISPOT.GT.NUM_SPOTS) THEN
C
C======  O.K. now we have finished reading...
C======  ...so sort the spots within the strip
C
            IF (N.EQ.0) GO TO 70
C
C====== Loop over each spot within the strip
C
            IF (N.GT.1) THEN
               DO 60 I = 1 , N - 1
                  DO 50 J = I + 1 , N
                     ISPOT1 = IWORK1(I)
                     ISPOT2 = IWORK1(J)
                     X1 = XF_D(ISPOT1)
                     Y1 = YF_D(ISPOT1)
                     X2 = XF_D(ISPOT2)
                     Y2 = YF_D(ISPOT2)
C
C====== If y difference exceeds delta then finish spot i.
C
                     Y_DIFF = ABS(Y1-Y2)
                     IF (Y_DIFF.GT.DELTAM) GO TO 60
C
C====== If y difference is possible then check x difference.
C
                     X_DIFF = ABS(X1-X2)
                     IF (X_DIFF.LE.DELTAM) THEN
C
C====== If both are possible then check precise difference.
C====== (Further approximation just if radial masking in use)
C
                        D2 = X_DIFF**2 + Y_DIFF**2
                        IF (D2.LE.DELTA2) THEN
C
C====== An overlap of some sort ... set both non-measurable and count.
C
                           IF (D2.LE.EPSILON2) THEN
C
C====== Some spot with spatial overlap, but too close to separate.
C====== Do not count twice! (set spatial overlap & too-close flags)
C
                              IF (NBITST(FLAGS(ISPOT1),1,1).EQ.0) THEN
                                 CALL STBITS (FLAGS(ISPOT1),0,2,3)
                                 NUM_CLOS = NUM_CLOS + 1
                              END IF

                              IF (NBITST(FLAGS(ISPOT2),1,1).EQ.0) THEN
                                 CALL STBITS (FLAGS(ISPOT2),0,2,3)
                                 NUM_CLOS = NUM_CLOS + 1
                              END IF

                           ELSE
C
C====== Some spot with spatial overlap, suitable for separation.
C====== (Further check needed if radial masking in use)
C
C====== Radial masking code makes approximation that overlapping
C====== Spots have the same size of ellipse and are in the same
C====== orientation
C
                              OVLP = .TRUE.
                              IF (RAD_MASK) THEN
                                 IF (IOVTYP.LE.3) THEN
                                    OVLP = LRL_OVQ (SPOT_L, SPOT_W,
     +                                     SPOTFAC, SPOT_B, CTOF,
     +                                     RASTOMM, X1, Y1, X2, Y2,
     +                                     IOVTYP)
                                 ELSE
                                    OVLP = LRL_OVM (SPOT_L, SPOT_W,
     +                                     SPOTFAC, SPOT_B, CTOF,
     +                                     MM_RAST_X, MM_RAST_Y, 
     +                                     X1, Y1, X2, Y2, IOVTYP-3)
                                 END IF
                              END IF

                              IF (OVLP) THEN
                                 IF (NBITST(FLAGS(ISPOT1),0,1).EQ.0) 
     +                              THEN
                                    CALL STBITS (FLAGS(ISPOT1),0,1,1)
                                    NUM_SPAT = NUM_SPAT + 1
                                 END IF
   
                                 IF (NBITST(FLAGS(ISPOT2),0,1).EQ.0) 
     +                              THEN
                                    CALL STBITS (FLAGS(ISPOT2),0,1,1)
                                    NUM_SPAT = NUM_SPAT + 1
                                 END IF
                              END IF
                           END IF
                         END IF
                     END IF

 50               CONTINUE
C
C====== End of loop through the strip of spots
C
 60            CONTINUE
            END IF
         ELSE
            II = INDX(ISPOT)
            XFD = XF_D(II)
            YFD = YF_D(II)
            IF (XFD.GE.XMIN .AND. XFD.LE.XMAX) THEN
C
C====== Add spot to strip, whilst checking for too many spots in strip
C
               N = N + 1
               IF (N.GT.MAX_SPOTS_IN_STRIP) THEN
                  IERR = 2
                  ERRSTR = 'Too many spots in a strip for overlaps -'
                  CALL STR_ADDS (ERRSTR,'program limit of',1)
                  CALL STR_ADDI (ERRSTR,MAX_SPOTS_IN_STRIP,1)
                  CALL STR_ADDS (ERRSTR,'exceeded',1)
                  RETURN
               END IF

               IWORK1(N) = II
            END IF
C
C====== Loop for more reflexions
C
            GO TO 30
         END IF
C
C====== End of loop over the strips
C
 70   CONTINUE
C
C====== Mark non-overlapped nodals
C
      IF (NUM_SPOTS.GT.0) THEN
         DO 80 ISPOT = 1 , NUM_SPOTS
            II = INDX(ISPOT)
            NOD_IDX = 0
            CALL LRL_GETHKL (NOD_HKL(II), NH, NK, NL)
            IF (IABS(NH).GT.NOD_IDX) NOD_IDX = ABS(NH)
            IF (IABS(NK).GT.NOD_IDX) NOD_IDX = ABS(NK)
            IF (IABS(NL).GT.NOD_IDX) NOD_IDX = ABS(NL)
            IF (NOD_IDX.LE.NODHKL_MAX) THEN
C
C====== Do not mark if overlapped
C
               IF (NBITST(FLAGS(II),0,1).EQ.0) 
     +            CALL STBITS (FLAGS(II),4,1,1)
            END IF
 80      CONTINUE
      END IF
      RETURN
      END
