C LAST UPDATE:
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C***************************
C**     LSC_RFACTS        **
C***************************
C
C **** LSC SERVICE ROUTINE ONLY ****
C
C Routine: Calculate merging R-factors for LSM - LSC_RFACTS
C This routine is called by LSC_RFACTS_LSP and calculates various weighted 
C and unweighted R-factors according to input in- or exclusion criterias
C and prepares the data for the R-factor tables to be output by LSC_RFACTS_OUT 
C
C Note:  OVERLOADS excluded
C
      SUBROUTINE LSC_RFACTS (LINDX, NBIN, SIG_CUT, SPOV, 
     +                       PLATEMERGE, OVLD, IERR)
C
      INCLUDE 'implicit_none.finc'
      INCLUDE 'limits.finc'
C
C Purpose: Calculate Merging R-factors for LSM
C
C Author:  Steffi Arzt, John Campbell, July 1997
C
C Arguments:
C
      INTEGER LINDX 
      INTEGER NBIN(MAXRANGE)
      REAL SIG_CUT
      LOGICAL SPOV
      LOGICAL OVLD
      LOGICAL PLATEMERGE
      INTEGER IERR
C
C Parameters:
C  LINDX   i  (R)   The index of the Laue Integrated Reflections List
C  NBIN()  i  (R)   Number of bins per wavelength range
C  SIG_CUT i  (R)   Reflections with I<SIG_CUT*sig(I) will be
C                   excluded from the analyses
C  SPOV    l  (R)   Include spatial overlaps in analyses flag 
C                   = .TRUE. yes, = .FALSE. no
C  PLATEMERGE l(R)  .true. = Merge plate data prior to statistics
C                   .false. = keep all data seperate
C  OVLD    l  (R)   Includeoverloads in analyses flag 
C                   = .TRUE. yes, = .FALSE. no
C  IERR    i  (W)   Error flag =0 OK
C                              =1 Invalid list index
C end
C
C====== PARAMETERS:
C
C
C====== GLOBALS:
C
      INCLUDE 'scp.finc'
      INCLUDE 'scps_counters.finc'
      INCLUDE 'next_refl.finc'
      INCLUDE 'next_refobs.finc'
      INCLUDE 'x_array.finc'
      INCLUDE 'rout_arrays.finc'
      INCLUDE 'scale_presets.finc'
      INCLUDE 'other_scale.finc'
      INCLUDE 'refine_lsm.finc'
      INCLUDE 'ref_flags.finc'
C
C====== FUNCTION TYPES:
C
      INTEGER LSM_NEXT_REFLN
      INTEGER LSM_INDX
C
C====== LOCALS:
C
      LOGICAL ABSC_USED
      INTEGER I, J, K, N, MM, JM, IM, IBIN, JBIN, IPACK, JPACK
      INTEGER NSUM, NXPAR
      INTEGER IH, IK, IL
      INTEGER KPLS
      INTEGER NBTOT
      INTEGER NR
      INTEGER NI(3)
      INTEGER IPLATE, KPLATE(MAX_SYMEQ)
      INTEGER KSM(MAX_SYMEQ)
      INTEGER KBIN(MAX_SYMEQ), KPACK(MAX_SYMEQ), MODSYM(MAX_SYMEQ)
      INTEGER KERR
      INTEGER KNORM, IUSE
      INTEGER KSYM
      REAL SCAL_INT(MAX_SYMEQ), SCAL_SIG(MAX_SYMEQ) 
      REAL ALAMS(MAX_SYMEQ)
      REAL OB_LPFAC
      REAL AI, SIGMA
      REAL SC
      REAL ABSCOR_NU, PSI_LOC, PSI_GLO
      REAL STOL2, AL3C2T, ABSCOR_LAMFN  
      REAL ALAM
      REAL REFINT
      REAL AMINL(MAXRANGE), AMAXL(MAXRANGE) 
      REAL BMINL(MAXRANGE,MAXRFR), BMAXL(MAXRANGE,MAXRFR)
      REAL BINW( MAXRANGE)
      REAL RSUM1(MAXPACK,MAXRFR)
      REAL RSUM2(MAXPACK,MAXRFR)
      REAL RSUM3(MAXPACK,MAXRFR)
      REAL RSUMP(MAXPACK,MAXRFR)
      REAL SIGI(3), RRD(3), AV(3), DL(3), RGOF(3)
      REAL SIGWI(3), SIGW(3),  RD(3), AVE(3), DEL(3)
      REAL WWI
      LOGICAL RF(3)
