CD-Title: Laue Reflection List (LRL) Routines
CD-end

CD-Intro:
C These are a set of routines for generating and accessing a Laue Reflection
C List (LRL) for the current pack and plate based on the current values in
C the Laue Data Module (LDM). The list is actually a list of spots rather
C than individual reflections and each spot has a slot in which an integrated
C intensity may be stored. Spot details include multiplicity information
C which enables generatation of the individual reflection data if required.
C Functions are available to sort the spot list in a number of ways.
CD-end

CD-Section_order:
C 1 3 2 4 7 9 8 5 6 10 11
CD-end

CD-Section: Generating the Laue Reflection List
C A high level routine is available which checks whether a list needs to
C be generated or whether the relevant LDM parameters are unchanged since the
C last call to the function; this also calls routines to calculate the
C orientation and 'A' matrices prior to the actual generation of the
C list. The lower level functions for these tasks may also be accessed
C directly if desired.
CD-end

C LAST UPDATE: 
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C**********************
C**     LRL_GEN      **
C**********************
C
CD-Routine: Generate Laue Reflection List (High Level) - LRL_GEN
C This routine checks whether a new reflection list needs to
C be generated or whether the relevant LDM parameters are unchanged since the
C last call to the function. If required it calls routines to calculate the
C orientation and 'A' matrices prior to the actual generation of the
C list. Normally distortion corrected coordinates are calculated for the
C spot positions but there is an option to omit the distortion correction.
C Various ways of determining the criterion for spatially overlapped spots
C are available.
CD-end
CD-Fortran:      
      SUBROUTINE LRL_GEN (IPACK, IPLATE, DISTOR, IOVTYP, IERR, ERRSTR)
CD-end
      INCLUDE 'implicit_none.finc'
C
C Purpose: Generate a predicted Laue Reflection List for a given
C          pack and plate
C
C Author:  John W. Campbell, July 1993
C
C Arguments:
C
       INTEGER IPACK
       INTEGER IPLATE
       LOGICAL DISTOR
       INTEGER IOVTYP
       INTEGER IERR
       CHARACTER*(*) ERRSTR
CD-Parameters:
C IPACK     i (R)   Pack number
C IPLATE    i (R)   Plate number
C DISTOR    l (R)   Flag = .true.  Include distortion corrections 
C                                  (normal case)
C                        = .false. Do not include distortion corrections
C IOVTYP    i (R)   Flag for calculation of spatial overlaps
C                        = 0, do not calculate these
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                    Note: If radial masking is not used (SPOT_WIDTH = 0.0)
C                          then if IOVTYP > 0, the overlaps are calculated
C                          on the basis of SPOT_DELTA
C
C IERR      i (W)    Returned error flag = 0, OK
C                                        =-1, Invalid pack (outside 
C                                             program limits)
C                                        =-2, Invalid plate (outside 
C                                             program limits)
C                                        = 1, Too many reflections to store    
C                                        = 2, Invalid overlap type flag
C                                        = 3, Too many spots in a strip for
C                                             overlaps calculation.
C
C ERRSTR    c (W)    Returned error message string (max 80 chars)
CD-end
C====== Parameter statements
C
C
C====== GLOBALS:
C
      INCLUDE 'lrl_data.finc'
C
C====== FUNCTION TYPES:
C
      EXTERNAL LRL_BDAT
C
C====== LOCALS:
C
      LOGICAL MTX
      LOGICAL GEN
      LOGICAL OVL
      LOGICAL DISTOR_CUR
      INTEGER KERR
      REAL ORIENT(3,3)
      REAL A(3,3)

C     LOGICAL MTX          !Matrices need to be recalculated flag
C     LOGICAL GEN          !LRL needs to be regenerated flag
C     LOGICAL OVL          !Overlaps need to be recalculated flag
C     LOGICAL DISTOR_CUR   !Distortion flag for current generation when
C                           routine was called
C     INTEGER KERR         !Error flag
C     REAL ORIENT(3,3)     !Orientation matrix
C     REAL A(3,3)          !'A' matrix
C
C-------------------------------------------------------------------------------
C
C====== Check parameters
C
      IERR = 0
      ERRSTR = ' '
      IF (IOVTYP.LT.0.OR.IOVTYP.GT.6) THEN
         IERR = 2
         ERRSTR = 'Invalid overlaps calculation flag'
         RETURN
      END IF
C
C====== Check pack/plate number
C
      CALL LDM_CHK_GPL (IPACK, IPLATE, IERR, ERRSTR)
      IF (IERR.NE.0) RETURN
C
C====== Check what needs to be done
C
      DISTOR_CUR = DISTOR_CORR
      CALL LDM_CHK_GEN (IPACK, IPLATE, MTX, GEN, OVL)
C
C====== Calculate orientation and 'A' matrices
C
      IF (MTX) THEN
         CALL LRL_MK_ORIENT (ORIENT)
         CALL LRL_A_CALC (IPACK, ORIENT, A)
         CALL LDM_LRL_CALC (IPACK, IPLATE, 1)
      END IF
C
C====== Generate reflection list
C
      IF (GEN.OR.(DISTOR.NEQV.DISTOR_CUR)) THEN
         CALL LRL_GNR (IPACK, IPLATE, A, DISTOR, KERR, ERRSTR)
         IF (KERR.NE.0) THEN
            IERR = 1
            RETURN
         END IF
         CALL LDM_LRL_CALC (IPACK, IPLATE, 2)
      END IF
C
C====== Find overlaps
C
      IF (OVL.OR.(DISTOR.NEQV.DISTOR_CUR)) THEN
         CALL LRL_OVLP (IOVTYP, KERR, ERRSTR)
         IF (KERR.NE.0) THEN
            IF (KERR.EQ.1) IERR = 2
            IF (KERR.EQ.2) IERR = 3
            RETURN
         END IF
         CALL LDM_LRL_CALC (IPACK, IPLATE, 3)
      END IF
      RETURN
      END
C LAST UPDATE: 
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C**************************
C**     LRL_GNR          **
C**************************
C
CD-Routine: Generate Laue Reflection List (Low Level) - LRL_GNR
C This routine performs the actual generation of a Laue Reflection List
C using the current LDM parameters and Orientation and 'A' matrices.
C Normally distortion corrected coordinates are calculated for the
C spot positions but there is an option to omit the distortion correction.
C This is used as a service routine to LRL_GEN but may also be called 
C directly if required.
CD-end
CD-Fortran:
      SUBROUTINE LRL_GNR   (IPACK, IPLATE, A, DISTOR, IERR, ERRSTR)
CD-end
      INCLUDE 'implicit_none.finc'
C
C Purpose: Generate reflections list for a given pack and plate 
C
C          (Used as a service routine for LRL_GEN but may be
C           called directly if desired)
C
C Author:  Based on code from GENLAUE by Clifton et al.
C
C Arguments:
C
      INTEGER IPACK
      INTEGER IPLATE
      REAL A(3,3)
      LOGICAL DISTOR
      INTEGER IERR
      CHARACTER*(*) ERRSTR

CD-Parameters:
C  IPACK    i (R)   Pack number (pre-check before calling this routine)
C  IPLATE   i (R)   Plate number (pre-check before calling this routine)
C  A(3,3)   r (R)   'A' matrix
C  DISTOR   l (R)   Flag = .true. apply distortion correction to stored
C                   coordinates (normal case), =.false. do not.
C  IERR     i (W)   Error flag, =0 OK, =1 too many reflections for arrays
C                   (list truncated)
C  ERRSTR   c (W)   Error message string (max 80 chars)
CD-end
C
C====== GLOBALS:
C
      INCLUDE 'lrl_data.finc'
C
C====== FUNCTION TYPES:
C
      EXTERNAL LRL_BDAT
      LOGICAL LDM_RTOD
      INTEGER LRL_HCF
C
C====== LOCALS:
C
      LOGICAL FOUND_ANY_L
      LOGICAL FOUND_ANY_K
      LOGICAL ABSENT
      LOGICAL EOR
      LOGICAL WITHIN
      INTEGER I, J
      INTEGER II
      INTEGER H, K, L
      INTEGER HN, KN, LN
      INTEGER NEXT_H, NEXT_K, NEXT_L
      INTEGER NTWO20
      INTEGER NTWO10
      INTEGER NTWO12
      INTEGER NTWO4
      INTEGER HCF1
      INTEGER HCF2
      INTEGER IHRM
      INTEGER MINHRM
      INTEGER MAXHRM
      INTEGER MULT
      INTEGER LASTNOD
      INTEGER KERR
      INTEGER IFLAG
      INTEGER ND
      INTEGER NXRAST
      INTEGER NYRAST
      REAL    WSTD
      REAL    WMIN
      REAL    WMAX
      REAL    DMIN
      REAL    FILMR
      REAL    FILM_RAD2
      REAL    X_CEN_F
      REAL    Y_CEN_F
      REAL    X_C
      REAL    Y_C
      REAL    W_C
      REAL    SPINDLE
      REAL    CTOF
      REAL    XLOW
      REAL    XHIGH
      REAL    YLOW
      REAL    YHIGH
      REAL    RREJ
      REAL    RREJSQ
      REAL    XMIN
      REAL    XMAX
      REAL    YMIN
      REAL    YMAX
      REAL    XFD
      REAL    YFD
      REAL    RSQ
      REAL    DSTAR2
      REAL    DSTAR_MAX2
      REAL    DSTAR2_K_MAX
      REAL    DSTAR2_L_MAX
      REAL    SPIN(3,3)
      REAL    X, Y, Z
      REAL    X1, Y1, Z1
      REAL    XF, YF
      REAL    RMIN, RMAX
      REAL    RMIN2, RMAX2
      REAL    RF
      REAL    DEG_TO_RAD
      REAL    ALM
      REAL    DM2

C     LOGICAL FOUND_ANY_L             !Reflections found along l row
C     LOGICAL FOUND_ANY_K             !Reflections found along k row
C     LOGICAL ABSENT                  !Systematic absence flag
C     LOGICAL EOR                     !End of reflections flag
C     LOGICAL WITHIN                  !Reflection within bounds flag
C     INTEGER I, J                    !Array index variables
C     INTEGER II                      !Do loop variable 
C     INTEGER H, K, L                 !Reflection indices h, k, l
C     INTEGER HN, KN, LN              !Nodal indices h, k, l
C     INTEGER NEXT_H, NEXT_K, NEXT_L  !Increment in h or k or l for search
C     INTEGER NTWO20                  !2**20
C     INTEGER NTWO10                  !2**10
C     INTEGER NTWO12                  !2**12
C     INTEGER NTWO4                   !2**4
C     INTEGER HCF1                    !Highest common factor (1)
C     INTEGER HCF2                    !Highest common factor (2)
C     INTEGER IHRM                    !Harmonic number
C     INTEGER MINHRM                  !Minimum harmonic number
C     INTEGER MAXHRM                  !Maximum harmonic number
C     INTEGER MULT                    !Multiplicity (count)
C     INTEGER KERR                    !Error flag
C     INTEGER IFLAG                   !Flag from LDM_GET.. routines
C     INTEGER ND                      !No. decimal places from LDM_GET.. 
C     INTEGER NXRAST                  !No. of x-rasters
C     INTEGER NYRAST                  !No. of y-rasters
C     REAL    WSTD                    !Standard wavelength
C     REAL    WMIN                    !Minimum wavelength
C     REAL    WMAX                    !Maximum wavelength
C     REAL    DMIN                    !'dmin'
C     REAL    FILMR                   !Max radius to use (mm)
C     REAL    FILM_RAD2               !FILMR**2
C     REAL    X_CEN_F                 !Uncorrected x-centre in rasters
C     REAL    Y_CEN_F                 !Uncorrected y-centre in rasters
C     REAL    X_C                     !x_c in mm
C     REAL    Y_C                     !y_c in mm
C     REAL    W_C                     !w_c in degrees
C     REAL    SPINDLE                 !Spindle angle
C     INTEGER LASTNOD                 !Saves value of last packed nodal indices 
C     REAL    CTOF                    !Crystal to film distance
C     REAL    XLOW                    !Low x limit in rasters
C     REAL    XHIGH                   !High x limit in rasters
C     REAL    YLOW                    !Low y limit in raster
C     REAL    YHIGH                   !High y limit in rasters
C     REAL    RREJ                    !Inner rejection radius (mm)
C     REAL    RREJSQ                  !The above squared
C     REAL    XMIN                    !Min allowed x (mm from centre)
C     REAL    XMAX                    !Max allowed x (mm from centre)
C     REAL    YMIN                    !Min allowed y (mm from centre)
C     REAL    YMAX                    !Max allowed y (mm from centre)
C     REAL    XFD                     !Predicted x coodinate (mm from centre)
C                                      (may be distortion corrected)
C     REAL    YFD                     !Predicted y coordinate (mm from centre)
C                                      (may be distortion corrected)
C     REAL    RSQ                     !Distance squared of predicted spot
C                                      from centre (in mm**2)
C     REAL    DSTAR2                  !d* squared
C     REAL    DSTAR_MAX2              !Maximum d* squared
C     REAL    DSTAR2_K_MAX            !Maximum d* squared along k
C     REAL    DSTAR2_L_MAX            !Maximum d* squared along l
C     REAL    SPIN(3,3)               !Spindle rotation matrix
C     REAL    X, Y, Z                 !Reciprocal lattice coordinates
C     REAL    X1, Y1, Z1              !Recip. lattice coords after spindle rotn.
C     REAL    XF, YF                  !Film coordinates
C     REAL    RMIN, RMAX              !Min, max reflecting sphere radii
C     REAL    RMIN2, RMAX2            !RMIN squared, RMAX squared
C     REAL    RF                      !Radius of reflecting sphere for lambda
C                                      = 1/|RF|
C     REAL    DEG_TO_RAD              !Conversion factor, degrees to radians
C     REAL    ALM                     !Saved lambda value
C     REAL    DM2                     !Saved dmin_threshold2 value
C
      DATA NTWO20/1048576/, NTWO10/1024/
      DATA NTWO12/4096/, NTWO4/16/
      DATA DEG_TO_RAD /0.0174533/
C
C-------------------------------------------------------------------------------
C
C====== LRL initialisations
C
      NUM_SPOTS = 0
      NSING = 0
      NMULT = 0
      NUM_REFLEX = 0
      NTRIES = 0
      NUM_SPAT = 0
      NUM_CLOS = 0
      LASTHKL(1) = 0
      LASTHKL(2) = 0
      LASTHKL(3) = 0
      XSTART = 1.0E09
      XEND = -1.0E09
      INTEGRATED = .FALSE.
      OVLAP_STATUS = 0
      ISORT = 0
      DISTOR_CORR = DISTOR
      INTCODE_TYPE = 0
      IPK = IPACK
      IPL = IPLATE
C
C====== Other initialisations
C
      IERR = 0
      WSTD = 1.0
      CALL LDM_GET_SPIN (IPACK, SPINDLE, IFLAG, ND, KERR, ERRSTR)
      CALL LDM_GET_LMIN (IPACK, WMIN, IFLAG, ND, KERR, ERRSTR)
      CALL LDM_GET_LMAX (IPACK, WMAX, IFLAG, ND, KERR, ERRSTR)
      CALL LDM_GET_DMIN (IPACK, DMIN, IFLAG, ND, KERR, ERRSTR)
      CALL LDM_GET_RMAX (FILMR, IFLAG, ND)
      CALL LDM_GET_CTOF (IPACK, IPLATE, CTOF, IFLAG, ND, KERR, ERRSTR)
      GCTOF = CTOF
      GLMIN = WMIN
      GLMAX = WMAX
      GDMIN = DMIN
      LMIN_NEW = 0.0
      LMAX_NEW = 0.0
      DMIN_NEW = 0.0
      DTHR2_NEW = 0.0
      CALL LDM_GET_NXRAS (NXRAST, IFLAG)
      CALL LDM_GET_NYRAS (NYRAST, IFLAG)
      CALL LDM_MMRAST (IPK, IPL, GMM_RAST_X, GMM_RAST_Y)
      GXCEN = 0.0
      GYCEN = 0.0
      CALL LDM_GET_XCENF (IPK, IPL, X_CEN_F, IFLAG, ND, KERR, ERRSTR)
      CALL LDM_GET_YCENF (IPK, IPL, Y_CEN_F, IFLAG, ND, KERR, ERRSTR)
      CALL LDM_GET_XC (IPK, IPL, X_C, IFLAG, ND, KERR, ERRSTR)
      CALL LDM_GET_YC (IPK, IPL, Y_C, IFLAG, ND, KERR, ERRSTR)
      CALL LDM_GET_WC (IPK, IPL, W_C, IFLAG, ND, KERR, ERRSTR)
      GCOS_WC = COSD(W_C)
      GSIN_WC = SIND(W_C)
      IF (X_CEN_F.GT.0.0.AND.Y_CEN_F.GT.0.0) THEN 
         GXCEN = X_CEN_F + X_C*GMM_RAST_X
         GYCEN = Y_CEN_F + Y_C*GMM_RAST_Y
      ELSE
         GXCEN = NXRAST/2
         GYCEN = NYRAST/2
      END IF
      CALL LDM_GET_XLOW (XLOW, IFLAG, ND)
      CALL LDM_GET_XHIGH (XHIGH, IFLAG, ND)
      CALL LDM_GET_YLOW (YLOW, IFLAG, ND)
      CALL LDM_GET_YHIGH (YHIGH, IFLAG, ND)
      CALL LDM_GET_RMIN (RREJ,IFLAG, ND)
C
C====== Calculate allowed x, y limits in mm from predicted centre based on
C       xlow, xhigh, ylow, yhigh
C
      XMIN = (1.0 - GXCEN)/GMM_RAST_X
      XMAX = (NXRAST - GXCEN)/GMM_RAST_X
      YMIN = (1.0 - GYCEN)/GMM_RAST_Y
      YMAX = (NYRAST - GYCEN)/GMM_RAST_Y
      IF (XLOW.GT.0.0) XMIN = (XLOW - GXCEN)/GMM_RAST_X
      IF (XHIGH.GT.0.0) XMAX = (XHIGH - GXCEN)/GMM_RAST_X
      IF (YLOW.GT.0.0) YMIN = (YLOW - GYCEN)/GMM_RAST_Y
      IF (YHIGH.GT.0.0) YMAX = (YHIGH - GYCEN)/GMM_RAST_Y
      RREJSQ = RREJ*RREJ
C
C====== Set up spindle rotation matrix
C
      SPIN(1,1) = COS(SPINDLE*DEG_TO_RAD)
      SPIN(1,2) = -SIN(SPINDLE*DEG_TO_RAD)
      SPIN(2,1) = -SPIN(1,2)
      SPIN(2,2) = SPIN(1,1)
