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
      REAL*8 D2RAD
      PARAMETER ( D2RAD = 3.14159265358979/180 )
      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 = COS(W_C*D2RAD)
      GSIN_WC = SIN(W_C*D2RAD)
      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