C     INTEGER I, MM         ! Temp/loop variables 
C     INTEGER NR            ! No. of reflections in set of equivalents
C     INTEGER KERR          ! Error flag
C     REAL SUMD, SUMI       ! Temp/loop variables 
C     REAL RFACT            ! Calculated R-factor
C     REAL OB_LPFAC         ! Lorenz-polarisation/obliquity factor
C     LOGICAL REFERENCE     ! Flag = .true.  return reference intensities
C                           ! = .false. return weighted mean scaled intensities
C
C
C====== Initialisations
C
      MM = 0
      NSUM = 0
      NRANGE = 0
      NBTOT = 0
C
C====== Calculate No. of lambda-bins NBTOT
C
      DO 200 I = 1, MAXRANGE
         IF (NORMLAMMIN(I).GT.(0.0).AND.
     +       NORMLAMMAX(I).GT.(0.0)) THEN
            NBINR(I) = NBIN(I)
            NRANGE = NRANGE + 1
         ELSE
            NBINR(I) = 0
         END IF
 200  CONTINUE   

      DO 260 N = 1, MAXRANGE
         NBTOT = NBTOT + NBINR(N)
 260  CONTINUE

      DO 30 I = 1 , NBTOT
         DO 20 J = 1 , NBTOT
            RFAC1(I,J) = 0.0
            RSUM1(I,J) = 0.0
            NREF1(I,J) = 0
            RFAC2(I,J) = 0.0
            RSUM2(I,J) = 0.0
            NREF2(I,J) = 0
            RFAC3(I,J) = 0.0
            RSUM3(I,J) = 0.0
            NREF3(I,J) = 0
 20      CONTINUE
 30   CONTINUE
      DO 40 I = 1 , 3
         RN(I) = 0.0
         RRN(I) = 0.0
         RD(I) = 0.0
         RRD(I) = 0.0
         RGOF(I) = 0.0
         NRR(I) = 0
 40   CONTINUE

      DO 60 I = 1 , NPK
         DO 50 J = 1 , NPK
            RFACP(I,J) = 0.0
            RSUMP(I,J) = 0.0
            NREFP(I,J) = 0
 50      CONTINUE
 60   CONTINUE

      DO 180 I = 1, MAXRANGE
         DO 190 J = 1, NBINR(I)
            AMINL(I) = NORMLAMMIN(I)
            AMAXL(I) = NORMLAMMAX(I)
            BINW(I) = (AMAXL(I)-AMINL(I))/NBINR(I)
            BMINL(I,J) = AMINL(I) + (J-1)*BINW(I)
            BMAXL(I,J) = AMINL(I) + J*BINW(I)              
 190     CONTINUE
 180  CONTINUE
C
C====== Save  current selection flags
C
      CALL LSM_SEL_SAVE 
C
C====== Select all (except omitted) packs
C
      CALL LSM_SELECT_PKPL (0, 0) 
C
C====== Select required reflection categories for use in the analyses
C
      CALL LSM_SELECT_SCREFS (.FALSE., .FALSE., SIG_CUT, SPOV, 
     +                        OVLD) 
    
      CALL LSM_TOP_LIRL(LINDX, KERR)

      IF (KERR.NE.0) THEN
         IERR = 1
         RETURN
      END IF
 
      CALL LSM_CALC_PRESETS 

 10   NR = LSM_NEXT_REFLN(KERR)

      IF (NR.LE.-1) THEN