C
C====== Generate reflections list
C
      RMIN = WSTD/WMAX
      RMAX = WSTD/WMIN
      RMIN2 = RMIN*RMIN
      RMAX2 = RMAX*RMAX

      FILM_RAD2 = FILMR*FILMR 
      DSTAR_MAX2 = (WSTD/DMIN)**2
C
C====== Start h loop
C
      H = 0
      NEXT_H = 1
C
C====== Start k loop
C
 20   K = 0
      NEXT_K = 1
      FOUND_ANY_K = .FALSE.
      DSTAR2_K_MAX = 1000.0
C
C====== Start l loop - skip 000
C
C
C====== Check h & k for systematic absences
C
 30   CALL LDM_SYSAHK(H,K,ABSENT)
      IF (ABSENT) THEN
         K = K + NEXT_K
         GO TO 30
      END IF
C
      L = 0
      NEXT_L = 1
      IF (H.EQ.0 .AND. K.EQ.0) L = 1
      FOUND_ANY_L = .FALSE.
      DSTAR2_L_MAX = 1000.0
C
C====== Next l
C
 40   NTRIES = NTRIES + 1
C
C====== Check h, k & l for systematic absences
C
      CALL LDM_SYSABS(H,K,L,ABSENT)
      IF (ABSENT) THEN
         L = L + NEXT_L
         GO TO 40
      END IF
C
C====== Check if within resolution
C====== Calculate reciprocal lattice co-ords from hkl values
C
      X1 = A(1,1)*H + A(1,2)*K + A(1,3)*L
      Y1 = A(2,1)*H + A(2,2)*K + A(2,3)*L
      Z1 = A(3,1)*H + A(3,2)*K + A(3,3)*L

      DSTAR2 = X1*X1 + Y1*Y1 + Z1*Z1

      IF (DSTAR2.LE.DSTAR_MAX2) THEN
C
C====== Apply spindle rotation
C
         X = SPIN(1,1)*X1 + SPIN(1,2)*Y1
         Y = SPIN(2,1)*X1 + SPIN(2,2)*Y1
         Z = Z1
C
C====== Convert x,y,z to xf,yf
C
C====== Is r.l. point between min and max spheres ?
C
         IF ((X+RMIN)**2+Y*Y+Z*Z.GE.RMIN2) THEN
            IF ((X+RMAX)**2+Y*Y+Z*Z.LE.RMAX2) THEN
               RF = DSTAR2/(2.0*X)
C
C====== Get detector coordinates  (xf yf)
C
               IF (LDM_RTOD (X, Y, Z, DSTAR2, CTOF, XF, YF)) THEN
C
C====== See if spot within allowed limits (Check before & after distortion
C       correction (if done). Done like this to avoid problems esp. with
C       polynomial spd correction where large coordinates outside the range
C       used to fit the polynomials may get large and bad corrections and 
C       spuriously appear within the image area
C
                  WITHIN = .TRUE.
                  RSQ = XF*XF + YF*YF
                  IF (RSQ.GT.FILM_RAD2) WITHIN = .FALSE.
                  IF (WITHIN.AND.RREJSQ.GT.0.0) THEN
                     IF (RSQ.LT.RREJSQ) WITHIN = .FALSE.
                  END IF
                  IF (WITHIN) THEN
                     IF (XF.LT.XMIN.OR.XF.GT.XMAX) WITHIN = .FALSE.
                     IF (YF.LT.YMIN.OR.YF.GT.YMAX) WITHIN = .FALSE.
                  END IF
                  IF (WITHIN) THEN
                     IF (DISTOR) THEN
                        CALL LDM_CORR (IPK, IPL, GCOS_WC, GSIN_WC,
     +                                 GCTOF, XF, YF, XFD, YFD)
                        RSQ = XFD*XFD + YFD*YFD
                        IF (RSQ.GT.FILM_RAD2) WITHIN = .FALSE.
                        IF (WITHIN.AND.RREJSQ.GT.0.0) THEN
                           IF (RSQ.LT.RREJSQ) WITHIN = .FALSE.
                        END IF
                        IF (WITHIN) THEN
                           IF (XFD.LT.XMIN.OR.XFD.GT.XMAX) 
     +                        WITHIN = .FALSE.
                           IF (YFD.LT.YMIN.OR.YFD.GT.YMAX) 
     +                     WITHIN = .FALSE.
                        END IF
                     ELSE
                        XFD = XF
                        YFD = YF
                     END IF
                  END IF
                  IF (WITHIN) THEN
C
C====== O.K. so present on detector
C
                     NUM_REFLEX = NUM_REFLEX + 1
                     IF (NUM_REFLEX.GT.MAX_REF) THEN
                        NUM_REFLEX = MAX_REF
                        IERR = 1
                        ERRSTR = 'Too many reflections generated -'
                        CALL STR_ADDS (ERRSTR,'program limit of',1)
                        CALL STR_ADDI (ERRSTR,MAX_REF,1)
                        CALL STR_ADDS (ERRSTR,'exceeded',1)
                        GO TO 50
                     END IF

                     HKL(NUM_REFLEX) = NTWO20*(H+512)
     +                               + NTWO10*(K+512)
     +                               + L + 512
                     XF_D(NUM_REFLEX) = XFD
                     YF_D(NUM_REFLEX) = YFD
                     LAMBDAS(NUM_REFLEX) = -WSTD/RF
                     DMIN_THR2(NUM_REFLEX) = DSTAR2
                     HCF1 = LRL_HCF(H,K)
                     HCF2 = LRL_HCF(HCF1,L)
                     HN = ISIGN(H/HCF2,H)
                     KN = ISIGN(K/HCF2,K)
                     LN = ISIGN(L/HCF2,L)
                     NOD_HKL(NUM_REFLEX) = NTWO20*(HN+512)
     +                                   + NTWO10*(KN+512)
     +                                   + LN + 512
                     MULTS(NUM_REFLEX) = 0
                     AII(NUM_REFLEX) = 0.0
                     SIGII(NUM_REFLEX) = 0.0
                     FLAGS(NUM_REFLEX) = 0
                     INDX(NUM_REFLEX) = 0
                     IF (XFD.LT.XSTART) XSTART = XFD
                     IF (XFD.GT.XEND) XEND = XFD
                  END IF
               END IF
            END IF
         END IF
C
         FOUND_ANY_L = .TRUE.
         FOUND_ANY_K = .TRUE.
C
C====== Increment l and continue
C
         L = L + NEXT_L
         GO TO 40
C
      END IF
C
C====== hkl outside resolution ... test if any hits for this l so far
C
      IF (FOUND_ANY_L) THEN
C
C====== yes ... reverse l direction
C
         IF (NEXT_L.EQ.1) THEN
            NEXT_L = -1
            L = -1
            GO TO 40
C
C====== or increment k if already reversed l
C
         ELSE
            K = K + NEXT_K
            GO TO 30
         END IF
C
      END IF
C
C====== no  ... check if getting closer ... save and continue if so
C
      IF (DSTAR2.LT.DSTAR2_L_MAX) THEN
         DSTAR2_L_MAX = DSTAR2
         L = L + NEXT_L
         GO TO 40
      END IF
C
C====== getting further away ... reverse l and continue
C
      IF (NEXT_L.EQ.1) THEN
         NEXT_L = -1
         L = -1
         GO TO 40
      END IF
C
C====== no hits for this l so far ... test if any for this k so far
C
      IF (FOUND_ANY_K) THEN
C
C====== yes ... reverse k direction
C
         IF (NEXT_K.EQ.1) THEN
            NEXT_K = -1
            K = -1
            GO TO 30
C
C====== or bump h if already reversed k
C
         ELSE
            H = H + NEXT_H
            GO TO 20
         END IF
C
      END IF
C
C====== no  ... check if getting closer ... save and continue if so
C
      IF (DSTAR2.LT.DSTAR2_K_MAX) THEN
         DSTAR2_K_MAX = DSTAR2
         K = K + NEXT_K
         GO TO 30
      END IF
C
C====== getting further away ... reverse k and continue
C
      IF (NEXT_K.EQ.1) THEN
         NEXT_K = -1
         K = -1
         GO TO 30
      END IF
C
C====== check if we can go in other h direction
C
      IF (NEXT_H.EQ.1) THEN
         NEXT_H = -1
         H = -1
         GO TO 20
      END IF
C
C====== no hits at all for this h value so stop searching further
C
50    LASTHKL(1) = H
      LASTHKL(2) = K
      LASTHKL(3) = L
C
C====== End of generation - now find multiples
C
      IF (NUM_REFLEX.EQ.0) RETURN
      CALL INDXIV (NOD_HKL, NUM_REFLEX, INDX)
      ISORT = 2
      EOR = .FALSE.
      MULT = 1
      I = 1
      J = INDX(1)
      LASTNOD = NOD_HKL(I)
60    I = I + 1
      IF (I.GT.NUM_REFLEX) THEN
         EOR = .TRUE.
      ELSE
         J = INDX(I)
         IF (NOD_HKL(J).EQ.LASTNOD) THEN
            MULT = MULT + 1
            GO TO 60
         END IF
      END IF
      IF (MULT.GT.1) THEN
         MINHRM = 100000
         MAXHRM = 0
         J = INDX(I-MULT)
         CALL LRL_GETHKL (NOD_HKL(J), HN, KN, LN)
         DO 70 II = I-MULT, I-1
            J = INDX(II)
            CALL LRL_GETHKL (HKL(J), H, K, L)
            IHRM = 1
            IF (HN.NE.0) THEN
               IHRM = H/HN
            ELSE IF (KN.NE.0) THEN
               IHRM = K/KN
            ELSE IF (LN.NE.0) THEN
               IHRM = L/LN
            END IF
            IF (IHRM.LT.MINHRM) THEN
               MINHRM = IHRM
               ALM = LAMBDAS(J)
               DM2 = DMIN_THR2(J)
            END IF
            IF (IHRM.GT.MAXHRM) MAXHRM = IHRM
70       CONTINUE
         J = INDX(I-MULT)
         H = MINHRM*HN
         K = MINHRM*KN
         L = MINHRM*LN
         HKL(J) = NTWO20*(H+512) + NTWO10*(K+512) + L + 512
         LAMBDAS(J) = ALM
         DMIN_THR2(J) = DM2
         MULTS(J) = NTWO20*MULT + NTWO12*MINHRM + NTWO4*MAXHRM
     +              + (MAXHRM-MINHRM)/(MULT-1)
         NMULT = NMULT + 1
      ELSE
         J = INDX(I-1)
         MULTS(J) = 1
         NSING = NSING + 1
      END IF
      IF (.NOT.EOR) THEN
         MULT = 1
         J = INDX(I)
         LASTNOD = NOD_HKL(J)
         GO TO 60
      END IF
C
C====== Compress reflections list to spots list (by eliminating entries
C       with MULTS(I) = 0)
C
      NUM_SPOTS = 0
      DO 80 I = 1, NUM_REFLEX
         IF (MULTS(I).EQ.0) GO TO 80
         NUM_SPOTS = NUM_SPOTS + 1
         IF (NUM_SPOTS.LT.I) THEN
            HKL(NUM_SPOTS) = HKL(I)
            XF_D(NUM_SPOTS) = XF_D(I)
            YF_D(NUM_SPOTS) = YF_D(I)
            LAMBDAS(NUM_SPOTS) = LAMBDAS(I)
            DMIN_THR2(NUM_SPOTS) = DMIN_THR2(I)
            MULTS(NUM_SPOTS) = MULTS(I)
            NOD_HKL(NUM_SPOTS) = NOD_HKL(I)
C
C====== Do not need to copy FLAGS, AII, SIGII as all zero at this stage
C
         END IF
80    CONTINUE      
      RETURN
      END
C LAST UPDATE: 
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C**************************
C**     LRL_MK_ORIENT    **
C**************************
C
CD-Routine: Calculate Orientation Matrix - LDM_MK_ORIENT
C This routine calculates the orientation matrix based on the current
C LDM parameters. This is normally used as a service routine to the
C subroutine LRL_GEN but may be called directly if required.
CD-end
CD-Fortran:
      SUBROUTINE LRL_MK_ORIENT (ORIENT)
CD-end
      INCLUDE 'implicit_none.finc'
C
C Purpose: Make the orientation matrix
C
C          (Used as a service routine for LRL_GEN but may be
C           called directly if desired)
C
C Author:  (Derived from make_orient in GENLAUE)
C
C Arguments:
C
      REAL ORIENT(3,3)
CD-Parameters:
C  ORIENT(3,3) r (W)    Returns the orientation matrix
CD-end
C Notes: Generalised matrix builder
C
C    The reciprocal cell parameters are calculated from the real cell
C    parameters.
C
C    ORIENT returns the orientation matrix
C      ... row 1 is components of reciprocal cell along X,
C      ... row 2 along Y, row 3 along Z.
C    Convention for crystal system:
C
C         1 - triclinic    2 - monoclinic    3 - orthorhombic
C         4 - tetragonal   5 - hexagonal     6 - rhombohedral
C         7 - cubic
C
C    Conventions for alignment are:
C       1..astar,   2..bstar,   3..cstar
C      -1..-astar, -2..-bstar, -3..-cstar
C
C    Thus if the astar axis lies along the X-axis in the direction of +x
C    then IBEAM = 1, and similarly for IROT.
C
C
C====== FUNCTION TYPES:
C
      EXTERNAL LRL_BDAT
C
C====== LOCALS:
C
      CHARACTER*1 VALSTR
      INTEGER I, J, K
      INTEGER IFLAG, ND
      INTEGER SI, SJ, SK
      INTEGER ISYST, IBEAM, IROT
      REAL    CELLP(6)
      REAL    RCELL(6)
      REAL    RL(3)
      REAL    CC(3), SS(3)
      REAL    COS_PSI, COS_RHO, A, C, S, FUN
      REAL    DEG_TO_RAD
      REAL    ALPHA, BETA, GAMMA
      REAL    WSTD
C
C-------------------------------------------------------------------------------
C
C====== Initialise & get LDM parameters
C
      WSTD = 1.0
      CALL LDM_GET_A (CELLP(1), IFLAG, ND)
      CALL LDM_GET_B (CELLP(2), IFLAG, ND)
      CALL LDM_GET_C (CELLP(3), IFLAG, ND)
      CALL LDM_GET_ALPHA (CELLP(4), IFLAG, ND)
      CALL LDM_GET_BETA (CELLP(5), IFLAG, ND)
      CALL LDM_GET_GAMMA (CELLP(6), IFLAG, ND)
      CALL LDM_GET_SYST (VALSTR, ISYST, IFLAG)
      CALL LDM_GET_BEAM (VALSTR, IBEAM, IFLAG)
      CALL LDM_GET_ROT (VALSTR, IROT, IFLAG)
C
C====== Get reciprocal cell parameters
C
      CALL LRL_RECCEL (RCELL,CELLP,WSTD)
C
C====== Clear the orientation matrix
C
      DO 80 I = 1 , 3
         DO 70 J = 1 , 3
            ORIENT(I,J) = 0.
 70      CONTINUE
 80   CONTINUE
C
C====== i, j, k have special meanings in this code.
C====== i refers to the crystal rotation axis, Z
C
      I = IABS(IROT)
      SI = 1
      IF (IROT.LT.0) SI = -1
C
C====== j refers to IBEAM, the x-ray beam axis, X
C
      J = IABS(IBEAM)
      SJ = 1
      IF (IBEAM.LT.0) SJ = -1
C
C====== k is therefore the other axis, Y
C
      K = 6 - I - J
C
C====== The sign of k relies on i,j,k and X,Y,Z being right-handed
C
      IF (K.EQ.J+1 .OR. K.EQ.J-2) THEN
         SK = SI*SJ
      ELSE
         SK = -SI*SJ
      END IF
C
C====== Zero orientation matrix
C
      CALL MTX_ZERO(ORIENT,3,3)
C
C====== Convert cell angles to radians
C
      DEG_TO_RAD = 3.1415926536/180.0
      ALPHA = RCELL(4)*DEG_TO_RAD
      BETA = RCELL(5)*DEG_TO_RAD
      GAMMA = RCELL(6)*DEG_TO_RAD
C
C====== Set reciprocal cell lengths
C
      RL(1) = RCELL(1)
      RL(2) = RCELL(2)
      RL(3) = RCELL(3)
C
C====== Now go and do specific stuff for each crystal system
C
      IF (ISYST.EQ.2) THEN
C
C====== Monoclinic
C
         CC(1) = 0.0
         CC(2) = COS(BETA)
         CC(3) = 0.0
C
         SS(1) = 1.0
         SS(2) = SIN(BETA)
         SS(3) = 1.0
C
         ORIENT(1,J) = SJ*RL(J)*SS(K)
         ORIENT(1,K) = SJ*RL(K)*CC(I)
         ORIENT(2,K) = SK*RL(K)*SS(I)*SS(J)
         ORIENT(3,I) = SI*RL(I)
         ORIENT(3,J) = SI*RL(J)*CC(K)
         ORIENT(3,K) = SI*RL(K)*CC(J)

      ELSE IF (ISYST.EQ.3) THEN
C
C====== Orthorhombic
C
         ORIENT(1,J) = SJ*RL(J)
         ORIENT(2,K) = SK*RL(K)
         ORIENT(3,I) = SI*RL(I)

      ELSE IF (ISYST.EQ.4) THEN
C
C====== Tetragonal
C
         ORIENT(1,J) = SJ*RL(J)
         ORIENT(2,K) = SK*RL(K)
         ORIENT(3,I) = SI*RL(I)

      ELSE IF (ISYST.EQ.5) THEN
C
C====== Hexagonal
C
         CC(1) = 0.0
         CC(2) = 0.0
         CC(3) = 0.5
C
         SS(1) = 1.0
         SS(2) = 1.0
         SS(3) = SQRT(3.0)/2.0
C
         ORIENT(1,J) = SJ*RL(J)*SS(K)
         ORIENT(1,K) = SJ*RL(K)*CC(I)
         ORIENT(2,K) = SK*RL(K)*SS(I)*SS(J)
         ORIENT(3,I) = SI*RL(I)
         ORIENT(3,J) = SI*RL(J)*CC(K)
         ORIENT(3,K) = SI*RL(K)*CC(J)

      ELSE IF (ISYST.EQ.6) THEN
C
C====== Rhombohedral
C
         A = RCELL(1)
         C = COS(ALPHA)
         S = SIN(ALPHA)
         FUN = SQRT(1.-3.*C**2+2.*C**3)
C
         ORIENT(1,J) = SJ*A*S
         ORIENT(1,K) = SJ*A*C*(1.-C)/S
         ORIENT(2,K) = SK*A*FUN/S
         ORIENT(3,I) = SI*A
         ORIENT(3,J) = SI*A*C
         ORIENT(3,K) = SI*A*C

      ELSE IF (ISYST.EQ.7) THEN
C
C====== Cubic
C
         A = RCELL(1)
C
         ORIENT(1,J) = SJ*A
         ORIENT(2,K) = SK*A
         ORIENT(3,I) = SI*A

      ELSE
C
C====== Triclinic
C
         CC(1) = COS(ALPHA)
         CC(2) = COS(BETA)
         CC(3) = COS(GAMMA)
C
         SS(1) = SIN(ALPHA)
         SS(2) = SIN(BETA)
         SS(3) = SIN(GAMMA)
C
C====== See Buerger's Xray Crystallography for the next formulae, p348
C
         COS_PSI = (CC(I)-CC(K)*CC(J))/SS(K)
         COS_RHO = SQRT(1.0-CC(I)**2-CC(J)**2-CC(K)**2+2.0*CC(I)*CC(J)
     +             *CC(K))/SS(K)
C
         ORIENT(1,J) = SJ*RL(J)*SS(K)
         ORIENT(1,K) = SJ*RL(K)*COS_PSI
         ORIENT(2,K) = SK*RL(K)*COS_RHO
         ORIENT(3,I) = SI*RL(I)
         ORIENT(3,J) = SI*RL(J)*CC(K)
         ORIENT(3,K) = SI*RL(K)*CC(J)
      END IF
      RETURN
      END
C LAST UPDATE: 
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C**************************
C**     LRL_A_CALC       **
C**************************
C
CD-Routine: Calculate the Master 'A' Matrix - LRL_A_CALC
C This routine calculates the master matrix, a_mat, from the crystal system,
C crystal alignment and cell parameters. This is then premultiplied by the phi 
C matrices. It is used as a service routine for the subroutine LRL_GEN but may
C be called directly if required.
CD-end
CD-Fortran:
      SUBROUTINE  LRL_A_CALC (IPACK, ORIENT, A_MAT)
CD-end
      INCLUDE 'implicit_none.finc'
C
C
C Purpose: Calculates the master matrix, a_mat, from the crystal system,
C          crystal alignment and cell parameters.
C          This is then premultiplied by the phi matrices.
C
C          (Used as a service routine for LRL_GEN but may be
C           called directly if desired)
C
C Author:  Based on code from GENLAUE by Clifton et al.
C
C Arguments:
C
      INTEGER    IPACK
      REAL       ORIENT(3,3)
      REAL       A_MAT(3,3)
CD-Parameters:
C IPACK        i (R)   Pack number
C ORIENT(3,3)  r (R)   Orientation matrix from crystal setting & cell 
C                      parameters
C A_MAT(3,3)   r (W)   The master matrix
CD-end
C
C====== FUNCTION TYPES:
C
      EXTERNAL LRL_BDAT
C
C====== LOCALS:
C
      CHARACTER*1 ERRSTR
      INTEGER IFLAG
      INTEGER ND
      INTEGER IERR
      REAL    PHIS(3)
      REAL    PC(3,3)
C
C------------------------------------------------------------------------------
C
C====== Get Phis
C
      CALL LDM_GET_PHIX (IPACK, PHIS(1), IFLAG, ND, IERR, ERRSTR)
      CALL LDM_GET_PHIY (IPACK, PHIS(2), IFLAG, ND, IERR, ERRSTR)
      CALL LDM_GET_PHIZ (IPACK, PHIS(3), IFLAG, ND, IERR, ERRSTR)
C
C====== Define phi rotation matrix
C
      CALL MTX_ROTATE(PHIS,PC)
C
C====== Apply phi rotations
C
      CALL MTX_MULT(PC,ORIENT,A_MAT,3,3,3)
      RETURN
      END
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

C LAST UPDATE:
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C************************
C**      LRL_OVALL     **
C************************
C
CD-Routine: Find All Spatial Overlaps - LRL_OVALL
C This routine is used to find all the spots (up to a requested limit) which 
C overlap with a given spatially overlapped spot. A number of ways of 
C specifying the overlap criterion are available (as in LRL_OVLP). The
C Laue Reflection List (LRL) must be sorted on 'xfd' or 'yfd' before this
C routine is called. 
CD-end
CD-Fortran:
      SUBROUTINE  LRL_OVALL (IOVTYP, N, NSPOT, SPOT_L, SPOT_W, 
     +                       SPOTFAC, SPOT_B, SPOT_D, RMAX, CTOF, 
     +                       MM_RAST_X, MM_RAST_Y, MAX_OVLPS, 
     +                       NUM_OV, X_OV, Y_OV, IERR, ERRSTR)
CD-end
      INCLUDE 'implicit_none.finc'
C
C Purpose:  to find list of overlapping spots
C
C          (Used as a service routine for LFN_SAINT but may be
C           called directly if desired)
C
C Author:  Q. Hao. Based on LRL_OVLP from J.W.Campbell et al.
C
C Arguments:
C
      INTEGER IOVTYP
      INTEGER N
      INTEGER NSPOT
      REAL SPOT_L
      REAL SPOT_W
      REAL SPOTFAC
      REAL SPOT_B
      REAL SPOT_D
      REAL RMAX
      REAL CTOF
      REAL MM_RAST_X
      REAL MM_RAST_Y
      INTEGER MAX_OVLPS
      INTEGER NUM_OV
      REAL X_OV(MAX_OVLPS)
      REAL Y_OV(MAX_OVLPS)
      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 N         i (R)    Spot number in sorted LRL
C NSPOT     i (R)    Number of spots in the LRL
C SPOT_L    r (R)    Spot length (mm)
C SPOT_W    r (R)    Spot width (mm)
C SPOTFAC   r (R)    Spot expansion factor (See LDM parameter defns.)
C SPOT_B    r (R)    Spot border width (pixels)
C SPOT_D    r (R)    Spot delta (mm)
C RMAX      r (R)    Maximum radius (mm)
C CTOF      r (R)    Crystal to image diastance
C MM_RAST_X r (R)    Conversion factor mm to x-rasters
C MM_RAST_Y r (R)    Conversion factor mm to y-rasters
C MAX_OVLPS i (R)    Maximum number of overlaps to return
C NUM_OV    i (W)    Number of overlapping spots found
C X_OV()    r (W)    'x' coordinates of the NUM_OV overlapping spots (rasters)
C Y_OV()    r (W)    'y' coordinates of the NUM_OV overlapping spots (rasters)
C IERR      i (W)     Error flag = 0, OK
C                                = 1, Spot list not sorted on 'xfd' 
C                                     or 'yfd'
C                                = 2, Spot was not found to be overlapped
C                                  
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
      EXTERNAL LRL_BDAT
C
C====== LOCALS:
C 
      LOGICAL    RAD_MASK
      LOGICAL    OVLP
      INTEGER    I
      INTEGER    J
      INTEGER    JJ
      INTEGER    ISPOT1
      INTEGER    ISPOT2
      REAL       X1
      REAL       Y1
      REAL       X2
      REAL       Y2
      REAL       XD
      REAL       YD
      REAL       RASTOMM
      REAL       DELTAM
      REAL       DELTA2
      REAL       D2
      REAL       X_DIFF
      REAL       Y_DIFF
C
C------------------------------------------------------------------------------
C
C====== Get sort order in LRL
C
      IF (ISORT.NE.3.AND.ISORT.NE.4) THEN
         IERR = 1
         ERRSTR = 'Reflections not sorted on xfd or yfd'
         RETURN
      END IF
C
C====== Get delta value or ellipse max. dimension if radial masking
C
      RASTOMM = 1.0/MM_RAST_X
      RAD_MASK = .FALSE.
      IF (SPOT_W.GT.0.0) RAD_MASK = .TRUE.
      IF (RAD_MASK) THEN
         DELTAM = SPOT_L
         IF (SPOTFAC.GT.0.0) THEN
            DELTAM = DELTAM*(1.0+SPOTFAC*RMAX/CTOF)       
         END IF
         IF (SPOT_W.GT.DELTAM) DELTAM = SPOT_W
         DELTAM = DELTAM + (1.0+SPOT_B)*RASTOMM
      ELSE
         DELTAM = SPOT_D
      END IF
C
C====== Adjust delta to a minimum of 10 micron
C
      DELTAM = MAX(DELTAM,0.01)
      DELTA2 = DELTAM**2
C
C====== Other initialisations
C
      NUM_OV = 0
      ISPOT1 = INDX(N)
      IF (DISTOR_CORR) THEN
        X1 = XF_D(ISPOT1)
        Y1 = YF_D(ISPOT1)
      ELSE
        CALL LDM_CORR (IPK, IPL, GCOS_WC, GSIN_WC, GCTOF, 
     +                 XF_D(ISPOT1), YF_D(ISPOT1), X1, Y1) 
      END IF
C
C====== Treat depending on sort order
C
C
C====== Sorted on 'xfd' - search back then forward
C
      IF (ISORT.EQ.3) THEN
         DO 60 JJ = -1, 1, 2
            DO 50 J = 1, NSPOT
               I = N + JJ*J
               IF (I.LT.1.OR.I.GT.NSPOT) GOTO 60   
               ISPOT2 = INDX(I)
               IF (DISTOR_CORR) THEN
                  X2 = XF_D(ISPOT2)
                  Y2 = YF_D(ISPOT2)
               ELSE
                  CALL LDM_CORR (IPK, IPL, GCOS_WC, GSIN_WC, GCTOF, 
     +                           XF_D(ISPOT2), YF_D(ISPOT2), X2, Y2) 
               END IF
C
C====== If X difference exceeds delta then finish this loop
C
               X_DIFF = ABS(X1-X2)
               IF (X_DIFF.GT.DELTAM) GO TO 60
C
C====== If X difference is possible then check Y difference.
C
               Y_DIFF = ABS(Y1-Y2)
               IF (Y_DIFF.LE.DELTAM) THEN
                  D2 = X_DIFF**2 + Y_DIFF**2
                  IF (D2.LE.DELTA2) THEN
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
                        NUM_OV = NUM_OV + 1
                        IF (NUM_OV.LE.MAX_OVLPS) THEN
                           XD = GXCEN + X2*GMM_RAST_X
                           YD = GYCEN + Y2*GMM_RAST_Y
                           X_OV(NUM_OV) = XD
                           Y_OV(NUM_OV) = YD
                        END IF
                     END IF
                  END IF
               END IF
 50         CONTINUE
 60      CONTINUE
C
C====== Sorted on 'yfd' - search back then forward
C
      ELSE IF (ISORT.EQ.4) THEN
         DO 80 JJ = -1, 1, 2
            DO 70 J = 1, NSPOT
               I = N + JJ*J
               IF (I.LT.1.OR.I.GT.NSPOT) GOTO 80   
               ISPOT2 = INDX(I)
               IF (DISTOR_CORR) THEN
                  X2 = XF_D(ISPOT2)
                  Y2 = YF_D(ISPOT2)
               ELSE
                  CALL LDM_CORR (IPK, IPL, GCOS_WC, GSIN_WC, GCTOF, 
     +                           XF_D(ISPOT2), YF_D(ISPOT2), X2, Y2) 
               END IF
C
C====== If Y difference exceeds delta then finish this loop
C
               Y_DIFF = ABS(Y1-Y2)
               IF (Y_DIFF.GT.DELTAM) GO TO 80
C
C====== If Y difference is possible then check X difference.
C
               X_DIFF = ABS(X1-X2)
               IF (X_DIFF.LE.DELTAM) THEN
                  D2 = X_DIFF**2 + Y_DIFF**2
                  IF (D2.LE.DELTA2) THEN
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
                        NUM_OV = NUM_OV + 1
                        IF (NUM_OV.LE.MAX_OVLPS) THEN
                           XD = GXCEN + X2*GMM_RAST_X
                           YD = GYCEN + Y2*GMM_RAST_Y
                           X_OV(NUM_OV) = XD
                           Y_OV(NUM_OV) = YD
                        END IF
                     END IF
                  END IF
               END IF
 70         CONTINUE
 80      CONTINUE
      END IF
C
C====== Check that overlap(s) found
C
      IF (NUM_OV.EQ.0) THEN
         IERR = 2
         ERRSTR = 'Inconsistent overlapping criteria'
         RETURN
      END IF
      IF (NUM_OV.GT.MAX_OVLPS) NUM_OV = MAX_OVLPS
      RETURN
      END

CD-Section: Calculate an Individual Spot Position
C A routine is available to calculate the predicted detector position
C for an individual Laue reflection. 
CD-end

C LAST UPDATE: 
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C**************************
C**     LRL_CALC_POS     **
C**************************
C
CD-Routine: Calculate Individual Reflection Position - LRL_CALC_POS
C This routine calculates the predicted detector position for an
C individual reflection. Distortion corrections are applied.
CD-end
CD-Fortran:
      SUBROUTINE LRL_CALC_POS (IPACK, IPLATE, X_CEN, Y_CEN, 
     +                         COSOM, SINOM,
     +                         AMAT, H, K, L, XFD, YFD, IERR)
CD-end
      INCLUDE 'implicit_none.finc'
C
C Purpose: Calculate the position of a Laue spot on the detector
C
C Author:  John W. Campbell (based on code from GENLAUE)
C
C Arguments:
C
      INTEGER IPACK
      INTEGER IPLATE
      REAL X_CEN
      REAL Y_CEN
      REAL COSOM
      REAL SINOM
      REAL AMAT(3,3)
      INTEGER H
      INTEGER K
      INTEGER L
      REAL XFD
      REAL YFD
      INTEGER IERR
CD-Parameters:
C IPACK     i (R)  Pack number
C IPLATE    i (R)  Plate number
C X_CEN     r (R)  Corrected x centre value (raster units)
C Y_CEN     r (R)  Corrected y centre value (raster units)
C COSOM     r (R)  cos (omega)
C SINOM     r (R)  sin (omega)
C AMAT(3,3) r (R)  A matrix calculated using current LDM parameters
C H         i (R)  Reflection h index
C K         i (R)  Reflection k index
C L         i (R)  Reflection l index
C XFD       r (W)  Predicted x position (raster units) (distortion 
C                  corrected)
C YFD       r (W)  Predicted y position (raster units) (distortion 
C                  corrected)
C IERR      i (W)  = 0 OK
C                  =-1, Invalid pack (outside program limits)
C                  =-2, Invalid plate (outside program limits)
C                  = 1 could not be recorded with current detector position
CD-end
C====== Parameter statements
C 

C
C====== GLOBALS:
C
C
C====== FUNCTION TYPES:
C
      LOGICAL LDM_RTOD
C
C====== LOCALS:
C
      INTEGER IFLAG, ND, KERR
      REAL CTOF
      REAL SPINDLE
      REAL SPIN(2,2)
      REAL X, Y, Z
      REAL XF, YF
      REAL X1, Y1, Z1
      REAL DSTAR2
      REAL DEG_TO_RAD
      REAL MM_RAST_X, MM_RAST_Y
      CHARACTER*2 ERRSTR

      DATA DEG_TO_RAD/0.0174533/
C
C-------------------------------------------------------------------------------
C
      XFD = 0.0
      YFD = 0.0
C
C====== Check pack/plate number
C
      CALL LDM_CHK_GPL (IPACK, IPLATE, IERR, ERRSTR)
      IF (IERR.NE.0) RETURN
C
C====== Get LDM parameters
C
      CALL LDM_GET_CTOF (IPACK, IPLATE, CTOF, IFLAG, ND, KERR, ERRSTR)
      CALL LDM_GET_SPIN (IPACK, SPINDLE, IFLAG, ND, KERR, ERRSTR)
      CALL LDM_MMRAST (IPACK, IPLATE, MM_RAST_X, MM_RAST_Y)

C
C====== Set up spindle rotation matrix
C
      SPIN(1,1) = COS(SPINDLE*DEG_TO_RAD)
      SPIN(1,2) = -SIN(SPINDLE*DEG_TO_RAD)
      SPIN(2,1) = -SPIN(1,2)
      SPIN(2,2) = SPIN(1,1)
      X1 = AMAT(1,1)*H + AMAT(1,2)*K + AMAT(1,3)*L
      Y1 = AMAT(2,1)*H + AMAT(2,2)*K + AMAT(2,3)*L
      Z1 = AMAT(3,1)*H + AMAT(3,2)*K + AMAT(3,3)*L

      DSTAR2 = X1*X1 + Y1*Y1 + Z1*Z1
C
C====== Apply spindle rotation
C
      X = SPIN(1,1)*X1 + SPIN(1,2)*Y1
      Y = SPIN(2,1)*X1 + SPIN(2,2)*Y1
      Z = Z1
 
      IF (LDM_RTOD(X, Y, Z, DSTAR2, CTOF, XF, YF)) THEN
         CALL LDM_CORR (IPACK, IPLATE, COSOM, SINOM, CTOF, XF, YF,
     +                  XFD, YFD)
         IERR = 0
         XFD = X_CEN + XFD*MM_RAST_X
         YFD = Y_CEN + YFD*MM_RAST_Y
      ELSE
         IERR = 1
      END IF
 
      RETURN
      END

CD-Section: Sort the Laue Reflection List
C A routine is available to sort the Laue Reflection List in a number of
C ways. In actual fact an internal index list is generated through which
C the spots in the 'sorted' list are accessed though the actual generated
C spot data remains in situ.
CD-end

C LAST UPDATE: 
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C**********************
C**     LRL_SORT     **
C**********************
C
CD-Routine: Sort Reflections List - LRL_SORT
C This routine sorts the Laue Reflection List in a number of
C ways. In actual fact an internal index list is generated through which
C the spots in the 'sorted' list are accessed; the actual generated
C spot data remains in situ. Spots may be sorted by index, detector
C position, lambda or resolution.
CD-end
CD-Fortran:
      SUBROUTINE LRL_SORT (ISOR)
CD-end
      INCLUDE 'implicit_none.finc'
C
C Purpose: Sort a Laue Reflection List in a number of different orders
C
C Author:  John W. Campbell, July 1993
C
C Arguments:
C
      INTEGER ISOR
CD-Parameters:
C ISOR    i (R)   Sort order flag
C                      = 1,  sort on k, k, l
C                      = 2,  sort on nodal h, k, l
C                      = 3,  sort on XFD (stored x-coordinate)
C                      = 4,  sort on YFD (stored y-coordinate)
C                      = 5,  sort on lambda
C                      = 6,  sort on dmin_thr2
CD-end
C====== GLOBALS:
C
      INCLUDE 'lrl_data.finc'
C
C====== FUNCTION TYPES:
C
      EXTERNAL LRL_BDAT
C
C====== LOCALS:
C

C
C-------------------------------------------------------------------------------
C
      IF (ISOR.LE.0.OR.ISOR.GT.6) THEN
         ISORT = 0
         RETURN
      END IF
      ISORT = ISOR
      IF (ISORT.EQ.1) THEN
         CALL INDXIV (HKL, NUM_SPOTS, INDX)
      ELSE IF (ISORT.EQ.2) THEN
         CALL INDXIV (NOD_HKL, NUM_SPOTS, INDX)
      ELSE IF (ISORT.EQ.3) THEN
         CALL INDXRV (XF_D, NUM_SPOTS, INDX)
      ELSE IF (ISORT.EQ.4) THEN
         CALL INDXRV (YF_D, NUM_SPOTS, INDX)
      ELSE IF (ISORT.EQ.5) THEN
         CALL INDXRV (LAMBDAS, NUM_SPOTS, INDX)
      ELSE IF (ISORT.EQ.6) THEN
         CALL INDXRV (DMIN_THR2, NUM_SPOTS, INDX)
      ELSE
         ISORT = 0
      END IF
      RETURN
      END

CD-Section: Get Global Laue Reflection List Data Items
C A series of routines are available to access data from the Laue
C Reflection List which refer to the list as a whole.
CD-end

C LAST UPDATE: 
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C**************************
C**     LRL_NUMSPOTS     **
C**************************
C
CD-Routine: Return the Number of Spots - LRL_NUMSPOTS
C This routine returns the number of spots in a Laue Reflection List (LRL)
CD-end
CD-Fortran:
      SUBROUTINE LRL_NUMSPOTS (NUMSPOTS)
CD-end
      INCLUDE 'implicit_none.finc'
C
C Purpose: Return number of spots in Laue Reflection List (LRL)
C
C Author:  John W. Campbell, July 1993
C
C Arguments:
C
      INTEGER NUMSPOTS
CD-Parameters:
C NUMSPOTS  i (W)  The number of spots in the Laue Reflection List
CD-end
C====== GLOBALS:
C
      INCLUDE 'lrl_data.finc'
C
C====== FUNCTION TYPES:
C
      EXTERNAL LRL_BDAT
C
C-------------------------------------------------------------------------------
C
      NUMSPOTS = NUM_SPOTS
      RETURN
      END
C LAST UPDATE: 
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C**************************
C**     LRL_COUNTS       **
C**************************
C
CD-Routine: Get Spot Counts - LRL_COUNTS
C This routine returns counts of the number of spots/reflections in the
C Laue Reflection List in a number of different categories.
CD-end
CD-Fortran:
      SUBROUTINE LRL_COUNTS (NUMSPOTS, NUMREFL, NUMSING, NUMMULT,
     +                       NUMSPAT, NUMCLOS)
CD-end
      INCLUDE 'implicit_none.finc'
C
C Purpose: Return counts from the Laue Reflection List (LRL)
C
C Author:  John W. Campbell, July 1993
C
C Arguments:
C
      INTEGER NUMSPOTS
      INTEGER NUMREFL
      INTEGER NUMSING
      INTEGER NUMMULT
      INTEGER NUMSPAT
      INTEGER NUMCLOS
CD-Parameters:
C NUMSPOTS  i (W)  The number of spots in the Laue Reflection List
C NUMREFL   i (W)  The number of reflections in the Laue Reflection List
C NUMSING   i (W)  The number of single spots in the Laue Reflection List
C NUMMULT   i (W)  The number of multiple spots in the Laue Reflection List    
C NUMSPAT   i (W)  The number of spots marked as spatially overlapped in
C                  in the Laue Reflection List
C NUMCLOS   i  (W) The number of spots marked as too close to spatially 
C                  deconvolute at integration time in the Laue Reflection 
C                  List
CD-end
C====== GLOBALS:
C
      INCLUDE 'lrl_data.finc'
C
C====== FUNCTION TYPES:
C
      EXTERNAL LRL_BDAT
C
C-------------------------------------------------------------------------------
C
      NUMSPOTS = NUM_SPOTS
      NUMREFL = NUM_REFLEX
      NUMSING = NSING
      NUMMULT = NMULT
      NUMSPAT = NUM_SPAT
      NUMCLOS = NUM_CLOS
      RETURN
      END
C LAST UPDATE: 
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C**************************
C**     LRL_STATUS       **
C**************************
C
CD-Routine: Return Status of Laue Reflection List - LRL_STATUS
C This routine returns status information from the Laue Reflection
C List including the number of spots, the current sort order, the
C spatial overlap criterion, whether or not the coordinates were
C distortion corrected and whether or not intensity integration has
C been carried out.
CD-end
CD-Fortran:
      SUBROUTINE LRL_STATUS (NUMSPOTS, ISORTED, IOVTYP, DISTOR, 
     +                       INTEG, INTCOD)
CD-end
      INCLUDE 'implicit_none.finc'
C
C Purpose: Return status flags from the Laue Reflection List (LRL)
C
C Author:  John W. Campbell, July 1993
C
C Arguments:
C
      INTEGER NUMSPOTS
      INTEGER ISORTED
      INTEGER IOVTYP
      LOGICAL DISTOR
      LOGICAL INTEG
      INTEGER INTCOD
CD-Parameters:
C NUMSPOTS   i (W)  The number of spots in the Laue Reflection List
C ISORTED    i (W)  Sort order 0 = none
C                              1 = sorted on h, k, l
C                              2 = sorted on nodal h, k, l
C                              3 = sorted on 'xfd'
C                              4 = sorted on 'yfd'
C                              5 = sorted on lambda
C                              6 = sorted on dminthr2
C IOVTYP    i (W)    Overlaps calculation option used (0=none calculated)
C                    (Value as used in call to LRL_GEN)
C DISTOR    l (W)    Distortion correction flag 
C                    = .true.  distortion correction was applied to stored
C                              coordinates
C                    = .false. distortion correction was not applied
C                    Note: LRL_GET routine returns distortion corrected 
C                          coordinates in either case; LRL_GET_XFYFD will
C                          return stored coordinates.
C INTEG     l (W)   Integration done flag =.true. yes, =.false. no
C INTCOD    i (W)   Integration code type flag = 0 if not known to the
C                   LRL routines.
C                   =1 Integration using Q. Hao's routine LFN_SAINT as
C                      used in INTLDM/LAUEGEN
CD-end
C====== GLOBALS:
C
      INCLUDE 'lrl_data.finc'
C
C====== FUNCTION TYPES:
C
      EXTERNAL LRL_BDAT
C
C-------------------------------------------------------------------------------
C
      NUMSPOTS = NUM_SPOTS
      ISORTED = ISORT
      IOVTYP = OVLAP_STATUS
      DISTOR = DISTOR_CORR
      INTEG = INTEGRATED
      INTCOD = INTCODE_TYPE
      RETURN
      END
C LAST UPDATE: 
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C**************************
C**     LRL_GET_XYCEN    **
C**************************
C
CD-Routine: Get Corrected Centre Position - LRL_GET_XYCEN
C This routine returns the corrected centre position used when the Laue
C Reflection List was generated. Note that it is only valid after such
C a list has been generated and does not provide a means of calculating
C the corrected centre position.
CD-end
CD-Fortran:
      SUBROUTINE LRL_GET_XYCEN (XCEN, YCEN)
CD-end
      INCLUDE 'implicit_none.finc'
C
C Purpose: Get corrected centre position (of predicted pattern) on image
C     
C          Note: Only valid after generating the reflection list
C
C Author:  John W. Campbell, July 1993
C
C Arguments:
C
      REAL XCEN
      REAL YCEN
CD-Parameters:    
C XCEN   r (W)  Corrected 'x' centre position in rasters (wrt image origin)    
C YCEN   r (W)  Corrected 'y' centre position in rasters (wrt image origin)
CD-end
C
C====== GLOBALS:
C
      INCLUDE 'lrl_data.finc'
C
C====== FUNCTION TYPES:
C
      EXTERNAL LRL_BDAT
C
C-------------------------------------------------------------------------------
C
      XCEN = GXCEN
      YCEN = GYCEN
      RETURN
      END
C LAST UPDATE: 
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C**************************
C**     LRL_TRIES        **
C**************************
C
CD-Routine: Number of Reflection Tries - LRL_TRIES
C this routine returns the number of reflection tried when generating the
C Laue Reflection List and the indices of the last reflection tried.
CD-end
CD-Fortran:
      SUBROUTINE LRL_TRIES (NUMTRIES, LAST_HKL)
CD-end
      INCLUDE 'implicit_none.finc'
C
C Purpose: Return number of reflections tried in generation and last h, k, l
C          tried
C
C Author:  John W. Campbell, July 1993
C
C Arguments:
C
      INTEGER NUMTRIES
      INTEGER LAST_HKL(3)
CD-Parameters:
C NUMTRIES    i (W)  The number reflections tried during the reflection 
C                    list generation.
C LAST_HKL(3) i (W)  The indices of the last reflection tried
CD-end
C
C====== GLOBALS:
C
      INCLUDE 'lrl_data.finc'
C
C====== FUNCTION TYPES:
C
      EXTERNAL LRL_BDAT
C
C-------------------------------------------------------------------------------
C
      NUMTRIES = NTRIES
      LAST_HKL(1) = LASTHKL(1)
      LAST_HKL(2) = LASTHKL(2)
      LAST_HKL(3) = LASTHKL(3)
      RETURN
      END

C LAST UPDATE: 
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C**************************
C**     LRL_GET_PLATE    **
C**************************
C
CD-Routine: Get Pack and Plate Number - LRL_GET_PLATE
C This routine returns the pack and plate number used when the Laue
C Reflection List was generated. Note that it is only valid after such
C a list has been generated.
CD-end
CD-Fortran:
      SUBROUTINE LRL_GET_PLATE (IPACK, IPLATE)
CD-end
      INCLUDE 'implicit_none.finc'
C
C Purpose: Get pack and plate no. used when LRL generated
C     
C          Note: Only valid after generating the reflection list
C
C Author:  John W. Campbell, September 1994
C
C Arguments:
C
      INTEGER IPACK
      INTEGER IPLATE
CD-Parameters:    
C IPACK  i (W)  Pack number
C IPLATE i (W)  Plate number
CD-end
C
C====== GLOBALS:
C
      INCLUDE 'lrl_data.finc'
C
C====== FUNCTION TYPES:
C
      EXTERNAL LRL_BDAT
C
C-------------------------------------------------------------------------------
C
      IPACK = IPK
      IPLATE = IPL
      RETURN
      END

C LAST UPDATE: 
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C**************************
C**     LRL_NODMAX       **
C**************************
C
CD-Routine: Return the Nodal Index Limit - LRL_NODMAX
C This routine returns the maximum nodal index for which a spot may be
C considered to be a nodal spot. This value is used in flagging nodal
C spots in the LRL.
CD-end
CD-Fortran:
      SUBROUTINE LRL_NODMAX (NODMAX)
CD-end
      INCLUDE 'implicit_none.finc'
C
C Purpose: Return maximum nodal index used
C
C Author:  John W. Campbell, September 1994
C
C Arguments:
C
      INTEGER NODMAX
CD-Parameters:
C NODMAX  i (W)  The maximum nodal index for which a spot may be
C                considered to be a nodal spot. This value is used 
C                in flagging nodal spots in the LRL.
CD-end
C====== GLOBALS:
C
      INCLUDE 'lrl_nodmax.finc'
C
C====== FUNCTION TYPES:
C
      EXTERNAL LRL_BDAT
C
C-------------------------------------------------------------------------------
C
      NODMAX = NODHKL_MAX
      RETURN
      END

CD-Section: Interface for XDL_VIEW Laue Simulations
C Special routines are available for interfacing to the XDL_VIEW X-windows
C toolkit Laue simulations view-object.
CD-end

C LAST UPDATE: 
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C**************************
C**     LRL_XDLSIM       **
C**************************
C
CD-Routine: Interface to X-windows Laue Simulation - LRL_XDLSIM
C This routine provides an interface the the XDL_VIEW Laue simulation
C view-object routine xdlf_laue_sim. Its use enables the internal arrays 
C storing the reflection data to be accessed directly and hence saves the user
C having to duplicate such information.
CD-end
CD-Fortran:
      SUBROUTINE LRL_XDLSIM (IVH, IVHPARENT, IX, IY, ICSET, ITYPE,
     +                       ISIZE, IREPL, MINW, MINH, IERR)
CD-end
      INCLUDE 'implicit_none.finc'
C
C Purpose: Interface to xdl_laue_sim routine
C
C          Note: This routine is supplied with the 'LRL' routines in
C                order to allow access to the internal 'LRL' data
C                structures which may not be accessed directly by
C                users' programs.
C
C Author:  John W. Campbell, July 1993
C
C Arguments:
C
      INTEGER IVH
      INTEGER IVHPARENT
      INTEGER IX
      INTEGER IY
      INTEGER ICSET
      INTEGER ITYPE
      INTEGER ISIZE
      INTEGER IREPL
      INTEGER MINW
      INTEGER MINH
      INTEGER IERR
CD-Parameters:   
C IVH        i (R)  View object handle for laue simulations view-object
C IVHPARENT  i (R)  View-object handle for parent, 0 if none
C IX         i (R)  X position (may be -1 if no parent to give default)
C IY         i (R)  Y position (may be -1 if no parent to give default)
C ICSET      i (R)  Colorset number
C ITYPE      i (R)  Type of simulation = 1, colour
C                                      = 2, interactive (b/w with sliders)
C ISIZE      i (R)  Halfwidth size in pixels
C IREPL      i (R)  =0 new simulation, =1 replace an existing one
C MINW       i (R)  Minimum width for view-object or 0 if no minimum
C MINH       i (R)  Minimum height for view-object or 0 if no minimum
C IERR       i (W)  Error flag =  0, OK
C                              = -1, no spots in list
C                              >  0, error code from xdl_laue_sim function
CD-end
C====== Parameter statements
C 
C
C====== GLOBALS:
C
      INCLUDE 'lrl_data.finc'
C
C====== LOCALS:
C
      CHARACTER*1 VALSTR
      CHARACTER*1 ERRSTR
      INTEGER I
      INTEGER LATT(4)
      INTEGER IFLAG
      INTEGER KERR
      INTEGER ND
      REAL  FILMR
      REAL  CELL(6)
      REAL  PHIS(4)
      REAL  CTOF
      REAL  WMIN
      REAL  WMAX
      REAL  DMIN
      REAL  D_LMIN
      REAL  D_LMAX
      REAL  D_DMIN
C
C-------------------------------------------------------------------------------
C
      IERR = 0
      IF (NUM_SPOTS.EQ.0) THEN
         IERR = -1
         RETURN
      END IF
      D_LMIN = 0.05
      D_LMAX = 0.05
      D_DMIN = 0.05
C
C====== Get sort index arrays
C
      IF (ISORT.NE.2) CALL LRL_SORT(2)
      CALL INDXRV (DMIN_THR2, NUM_SPOTS, IWORK1)
      DO 10 I = 1, NUM_SPOTS
         IWORK2(I) = IWORK1(NUM_SPOTS-I+1)
10    CONTINUE
      CALL INDXRV (LAMBDAS, NUM_SPOTS, IWORK1)
C
C====== Get LDM data values
C         
      CALL LDM_GET_RMAX (FILMR, IFLAG, ND)
      CALL LDM_GET_LMIN (IPK, WMIN, IFLAG, ND, KERR, ERRSTR)
      CALL LDM_GET_LMAX (IPK, WMAX, IFLAG, ND, KERR, ERRSTR)
      CALL LDM_GET_DMIN (IPK, DMIN, IFLAG, ND, KERR, ERRSTR)
      CALL LDM_GET_A (CELL(1), IFLAG, ND)
      CALL LDM_GET_B (CELL(2), IFLAG, ND)
      CALL LDM_GET_C (CELL(3), IFLAG, ND)
      CALL LDM_GET_ALPHA (CELL(4), IFLAG, ND)
      CALL LDM_GET_BETA (CELL(5), IFLAG, ND)
      CALL LDM_GET_GAMMA (CELL(6), IFLAG, ND)
      CALL LDM_GET_SYST (VALSTR, LATT(1), IFLAG)
      CALL LDM_GET_LATT (VALSTR, LATT(2), IFLAG)
      CALL LDM_GET_BEAM (VALSTR, LATT(3), IFLAG)
      CALL LDM_GET_ROT (VALSTR, LATT(4), IFLAG)
      CALL LDM_GET_PHIX (IPK, PHIS(1), IFLAG, ND, KERR, ERRSTR)
      CALL LDM_GET_PHIY (IPK, PHIS(2), IFLAG, ND, KERR, ERRSTR)
      CALL LDM_GET_PHIZ (IPK, PHIS(3), IFLAG, ND, KERR, ERRSTR)
      CALL LDM_GET_SPIN (IPK, PHIS(4), IFLAG, ND, KERR, ERRSTR)
      CALL LDM_GET_CTOF (IPK, IPL, CTOF, IFLAG, ND, KERR, ERRSTR)
C
C====== Call Laue simulation routine
C
      IF (IREPL.EQ.0) THEN
         CALL XDLF_LAUE_SIM (IVH, IVHPARENT, IX, IY, ICSET, ITYPE,
     +     ISIZE, NUM_SPOTS, HKL, XF_D, YF_D, LAMBDAS, DMIN_THR2, 
     +     MULTS, FLAGS, NOD_HKL, IWORK1, IWORK2, INDX, FILMR, 0.0,
     +     WMIN, WMAX, DMIN, D_LMIN, D_LMAX, D_DMIN, 
     +     CELL, PHIS, LATT, CTOF,
     +     MINW, MINH, IERR)
      ELSE
         CALL XDLF_LAUE_SIM_NEW (IVH, NUM_SPOTS, HKL, XF_D, YF_D, 
     +     LAMBDAS, DMIN_THR2, MULTS, FLAGS, NOD_HKL, 
     +     IWORK1, IWORK2, INDX, FILMR, 0.0, WMIN, WMAX, DMIN, 
     +     CELL, PHIS, LATT, CTOF, IERR)
      END IF
      RETURN
      END

C LAST UPDATE: 
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C**************************
C**     LRL_XDLGNM       **
C**************************
C
CD-Routine: Interface to X-windows Gnomonic Simulation - LRL_XDLGNM
C This routine provides an interface the the XDL_VIEW Laue simulation
C view-object routine xdlf_laue_sim in its mode to display a gnomonic
C projection of the Laue patters. Its use enables the internal arrays 
C storing the reflection data to be accessed directly and hence saves the user
C having to duplicate such information.
CD-end
CD-Fortran:
      SUBROUTINE LRL_XDLGNM (IVH, IVHPARENT, IX, IY, ICSET, ITYPE,
     +                       ISIZE, IREPL, GNOMR, MINW, MINH, IERR)
CD-end
      INCLUDE 'implicit_none.finc'
C
C Purpose: Interface to xdl_laue_sim routine for a gnomonic projection
C
C          Note: This routine is supplied with the 'LRL' routines in
C                order to allow access to the internal 'LRL' data
C                structures which may not be accessed directly by
C                users' programs.
C
C Author:  John W. Campbell, July 1993
C
C Arguments:
C
      INTEGER IVH
      INTEGER IVHPARENT
      INTEGER IX
      INTEGER IY
      INTEGER ICSET
      INTEGER ITYPE
      INTEGER ISIZE
      INTEGER IREPL
      REAL    GNOMR
      INTEGER MINW
      INTEGER MINH
      INTEGER IERR
CD-Parameters:   
C IVH        i (R)  View object handle for laue simulations view-object
C IVHPARENT  i (R)  View-object handle for parent, 0 if none
C IX         i (R)  X position (may be -1 if no parent to give default)
C IY         i (R)  Y position (may be -1 if no parent to give default)
C ICSET      i (R)  Colorset number
C ITYPE      i (R)  Type of simulation = 1, colour
C                                      = 2, interactive (b/w with sliders)
C ISIZE      i (R)  Halfwidth size in pixels
C IREPL      i (R)  =0 new simulation, =1 replace an existing one
C GNOMR      r (R)  Minimum radius in mm (>0.1) from image centre for 
C                   conversion of plot to a gnomonic projection.
C                   (If a value <= 0.1 is given a normal Laue 
C                   simulation will be given) 
C MINW       i (R)  Minimum width for view-object or 0 if no minimum
C MINH       i (R)  Minimum height for view-object or 0 if no minimum
C IERR       i (W)  Error flag =  0, OK
C                              = -1, no spots in list
C                              >  0, error code from xdl_laue_sim function
CD-end
C====== Parameter statements
C 
C
C====== GLOBALS:
C
      INCLUDE 'lrl_data.finc'
C
C====== LOCALS:
C
      CHARACTER*1 VALSTR
      CHARACTER*1 ERRSTR
      INTEGER I
      INTEGER LATT(4)
      INTEGER IFLAG
      INTEGER KERR
      INTEGER ND
      REAL  FILMR
      REAL  CELL(6)
      REAL  PHIS(4)
      REAL  CTOF
      REAL  WMIN
      REAL  WMAX
      REAL  DMIN
      REAL  D_LMIN
      REAL  D_LMAX
      REAL  D_DMIN
C
C-------------------------------------------------------------------------------
C
      IERR = 0
      IF (NUM_SPOTS.EQ.0) THEN
         IERR = -1
         RETURN
      END IF
      D_LMIN = 0.05
      D_LMAX = 0.05
      D_DMIN = 0.05
C
C====== Get sort index arrays
C
      IF (ISORT.NE.2) CALL LRL_SORT(2)
      CALL INDXRV (DMIN_THR2, NUM_SPOTS, IWORK1)
      DO 10 I = 1, NUM_SPOTS
         IWORK2(I) = IWORK1(NUM_SPOTS-I+1)
10    CONTINUE
      CALL INDXRV (LAMBDAS, NUM_SPOTS, IWORK1)
C
C====== Get LDM data values
C         
      CALL LDM_GET_RMAX (FILMR, IFLAG, ND)
      CALL LDM_GET_LMIN (IPK, WMIN, IFLAG, ND, KERR, ERRSTR)
      CALL LDM_GET_LMAX (IPK, WMAX, IFLAG, ND, KERR, ERRSTR)
      CALL LDM_GET_DMIN (IPK, DMIN, IFLAG, ND, KERR, ERRSTR)
      CALL LDM_GET_A (CELL(1), IFLAG, ND)
      CALL LDM_GET_B (CELL(2), IFLAG, ND)
      CALL LDM_GET_C (CELL(3), IFLAG, ND)
      CALL LDM_GET_ALPHA (CELL(4), IFLAG, ND)
      CALL LDM_GET_BETA (CELL(5), IFLAG, ND)
      CALL LDM_GET_GAMMA (CELL(6), IFLAG, ND)
      CALL LDM_GET_SYST (VALSTR, LATT(1), IFLAG)
      CALL LDM_GET_LATT (VALSTR, LATT(2), IFLAG)
      CALL LDM_GET_BEAM (VALSTR, LATT(3), IFLAG)
      CALL LDM_GET_ROT (VALSTR, LATT(4), IFLAG)
      CALL LDM_GET_PHIX (IPK, PHIS(1), IFLAG, ND, KERR, ERRSTR)
      CALL LDM_GET_PHIY (IPK, PHIS(2), IFLAG, ND, KERR, ERRSTR)
      CALL LDM_GET_PHIZ (IPK, PHIS(3), IFLAG, ND, KERR, ERRSTR)
      CALL LDM_GET_SPIN (IPK, PHIS(4), IFLAG, ND, KERR, ERRSTR)
      CALL LDM_GET_CTOF (IPK, IPL, CTOF, IFLAG, ND, KERR, ERRSTR)
C
C====== Call Laue simulation routine
C
      IF (IREPL.EQ.0) THEN
         CALL XDLF_LAUE_SIM (IVH, IVHPARENT, IX, IY, ICSET, ITYPE,
     +     ISIZE, NUM_SPOTS, HKL, XF_D, YF_D, LAMBDAS, DMIN_THR2, 
     +     MULTS, FLAGS, NOD_HKL, IWORK1, IWORK2, INDX, FILMR, GNOMR,
     +     WMIN, WMAX, DMIN, D_LMIN, D_LMAX, D_DMIN, 
     +     CELL, PHIS, LATT, CTOF,
     +     MINW, MINH, IERR)
      ELSE
         CALL XDLF_LAUE_SIM_NEW (IVH, NUM_SPOTS, HKL, XF_D, YF_D, 
     +     LAMBDAS, DMIN_THR2, MULTS, FLAGS, NOD_HKL, 
     +     IWORK1, IWORK2, INDX, FILMR, GNOMR, WMIN, WMAX, DMIN, 
     +     CELL, PHIS, LATT, CTOF, IERR)
      END IF
      RETURN
      END

CD-Section: Re-assign Spatial Overlaps and Multiplicities
C A set of routines have been written to enable an alternative assignment
C of spatially overlapped and multiple spots using revised values of the
C Laue soft limits. This effectively enables some re-processing of the
C data without having to re-integrate spot intensities provided that
C conservative values have been used for the soft limits when the reflection
C list was generated. The re-assigned flags are held separately from the
C original flags and are accessed via the routines described in this section.
CD-end

C LAST UPDATE:
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C****************************
C**      LRL_REASSIGN      **
C****************************
C
CD-Routine: Reassigning spatial overlaps - LRL_REASSIGN
C This routine enables the spatial overlap/multiplicity flags to be 
C re-assigned for different values of the Laue soft limits (lambda-min, 
C lambda-max, dmin) without regenerating the Laue Reflection List. The 
C re-assigned flags are saved separately from the original flags and
C may be accessed via the LRL_GET_REASS routine. The
C LDM values for the soft limits are not reset.
CD-end
CD-Fortran:
      SUBROUTINE  LRL_REASSIGN (IOVTYP, WMIN, WMAX, DMIN,
     +                          IERR, ERRSTR)
CD-end
      INCLUDE 'implicit_none.finc'
C
C Purpose: Sort reflexions in LRL to find and mark spatial overlaps with
C          re-assigned soft limits
C
C
C Author:  John W. Campbell, October 1993 (see also LRL_OVLP)
C
C Arguments:
C
      INTEGER IOVTYP
      REAL WMIN
      REAL WMAX
      REAL DMIN
      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 WMIN      r (R)   New value for lambda-min
C WMAX      r (R)   New value for lambda-max
C DMIN      r (R)   New value for 'dmin'
C
C                   Note: New lambda soft limits must be within range
C                         of those used in generating the reflection list 
C                         and new 'dmin' must not be lower than that used
C                         generating the reflection list; otherwise the
C                         limits will be reset to meet these conditions
C                         and IERR will be returned with a value of -1
C
C IERR      i (W)  Error flag = 0, OK
C                             =-1, New limits reset to be within range
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
      LOGICAL LRL_STILLOK
      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    N
      INTEGER    NUM_STRIPS
      INTEGER    STRIP
      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),2,2,0)