C
C====== End of reflections reached, prepare tables for R-factor
C
C
C====== Calculate R-factors
C
         DO 140 I = 1 , NBTOT
            DO 130 J = 1 , NBTOT
               IF (NREF1(I,J).GT.0) RFAC1(I,J) = RFAC1(I,J)/RSUM1(I,J)
               IF (NREF2(I,J).GT.0) RFAC2(I,J) = RFAC2(I,J)/RSUM2(I,J)
               IF (NREF3(I,J).GT.0) RFAC3(I,J) = RFAC3(I,J)/RSUM3(I,J)
 130        CONTINUE
 140     CONTINUE
         DO 160 I = 1 , NPK
            DO 150 J = 1 , NPK
               IF (NREFP(I,J).GT.0) RFACP(I,J) = RFACP(I,J)/RSUMP(I,J)
 150        CONTINUE
 160     CONTINUE
C
C====== Calculate number of non-default parameters
C
         NXPAR = 0
         DO 300 I = 1, MAXRANGE
            IF (PLSCRAN(I)) THEN
               DO 310 J = 1, NPK
                  NXPAR = NXPAR + NPL(J)-1
 310           CONTINUE
            END IF
 300     CONTINUE
         DO 70 I = 2, NPK
            IF(SCP(NCSCP(2)+I).NE.1.0) NXPAR = NXPAR + 1
            IF(SCP(NCSCP(3)+I).NE.0.0) NXPAR = NXPAR + 1
 70      CONTINUE
         DO 330 J = 1, MAXRANGE
            IF (NORMRAN(J)) NXPAR = NXPAR + NNORM(J)
 330     CONTINUE
         IF (REFSETFLAG) NXPAR = NXPAR + 2
         IF ((NABSG-1.NE.0.OR.MABSG-1.NE.0).AND.ABSCORGLOB) THEN
            ABSC_USED = .FALSE.
            DO 335 I = 1, NABSG*MABSG-1
               IF (SCP(NCSCP(7)+I).NE.0.0) ABSC_USED = .TRUE.
 335        CONTINUE
            IF (ABSC_USED) THEN
               IF (IOPTABSG.EQ.1) NXPAR = NXPAR + NABSG*MABSG-1
               IF (IOPTABSG.EQ.2) NXPAR = NXPAR + 2*NABSG*MABSG-2
            END IF
         END IF
         DO 340 I = 1, NPK
            IF ((NABSL(IPACK)-1.NE.0.OR.MABSL(IPACK)-1.NE.0).AND.
     +         ABSCORLOC) THEN
               ABSC_USED = .FALSE.
               DO 338 J = 1, NABSL(I)*MABSL(I)-1
                   IF (SCPLOC(J+(I-1)*MAXCOEFFS).NE.0.0) 
     +                 ABSC_USED = .TRUE.
 338           CONTINUE
               IF (ABSC_USED) NXPAR = NXPAR + NABSL(I)*MABSL(I)-1
            END IF
 340     CONTINUE

         DO 170 I = 1 , 3
            IF (NRR(I).GT.0) THEN
C
C====== RN() weighted R-factors, RRN() unweighted R-factors
C  
               RN(I) = RN(I)/RD(I)
               RRN(I) = RRN(I)/RRD(I)
               GOFR(I) = SQRT(RGOF(I)/(MM-NXPAR))
            END IF
 170     CONTINUE

         NPK_RF = NPK

         CALL LSM_SEL_RESTORE

         RETURN

      END IF

C
C====== Get reflection details for set
C      
      IF (NR.LT.2) GO TO 10

      DO 240 I = 1, NR
         CALL LIRLF_LSM_SCALOBS (LINDX,LSM_INDX(I), 
     +             IPACK, IPLATE, AI, SIGMA, ALAM, REFINT,
     +             ABSCOR_NU, PSI_LOC, PSI_GLO, STOL2,
     +             AL3C2T, ABSCOR_LAMFN, OB_LPFAC,
     +             KNORM, KPLS, IUSE, KERR)
         CALL LSM_CALCSCFAC (SC, IPACK, IPLATE, KPLS, KNORM,
     +             ABSCOR_NU, PSI_LOC, PSI_GLO, ABSCOR_LAMFN,  
     +             STOL2, ALAM, AL3C2T, OB_LPFAC,  .TRUE.,
     +             PQ00_GLOB_STORE, P00_LOC_STORE, AKNORM_STORE,
     +             ANORMREF_STORE, KERR)
         SCAL_INT(I) = SC*AI
         SCAL_SIG(I) = SC*SIGMA
         KPLATE(I) = IPLATE
         KPACK(I) = IPACK
         ALAMS(I) = ALAM
         CALL LIRLF_GETHKL (LINDX, LSM_INDX(I), IH, IK, IL, KSYM, KERR)
         MODSYM(I) = MOD(KSYM,2)
         KSM(I) = KSYM
 240  CONTINUE

      IF (PLATEMERGE.AND.PLFLAG) THEN
         CALL LSM_MERGE_PLATES (NR, KPACK, KPLATE, KSM, MODSYM,
     +                             SCAL_INT, SCAL_SIG, ALAMS)
      END IF   

      DO 210 K = 1, NR
         DO 220 I = 1, MAXRANGE
            DO 230 J = 1, NBINR(I)
               IF (ALAMS(K).GE.BMINL(I,J).AND. 
     +             ALAMS(K).LE.BMAXL(I,J)) THEN
                  NSUM = 0
                  IF (I.GT.1) THEN
                     DO 250 N = 1, I-1 
                        NSUM = NSUM + NBINR(N)
 250                 CONTINUE
                  END IF
                  KBIN(K) = NSUM + J
                  GO TO 210
               END IF
 230        CONTINUE
  220     CONTINUE
 210  CONTINUE
      
C
C====== Step through reflections
C
      DO 120 IM = 1, NR
         MM = MM + 1
         IBIN = KBIN(IM)
         IPACK = KPACK(IM)
         DO 90 I = 1 , 3
            SIGWI(I) = 0.0
            SIGW(I) = 0.0
            SIGI(I) = 0.0
            NI(I) = 0
            RF(I) = .FALSE.
 90      CONTINUE
C
C====== Find the mean intensities for the reflection required for calculating
C====== the R-factors. the means are (1) the mean of all measurements,
C====== (2) the mean of all measurements of the same sign as the current
C====== reflection (i.e. i+ or i-) and (3) the mean of measurements with the
C====== same sign and measured withinn 0.1 Angstrom of the same wavelength.
C
         DO 100 JM = 1, NR
            IF (IM.NE.JM) RF(1) = .TRUE.
            JBIN = KBIN(JM)
            JPACK = KPACK(JM)
            WWI = 1/(SCAL_SIG(JM)*SCAL_SIG(JM))
            SIGWI(1) = SIGWI(1) + WWI*SCAL_INT(JM)
            SIGW(1) = SIGW(1) + WWI
            SIGI(1) = SIGI(1) + SCAL_INT(JM)
            NI(1) = NI(1) + 1
            IF (IM.LT.JM) THEN
               RFAC1(IBIN,JBIN) = RFAC1(IBIN,JBIN) + 
     +                            ABS(SCAL_INT(IM)-SCAL_INT(JM))
               RSUM1(IBIN,JBIN) = RSUM1(IBIN,JBIN) + 
     +                            SCAL_INT(IM) + SCAL_INT(JM)
               NREF1(IBIN,JBIN) = NREF1(IBIN,JBIN) + 1
               RFAC1(JBIN,IBIN) = RFAC1(IBIN,JBIN)
               RSUM1(JBIN,IBIN) = RSUM1(IBIN,JBIN)
               NREF1(JBIN,IBIN) = NREF1(IBIN,JBIN)
               RFACP(IPACK,JPACK) = RFACP(IPACK,JPACK) + 
     +                            ABS(SCAL_INT(IM)-SCAL_INT(JM))
               RSUMP(IPACK,JPACK) = RSUMP(IPACK,JPACK) +
     +                            SCAL_INT(IM) + SCAL_INT(JM)
               NREFP(IPACK,JPACK) = NREFP(IPACK,JPACK) + 1
               RFACP(JPACK,IPACK) = RFACP(IPACK,JPACK)
               RSUMP(JPACK,IPACK) = RSUMP(IPACK,JPACK)
               NREFP(JPACK,IPACK) = NREFP(IPACK,JPACK)
            END IF