10    CONTINUE
      IF (IOVTYP.LT.0.OR.IOVTYP.GT.6) THEN
         IERR = 1
         ERRSTR = 'Invalid overlaps calculation flag'
         RETURN
      END IF
      LMIN_NEW = WMIN
      IF (LMIN_NEW.LT.GLMIN) THEN
         IERR = -1
         LMIN_NEW = GLMIN
      END IF
      IF (LMIN_NEW.GT.GLMAX) THEN
         IERR = -1
         LMIN_NEW = GLMAX
      END IF
      LMAX_NEW = WMAX
      IF (LMAX_NEW.GT.GLMAX) THEN
         IERR = -1
         LMAX_NEW = GLMAX
      END IF
      IF (LMAX_NEW.LT.GLMIN) THEN
         IERR = -1
         LMAX_NEW = GLMIN
      END IF
      DMIN_NEW = WMIN
      IF (DMIN_NEW.LT.GDMIN) THEN
         IERR = -1
         DMIN_NEW = GDMIN
      END IF
      DTHR2_NEW = 1.0/ ((DMIN_NEW)**2)

      OVLAP_STATUS = IOVTYP
C
C====== Get delta value or ellipse max. dimension if radial masking
C
      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

      CALL LDM_MMRAST (IPK, IPL, MM_RAST_X, MM_RAST_Y)
      RASTOMM = 1.0/MM_RAST_X

      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),3,1).EQ.0) THEN
                                 CALL STBITS (FLAGS(ISPOT1),2,2,3)
                                 NUM_CLOS = NUM_CLOS + 1
                              END IF

                              IF (NBITST(FLAGS(ISPOT2),3,1).EQ.0) THEN
                                 CALL STBITS (FLAGS(ISPOT2),2,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),2,1).EQ.0) 
     +                              THEN
                                    CALL STBITS (FLAGS(ISPOT1),2,1,1)
                                    NUM_SPAT = NUM_SPAT + 1
                                 END IF
   
                                 IF (NBITST(FLAGS(ISPOT2),2,1).EQ.0) 
     +                              THEN
                                    CALL STBITS (FLAGS(ISPOT2),2,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
 65         II = INDX(ISPOT)
            IF (.NOT.LRL_STILLOK(II)) THEN
               ISPOT = ISPOT + 1
               IF (ISPOT.LE.NUM_SPOTS) GO TO 65
               ISPOT = NUM_SPOTS
               GO TO 30
            END IF
            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
      RETURN
      END
C LAST UPDATE: 
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C**************************
C**     LRL_GET_REASS    **
C**************************
C
CD-Routine: Get Re-assigned Spot Flags - LRL_GET_REASS
C This routine enables the reassigned spot flags to be retrieved for
C a given spot in the Laue Reflections List after such re-assignment has
C been done via a call to the routine LRL_REASSIGN
CD-end
CD-Fortran:
      SUBROUTINE LRL_GET_REASS (ISPOT, MULT, MINHARM, MAXHARM, INCHARM,
     +                         SPOV, CLOS, OKNOD, IERR)
CD-end
      INCLUDE 'implicit_none.finc'
C
C Purpose: Get re-assigned data for a spot (after using LRL_REASSIGN)
C
C Author:  John W. Campbell, October 1993
C
C Arguments:
C
      INTEGER ISPOT
      INTEGER MULT
      INTEGER MINHARM
      INTEGER MAXHARM
      INTEGER INCHARM
      LOGICAL SPOV
      LOGICAL CLOS
      LOGICAL OKNOD
      INTEGER IERR
CD-Parameters:
C ISPOT   i (R)   Spot index number (in current sort order)
C MULT    i (W)   Multiplicity  (Note: may be 0 after re-assigning limits
C                 i.e. spot may no longer be recorded
C MINHARM i (W)   Minimum harmonic number
C MAXHARM i (W)   Maximum harmonic number
C INCHARM i (W)   Harmonics increment (returns 1 for a single)
C SPOV    l (W)   Spatially overlapped = .TRUE. yes, = .FALSE. no
C CLOS    l (W)   Too close to deconvolute = .TRUE. yes, = .FALSE. no
C OKNOD   l (W)   Non-spatially overlapped nodal = .TRUE. yes, = .FALSE. no
C IERR    i (W)   Error flag =  0,  OK
C                            =  1,  Invalid spot number
CD-end
C
C====== GLOBALS:
C
      INCLUDE 'lrl_data.finc'
C
C====== FUNCTION TYPES:
C
      INTEGER NBITST
      EXTERNAL LRL_BDAT
C
C====== LOCALS:
C   
      LOGICAL OK
      INTEGER I
      INTEGER IH
      INTEGER IK
      INTEGER IL
      INTEGER H_N
      INTEGER K_N
      INTEGER L_N
      INTEGER IHARM
      INTEGER MINH
      INTEGER MAXH
      REAL LAMBDA
      REAL DTHR2
C
C     LOGICAL OK      !Harmonic still OK flag
C     INTEGER I       !Index to spot
C     INTEGER IH      !'h' index
C     INTEGER IK      !'k' index
C     INTEGER IL      !'l' index
C     INTEGER H_N     !'h' nodal index
C     INTEGER K_N     !'k' nodal index
C     INTEGER L_N     !'l' nodal index
C     INTEGER IHARM   !Loop counter for harmonics
C     INTEGER MINH    !Minimum harmonic (before re-assign)
C     INTEGER MAXH    !Maximum harmonic (before re-assign)
C     REAL LAMBDA     !Lambda
C     REAL DTHR2      !1/dmin**2
C
C-------------------------------------------------------------------------------
C
      IERR = 0
      MINHARM = 0
      MAXHARM = 0
      INCHARM = 1
      SPOV = .FALSE.
      CLOS = .FALSE.
      OKNOD = .FALSE.
      IF (ISPOT.LT.1.OR.ISPOT.GT.NUM_SPOTS) THEN
         IERR = 1
         RETURN
      END IF
      I = INDX(ISPOT)
      CALL LRL_GETHKL (HKL(I), IH, IK, IL)
      CALL LRL_GETHKL (NOD_HKL(I), H_N, K_N, L_N)
      IF (MULTS(I).EQ.1) THEN
         MULT = 1
         IF (LAMBDAS(I).LT.LMIN_NEW.OR.LAMBDAS(I).GT.LMAX_NEW) THEN
            MULT = 0
            RETURN
         END IF
         IF (DMIN_THR2(I).GT.DTHR2_NEW) THEN
            MULT = 0
            RETURN
         END IF
         MINHARM = 0
         IF (H_N.NE.0) THEN
            MINHARM = IH/H_N
         ELSE IF (K_N.NE.0) THEN
            MINHARM = IK/K_N
         ELSE IF (L_N.NE.0) THEN
            MINHARM = IL/L_N
         END IF
         MAXHARM = MINHARM
         INCHARM = 1
      ELSE
         CALL LRL_GETMULT (MULTS(I), MULT, MINH, MAXH, INCHARM)
         MINHARM = MAXH
         MAXHARM = MINH
         MULT = 0
         DO 10 IHARM = MINH, MAXH, INCHARM
            OK = .TRUE.
            LAMBDA = LAMBDAS(I) * FLOAT(MINH) / FLOAT(IHARM)
            DTHR2 = DMIN_THR2(I) * FLOAT(IHARM) / FLOAT(MINH)
            IF (LAMBDA.LT.LMIN_NEW.OR.LAMBDA.GT.LMAX_NEW) OK = .FALSE.
            IF (DTHR2.GT.DTHR2_NEW) OK = .FALSE.
            IF (OK) THEN
               MULT = MULT + 1
               IF (IHARM.LT.MINHARM) MINHARM = IHARM
               IF (IHARM.GT.MAXHARM) MAXHARM = IHARM
            END IF
10       CONTINUE
         IF (MULT.EQ.0) THEN
            MINHARM = 0
            MAXHARM = 0
         END IF
      OKNOD = .FALSE.
      IF (NBITST(FLAGS(I),2,1).NE.0) SPOV = .TRUE.
      IF (NBITST(FLAGS(I),4,1).NE.0) OKNOD = .TRUE.
      IF (SPOV) OKNOD = .FALSE.
      IF (NBITST(FLAGS(I),3,1).NE.0) CLOS = .TRUE.
      END IF
      RETURN
      END
C LAST UPDATE: 
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C**************************
C**     LRL_REASS_OV     **
C**************************
C
C
      SUBROUTINE LRL_REASS_OV (ISPOT)
C
C
      INCLUDE 'implicit_none.finc'
C
C Purpose: Set reassigned spatially overlapped flag for a spot
C
C Author:  John W. Campbell, July 1993
C
C **** Service Routine Only ****
C
C Arguments:
C
      INTEGER ISPOT
C
C ISPOT   i  (R)   Spot index number (in current sort order)
C
C====== GLOBALS:
C
      INCLUDE 'lrl_data.finc'
C
C====== FUNCTION TYPES:
C
      EXTERNAL LRL_BDAT
C
C====== LOCALS:
C
      INTEGER I
C
C-------------------------------------------------------------------------------
C 
      IF (ISPOT.LT.1.OR.ISPOT.GT.NUM_SPOTS) RETURN
      I = INDX(ISPOT)
      CALL STBITS (FLAGS(I),3,1,1)
      RETURN
      END
C LAST UPDATE: 
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C**************************
C**     LRL_REASS_CL     **
C**************************
C
C
      SUBROUTINE LRL_REASS_CL (ISPOT)
C
C
      INCLUDE 'implicit_none.finc'
C
C Purpose: Set re-assigned 'too close' flag for a spot
C
C Author:  John W. Campbell, July 1993
C
C **** Service Routine Only ****
C
C Arguments:
C
      INTEGER ISPOT
C
C ISPOT   i (R)   Spot index number (in current sort order)
C
C====== GLOBALS:
C
      INCLUDE 'lrl_data.finc'
C
C====== FUNCTION TYPES:
C
      EXTERNAL LRL_BDAT
C
C====== LOCALS:
C
      INTEGER I
C
C-------------------------------------------------------------------------------
C 
      IF (ISPOT.LT.1.OR.ISPOT.GT.NUM_SPOTS) RETURN
      I = INDX(ISPOT)
      CALL STBITS (FLAGS(I),4,1,1)
      RETURN
      END

CD-Section: Get Individual Spot Data from Laue Reflection List
C This section describes a set of routines to get information about
C individual spots from the Laue Reflection List.
CD-end

C LAST UPDATE: 
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C**************************
C**     LRL_GET          **
C**************************
C
CD-Routine: Get Spot Details - LRL_GET
C This routine returns the basic details of a spot from the Laue Reflections 
C List. Additional routines are available to get other information as 
C described below.
CD-end
CD-Fortran:
      SUBROUTINE LRL_GET (ISPOT, IH, IK, IL, XD, YD, ALAM, DMTHR2,
     +                    MULT, OKNOD, SPOV, IERR)
CD-end
      INCLUDE 'implicit_none.finc'
C
C Purpose: Return spot details
C
C Author:  John W. Campbell, July 1993
C
C Arguments:
C
      INTEGER ISPOT
      INTEGER IH
      INTEGER IK
      INTEGER IL
      REAL XD
      REAL YD
      REAL ALAM
      REAL DMTHR2
      INTEGER MULT
      LOGICAL OKNOD
      LOGICAL SPOV
      INTEGER IERR
CD-Parameters:
C ISPOT  i (R)   Spot index number (in current sort order)
C IH     i (W)   Reflection 'h' index (lowest harmonic present for 
C                a multiple)
C IK     i (W)   Reflection 'k' index (lowest harmonic present 
C                for a multiple)
C IL     i (W)   Reflection 'l' index (lowest harmonic present for 
C                a multiple)
C XD     r (W)   'x' position on detector in raster units (corrected for
C                distortion and centre) wrt image origin
C YD     r (W)   'y' position on detector in raster units (corrected for
C                distortion and centre) wrt image origin
C ALAM   r (W)   Reflection lambda value (lowest harmonic for a multiple)
C DMTHR2 r (W)   Reflection dmin_thresold2 (lowest harmonic for a multiple)
C MULT   i (W)   Multiplicity
C OKNOD  l (W)   Flag = .true.  non-spatially overlapped nodal spot
C                     = .false. not
C SPOV   l (W)   Flag = .true. spatially overlapped spot, = .false. not
C IERR   i  (W)   Error flag =  0,  OK
C                            =  1,  Invalid spot number
C                            = -1,  Plate centre was not set (coordinates 
C                                   will be returned wrt. predicted pattern 
C                                   centre)
CD-end
C
C====== GLOBALS:
C
      INCLUDE 'lrl_data.finc'
C
C====== FUNCTION TYPES:
C
      INTEGER NBITST
      EXTERNAL LRL_BDAT
C
C====== LOCALS:
C   
      INTEGER J
      INTEGER NTWO20
      REAL XFD
      REAL YFD

      DATA NTWO20/1048576/
C
C-------------------------------------------------------------------------------
C
      IERR = 0
      IF (ISPOT.LT.1.OR.ISPOT.GT.NUM_SPOTS) THEN
         IERR = 1
         RETURN
      END IF
      J = INDX(ISPOT)
      CALL LRL_GETHKL (HKL(J), IH, IK, IL)
      IF (DISTOR_CORR) THEN
         XD = GXCEN + XF_D(J)*GMM_RAST_X
         YD = GYCEN + YF_D(J)*GMM_RAST_Y
      ELSE
         CALL LDM_CORR (IPK, IPL, GCOS_WC, GSIN_WC, GCTOF,
     +                  XF_D(J), YF_D(J), XFD, YFD)
         XD = GXCEN + XFD*GMM_RAST_X
         YD = GYCEN + YFD*GMM_RAST_Y
      END IF
      IF (GXCEN.LE.0.0.OR.GYCEN.LE.0.0) IERR = -1
      ALAM = LAMBDAS(J)
      DMTHR2 = DMIN_THR2(J)
      IF (MULTS(J).EQ.1) THEN
         MULT = 1
      ELSE
         MULT = MULTS(J)/NTWO20
      END IF
      OKNOD = .FALSE.
      IF (NBITST(FLAGS(J),4,1).NE.0) OKNOD = .TRUE.
      SPOV = .FALSE.
      IF (NBITST(FLAGS(J),0,1).NE.0) SPOV = .TRUE.
      RETURN
      END
C LAST UPDATE: 
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C**************************
C**     LRL_GET_HARM     **
C**************************
C
CD-Routine: Get Harmonics Data for a Spot - LRL_GET_HARM
C This routine returns the harmonics data for a spot in the Laue Reflections
C List. All the component reflections of a multiple may easily be generated
C from this information if required.
CD-end
CD-Fortran:
      SUBROUTINE LRL_GET_HARM (ISPOT, IH, IK, IL, H_N, K_N, L_N,
     +                         MULT, MINHARM, MAXHARM, INCHARM,
     +                         ALNOD, DTHR2NOD, IERR)
CD-end
      INCLUDE 'implicit_none.finc'
C
C Purpose: Get harmonics data for a spot
C
C Author:  John W. Campbell, July 1993
C
C Arguments:
C
      INTEGER ISPOT
      INTEGER IH
      INTEGER IK
      INTEGER IL
      INTEGER H_N
      INTEGER K_N
      INTEGER L_N
      INTEGER MULT
      INTEGER MINHARM
      INTEGER MAXHARM
      INTEGER INCHARM
      REAL    ALNOD
      REAL    DTHR2NOD
      INTEGER IERR
CD-Parameters:
C ISPOT   i (R)   Spot index number (in current sort order)
C IH      i (W)   Reflection 'h' index (lowest harmonic present 
C                 for a multiple)
C IK      i (W)   Reflection 'k' index (lowest harmonic present 
C                 for a multiple)
C IL      i (W)   Reflection 'l' index (lowest harmonic present 
C                 for a multiple)
C H_N     i (W)   Nodal 'h' index
C K_N     i (W)   Nodal 'k' index
C L_N     i (W)   Nodal 'l' index
C MULT    i (W)   Multiplicity
C MINHARM i (W)   Minimum harmonic number
C MAXHARM i (W)   Maximum harmonic number
C INCHARM i (W)   Harmonics increment (returns 1 for a single)
C IERR    i (W)   Error flag =  0,  OK
C                            =  1,  Invalid spot number
CD-end
C
C====== GLOBALS:
C
      INCLUDE 'lrl_data.finc'
C
C====== FUNCTION TYPES:
C
      EXTERNAL LRL_BDAT
C
C====== LOCALS:
C   
      INTEGER I
C
C     INTEGER I       !Index to spot
C
C-------------------------------------------------------------------------------
C
      IERR = 0
      IF (ISPOT.LT.1.OR.ISPOT.GT.NUM_SPOTS) THEN
         IERR = 1
         RETURN
      END IF
      I = INDX(ISPOT)
      CALL LRL_GETHKL (HKL(I), IH, IK, IL)
      CALL LRL_GETHKL (NOD_HKL(I), H_N, K_N, L_N)
      IF (MULTS(I).EQ.1) THEN
         MULT = 1
         MINHARM = 0
         IF (H_N.NE.0) THEN
            MINHARM = IH/H_N
         ELSE IF (K_N.NE.0) THEN
            MINHARM = IK/K_N
         ELSE IF (L_N.NE.0) THEN
            MINHARM = IL/L_N
         END IF
         MAXHARM = MINHARM
         INCHARM = 1
      ELSE
         CALL LRL_GETMULT (MULTS(I), MULT, MINHARM, MAXHARM, INCHARM)
      END IF
      ALNOD = LAMBDAS(I)*MINHARM
      DTHR2NOD = DMIN_THR2(I)/MINHARM
      RETURN
      END
C LAST UPDATE: 
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C**************************
C**     LRL_GET_XFYFD    **
C**************************
C
CD-Routine: Return Coordinates for a Spot - LRL_GET_XFYFD
C This routine returns the predicted coordinates for a spot in the Laue
C Reflections List. The spot position may or may not be distortion
C corrected depending on the option selected when the Laue Reflection List
C was generated.
CD-end
CD-Fortran:
      SUBROUTINE LRL_GET_XFYFD (ISPOT, XFD, YFD, IERR)
CD-end
      INCLUDE 'implicit_none.finc'
C
C Purpose: Get predicted coordinates for a spot 
C
C Author:  John W. Campbell, July 1993
C
C Arguments:
C
      INTEGER ISPOT
      REAL  XFD
      REAL  YFD
      INTEGER IERR
CD-Parameters:
C ISPOT   i (R)   Spot index number (in current sort order)
C XFD     r (W)   Predicted 'x' coordinate of stored spot along 'xd' ('xf')
C                 wrt. predicted pattern centre
C YFD     r (W)   Predicted 'y' coordinate of stored spot along 'yd' ('yf')
C                 wrt. predicted pattern centre
C                 Note: These will be distortion corrected if this option
C                       was selected when the reflection list was generated
C                       (normal) case (LRL_GEN, LRL_GNR). The distortion
C                       correction status can be found by calling the 
C                       routine LRL_STATUS
C IERR    i (W)   Error flag =  0,  OK
C                            =  1,  Invalid spot number
CD-end
C
C====== GLOBALS:
C
      INCLUDE 'lrl_data.finc'
C
C====== FUNCTION TYPES:
C
      EXTERNAL LRL_BDAT
C
C====== LOCALS:
C
      INTEGER I
C
C-------------------------------------------------------------------------------
C
      IERR = 0
      IF (ISPOT.LT.1.OR.ISPOT.GT.NUM_SPOTS) THEN
         IERR = 1
         RETURN
      END IF
      I = INDX(ISPOT)
      XFD = XF_D(I)
      YFD = YF_D(I)
      RETURN
      END
C LAST UPDATE: 
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C**************************
C**     LRL_GET_INT     **
C**************************
C
CD-Routine: Return Integrated Intensity of a Spot - LRL_GET_INT
C This routine will return the integrated intensity of a spot in the
C Laue Reflection List if this has been measured and stored (LRL_SET_INT).
CD-end
CD-Fortran:
      SUBROUTINE LRL_GET_INT (ISPOT, AI, SIGI, MEASURED, BAD,
     +                        OVERLOAD, ICODE, ICODE_TYPE, IERR)
CD-end
      INCLUDE 'implicit_none.finc'
C
C Purpose: Return integrated intensity data for a spot
C
C Author:  John W. Campbell, July 1993
C
C Arguments:
C
      INTEGER ISPOT
      REAL    AI
      REAL    SIGI
      LOGICAL MEASURED
      LOGICAL BAD
      LOGICAL OVERLOAD
      INTEGER ICODE
      INTEGER ICODE_TYPE
      INTEGER IERR
CD-Parameters:
C ISPOT      i  (R)   Spot index number (in current sort order)
C AI         r  (W)   The integrated intensity
C SIG I      r  (W)   sig(I)
C MEASURED   l  (W)   Spot has been integrated flag .TRUE. yes, .FALSE. no
C BAD        l  (W)   Bad spot flagged .TRUE. bad, .FALSE. good
C OVERLOAD   l  (W)   Overload flag .TRUE. yes, .FALSE. no
C ICODE      i  (W)   Status code fom integration program
C                     For ICODE_TYPE=1   bit 0 no non-overlapped peak pixels
C                                        bit 1 < 20% of peak pixels are
C                                              non-overlapped.
C                                        bit 2 <3 non-overlapped background
C                                              pixels
C                                        bit 4 < 20% of background pixels are
C                                              non-overlapped.
C ICODE_TYPE i  (W)   Integration code type; should be 0 if of a type 
C                     unknown to the 'LRL' routines
C                     = 1 Integration using Q. Hao's routines LFN_SAINT
C                         as used in INTLDM and LAUEGEN
C IERR       i  (W)   Error flag =  0,  OK
C                                =  1,  Invalid spot number
C                                =  2,  No integration done
CD-end
C
C
C====== GLOBALS:
C
      INCLUDE 'lrl_data.finc'
C
C====== FUNCTION TYPES:
C
      INTEGER NBITST
      EXTERNAL LRL_BDAT
C
C====== LOCALS:
C   
      INTEGER I
      INTEGER IFLAG
C
C     INTEGER I        !Index to spot
C     INTEGER IFLAG    !Spot flags word
C
C-------------------------------------------------------------------------------
C
      IERR = 0
      AI = 0.0
      SIGI = 0.0
      MEASURED = .FALSE.
      IF (ISPOT.LT.1.OR.ISPOT.GT.NUM_SPOTS) THEN
         IERR = 1
         RETURN
      END IF
      IF (.NOT.INTEGRATED) THEN
         IERR = 2
         RETURN
      END IF
      I = INDX(ISPOT)
      AI = AII(I)
      SIGI = SIGII(I)
      IFLAG = FLAGS(I)
      IF (NBITST(IFLAG,8,1).NE.0) MEASURED =.TRUE.
      BAD = .FALSE.
      IF (NBITST(IFLAG,9,1).NE.0) BAD = .TRUE.
      OVERLOAD = .FALSE.
      IF (NBITST(IFLAG,10,1).NE.0) OVERLOAD = .TRUE.
      ICODE = NBITST(IFLAG,12,8)
      ICODE_TYPE = INTCODE_TYPE
      RETURN
      END
C LAST UPDATE: 
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C**************************
C**     LRL_GET_NIDX     **
C**************************
C
CD-Routine: Get Nodal Index for a Spot - LRL_GET_NIDX
C This routine returns the 'nodal index' for a spot from the Laue Reflections
C List (i.e. the maximum absolute nodal index for the three indices).
CD-end
CD-Fortran:
      SUBROUTINE LRL_GET_NIDX (ISPOT, NOD_IDX, IERR)
CD-end
       INCLUDE 'implicit_none.finc'
C
C Purpose: Get 'nodal index' for a spot (i.e. max absolute nodal index)
C
C Author:  John W. Campbell, July 1993
C
C
      INTEGER ISPOT
      INTEGER NOD_IDX
      INTEGER IERR
CD-Parameters:
C ISPOT   i (R)   Spot index number (in current sort order)
C NOD_IDX i (W)   The nodal index
C IERR    i (W)   Error flag =  0,  OK
C                            =  1,  Invalid spot number
CD-end
C
C====== GLOBALS:
C
      INCLUDE 'lrl_data.finc'
C
C====== FUNCTION TYPES:
C
      EXTERNAL LRL_BDAT
C
C====== LOCALS:
C
      INTEGER I
      INTEGER NH, NK, NL
C
C-------------------------------------------------------------------------------
C
      NOD_IDX = 0
      IERR = 0
      IF (ISPOT.LT.1.OR.ISPOT.GT.NUM_SPOTS) THEN
         IERR = 1
         RETURN
      END IF
      I = INDX(ISPOT)
      CALL LRL_GETHKL (NOD_HKL(I), 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)
      RETURN
      END 
C LAST UPDATE: 
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C**************************
C**     LRL_GET_OVLP     **
C**************************
C
CD-Routine: Return Spatial Overlaps Flags for a Spot - LRL_GET_OVLP
C This routine returns the spatial overlaps flags for a spot in the
C Laue Reflections List.
CD-end
CD-Fortran:
      SUBROUTINE LRL_GET_OVLP (ISPOT, SPOV, CLOS, IERR)
CD-end
       INCLUDE 'implicit_none.finc'
C
C Purpose: Get spatial overlap flags for a spot 
C
C Author:  John W. Campbell, July 1993
C
C
      INTEGER ISPOT
      LOGICAL SPOV
      LOGICAL CLOS
      INTEGER IERR
CD-Parameters:
C ISPOT   i (R)   Spot index number (in current sort order)
C SPOV    l (W)   Spatially overlapped flag .TRUE. = yes, .FALSE. = no
C CLOS    l (W)   Too close for spatial deconvolution flag = .TRUE. yes,
C                 =. FALSE. no
C IERR    i (W)   Error flag =  0,  OK
C                            =  1,  Invalid spot number
CD-end
C====== GLOBALS:
C
      INCLUDE 'lrl_data.finc'
C
C====== FUNCTION TYPES:
C
      INTEGER NBITST
      EXTERNAL LRL_BDAT
C
C====== LOCALS:
C
      INTEGER I
C
C-------------------------------------------------------------------------------
C
      SPOV = .FALSE.
      CLOS = .FALSE.
      IERR = 0
      IF (ISPOT.LT.1.OR.ISPOT.GT.NUM_SPOTS) THEN
         IERR = 1
         RETURN
      END IF
      I = INDX(ISPOT)
      IF (NBITST(FLAGS(I),0,1).NE.0) SPOV = .TRUE.
      IF (NBITST(FLAGS(I),1,1).NE.0) CLOS = .TRUE.
      RETURN
      END
 
C LAST UPDATE: 
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C**************************
C**     LRL_GET_OVPIX    **
C**************************
C
CD-Routine: Return No. of Overload Pixels for a Spot - LRL_GET_OVPIX
C This routine returns the number of overloaded pixels in a spot if
C set by the integration routine for a spot in the Laue Reflections List.
CD-end
CD-Fortran:
      SUBROUTINE LRL_GET_OVPIX (ISPOT, NOVLD, IERR)
CD-end
       INCLUDE 'implicit_none.finc'
C
C Purpose: Returns the  no. of overloaded pixels if set 
C
C Author:  John W. Campbell, June 1997
C
C
      INTEGER ISPOT
      INTEGER NOVLD
      INTEGER IERR
CD-Parameters:
C ISPOT   i (R)   Spot index number (in current sort order)
C NOVLD   i (W)   Returns the number of overloaded pixels in the
C                 spot if set by the integration routine (0-255). 
C                 -1 indicates not set; 255 means 255 or more
C                 are overloaded. 
C IERR    i (W)   Error flag =  0,  OK
C                            =  1,  Invalid spot number
CD-end
C====== GLOBALS:
C
      INCLUDE 'lrl_data.finc'
C
C====== FUNCTION TYPES:
C
      INTEGER NBITST
      EXTERNAL LRL_BDAT
C
C====== LOCALS:
C
      INTEGER I
C
C-------------------------------------------------------------------------------
C
      NOVLD = -1
      IERR = 0
      IF (ISPOT.LT.1.OR.ISPOT.GT.NUM_SPOTS) THEN
         IERR = 1
         RETURN
      END IF
      I = INDX(ISPOT)
      IF (NBITST(FLAGS(I),20,1).EQ.0) RETURN
      NOVLD = (NBITST(FLAGS(I),21,8))
      RETURN
      END 

CD-Section: Set Individual Spot Data in Laue Reflection List
C The routines described in this section enable the setting of values
C associated with a particular spot in the Laue Reflection List.
CD-end

C LAST UPDATE: 
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C**************************
C**     LRL_SET_INT      **
C**************************
C
CD-Routine: Store Integrated Intensity - LRL_SET_INT
C This routine is used to store an integrated intensity value for a spot
C in the Laue Reflections List.
CD-end
CD-Fortran:
      SUBROUTINE LRL_SET_INT (ISPOT, AI, SIGI, BAD, OVERLOAD, 
     +                        ICODE, IERR)
CD-end
      INCLUDE 'implicit_none.finc'
C
C Purpose: Set integrated intensity data for a spot
C
C Author:  John W. Campbell, July 1993
C
C Arguments:
C
      INTEGER ISPOT
      REAL    AI
      REAL    SIGI
      LOGICAL BAD
      LOGICAL OVERLOAD
      INTEGER ICODE
      INTEGER IERR
CD-Parameters:
C ISPOT      i  (R)   Spot index number (in current sort order)
C AI         r  (R)   The integrated intensity
C SIG I      r  (R)   sig(I)
C BAD        l  (R)   Bad spot flag .TRUE. bad, .FALSE. good
C OVERLOAD   l  (R)   Overload flag .TRUE. yes, .FALSE. no
C ICODE      i  (R)   Status code from integration program (0 to 255)
C                     See for example LFN_SAINT or LRL_GET_INT
C IERR       i  (W)   Error flag =  0,  OK
C                                =  1,  Invalid spot number
CD-end
C
C====== GLOBALS:
C
      INCLUDE 'lrl_data.finc'
C
C====== FUNCTION TYPES:
C
      EXTERNAL LRL_BDAT
C
C====== LOCALS:
C   
      INTEGER I
      INTEGER J
C
C     INTEGER I        !Index to spot
C     INTEGER J        !Temp value for flag
C
C-------------------------------------------------------------------------------
C
      IERR = 0
      IF (ISPOT.LT.1.OR.ISPOT.GT.NUM_SPOTS) THEN
         IERR = 1
         RETURN
      END IF
      I = INDX(ISPOT)
      AII(I) = AI
      SIGII(I) = SIGI
      CALL STBITS (FLAGS(I),8,1,1)
      J = 0
      IF (BAD) J = 1
      CALL STBITS (FLAGS(I),9,1,J)
      J = 0
      IF (OVERLOAD) J = 1
      CALL STBITS (FLAGS(I),10,1,J)
      CALL STBITS (FLAGS(I),12,8,ICODE)
      RETURN
      END
C LAST UPDATE: 
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C**************************
C**     LRL_CLR_INT      **
C**************************
C
CD-Routine: Clear Intensity Data for a Spot - LRL_CLR_INT
C This routine will clear (re-initialise) the integrated intensity data for
C a spot in the Laue Reflections List.
CD-end
CD-Fortran:
      SUBROUTINE LRL_CLR_INT (ISPOT)
CD-end
      INCLUDE 'implicit_none.finc'
C
C Purpose: Clear (re-initialise) integrated intensity data for a spot
C
C Author:  John W. Campbell, July 1993
C
C Arguments:
C
      INTEGER ISPOT
CD-Parameters:
C ISPOT     i  (R)   Spot index number (in current sort order)
CD-end
C====== GLOBALS:
C
      INCLUDE 'lrl_data.finc'
C
C====== FUNCTION TYPES:
C
      EXTERNAL LRL_BDAT
C
C====== LOCALS:
C   
      INTEGER I
C
C     INTEGER I        !Index to spot
C
C-------------------------------------------------------------------------------
C
      IF (ISPOT.LT.1.OR.ISPOT.GT.NUM_SPOTS) RETURN
      I = INDX(ISPOT)
      AII(I) = 0.0
      SIGII(I) = 0.0
      CALL STBITS (FLAGS(I),8,1,0)
      CALL STBITS (FLAGS(I),9,1,0)
      CALL STBITS (FLAGS(I),10,1,0)
      CALL STBITS (FLAGS(I),12,8,0)
      CALL STBITS (FLAGS(I),20,1,0)
      CALL STBITS (FLAGS(I),21,8,0)
      RETURN
      END
C LAST UPDATE: 
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C**************************
C**     LRL_SET_SPOV     **
C**************************
C
CD-Routine: Set Spatially Overlapped Flag for a Spot - LRL_SET_SPOV
C This routine sets the spatially overlapped flag for a spot in the
C Laue Reflection List. It is basically written for internal use but
C may be called by the user if desired.
CD-end
CD-Fortran:
      SUBROUTINE LRL_SET_SPOV (ISPOT)
CD-end
      INCLUDE 'implicit_none.finc'
C
C Purpose: Set spatially overlapped flag for a spot
C
C Author:  John W. Campbell, July 1993
C
C Arguments:
C
      INTEGER ISPOT
CD-Parameters:
C ISPOT   i (R)   Spot index number (in current sort order)
CD-end
C====== GLOBALS:
C
      INCLUDE 'lrl_data.finc'
C
C====== FUNCTION TYPES:
C
      EXTERNAL LRL_BDAT
C
C====== LOCALS:
C
      INTEGER I
C
C-------------------------------------------------------------------------------
C 
      IF (ISPOT.LT.1.OR.ISPOT.GT.NUM_SPOTS) RETURN
      I = INDX(ISPOT)
      CALL STBITS (FLAGS(I),0,1,1)
      RETURN
      END
C LAST UPDATE: 
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C**************************
C**     LRL_SET_CLOS     **
C**************************
C
CD-Routine: Set 'Too Close' Flag for a Spot - LRL_SET_CLOS
C This routine sets the 'spots too close for spatial deconvolution' flag 
C for a spot. It is basically written for internal use but
C may be called by the user if desired.
CD-end
CD-Fortran:
      SUBROUTINE LRL_SET_CLOS (ISPOT)
CD-end
      INCLUDE 'implicit_none.finc'
C
C Purpose: Set 'spots too close for spatial deconvolution' flag for a spot
C
C Author:  John W. Campbell, July 1993
C
C Arguments:
C
      INTEGER ISPOT
CD-Parameters:
C ISPOT   i (R)   Spot index number (in current sort order)
CD-end
C====== GLOBALS:
C
      INCLUDE 'lrl_data.finc'
C
C====== FUNCTION TYPES:
C
      EXTERNAL LRL_BDAT
C
C====== LOCALS:
C
      INTEGER I
C
C-------------------------------------------------------------------------------
C 
      IF (ISPOT.LT.1.OR.ISPOT.GT.NUM_SPOTS) RETURN
      I = INDX(ISPOT)
      CALL STBITS (FLAGS(I),1,1,1)
      RETURN
      END
 
C LAST UPDATE: 
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C**************************
C**     LRL_SET_OVPIX    **
C**************************
C
CD-Routine: Set No. of Overload Pixels for a Spot - LRL_SET_OVPIX
C This routine sets the number of overloaded pixels in a spot for
C use within an integration routine for a spot in the Laue Reflections List.
CD-end
CD-Fortran:
      SUBROUTINE LRL_SET_OVPIX (ISPOT, NOVLD, IERR)
CD-end
       INCLUDE 'implicit_none.finc'
C
C Purpose: Sets the  no. of overloaded pixels if set 
C
C Author:  John W. Campbell, June 1997
C
C
      INTEGER ISPOT
      INTEGER NOVLD
      INTEGER IERR
CD-Parameters:
C ISPOT   i (R)   Spot index number (in current sort order)
C NOVLD   i (R)   The number of overloaded pixels in the spot. 
C                 If there are more than 255 pixels overloaded,
C                 a value of 255 will be stored. -1 clears the 
C                 count in the list.
C IERR    i (W)   Error flag =  0,  OK
C                            =  1,  Invalid spot number
CD-end
C====== GLOBALS:
C
      INCLUDE 'lrl_data.finc'
C
C====== FUNCTION TYPES:
C
      EXTERNAL LRL_BDAT
C
C====== LOCALS:
C
      INTEGER I
      INTEGER NOV
C
C-------------------------------------------------------------------------------
C
      IERR = 0
      IF (ISPOT.LT.1.OR.ISPOT.GT.NUM_SPOTS) THEN
         IERR = 1
         RETURN
      END IF
      I = INDX(ISPOT)
      IF (NOVLD.LT.0) THEN
         CALL STBITS (FLAGS(I),20,1,0)
         CALL STBITS (FLAGS(I),21,8,0)
      ELSE
         NOV = NOVLD
         IF (NOV.GT.255) NOV = 255
         CALL STBITS (FLAGS(I),20,1,1)
         CALL STBITS (FLAGS(I),21,8,NOV)
      END IF
      RETURN
      END 

C LAST UPDATE: 
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C**************************
C**     LRL_SET_OKNOD    **
C**************************
C
CD-Routine: Set Non-spatially Overlapped Nodal Flag - LRL_SET_OKNOD
C This routine sets the 'non-spatially overlapped nodal spot' flag for a spot.
C It is basically written for internal use but
C may be called by the user if desired.
CD-end
CD-Fortran:
      SUBROUTINE LRL_SET_OKNOD (ISPOT)
CD-end
      INCLUDE 'implicit_none.finc'
C
C Purpose: Set 'non-spatially overlapped nodal spot' flag for a spot
C
C Author:  John W. Campbell, July 1993
C
C Arguments:
C
      INTEGER ISPOT
CD-Parameters:
C ISPOT   i (R)   Spot index number (in current sort order)
CD-end
C====== GLOBALS:
C
      INCLUDE 'lrl_data.finc'
C
C====== FUNCTION TYPES:
C
      EXTERNAL LRL_BDAT
C
C====== LOCALS:
C
      INTEGER I
C
C-------------------------------------------------------------------------------
C 
      IF (ISPOT.LT.1.OR.ISPOT.GT.NUM_SPOTS) RETURN
      I = INDX(ISPOT)
      CALL STBITS (FLAGS(I),2,1,1)
      RETURN
      END

CD-Section: Set Global Laue Reflection List Data Items
C The routines described in this section are for setting parameters
C which apply to the Laue Reflection List as a whole.
CD-end

C LAST UPDATE: 
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C**************************
C**     LRL_CLR_ALL      **
C**************************
C
CD-Routine: Clear Intensity Data - LRL_CLR_ALL
C This routine clears (re-initialises) the intensities data for all the
C spots in the Laue Reflections List. The 'integration done' flag is
C also cleared.
CD-end
CD-Fortran:
      SUBROUTINE LRL_CLR_ALL
CD-end
      INCLUDE 'implicit_none.finc'
C
C Purpose: Clear (re-initialise) integrated intensity data for all spots
C          and clear 'integration done flag'
C
C Author:  John W. Campbell, July 1993
C
C Arguments:
CD-Parameters:
C     none
CD-end
C====== GLOBALS:
C
      INCLUDE 'lrl_data.finc'
C
C====== FUNCTION TYPES:
C
      EXTERNAL LRL_BDAT
C
C====== LOCALS:
C   
      INTEGER I
C
C     INTEGER I        !Index to spot
C
C-------------------------------------------------------------------------------
C
      INTEGRATED = .FALSE.
      INTCODE_TYPE = 0
      IF (NUM_SPOTS.LT.1) RETURN
      DO 100 I = 1, NUM_SPOTS
         AII(I) = 0.0
         SIGII(I) = 0.0
         CALL STBITS (FLAGS(I),8,1,0)
         CALL STBITS (FLAGS(I),9,1,0)
         CALL STBITS (FLAGS(I),10,1,0)
         CALL STBITS (FLAGS(I),12,8,0)
         CALL STBITS (FLAGS(I),20,1,0)
         CALL STBITS (FLAGS(I),21,8,0)
100   CONTINUE
      RETURN
      END
C LAST UPDATE: 
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C**************************
C**     LRL_SET_IDONE    **
C**************************
C
CD-Routine: Set the 'Integration Done' Flag - LRL_SET_IDONE
C This routine sets the 'integration done' flag for the Laue Reflections
C List and also stores the integration code type.
CD-end
CD-Fortran:
      SUBROUTINE LRL_SET_IDONE (ICODE_TYPE)
CD-end
      INCLUDE 'implicit_none.finc'
C
C Purpose: Set 'integration done' flag (and integration code type)
C
C Author:  John W. Campbell, July 1993
C
C Arguments:
C
      INTEGER ICODE_TYPE
CD-Parameters:
C ICODE_TYPE  i (R)  Integration code type 0 if not known to LRL routines
C                    (Used for interpreting the meanings of the 'ICODE'
C                    flags assigned when an intensity is integrated)
C                    =1 Integration using Q. Hao's routine LFN_SAINT
C                       as used in INTLDM/LAUEGEN
CD-end
C====== GLOBALS:
C
      INCLUDE 'lrl_data.finc'
C
C====== FUNCTION TYPES:
C
      EXTERNAL LRL_BDAT
C
C-------------------------------------------------------------------------------
C
      INTEGRATED = .TRUE.
      INTCODE_TYPE = ICODE_TYPE
      RETURN
      END

CD-Section: Determine Spatial Overlap Between Two Spots
C These routines are to determine whether two spots are spatially 
C overlapped. The elliptical LDM spot shapes are used and two different
C methods are available.
CD-end

C LAST UPDATE: 
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C**************************
C**     LRL_OVQ          **
C**************************
C
CD-Routine: Determine Spot Ellipse Overlap (analytical) - LRL_OVQ
C This routine determines whether two elliptical spots, with or without
C a background border, overlap using an anlytical procedure. The
C routine makes the approximation that the two ellipses are the
C same size and have the same orientation (reasonable as close
C spots are being considered).
CD-end
CD-Fortran:
      LOGICAL FUNCTION LRL_OVQ (SPOT_L, SPOT_W, SPOTFAC, SPOT_B, CTOF,
     +                          RASTOMM, X1, Y1, X2, Y2, ITYP)
CD-end
      INCLUDE 'implicit_none.finc'
C
C Purpose: Determine if two elliptical reflection masks overlap using
C          analytical procedure. 
C
C          Note: Makes approximations that the two ellipses are the
C                same size and orientation (reasonable as close
C                spots are being considered)
C
C Author:  Quan Hao, John W. Campbell, July 1993
C
C Return:  .TRUE. if spots overlap, otherwise .FALSE.
C
C Arguments:
C
      REAL SPOT_L
      REAL SPOT_W
      REAL SPOTFAC
      REAL SPOT_B
      REAL CTOF
      REAL RASTOMM
      REAL X1
      REAL Y1
      REAL X2
      REAL Y2
      INTEGER ITYP
CD-Parameters:   
C SPOT_L    r (R)  Spot length in mm.
C SPOT_W    r (R)  Spot width in mm.
C SPOTFAC   r (R)  Spot expansion factor
C SPOT_B    r (R)  Background border width in rasters
C CTOF      r (R)  Crystal to plate distance in mm.
C RASTOMM   r (R)  Conversion factor rasters to mm (x)
C X1        r (R)  'x' position of reference spot in mm. 
C Y1        r (R)  'y' position of reference spot in mm.
C X2        r (R)  'x' position of second spot in mm.
C Y2        r (R)  'y' position of second spot in mm.
C ITYP      i (R)  Flag = 1, Check that (spot + unused pixel + background 
C                            border) does not overlap with next spot 
C                       = 2, Check that (spot + background border) does not    
C                            overlap with next spot
C                       = 3, Check overlap strictly for spot area
C
C           Return is .TRUE. if spots overlap, otherwise .FALSE.
C
CD-end
C
C====== LOCALS:
C
      REAL CMASK
      REAL RADIUS
      REAL AMAJOR
      REAL BMINOR
      REAL CTH, STH
      REAL EXTRA
      REAL PI
      REAL XO
      REAL YO
      REAL XO2
      REAL YO2
      REAL A2
      REAL B2
      REAL CC
      REAL QTEST

      DATA PI/3.1415927/
C
C-------------------------------------------------------------------------------
C
      LRL_OVQ = .FALSE.
      CMASK = SPOT_L/2.0
      RADIUS = SQRT(X1*X1+Y1*Y1)
      AMAJOR = CMASK*(1.0+SPOTFAC*RADIUS/CTOF)
      BMINOR = SPOT_W/2.0
      EXTRA = 0.0
      IF (ITYP.EQ.1) THEN
         EXTRA = RASTOMM*(1.0 + SPOT_B)/2.0
      ELSE IF (ITYP.EQ.2) THEN
         EXTRA = RASTOMM*SPOT_B/2.0
      END IF
      CTH = X1/RADIUS
      STH = Y1/RADIUS
      XO = (X2-X1)*CTH + (Y2-Y1)*STH
      YO = (Y2-Y1)*CTH - (X2-X1)*STH
      AMAJOR = AMAJOR + EXTRA
      BMINOR = BMINOR + EXTRA
      A2 = AMAJOR*AMAJOR
      B2 = BMINOR*BMINOR
      XO2 = XO*XO
      YO2 = YO*YO
      CC = (B2*XO2 + A2*YO2)
      CC = CC*CC/(4.0*B2)
      QTEST = A2*B2*XO2 + A2*A2*YO2 - CC
      IF (QTEST.GT.0.0) LRL_OVQ = .TRUE.
      RETURN
      END
C LAST UPDATE: 
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C**************************
C**     LRL_OVM          **
C**************************
C
CD-Routine: Determine Spot Ellipse Overlap (using mask) - LRL_OVM
C This routine determines whether two elliptical spots, with or without
C a background border, overlap using a mask procedure based on code from
C the UPDATE program by T.J. Greenhough and A.K. Shrive.
CD-end
CD-Fortran:
      LOGICAL FUNCTION LRL_OVM (SPOT_L, SPOT_W, SPOTFAC, SPOT_B, CTOF,
     +                      MM_RAST_X, MM_RAST_Y, X1, Y1, X2, Y2, ITYP)
CD-end
      INCLUDE 'implicit_none.finc'
C
C Purpose: Determine if two elliptical reflection masks overlap using
C          mask procedure. 
C
C Author:  John W. Campbell, July 1993 (based on code from the program UPDATE
C                                       of T.J. Greenhough and A.K. Shrive)
C
C Return:  .TRUE. if spots overlap, otherwise .FALSE.
C
C Arguments:
C
      REAL SPOT_L
      REAL SPOT_W
      REAL SPOTFAC
      REAL SPOT_B
      REAL CTOF
      REAL MM_RAST_X
      REAL MM_RAST_Y
      REAL X1
      REAL Y1
      REAL X2
      REAL Y2
      INTEGER ITYP
CD-Parameters:   
C SPOT_L    r (R)  Spot length in mm.
C SPOT_W    r (R)  Spot width in mm.
C SPOTFAC   r (R)  Spot expansion factor
C SPOT_B    r (R)  Background border width in rasters
C CTOF      r (R)  Crystal to plate distance in mm.
C MM_RAST_X r (R)  Conversion factor mm to x-rasters 
C MM_RAST_Y r (R)  Conversion factor mm to y-rasters
C X1        r (R)  'x' position of reference spot in mm. 
C Y1        r (R)  'y' position of reference spot in mm.
C X2        r (R)  'x' position of second spot in mm.
C Y2        r (R)  'y' position of second spot in mm.
C ITYP      i (R)  Flag = 1, Check that (spot + unused pixel + background 
C                            border) does not overlap with next spot 
C                       = 2, Check that (spot + background border) does not
C                            overlap with next spot
C                       = 3, Check overlap strictly for spot area
C
C           Returns .TRUE. if spots overlap, otherwise .FALSE.
C
CD-end
C
C====== LOCALS:
C
      INTEGER IX1
      INTEGER IY1
      INTEGER IX2
      INTEGER IY2
      INTEGER IXS
      INTEGER IXE
      INTEGER IYS
      INTEGER IYE
      INTEGER IOFF
      INTEGER IX
      INTEGER IY
      REAL CMASK
      REAL BMINOR
      REAL AMAJOR
      REAL RADIUS
      REAL AM1
      REAL BM1
      REAL AM2
      REAL BM2
      REAL STH1
      REAL CTH1
      REAL STH2
      REAL CTH2
      REAL XPRIME
      REAL YPRIME
      REAL ELLIP

C     INTEGER IX1       !X ccordinate of 1'st spot centre (rasters)
C     INTEGER IY1       !Y coordinate of 1'st spot centre (rasters)
C     INTEGER IX2       !X coordinate of 2'nd spot centre (rasters)
C     INTEGER IY2       !Y coordinate of 2'nd spot centre (rasters)
C     INTEGER IXS       !X start of box surrounding 1'st spot (rasters)
C     INTEGER IXE       !Y start of box surrounding 1'st spot (rasters)
C     INTEGER IYS       !X start of box surrounding 2'nd spot (rasters)
C     INTEGER IYE       !Y start of box surrounding 2'nd spot (rasters)
C     INTEGER IOFF      !Max. offset from centre of box (rasters)
C     INTEGER IX        !Loop x counter for box
C     INTEGER IY        !Loop y counter for box
C     REAL CMASK        !Major axis of ellipse uncorrected for radius (rasters)
C     REAL BMINOR       !Minor axis of ellipse (rasters)
C     REAL AMAJOR       !Corrected major axis of ellipse (rasters)
C     REAL RADIUS       !Radius of spot from centre (rasters)
C     REAL AM1          !1/(amajor1*amajor1)
C     REAL BM1          !1/(bminor1*bminor1)
C     REAL AM2          !1/(amajor2*amajor2)
C     REAL BM2          !1/(bminor2*bminor2)
C     REAL STH1         !sin(theta) for 1'st spot
C     REAL CTH1         !cos(theta) for 1'st spot
C     REAL STH2         !sin(theta) for 2'nd spot
C     REAL CTH2         !sin(theta) for 2'nd spot
C     REAL XPRIME       !x coordinate wrt local ellipse axes (rasters)
C     REAL YPRIME       !y coordinate wrt local ellipse axes (rasters)
C     REAL ELLIP        !Ellipse function test value
C
C
C-------------------------------------------------------------------------------
C
      LRL_OVM = .FALSE.
      CMASK = SPOT_L*MM_RAST_X/2.0
C
C====== Get basic parameters for 1'st spot
C
      IX1 = NINT(X1*MM_RAST_X)
      IY1 = NINT(Y1*MM_RAST_Y)
      RADIUS = SQRT(X1*X1 + Y1*Y1)
      AMAJOR = CMASK*(1.0+SPOTFAC*RADIUS/
     +         (CTOF))
      BMINOR = SPOT_W*MM_RAST_X/2.0
      IF (ITYP.EQ.1) THEN
         AMAJOR = AMAJOR + 1.0 + SPOT_B
         BMINOR = BMINOR + 1.0 + SPOT_B
      ELSE IF (ITYP.EQ.2) THEN
         AMAJOR = AMAJOR + SPOT_B
         BMINOR = BMINOR + SPOT_B
      END IF
      AM1 = 1.0/(AMAJOR*AMAJOR)
      BM1 = 1.0/(BMINOR*BMINOR)
      IOFF = NINT(AMAJOR+1.0)
      IXS = IX1 - IOFF
      IXE = IX1 + IOFF
      IYS = IY1 - IOFF
      IYE = IY1 + IOFF
      STH1 = Y1/RADIUS
      CTH1 = X1/RADIUS
C
C====== Get basic parameters for 2'nd spot
C
      IX2 = NINT(X2*MM_RAST_X)
      IY2 = NINT(Y2*MM_RAST_Y)
      RADIUS = SQRT(X2*X2 + Y2*Y2)
      AMAJOR = CMASK*(1.0+SPOTFAC*RADIUS/
     +         (CTOF))
      BMINOR = SPOT_W*MM_RAST_X/2.0
      AM2 = 1.0/(AMAJOR*AMAJOR)
      BM2 = 1.0/(BMINOR*BMINOR)
      STH2 = Y2/RADIUS
      CTH2 = X2/RADIUS
C
C====== Check for overlap
C
      DO 150 IX = IXS, IXE
         DO 100 IY = IYS, IYE
            XPRIME = (IX-IX1)*CTH1 + (IY-IY1)*STH1
            YPRIME = (IY-IY1)*CTH1 - (IX-IX1)*STH1
            ELLIP = XPRIME*XPRIME*AM1 + YPRIME*YPRIME*BM1
            IF (ELLIP.LE.1.0) THEN
               XPRIME = (IX-IX2)*CTH2 + (IY-IY2)*STH2
               YPRIME = (IY-IY2)*CTH2 - (IX-IX2)*STH2
               ELLIP = XPRIME*XPRIME*AM2 + YPRIME*YPRIME*BM2
               IF (ELLIP.LE.1.0) THEN
                  LRL_OVM = .TRUE.
                  RETURN
               END IF
            END IF
100      CONTINUE
150   CONTINUE
      RETURN
      END
C LAST UPDATE: 
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C**************************
C**     LRL_GETHKL       **
C**************************
C
C
      SUBROUTINE LRL_GETHKL (HKL, IH, IK, IL)
C
      INCLUDE 'implicit_none.finc'
C
C Purpose: Unpack packed hkl or nodal hkl as stored in reflection list
C
C Author:  John W. Campbell, July 1993
C
C **** Service Routine Only ****
C
C Arguments:
C
      INTEGER HKL
      INTEGER IH
      INTEGER IK
      INTEGER IL
C
C HKL     (R)   Packed indices
C IH      (W)   Returns 'h' index
C IK      (W)   Returns 'k' index
C IL      (W)   Returns 'l' index
C
C
C====== LOCALS:
C
      INTEGER I
      INTEGER NTWO20
      INTEGER NTWO10

C     INTEGER I           !Temp. integer value
C     INTEGER NTWO20      !2**20
C     INTEGER NTWO10      !2**10

      DATA NTWO20/1048576/
      DATA NTWO10/1024/
C
C-------------------------------------------------------------------------------
C
      I = HKL
      IH = I/NTWO20
      I = I - NTWO20*IH
      IH = IH - 512
      IK = I/NTWO10
      IL = I - NTWO10*IK - 512
      IK = IK - 512
      RETURN
      END
C LAST UPDATE: 
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C**************************
C**     LRL_GETMULT      **
C**************************
C
C
      SUBROUTINE LRL_GETMULT (MULTFLAG, MULT, MINHARM, MAXHARM, INCHARM)      
C
      INCLUDE 'implicit_none.finc'
C
C Purpose: Unpack multiplicities flag as stored in reflection list
C
C Author:  John W. Campbell, July 1993
C
C **** Service Routine Only ****
C
C Arguments:
C
      INTEGER MULTFLAG
      INTEGER MULT
      INTEGER MINHARM
      INTEGER MAXHARM
      INTEGER INCHARM
C
C MULTFLAG  (R)   Packed 'LRL' multiplicities flag
C MULT      (W)   Returns multiplicity
C MINHARM   (W)   Returns minimum harmonic number (0 if single)
C MAXHARM   (W)   Returns maximum harmonic number (0 if single)
C INCHARM   (W)   Returns harmonic increment (0 if single)
C
C
C====== LOCALS:
C
      INTEGER I
      INTEGER NTWO20
      INTEGER NTWO12
      INTEGER NTWO4

C     INTEGER I           !Temp. integer value
C     INTEGER NTWO20      !2**20
C     INTEGER NTWO12      !2**12
C     INTEGER NTWO4       !2**4

      DATA NTWO20/1048576/
      DATA NTWO12/4096/
      DATA NTWO4/16/
C
C-------------------------------------------------------------------------------
C
      MULT = 0
      MINHARM = 0
      MAXHARM = 0
      INCHARM = 0
      IF (MULTFLAG.EQ.0) RETURN
      IF (MULTFLAG.EQ.1) THEN
         MULT = 1
         RETURN
      END IF
      I = MULTFLAG
      MULT = I/NTWO20
      I = I - NTWO20*MULT
      MINHARM = I/NTWO12
      I = I - NTWO12*MINHARM
      MAXHARM = I/NTWO4
      INCHARM = I - NTWO4*MAXHARM
      RETURN
      END
C LAST UPDATE: 
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C**************************
C**     LRL_BDAT         **
C**************************
C
C
      BLOCK DATA LRL_BDAT
C
      INCLUDE 'implicit_none.finc'
C
C Purpose: Initialise LRL data
C
C Author:  John W. Campbell, July 1993
C
C **** Service Routine Only ****
C
C====== GLOBALS:
C
      INCLUDE 'lrl_data.finc'

      DATA NUM_SPOTS/0/, NSING/0/, NMULT/0/, NUM_REFLEX/0/
      DATA NTRIES/0/, LASTHKL/3*0/, XSTART/-1.0E09/, XEND/1.0E09/
      DATA INTEGRATED /.FALSE./ DISTOR_CORR/.TRUE./
      DATA OVLAP_STATUS/0/, ISORT /0/, INTCODE_TYPE/0/
C
C-------------------------------------------------------------------------------
C
      END

CD-Section: Miscellaneous Routines
C This section describes some routines which are basically for internal use
C but may also be called by the user if desired.
CD-end

C LAST UPDATE:
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C*************************
C**       LRL_HCF       **
C*************************
C
CD-Routine: Highest Common Factor of Two Numbers - LRL_HCF
C This routine, based on code by I.J. Clifton, returns the highest common
C factor of two integers.
CD-end
CD-Fortran:
      INTEGER FUNCTION  LRL_HCF (II,JJ)
CD-end
      INCLUDE 'implicit_none.finc'
C
C Purpose: To calculate h.c.f. of 2 integers
C
C Author:  Based on code from GENLAUE by Clifton et al.
C
C Arguments:
C
      INTEGER    II
      INTEGER    JJ
CD-Parameters: 
C II    i (R)   First integer
C JJ    i (R)   Second integer
C
C       Returns the highest common factor of the two integers
C
CD-end
C
C====== LOCALS:
C
      INTEGER    I, J, ISAV, K
C
C------------------------------------------------------------------------------
C
      I = II
      J = JJ
      IF (I.EQ.J) THEN
         LRL_HCF = I
         RETURN
      END IF

      IF (I.LT.J) THEN
         ISAV = I
         I = J
         J = ISAV
      END IF

      IF (J.EQ.0) THEN
         LRL_HCF = I
         RETURN
      END IF

 20   K = MOD(I,J)

      IF (K.EQ.0) THEN
         LRL_HCF = J
         RETURN
      END IF

      I = J
      J = K
      GO TO 20

      END
C LAST UPDATE:
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C*************************
C**       LRL_RECCEL    **
C*************************
C
CD-Routine: Convert Between Real and Reciprocal Cells - LRL_RECCEL
C This routine converts between real and reciprocal cell parameters.
CD-end
CD-Fortran:
      SUBROUTINE LRL_RECCEL(RX ,CX ,WAVE)
CD-end
C Purpose: Convert between real and reciprocal cell parameters. CX input,
C          RX output; angles in degrees
C
CD-Parameters:
C RX(6)   r  (W) Converted cell parameters reciprocal or real
C CX(6)   r  (R) Input cell parameters real or reciprocal
C WAVE    r  (R) Standard wavelength
CD-end

      DIMENSION CX(6),RX(6)

      REAL C(6)
      REAL DTOR

C
C-------------------------------------------------------------------------------
C
      
       DTOR = 0.0174533

      DO 10 I = 1, 3
         C(I) = CX(I)
 10   CONTINUE

      SUM = 0.
      DO 20 I = 4 , 6
         C(I) = DTOR*CX(I)
         SUM = SUM + 0.5*C(I)
 20   CONTINUE
 
      V = 1.
      DO 30 I = 4 , 6
         V = V*SIN(SUM-C(I))
 30   CONTINUE
 
      WV = WAVE/(2.*C(1)*C(2)*C(3)*SQRT(SIN(SUM)*V))

      RX(1) = C(2)*C(3)*SIN(C(4))*WV
      RX(2) = C(1)*C(3)*SIN(C(5))*WV
      RX(3) = C(1)*C(2)*SIN(C(6))*WV
 
      RX(4) = ACOS((COS(C(5))*COS(C(6))-COS(C(4)))
     +       /(SIN(C(5))*SIN(C(6))))
      RX(4) = RX(4)/DTOR

      RX(5) = ACOS((COS(C(6))*COS(C(4))-COS(C(5)))
     +       /(SIN(C(6))*SIN(C(4))))
      RX(5) = RX(5)/DTOR

      RX(6) = ACOS((COS(C(4))*COS(C(5))-COS(C(6)))
     +       /(SIN(C(4))*SIN(C(5))))
      RX(6) = RX(6)/DTOR
 
      RETURN
      END
C LAST UPDATE:
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C***************************
C**       LRL_NODIDX      **
C***************************
C
CD-Routine: Get Nodal Index from Indices - LRL_NODIDX
C This routine determines the nodal index from a set of reflection indices.
CD-end
CD-Fortran:
      SUBROUTINE LRL_NODIDX (IH, IK, IL, NODIDX)
CD-end
C
C      INCLUDE 'implicit_none.finc'
C
C Purpose: To calculate the "nodal spot selection index" at which a
C          nodal would be selected.
C
C Author:  Based on code from GENLAUE by Clifton et al.
C
C Arguments:
C
      INTEGER IH
      INTEGER IK
      INTEGER IL
      INTEGER NODIDX
C
CD-Parameters:
C IH      i  (R) 'h' index
C IK      i  (R) 'k' index
C IL      i  (R) 'l' index
C NODIDX  i  (W) Returns nodal index
CD-end
C
C====== FUNCTION TYPES:
C
      INTEGER LRL_HCF
C
C====== LOCALS:
C
      INTEGER   HCF1, HCF2
      INTEGER   H1, K1, L1
      INTEGER   HI, LO
C
C------------------------------------------------------------------------------
C
      H1 = IH
      K1 = IK
      L1 = IL
      HCF1 = LRL_HCF(H1,K1)
      HCF2 = LRL_HCF(HCF1,L1)
      H1 = H1/HCF2
      K1 = K1/HCF2
      L1 = L1/HCF2
C
      HI = MAX0(H1,K1,L1)
      LO = MIN0(H1,K1,L1)
      NODIDX = MAX0(HI,-LO)
C
      RETURN
      END
C LAST UPDATE: 
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C**************************
C**     LRL_STILLOK      **
C**************************
C
C
      LOGICAL FUNCTION LRL_STILLOK (II)
C
      INCLUDE 'implicit_none.finc'
C
C Purpose: See if spot still predicted for revised soft limits
C
C Author:  John W. Campbell, October 1993
C
C  **** LRL Service Routine Only ***
C
C Arguments:
C
      INTEGER II
C
C II    (R)   Index to reflection list arrays.
C
C====== Parameter statements
C 

C
C====== GLOBALS:
C
      INCLUDE 'lrl_data.finc'
C
C====== FUNCTION TYPES:
C

C
C====== LOCALS:
C
      LOGICAL OK
      INTEGER IHARM
      INTEGER MULT
      INTEGER MINHARM
      INTEGER MAXHARM
      INTEGER INCHARM
      REAL LAMBDA
      REAL DTHR2

C     LOGICAL OK           !Current harmonic OK flag
C     INTEGER IHARM        !Harmonics loop counter
C     INTEGER MULT         !Multiplicity
C     INTEGER MINHARM      !Minimum harmonic
C     INTEGER MAXHARM      !Maximum harmonic
C     INTEGER INCHARM      !Harmonic increment
C     REAL LAMBDA          !Lambda
C     REAL DTHR2           !1/dmin**2
C
C-------------------------------------------------------------------------------
C
      LRL_STILLOK = .TRUE.
      IF (MULTS(II).EQ.1) THEN
         IF (LAMBDAS(II).LT.LMIN_NEW.OR.LAMBDAS(II).GT.LMAX_NEW) THEN
            LRL_STILLOK = .FALSE.
            RETURN
         END IF
         IF (DMIN_THR2(II).GT.DTHR2_NEW) THEN
            LRL_STILLOK = .FALSE.
            RETURN
         END IF
      ELSE
         CALL LRL_GETMULT (MULTS(II), MULT, MINHARM, MAXHARM, INCHARM)
         DO 10 IHARM = MINHARM, MAXHARM, INCHARM
            OK = .TRUE.
            LAMBDA = LAMBDAS(II) * FLOAT(MINHARM) / FLOAT(IHARM)
            DTHR2 = DMIN_THR2(II) * FLOAT(IHARM) / FLOAT(MINHARM)
            IF (LAMBDA.LT.LMIN_NEW.OR.LAMBDA.GT.LMAX_NEW) OK = .FALSE.
            IF (DTHR2.GT.DTHR2_NEW) OK = .FALSE.
            IF (OK) RETURN
10       CONTINUE
         LRL_STILLOK = .FALSE.
      END IF
      RETURN
      END