C
C====== Same sign
C
            IF (MODSYM(IM).EQ.MODSYM(JM)) THEN
               SIGWI(2) = SIGWI(2) + WWI*SCAL_INT(JM)
               SIGW(2) = SIGW(2) + WWI
               SIGI(2) = SIGI(2) + SCAL_INT(JM)
               NI(2) = NI(2) + 1
               IF (IM.NE.JM) RF(2) = .TRUE.
               IF (IM.LT.JM) THEN
                  RFAC2(IBIN,JBIN) = RFAC2(IBIN,JBIN) + 
     +                               ABS(SCAL_INT(IM)-SCAL_INT(JM))
                  RSUM2(IBIN,JBIN) = RSUM2(IBIN,JBIN) + 
     +                               SCAL_INT(IM) + SCAL_INT(JM)
                  NREF2(IBIN,JBIN) = NREF2(IBIN,JBIN) + 1
                  RFAC2(JBIN,IBIN) = RFAC2(IBIN,JBIN)
                  RSUM2(JBIN,IBIN) = RSUM2(IBIN,JBIN)
                  NREF2(JBIN,IBIN) = NREF2(IBIN,JBIN)
               END IF
C
C====== Same sign and lambdas within 0.1 Angstrom
C
               IF (ABS(ALAMS(JM)-ALAMS(IM)).LE.0.1) THEN
                  SIGWI(3) = SIGWI(3) + WWI*SCAL_INT(JM)
                  SIGW(3) = SIGW(3) + WWI
                  SIGI(3) = SIGI(3) + SCAL_INT(JM)
                  NI(3) = NI(3) + 1
                  IF (IM.NE.JM) RF(3) = .TRUE.
                  IF (IM.LT.JM) THEN
                     RFAC3(IBIN,JBIN) = RFAC3(IBIN,JBIN) +
     +                    ABS(SCAL_INT(IM)-SCAL_INT(JM))
                     RSUM3(IBIN,JBIN) = RSUM3(IBIN,JBIN) + 
     +                    SCAL_INT(IM) + SCAL_INT(JM)
                     NREF3(IBIN,JBIN) = NREF3(IBIN,JBIN) + 1
                     RFAC3(JBIN,IBIN) = RFAC3(IBIN,JBIN)
                     RSUM3(JBIN,IBIN) = RSUM3(IBIN,JBIN)
                     NREF3(JBIN,IBIN) = NREF3(IBIN,JBIN)
                  END IF
               END IF
            END IF
 100     CONTINUE
C
C====== Accumulate values for R-factors
C
         DO 110 I = 1 , 3
            AVE(I) = SIGWI(I)/SIGW(I)
            DEL(I) = ABS(SCAL_INT(IM)-AVE(I))
            AV(I) = SIGI(I)/NI(I)
            DL(I) = ABS(SCAL_INT(IM)-AV(I))
            IF (RF(I)) THEN
               RN(I) = RN(I) + DEL(I)
               RD(I) = RD(I) + AVE(I)
               RGOF(I) = RGOF(I)+DEL(I)**2
               RRN(I) = RRN(I) + DL(I)
               RRD(I) = RRD(I) + AV(I)
               NRR(I) = NRR(I) + 1
            END IF
 110     CONTINUE
 120  CONTINUE
      GO TO 10
  999 FORMAT (3I4,2F9.3,I3,I4,2F9.4)
      END