**==aa0001.f  
C
C PROGRAM 'LAUESCALE'
C ===================
C
C FORTRAN 77
C
C JOHN W. CAMPBELL, APRIL 1988
C
C THIS PROGRAM IS A MODIFIED VERSION OF THE ORIGINAL 'LAUESCAL' USED TO
C NORMALISE A SET OF LAUE INTENSITY DATA (AS A FUNCTION OF WAVELENGTH)
C AGAINST A SET OF MONOCHROMATIC DATA AND TO OUTPUT A SET OF SCALED
C AND MERGED DATA. THIS UPDATED VERSION ALLOWS FOR THE CALCULATION OF
C AN ABSORPTION CORRECTION FUNCTION. THE WAVELENGTH NORMALISATION AND
C ABSORPTION CORRECTION STEPS MAY BE ITERATED.
C
C PARAMETER STATEMENT. (MXR = MAX. NO. OF LAUE REFLNS. ALLOWED)
C                      (MXL = MAX. NO. OF LAMBDA RANGES)
C
      PARAMETER(MXR=25000,MXL=8)
C
C IF VALUES ARE CHANGED, THEN CHANGE PARAMETER STATEMENTS
C THROUGHOUT THE PROGRAM
C
C
C
C FORMAT OF THE ARRAY 'IA' CONTAINING THE LAUE DATA
C
C 7 WORDS/REFLECTION
C
C        WORD 1  H,K,L,ICS   PACKED  2**18(H+128)+2**10(K+128)
C                                   +4*(L+128) + IC + 1
C        WORD 2  10000*LAMBDA
C        WORD 3  INTENSITY  (REAL)
C        WORD 4  SIG(I) + 2**16(ISYM)    ISYM=SYMMETRY NO.
C        WORD 5  10000*THETA             THETA IN RADIANS.
C        WORD 6  FILM X COORDINATE IN MM (REAL)
C        WORD 7  FILM Y COORDINATE IN MM (REAL)
C
C   ICS = 1 odd symmetry numbers, -1 even symmetry numbers if non-centro
C   symmetric data separated (otherwise 0)
C
C SPECIFICATION STATEMENTS
C ------------------------
C
      DIMENSION IA(MXR*7)
      CHARACTER*10 PGNAME
      LOGICAL NORM,STFIL,FILE1,FILE2,FILAGR,AFSCAL,SHLX1,SHLX2,SHLX3
      COMMON /FLAGS/NORM,STFIL,FILE1,FILE2,FILAGR,AFSCAL,SHLX1,SHLX2,
     +              SHLX3
      LOGICAL UPDATB,NRMCAL,ABSCAL
      COMMON /SCFLAG/UPDATB,NRMCAL,ABSCAL
      COMMON /CDATA/AMINL,AMAXL,RES1,RES2,TH1,TH2,BFAC,NRANGE,IBATCH,
     *              NITER,IABS,LASTIT,KBATCH
      COMMON /WGHTS/IWTS,AWTS,IWTL,AWTL,IWMEAN
      COMMON /RANGES/BMINL(MXL),BMAXL(MXL),BINS(MXL),NBIN(MXL),
     *               NFIT(MXL),NP(MXL),IFEX(MXL)
      COMMON /CURVE/CMIN(MXL),CMAX(MXL),NCOEFF(MXL),C(MXL,11)
      DOUBLE PRECISION PS
      COMMON /ABSURF/PS(121),NQORD,NRORD,QXMIN,QXMAX,RYMIN,RYMAX
      COMMON /BINS/NBTOT,SWL(100),SWIL(100),SWS(100),SWIS(100),W(100),
     *             R(100),ALAM(100),NVAL(100),AIS(100),AIL(100)
      COMMON /RFACTS/RNS,RDS,NRS,RN(3),RD(3),NRF(3),RNB(100),RDB(100),
     *               NRB(100)
      COMMON /RFACT2/THEMIN,THEMAX,RNTHE(10),RDTHE(10),NRTHE(10),
     *SISTHE(10),SILTHE(10)
      COMMON /RFACT3/STLMIN,STLMAX,RNSTL(10),RDSTL(10),NRSTL(10),
     *SISSTL(10),SILSTL(10)
      COMMON /RFACT4/FSQMIN,FSQMAX,RNFSQ(10),RDFSQ(10),NRFSQ(10),
     *SISFSQ(10),SILFSQ(10),ALOGI(10)
      COMMON /SCAL/SCALEI,SCAL2,SCALIN,SIGREJ
      COMMON /DIAGS/NDIAG1,NDIAG2
      COMMON /SYMM1/PGNAME
      COMMON /NWR/NWR
      COMMON /ABSCUT/ISMIN,ILMIN,DIAGR
      LOGICAL LAMOPT
      COMMON /LAMO/ LAMOPT, DLAM
      LOGICAL CENTYP
      COMMON /CENT/ CENTYP
      COMMON /GEFL/ INTYP, ISPAT
C
C INITIALISATIONS
C ---------------
C
      MAXREF = MXR
      NWR = 7
      IMAX = NWR*MXR
      SCALEI = 1000.0
      NRMCAL = .FALSE.
      ABSCAL = .FALSE.
      CENTYP = .FALSE.
      IFAIL = 0
      CALL CCPOPN(5,'DATA',5,1,LDUM,IFAIL)
      CALL CCPOPN(6,'PRINTER',6,1,LDUM,IFAIL)
      CALL CCPVRS(6,'LAUESCALE','15/10/92')
      CALL MTZINI
C
C READ IN CONTROL DATA
C --------------------
C
      CALL DATAIN
C
C READ IN LAUE DATA AND SORT ON HKL
C ---------------------------------
C
      IF (AFSCAL) THEN
         IFAIL = 0
         CALL CCPOPN(1,'LAUE',5,1,LDUM,IFAIL)
         CALL RDLAUE(IA,IMAX,MAXREF,NREF)
      ELSE
         IFAIL = 0
         CALL CCPOPN(1,'LAUE',5,4,19,IFAIL)
         CALL RDGE(1,INTYP,ISPAT,IA,IMAX,MAXREF,NREF)
         CLOSE (UNIT=1)
      END IF
      CALL REFSRT(IA,NREF,NWR,IMAX)
      IF (LAMOPT) CALL LAMSRT(IA,NREF)
C
C CALCULATE BIN REQUIREMENTS
C --------------------------
C
      NBTOT = 0
      DO 20 NN = 1 , NRANGE
         IF (BINS(NN).LE.0.0) GO TO 60
         IF (BINS(NN).GT.1.0) NBIN(NN) = NINT(BINS(NN))
         IF (BINS(NN).LE.1.0) NBIN(NN) = (BMAXL(NN)-BMINL(NN))/BINS(NN)
     +       + 0.999
         NBTOT = NBTOT + NBIN(NN)
         IF (NBIN(NN).LT.2) GO TO 70
         IF (NBTOT.GT.100) GO TO 80
         WRITE (6,6010) NBIN(NN) , NN
 20   CONTINUE
      NPP = 0
      DO 30 NN = 1 , NRANGE
         NP(NN) = NPP
         NPP = NPP + NBIN(NN)
 30   CONTINUE
      IF (.NOT.NORM) THEN
         IF (STFIL) CALL COLECT(IA,IMAX,NREF)
         NRMCAL = .TRUE.
         GO TO 50
      END IF
C
C ITERATION LOOP
C --------------
C
      DO 40 NIT = 1 , NITER
C
C COLLECT DATA IN BINS AND FIT NORMALISATION CURVE
C
         CALL COLECT(IA,IMAX,NREF)
         CALL CURFIT
         NRMCAL = .TRUE.
C
C OUTPUT STATISTICS FOR THE NORMALISATION
C
         CALL TABLE1
C
C OUTPUT SCALING STATISTICS AND CALCULATE B FACTOR
C
         IF (IABS.GT.0 .OR. UPDATB .OR. NIT.LT.NITER)
     +       CALL RSCALB(IA,IMAX,NREF)
C
C CALCULATE ASORPTION SURFACE IF REQUIRED
C
         IF (IABS.GT.0 .AND. NIT.LE.NITER-LASTIT) THEN
            CALL COLABS(IA,IMAX,NREF)
            ABSCAL = .TRUE.
            IF (UPDATB .OR. NIT.LT.NITER) CALL RSCALB(IA,IMAX,NREF)
         END IF
C
C END OF ITERATION LOOP
C ---------------------
C
 40   CONTINUE
C
C SCALE AND OUTPUT THE REFLECTIONS
C --------------------------------
C
 50   CALL OPSCAL(IA,IMAX,NREF)
C
C OUTPUT STATISTICS AFTER SCALING
C -------------------------------
C
      CALL TABLE2(IA,IMAX,NREF)
      STOP
C
C ERROR CONDITIONS
C ----------------
C
 60   WRITE (6,6020) BINS(NN) , NN
      STOP
 70   WRITE (6,6030)
      STOP
 80   WRITE (6,6040)
      STOP
C
C FORMAT STATEMENTS
C -----------------
C
 6010 FORMAT (//,' NUMBER OF LAMBDA BINS TO BE USED FOR NORMALISATION ='
     +        ,I4,' FOR RANGE',I2,/)
 6020 FORMAT (/,' **BINS PARAMETER MUST BE GREATER THAN 0, BINS =',
     +        F10.5,' FOR RANGE',I2,' **')
 6030 FORMAT (/,' **MINIMUM OF TWO BINS ALLOWED IN EACH RANGE**')
 6040 FORMAT (/,' **MAXIMUM OF 100 BINS ALLOWED**')
      END
**==datain.f
C
C SUBROUTINE 'DATAIN'
C ===================
C
C THIS SUBROUTINE READS IN THE CONTROL DATA
C
      SUBROUTINE DATAIN
C
C PARAMETERS
C
C NONE
C
      PARAMETER(MXR=25000,MXL=8)
C
C SPECIFICATION STATEMENTS
C ------------------------
C
      LOGICAL EOF
      INTEGER MTOK
      PARAMETER (MTOK=20)
      CHARACTER*4 KEY, CCVAL(MTOK)
      CHARACTER*400 LINE
      INTEGER NTOK, IBEG(MTOK), IEND(MTOK), ITYP(MTOK), IDEC(MTOK)
      REAL FVAL(MTOK)
      INTEGER IPOINT(6)
      CHARACTER*70 TITLE
      CHARACTER*80 ITEXT
      CHARACTER*30 PRGLAB(7)
      CHARACTER*1 CTYP(7)
      CHARACTER*30 LABL1(9)
      CHARACTER*1 CTYP1(9)
      CHARACTER*30 LABL2(7)
      CHARACTER*1 CTYP2(7)
      CHARACTER*30 LABO2(9)
      CHARACTER*1  CTYPO2(9)
      CHARACTER*30 LABO3(7)
      CHARACTER*1 CTYPO3(7)
      CHARACTER*10 SPGNAM
      CHARACTER*6 PGN
      CHARACTER*1 LTYPE
      CHARACTER*10 PGNAME
      CHARACTER*80 HSTR(2)
      INTEGER ISORT(5)
      REAL RSYMM(4,4,96)
      CHARACTER*15 INTT(2)
      DIMENSION LOOKUP(10),CELL(6)
      LOGICAL NORM,STFIL,FILE1,FILE2,FILAGR,AFSCAL,SHLX1,SHLX2,SHLX3
      COMMON /FLAGS/NORM,STFIL,FILE1,FILE2,FILAGR,AFSCAL,SHLX1,SHLX2,
     +              SHLX3
      LOGICAL UPDATB,NRMCAL,ABSCAL
      COMMON /SCFLAG/UPDATB,NRMCAL,ABSCAL
      COMMON /CDATA/AMINL,AMAXL,RES1,RES2,TH1,TH2,BFAC,NRANGE,IBATCH,
     *              NITER,IABS,LASTIT,KBATCH
      COMMON /WGHTS/IWTS,AWTS,IWTL,AWTL,IWMEAN
      COMMON /RANGES/BMINL(MXL),BMAXL(MXL),BINS(MXL),NBIN(MXL),
     *               NFIT(MXL),NP(MXL),IFEX(MXL)
      COMMON /CURVE/CMIN(MXL),CMAX(MXL),NCOEFF(MXL),C(MXL,11)
      DOUBLE PRECISION PS
      COMMON /ABSURF/PS(121),NQORD,NRORD,QXMIN,QXMAX,RYMIN,RYMAX
      COMMON /SCAL/SCALEI,SCAL2,SCALIN,SIGREJ
      COMMON /DIAGS/NDIAG1,NDIAG2
      COMMON /SYMM1/PGNAME
      COMMON /ABSCUT/ISMIN,ILMIN,DIAGR
      LOGICAL LAMOPT
      COMMON /LAMO/ LAMOPT, DLAM
      LOGICAL CENTYP
      COMMON /CENT/ CENTYP
      COMMON /GEFL/ INTYP, ISPAT
      REAL CVAL(11)
      DATA CELL/6*0.0/
      DATA PRGLAB/'H','K','L','FP','SIGFP','    ','    '/
      DATA CTYP/'H','H','H','F','Q',' ',' '/
      DATA LABL1/'H','K','L','FP','SIGFP','FL','SIGFL',' ',' '/
      DATA CTYP1/'H','H','H','F','Q','F','Q',' ',' '/
      DATA LABL2/'H','K','L','FL','SIGFL',' ',' '/
      DATA CTYP2/'H','H','H','F','Q',' ',' '/
      DATA LABO2/'H','K','L','LAMBDA','THETA','FP',
     +            'SIGFP','DANO','SIGDANO'/
      DATA CTYPO2/'H','H','H','R','R','F','Q','D','Q'/
      DATA LABO3/'H','K','L','M/ISYM','BATCH','I','SIGI'/
      DATA CTYPO3/'H','H','H','Y','B','J','Q'/
      DATA INTT/'BOX','PROFILE FITTED'/
      DATA NRMAX/MXL/
      DATA HSTR/'Created by the program LAUESCALE',' '/
C
C OUTPUT PRINTER HEADING
C ----------------------
C
      CALL HEADXX(6,'CONTROL DATA ITEMS$','$','0',2)
C
C DATA CARD 1     TITLE
C -----------
C
      READ (5,6010) TITLE
C
C DATA CARD 2      NORMALISE OR SCALE OR RSCALE ['GE', INTYP, ISPAT]
C -----------
C
      NORM = .FALSE.
      STFIL = .FALSE.
      AFSCAL = .TRUE.
      CALL FFREAD(5,ITEXT,80,6,.TRUE.)
      CALL TEXLIM(M1,M2)
      IF (M1.EQ.0) THEN
         WRITE (6,6420)
         STOP
      ELSE
         IF (ITEXT(M1:M2).EQ.'NORMALISE') THEN
            NORM = .TRUE.
            STFIL = .TRUE.
            GO TO 20
         END IF
         IF (ITEXT(M1:M2).EQ.'RSCALE') THEN
            STFIL = .TRUE.
            GO TO 20
         END IF
         IF (ITEXT(M1:M2).NE.'SCALE') THEN
            WRITE (6,6430) ITEXT(M1:M2)
            STOP
         END IF
      END IF
 20   WRITE (6,6020) ITEXT(M1:M2)
      CALL TEXLIM(M1,M2)
      IF (M1.GT.0) THEN
         IF (ITEXT(M1:M1).NE.'G' .AND. ITEXT(M1:M1).NE.'g') THEN
            WRITE (6,6560) ITEXT(M1:M2)
            STOP
         ELSE
            AFSCAL = .FALSE.
            CALL IVALUE(INTYP)
            IF (INTYP.LE.0 .OR. INTYP.GT.2) INTYP = 2
            WRITE (6,6530) INTT(INTYP)
            CALL IVALUE(ISPAT)
            IF (ISPAT.NE.1) ISPAT = 0
            IF (ISPAT.EQ.0) WRITE (6,6540)
            IF (ISPAT.EQ.1) WRITE (6,6550)
         END IF
      END IF
C
C DATA CARD 3     NSPGRP AMINL AMAXL RES1 RES2 TH1 TH2 BFAC NRANGE
C -----------
C
      CALL FFREAD(5,ITEXT,80,6,.TRUE.)
      CALL IVALUE(NSPGRP)
      CALL FVALUE(AMINL)
      CALL FVALUE(AMAXL)
      CALL FVALUE(RES1)
      CALL FVALUE(RES2)
      CALL FVALUE(TH1)
      CALL FVALUE(TH2)
      CALL FVALUE(BFAC)
      CALL IVALUE(NRANGE)
      CALL MSYMLB (10,NSPGRP,SPGNAM,PGN,NSYMP,NSYMM,RSYMM)
      CALL CENTRIC (NSYMM,RSYMM,0)
      LTYPE = SPGNAM(1:1)
      PGNAME = PGN(3:)
      CALL SYMMAT
      IF (AMINL.GE.AMAXL) THEN
         TEMP = AMINL
         AMINL = AMAXL
         AMAXL = TEMP
      END IF
      IF (AMINL.EQ.0.0 .AND. AMAXL.EQ.0.0) THEN
         WRITE (6,6040)
      ELSE IF (AMINL.LT.0.0 .OR. AMINL.EQ.AMAXL) THEN
         WRITE (6,6460)
         STOP
      ELSE
         WRITE (6,6030) AMINL , AMAXL
      END IF
      IF (RES1.EQ.0.0 .AND. RES2.EQ.0.0) THEN
         WRITE (6,6120)
      ELSE
         IF (RES1.LE.RES2) THEN
            TEMP = RES2
            RES2 = RES1
            RES1 = TEMP
         END IF
         WRITE (6,6110) RES1 , RES2
      END IF
      IF (TH1.EQ.0.0 .AND. TH2.EQ.0.0) THEN
         WRITE (6,6140)
      ELSE
         IF (TH1.GE.TH2) THEN
            TEMP = TH2
            TH2 = TH1
            TH1 = TEMP
         END IF
         WRITE (6,6130) TH1 , TH2
         TH1 = 0.017453*TH1
         TH2 = 0.017453*TH2
      END IF
      IF (NRANGE.LT.1 .OR. NRANGE.GT.NRMAX) THEN
         WRITE (6,6490) NRMAX
         STOP
      ELSE
         WRITE (6,6150) BFAC
         WRITE (6,6160) NRANGE
C
C DATA CARDS 4    BMINL BMAXL BINS NFIT IFEX [IORD P(0) .. P(IORD)]
C ------------
C
C  * NRANGE
C
         DO 60 NN = 1 , NRANGE
            CALL FFREAD(5,ITEXT,80,6,.TRUE.)
            CALL FVALUE(BMINL(NN))
            CALL FVALUE(BMAXL(NN))
            CALL FVALUE(BINS(NN))
            CALL IVALUE(NFIT(NN))
            CALL IVALUE(IFEX(NN))
            IF (IFEX(NN).LT.0 .OR. IFEX(NN).GT.3) GO TO 100
            IF (NFIT(NN).LT.0 .OR. NFIT(NN).GT.10) GO TO 80
            WRITE (6,6170) NN
            WRITE (6,6030) BMINL(NN) , BMAXL(NN)
            WRITE (6,6050) BINS(NN)
            WRITE (6,6060) NFIT(NN)
            IF (MOD(IFEX(NN),2).NE.0) WRITE (6,6290)
            IF (IFEX(NN).GT.1) WRITE (6,6300)
            IF (.NOT.NORM) THEN
               CMIN(NN) = BMINL(NN)
               CMAX(NN) = BMAXL(NN)
               CALL TEXLIM(M1,M2)
               IF (M1.EQ.0) GO TO 110
               CALL DECNUM(ITEXT,M1,M2,1,XX,IERR)
               IF (IERR.GT.0) GO TO 110
               IORD = NINT(XX)
               NCOEFF(NN) = IORD + 1
               NFIT(NN) = IORD
               DO 40 N = 1 , NCOEFF(NN)
                  CALL TEXLIM(M1,M2)
                  IF (M1.EQ.0) THEN
                     CALL GTVALS(5,ITEXT,80,IORD-N+2,CVAL,6,.TRUE.)
                     K = 0
                     DO 30 J = N , NCOEFF(NN)
                        K = K + 1
                        C(NN,J) = CVAL(K)
 30                  CONTINUE
                     GO TO 50
                  END IF
                  CALL DECNUM(ITEXT,M1,M2,2,C(NN,N),IERR)
                  IF (IERR.GT.0) GO TO 110
 40            CONTINUE
 50            WRITE (6,6070) IORD , (C(NN,N),N=1,NCOEFF(NN))
            END IF
 60      CONTINUE
         IF (NRANGE.NE.1 .OR. BMINL(1).NE.0.0 .OR. BMAXL(1).NE.0.0) THEN
            DO 70 NN = 1 , NRANGE
               IF (BMINL(NN).LT.0.0) GO TO 90
               IF (BMAXL(NN).LT.0.0) GO TO 90
               IF (BMAXL(NN).LE.BMINL(NN)) GO TO 90
               IF (NN.NE.NRANGE) THEN
                  IF (BMAXL(NN).GE.BMINL(NN+1)) GO TO 90
               END IF
 70         CONTINUE
         END IF
C
C DATA CARD 5     IFILE1 IFILE2 SCAL2 IWTS AWTS ITWL AWTL IWMEAN
C -----------     SCALIN SIGREJ NDIAG1 NDIAG2
C
         CALL FFREAD(5,ITEXT,80,6,.TRUE.)
         CALL IVALUE(IFILE1)
         CALL IVALUE(IFILE2)
         CALL FVALUE(SCAL2)
         CALL IVALUE(IWTS)
         CALL FVALUE(AWTS)
         CALL IVALUE(IWTL)
         CALL FVALUE(AWTL)
         CALL IVALUE(IWMEAN)
         CALL FVALUE(SCALIN)
         CALL FVALUE(SIGREJ)
         CALL IVALUE(NDIAG1)
         CALL IVALUE(NDIAG2)
         IF (SCAL2.EQ.0.0) SCAL2 = 20000.0
         IF (SCALIN.EQ.0.0) SCALIN = 1.0
         FILE1 = .FALSE.
         IF (IFILE1.GT.0) FILE1 = .TRUE.
         SHLX1 = .FALSE.
         SHLX2 = .FALSE.
         SHLX3 = .FALSE.
         IF (IFILE1.EQ.-1) SHLX1 = .TRUE.
         IF (IFILE1.EQ.-2) SHLX2 = .TRUE.
         IF (IFILE1.LE.-1000) THEN
            SHLX3 = .TRUE.
            KBATCH = - IFILE1 - 1000
         END IF
         IF (FILE1) WRITE (6,6270)
         IF (SHLX1) WRITE (6,6570)
         IF (SHLX2) WRITE (6,6580)
         IF (SHLX3) WRITE (6,6610)
         FILE2 = .FALSE.
         IF (IFILE2.LT.0) FILE2 = .TRUE.
         FILAGR = .FALSE.
         IF (IFILE2.GT.0) THEN
            FILAGR = .TRUE.
            IBATCH = IFILE2
         END IF
         IF (FILE2) WRITE (6,6080)
         IF (FILAGR) WRITE (6,6260)
         IF (FILAGR) WRITE (6,6280) SCAL2
         IFAIL = 0
         IF (NDIAG1.GT.0 .OR. NDIAG2.GT.0)
     +       CALL CCPOPN(9,'DIAGOUT',6,1,LDUM,IFAIL)
         IF (IWTS.EQ.0) WRITE (6,6180)
         IF (IWTS.EQ.1) WRITE (6,6190)
         IF (IWTS.GT.1) WRITE (6,6200) AWTS , SCALEI
         IF (IWTL.EQ.0) WRITE (6,6210)
         IF (IWTL.EQ.1) WRITE (6,6220)
         IF (IWTL.GT.1) WRITE (6,6230) AWTL
         WRITE (6,6240) SCALIN
         IF (SIGREJ.EQ.0.0) WRITE (6,6250)
         IF (SIGREJ.GT.0.0) WRITE (6,6250) SIGREJ
         IF (IWMEAN.EQ.0) WRITE (6,6090)
         IF (IWMEAN.NE.0) WRITE (6,6100)
C
C DATA CARD 6    NITER IUPDB IABS NQORD NRORD ISMIN ILMIN LASTIT DIAGR
C -----------    DLAM ICTYP
C
         CALL FFREAD(5,ITEXT,80,6,.TRUE.)
         CALL IVALUE(NITER)
         UPDATB = .FALSE.
         CALL IVALUE(IUPDB)
         IF (IUPDB.GT.0) UPDATB = .TRUE.
         CALL IVALUE(IABS)
         CALL IVALUE(NQORD)
         CALL IVALUE(NRORD)
         CALL IVALUE(ISMIN)
         CALL IVALUE(ILMIN)
         CALL IVALUE(LASTIT)
         IF (LASTIT.NE.0) LASTIT = 1
         CALL FVALUE(DIAGR)
         CALL FVALUE(DLAM)
         CALL IVALUE(ICTYP)
         
         IF (DIAGR.LE.0.0) DIAGR = 10.0
         IF (NQORD.GT.10 .OR. NRORD.GT.10) THEN
            WRITE (6,6480)
            STOP
         ELSE
            WRITE (6,6310) NITER
            IF (UPDATB) THEN
               WRITE (6,6320)
            ELSE
               WRITE (6,6330)
            END IF
            IF (IABS.GT.2) THEN
               WRITE (6,6520) IABS
               STOP
            ELSE
               IF (IABS.LE.0) THEN
                  IABS = 0
                  WRITE (6,6340)
               ELSE
                  WRITE (6,6350) IABS , NQORD , NRORD
               END IF
               IF (ISMIN.GT.0) THEN
                  WRITE (6,6360) ISMIN
               ELSE
                  WRITE (6,6370)
               END IF
               IF (ILMIN.GT.0) THEN
                  WRITE (6,6380) ILMIN
               ELSE
                  WRITE (6,6390)
               END IF
               IF (IABS.GT.0 .AND. LASTIT.EQ.1) WRITE (6,6400)
               IF (IABS.GT.0) WRITE (6,6410) DIAGR
               LAMOPT = .FALSE.
               IF (DLAM.GT.0.0) THEN
                  LAMOPT = .TRUE.
               ELSE
                  DLAM = 0.0
               END IF
               IF (LAMOPT) WRITE(6,6600) DLAM
               CENTYP = .FALSE.
               IF (ICTYP.GT.0) CENTYP = .TRUE.
               IF (CENTYP) WRITE(6,6620) 
C
C OPEN INPUT MTZ FILE IF USING 'NORMALISE' OR 'RSCALE'
C ----------------------------------------------------
C
               IF (STFIL) THEN
                  LINE = ' '
                  NTOK = MTOK
                  CALL PARSER (KEY,LINE,IBEG,IEND,ITYP,FVAL,CCVAL,
     +                         IDEC,NTOK,EOF,.FALSE.)
                  CALL CCPUPC(KEY)
                  IF (KEY.NE.'LABI') GO TO 80
                  NL = 5
                  IF (LAMOPT) THEN
                      NL = 6
                      PRGLAB(6) = 'LAM'
                      CTYP(6) = 'R'
                  ENDIF
                  IF (CENTYP) THEN
                     NL = NL + 1
                     PRGLAB(NL) = 'SIGN'
                     CTYP(NL) = 'I'
                  ENDIF
                  CALL LKYIN(1,PRGLAB,NL,NTOK,LINE,IBEG,IEND)
                  DO 75 I = 1 , NL
                     IPOINT(I) = -1
 75               CONTINUE
                  CALL LROPEN (1,'HKLIN',1,IFAIL)
                  IF (IFAIL.EQ.1) STOP
                  CALL LRASSN (1,PRGLAB,NL,IPOINT,CTYP)
                  CALL LRCELL (1,CELL)
               END IF
C
C OPEN OUTPUT MTZ(/SHELX) FILE(S)
C -------------------------------
C
               ISORT(1) = 1
               ISORT(2) = 2
               ISORT(3) = 3
               ISORT(4) = 0
               ISORT(5) = 0
               CALL LWOPEN (1,'HKLOUT')
               CALL LWTITL (1,TITLE,0)
               IF (STFIL) THEN
                  NLOUT = 7
                  IF (CENTYP) THEN
                     NLOUT = NLOUT + 1
                     LABL1(NLOUT) = 'SIGN'
                     CTYP1(NLOUT) = 'I'
                  ENDIF
                  IF (LAMOPT) THEN
                     NLOUT = NLOUT + 1
                     LABL1(NLOUT) = 'LAM'
                     CTYP1(NLOUT) = 'R'
                  ENDIF
                  CALL LWCLAB (1,LABL1,NLOUT,CTYP1,0)
               ELSE
                  NLOUT = 5
                  IF (CENTYP) THEN
                     NLOUT = NLOUT + 1
                     LABL2(NLOUT) = 'SIGN'
                     CTYP2(NLOUT) = 'I'
                  ENDIF
                  IF (LAMOPT) THEN
                     NLOUT = NLOUT + 1
                     LABL2(NLOUT) = 'LAM'
                     CTYP2(NLOUT) ='R'
                  ENDIF
                  CALL LWCLAB (1,LABL2,NLOUT,CTYP2,0)
               END IF
               CALL LWHIST (1,HSTR,1)
               CALL LWCELL (1,CELL)
               CALL LWSYMM (1,NSYMM,NSYMP,RSYMM,LTYPE,NSPGRP,
     +                      SPGNAM,PGNAME)        
               CALL LWSORT (1,ISORT)
               IF (FILE2) THEN
                  ISORT(1) = 1
                  ISORT(2) = 2
                  ISORT(3) = 3
                  ISORT(4) = 0
                  ISORT(5) = 0
                  CALL LWOPEN (2,'HKLOUT2')
                  CALL LWTITL (2,TITLE,0)
                  CALL LWCLAB (2,LABO2,9,CTYPO2,0)
                  CALL LWHIST (2,HSTR,1)
                  CALL LWCELL (2,CELL)
                  CALL LWSYMM (2,NSYMM,NSYMP,RSYMM,LTYPE,NSPGRP,
     +                      SPGNAM,PGNAME)        
                  CALL LWSORT (2,ISORT)
               END IF
               IF (FILAGR) THEN
                  ISORT(1) = 1
                  ISORT(2) = 2
                  ISORT(3) = 3
                  ISORT(4) = 0
                  ISORT(5) = 0
                  CALL LWOPEN (2,'HKLOUT2')
                  CALL LWTITL (2,TITLE,0)
                  CALL LWCLAB (2,LABO3,7,CTYPO3,0)
                  CALL LWHIST (2,HSTR,1)
                  CALL LWCELL (2,CELL)
                  CALL LWSYMM (2,NSYMM,NSYMP,RSYMM,LTYPE,NSPGRP,
     +                      SPGNAM,PGNAME)        
                  CALL LWSORT (2,ISORT)
                  CALL LWBTIT (2,IBATCH,TITLE)
               END IF
               IF (SHLX1 .OR. SHLX2 .OR. SHLX3) THEN
                  IFAIL = 0
                  CALL CCPOPN(2,'SHELX',4,1,LDUM,IFAIL)
               END IF
               RETURN
            END IF
         END IF
      END IF
 80   WRITE (6,6470)
      STOP
 90   WRITE (6,6500) (BMINL(NN),BMAXL(NN),NN=1,NRANGE)
      STOP
 100  WRITE (6,6510) NN
      STOP
 110  WRITE (6,6590)
      STOP
C
C FORMAT STATEMENTS
C -----------------
C
 6010 FORMAT (A80)
 6020 FORMAT (//,' OPTION SELECTED: ',A)
 6030 FORMAT (//,' MINIMUM LAMBDA =',F5.3,'A  MAXIMUM LAMBDA =',F5.3,
     +        'A')
 6040 FORMAT (//,' LAMBDA LIMITS WILL BE TAKEN FROM INPUT LAUE FILE')
 6050 FORMAT (//,' BINS PARAMETER =',F7.3)
 6060 FORMAT (//,' ORDER FOR CURVE FITTING =',I3)
 6070 FORMAT (//,' POLYNOMIAL COEFFICIENTS FOR SCALING OPTIONS  P 0',
     +        ' TO  P',I2,//,(' ',6E20.8))
 6080 FORMAT (//,' WRITE FILE OF UNMERGED DATA WITH ANOMALOUS',
     +        ' DIFFERENCES')
 6090 FORMAT (//,' USE WEIGHTS CALCULATED FROM SIGMA VALUES',
     +        ' WHEN CALCULATING LAUE MEANS')
 6100 FORMAT (//,' USE UNIT WEIGHTS WHEN CALCULATING LAUE MEANS')
 6110 FORMAT (//,' EXCLUDE REFLECTIONS OUTSIDE RESOLUTION LIMITS OF',
     +        F8.1,'A  TO',F8.1,'A')
 6120 FORMAT (//,' NO RESOLUTION LIMITS APPLIED')
 6130 FORMAT (//,' EXCLUDE REFLECTIONS OUTSIDE THETA LIMITS OF',F8.1,
     +        ' DEG  TO',F8.1,' DEG')
 6140 FORMAT (//,' NO THETA LIMITS APPLIED')
 6150 FORMAT (//,' TEMPERATURE FACTOR TO BE APPLIED TO LAUE',
     +        ' REFLECTIONS ON INPUT =',F8.2)
 6160 FORMAT (//,' NUMBER OF LAMBDA RANGES TO BE USED IN SCALING =',I2)
 6170 FORMAT (////,' CURVE FITTING RANGE',I2)
 6180 FORMAT (//,' USE SIGMAS IN WEIGHTING WHEN CALCULATING MEAN',
     +        ' STANDARD INTENSITY VALUES FOR LAMBDA BINS')
 6190 FORMAT (//,' USE UNIT WEIGHTS WHEN CALCULATING MEAN STANDARD',
     +        ' INTENSITY VALUES FOR LAMBDA BINS')
 6200 FORMAT (//,' USE FOLLOWING WEIGHTING SCHEME WHEN CALCULATING',
     +        ' STANDARD INTENSITIES FOR LAMBDA BINS:  W = 1.0 + ',F8.5,
     +        '*IS/(',F7.1,')')
 6210 FORMAT (//,' USE SIGMAS IN WEIGHTING WHEN CALCULATING MEAN',
     +        ' LAUE INTENSITY VALUES FOR LAMBDA BINS')
 6220 FORMAT (//,' USE UNIT WEIGHTS WHEN CALCULATING MEAN LAUE',
     +        ' INTENSITY VALUES FOR LAMBDA BINS')
 6230 FORMAT (//,' USE FOLLOWING WEIGHTING SCHEME WHEN CALCULATING',
     +        ' LAUE INTENSITIES FOR LAMBDA BINS:  W = 1.0 + ',F8.5,
     +        '*IL')
 6240 FORMAT (//,' INTENSITIES AND SIGMA(INTENSITIES) MULTIPLIED BY',
     +        F7.3,' ON INPUT')
 6250 FORMAT (//,' REJECT REFLECTIONS WITH SIGMA VALUES OF ZERO',:,
     +        ' OR WITH I LESS THAN',F4.1,' * SIG(I)')
 6260 FORMAT (//,' WRITE FILE OF UNMERGED DATA FOR ROTAVATA/AGROVATA',
     +        ' INPUT')
 6270 FORMAT (//,' WRITE FILE OF MERGED AND SCALED LAUE DATA')
 6280 FORMAT (/,' OUTPUT INTENSITIES DIVIDED BY',F12.1)
 6290 FORMAT (/' ADD EXTRA EXTRAPOLATED POINT AT LOW WAVELENGTH',
     +        ' END OF RANGE BEFORE CURVE FITTING')
 6300 FORMAT (/' ADD EXTRA EXTRAPOLATED POINT AT HIGH WAVELENGTH',
     +        ' END OF RANGE BEFORE CURVE FITTING')
 6310 FORMAT (//,' NUMBER OF ITERATIONS REQUIRED =',I3)
 6320 FORMAT (/,' UPDATE TEMPERATURE FACTOR USING CALCULATED CHANGE')
 6330 FORMAT (/,' DO NOT UPDATE TEMPERATURE FACTOR AT THE END',
     +        ' OF THE ITERATION')
 6340 FORMAT (/,' NO ABSORPTION CORRECTION TO BE APPLIED')
 6350 FORMAT (/,' APPLY ABSORPTION CORRECTION TYPE',I2,//,
     +        ' POLYNOMIAL ORDERS ARE',I3,' FOR Q AND',I3,' FOR R')
 6360 FORMAT (/,
     +' REFLECTIONS EXCLUDED FROM ABSORPTION SURFACE ''SURFACE FOR STAND
     +        ARD INTENSITIES (UNSCALED) LESS THAN',I12)
 6370 FORMAT (/,
     +' NO LOW INTENSITY CUTOFF APPLIED TO STANDARD REFLECTION'' INTENSI
     +        TIES WHEN CALCULATING ABSORPTION SURFACE')
 6380 FORMAT (/,
     +' REFLECTIONS EXCLUDED FROM ABSORPTION SURFACE ''SURFACE FOR LAUEI
     +        NTENSITIES (UNSCALED) LESS THAN',I12)
 6390 FORMAT (/,
     +' NO LOW INTENSITY CUTOFF APPLIED TO LAUE REFLECTION'' INTENSITIES
     +        WHEN CALCULATING ABSORPTION SURFACE')
 6400 FORMAT (/,' OMIT ABSORPTION CORRECTION ON FINAL ITERATION')
 6410 FORMAT (/,' PRINT ABSORPTION SURFACE DIAGNOSTIC DETAILS FOR',
     +        ' REFLECTIONS WITH IL/IS OR IS/IL >= ',F8.2)
C
 6420 FORMAT (/,' **NORMALISE,SCALE OR RSCALE MUST BE SPECIFIED**')
 6430 FORMAT (/,' **INVALID CODE**: ',A)
 6440 FORMAT (/,' **POINT GROUP NAME MUST BE GIVEN**')
 6450 FORMAT (/,' **POINT GROUP NAME MUST HAVE MAX. OF 4 CHARACTERS**: '
     +        ,A)
 6460 FORMAT (/,' **INVALID LAMBDA LIMITS**')
 6470 FORMAT (/,' **ORDER FOR CURVE FITTING MUST BE 0-10**')
 6480 FORMAT (/,' **MAXIMUM ORDER ALLOWED IS 10 FOR POLYNOMIAL',
     +        ' ABSORPTION SURFACE**')
 6490 FORMAT (/,' **NUMBER OF LAMBDA RANGES MUST BE 1 TO',I2,' **')
 6500 FORMAT (/,' **LAMBDA RANGE LIMITS MUST BE IN INCREASING POSITIVE',
     +        ' VALUES OF LAMBDA**:',12F7.3)
 6510 FORMAT (/,' **INVALID VALUE OF IFEX FOR RANGE',I2,' **')
 6520 FORMAT (/,' **ABSORPTION TYPE',I3,' NOT AVAILABLE**')
C
 6530 FORMAT (/,' ''A'' FILM DATA INPUT FROM GE1 FILE - INTEGRATION',
     +        ' TYPE: ',A)
 6540 FORMAT (/,' SPATIAL OVERLAP DATA NOT INCLUDED')
 6550 FORMAT (/,' SPATIAL OVERLAP DATA INCLUDED')
 6560 FORMAT (/,' **UNKNOWN OPTION: ',A,' **')
 6570 FORMAT (/,' OUTPUT RELECTION FILE WILL CONTAIN MERGED REFLECTION',
     +        ' DATA (F AND SIGF) IN SHELX FORMAT')
 6580 FORMAT (/,' OUTPUT REFLECTION FILE WILL CONTAIN UNMERGED',
     +        ' REFLECTION DATA (F, SIGF, LAMBDA) IN SHELX FORMAT')
 6590 FORMAT (/,' **INVALID CURVE FITTING PARAMETER SPECIFICATION**')
 6600 FORMAT (/,' USING WAVELENGTH TAGGED DATA WITH DLAM CHECK OF',
     +        F7.4)
 6610 FORMAT (/,' OUTPUT REFLECTION FILE WILL CONTAIN UNMERGED',
     +        ' REFLECTION DATA (F**2 SIG(F**2) BATCH LAMDA) IN',
     +        ' SHELX FORMAT')
 6620 FORMAT (/,' USING SEPARATED I+ AND I- REFERENCE DATA')
      END
**==gecell.f  
C
C SUBROUTINE 'GECELL'
C ===================
C
C GET REAL CELL PARAMETERS FROM A 'GE' FILE
C
      SUBROUTINE GECELL(LUN,CP)
C
C PARAMETERS
C
C         LUN (I)   STREAM NUMBER FOR 'GE' FILE
C          CP (O)   6 WORD ARRAY RETURNING CELL PARAMETERS
C
C
C ***WARNING - THIS ROUTINE CONTAINS NON-STANDARD FORTRAN77 CODE***
C
C SPECIFICATION STATEMENTS
C
      REAL CP(6),BUFF(9),BUFF2(3)
      REAL PARAMS(6)
      INTEGER XALNO
      INTEGER ISYLAT(2)
      CHARACTER*40 BUF
C
C CODE
C ----
C
      READ (LUN,REC=3) WFS , BUFF , CP , ISL
      CALL CCPI2I(ISYLAT,ISL,2,.TRUE.,.FALSE.)
      XALNO = ISYLAT(1)
      READ (LUN,REC=4) BUF , BUFF2 , WAVE
      IF (XALNO.EQ.1) THEN
         CALL RECCEL(PARAMS,CP,WAVE)
      ELSE IF (XALNO.EQ.2) THEN
         CP(5) = CP(4)
         CP(4) = 90
         CP(6) = 90
         CALL RECCEL(PARAMS,CP,WAVE)
      ELSE IF (XALNO.EQ.3) THEN
         CP(4) = 90
         CP(5) = 90
         CP(6) = 90
         CALL RECCEL(PARAMS,CP,WAVE)
      ELSE IF (XALNO.EQ.4) THEN
         CP(3) = CP(2)
         CP(2) = CP(1)
         CP(4) = 90
         CP(5) = 90
         CP(6) = 90
         CALL RECCEL(PARAMS,CP,WAVE)
      ELSE IF (XALNO.EQ.5) THEN
         CP(3) = CP(2)
         CP(2) = CP(1)
         CP(4) = 90
         CP(5) = 90
         CP(6) = 60
         CALL RECCEL(PARAMS,CP,WAVE)
      ELSE IF (XALNO.EQ.6) THEN
         CP(4) = CP(2)
         CP(5) = CP(2)
         CP(6) = CP(2)
         CP(2) = CP(1)
         CP(3) = CP(1)
         CALL RECCEL(PARAMS,CP,WAVE)
      ELSE IF (XALNO.EQ.7) THEN
         CP(2) = CP(1)
         CP(3) = CP(1)
         CP(4) = 90
         CP(5) = 90
         CP(6) = 90
         CALL RECCEL(PARAMS,CP,WAVE)
      ELSE
         WRITE (6,6010)
         STOP
      END IF
      DO 20 I = 1 , 6
         CP(I) = PARAMS(I)
 20   CONTINUE
      RETURN
C
C FORMAT STATEMENTS
C -----------------
C
 6010 FORMAT (/,' **ERROR IN GETTING CELL PARAMETERS FROM GE FILE**')
      END
**==getint.f  
C
C SUBROUTINE 'GETINT'
C ===================
C
C READ A REFLECTION FROM A LAUE 'GE' FILE
C
      SUBROUTINE GETINT(IUN,IREC,IFILM,INTYP,ISPAT,IH,IK,IL,ALBDA,
     *                  II,ISIG,X,Y,REJ)
C
C PARAMETERS
C
C         IUN (I)   UNIT NUMBER OF 'GE' FILE
C        IREC (I)   RECORD NUMBER OF REFLECTION
C       IFILM (I)   NO. OF THE FILM FOR WHICH THE INTENSITY IS TO
C                   BE RETURNED
C       INTYP (I)   INTEGRATION TYPE 1=BOX, 2=PROFILE FITTED
C       ISPAT (I)   =1 INCLUDE SPATIAL OVERLAP DATA, =0 DO NOT
C          IH (O)   INDICES
C          IK (O)   "
C          IL (O)   "
C       ALBDA (O)   WAVELENGTH
C          II (O)   INTENSITY
C        ISIG (O)   SIGMA(INTENSITY)
C         REJ (O)   REJECT FLAG = .TRUE. REJECT, = .FALSE. DO NOT
C
C
C SPECIFICATION STATEMENTS
C ------------------------
C
      LOGICAL REJ
      INTEGER MHKL(4), MULTF(1)
      INTEGER IBUFF(12)
      INTEGER INTB(6),ISIGB(6),INTP(6),ISIGP(6)
C
C READ AND INTERPRET RECORD
C -------------------------
C
      REJ = .FALSE.
      READ (IUN,REC=IREC) X , Y , IBUFF , IMULT , ALBDA
C
C REJECT MULTIPLE RELECTIONS, NODALS, SPATIAL OVERLAPS
C
      CALL CCPBYI(MULTF,IMULT,1)
      CALL CCPBYI(MHKL,IMULT,-4)
      IF (NBITST(MULTF(1),0,5).GT.1) REJ = .TRUE.
      IF (ISPAT.EQ.0 .AND. NBITST(MULTF(1),6,1).NE.0) REJ = .TRUE.
      IH = MHKL(2)
      IK = MHKL(3)
      IL = MHKL(4)
      IF (INTYP.EQ.1) THEN
         CALL CCPI2I(INTB,IBUFF(1),6,.TRUE.,.FALSE.)
         CALL CCPI2I(ISIGB,IBUFF(4),6,.TRUE.,.FALSE.)
         II = INTB(IFILM)
         ISIG = ISIGB(IFILM)
      ELSE
         CALL CCPI2I(INTP,IBUFF(7),6,.TRUE.,.FALSE.)
         CALL CCPI2I(ISIGP,IBUFF(10),6,.TRUE.,.FALSE.)
         II = INTP(IFILM)
         ISIG = ISIGP(IFILM)
      END IF
      IF (ISIG.EQ.-9999) REJ = .TRUE.
      RETURN
      END
**==rdge.f  
C
C SUBROUTINE 'RDGE'
C =================
C
C THIS SUBROUTINE READS IN THE LAUE DATA FROM A .GE1 FILE  AND STORES IT IN THE
C ARRAY 'IA'
C
      SUBROUTINE RDGE(LUN,INTYP,ISPAT,IA,IMAX,MAXREF,NREF)
C
C PARAMETERS
C
C        LUN (I)   LOGICAL UNIT NUMBER FOR READING .GE1 FILE
C      INTYP (I)   INTEGRATION TYPE 1=BOX, 2=PROFILE-FITTED
C      ISPAT (I)   SPATIALS FLAG 1=INCLUDE SPATIALS, 0=DO NOT
C         IA (I)   ARRAY TO HOLD LAUE DATA (NWR WORDS/REFLN)
C       IMAX (I)   NO. OF WORDS IN 'IA'
C     MAXREF (I)   MAXIMUM NO. OF REFLNS ALLOWED
C       NREF (I)   NUMBER OF REFLECTIONS READ AND ACCEPTED.
C
      PARAMETER(MXR=25000,MXL=8)
C
C SPECIFICATION STATEMENTS
C ------------------------
C
      DIMENSION IA(IMAX)
      LOGICAL NORM,STFIL,FILE1,FILE2,FILAGR,LAMLIM,REJ,AFSCAL
      LOGICAL SHLX1,SHLX2,SHLX3
      COMMON /FLAGS/NORM,STFIL,FILE1,FILE2,FILAGR,AFSCAL,SHLX1,SHLX2,
     +              SHLX3
      COMMON /CDATA/AMINL,AMAXL,RES1,RES2,TH1,TH2,BFAC,NRANGE,IBATCH,
     *              NITER,IABS,LASTIT,KBATCH
      COMMON /WGHTS/IWTS,AWTS,IWTL,AWTL,IWMEAN
      COMMON /RANGES/BMINL(MXL),BMAXL(MXL),BINS(MXL),NBIN(MXL),
     *               NFIT(MXL),NP(MXL),IFEX(MXL)
      COMMON /RFACT2/THEMIN,THEMAX,RNTHE(10),RDTHE(10),NRTHE(10),
     *SISTHE(10),SILTHE(10)
      COMMON /RFACT3/STLMIN,STLMAX,RNSTL(10),RDSTL(10),NRSTL(10),
     *SISSTL(10),SILSTL(10)
      COMMON /RFACT4/FSQMIN,FSQMAX,RNFSQ(10),RDFSQ(10),NRFSQ(10),
     *SISFSQ(10),SILFSQ(10),ALOGI(10)
      COMMON /SCAL/SCALEI,SCAL2,SCALIN,SIGREJ
      COMMON /CURVE/CMIN(MXL),CMAX(MXL),NCOEFF(MXL),C(MXL,11)
      DOUBLE PRECISION PS
      COMMON /ABSURF/PS(121),NQORD,NRORD,QXMIN,QXMAX,RYMIN,RYMAX
      CHARACTER*10 PGNAME
      COMMON /SYMM1/PGNAME
      COMMON /NWR/NWR
      LOGICAL CENTYP
      COMMON /CENT/CENTYP
      DOUBLE PRECISION Q,R
      REAL CELL(6)
      DIMENSION IHKL(3)
      EQUIVALENCE (IAA,AAA)
      DATA NTWO18/262144/,NTWO10/1024/,NTWO16/65536/
C
C INITIALISATIONS
C ---------------
C
      LAMLIM = .TRUE.
      IF (NRANGE.EQ.1 .AND. BMINL(1).EQ.0.0 .AND. BMAXL(1).EQ.0.0)
     +    LAMLIM = .FALSE.
      NREF = 0
      NEXCL = 0
      NREJ = 0
      NNEG = 0
      NSIGRJ = 0
      NREJL = 0
      NREJRS = 0
      NREJTH = 0
      QMINL = 10000.0
      QMAXL = 0.0
      THEMIN = 10E20
      THEMAX = 0.0
      STLMIN = 10E20
      STLMAX = 0.0
      QXMIN = 10E20
      QXMAX = -10E20
      RYMIN = 10E20
      RYMAX = -10E20
C
C READ IN LAUE DATA, CONVERTING INDICES TO UNIQUE SET
C ---------------------------------------------------
C
      AKP = 0.105
      ATAU = 1.0
      SCIN = 100.0*SCALIN
      CALL GECELL(LUN,CELL)
      CALL SETRSL(CELL(1),CELL(2),CELL(3),CELL(4),CELL(5),CELL(6))
      READ (LUN,REC=10) IPACK , NSPOT
      IF (NSPOT.LE.0) READ (LUN,REC=9) JJJ , KKK , NSPOT
      NREC = NSPOT + 10
      IFILM = 1
      IREC = 10
 20   IREC = IREC + 1
      IF (IREC.GT.NREC) GO TO 50
      CALL GETINT(LUN,IREC,IFILM,INTYP,ISPAT,IH,IK,IL,ALBDA,II,ISIG,X,Y,
     +            REJ)
      IF (ALBDA.GT.QMAXL) QMAXL = ALBDA
      IF (ALBDA.LT.QMINL) QMINL = ALBDA
      If (REJ) THEN
         NEXCL = NEXCL + 1
         NREJ = NREJ + 1
         GO TO 20
      END IF
      STHL2 = STHLSQ(IH,IK,IL)
      STH = ALBDA*SQRT(STHL2)
      THETA = ASIN(STH)
      IF (SIGREJ.GT.0.0) THEN
         IF (II.LT.SIGREJ*ISIG) THEN
            NREJ = NREJ + 1
            NSIGRJ = NSIGRJ + 1
            GO TO 20
         END IF
      END IF
      IF (AMAXL.GT.0) THEN
         IF (ALBDA.LT.AMINL .OR. ALBDA.GT.AMAXL) THEN
            NREJ = NREJ + 1
            NREJL = NREJL + 1
            GO TO 20
         END IF
      END IF
      IF (LAMLIM) THEN
         DO 30 NN = 1 , NRANGE
            IF (ALBDA.GE.BMINL(NN) .AND. ALBDA.LE.BMAXL(NN)) GO TO 40
 30      CONTINUE
         NREJ = NREJ + 1
         NREJL = NREJL + 1
         GO TO 20
      END IF
 40   IF (RES2.GT.0.0) THEN
         RES = ALBDA/(2*SIN(THETA))
         IF (RES.GT.RES1 .OR. RES.LT.RES2) THEN
            NREJ = NREJ + 1
            NREJRS = NREJRS + 1
            GO TO 20
         END IF
      END IF
      IF (TH2.GT.0.0) THEN
         IF (THETA.LT.TH1 .OR. THETA.GT.TH2) THEN
            NREJ = NREJ + 1
            NREJTH = NREJTH + 1
            GO TO 20
         END IF
      END IF
      IF (II.LE.0) THEN
         NREJ = NREJ + 1
         NNEG = NNEG + 1
         GO TO 20
      END IF
      NREF = NREF + 1
      IF (NREF.GT.MAXREF) THEN
         WRITE (6,6120) MAXREF
         STOP
      ELSE
         IF (IH.LT.-127 .OR. IH.GT.127) GO TO 60
         IF (IK.LT.-127 .OR. IK.GT.127) GO TO 60
         IF (IL.LT.-127 .OR. IL.GT.127) GO TO 60
         IF (PGNAME(1:1).EQ.'3' .OR. PGNAME(1:1).EQ.'6'
     +       .OR. PGNAME(1:1).EQ.'4') THEN
            CALL TURN3(IH,IK,IL,ISYM)
         ELSE
            CALL TURNIP(IH,IK,IL,ISYM)
         END IF
         ICS = 0
         IF (CENTYP) THEN
            ICS = 1
            IF (MOD(ISYM,2).EQ.0) ICS = -1
            IHKL(1) = IH
            IHKL(2) = IK
            IHKL(3) = IL
            CALL CENTR(IHKL,IIC)
            IF (IIC.EQ.1) ICS = 0
         ENDIF
         I = NWR*(NREF-1)
         IA(I+1) = NTWO18*(IH+128) + NTWO10*(IK+128) + 
     +     4*(IL+128) + ICS + 1
         IA(I+2) = 10000*ALBDA
         STL = SIN(THETA)/ALBDA
         STL = STL*STL
C
C APPLY LORENTZ, POLARISATION AND OBLIQUITY CORRECTIONS
C
         ALOR = 1.0/(STH*STH)
C
         IF (Y.EQ.0.0) THEN
            APOL = 1.0
         ELSE
            PSI = ATAN(X/Y)
            C2TH2 = COS(2.0*THETA)**2
            P0 = (1+C2TH2)/2.0
            ALPHA = COS(2.0*PSI)*(SIN(2.0*THETA)**2)
            ALPHA = ALPHA/(1.0+C2TH2)
            APOL = P0*(1.0-ATAU*ALPHA)
         END IF
C
         C2TH = COS(2.0*THETA)
         OBLIQ = EXP(AKP*ATAU/C2TH)
C
         AAA = II
         AAA = AAA*SCIN*OBLIQ/(ALOR*APOL)
         ISIG = NINT(ISIG*SCIN*OBLIQ/(ALOR*APOL))
         IF (ISIG.LE.0) ISIG = 1
         IA(I+3) = IAA
         IF (ISIG.GT.32767) THEN
            WRITE (6,6140) IH , IK , IL , ISIG
            STOP
         ELSE
            IA(I+4) = NTWO16*ISYM + ISIG
            IA(I+5) = 10000*THETA
            IF (THETA.GT.THEMAX) THEMAX = THETA
            IF (THETA.LT.THEMIN) THEMIN = THETA
            IF (STL.GT.STLMAX) STLMAX = STL
            IF (STL.LT.STLMIN) STLMIN = STL
            AAA = X
            IA(I+6) = IAA
            AAA = Y
            IA(I+7) = IAA
            IF (IABS.GT.0) THEN
               CALL QRVAL(IA(I+1),IABS,Q,R)
               IF (Q.GT.QXMAX) QXMAX = Q
               IF (Q.LT.QXMIN) QXMIN = Q
               IF (R.GT.RYMAX) RYMAX = R
               IF (R.LT.RYMIN) RYMIN = R
            END IF
            GO TO 20
         END IF
      END IF
C
C END OF REFLECTIONS REACHED
C --------------------------
C
 50   IF (NREF.EQ.0) THEN
         WRITE (6,6150)
         STOP
      ELSE
         WRITE (6,6010) NSPOT
         WRITE (6,6020) NREF
         WRITE (6,6030) NREJ
         WRITE (6,6040) NEXCL
         WRITE (6,6050) NNEG
         WRITE (6,6060) SIGREJ , NSIGRJ
         WRITE (6,6070) NREJL
         WRITE (6,6080) NREJTH
         WRITE (6,6090) NREJRS
         WRITE (6,6100) QMINL , QMAXL
         IF (NORM) WRITE (6,6110) SCALEI
         AMINL = QMINL
         AMAXL = QMAXL
         IF (NRANGE.EQ.1 .AND. BMINL(1).EQ.0.0 .AND. BMAXL(1).EQ.0.0)
     +       THEN
            BMINL(1) = QMINL
            BMAXL(1) = QMAXL
            IF (.NOT.NORM) THEN
               CMIN(1) = QMINL
               CMAX(1) = QMAXL
            END IF
         END IF
         RETURN
      END IF
 60   WRITE (6,6130) IH , IK , IL
      STOP
C
C FORMAT STATEMENTS
C -----------------
C
 6010 FORMAT (//,' NUMBER OF RELECTIONS IN .GE1 FILE          =',I6)
 6020 FORMAT (' NUMBER OF REFLECTIONS ACCEPTED             =',I6)
 6030 FORMAT (' NUMBER REJECTED (TOTAL)                    =',I6)
 6040 FORMAT (' NUMBER REJECTED (NOT SELECTED/MEASURED)    =',I6)
 6050 FORMAT (' NUMBER REJECTED (INTENSITY <= 0)           =',I6)
 6060 FORMAT (' NUMBER REJECTED (LESS THAN',F6.2,' SIGMA)    =',I6)
 6070 FORMAT (' NUMBER REJECTED (OUTSIDE LAMBDA LIMITS)    =',I6)
 6080 FORMAT (' NUMBER REJECTED (OUTSIDE THETA LIMITS)     =',I6)
 6090 FORMAT (' NUMBER REJECTED (OUTSIDE RESOLUTION LIMITS)=',I6)
 6100 FORMAT (//,' LAMBDA LIMITS FROM LAUE FILE ARE',F7.3,',',F7.3,
     +        ' ANGSTROMS')
 6110 FORMAT (//,' INPUT STANDARD INTENSITIES ARE DIVIDED BY',F7.1,
     +        ' DURING SCALING PROCESS')
C
 6120 FORMAT (/,' **NUMBER OF REFLECTIONS MAY NOT EXCEED',I6,' **')
 6130 FORMAT (/,' **INDICES OUT OF RANGE -127 TO 127, HKL=',3I4,' **')
 6140 FORMAT (/,' **SIGMA OUT OF RANGE 0-32767 FOR REFLECTION HKL=',3I4,
     +        '  SIGMA =',I8,' **')
 6150 FORMAT (/,' **NO LAUE REFLECTIONS ACCEPTED**')
      END
**==rdlaue.f  
C
C SUBROUTINE 'RDLAUE'
C ===================
C
C THIS SUBROUTINE READS IN THE LAUE DATA AND STORES IT IN THE
C ARRAY 'IA'
C
      SUBROUTINE RDLAUE(IA,IMAX,MAXREF,NREF)
C
C PARAMETERS
C
C         IA (I)   ARRAY TO HOLD LAUE DATA (NWR WORDS/REFLN)
C       IMAX (I)   NO. OF WORDS IN 'IA'
C     MAXREF (I)   MAXIMUM NO. OF REFLNS ALLOWED
C       NREF (I)   NUMBER OF REFLECTIONS READ AND ACCEPTED.
C
      PARAMETER(MXR=25000,MXL=8)
C
C SPECIFICATION STATEMENTS
C ------------------------
C
      DIMENSION IA(IMAX)
      DIMENSION IHKL(3)
      LOGICAL NORM,STFIL,FILE1,FILE2,FILAGR,LAMLIM,AFSCAL,SHLX1,SHLX2,
     +        SHLX3
      COMMON /FLAGS/NORM,STFIL,FILE1,FILE2,FILAGR,AFSCAL,SHLX1,SHLX2,
     +              SHLX3
      COMMON /CDATA/AMINL,AMAXL,RES1,RES2,TH1,TH2,BFAC,NRANGE,IBATCH,
     *              NITER,IABS,LASTIT,KBATCH
      COMMON /WGHTS/IWTS,AWTS,IWTL,AWTL,IWMEAN
      COMMON /RANGES/BMINL(MXL),BMAXL(MXL),BINS(MXL),NBIN(MXL),
     *               NFIT(MXL),NP(MXL),IFEX(MXL)
      COMMON /RFACT2/THEMIN,THEMAX,RNTHE(10),RDTHE(10),NRTHE(10),
     *SISTHE(10),SILTHE(10)
      COMMON /RFACT3/STLMIN,STLMAX,RNSTL(10),RDSTL(10),NRSTL(10),
     *SISSTL(10),SILSTL(10)
      COMMON /RFACT4/FSQMIN,FSQMAX,RNFSQ(10),RDFSQ(10),NRFSQ(10),
     *SISFSQ(10),SILFSQ(10),ALOGI(10)
      COMMON /SCAL/SCALEI,SCAL2,SCALIN,SIGREJ
      COMMON /CURVE/CMIN(MXL),CMAX(MXL),NCOEFF(MXL),C(MXL,11)
      DOUBLE PRECISION PS
      COMMON /ABSURF/PS(121),NQORD,NRORD,QXMIN,QXMAX,RYMIN,RYMAX
      CHARACTER*10 PGNAME
      COMMON /SYMM1/PGNAME
      COMMON /NWR/NWR
      LOGICAL CENTYP
      COMMON /CENT/CENTYP
      DOUBLE PRECISION Q,R
      EQUIVALENCE (IAA,AAA)
      DATA NTWO16/65536/,NTWO8/256/,NTWO18/262144/,NTWO10/1024/
C
C INITIALISATIONS
C ---------------
C
      LAMLIM = .TRUE.
      IF (NRANGE.EQ.1 .AND. BMINL(1).EQ.0.0 .AND. BMAXL(1).EQ.0.0)
     +    LAMLIM = .FALSE.
      NREF = 0
      NTOT = 0
      NREJ = 0
      NNEG = 0
      NSIGZ = 0
      NSIGRJ = 0
      NREJL = 0
      NREJRS = 0
      NREJTH = 0
      QMINL = 10000.0
      QMAXL = 0.0
      THEMIN = 10E20
      THEMAX = 0.0
      STLMIN = 10E20
      STLMAX = 0.0
      QXMIN = 10E20
      QXMAX = -10E20
      RYMIN = 10E20
      RYMAX = -10E20
C
C READ IN LAUE DATA, CONVERTING INDICES TO UNIQUE SET
C ---------------------------------------------------
C
 20   READ (1,6010,END=50) IH , IK , IL , ALBDA , THETA , II , ISIG , 
     +                     X , Y
      NTOT = NTOT + 1
      II = NINT(SCALIN*II)
      ISIG = NINT(SCALIN*ISIG)
      IF (ISIG.EQ.0) THEN
         NREJ = NREJ + 1
         NSIGZ = NSIGZ + 1
         GO TO 20
      END IF
      IF (SIGREJ.GT.0.0) THEN
         IF (II.LT.SIGREJ*ISIG) THEN
            NREJ = NREJ + 1
            NSIGRJ = NSIGRJ + 1
            GO TO 20
         END IF
      END IF
      IF (AMAXL.GT.0) THEN
         IF (ALBDA.LT.AMINL .OR. ALBDA.GT.AMAXL) THEN
            NREJ = NREJ + 1
            NREJL = NREJL + 1
            GO TO 20
         END IF
      END IF
      IF (LAMLIM) THEN
         DO 30 NN = 1 , NRANGE
            IF (ALBDA.GE.BMINL(NN) .AND. ALBDA.LE.BMAXL(NN)) GO TO 40
 30      CONTINUE
         NREJ = NREJ + 1
         NREJL = NREJL + 1
         GO TO 20
      END IF
 40   IF (RES2.GT.0.0) THEN
         RES = ALBDA/(2*SIN(THETA))
         IF (RES.GT.RES1 .OR. RES.LT.RES2) THEN
            NREJ = NREJ + 1
            NREJRS = NREJRS + 1
            GO TO 20
         END IF
      END IF
      IF (TH2.GT.0.0) THEN
         IF (THETA.LT.TH1 .OR. THETA.GT.TH2) THEN
            NREJ = NREJ + 1
            NREJTH = NREJTH + 1
            GO TO 20
         END IF
      END IF
      IF (II.LE.0) THEN
         NREJ = NREJ + 1
         NNEG = NNEG + 1
         GO TO 20
      END IF
      NREF = NREF + 1
      IF (NREF.GT.MAXREF) THEN
C
C ERROR CONDITIONS
C ----------------
C
         WRITE (6,6120) MAXREF
         STOP
      ELSE
         IF (ALBDA.GT.QMAXL) QMAXL = ALBDA
         IF (ALBDA.LT.QMINL) QMINL = ALBDA
         IF (IH.LT.-127 .OR. IH.GT.127) GO TO 60
         IF (IK.LT.-127 .OR. IK.GT.127) GO TO 60
         IF (IL.LT.-127 .OR. IL.GT.127) GO TO 60
         IF (PGNAME(1:1).EQ.'3' .OR. PGNAME(1:1).EQ.'6'
     +       .OR. PGNAME(1:1).EQ.'4') THEN
            CALL TURN3(IH,IK,IL,ISYM)
         ELSE
            CALL TURNIP(IH,IK,IL,ISYM)
         END IF
         ICS = 0
         IF (CENTYP) THEN
            ICS = 1
            IF (MOD(ISYM,2).EQ.0) ICS = -1
            IHKL(1) = IH
            IHKL(2) = IK
            IHKL(3) = IL
            CALL CENTR(IHKL,IIC)
            IF (IIC.EQ.1) ICS = 0
         ENDIF
         I = NWR*(NREF-1)
         IA(I+1) = NTWO18*(IH+128) + NTWO10*(IK+128) + 
     +     4*(IL+128) + ICS + 1
         IA(I+2) = 10000*ALBDA
         STL = SIN(THETA)/ALBDA
         STL = STL*STL
         AAA = II
         IA(I+3) = IAA
         IF (ISIG.LT.0 .OR. ISIG.GT.32767) THEN
            WRITE (6,6140) IH , IK , IL , ISIG
            STOP
         ELSE
            IA(I+4) = NTWO16*ISYM + ISIG
            IA(I+5) = 10000*THETA
            IF (THETA.GT.THEMAX) THEMAX = THETA
            IF (THETA.LT.THEMIN) THEMIN = THETA
            IF (STL.GT.STLMAX) STLMAX = STL
            IF (STL.LT.STLMIN) STLMIN = STL
            AAA = X
            IA(I+6) = IAA
            AAA = Y
            IA(I+7) = IAA
            IF (IABS.GT.0) THEN
               CALL QRVAL(IA(I+1),IABS,Q,R)
               IF (Q.GT.QXMAX) QXMAX = Q
               IF (Q.LT.QXMIN) QXMIN = Q
               IF (R.GT.RYMAX) RYMAX = R
               IF (R.LT.RYMIN) RYMIN = R
            END IF
            GO TO 20
         END IF
      END IF
C
C END OF REFLECTIONS REACHED
C --------------------------
C
 50   IF (NREF.EQ.0) THEN
         WRITE (6,6150)
         STOP
      ELSE
         WRITE (6,6020) NTOT
         WRITE (6,6030) NREF
         WRITE (6,6040) NREJ
         WRITE (6,6050) NNEG
         WRITE (6,6060) NSIGZ , SIGREJ , NSIGRJ
         WRITE (6,6070) NREJL
         WRITE (6,6080) NREJTH
         WRITE (6,6090) NREJRS
         WRITE (6,6100) QMINL , QMAXL
         IF (NORM) WRITE (6,6110) SCALEI
         AMINL = QMINL
         AMAXL = QMAXL
         IF (NRANGE.EQ.1 .AND. BMINL(1).EQ.0.0 .AND. BMAXL(1).EQ.0.0)
     +       THEN
            BMINL(1) = QMINL
            BMAXL(1) = QMAXL
            IF (.NOT.NORM) THEN
               CMIN(1) = QMINL
               CMAX(1) = QMAXL
            END IF
         END IF
         RETURN
      END IF
 60   WRITE (6,6130) IH , IK , IL
      STOP
C
C FORMAT STATEMENTS
C -----------------
C
 6010 FORMAT (3I5,2F10.5,2I10,2F8.3)
 6020 FORMAT (//,' NUMBER OF LAUE REFLECTIONS READ            =',I6)
 6030 FORMAT (' NUMBER OF REFLECTIONS ACCEPTED             =',I6)
 6040 FORMAT (' NUMBER REJECTED (TOTAL)                    =',I6)
 6050 FORMAT (' NUMBER REJECTED (INTENSITY <= 0)           =',I6)
 6060 FORMAT (' NUMBER REJECTED (SIGMA = 0)                =',I6,/,
     +        ' NUMBER REJECTED (LESS THAN',F6.2,' SIGMA)    =',I6)
 6070 FORMAT (' NUMBER REJECTED (OUTSIDE LAMBDA LIMITS)    =',I6)
 6080 FORMAT (' NUMBER REJECTED (OUTSIDE THETA LIMITS)     =',I6)
 6090 FORMAT (' NUMBER REJECTED (OUTSIDE RESOLUTION LIMITS)=',I6)
 6100 FORMAT (//,' LAMBDA LIMITS FROM LAUE FILE ARE',F7.3,',',F7.3,
     +        ' ANGSTROMS')
 6110 FORMAT (//,' INPUT STANDARD INTENSITIES ARE DIVIDED BY',F7.1,
     +        ' DURING SCALING PROCESS')
C
 6120 FORMAT (/,' **NUMBER OF REFLECTIONS MAY NOT EXCEED',I6,' **')
 6130 FORMAT (/,' **INDICES OUT OF RANGE -127 TO 127, HKL=',3I4,' **')
 6140 FORMAT (/,' **SIGMA OUT OF RANGE 0-32767 FOR REFLECTION HKL=',3I4,
     +        '  SIGMA =',I8,' **')
 6150 FORMAT (/,' **NO LAUE REFLECTIONS ACCEPTED**')
      END
**==refsrt.f  
C
C SUBROUTINE 'REFSRT'
C ===================
C
C THIS SUBROUTINE IS USED TO SORT THE REFLECTION DATA HELD IN CORE. IT
C IS A MODIFIED VERSION OF THE SUBROUTINE 'AHVSOR' FROM ELEANOR DODSON'S
C VERSION OF THE PROGRAM 'CAD'. THIS VERSION SORTS ON A WORD CONTAINING
C THE PACKED INDICES.
C
      SUBROUTINE REFSRT(IA,NREF,NWRD,IDIM)
C
C PARAMETERS
C
C         IA (I/O)  INTEGER ARRAY CONTAINING THE REFLECTION DATA WITH
C                   NWRD WORDS FOR REFLN 1 FOLLOWED BY NWRD WORDS FOR
C                   EACH REFLECTION UP TO NREF REFLECTIONS  (NWRD <=10)
C       NREF (I)    THE NO. OF REFLECTIONS TO BE SORTED
C       NWRD (I)    NUMBER OF WORDS/REFLECTION (INDICES PACKED IN FIRST
C                   WORD FOR EACH REFLECTION)
C       IDIM (I)    THE NUMBER OF ELEMENTS IN THE ARRAY IA
C
C SPECIFICATION STATEMENTS
C ------------------------
C
      IMPLICIT INTEGER(A-Z)
      INTEGER IA(IDIM),ITEMP,CUT(10)
      INTEGER SAVE(2,16)
      DATA CUTOFF/20/
C
C SORT DATA
C ---------
C
      STACK = 0
      LOW = 1
      HIGH = NREF
 20   CHUNK = HIGH - LOW
      IF (CHUNK.GT.CUTOFF) THEN
         MIDDLE = (LOW+HIGH)/2
         IF (ICMP(IA,LOW,IA,MIDDLE,NWRD).GT.0) THEN
            J1 = NWRD*(LOW-1)
            J2 = NWRD*(MIDDLE-1)
            DO 30 I = 1 , NWRD
               J1 = J1 + 1
               J2 = J2 + 1
               ITEMP = IA(J1)
               IA(J1) = IA(J2)
               IA(J2) = ITEMP
 30         CONTINUE
         END IF
         IF (ICMP(IA,MIDDLE,IA,HIGH,NWRD).GT.0) THEN
            J1 = NWRD*(MIDDLE-1)
            J2 = NWRD*(HIGH-1)
            DO 40 I = 1 , NWRD
               J1 = J1 + 1
               J2 = J2 + 1
               ITEMP = IA(J1)
               IA(J1) = IA(J2)
               IA(J2) = ITEMP
 40         CONTINUE
            IF (ICMP(IA,LOW,IA,MIDDLE,NWRD).GT.0) THEN
               J1 = NWRD*(LOW-1)
               J2 = NWRD*(MIDDLE-1)
               DO 50 I = 1 , NWRD
                  J1 = J1 + 1
                  J2 = J2 + 1
                  ITEMP = IA(J1)
                  IA(J1) = IA(J2)
                  IA(J2) = ITEMP
 50            CONTINUE
            END IF
         END IF
         J = NWRD*(MIDDLE-1)
         DO 60 I = 1 , NWRD
            J = J + 1
            CUT(I) = IA(J)
 60      CONTINUE
         UP = LOW
         DOWN = HIGH
 70      DOWN = DOWN - 1
         IF (ICMP(CUT,1,IA,DOWN,NWRD).LT.0) GO TO 70
 80      UP = UP + 1
         IF (ICMP(CUT,1,IA,UP,NWRD).GT.0) GO TO 80
         IF (UP.LT.DOWN) THEN
            J1 = NWRD*(UP-1)
            J2 = NWRD*(DOWN-1)
            DO 90 I = 1 , NWRD
               J1 = J1 + 1
               J2 = J2 + 1
               ITEMP = IA(J1)
               IA(J1) = IA(J2)
               IA(J2) = ITEMP
 90         CONTINUE
            GO TO 70
         ELSE IF (UP.EQ.DOWN) THEN
            DOWN = DOWN - 1
            UP = UP + 1
         END IF
      ELSE
         IF (CHUNK.GE.2) THEN
            MM = HIGH - 1
            DO 120 M = LOW , MM
               MIN = M
               MMM = M + 1
               DO 100 N = MMM , HIGH
                  IF (ICMP(IA,MIN,IA,N,NWRD).GT.0) THEN
                     MIN = N
                  END IF
 100           CONTINUE
               IF (MIN.NE.M) THEN
                  J1 = NWRD*(M-1)
                  J2 = NWRD*(MIN-1)
                  DO 110 I = 1 , NWRD
                     J1 = J1 + 1
                     J2 = J2 + 1
                     ITEMP = IA(J1)
                     IA(J1) = IA(J2)
                     IA(J2) = ITEMP
 110              CONTINUE
               END IF
 120        CONTINUE
         END IF
         IF (STACK.EQ.0) RETURN
         LOW = SAVE(1,STACK)
         HIGH = SAVE(2,STACK)
         STACK = STACK - 1
         GO TO 20
      END IF
      STACK = STACK + 1
      IF (STACK.GT.16) THEN
C
C ERROR CONDITION
C ----------------
C
         WRITE (6,6010) STACK
         STOP
      ELSE
         BOTTOM = DOWN - LOW
         TOP = HIGH - UP
         IF (BOTTOM.GT.TOP) THEN
            SAVE(1,STACK) = LOW
            SAVE(2,STACK) = DOWN
            LOW = UP
         ELSE
            SAVE(1,STACK) = UP
            SAVE(2,STACK) = HIGH
            HIGH = DOWN
         END IF
         GO TO 20
      END IF
C
C FORMAT STATEMENTS
C -----------------
C
 6010 FORMAT (' **STACK LIMIT OF 16 EXCEEDED IN REFSOR**')
      END
**==icmp.f  
C
C FUNCTION 'ICMP'
C ===============
C
C THIS FUNCTION IS USED TO COMPARE THE VALUES OF THE SORT KEYS
C FOR 2 REFLECTIONS
C
      FUNCTION ICMP(IARR1,I1,IARR2,I2,NWRD)
C
C PARAMETERS
C
C        IARR1 (I)   ARRAY CONTAINING 1ST REFLN
C           I1 (I)   POINTER TO REFLECTION IN IARR1
C        IARR2 (I)   ARRAY CONTAINING 2ND REFLN
C           I2 (I)   POINTER TO REFLECTIONIIN IARR2
C         NWRD (I)   NO. OF WORDS/REFLN
C
C FUNCTION RETURN VALUES  -1  REFLN 1 LT REFLN 2
C                          0  REFLN 1 EQ REFLN 2
C                          1  REFLN 1 GT REFLN 2
C
C SPECIFICATION STATEMENTS
C ------------------------
C
      INTEGER IARR1(*),IARR2(*)
C
C COMPARE KEYS
C ------------
C
      J1 = NWRD*I1 - NWRD + 1
      J2 = NWRD*I2 - NWRD + 1
      IF (IARR1(J1).LT.IARR2(J2)) THEN
         ICMP = -1
         RETURN
      ELSE IF (IARR1(J1).EQ.IARR2(J2)) THEN
         ICMP = 0
         RETURN
      ELSE
         ICMP = 1
         RETURN
      END IF
      END
**==lamsrt.f  
C
C SUBROUTINE 'LAMSRT'
C ===================
C
C SORT EQUIVALENTS INTO LAMBDA ORDER WITHIN LAUESCALE REFLECTION ARRAY
C
      SUBROUTINE LAMSRT(IA,NREF)
C
C PARAMETERS
C
C          IA (I/O) INTEGER ARRAY CONTAINING THE REFLECTION DATA.
C        NREF (I)   THE NUMBER OF REFLECTIONS
C
C SPECIFICATION STATEMENTS
C ------------------------
C
      DIMENSION IA(*)
      LOGICAL END,SWAP
C
C SORT REFLECTIONS ON LAMBDA
C --------------------------
C
      END = .FALSE.
      M1 = 1
      ILAST = IA(1)
      K = 0
C
C STEP THROUGH THE REFLECTIONS
C
 20   K = K + 1
      IF (K.GT.NREF) THEN
         END = .TRUE.
         GO TO 30
      END IF
      JJ = 7*K - 6
      IF (IA(JJ).EQ.ILAST) GO TO 20
 30   M2 = K - 1
      MM = 7*M1 - 13
      NVAL = M2 - M1 + 1
C
C SORT EACH REFLECTION ON LAMBDA
C
      IF (NVAL.GE.2) THEN
         NMAX = NVAL - 1
         DO 60 N = 1 , NMAX
            LMAX = NVAL - N
            SWAP = .FALSE.
            DO 50 L = 1 , LMAX
               L1 = MM + 7*L
               L2 = L1 + 7
               IF (IA(L1+1).GT.IA(L2+1)) THEN
                  SWAP = .TRUE.
                  DO 40 II = 0 , 6
                     ITEMP = IA(L1+II)
                     IA(L1+II) = IA(L2+II)
                     IA(L2+II) = ITEMP
 40               CONTINUE
               END IF
 50         CONTINUE
            IF (.NOT.SWAP) GO TO 70
 60      CONTINUE
      END IF
 70   IF (END) THEN
         RETURN
      ELSE
         M1 = K
         ILAST = IA(JJ)
         GO TO 20
      END IF
      END
**==irange.f  
C
C FUNCTION 'IRANGE'
C =================
C
C FIND A RANGE NUMBER
C
      FUNCTION IRANGE(X,XMIN,XMAX,NX)
C
C PARAMETERS
C
C           X (I)   THE VALUE FOR WHICH THE RANGE NO. IS TO BE FOUND
C        XMIN (I)   THE MINIMUM X VALUE
C        XMAX (I)   THE MAXIMUM X VALUE
C          NX (I)   THE NUMBER OR RANGES (FROM XMIN TO XMAX)
C
C  RETURNS  IRANGE = THE RANGE NUMBER
C
C CALCULATE THE RANGE NUMBER
C --------------------------
C
      IRANGE = 1
      IF (XMAX.LE.XMIN) RETURN
      IRANGE = NX*(X-XMIN)/(XMAX-XMIN) + 1
      IF (IRANGE.LT.1) IRANGE = 1
      IF (IRANGE.GT.NX) IRANGE = NX
      RETURN
      END
**==lambin.f  
C
C FUNCTION 'LAMBIN'
C =================
C
C FIND THE BIN NUMBER FOR A GIVEN LAMBDA VALUE
C
       FUNCTION LAMBIN(BLAM)
C
C PARAMETERS
C
C       BLAM (I)   THE LAMBDA VALUE FOR WHICH THE BIN NO. IS REQUIRED
C
C RETURNS  LAMBIN = THE BIN NUMBER (0 IF OUTSIDE RANGES)
C
C
      PARAMETER(MXR=25000,MXL=8)
C
C
C SPECIFICATION STATEMENTS
C ------------------------
C
      COMMON /CDATA/AMINL,AMAXL,RES1,RES2,TH1,TH2,BFAC,NRANGE,IBATCH,
     *              NITER,IABS,LASTIT,KBATCH
      COMMON /WGHTS/IWTS,AWTS,IWTL,AWTL,IWMEAN
      COMMON /RANGES/BMINL(MXL),BMAXL(MXL),BINS(MXL),NBIN(MXL),
     *               NFIT(MXL),NP(MXL),IFEX(MXL)
C
C CALCULATE BIN NUMBER
C --------------------
C
      DO 20 NN = 1 , NRANGE
         IF (BLAM.GE.BMINL(NN) .AND. BLAM.LE.BMAXL(NN)) GO TO 30
 20   CONTINUE
      LAMBIN = 0
      RETURN
 30   IBIN = NBIN(NN)*(BLAM-BMINL(NN))/(BMAXL(NN)-BMINL(NN)) + 1
      IF (IBIN.LT.1) IBIN = 1
      IF (IBIN.GT.NBIN(NN)) IBIN = NBIN(NN)
      LAMBIN = IBIN + NP(NN)
      RETURN
      END
**==colect.f  
C
C SUBROUTINE 'COLECT'
C ===================
C
C THIS SUBROUTINE COLLECTS THE DATA TO BE USED FOR SCALING FROM THE
C LAUE DATA AND THE INPUT MTZ FILE INTO WAVELENGTH BINS.
C
      SUBROUTINE COLECT(IA,IMAX,NREF)
C
C PARAMETERS
C
C         IA (I)   ARRAY TO HOLD LAUE DATA (5 WORDS/REFLN)
C       IMAX (I)   NO. OF WORDS IN 'IA'
C       NREF (I)   NUMBER OF LAUE REFLECTIONS
C
      PARAMETER(MXR=25000,MXL=8)
C
C SPECIFICATION STATEMENTS
C ------------------------
C
      LOGICAL EOF
      DIMENSION IA(IMAX),IDXL(3),ADATA(10),IDATA(3)
      LOGICAL NORM,STFIL,FILE1,FILE2,FILAGR,AFSCAL,SHLX1,SHLX2,SHLX3
      COMMON /FLAGS/NORM,STFIL,FILE1,FILE2,FILAGR,AFSCAL,SHLX1,SHLX2,
     +              SHLX3
      LOGICAL UPDATB,NRMCAL,ABSCAL
      COMMON /SCFLAG/UPDATB,NRMCAL,ABSCAL
      COMMON /CDATA/AMINL,AMAXL,RES1,RES2,TH1,TH2,BFAC,NRANGE,IBATCH,
     *              NITER,IABS,LASTIT,KBATCH
      COMMON /WGHTS/IWTS,AWTS,IWTL,AWTL,IWMEAN
      COMMON /RANGES/BMINL(MXL),BMAXL(MXL),BINS(MXL),NBIN(MXL),
     *               NFIT(MXL),NP(MXL),IFEX(MXL)
      COMMON /BINS/NBTOT,SWL(100),SWIL(100),SWS(100),SWIS(100),W(100),
     *             R(100),ALAM(100),NVAL(100),AIS(100),AIL(100)
      COMMON /RFACT4/FSQMIN,FSQMAX,RNFSQ(10),RDFSQ(10),NRFSQ(10),
     *SISFSQ(10),SILFSQ(10),ALOGI(10)
      COMMON /SCAL/SCALEI,SCAL2,SCALIN,SIGREJ
      COMMON /NWR/NWR
      LOGICAL LAMOPT
      COMMON /LAMO/ LAMOPT, DLAM
      LOGICAL CENTYP
      COMMON /CENT/CENTYP
      EQUIVALENCE (IAA,AAA)
      DATA NTWO16/65536/
C
C INITIALISE BIN ARRAYS
C ---------------------
C
      NMATCH = 0
      DO 20 I = 1 , 100
         NVAL(I) = 0
         SWL(I) = 0.0
         SWIL(I) = 0.0
         SWS(I) = 0.0
         SWIS(I) = 0.0
         W(I) = 0.0
         R(I) = 0.0
         ALAM(I) = 0.0
         AIS(I) = 0.0
         AIL(I) = 0.0
 20   CONTINUE
      FSQMIN = 10E20
      FSQMAX = 0.0
      NZEROW = 0
      NR = 1
      CALL LRREWD(1)
C
C GET FIRST LAUE AND FIRST MTZ REFLECTION
C ---------------------------------------
C
      CALL GETHKL(IA(1),IDXL,ICS)
      BLMBDA = IA(2)/10000.0
      CALL LRREFF(1,RESOL,ADATA,EOF)
      IF (EOF) GO TO 50
      IDATA(1) = NINT(ADATA(1))
      IDATA(2) = NINT(ADATA(2))
      IDATA(3) = NINT(ADATA(3))
      ALMBDA = 0.0
      KCS = 0
      IF (LAMOPT) THEN
          ALMBDA = ADATA(6)
          IF (CENTYP) KCS = NINT(ADATA(7))
      ELSE
          IF (CENTYP) KCS = NINT(ADATA(6))
      ENDIF
C
C COMPARE LAUE AND MTZ REFLECTIONS AND READ FURTHER REFLECTIONS FROM
C ------------------------------------------------------------------
C
C WHICHEVER SET IS REQUIRED TILL A MATCH IS FOUND
C -----------------------------------------------
C
 30   IF (ICMPR(IDXL,IDATA,ALMBDA,BLMBDA,KCS,ICS,LAMOPT,DLAM).LT.0) THEN 
      ELSE IF (ICMPR(IDXL,IDATA,ALMBDA,BLMBDA,KCS,ICS,LAMOPT,DLAM)
     +        .EQ.0) THEN
C
C MATCHING REFLECTION FOUND, COLLECT DATA FOR ANALYSIS AND CURVE FIT
C ------------------------------------------------------------------
C
         IB = NWR*(NR-1)
         NMATCH = NMATCH + 1
         SC = SCFAC(IA(IB+1),.FALSE.,ABSCAL)
         SIGL = MOD(IA(IB+4),NTWO16)*SC
         IAA = IA(IB+3)
         FSQL = AAA*SC
         BLAM = IA(IB+2)/10000.0
         FS = ADATA(4)
         FSQS = ADATA(4)*ADATA(4)
         IF (FSQS.GT.FSQMAX) FSQMAX = FSQS
         IF (FSQS.LT.FSQMIN) FSQMIN = FSQS
         FSQS = FSQS/SCALEI
         IF (SIGL.NE.0.0 .AND. FSQL.NE.0.0) THEN
            IF (FSQS.NE.0.0 .AND. ADATA(5).NE.0.0) THEN
               SIGS = 2*FS*ADATA(5)/SCALEI
               WL = 1/(SIGL*SIGL)
               WS = 1/(SIGS*SIGS)
               IF (IWTL.EQ.1) WL = 1.0
               IF (IWTL.GT.1) WL = 1.0 + AWTL*FSQL
               IF (WL.LT.0.0) WL = 0.0
               IF (IWTS.EQ.1) WS = 1.0
               IF (IWTS.GT.1) WS = 1.0 + AWTS*FSQS
               IF (WS.LT.0.0) WS = 0.0
               IF (WS.LE.0.0 .OR. WL.LE.0.0) THEN
                  NZEROW = NZEROW + 1
                  GO TO 40
               END IF
               IBIN = LAMBIN(BLAM)
               IF (IBIN.NE.0) THEN
                  NVAL(IBIN) = NVAL(IBIN) + 1
                  SWIL(IBIN) = SWIL(IBIN) + FSQL*WL
                  SWL(IBIN) = SWL(IBIN) + WL
                  SWIS(IBIN) = SWIS(IBIN) + FSQS*WS
                  SWS(IBIN) = SWS(IBIN) + WS
                  AIS(IBIN) = AIS(IBIN) + FSQS
                  AIL(IBIN) = AIL(IBIN) + FSQL
                  CALL DIAG1(IDXL,FSQL,FSQS,BLAM,IBIN)
               END IF
            END IF
         END IF
      ELSE
C
C NO MATCH, LAUE>STANDARD, GET NEXT STANDARD
C
         CALL LRREFF(1,RESOL,ADATA,EOF)
         IF (EOF) GO TO 50
         IDATA(1) = NINT(ADATA(1))
         IDATA(2) = NINT(ADATA(2))
         IDATA(3) = NINT(ADATA(3))
         ALMBDA = 0.0
         KCS = 0
         IF (LAMOPT) THEN
             ALMBDA = ADATA(6)
             IF (CENTYP) KCS = NINT(ADATA(7))
         ELSE
             IF (CENTYP) KCS = NINT(ADATA(6))
         ENDIF
         GO TO 30
      END IF
C
C NO MATCH, LAUE<STANDARD, GET NEXT LAUE
C
 40   NR = NR + 1
      IF (NR.LE.NREF) THEN
         I = NWR*NR - NWR + 1
         CALL GETHKL(IA(I),IDXL,ICS)
         BLMBDA = IA(I+1)/10000.0
         GO TO 30
      END IF
C
C END OF DATA REACHED, GET POINTS FOR CURVE FITTING
C -------------------------------------------------
C
 50   IF (.NOT.NORM) RETURN
      DO 70 NN = 1 , NRANGE
         NB = 0
         DO 60 NA = 1 , NBIN(NN)
            N = NP(NN) + NA
            ALAM(N) = BMINL(NN) + (NA-0.5)*(BMAXL(NN)-BMINL(NN))
     +                /NBIN(NN)
            IF (NVAL(N).GT.0) THEN
               AIS(N) = AIS(N)/NVAL(N)
               AIL(N) = AIL(N)/NVAL(N)
               NB = NB + 1
            END IF
 60      CONTINUE
         IF (NB.LT.2) GO TO 110
         IF (NB.LT.NBIN(NN)) WRITE (6,6010)
 70   CONTINUE
      NW = 0
      AW = 0
      W1 = 1.0
      DO 80 N = 1 , NBTOT
         IF (SWL(N).GT.0.0) THEN
            FSQLM = SWIL(N)/SWL(N)
            IF (SWS(N).GT.0.0) THEN
               FSQSM = SWIS(N)/SWS(N)
               IF (FSQSM.NE.0.0) THEN
                  R(N) = FSQLM/FSQSM
                  SIGM2 = R(N)*(1/(SWL(N)*FSQLM)+1/(SWS(N)*FSQLM))
                  W(N) = 1/SQRT(SIGM2)
C NOTE: WEIGHT FOR CURVE FIT ROUTINE IS 1/SIGMA NOT 1/SIGMA**2
                  NW = NW + NVAL(N)
                  AW = AW + W(N)
               END IF
            END IF
         END IF
 80   CONTINUE
C
C SET WEIGHTS FOR SINGLE REFLECTION BINS
C --------------------------------------
C
      IF (NW.GT.0) W1 = AW/NW
      DO 90 N = 1 , NBTOT
         IF (NVAL(N).EQ.1) THEN
            W(N) = W1
         END IF
 90   CONTINUE
      WRITE (6,6020) NMATCH
      IF (NZEROW.GT.0) WRITE (6,6030) NZEROW
      RETURN
C
C ERROR CONDITIONS
C ----------------
C
 110  WRITE (6,6050) NN
      STOP
C
C FORMAT STATEMENTS
C -----------------
C
 6010 FORMAT (/,' *WARNING* EMPTY BIN(S) PRESENT')
 6020 FORMAT (//,' NUMBER OF MATCHING REFLECTIONS =',I6,////)
 6030 FORMAT (' NUMBER OF REFLECTIONS EXCLUDED FROM CALCULATION ',
     +        ' OF SCALING POINTS BECAUSE OF ZERO WEIGHTS =',I5,///)
C
 6050 FORMAT (/,' **NUMBER OF BINS CONTAINING DATA IS LESS THAN TWO',
     +        ' FOR RANGE',I2,' **')
      END
**==gethkl.f  
C
C SUBROUTINE 'GETHKL'
C ===================
C
C THIS SUBROUTINE UNPACKS THE INDICES FOR A LAUE REFLECTION
C
      SUBROUTINE GETHKL(IDX,IDXL,ICS)
C
C PARAMETERS
C
C      IDX (I)   WORD CONTAINING PACKED INDICES
C     IDXL (O)   3 WORD ARRAY CONTAINING H, K, L.
C      ICS (O)   CENTERIG SIGN FLAG
C
C SPECIFICATION STATEMENTS
C ------------------------
C
      DIMENSION IDXL(3)
      DATA NTWO18/262144/,NTWO10/1024/,NTWO16/65536/
C
C UNPACK INDICES
C --------------
C
      I = IDX
      IDXL(1) = I/NTWO18
      I = I - NTWO18*IDXL(1)
      IDXL(1) = IDXL(1) - 128
      IDXL(2) = I/NTWO10
      I = I - NTWO10*IDXL(2)
      IDXL(2) = IDXL(2) - 128
      IDXL(3) = I/4
      I = I - 4*IDXL(3)
      IDXL(3) = IDXL(3) -128
      ICS = I - 1
      RETURN
      END
**==icmpr.f  
C
C FUNCTION 'ICMPR'
C ================
C
C THIS FUNCTION COMPARES THE INDICES OF A LAUE REFLECTION WITH THE
C INDICES OF A STANDARD REFLECTION (WITH OPTIONAL LAMBDA CHECK)
C
      FUNCTION ICMPR(IHKL1,IHKL2,ALMBDA,BLMBDA,KCS,ICS,LAMOPT,DLAM)
C
C PARAMETERS
C
C       IHKL1 (I)   3 WORD ARRAY HOLDING LAUE REFLN INDICES
C       IHKL2 (I)   3 WORD ARRAY HOLDING STANDARD REFLN INDICES
C      ALMBDA (I)   LAMBDA VALUE OF STANDARD REFLN
C      BLMBDA (I)   LAMBDA VALUE OF LAUE REFLN
C         KCS (I)   SIGN FLAG FOR STANDARD REFLN
C         ICS (I)   SIGN FLAG FOR LAUE REFLECTION
C      LAMOPT (I)   FLAG = .TRUE. INCLUDE LAMBDA CHECK, = .FALSE. DO NOT
C        DLAM (I)   TOLERANCE FOR LAMBDA CHECK
C
C RETURNS  -1  LAUE < STANDARD
C           0  LAUE = STANDARD
C           1  LAUE > STANDARD
C
C SPECIFICATION STATEMENTS
C ------------------------
C
      DIMENSION IHKL1(3),IHKL2(3)
      LOGICAL LAMOPT
C
C COMPARE INDICES
C ---------------
C
      ICMPR = 0
      DO 20 I = 1 , 3
         IF (IHKL1(I).LT.IHKL2(I)) GO TO 30
         IF (IHKL1(I).NE.IHKL2(I)) GO TO 40
 20   CONTINUE
      IF (ICS.LT.KCS) GO TO 30
      IF (ICS.GT.KCS) GO TO 40
      IF (.NOT.LAMOPT) RETURN
      IF (ABS(ALMBDA-BLMBDA).LE.DLAM) RETURN
      IF (BLMBDA.GT.ALMBDA) GO TO 40
C
C LAUE < STANDARD
C
 30   ICMPR = -1
      RETURN
C
C LAUE > STANDARD
C
 40   ICMPR = 1
      RETURN
      END
**==curfit.f  
C
C SUBROUTINE 'CURFIT'
C ===================
C
C CALCULATE CURVE FOR NORMALISATION OF THE DATA
C
C
      SUBROUTINE CURFIT
C
C PARAMETERS
C
C NONE
C
      PARAMETER(MXR=25000,MXL=8)
C
C SPECIFICATION STATEMENTS
C ------------------------
C
      CHARACTER*1 SEL
      DOUBLE PRECISION COEFFS(11),RMS(11),CW(100),CR(100),CLAM(100)
      COMMON /CDATA/AMINL,AMAXL,RES1,RES2,TH1,TH2,BFAC,NRANGE,IBATCH,
     *              NITER,IABS,LASTIT,KBATCH
      COMMON /WGHTS/IWTS,AWTS,IWTL,AWTL,IWMEAN
      COMMON /RANGES/BMINL(MXL),BMAXL(MXL),BINS(MXL),NBIN(MXL),
     *               NFIT(MXL),NP(MXL),IFEX(MXL)
      COMMON /CURVE/CMIN(MXL),CMAX(MXL),NCOEFF(MXL),C(MXL,11)
      COMMON /BINS/NBTOT,SWL(100),SWIL(100),SWS(100),SWIS(100),W(100),
     *             R(100),ALAM(100),NVAL(100),AIS(100),AIL(100)
C
C PREPARE POINTS FOR CURVE FITTING
C --------------------------------
C
      CALL HEADXX(6,'CURVE FITTING RESULTS$','$','0',2)
      DO 60 NN = 1 , NRANGE
         WRITE (6,6010) NN
         DO 20 I = 1 , 11
            C(NN,I) = 0.0
 20      CONTINUE
         KMX = 10
         NPOINT = 0
         NB1 = NP(NN) + 1
         NB2 = NP(NN) + NBIN(NN)
         IF ((MOD(IFEX(NN),2).NE.0) .AND. NB2-NB1.GT.1) THEN
            AL = 1.5*ALAM(NB1) - 0.5*ALAM(NB1+1)
            AL0 = ALAM(NB1+2)
            SC0 = R(NB1+2)
            AL1 = ALAM(NB1+1)
            SC1 = R(NB1+1)
            AL2 = ALAM(NB1)
            SC2 = R(NB1)
            NPOINT = 1
            CW(NPOINT) = W(NB1)
            CR(NPOINT) = EXTRAP(AL,AL0,SC0,AL1,SC1,AL2,SC2,6)
            CLAM(NPOINT) = AL
            WRITE (6,6050) CR(1)
         END IF
         DO 30 N = NB1 , NB2
            NPOINT = NPOINT + 1
            CW(NPOINT) = W(N)
            CR(NPOINT) = R(N)
            CLAM(NPOINT) = ALAM(N)
 30      CONTINUE
         IF (IFEX(NN).GT.1 .AND. NB2-NB1.GT.1) THEN
            AL = 1.5*ALAM(NB2) - 0.5*ALAM(NB2-1)
            AL0 = ALAM(NB2-2)
            SC0 = R(NB2-2)
            AL1 = ALAM(NB2-1)
            SC1 = R(NB2-1)
            AL2 = ALAM(NB2)
            SC2 = R(NB2)
            NPOINT = NPOINT + 1
            CW(NPOINT) = W(NB2)
            CR(NPOINT) = EXTRAP(AL,AL0,SC0,AL1,SC1,AL2,SC2,6)
            CLAM(NPOINT) = AL
            WRITE (6,6060) CR(NPOINT)
         END IF
         WRITE (6,6020) NPOINT
         IF (NPOINT.LE.KMX) KMX = NPOINT - 1
         IF (KMX.LT.NFIT(NN)) THEN
            NFIT(NN) = KMX
            WRITE (6,6030) NFIT(NN)
         END IF
         IF (NPOINT.LT.2) GO TO 70
         CMIN(NN) = CLAM(1)
         CMAX(NN) = CLAM(NPOINT)
C#### Invert curve
         DO 35 I = 1, NPOINT
            CR(I) = 1.0/CR(I)
35       CONTINUE
C####
         CALL POLYFT(NPOINT,CLAM,CR,CW,KMX,NFIT(NN),RMS,COEFFS)
         WRITE (6,6040) NFIT(NN) , (COEFFS(I),I=1,NFIT(NN)+1)
         WRITE (6,6070)
         DO 40 N = 1 , KMX + 1
            SEL = ' '
            IF (N.EQ.NFIT(NN)+1) SEL = '#'
            WRITE (6,6080) N - 1 , SEL , RMS(N)
 40      CONTINUE
         NCOEFF(NN) = NFIT(NN) + 1
         DO 50 I = 1 , NCOEFF(NN)
            C(NN,I) = COEFFS(I)
 50      CONTINUE
 60   CONTINUE
      RETURN
C
C ERROR CONDITIONS
C ----------------
C
 70   WRITE (6,6090)
      STOP
C
C FORMAT STATEMENTS
C -----------------
C
 6010 FORMAT (///,' RANGE',I2)
 6020 FORMAT (//,' NUMBER OF POINTS TO WHICH CURVE IS TO BE FITTED =',
     +        I3)
 6030 FORMAT (/,' *WARNING* MAXIMUM ORDER FOR CURVE FITTING HAS BEEN',
     +        ' REDUCED TO ',I2)
 6040 FORMAT (//,' POLYNOMIAL COEFFICIENTS   P 0  TO  P',I2,//,
     +        (' ',6E20.8))
 6050 FORMAT (/,' EXTRAPOLATED POINT ADDED AT  LOW WAVELENGTH END',
     +        ' OF RANGE:',F10.3)
 6060 FORMAT (/,' EXTRAPOLATED POINT ADDED AT HIGH WAVELENGTH END',
     +        ' OF RANGE:',F10.3)
 6070 FORMAT (//,' MEAN SQUARE DEVIATIONS',//,' ORDER    DEVIATION',/)
 6080 FORMAT (' ',I3,A1,E11.3)
 6090 FORMAT (/,' **LESS THAN TWO POINTS DETERMINED FOR CURVE FITTING**'
     +        )
      END
**==curval.f  
C
C FUNCTION 'CURVAL'
C =================
C
C THIS FUNCTION RETURNS THE SCALING FUNCTION VALUE FOR A GIVEN
C LAMBDA. FOR POINTS OUTSIDE THE RANGE OF THE FITTED CURVE, THE
C SUBROUTINE ALLOWS POLYNOMIAL EVALUATION OUTSIDE THE RANGE OF THE
C CURVE FITTED POINTS.
C
      FUNCTION CURVAL(X)
C
C PARAMETERS
C
C           X (I)   THE LAMBDA VALUE
C
C RETURNS   CURVAL, THE SCALING FUNCTION VALUE (0 IF OUTSIDE RANGES)
C
      PARAMETER(MXL=8)
C
C SPECIFICATION STATEMENTS
C ------------------------
C
      LOGICAL CFLAG
      COMMON /CDATA/AMINL,AMAXL,RES1,RES2,TH1,TH2,BFAC,NRANGE,IBATCH,
     *              NITER,IABS,LASTIT,KBATCH
      COMMON /RANGES/BMINL(MXL),BMAXL(MXL),BINS(MXL),NBIN(MXL),
     *               NFIT(MXL),NP(MXL),IFEX(MXL)
      COMMON /CURVE/CMIN(MXL),CMAX(MXL),NCOEFF(MXL),C(MXL,11)
C
C CALCULATE REQUIRED SCALING VALUE
C --------------------------------
C
      DO 20 N = 1 , NRANGE
         NN = N
         IF (X.GE.BMINL(NN) .AND. X.LE.BMAXL(NN)) GO TO 30
 20   CONTINUE
      CURVAL = 0.0
      RETURN
 30   CURVAL = C(NN,NCOEFF(NN))
      DO 40 J = NCOEFF(NN) - 1 , 1 , -1
         CURVAL = CURVAL*X + C(NN,J)
 40   CONTINUE
C#### Re-invert curve 
      CURVAL = 1.0/CURVAL
C####
      RETURN
      END
**==extrap.f  
C
C FUNCTION 'EXTRAP'
C =================
C
C THIS FUNCTION PERFORMS A QUADRATIC EXTRAPOLATION FROM THREE POINTS
C OF A GIVEN FUNCTION
C
      FUNCTION EXTRAP(X,X0,Y0,X1,Y1,X2,Y2,MSG)
C
C PARAMETERS
C
C      X (I)   THE VALUE AT WHICH THE EXTRAPOLATED FUNCTION IS TO BE
C              CALCULATED. X MUST BE >= X1. FOR VALUES BETWEEN X1 AND
C              X2 A QUADRATIC INTERPOLATION WILL TAKE PLACE.
C
C      X0,Y0,X1,Y1,X2,Y2 (I)  THE X AND Y COORDINATES OF THREE POINTS IN
C                             THE FUNCTION. X2>X1>X0
C
C    MSG (I)   THE UNIT NUMBER FOR PRINTING AN ERROR MESSAGE. IF THIS
C              IS GIVEN THE PROGRAM WILL STOP AFTER AN ERROR. IF A VALUE
C              OF 0 IS GIVEN, NO MESSAGE WILL BE PRINTED AND A FUNCTION
C              VALUE OF 0 WILL BE RETURNED.
C
C  RETURNS  EXTRAP, THE EXTRAPOLATED (OR INTERPOLATED) VALUE.
C
C NOTE: CARE SHOULD BE TAKEN WHEN USING EXTRAPOLATION AS VALUES MAY
C       BECOME VERY INACCURATE ESPECIALLY IF X IS MUCH GREATER THAN X2
C
C PERFORM EXTRAPOLATION
C ---------------------
C
      EXTRAP = 0.0
      IF (X.LT.X1 .OR. X1.LE.X0 .OR. X2.LE.X1) THEN
         IF (X.GT.X1 .OR. X1.GE.X0 .OR. X2.GE.X1) THEN
C
C ERROR CONDITION
C ---------------
C
            IF (MSG.LE.0) RETURN
            WRITE (MSG,6010) X , X0 , X1 , X2
            STOP
         END IF
      END IF
      YY0 = (Y0*(X-X1)-Y1*(X-X0))/(X0-X1)
      YY1 = (Y1*(X-X2)-Y2*(X-X1))/(X1-X2)
      EXTRAP = (YY0*(X-X2)-YY1*(X-X0))/(X0-X2)
      RETURN
C
C FORMAT STATEMENTS
C -----------------
C
 6010 FORMAT (/,' **INVALID PARAMETERS FOR EXTRAPOLATION**',4E12.3)
      END
**==table1.f  
 
C
C SUBROUTINE 'TABLE1'
C ===================
C
C OUTPUT THE STATISTICS FOR THE STANDARD AND LAUE DATA USED IN THE
C CALCULATION OF THE NORMALISATION CURVE
C
      SUBROUTINE TABLE1
C
C PARAMETERS
C
C NONE
C
      PARAMETER(MXR=25000,MXL=8)
C
C SPECIFICATION STATEMENTS
C ------------------------
C
      CHARACTER*32 GRLIN,GR0,GR1
      CHARACTER*1 STAR
      COMMON /CDATA/AMINL,AMAXL,RES1,RES2,TH1,TH2,BFAC,NRANGE,IBATCH,
     *              NITER,IABS,LASTIT,KBATCH
      COMMON /WGHTS/IWTS,AWTS,IWTL,AWTL,IWMEAN
      COMMON /RANGES/BMINL(MXL),BMAXL(MXL),BINS(MXL),NBIN(MXL),
     *               NFIT(MXL),NP(MXL),IFEX(MXL)
      COMMON /CURVE/CMIN(MXL),CMAX(MXL),NCOEFF(MXL),C(MXL,11)
      COMMON /BINS/NBTOT,SWL(100),SWIL(100),SWS(100),SWIS(100),W(100),
     *             R(100),ALAM(100),NVAL(100),AIS(100),AIL(100)
      DATA GR0/'................................'/
      DATA GR1/'.                              .'/
C
C OUTPUT STATISTICS
C -----------------
C
      CALL HEADXX(6,'NORMALISATION ANALYSIS$','$','1',2)
      WRITE (6,6010) GR0
      RMAX = 1E-10
      RMIN = 1E+10
      DO 20 N = 1 , NBTOT
         CV = CURVAL(ALAM(N))
         IF (CV.GT.RMAX) RMAX = CV
         IF (R(N).GT.RMAX) RMAX = R(N)
         IF (CV.LT.RMIN) RMIN = CV
         IF (R(N).LT.RMIN) RMIN = R(N)
 20   CONTINUE
      DO 30 NN = 1 , NRANGE
         IF (CURVAL(BMAXL(NN)).GT.RMAX) RMAX = CURVAL(BMAXL(NN))
         IF (CURVAL(BMINL(NN)).GT.RMAX) RMAX = CURVAL(BMINL(NN))
         IF (CURVAL(BMAXL(NN)).LT.RMIN) RMIN = CURVAL(BMAXL(NN))
         IF (CURVAL(BMINL(NN)).LT.RMIN) RMIN = CURVAL(BMINL(NN))
 30   CONTINUE
      GRLIN = GR1
      CV = CURVAL(BMINL(1))
      EXT1 = CV
      IP = 15
      IF (RMAX.GT.RMIN) IP = NINT(29*(CV-RMIN)/(RMAX-RMIN)) + 2
      GRLIN(IP:IP) = 'E'
      WRITE (6,6040) GRLIN
      DO 60 NN = 1 , NRANGE
         DO 50 NA = 1 , NBIN(NN)
            N = NP(NN) + NA
            STAR = ' '
            X = 0.0
            IF (ALAM(N).GE.CMIN(NN) .AND. ALAM(N).LE.CMAX(NN)) THEN
               X = CURVAL(ALAM(N))
            END IF
            IF (W(N).LE.0.0) STAR = '*'
            GRLIN = GR1
            IP1 = 15
            IP2 = 15
            IF (RMAX.GT.RMIN) THEN
               IP1 = NINT(29*(R(N)-RMIN)/(RMAX-RMIN)) + 2
               CV = CURVAL(ALAM(N))
               IP2 = NINT(29*(CV-RMIN)/(RMAX-RMIN)) + 2
            END IF
            GRLIN(IP1:IP1) = 'X'
            GRLIN(IP2:IP2) = '+'
            IF (IP1.EQ.IP2) GRLIN(IP1:IP1) = '*'
            
            WRITE (6,6020) N , STAR , ALAM(N) , NVAL(N) , AIS(N) , 
     +                     AIL(N) , R(N) , W(N) , X , GRLIN
            IF (NA.NE.NBIN(NN)) WRITE (6,6040) GR1
 50      CONTINUE
         GRLIN = GR1
         CV = CURVAL(BMAXL(NN))
         EXT2 = CV
         IP = 15
         IF (RMAX.GT.RMIN) IP = NINT(29*(CV-RMIN)/(RMAX-RMIN)) + 2
         GRLIN(IP:IP) = 'E'
         WRITE (6,6040) GRLIN
         WRITE (6,6040) GR0
         WRITE (6,6070) EXT1 , EXT2
         IF (NN.NE.NRANGE) THEN
            WRITE (6,6060)
            WRITE (6,6040) GR0
            GRLIN = GR1
            CV = CURVAL(BMINL(NN+1))
            EXT1 = CV
            IP = 15
            IF (RMAX.GT.RMIN) IP = NINT(29*(CV-RMIN)/(RMAX-RMIN)) + 2
            GRLIN(IP:IP) = 'E'
            WRITE (6,6040) GRLIN
         END IF
 60   CONTINUE
      WRITE (6,6030)
      RETURN
C
C FORMAT STATEMENTS
C -----------------
C
 6010 FORMAT (//,'  BIN  AVERAGE  NO. OF     MEAN        MEAN    ',
     +        'WEIGHTED    MEAN      CURVE FITTED',/,
     +        '  NO.   LAMBDA  REFLNS',
     +  '      IS          IL    MEAN IL/IS  WEIGHT   VALUE(MEAN IL/IS)'
     +  ,6X,A)
 6020 FORMAT (' ',I4,A1,F7.3,I7,F12.1,F12.1,F10.4,F10.1,F14.4,12X,A)
 6030 FORMAT (//,' KEY:  *  = BIN NOT USED IN CURVE FITTING',/,
     +        '       IL = LAUE INTENSITY',/,
     +        '       IS = STANDARD (MONOCHROMATIC) INTENSITY',/,/,
     +        '        E = EXTRAPOLATED POINT ON GRAPH',/,
     +        '        X = POINT USED FOR CURVE FITTING',/,
     +        '        + = CALCULATED POINT FROM FITTED CURVE',/,
     +        '        * = POINTS OVERLAP (OR X AND + OVERPRINTED)',/)
 6040 FORMAT (' ',89X,A)
 6050 FORMAT ('+',89X,A)
 6060 FORMAT (//)
 6070 FORMAT (/,' CURVE FITTED MEAN IL/IS AT FIRST POINT =',F8.4,/,
     +        ' CURVE FITTED MEAN IL/IS AT FINAL POINT =',F8.4)
      END
**==rscalb.f  
C
C SUBROUTINE 'RSCALB'
C ===================
C
C THIS SUBROUTINE PERFORMS 'RSCALE' ANALYSIS AND B FACTOR CALCULATION
C
      SUBROUTINE RSCALB(IA,IMAX,NREF)
C
C PARAMETERS
C
C         IA (I)   ARRAY TO HOLD LAUE DATA (5 WORDS/REFLN)
C       IMAX (I)   NO. OF WORDS IN 'IA'
C       NREF (I)   NUMBER OF LAUE REFLECTIONS
C
C
      PARAMETER(MXR=25000,MXL=8)
C
C SPECIFICATION STATEMENTS
C ------------------------
C
      LOGICAL EOF
      LOGICAL NORM,STFIL,FILE1,FILE2,FILAGR,END,AFSCAL,SHLX1,SHLX2,SHLX3       
      DIMENSION IA(IMAX)
      DIMENSION ADATA(10),IDATA(3),IDXOLD(3),IDXL(3)
      COMMON /FLAGS/NORM,STFIL,FILE1,FILE2,FILAGR,AFSCAL,SHLX1,SHLX2,
     +              SHLX3
      LOGICAL UPDATB,NRMCAL,ABSCAL
      COMMON /SCFLAG/UPDATB,NRMCAL,ABSCAL
      COMMON /CDATA/AMINL,AMAXL,RES1,RES2,TH1,TH2,BFAC,NRANGE,IBATCH,
     *              NITER,IABS,LASTIT,KBATCH
      COMMON /WGHTS/IWTS,AWTS,IWTL,AWTL,IWMEAN
      COMMON /RANGES/BMINL(MXL),BMAXL(MXL),BINS(MXL),NBIN(MXL),
     *               NFIT(MXL),NP(MXL),IFEX(MXL)
      COMMON /BINS/NBTOT,SWL(100),SWIL(100),SWS(100),SWIS(100),W(100),
     *             R(100),ALAM(100),NVAL(100),AIS(100),AIL(100)
      COMMON /RFACTS/RNS,RDS,NRS,RN(3),RD(3),NRF(3),RNB(100),RDB(100),
     *               NRB(100)
      COMMON /RFACT3/STLMIN,STLMAX,RNSTL(10),RDSTL(10),NRSTL(10),
     *SISSTL(10),SILSTL(10)
      COMMON /RFACT4/FSQMIN,FSQMAX,RNFSQ(10),RDFSQ(10),NRFSQ(10),
     *SISFSQ(10),SILFSQ(10),ALOGI(10)
      COMMON /SCAL/SCALEI,SCAL2,SCALIN,SIGREJ
      COMMON /NWR/NWR
      LOGICAL LAMOPT
      COMMON /LAMO/ LAMOPT, DLAM
      LOGICAL CENTYP
      COMMON /CENT/ CENTYP
      EQUIVALENCE (IAA,AAA)
      DATA NTWO16/65536/
C
C INITIALISATIONS
C ---------------
C
      DO 20 N = 1 , NBTOT
         RNB(N) = 0.0
         RDB(N) = 0.0
         NRB(N) = 0
 20   CONTINUE
      END = .TRUE.
      CALL LRREWD(1)
      CALL LRREFF(1,RESOL,ADATA,EOF)
      END = .FALSE.
      IDATA(1) = NINT(ADATA(1))
      IDATA(2) = NINT(ADATA(2))
      IDATA(3) = NINT(ADATA(3))
      ALMBDA = 0.0
      KCS = 0
      IF (LAMOPT) THEN
          ALMBDA = ADATA(6)
          IF (CENTYP) KCS = NINT(ADATA(7))
      ELSE
          IF (CENTYP) KCS = NINT(ADATA(6))
      ENDIF
 30   DO 40 I = 1 , 3
         IDXOLD(I) = 0
 40   CONTINUE
      BLMOLD = -1.0
      ICSOLD = -100
      DO 50 I = 1 , 10
         RNSTL(I) = 0.0
         RDSTL(I) = 0.0
         NRSTL(I) = 0
         SISSTL(I) = 0.0
         SILSTL(I) = 0.0
 50   CONTINUE
      RNS = 0.0
      RDS = 0.0
      NRS = 0
      ICOUNT = 0
      NR = 0
C
C STEP THROUGH LAUE REFLECTIONS
C -----------------------------
C
 60   NR = NR + 1
      IF (NR.LE.NREF) THEN
         IB = NWR*(NR-1)
         CALL GETHKL(IA(IB+1),IDXL,ICS)
         BLMBDA = IA(IB+2)/10000.0
         IF (ICMPR(IDXL,IDXOLD,BLMOLD,BLMBDA,KCS,ICSOLD,LAMOPT,DLAM)
     +      .EQ.0) THEN       
            ICOUNT = ICOUNT + 1
            GO TO 60
         END IF
      END IF
C
C NEW REFLECTION FOUND, PROCESS THE DATA FOR THE PREVIOUS ONE
C -----------------------------------------------------------
C
      IF (ICOUNT.EQ.0) GO TO 110
C
C LOOP THROUGH THE MEASUREMENTS FOR THIS REFLECTION
C
      SIGWI = 0.0
      SIGW = 0.0
      DO 70 IC = 1 , ICOUNT
         IB = NWR*(NR-ICOUNT+IC-2)
         ALM = IA(IB+2)/10000.0
         TH = IA(IB+5)/10000.0
         IAA = IA(IB+3)
         SC = SCFAC(IA(IB+1),NRMCAL,ABSCAL)
         AI = AAA*SC
         KSYM = IA(IB+4)/NTWO16
         SI = (IA(IB+4)-NTWO16*KSYM)*SC
         WI = 1/(SI*SI)
         IF (IWMEAN.NE.0) WI = 1.0
         SIGWI = SIGWI + WI*AI
         SIGW = SIGW + WI
 70   CONTINUE
      IF (SIGW.LE.0.0) GO TO 110
      FF = SCALEI*SIGWI/SIGW
      SS = SCALEI*SQRT(1/SIGW)
      IF (IWMEAN.NE.0) SIG = 1.0
      IF (END) GO TO 110
 80   IF (ICMPR(IDXOLD,IDATA,ALMBDA,BLMOLD,KCS,ICSOLD,LAMOPT,DLAM)
     +   .LT.0) GO TO 110       
      IF (ICMPR(IDXOLD,IDATA,ALMBDA,BLMOLD,KCS,ICSOLD,LAMOPT,DLAM)
     +   .EQ.0) THEN
C
C MATCH FOUND, ADD CONTRIBUTIONS TO THE R-FACTOR
C
         IF (ADATA(4).NE.0.0 .AND. ADATA(5).NE.0.0) THEN
            FSQ = ADATA(4)*ADATA(4)
            RNS = RNS + ABS(FSQ-FF)
            RDS = RDS + (FSQ+FF)
            NRS = NRS + 1
            DO 90 IC = 1 , ICOUNT
               IB = NWR*(NR-ICOUNT+IC-2)
               ALM = IA(IB+2)/10000.0
               IBIN = LAMBIN(ALM)
               IAA = IA(IB+3)
               AI = AAA*SCALEI*SCFAC(IA(IB+1),NRMCAL,ABSCAL)
               ADEL = ABS(FSQ-AI)
               RNB(IBIN) = RNB(IBIN) + ADEL
               RDB(IBIN) = RDB(IBIN) + FSQ + AI
               NRB(IBIN) = NRB(IBIN) + 1
               TH = IA(IB+5)/10000.0
               STL = (SIN(TH)/ALM)**2
               K = IRANGE(STL,STLMIN,STLMAX,10)
               RNSTL(K) = RNSTL(K) + ADEL
               RDSTL(K) = RDSTL(K) + FSQ + AI
               NRSTL(K) = NRSTL(K) + 1
               SILSTL(K) = SILSTL(K) + AI
               SISSTL(K) = SISSTL(K) + FSQ
 90         CONTINUE
         END IF
         GO TO 110
      ELSE
         CALL LRREFF(1,RESOL,ADATA,EOF)
         IF (EOF) GO TO 100
         IDATA(1) = NINT(ADATA(1))
         IDATA(2) = NINT(ADATA(2))
         IDATA(3) = NINT(ADATA(3))
         ALMBDA = 0.0
         KCS = 0
         IF (LAMOPT) THEN
             ALMBDA = ADATA(6)
             IF (CENTYP) KCS = NINT(ADATA(7))
         ELSE
             IF (CENTYP) KCS = NINT(ADATA(6))
         ENDIF
         GO TO 80
      END IF
 100  END = .TRUE.
 110  ICOUNT = 1
      IF (NR.GT.NREF) THEN
C
C END OF DATA REACHED
C --------------------
C
         RSCALE = -1.0
         DO 130 NN = 1 , NRANGE
            DO 120 NA = 1 , NBIN(NN)
               N = NP(NN) + NA
               AIL(N) = 0.0
               NVAL(N) = 0
               ALAM(N) = BMINL(NN) + (NA-0.5)*(BMAXL(NN)-BMINL(NN))
     +                   /NBIN(NN)
 120        CONTINUE
 130     CONTINUE
         NR = 0
 140     NR = NR + 1
         IF (NR.GT.NREF) THEN
C
C PRINT INTENSITY ANALYSIS
C ------------------------
C
            CALL HEADXX(6,'ANALYSIS OF SCALED LAUE INTENSITIES$','$',
     +                  '1',2)
            DO 150 N = 1 , NBTOT
               IF (NVAL(N).GT.0) AIL(N) = AIL(N)/NVAL(N)
 150        CONTINUE
            WRITE (6,6010)
            DO 160 N = 1 , NBTOT
               RFAC = -1.0
               IF (NRB(N).GT.0) RFAC = RNB(N)/RDB(N)
               WRITE (6,6020) N , ALAM(N) , NVAL(N) , AIS(N) , AIL(N) , 
     +                        RFAC , NRB(N)
 160        CONTINUE
C
C PRINT RSCALE
C ------------
            IF (RDS.GT.0) RSCALE = RNS/RDS
            WRITE (6,6030)
            WRITE (6,6050)
            WRITE (6,6060) RSCALE
            WRITE (6,6040) NRS
C
C SIN**2THETA/LAMBDA**2 ANALYSIS
C ------------------------------
C
            WRITE (6,6070) STLMIN , STLMAX
            WRITE (6,6080)
            SX = 0
            SY = 0
            SXX = 0
            SXY = 0
            DO 170 K = 1 , 10
               RFAC = -1.0
               RATIO = 0.0
               ALOGI(K) = 0.0
               IF (NRSTL(K).NE.0) THEN
                  RFAC = RNSTL(K)/RDSTL(K)
                  RATIO = SILSTL(K)/SISSTL(K)
                  ALOGI(K) = ALOG(RATIO)
               END IF
               AV = STLMIN + (K-0.5)*(STLMAX-STLMIN)/10
               WRITE (6,6090) K , AV , NRSTL(K) , RATIO , RFAC , 
     +                        ALOGI(K)
               SX = SX + AV
               SY = SY + ALOGI(K)
               SXY = SXY + AV*ALOGI(K)
               SXX = SXX + AV*AV
 170        CONTINUE
            BFCT = (10*SXY-SX*SY)/(10*SXX-SX*SX)
            WRITE (6,6100) BFCT
C
C UPDATE B FACTOR IF REQUIRED
C
            IF (UPDATB) THEN
               BFAC = BFAC + BFCT
               WRITE (6,6110) BFAC
            END IF
            RETURN
         ELSE
            I = NWR*(NR-1)
            BLAM = IA(I+2)/10000.0
            IAA = IA(I+3)
            FSQL = AAA
            IBIN = LAMBIN(BLAM)
            IF (IBIN.NE.0) THEN
               AIL(IBIN) = AIL(IBIN) + FSQL*SCFAC(IA(I+1),NRMCAL,ABSCAL)
               NVAL(IBIN) = NVAL(IBIN) + 1
            END IF
            GO TO 140
         END IF
      ELSE
         IDXOLD(1) = IDXL(1)
         IDXOLD(2) = IDXL(2)
         IDXOLD(3) = IDXL(3)
         BLMOLD = BLMBDA
         ICSOLD = ICS
         GO TO 60
      END IF
C
C FORMAT STATEMENTS
C -----------------
C
 6010 FORMAT (//,'  BIN  AVERAGE NO. LAUE    MEAN        MEAN',
     +        '      R-FACTOR',/,
     +        '  NO.   LAMBDA  REFLNS      IS          IL',/)
 6020 FORMAT (' ',I4,F8.3,I7,F13.1,F12.1,F8.3,' (',I5,')')
 6030 FORMAT (////)
 6040 FORMAT (/,' N =',I6)
 6050 FORMAT (///,' SCALING R-FACTOR = SUM(ABS(IL-IS)/SUM(IL+IS)',//,
     +        '                    (SUM FOR EACH MERGED REFLECTION',/,
     +        '                     PRESENT IN BOTH DATA SETS)')
 6060 FORMAT (//,' RSCALE =',F6.3)
 6070 FORMAT (///,' ANALYSIS OF DATA IN 10 EQUAL BINS OF ',
     +        '(SIN(THETA)/LAMBDA)**2 FROM',F7.4,' TO',F7.4)
 6080 FORMAT (//,
     +        '  BIN  AVERAGE NO. LAUE     MEAN  R-FACTOR  LOG(IL/IS)'/,
     +        '  NO.    STL2   REFLNS     IL/IS',/)
 6090 FORMAT (' ',I4,F8.3,I7,F12.3,F8.3,F11.3)
 6100 FORMAT (//,' B-FACTOR DIFFERENCE = ',F6.2)
 6110 FORMAT (//,' UPDATED B-FACTOR =',F6.2)
C
      END
**==opscal.f  
C
C SUBROUTINE 'OPSCAL'
C ===================
C
C THIS SUBROUTINE OUTPUTS THE SCALED DATA AND COLLECTS DATA FOR THE
C R-FACTOR ANALYSES
C
      SUBROUTINE OPSCAL(IA,IMAX,NREF)
C
C PARAMETERS
C
C         IA (I)   ARRAY TO HOLD LAUE DATA (5 WORDS/REFLN)
C       IMAX (I)   NO. OF WORDS IN 'IA'
C       NREF (I)   NUMBER OF LAUE REFLECTIONS
C
C
      PARAMETER(MXR=25000,MXL=8)
C
C SPECIFICATION STATEMENTS
C ------------------------
C
      LOGICAL EOF
      LOGICAL NORM,STFIL,FILE1,FILE2,FILAGR,END,PAIR,RF(3),AFSCAL
      LOGICAL SHLX1,SHLX2,SHLX3
      DIMENSION IA(IMAX)
      DIMENSION IDATA(3)
      DIMENSION ADATA(10),BDATA(10),CDATA(10),IDXOLD(3),IDXL(3)
      DIMENSION SIGWI(3),SIGW(3),AVEI(3),DEL(3),HKL(3)
      COMMON /FLAGS/NORM,STFIL,FILE1,FILE2,FILAGR,AFSCAL,SHLX1,SHLX2,
     +              SHLX3
      LOGICAL UPDATB,NRMCAL,ABSCAL
      COMMON /SCFLAG/UPDATB,NRMCAL,ABSCAL
      COMMON /CDATA/AMINL,AMAXL,RES1,RES2,TH1,TH2,BFAC,NRANGE,IBATCH,
     *              NITER,IABS,LASTIT,KBATCH
      COMMON /WGHTS/IWTS,AWTS,IWTL,AWTL,IWMEAN
      COMMON /RANGES/BMINL(MXL),BMAXL(MXL),BINS(MXL),NBIN(MXL),
     *               NFIT(MXL),NP(MXL),IFEX(MXL)
      COMMON /BINS/NBTOT,SWL(100),SWIL(100),SWS(100),SWIS(100),W(100),
     *             R(100),ALAM(100),NVAL(100),AIS(100),AIL(100)
      COMMON /RFACTS/RNS,RDS,NRS,RN(3),RD(3),NRF(3),RNB(100),RDB(100),
     *               NRB(100)
      COMMON /RFACT2/THEMIN,THEMAX,RNTHE(10),RDTHE(10),NRTHE(10),
     *SISTHE(10),SILTHE(10)
      COMMON /RFACT3/STLMIN,STLMAX,RNSTL(10),RDSTL(10),NRSTL(10),
     *SISSTL(10),SILSTL(10)
      COMMON /RFACT4/FSQMIN,FSQMAX,RNFSQ(10),RDFSQ(10),NRFSQ(10),
     *SISFSQ(10),SILFSQ(10),ALOGI(10)
      COMMON /SCAL/SCALEI,SCAL2,SCALIN,SIGREJ
      COMMON /NWR/NWR
      LOGICAL LAMOPT
      COMMON /LAMO/ LAMOPT, DLAM
      LOGICAL CENTYP
      COMMON /CENT/ CENTYP
      EQUIVALENCE (IAA,AAA)
      DATA NTWO16/65536/
C
C INITIALISATIONS
C ---------------
C
      DO 20 N = 1 , NBTOT
         RNB(N) = 0.0
         RDB(N) = 0.0
         NRB(N) = 0
 20   CONTINUE
      NC1 = 0
      NC2 = 0
      NCLO = 5
      IF (STFIL) NCLO = 7
      IF (CENTYP) THEN
         NCLO = NCLO + 1
         NC1 = NCLO
      ENDIF
      IF (LAMOPT) THEN
         NCLO = NCLO + 1
         NC2 = NCLO
      ENDIF
      IF (STFIL) THEN
         END = .TRUE.
         CALL LRREWD(1)
         CALL LRREFF(1,RESOL,ADATA,EOF)
         IF (EOF) GO TO 30
         IDATA(1) = NINT(ADATA(1))
         IDATA(2) = NINT(ADATA(2))
         IDATA(3) = NINT(ADATA(3))
         ALMBDA = 0.0
         KCS = 0
         IF (LAMOPT) THEN
             ALMBDA = ADATA(6)
             IF (CENTYP) KCS = NINT(ADATA(7))
         ELSE
             IF (CENTYP) KCS = NINT(ADATA(6))
         ENDIF
         END = .FALSE.
      END IF
 30   DO 40 I = 1 , 3
         IDXOLD(I) = 0
         RN(I) = 0.0
         RD(I) = 0.0
         NRF(I) = 0
 40   CONTINUE
      BLMOLD = -1.0
      ICSOLD = -100
      DO 50 I = 1 , 10
         RNTHE(I) = 0.0
         RDTHE(I) = 0.0
         NRTHE(I) = 0
         SISTHE(I) = 0.0
         SILTHE(I) = 0.0
         RNSTL(I) = 0.0
         RDSTL(I) = 0.0
         NRSTL(I) = 0
         SISSTL(I) = 0.0
         SILSTL(I) = 0.0
         RNFSQ(I) = 0.0
         RDFSQ(I) = 0.0
         NRFSQ(I) = 0
         SISFSQ(I) = 0.0
         SILFSQ(I) = 0.0
 50   CONTINUE
      RNS = 0.0
      RDS = 0.0
      NRS = 0
      CALL SYMINV
      ICOUNT = 0
      NR = 0
C
C STEP THROUGH LAUE REFLECTIONS
C -----------------------------
C
 60   NR = NR + 1
      IF (NR.LE.NREF) THEN
         IB = NWR*(NR-1)
         CALL GETHKL(IA(IB+1),IDXL,ICS)
         BLMBDA = IA(IB+2)/10000.0
         IF (ICMPR(IDXL,IDXOLD,BLMOLD,BLMBDA,ICSOLD,ICS,LAMOPT,DLAM)
     +      .EQ.0) THEN
            ICOUNT = ICOUNT + 1
            GO TO 60
         END IF
      END IF
C
C NEW REFLECTION FOUND, PROCESS THE DATA FOR THE PREVIOUS ONE
C -----------------------------------------------------------
C
      IF (ICOUNT.EQ.0) THEN
         GO TO 140
      ELSE
C
C LOOP THROUGH THE MEASUREMENTS FOR THIS REFLECTION
C
         DO 100 IC = 1 , ICOUNT
            PAIR = .FALSE.
            IB = NWR*(NR-ICOUNT+IC-2)
            ALM = IA(IB+2)/10000.0
            TH = IA(IB+5)/10000.0
            IF (IC.EQ.1) THEN
               STH = SIN(TH)
               STL4 = 4*STH*STH/(ALM*ALM)
            END IF
            SC = SCFAC(IA(IB+1),NRMCAL,ABSCAL)
            IAA = IA(IB+3)
            AI = AAA*SC
            KSYM = IA(IB+4)/NTWO16
            SI = (IA(IB+4)-NTWO16*KSYM)*SC
            WI = 1/(SI*SI)
            IF (IWMEAN.NE.0) WI = 1.0
            DO 70 I = 1 , 3
               SIGWI(I) = 0.0
               SIGW(I) = 0.0
               RF(I) = .FALSE.
 70         CONTINUE
C
C FIND THE MEAN INTENSITIES FOR THE REFLECTION REQUIRED FOR CALCULATING
C THE STATISTICS. 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 MEANS OF MEASUREMENTS WITH THE
C SAME SIGN AND MEASURED WITHIN 0.1 ANGSTROMS OF THE SAME WAVELENGTH.
C
            DO 80 JC = 1 , ICOUNT
               IF (IC.NE.JC) RF(1) = .TRUE.
               JB = NWR*(NR-ICOUNT+JC-2)
               AALM = IA(JB+2)/10000.0
               TTH = IA(IB+5)/10000.0
               SSC = SCFAC(IA(JB+1),NRMCAL,ABSCAL)
               IAA = IA(JB+3)
               AAI = AAA*SSC
               KKSYM = IA(JB+4)/NTWO16
               SSI = (IA(JB+4)-NTWO16*KKSYM)*SSC
               WWI = 1/(SSI*SSI)
               IF (IWMEAN.NE.0) WWI = 1.0
               SIGWI(1) = SIGWI(1) + WWI*AAI
               SIGW(1) = SIGW(1) + WWI
               IF (FILE2) THEN
                  IF (MOD(KSYM,2).EQ.1 .AND. MOD(KKSYM,2).EQ.0 .AND. 
     +                ABS(ALM-AALM).LT.0.1) THEN
                     PAIR = .TRUE.
                     CDATA(4) = 0.5*(ALM+AALM)
                     CDATA(5) = 0.5*(TH+TTH)
                     F1 = SQRT(AI*SCALEI)
                     S1 = SQRT(SCALEI*(AI+SI)) - SQRT(SCALEI*AI)
                     F2 = SQRT(AAI*SCALEI)
                     S2 = SQRT(SCALEI*(AAI+SSI)) - SQRT(SCALEI*AAI)
                     CDATA(6) = 0.5*(F1+F2)
                     CDATA(7) = 0.5*SQRT(S1*S1+S2*S2)
                     CDATA(8) = F1-F2
                     CDATA(9) = SQRT(S1*S1+S2*S2)
                  END IF
               END IF
               IF (MOD(KKSYM,2).EQ.MOD(KSYM,2)) THEN
C
C SAME SIGN
C
                  SIGWI(2) = SIGWI(2) + WWI*AAI
                  SIGW(2) = SIGW(2) + WWI
                  IF (IC.NE.JC) RF(2) = .TRUE.
                  IF (ABS(AALM-ALM).LE.0.1) THEN
C
C SAME SIGN AND LAMBDAS WITHIN 0.1 ANGSTROM
C
                     SIGWI(3) = SIGWI(3) + WWI*AAI
                     SIGW(3) = SIGW(3) + WWI
                     IF (IC.NE.JC) RF(3) = .TRUE.
                  END IF
               END IF
 80         CONTINUE
C
C ACCUMULATE VALUES FOR R-FACTORS
C
            DO 90 I = 1 , 3
               AVEI(I) = SIGWI(I)/SIGW(I)
               DEL(I) = ABS(AI-AVEI(I))
               IF (RF(I)) THEN
                  RN(I) = RN(I) + DEL(I)
                  RD(I) = RD(I) + AVEI(I)
                  NRF(I) = NRF(I) + 1
               END IF
 90         CONTINUE
            CALL DIAG2(IDXOLD,SC,STL4,ALM,AI,SI,AVEI,KSYM,F1,F2,PAIR)
C
C WRITE REFLECTION TO OUTPUT FILE 2 IF REQUIRED (ANOM DIFF FILE)
C
            IF (FILE2) THEN
               IF (PAIR) THEN
                  CALL DECSYM(IDXOLD,HKL,KSYM)
                  CDATA(1) = HKL(1)
                  CDATA(2) = HKL(2)
                  CDATA(3) = HKL(3)
                  CALL LWREFL(2,CDATA)
               END IF
            END IF
C
C WRITE REFLECTION DATA TO OUTPUT FILE 2 IF REQUIRED (ROTAVATA/AGROVATA)
C
            IF (FILAGR) THEN
               CDATA(1) = IDXOLD(1)
               CDATA(2) = IDXOLD(2)
               CDATA(3) = IDXOLD(3)
               CDATA(4) = KSYM
               CDATA(5) = IBATCH
               CDATA(6) = SCALEI*AI/SCAL2
               CDATA(7) = SCALEI*SI/SCAL2
               CALL LWREFL(2,CDATA)
            END IF
            IF (SHLX2) THEN
               FF = SCALEI*AI
               SS = SCALEI*SI
               F = SQRT(FF)
               SIGF = SQRT(FF+SS) - SQRT(FF)
               CALL DECSYM(IDXOLD,HKL,KSYM)
               WRITE (2,6030) NINT(HKL(1)), NINT(HKL(2)), NINT(HKL(3)),
     +                        F, SIGF, 
     +                        ALM
            END IF
            IF (SHLX3) THEN
               FF = SCALEI*AI
               SS = SCALEI*SI
               CALL DECSYM(IDXOLD,HKL,KSYM)
               WRITE (2,6040) NINT(HKL(1)), NINT(HKL(2)), NINT(HKL(3)),
     +                        FF/100.0, SS/100.0, KBATCH,
     +                        ALM
            END IF
C
C END OF LOOP FOR THIS REFLECTION
C
 100     CONTINUE
C
C WRITE REFLECTION TO OUTPUT FILE 1
C ---------------------------------
C
         BDATA(4) = 0.0
         BDATA(5) = 0.0
         IF (SIGW(1).LE.0.0) GO TO 140
         FF = SCALEI*SIGWI(1)/SIGW(1)
         SS = SCALEI*SQRT(1/SIGW(1))
         F = SQRT(FF)
         SIG = SQRT(FF+SS) - SQRT(FF)
         IF (IWMEAN.NE.0) SIG = 1.0
         BDATA(4) = F
         BDATA(5) = SIG
         BDATA(1) = IDXOLD(1)
         BDATA(2) = IDXOLD(2)
         BDATA(3) = IDXOLD(3)
C
C COMPARE REFLECTION WITH STANDARD IF REQUIRED
C ---------------------------------------------
C
         IF (.NOT.STFIL) THEN
            IF (NC1.GT.0) BDATA(NC1) = ICSOLD
            IF (NC2.GT.0) BDATA(NC2) = BLMOLD
            IF (ICOUNT.GT.0 .AND. FILE1) CALL LWREFL(1,BDATA)
            IF (ICOUNT.GT.0 .AND. SHLX1) THEN
               IF (CENTYP.OR.LAMOPT) THEN
                  IF (LAMOPT) THEN
                     IF (CENTYP) THEN
                        WRITE (2,6050) 
     +                     IDXOLD(1), IDXOLD(2),IDXOLD(3),
     +                     F, SIG, BDATA(NC2), NINT(BDATA(NC1))
                      ELSE
                        WRITE (2,6050) 
     +                     IDXOLD(1), IDXOLD(2),IDXOLD(3),
     +                     F, SIG, BDATA(NC2)
                      END IF
                   ELSE
                      WRITE (2,6060) 
     +                  IDXOLD(1), IDXOLD(2),IDXOLD(3),
     +                  F, SIG, NINT(BDATA(NC1))
                   END IF
               ELSE
                  WRITE (2,6030)
     +            IDXOLD(1), IDXOLD(2), IDXOLD(3), F, SIG
               END IF
            END IF
            GO TO 140
         ELSE
            BDATA(7) = BDATA(5)
            BDATA(6) = BDATA(4)
            BDATA(5) = 0.0
            BDATA(4) = 0.0
            IF (NC1.GT.0) BDATA(NC1) = KCS
            IF (NC2.GT.0) BDATA(NC2) = ALMBDA
            IF (END) GO TO 140
 110        IF (ICMPR(IDXOLD,IDATA,ALMBDA,BLMOLD,KCS,ICSOLD,
     +          LAMOPT,DLAM).LT.0)
     +         GO TO 140
            IF (ICMPR(IDXOLD,IDATA,ALMBDA,BLMOLD,KCS,ICSOLD,
     +          LAMOPT,DLAM).EQ.0) THEN
C
C MATCH FOUND, ADD CONTRIBUTIONS TO THE R-FACTOR
C
               BDATA(4) = ADATA(4)
               BDATA(5) = ADATA(5)
               IF (ADATA(4).NE.0.0 .AND. ADATA(5).NE.0.0) THEN
                  FSQ = ADATA(4)*ADATA(4)
                  RNS = RNS + ABS(FSQ-FF)
                  RDS = RDS + (FSQ+FF)
                  NRS = NRS + 1
                  DO 120 IC = 1 , ICOUNT
                     IB = NWR*(NR-ICOUNT+IC-2)
                     ALM = IA(IB+2)/10000.0
                     IBIN = LAMBIN(ALM)
                     IF (IBIN.NE.0) THEN
                        IAA = IA(IB+3)
                        AI = AAA*SCALEI*SCFAC(IA(IB+1),NRMCAL,ABSCAL)
                        ADEL = ABS(FSQ-AI)
                        RNB(IBIN) = RNB(IBIN) + ADEL
                        RDB(IBIN) = RDB(IBIN) + FSQ + AI
                        NRB(IBIN) = NRB(IBIN) + 1
                        TH = IA(IB+5)/10000.0
                        STL = (SIN(TH)/ALM)**2
                        K = IRANGE(TH,THEMIN,THEMAX,10)
                        RNTHE(K) = RNTHE(K) + ADEL
                        RDTHE(K) = RDTHE(K) + FSQ + AI
                        NRTHE(K) = NRTHE(K) + 1
                        SILTHE(K) = SILTHE(K) + AI
                        SISTHE(K) = SISTHE(K) + FSQ
                        K = IRANGE(STL,STLMIN,STLMAX,10)
                        RNSTL(K) = RNSTL(K) + ADEL
                        RDSTL(K) = RDSTL(K) + FSQ + AI
                        NRSTL(K) = NRSTL(K) + 1
                        SILSTL(K) = SILSTL(K) + AI
                        SISSTL(K) = SISSTL(K) + FSQ
                        K = IRANGE(FSQ,FSQMIN,FSQMAX,10)
                        RNFSQ(K) = RNFSQ(K) + ADEL
                        RDFSQ(K) = RDFSQ(K) + FSQ + AI
                        NRFSQ(K) = NRFSQ(K) + 1
                        SILFSQ(K) = SILFSQ(K) + AI
                        SISFSQ(K) = SISFSQ(K) + FSQ
                     END IF
 120              CONTINUE
                  IF (ICOUNT.GT.0 .AND. FILE1) CALL LWREFL(1,BDATA)
                  IF (ICOUNT.GT.0 .AND. SHLX1) THEN
                     IF (CENTYP.OR.LAMOPT) THEN
                        IF (LAMOPT) THEN
                           IF (CENTYP) THEN
                              WRITE (2,6050) 
     +                           IDXOLD(1), IDXOLD(2),IDXOLD(3),
     +                           F, SIG, BDATA(NC2), NINT(BDATA(NC1))
                            ELSE
                              WRITE (2,6050) 
     +                           IDXOLD(1), IDXOLD(2),IDXOLD(3),
     +                           F, SIG, BDATA(NC2)
                            END IF
                         ELSE
                            WRITE (2,6060) 
     +                        IDXOLD(1), IDXOLD(2),IDXOLD(3),
     +                        F, SIG, NINT(BDATA(NC1))
                         END IF
                     ELSE
                        WRITE (2,6030)
     +                  IDXOLD(1), IDXOLD(2), IDXOLD(3), F, SIG
                     END IF
                  END IF
               END IF
               GO TO 140
            ELSE
               CALL LRREFF(1,RESOL,ADATA,EOF)
               IF (EOF) GO TO 130
               IDATA(1) = NINT(ADATA(1))
               IDATA(2) = NINT(ADATA(2))
               IDATA(3) = NINT(ADATA(3))
               ALMBDA = 0.0
               KCS = 0
               IF (LAMOPT) THEN
                   ALMBDA = ADATA(6)
                   IF (CENTYP) KCS = NINT(ADATA(7))
               ELSE
                   IF (CENTYP) KCS = NINT(ADATA(6))
               ENDIF
               GO TO 110
            END IF
         END IF
      END IF
 130  END = .TRUE.
 140  ICOUNT = 1
      IF (NR.GT.NREF) THEN
C
C END OF DATA REACHED
C --------------------
C
         IF (FILE1) CALL LWCLOS(1,1)
         IF (FILE2 .OR. FILAGR) CALL LWCLOS(2,1)
         RETURN
      ELSE
         IDXOLD(1) = IDXL(1)
         IDXOLD(2) = IDXL(2)
         IDXOLD(3) = IDXL(3)
         BLMOLD = BLMBDA
         ICSOLD = ICS
         GO TO 60
      END IF
C
C ERROR CONDITIONS
C ----------------
C
 160  WRITE (6,6020) IDXOLD , CDATA(6) , CDATA(7)
      STOP
C
C FORMAT STATEMENTS
C -----------------
C
 6020 FORMAT (/,' **OUTPUT INTENSITY OR SIGI OUT OF RANGE** HKL =',3I5,
     +        ', I =',G12.2,', SIGI =',G12.2)
 6030 FORMAT (3I4,2F8.2,F8.4)
 6040 FORMAT (3I4,2F8.2,I4,F8.4)
 6050 FORMAT (3I4,2F8.2,F8.4,I8)
 6060 FORMAT (3I4,2F8.2,I8)
      END
**==table2.f  
C
C SUBROUTINE 'TABLE2'
C ===================
C
C PRINT INTENSITY AND R-FACTOR ANALYSES
C
      SUBROUTINE TABLE2(IA,IMAX,NREF)
C
C PARAMETERS
C
C         IA (I)   ARRAY TO HOLD LAUE DATA (5 WORDS/REFLN)
C       IMAX (I)   NO. OF WORDS IN 'IA'
C       NREF (I)   NUMBER OF LAUE REFLECTIONS
C
      PARAMETER(MXR=25000,MXL=8)
C
C SPECIFICATION STATEMENTS
C ------------------------
C
      DIMENSION IA(IMAX)
      DIMENSION RFACT(3)
      LOGICAL NORM,STFIL,FILE1,FILE2,FILAGR,AFSCAL,SHLX1,SHLX2,SHLX3
      COMMON /FLAGS/NORM,STFIL,FILE1,FILE2,FILAGR,AFSCAL,SHLX1,SHLX2,
     +              SHLX3
      LOGICAL UPDATB,NRMCAL,ABSCAL
      COMMON /SCFLAG/UPDATB,NRMCAL,ABSCAL
      COMMON /CDATA/AMINL,AMAXL,RES1,RES2,TH1,TH2,BFAC,NRANGE,IBATCH,
     *              NITER,IABS,LASTIT,KBATCH
      COMMON /WGHTS/IWTS,AWTS,IWTL,AWTL,IWMEAN
      COMMON /RANGES/BMINL(MXL),BMAXL(MXL),BINS(MXL),NBIN(MXL),
     *               NFIT(MXL),NP(MXL),IFEX(MXL)
      COMMON /CURVE/CMIN(MXL),CMAX(MXL),NCOEFF(MXL),C(MXL,11)
      COMMON /BINS/NBTOT,SWL(100),SWIL(100),SWS(100),SWIS(100),W(100),
     *             R(100),ALAM(100),NVAL(100),AIS(100),AIL(100)
      COMMON /RFACTS/RNS,RDS,NRS,RN(3),RD(3),NRF(3),RNB(100),RDB(100),
     *               NRB(100)
      COMMON /RFACT2/THEMIN,THEMAX,RNTHE(10),RDTHE(10),NRTHE(10),
     *SISTHE(10),SILTHE(10)
      COMMON /RFACT3/STLMIN,STLMAX,RNSTL(10),RDSTL(10),NRSTL(10),
     *SISSTL(10),SILSTL(10)
      COMMON /RFACT4/FSQMIN,FSQMAX,RNFSQ(10),RDFSQ(10),NRFSQ(10),
     *SISFSQ(10),SILFSQ(10),ALOGI(10)
      COMMON /NWR/NWR
      EQUIVALENCE (IAA,AAA)
C
C COLLECT SCALED LAUE DATA IN LAMBDA BINS
C ---------------------------------------
C
      DO 20 I = 1 , 3
         RFACT(I) = -1.0
 20   CONTINUE
      RSCALE = -1.0
      DO 40 NN = 1 , NRANGE
         DO 30 NA = 1 , NBIN(NN)
            N = NP(NN) + NA
            AIL(N) = 0.0
            NVAL(N) = 0
            ALAM(N) = BMINL(NN) + (NA-0.5)*(BMAXL(NN)-BMINL(NN))
     +                /NBIN(NN)
 30      CONTINUE
 40   CONTINUE
      NR = 0
 50   NR = NR + 1
      IF (NR.GT.NREF) THEN
C
C PRINT INTENSITY ANALYSIS
C ------------------------
C
         CALL HEADXX(6,'ANALYSIS OF SCALED LAUE INTENSITIES$','$','1',2)
         DO 60 N = 1 , NBTOT
            IF (NVAL(N).GT.0) AIL(N) = AIL(N)/NVAL(N)
 60      CONTINUE
         IF (STFIL) THEN
            WRITE (6,6030)
            DO 70 N = 1 , NBTOT
               RFAC = -1.0
               IF (NRB(N).GT.0) RFAC = RNB(N)/RDB(N)
               WRITE (6,6040) N , ALAM(N) , NVAL(N) , AIS(N) , AIL(N) , 
     +                        RFAC , NRB(N)
 70         CONTINUE
         ELSE
            WRITE (6,6010)
            DO 80 N = 1 , NBTOT
               WRITE (6,6020) N , ALAM(N) , NVAL(N) , AIL(N)
 80         CONTINUE
         END IF
C
C PRINT R-FACTORS
C ---------------
C
         DO 90 I = 1 , 3
            IF (RD(I).GT.0 .AND. NRF(I).GT.0) RFACT(I) = RN(I)/RD(I)
 90      CONTINUE
         IF (RDS.GT.0) RSCALE = RNS/RDS
         WRITE (6,6050)
         CALL HEADXX(6,'R-FACTOR ANALYSES$','$','0',2)
         WRITE (6,6060)
         WRITE (6,6070)
         WRITE (6,6080) RFACT(1)
         WRITE (6,6110) NRF(1)
         WRITE (6,6090) RFACT(2)
         WRITE (6,6110) NRF(2)
         WRITE (6,6100) RFACT(3)
         WRITE (6,6110) NRF(3)
         IF (.NOT.STFIL) RETURN
         WRITE (6,6120)
         WRITE (6,6130) RSCALE
         WRITE (6,6110) NRS
C
C THETA ANALYSIS
C --------------
C
         WRITE (6,6140) THEMIN , THEMAX
         WRITE (6,6150)
         DO 100 K = 1 , 10
            RFAC = -1.0
            RATIO = 0.0
            IF (NRTHE(K).NE.0) THEN
               RFAC = RNTHE(K)/RDTHE(K)
               RATIO = SILTHE(K)/SISTHE(K)
            END IF
            AV = THEMIN + (K-0.5)*(THEMAX-THEMIN)/10
            WRITE (6,6160) K , AV , NRTHE(K) , RATIO , RFAC
 100     CONTINUE
C
C SIN**2THETA/LAMBDA**2 ANALYSIS
C ------------------------------
C
         WRITE (6,6170) STLMIN , STLMAX
         WRITE (6,6180)
         SX = 0
         SY = 0
         SXX = 0
         SXY = 0
         DO 110 K = 1 , 10
            RFAC = -1.0
            RATIO = 0.0
            ALOGI(K) = 0.0
            IF (NRSTL(K).NE.0) THEN
               RFAC = RNSTL(K)/RDSTL(K)
               RATIO = SILSTL(K)/SISSTL(K)
               ALOGI(K) = ALOG(RATIO)
            END IF
            AV = STLMIN + (K-0.5)*(STLMAX-STLMIN)/10
            WRITE (6,6190) K , AV , NRSTL(K) , RATIO , RFAC , ALOGI(K)
            SX = SX + AV
            SY = SY + ALOGI(K)
            SXY = SXY + AV*ALOGI(K)
            SXX = SXX + AV*AV
 110     CONTINUE
         BFCT = (10*SXY-SX*SY)/(10*SXX-SX*SX)
         WRITE (6,6230) BFCT
C
C INTENSITY ANALYSIS
C ------------------
C
         WRITE (6,6200) FSQMIN , FSQMAX
         WRITE (6,6210)
         DO 120 K = 1 , 10
            RFAC = -1.0
            RATIO = 0.0
            IF (NRFSQ(K).NE.0) THEN
               RFAC = RNFSQ(K)/RDFSQ(K)
               RATIO = SILFSQ(K)/SISFSQ(K)
            END IF
            AV = FSQMIN + (K-0.5)*(FSQMAX-FSQMIN)/10
            WRITE (6,6220) K , AV , NRFSQ(K) , RATIO , RFAC
 120     CONTINUE
         RETURN
      ELSE
         I = NWR*(NR-1)
         BLAM = IA(I+2)/10000.0
         IAA = IA(I+3)
         FSQL = AAA
         SC = SCFAC(IA(I+1),NRMCAL,ABSCAL)
         IBIN = LAMBIN(BLAM)
         IF (IBIN.NE.0) THEN
            AIL(IBIN) = AIL(IBIN) + FSQL*SC
            NVAL(IBIN) = NVAL(IBIN) + 1
         END IF
         GO TO 50
      END IF
C
C FORMAT STATEMENTS
C -----------------
C
 6010 FORMAT (//,'  BIN  AVERAGE NO. LAUE    MEAN',/,
     +        '  NO.   LAMBDA  REFLNS      IL',/)
 6020 FORMAT (' ',I4,F8.3,I7,F13.1)
 6030 FORMAT (//,'  BIN  AVERAGE NO. LAUE    MEAN        MEAN',
     +        '      R-FACTOR',/,
     +        '  NO.   LAMBDA  REFLNS      IS          IL',/)
 6040 FORMAT (' ',I4,F8.3,I7,F13.1,F12.1,F8.3,' (',I5,')')
 6050 FORMAT (////)
 6060 FORMAT (/' R VALUES OF -1.0 INDICATE THAT NO DATA WERE PRESENT')
 6070 FORMAT (//,' LAUE MERGING R-FACTORS = SUM(ABS(I-IMEAN)/SUM(I)',//,
     +        '                          (SUM OVER ALL MEASUREMENTS)')
 6080 FORMAT (//,' RFACT1 =',F6.3,
     +        '     IMEAN FROM ALL MEASUREMENTS FOR',
     +        ' THE REFLECTION (I+ AND I-)')
 6090 FORMAT (//,' RFACT2 =',F6.3,'     IMEAN FROM ALL MEASUREMENTS OF',
     +        ' THE SAME SIGN (I+ OR I-)')
 6100 FORMAT (//,' RFACT3 =',F6.3,'     IMEAN FROM ALL MEASUREMENTS OF',
     +        ' THE SAME SIGN (I+ OR I-) AND WITH LAMBDA WITHIN 0.1A')
 6110 FORMAT (/,' N =',I6)
 6120 FORMAT (///,' SCALING R-FACTOR = SUM(ABS(IL-IS)/SUM(IL+IS)',//,
     +        '                    (SUM FOR EACH MERGED REFLECTION',/,
     +        '                     PRESENT IN BOTH DATA SETS)')
 6130 FORMAT (//,' RSCALE =',F6.3)
 6140 FORMAT (///,' ANALYSIS OF DATA IN 10 EQUAL BINS OF THETA FROM',
     +        ' THETA =',F6.3,' TO',F6.3)
 6150 FORMAT (//,'  BIN  AVERAGE NO. LAUE     MEAN  R-FACTOR',/,
     +        '  NO.   THETA   REFLNS     IL/IS',/)
 6160 FORMAT (' ',I4,F8.3,I7,F12.3,F8.3)
 6170 FORMAT (///,' ANALYSIS OF DATA IN 10 EQUAL BINS OF ',
     +        '(SIN(THETA)/LAMBDA)**2 FROM',F7.4,' TO',F7.4)
 6180 FORMAT (//,
     +        '  BIN  AVERAGE NO. LAUE     MEAN  R-FACTOR  LOG(IL/IS)'/,
     +        '  NO.    STL2   REFLNS     IL/IS',/)
 6190 FORMAT (' ',I4,F8.3,I7,F12.3,F8.3,F11.3)
 6200 FORMAT (///,' ANALYSIS OF DATA IN 10 EQUAL BINS OF INTENSITY FROM'
     +        ,F12.1,'  TO',F12.1,' (UNSCALED)')
 6210 FORMAT (//,'  BIN     AVERAGE    NO. LAUE     MEAN  R-FACTOR',/,
     +        '  NO.    INTENSITY    REFLNS     IL/IS',/)
 6220 FORMAT (' ',I4,F14.1,I7,F12.3,F8.3)
 6230 FORMAT (//,' B-FACTOR DIFFERENCE = ',F6.2)
      END
**==diag1.f  
C
C SUBROUTINE 'DIAG1'
C ==================
C
C PRINT DIAGNOSTICS FROM THE SUBROUTINE 'COLECT' IF REQUIRED
C
      SUBROUTINE DIAG1(IDXL,FSQL,FSQS,BLAM,IBIN)
C
C PARAMETERS
C
C        IDXL (I)   3 WORD ARRAY HOLDING LAUE INDICES
C        FSQL (I)   INTENSITY OF LAUE REFLECTION
C        FSQS (I)   INTENSITY OF STANDARD REFLECTION
C        BLAM (I)   WAVELENGTH
C        IBIN (I)   BIN NUMBER
C
C
C SPECIFICATION STATEMENTS
C ------------------------
C
      LOGICAL FIRST
      DIMENSION IDXL(3)
      COMMON /DIAGS/NDIAG1,NDIAG2
      SAVE IFIRST,JCOUNT
      DATA FIRST/.TRUE./
C
C PRINT DIAGNOSTICS (FIRST NDIAG1)
C --------------------------------
C
      IF (NDIAG1.LE.0) RETURN
      IF (FIRST) THEN
         FIRST = .FALSE.
         WRITE (9,6010) NDIAG1
         JCOUNT = 0
      END IF
      JCOUNT = JCOUNT + 1
      IF (JCOUNT.GT.NDIAG1) RETURN
      INTL = NINT(FSQL)
      INTS = NINT(FSQS)
      WRITE (9,6020) IDXL , INTL , INTS , BLAM , IBIN
      RETURN
C
C FORMAT STATEMENTS
C -----------------
C
 6010 FORMAT ('1DIAGNOSTICS FROM SUBROUTINE COLECT',
     +        '          NUMBER REQUESTED =',I10,///,
     +        '    H    K    L    I(LAUE) I(STANDARD)  LAMBDA   BIN',/)
 6020 FORMAT (' ',I4,2I5,I11,I12,F8.3,I5)
      END
**==diag2.f  
C
C SUBROUTINE 'DIAG2'
C ==================
C
C PRINT DIAGNOSTICS FROM THE SUBROUTINE 'OPSCAL' IF REQUIRED
C
      SUBROUTINE DIAG2(IDXOLD,SC,STL,ALM,AI,SI,AVEI,KSYM,F1,F2,PAIR)
C
C PARAMETERS
C
C      IDXOLD (I)   3 WORD ARRAY HOLDING LAUE INDICES
C          SC (I)   SCALE FACTOR
C         STL (I)   'S' VALUE
C         ALM (I)   WAVELENGTH
C          AI (I)   SCALED LAUE INTENSITY
C          WI (I)   WEIGHT OF LAUE INTENSITY
C        AVEI (I)   3 WORD ARRAY HOLDING THE 3 AVERAGE (SCALED)
C                   INTENSITIES FOR THE REFLECTION
C        KSYM (I)   THE SYMMETRY NUMBER
C          F1 (I)   THE 'F' VALUE FOR I+
C          F2 (I)   THE 'F' VALUE FOR I-
C        PAIR (I)   FLAG = .TRUE. IF F+ AND F- VALUES CALCULATED
C                   FOR THIS REFLECTION, OTHERWISE = .FALSE.
C
C
C SPECIFICATION STATEMENTS
C ------------------------
C
      LOGICAL FIRST,PAIR
      DIMENSION IDXOLD(3),AVEI(3),IAVE(3)
      COMMON /DIAGS/NDIAG1,NDIAG2
      SAVE IFIRST,JCOUNT
      DATA FIRST/.TRUE./
C
C PRINT DIAGNOSTICS (FIRST NDIAG2)
C --------------------------------
C
      IF (NDIAG2.LE.0) RETURN
      IF (FIRST) THEN
         FIRST = .FALSE.
         WRITE (9,6010) NDIAG2
         JCOUNT = 0
      END IF
      JCOUNT = JCOUNT + 1
      IF (JCOUNT.GT.NDIAG2) RETURN
      DO 20 I = 1 , 3
         IAVE(I) = NINT(AVEI(I))
 20   CONTINUE
      IL = NINT(AI)
      ISIG = NINT(SI)
      IF (PAIR) THEN
         IF1 = NINT(F1)
         IF2 = NINT(F2)
         IDEL = NINT(F1-F2)
         WRITE (9,6020) IDXOLD , ALM , SC , STL , IL , ISIG , IAVE , 
     +                  KSYM , IF1 , IF2 , IDEL
      ELSE
         WRITE (9,6020) IDXOLD , ALM , SC , STL , IL , ISIG , IAVE , 
     +                  KSYM
      END IF
      RETURN
C
C FORMAT STATEMENTS
C -----------------
C
 6010 FORMAT ('1DIAGNOSTICS FROM SUBROUTINE OPSCAL',
     +        '          NUMBER REQUESTED =',I10,//,
     +        ' NOTE: AVERAGES ARE WEIGHTED',///,
     +  '    H    K    L  LAMBDA   SCALE    S        I(LAUE)      SIGMA'
     +  ,'     AVEI(1)     AVEI(2)     AVEI(3)  KSYM     F1',
     +  '      F2  DELTA',/)
 6020 FORMAT (' ',I4,2I5,F8.3,F8.3,F8.4,I12,I11,3I12,I5,2I9,I6)
      END
**==headxx.f  
C
C SUBROUTINE 'HEADXX'
C ===================
C
C THIS SUBROUTINE IS USED TO OUTPUT A HEADING
C (ORIGINAL VERSION ALLOWED FOR OVERPRINTING/UNDERLINING)
C
C SUBROUTINE CALL
C
C   CALL HEADXX(LUN,TEXT,ITERM,ICC,ITYPE)
C
C PARAMETERS
C             LUN (I) LOGICAL UNIT NUMBER FOR OUTPUT OF HEADING
C            TEXT (I) CHARACTER ARRAY CONTAINING HEADING TEXT
C           ITERM (I) CHARACTER TERMINATING HEADING TEXT STRING
C                     (NOT INCLUDED IN HEADING OUTPUT)
C             ICC (I) CARRIAGE CONTROL CHARACTER FOR HEADING
C           ITYPE (I) FLAG =1, NORMAL HEADING
C                          =2, DOUBLE SPACED HEADING
C
      SUBROUTINE HEADXX(LUN,TEXT,ITERM,ICC,ITYPE)
C
C SPECIFICATION STATEMENTS
C
      CHARACTER*1 TEXT(*),PLUS,UNDR,ITERM,ICC
      DATA PLUS/'+'/
      DATA UNDR/'_'/
C
C FIND NUMBER OF CHARACTERS IN HEADING
C
      I = 0
 20   I = I + 1
      IF (TEXT(I).NE.ITERM) GO TO 20
      NCHAR = I - 1
      IF (NCHAR.EQ.0) RETURN
C
C OUTPUT FIRST LINE OF HEADING
C
      IF (ITYPE.EQ.2) THEN
C
C DOUBLE SPACED
C
         WRITE (LUN,6020) ICC , (TEXT(N),N=1,NCHAR)
         RETURN
      ELSE
C
C NORMAL HEADING
C
         WRITE (LUN,6010) ICC , (TEXT(N),N=1,NCHAR)
         RETURN
      END IF
C
C FORMAT STATEMENTS
C
 6010 FORMAT (133A1)
 6020 FORMAT (A1,132(A1,1X))
      END
**==ffread.f  
C
C
C SUBROUTINE 'FFREAD'  ENTRIES: 'IVALUE' 'FVALUE' 'TEXLIM'
C ========================================================
C
C THIS SUBROUTINE ENABLES FREE FORMAT INPUT OF DATA. IT IS CALLED
C ONCE TO READ IN THE TEXT BUFFER  AND CALLED VIA AN ENTRY POINT
C FOR EACH INTEGER, FLOATING POINT OR TEXT FIELD IN TURN.
C
      SUBROUTINE FFREAD(IUN,LTEXT,LMAX,LLP,LSTOP)
C
C PARAMETERS
C
C       IUN (I)   UNIT NUMBER FOR READ (MAY BE SET TO ZERO IF THE TEXT
C                 TO BE INTERPRETED IS ALREADY IN THE ARRAY 'ITEXT')
C     LTEXT (I/O) CHARACTER ARRAY FOR TEXT BUFFER
C      LMAX (I)   THE MAXIMUM NO. OF CHARS IN 'ITEXT' (UP TO 150)
C       LLP (I)   UNIT NUMBER FOR ERROR MESSAGES (0 IF NONE REQUIRED)
C     LSTOP (I)   LOGICAL FLAG =TRUE, STOP IF ERROR, =FALSE CONTINUE.
C
C ENTRY IVALUE(IVAL)       RETURN INTEGER VALUE 'IVAL' FROM NEXT FIELD
C -----------------
C
C ENTRY FVALUE(FVAL)       RETURN F.P. VALUE 'FVAL' FROM NEXT FIELD
C -----------------
C
C ENTRY TEXLIM(MIN,MAX)    RETURN NOS. OF FIRST AND FINAL CHARACTERS IN
C ---------------------    FIELD IN 'MIN' AND 'MAX'
C
C NOTE: ENTRY PARAMETER VALUES ARE RETURNED AS ZERO FOR BLANK FIELDS
C       OR IF CONTINUING AFTER A SYNTAX ERROR.
C
C
C SPECIFICATION STATEMENTS
C
      CHARACTER*1 ITEXT
      CHARACTER*1 ISP,IPLUS,IMINUS,IDEC,ICOMMA
      CHARACTER*1 LTEXT(*),INUM(10)
      LOGICAL LSTOP,ISTOP,ERR
      COMMON /FFRDXX/ITEXT(150)
      COMMON /FFRDYY/NMAX,ILP,ISTOP,IFLD,M
      SAVE /FFRDXX/,/FFRDYY/
      DATA ISP/' '/,IPLUS/'+'/,IMINUS/'-'/,IDEC/'.'/,ICOMMA/','/
      DATA INUM/'0','1','2','3','4','5','6','7','8','9'/
C
C SAVE VARIABLES FOR ENTRY POINTS
C
      NMAX = LMAX
      IF (LMAX.GT.150) THEN
         IF (ILP.GT.0) WRITE (ILP,6030)
         STOP 810
      ELSE
         ILP = LLP
         ISTOP = LSTOP
C
C INITIAL ENTRY
C
         M = 0
         IFLD = 0
         IF (IUN.GT.0) READ (IUN,6010) (LTEXT(N),N=1,NMAX)
         DO 20 I = 1 , LMAX
            ITEXT(I) = LTEXT(I)
 20      CONTINUE
         RETURN
      END IF
C
C ENTRIES FOR INTEGER AND FLOATING POINT VALUES
C
      ENTRY IVALUE(IVAL)
C     ==================
C
      ITYP = 0
      IVAL = 0
      GO TO 30
      ENTRY FVALUE(FVAL)
C     ==================
C
      ITYP = 1
      FVAL = 0.0
 30   S = 1.0
      X = 0.0
      ERR = .FALSE.
      IFDEC = 0
      ID = 0
      IFLD = IFLD + 1
 40   M = M + 1
      IF (M.GT.NMAX) GO TO 100
      IF (ITEXT(M).EQ.ISP) GO TO 40
      IF (ITEXT(M).EQ.ICOMMA) GO TO 40
      IF (ITEXT(M).NE.IPLUS) THEN
         IF (ITEXT(M).NE.IMINUS) GO TO 60
         S = -1.0
      END IF
 50   M = M + 1
      IF (M.GT.NMAX) GO TO 90
      IF (ITEXT(M).EQ.ISP) GO TO 90
      IF (ITEXT(M).EQ.ICOMMA) GO TO 90
 60   IF (ITEXT(M).NE.IDEC) THEN
         DO 70 I = 1 , 10
            N = I - 1
            IF (ITEXT(M).EQ.INUM(I)) GO TO 80
 70      CONTINUE
         ERR = .TRUE.
      ELSE
         IF (IFDEC.EQ.1) ERR = .TRUE.
         IFDEC = 1
      END IF
      GO TO 50
 80   ID = ID + IFDEC
      X = 10.0*X + N
      GO TO 50
 90   X = S*X/(10.0**ID)
      IF (ERR) GO TO 130
 100  IF (ITYP.NE.0) THEN
         FVAL = X
         RETURN
      ELSE
         IF (IFDEC.NE.0) GO TO 130
         IVAL = X + 0.1*S
         RETURN
      END IF
C
C TEXT STRINGS ENTRY
C
      ENTRY TEXLIM(MIN,MAX)
C     =====================
C
      MIN = 0
      MAX = 0
      IFLD = IFLD + 1
 110  M = M + 1
      IF (M.GT.NMAX) RETURN
      IF (ITEXT(M).EQ.ISP) GO TO 110
      MIN = M
 120  MAX = M
      M = M + 1
      IF (M.GT.NMAX) RETURN
      IF (ITEXT(M).NE.ISP) GO TO 120
      RETURN
C
C ERROR FOUND
C
 130  IF (ILP.GT.0) WRITE (ILP,6020) IFLD , (ITEXT(I),I=1,NMAX)
      IF (ISTOP) STOP
      RETURN
C
C FORMAT STATEMENTS
C
 6010 FORMAT (255A1)
 6020 FORMAT (' **SYNTAX ERROR IN FIELD',I3,' **'/(1X,130A1/))
 6030 FORMAT (' **MAXIMUM LENGTH OF TEXT ALLOWED IN FFREAD IS 150',
     +        ' CHARACTERS**')
      END
**==polyft.f  
C
C SUBROUTINE 'POLYFT'
C ===================
C
C POLYNOMIAL LEAST SQUARES CURVE FITTING.
C
      SUBROUTINE POLYFT (M,X,Y,W,KMX,N,SI,P)
C
C PARAMETERS
C
C           M (I)   NUMBER OF POINTS TO BE FITTED
C           X (I)   ARRAY HOLDING ABSCISSA POINTS (DOUBLE PRECISION)
C           Y (I)   ARRAY HOLDING FUNCTION VALUES  (DOUBLE PRECISION)
C           W (I)   ARRAY HOLDING INPUT WEIGHTS FOR THE POINTS (DOUBLE
C                   PRECISION)
C         KMX (I)   MAXIMUM ORDER OF THE POLYNOMIAL FOR WHICH
C                   THE CURVE FITTING IS TO BE ASSESSED.
C           N (I)   ORDER OF POLYNOMIAL TO BE FITTED.
C          SI (O)   ARRAY RETURNING RESIDUALS FOR POLYNOMIAL ORDERS
C                   0 TO KMX  (KMX+1) VALUES   (DOUBLE PRECISION)
C           P (O)   ARRAY RETURNING POLYNOMIAL COEFFICIENTS P0..PN
C                   (N+1 VALUES)   (DOUBLE PRECISION)
C
C
C PARAMETER STATEMENT  MXP = MAXIMUM NUMBER OF POINTS ALLOWED
C                      MXO = MAXIMUM ORDER OF POLYNOMIAL ALLOWED
C
      PARAMETER(MXP=200,MXO=20)
C
C SPECIFICATION STATEMENTS
C ------------------------
C
      DOUBLE PRECISION X(*),Y(*),W(*),SI(*),P(*)
      DOUBLE PRECISION CTP(MXO+1),CPSAVE(MXO+1),CP(MXO+1),CLP(MXO+2)
      DOUBLE PRECISION AL(MXO+1),BE(MXO+1),S(MXO+1)
      DOUBLE PRECISION XLP(MXP),TP(MXP),XX(MXP)
      DOUBLE PRECISION A,B,DU,XTEMP,TW,OM,DELSQ,SIMIN,XLW,TPSAV
C
C INITIALISATIONS
C ---------------
C
      IF (N.GT.MXO .OR. KMX.GT.MXO)
     +    STOP '**MAXIMUM ORDER EXCEEDED FOR POLYFT**'
      IF (M.GT.MXP) STOP '**TOO MANY POINTS FOR POLYFT**'
      N1 = N + 1
      K1 = KMX + 1
      K = K1 - 1
      DO 20 I = 1 , K1
         CP(I) = 0.0
 20   CONTINUE
      BE(1) = 0.0
      CLP(2) = 0.0
      CLP(1) = 0.0
      DELSQ = 0.0
      OM = 0.0
      TW = 0.0
      CTP(1) = 1.0
C
C PERFORM THE CURVE FITTING
C -------------------------
C
C
C DELSQ=SUM(WF(X)**2)
C OM=SUM(WF(X))
C TW=SUM(W)
C
      DO 30 I = 1 , M
         XX(I) = X(I)
         DELSQ = DELSQ + W(I)*Y(I)**2
         TP(I) = 1.0
         XLP(I) = 0.0
         OM = OM + W(I)*Y(I)
         TW = TW + W(I)
 30   CONTINUE
      S(1) = OM/TW
      CP(1) = S(1)
C
C DELSQ=SUM(WF**2)-(SUM(WF))**2
C
      DELSQ = DELSQ - S(1)*OM
      XTEMP = M - 1
C
C SD=DELSQ/(M-1)
C
      SI(1) = ABS(DELSQ/XTEMP)
      A = 4.0/(XX(M)-XX(1))
      B = -2.0 - A*XX(1)
      DO 40 I = 1 , M
         XX(I) = A*XX(I) + B
 40   CONTINUE
      DO 100 I = 1 , K
         DU = 0.0
         DO 50 J = 1 , M
            DU = DU + W(J)*XX(J)*TP(J)**2
 50      CONTINUE
         AL(I+1) = DU/TW
         XLW = TW
         TW = 0.0
         OM = 0.0
         DO 60 J = 1 , M
            DU = BE(I)*XLP(J)
            TPSAV = TP(J)
            TP(J) = (XX(J)-AL(I+1))*TPSAV - DU
            TW = TW + W(J)*TP(J)**2
            OM = OM + W(J)*Y(J)*TP(J)
            XLP(J) = TPSAV
 60      CONTINUE
         BE(I+1) = TW/XLW
         S(I+1) = OM/TW
         DELSQ = DELSQ - S(I+1)*OM
         XTEMP = M - I - 1
         IF (M-I.GT.1) THEN
            SI(I+1) = ABS(DELSQ/XTEMP)
         ELSE
            SI(I+1) = 0.0
         END IF
         IF (N1.EQ.I) THEN
            DO 70 J = 1 , K1
               CPSAVE(J) = CP(J)
 70         CONTINUE
         END IF
         DO 80 J = 1 , I
            DU = CLP(J+1)*BE(I)
            CLP(J+1) = CTP(J)
            CTP(J) = CLP(J) - AL(I+1)*CTP(J) - DU
            CP(J) = CP(J) + S(I+1)*CTP(J)
 80      CONTINUE
         CP(I+1) = S(I+1)
         CTP(I+1) = 1.0
         CLP(I+2) = 0.0
         IF (I.EQ.K) THEN
            IF (N1.LT.K1) THEN
               DO 90 J = 1 , K1
                  CP(J) = CPSAVE(J)
 90            CONTINUE
            END IF
         END IF
 100  CONTINUE
C
C PREPARE POLYNOMIAL COEFFICIENTS
C -------------------------------
C
      CLP(1) = 1
      CPSAVE(1) = 1
      P(1) = CP(1)
      DO 110 I = 2 , K1
         CLP(I) = 1
         CPSAVE(I) = B*CPSAVE(I-1)
         P(1) = P(1) + CP(I)*CPSAVE(I)
 110  CONTINUE
      DO 130 J = 2 , K1
         CLP(1) = CLP(1)*A
         P(J) = CP(J)*CLP(1)
         KK = 2
         J1 = J + 1
         IF (J1.GT.K1) GO TO 140
         DO 120 I = J1 , K1
            CLP(KK) = A*CLP(KK) + CLP(KK-1)
            P(J) = P(J) + CP(I)*CLP(KK)*CPSAVE(KK)
            KK = KK + 1
 120     CONTINUE
 130  CONTINUE
 140  A = 1.0/A
      RETURN
      END
**==colabs.f  
C
C SUBROUTINE 'COLABS'
C ===================
C
C THIS SUBROUTINE COLLECTS THE DATA TO BE USED FOR THE ABSORPTION
C CORRECTION AND CALCULATES AN ABSORPTION SURFACE USING A
C BIVARIATE POLYNOMIAL DISTRIBUTION
C
      SUBROUTINE COLABS(IA,IMAX,NREF)
C
C PARAMETERS
C
C         IA (I)   ARRAY TO HOLD LAUE DATA (5 WORDS/REFLN)
C       IMAX (I)   NO. OF WORDS IN 'IA'
C       NREF (I)   NUMBER OF LAUE REFLECTIONS
C
      PARAMETER(MXR=25000,MXL=8)
C
C SPECIFICATION STATEMENTS
C ------------------------
C
      LOGICAL EOF
      DIMENSION IA(IMAX),IDXL(3),ADATA(10),IDATA(3)
      DOUBLE PRECISION AN(121,121),AI(121,121),A(121),QN(121),RN(121)
      DOUBLE PRECISION RHS(121),G,Q,R,DET,RMS,SCMAT
      DOUBLE PRECISION GCAL(16),QQ,RR
      DIMENSION GAVE(16,16),NAVE(16,16)
      LOGICAL DIAGHD
      LOGICAL UPDATB,NRMCAL,ABSCAL
      COMMON /SCFLAG/UPDATB,NRMCAL,ABSCAL
      COMMON /CDATA/AMINL,AMAXL,RES1,RES2,TH1,TH2,BFAC,NRANGE,IBATCH,
     *              NITER,IABS,LASTIT,KBATCH
      DOUBLE PRECISION PS
      COMMON /ABSURF/PS(121),NQORD,NRORD,QXMIN,QXMAX,RYMIN,RYMAX
      COMMON /SCAL/SCALEI,SCAL2,SCALIN,SIGREJ
      COMMON /NWR/NWR
      COMMON /ABSCUT/ISMIN,ILMIN,DIAGR
      LOGICAL LAMOPT
      COMMON /LAMO/ LAMOPT, DLAM
      LOGICAL CENTYP
      COMMON /CENT/ CENTYP
      EQUIVALENCE (IAA,AAA)
      DOUBLE PRECISION EVALSP
      DATA NTWO16/65536/
C
C INITIALISATIONS
C ---------------
C
      NRABS = 0
      NREJS = 0
      NREJL = 0
      DIAGHD = .FALSE.
      NR = 1
      NPAR = (NQORD+1)*(NRORD+1)
      CALL LRREWD(1)
      DO 30 I = 1 , NPAR
         RHS(I) = 0
         DO 20 J = 1 , NPAR
            AN(J,I) = 0
 20      CONTINUE
 30   CONTINUE
      DO 50 I = 1 , 16
         DO 40 J = 1 , 16
            GAVE(I,J) = 0.0
            NAVE(I,J) = 0
 40      CONTINUE
 50   CONTINUE
C
C GET FIRST LAUE AND FIRST MTZ REFLECTION
C ---------------------------------------
C
      CALL GETHKL(IA(1),IDXL,ICS)
      BLMBDA = IA(2)/10000.0
      CALL LRREFF(1,RESOL,ADATA,EOF)
      IF (EOF) GO TO 140
      IDATA(1) = NINT(ADATA(1))
      IDATA(2) = NINT(ADATA(2))
      IDATA(3) = NINT(ADATA(3))
      ALMBDA = 0.0
      KCS = 0
      IF (LAMOPT) THEN
          ALMBDA = ADATA(6)
          IF (CENTYP) KCS = NINT(ADATA(7))
      ELSE
          IF (CENTYP) KCS = NINT(ADATA(6))
      ENDIF
C
C COMPARE LAUE AND MTZ REFLECTIONS AND READ FURTHER REFLECTIONS FROM
C ------------------------------------------------------------------
C
C WHICHEVER SET IS REQUIRED TILL A MATCH IS FOUND
C -----------------------------------------------
C
 60   IF (ICMPR(IDXL,IDATA,ALMBDA,BLMBDA,KCS,ICS,LAMOPT,DLAM)
     +    .LT.0) THEN
      ELSE IF (ICMPR(IDXL,IDATA,ALMBDA,BLMBDA,KCS,ICS,LAMOPT,DLAM)
     +         .EQ.0) THEN
C
C MATCHING REFLECTION FOUND, COLLECT DATA FOR ANALYSIS AND CURVE FIT
C ------------------------------------------------------------------
C
         IB = NWR*(NR-1)
         SIGL = MOD(IA(IB+4),NTWO16)
         BLAM = IA(IB+2)/10000.0
         THETA = IA(IB+5)/10000.0
         IAA = IA(IB+3)
         FFUNSC = AAA
         FSQL = AAA*SCFAC(IA(IB+1),NRMCAL,.FALSE.)
         FSQS = ADATA(4)*ADATA(4)
         IF (ISMIN.GT.0 .AND. FSQS.LT.ISMIN) THEN
            NREJS = NREJS + 1
            GO TO 130
         END IF
         IF (ILMIN.GT.0 .AND. FFUNSC.LT.ILMIN) THEN
            NREJL = NREJL + 1
            GO TO 130
         END IF
         FSQS = FSQS/SCALEI
         IF (SIGL.NE.0.0 .AND. FSQL.NE.0.0) THEN
            IF (FSQS.NE.0.0 .AND. ADATA(5).NE.0.0) THEN
               NRABS = NRABS + 1
               WSQT = 1.0
               CALL QRNRM(IA(IB+1),IABS,Q,R)
C****      Q=(2*Q-QXMIN-QXMAX)/(QXMAX-QXMIN)
C****      R=(2*R-RYMIN-RYMAX)/(RYMAX-RYMIN)
               G = FSQS/FSQL
               IQ = IRANGE(REAL(Q),-1.0,1.0,16)
               IR = IRANGE(REAL(R),-1.0,1.0,16)
               GAVE(IR,IQ) = GAVE(IR,IQ) + G
               NAVE(IR,IQ) = NAVE(IR,IQ) + 1
               IF (G.GE.DIAGR .OR. G.LE.1.0/DIAGR) THEN
                  IF (.NOT.DIAGHD) THEN
                     WRITE (6,6140)
                     DIAGHD = .TRUE.
                  END IF
                  IAA = IA(IB+6)
                  XX = AAA
                  IAA = IA(IB+7)
                  YY = AAA
                  WRITE (6,6150) IDXL , FFUNSC , FSQS , FSQL , THETA , 
     +                           BLAM , G , XX , YY , Q , R , IQ , IR
               END IF
               QN(1) = 1.0
               IF (NQORD.GT.0) THEN
                  DO 70 I = 2 , NQORD + 1
                     QN(I) = QN(I-1)*Q
 70               CONTINUE
               END IF
               RN(1) = 1.0
               IF (NRORD.GT.0) THEN
                  DO 80 I = 2 , NRORD + 1
                     RN(I) = RN(I-1)*R
 80               CONTINUE
               END IF
               NP = 0
               DO 100 IQ = 1 , NQORD + 1
                  DO 90 IR = 1 , NRORD + 1
                     NP = NP + 1
                     A(NP) = WSQT*RN(IR)*QN(IQ)
                     RHS(NP) = RHS(NP) + WSQT*G*A(NP)
 90               CONTINUE
 100           CONTINUE
               DO 120 I = 1 , NP
                  DO 110 J = 1 , NP
                     AN(J,I) = AN(J,I) + A(I)*A(J)
 110              CONTINUE
 120           CONTINUE
            END IF
         END IF
      ELSE
C
C NO MATCH, LAUE>STANDARD, GET NEXT STANDARD
C
         CALL LRREFF(1,RESOL,ADATA,EOF)
         IF (EOF) GO TO 140
         IDATA(1) = NINT(ADATA(1))
         IDATA(2) = NINT(ADATA(2))
         IDATA(3) = NINT(ADATA(3))
         ALMBDA = 0.0
         KCS = 0
         IF (LAMOPT) THEN
             ALMBDA = ADATA(6)
             IF (CENTYP) KCS = NINT(ADATA(7))
         ELSE
             IF (CENTYP) KCS = NINT(ADATA(6))
         ENDIF
         GO TO 60
      END IF
C
C NO MATCH, LAUE<STANDARD, GET NEXT LAUE
C
 130  NR = NR + 1
      IF (NR.LE.NREF) THEN
         I = NWR*NR - NWR + 1
         CALL GETHKL(IA(I),IDXL,ICS)
         BLMBDA = IA(I+1)/10000.0
         GO TO 60
      END IF
C
C END OF DATA REACHED, INVERT NORMAL MATRIX AND CALC SURFACE POLYNOMIAL
C ---------------------------------------------------------------------
C
 140  WRITE (6,6010) NRABS
      WRITE (6,6020) NREJS
      WRITE (6,6030) NREJL
      IF (NRABS.LT.NPAR) THEN
         WRITE (6,6120)
         STOP
      ELSE
         SCMAT = 0.0
         DO 150 I = 1 , NP
            SCMAT = SCMAT + AN(I,I)
 150     CONTINUE
         SCMAT = SCMAT/NP
         DO 170 I = 1 , NP
            DO 160 J = 1 , NP
               AN(J,I) = AN(J,I)/SCMAT
 160        CONTINUE
 170     CONTINUE
         CALL INVMAT(AN,AI,NPAR,121,6,DET,IERR)
         IF (IERR.GT.0) THEN
            WRITE (6,6130)
            STOP
         ELSE
C
C EVALUATE POLYNOMIAL COEFFICIENTS
C
            DO 190 I = 1 , NPAR
               PS(I) = 0.0
               DO 180 J = 1 , NPAR
                  PS(I) = PS(I) + AI(I,J)*RHS(J)
 180           CONTINUE
               PS(I) = PS(I)/SCMAT
 190        CONTINUE
C
C PRINT TABLE OF AVERAGED BIN RATIOS IS/IL (BEFORE CORRECTION)
C ON A 16*16 GRID
C
            DO 210 I = 1 , 16
               DO 200 J = 1 , 16
                  IF (NAVE(I,J).GT.0) GAVE(I,J) = GAVE(I,J)/NAVE(I,J)
 200           CONTINUE
 210        CONTINUE
            WRITE (6,6060)
            DO 220 IR = 16 , 1 , -1
               WRITE (6,6070) (GAVE(IR,IQ),IQ=1,16)
               WRITE (6,6080) (NAVE(IR,IQ),IQ=1,16)
 220        CONTINUE
C
C PRINT POLYNOMIAL COEFFICIENTS
C
            WRITE (6,6040) NPAR
            WRITE (6,6050) (PS(I),I=1,NPAR)
C
C PLOT THE SURFACE ON A 16*16 GRID
C
            WRITE (6,6090)
            DO 240 IRR = 15 , -15 , -2
               RR = IRR/16.0
               DO 230 IQQ = -15 , 15 , 2
                  QQ = IQQ/16.0
                  III = (IQQ+17)/2
                  GCAL(III) = EVALSP(QQ,RR,PS,NQORD,NRORD,6,IERR)
 230           CONTINUE
               WRITE (6,6100) GCAL
 240        CONTINUE
            RETURN
         END IF
      END IF
C
C FORMAT STATEMENTS
C -----------------
C
 6010 FORMAT (//,' NUMBER OF OBSERVATIONS FOR ABSORPTION SURFACE',
     +        ' CALCULATION =',I6)
 6020 FORMAT (/,' NUMBER REJECTED BY ISMIN CUTOFF =',I6)
 6030 FORMAT (/,' NUMBER REJECTED BY ILMIN CUTOFF =',I6)
 6040 FORMAT (//,' POLYNOMIAL COEFFICIENTS FOR ABSORPTION SURFACE',
     +        ' PS 1  TO  PS',I3,/)
 6050 FORMAT (' ',10E12.4)
 6060 FORMAT (///,' TABLE OF AVERAGE IS/IL RATIOS ON A 16*16 GRID',//)
 6070 FORMAT (/,' ',16F7.2)
 6080 FORMAT (' ',16I7)
 6090 FORMAT ('1A B S O R P T I O N   S U R F A C E')
 6100 FORMAT (/,' ',17F7.2)
C
 6120 FORMAT (/,' **FEWER OBSERVATIONS THAN PARAMETERS FOR SURFACE',
     +        ' FITTING**')
 6130 FORMAT (/,' **UNABLE TO FIT ABSORPTION SURFACE**')
C
 6140 FORMAT (//,'    H   K   L  UNSCALED STANDARD    LAUE     THETA',
     +        ' LAMBDA  IS/IL FILM-X FILM-Y      Q      R IQ IR',/,
     +        '                LAUE I  INTENSITY INTENSITY',/)
 6150 FORMAT (' ',3I4,3F10.1,7F7.2,2I3)
      END
**==invmat.f  
C
C SUBROUTINE 'INVMAT'
C ===================
C
C INVERT A SQUARE MATRIX
C
C (INVERSION BY ELIMINATION WITH PARTIAL PIVOTING)
C
      SUBROUTINE INVMAT(A,AI,N,IDIM,MSG,DET,IERR)
C
C PARAMETERS
C
C       A (I/O) MATRIX TO BE INVERTED A(I,J), I= 1 TO N, AND J = 1 TO N
C               (DOUBLE PRECISION)   NOTE: THE ARRAY IS MODIFIED DURING
C               THE COMPUTATION.
C      AI (O)   INVERTED MATRIX (DOUBLE PRECISION)
C       N (I)   SIZE OF SQUARE MATRIX
C    IDIM (I)   FIRST DIMENSION OF ARRAYS A AND AI (MUST BE THE SAME FOR
C               EACH ARRAY)
C     MSG (I)   LOGICAL UNIT NUMBER FOR ERROR MESSAGES ,0 IF NONE DESIRED
C     DET (O)   RETURNS DETERMINANT VALUE (DOUBLE PRECISION)
C    IERR (O)   ERROR FLAG =0 OK
C                          =1 SINGULAR MATRIX
C
C
C SPECIFICATION STATEMENTS
C ------------------------
C
      DOUBLE PRECISION A(IDIM,N),AI(IDIM,N),DET
      DOUBLE PRECISION EPS,AMAX,ATMP,AITMP,DIV,AMULT
C
C PERFORM MATRIX INVERSION
C ------------------------
C
C INITIALISATIONS, CONSTRUCT IDENTITY MATRIX
C
      IERR = 0
      EPS = 1.0E-07
      DO 30 I = 1 , N
         DO 20 J = 1 , N
            AI(I,J) = 0.0
 20      CONTINUE
         AI(I,I) = 1.0
 30   CONTINUE
      DET = 1.0
      DO 90 K = 1 , N
C
C LOCATE MAX MAGNITUDE ON OR BELOW MAIN DIAGONAL
C
         IF (K.LT.N) THEN
            IMAX = K
            AMAX = ABS(A(K,K))
            DO 40 I = K + 1 , N
               IF (AMAX.LT.ABS(A(I,K))) THEN
                  IMAX = I
                  AMAX = ABS(A(I,K))
               END IF
 40         CONTINUE
C
C INTERCHANGE ROWS IMAX AND K IF IMAX.NE.K
C
            IF (IMAX.NE.K) THEN
               DO 50 J = 1 , N
                  ATMP = A(IMAX,J)
                  A(IMAX,J) = A(K,J)
                  A(K,J) = ATMP
                  AITMP = AI(IMAX,J)
                  AI(IMAX,J) = AI(K,J)
                  AI(K,J) = AITMP
 50            CONTINUE
               DET = -DET
            END IF
         END IF
C
C TEST FOR SINGULAR MATRIX
C
         IF (ABS(A(K,K)).LT.EPS) GO TO 100
         DET = A(K,K)*DET
C
C DIVIDE PIVOT ROW BY ITS MAIN DIAGONAL ELEMENT
C
         DIV = A(K,K)
         DO 60 J = 1 , N
            A(K,J) = A(K,J)/DIV
            AI(K,J) = AI(K,J)/DIV
 60      CONTINUE
C
C REPLACE EACH ROW BY LINEAR COMBINATION WITH PIVOTAL ROW
C
         DO 80 I = 1 , N
            AMULT = A(I,K)
            IF (I.NE.K) THEN
               DO 70 J = 1 , N
                  A(I,J) = A(I,J) - AMULT*A(K,J)
                  AI(I,J) = AI(I,J) - AMULT*AI(K,J)
 70            CONTINUE
            END IF
 80      CONTINUE
C
C END OF MAIN LOOP
C
 90   CONTINUE
      RETURN
C
C ERROR CONDITIONS
C ----------------
C
 100  IF (MSG.GT.0) WRITE (MSG,6010) K
      IERR = 1
      RETURN
C
C FORMAT STATEMENTS
C -----------------
C
 6010 FORMAT (/,'**SINGULAR MATRIX FOR ELEMENT',I5,' **')
      END
**==scfac.f  
C
C FUNCTION 'SCFAC'
C ================
C
C RETURN THE REFLECTION SCALING FACTOR APPLYING TEMPERATURE FACTOR, LAMBDA
C AND ABSORPTION FACTOR SCALING WHERE RELEVANT.
C
      FUNCTION SCFAC(IAREF,NRMC,ABSC)
C
C PARAMETERS
C
C       IAREF (I)   ARRAY CONTAINING THE WORDS FOR THE REFLECTION STARTING
C                   WITH THE WORD CONTAINING THE INDICES
C        NRMC (I)   FLAG = .TRUE.  INCLUDE LAMBDA SCALING FACTOR
C                        = .FALSE. DO NOT
C        ABSC (I)   FLAG = .TRUE.  INCLUDE ABSORPTION SCALING FACTOR
C                        = .FALSE. DO NOT
C
C SPECIFICATION STATEMENTS
C ------------------------
C
      DIMENSION IAREF(*)
      LOGICAL NRMC,ABSC
      COMMON /CDATA/AMINL,AMAXL,RES1,RES2,TH1,TH2,BFAC,NRANGE,IBATCH,
     *              NITER,IABS,LASTIT,KBATCH
      DOUBLE PRECISION PS
      COMMON /ABSURF/PS(121),NQORD,NRORD,QXMIN,QXMAX,RYMIN,RYMAX
      DOUBLE PRECISION Q,R
      DOUBLE PRECISION EVALSP
C
C GET SCALE FACTOR
C ----------------
C
      SCFAC = 1.0
      BLAM = IAREF(2)/10000.0
      THETA = IAREF(5)/10000.0
      IF (BFAC.NE.0.0) THEN
         STOL = SIN(THETA)/BLAM
         STOL2 = STOL*STOL
         SCFAC = EXP(-BFAC*STOL2)
      END IF
C
C GET LAMBDA SCALING FACTOR
C -------------------------
C
      IF (NRMC) SCFAC = SCFAC*1.0/CURVAL(BLAM)
C
C GET ABSORPTION SCALING FACTOR
C -----------------------------
C
      IF (.NOT.ABSC) RETURN
      CALL QRNRM(IAREF,IABS,Q,R)
      SCABS = EVALSP(Q,R,PS,NQORD,NRORD,6,IERR)
      SCFAC = SCFAC*SCABS
      RETURN
      END
**==qrnrm.f  
C
C SUBROUTINE 'QRNRM'
C ==================
C
C CALCULATE THE NORMALISED 'Q' AND 'R' ABSORPTION SURFACE COORDINATES
C FOR A REFLECTION
C
      SUBROUTINE QRNRM(IAREF,IABS,Q,R)
C
C PARAMETERS
C
C      IAREF (I)   ARRAY CONTAINING THE WORDS FOR THE REFLECTION
C       IABS (I)   ABSORPTION TYPE FLAG
C          Q (O)   NORMALISED Q COORDINATE (DOUBLE PRECISION)
C          R (O)   NORMALISED R COORDINATE (DOUBLE PRECISION)
C
C SPECIFICATION STATEMENTS
C ------------------------
C
      DIMENSION IAREF(*)
      DOUBLE PRECISION PS
      COMMON /ABSURF/PS(121),NQORD,NRORD,QXMIN,QXMAX,RYMIN,RYMAX
      DOUBLE PRECISION Q,R
C
C CALCULATE COORDINATES
C ---------------------
C
      CALL QRVAL(IAREF,IABS,Q,R)
      Q = (2*Q-QXMIN-QXMAX)/(QXMAX-QXMIN)
      R = (2*R-RYMIN-RYMAX)/(RYMAX-RYMIN)
      RETURN
      END
**==qrval.f  
C
C SUBROUTINE 'QRVAL'
C ==================
C
C CALCULATE THE 'Q' AND 'R' ABSORPTION SURFACE COORDINATES
C FOR A REFLECTION
C
      SUBROUTINE QRVAL(IAREF,IABS,Q,R)
C
C PARAMETERS
C
C      IAREF (I)   ARRAY CONTAINING THE WORDS FOR THE REFLECTION
C       IABS (I)   ABSORPTION TYPE FLAG
C          Q (O)   NORMALISED Q COORDINATE (DOUBLE PRECISION)
C          R (O)   NORMALISED R COORDINATE (DOUBLE PRECISION)
C
C SPECIFICATION STATEMENTS
C ------------------------
C
      DIMENSION IAREF(*)
      DOUBLE PRECISION Q,R
      EQUIVALENCE (IAA,AAA)
C
C CALCULATE COORDINATES
C ---------------------
C
      Q = 0.0
      R = 0.0
      IF (IABS.LE.0) RETURN
      IF (IABS.EQ.1) THEN
         THETA = IAREF(5)/10000.0
         IAA = IAREF(6)
         X = AAA
         IAA = IAREF(7)
         Y = AAA
         D = SQRT(X*X+Y*Y)/TAN(2.0*THETA)
         Q = ATAN(X/D)
         R = ATAN(Y/D)
      END IF
      IF (IABS.EQ.2) THEN
         IAA = IAREF(6)
         Q = AAA
         IAA = IAREF(7)
         R = AAA
      END IF
      RETURN
      END
**==evalsp.f  
C
C FUNCTION 'EVALSP'
C =================
C
C EVALUATE SURFACE POINT VALUE FROM THE POLYNOMIAL COEFFICIENTS
C GENERATED BY THE SUBROUTINE 'POLYSP' AND DEFINED THEREIN
C
      DOUBLE PRECISION FUNCTION EVALSP(X,Y,P,NXORD,NYORD,MSG,IERR)
C
C PARAMETERS
C
C           X (I)   X VALUE FOR THE POINT (DOUBLE PRECISION)
C           Y (I)   Y VALUE FOR THE POINT (DOUBLE PRECISION)
C           P (I)   ARRAY HOLDING THE POLYNOMIAL COEFFICIENTS
C                   (DOUBLE PRECISION)
C       NXORD (I)   THE POLYNOMIAL ORDER  TO BE FITTED FOR X
C       NYORD (I)   THE POLYNOMIAL ORDER  TO BE FITTED FOR Y
C         MSG (I)   UNIT NUMBER FOR ERROR MESSAGES (=0 IF NONE REQUIRED)
C        IERR (O)   ERROR FLAG = 0   OK
C                              = 1   TOO HIGH AN ORDER OF POLYNOMIAL FOR
C                                    THE ROUTINE
C
C PARAMETER STATEMENTS
C
C        MAXORD = MAXIMUM ORDER OF POLYNOMIAL IN X OR Y
C
      PARAMETER(MAXORD=10)
C
C SPECIFICATION STATEMENTS
C ------------------------
C
      DOUBLE PRECISION X,Y,P(*)
      DOUBLE PRECISION XN(MAXORD+1),YN(MAXORD+1)
      DOUBLE PRECISION GCAL
C
C EVALUATE VALUE FOR THE POINT X,Y
C --------------------------------
C
      IF (NXORD.GT.MAXORD .OR. NYORD.GT.MAXORD) THEN
C
C ERROR CONDITIONS
C ----------------
C
         IF (MSG.GT.0) WRITE (MSG,6010) MAXORD
         IERR = 1
         RETURN
      ELSE
         GCAL = 0.0
         IP = 0
         XN(1) = 1.0
         YN(1) = 1.0
         IF (NXORD.GT.0) THEN
            DO 20 I = 2 , NXORD + 1
               XN(I) = XN(I-1)*X
 20         CONTINUE
         END IF
         IF (NYORD.GT.0) THEN
            DO 30 I = 2 , NYORD + 1
               YN(I) = YN(I-1)*Y
 30         CONTINUE
         END IF
         DO 50 I = 1 , NXORD + 1
            DO 40 J = 1 , NYORD + 1
               IP = IP + 1
               GCAL = GCAL + P(IP)*YN(J)*XN(I)
 40         CONTINUE
 50      CONTINUE
         EVALSP = GCAL
         RETURN
      END IF
C
C FORMAT STATEMENTS
C -----------------
C
 6010 FORMAT (/,' **ROUTINE ONLY SET UP FOR A MAXIMUM ORDER OF',I3,'**')
      END
**==gtvals.f  
C
C SUBROUTINE 'GTVALS'
C ===================
C
C THIS SUBROUTINE READS A FIXED NUMBER OF REAL VALUES IN FREE
C FORMAT. THE VALUES MAY BE INPUT OVER SEVERAL CARDS
C
      SUBROUTINE GTVALS(IUN,ITEXT,IMAX,NVALS,A,ILP,ISTOP)
C
C PARAMETERS
C
C         IUN (I)   INPUT UNIT NO. FOR READ
C       ITEXT (I/O) CHARACTER ARRAY FOR READING INPUT
C        IMAX (I)   NO. OF CHARACTERS TO BE READ IN TO 'ITEXT'
C       NVALS (I)   THE NUMBER OF REAL VALUES TO BE INPUT
C           A (O)   ARRAY TO CONTAIN THE REAL VALUES
C         ILP (I)   UNIT NO. FOR ERROR MESSAGES (MAY BE 0 IF NONE
C                   ARE REQUIRED)
C       ISTOP (I)   FLAG =.TRUE. STOP AFTER SYNTAX ERROR
C                        =.FALSE. CONTINUE AFTER ERROR
C
C SUBROUTINES CALLED: FFREAD NFIELD
C
C
C SPECIFICATION STATEMENTS
C
      LOGICAL ISTOP
      CHARACTER*1 ITEXT(*)
      DIMENSION A(*)
C
C CODE
C
      K = 0
 20   CALL FFREAD(IUN,ITEXT,IMAX,ILP,ISTOP)
      N = NFIELD(ITEXT,IMAX)
      IF (N.EQ.0) GO TO 20
      DO 30 I = 1 , N
         K = K + 1
         IF (K.LE.NVALS) CALL FVALUE(A(K))
 30   CONTINUE
      IF (K.LT.NVALS) GO TO 20
      RETURN
      END
**==decnum.f  
C
C
C
C SUBROUTINE 'DECNUM'
C ===================
C
C THIS SUBROUTINE INTERPRETS A NUMERICAL FIELD WITHIN A TEXT STRING
C
C
C
C SUBROUTINE CALL    CALL DECNUM(ITEXT,MIN,MAX,ITYPE,X,IERR)
C
C PARAMETERS
C            ITEXT (I) CHARACTER ARRAY CONTAINING THE INPUT LINE
C                      OF TEXT
C              MIN (I) THE FIRST CHARACTER OF THE FIELD TO BE DECODED
C              MAX (I) THE FINAL CHARACTER OF THE FIELD TO BE DECODED
C            ITYPE (I) FLAG =1 INTEGER NUMBER TO BE INTERPRETED
C                           =2 REAL NUMBER TO BE INTERPRETED
C                X (O) THE VALUE FOUND FOR THE NUMBER (=0.0 IF ABSENT)
C             IERR (O) FLAG =0 VALID SYNTAX
C                           =1 ERROR IN SYNTAX
C
C
      SUBROUTINE DECNUM(ITEXT,MIN,MAX,ITYPE,X,IERR)
C
C DECLARATIONS
C
      CHARACTER*1 INUM(10),ITEXT(MAX),IPOINT,IPLUS,IMINUS
C
C DATA STATEMENTS
C
      DATA INUM/'1','2','3','4','5','6','7','8','9','0'/
      DATA IPOINT/'.'/
      DATA IPLUS/'+'/
      DATA IMINUS/'-'/
C
C INITIALISATION STATEMENTS
C
      X = 0.0
      A = 0.0
      IP = 0
      NP = 0
      S = 1.0
C
C CHECK VALIDITY OF MIN AND MAX
C
      IF (MIN.LT.1) GO TO 60
      IF (MAX.LT.MIN) GO TO 60
C
C INTERPRET NUMBER
C
      DO 40 I = MIN , MAX
         IF (I.LE.MIN) THEN
C FIRST CHARACTER MAY BE A SIGN
            IF (ITEXT(I).EQ.IPLUS .OR. ITEXT(I).EQ.IMINUS) THEN
               IF (ITEXT(I).EQ.IMINUS) S = -1.0
               GO TO 40
            END IF
         END IF
C CHECK FOR DECIMAL POINT AND ALSO THAT ONLY ONE IS PRESENT
         IF (ITEXT(I).NE.IPOINT) THEN
C CHECK THAT OTHER CHARACTERS ARE DIGITS
            DO 20 K = 1 , 10
               IF (ITEXT(I).EQ.INUM(K)) GO TO 30
 20         CONTINUE
            GO TO 60
         ELSE
            IP = IP + 1
            IF (IP.LE.1) GO TO 40
            GO TO 60
         END IF
C VALID DIGIT FOUND  ADD IT INTO NUMBER
 30      IF (K.EQ.10) K = 0
         IF (IP.EQ.1) NP = NP + 1
         A = 10.0*A + K
 40   CONTINUE
C EXTRA CHECK FOR INTEGER IF REQUIRED
      IF (ITYPE.EQ.1) THEN
         IF (IP.EQ.1) GO TO 60
      END IF
C FORM RESULT
      IERR = 0
      AA = S*A
      IF (NP.NE.0) THEN
         DO 50 I = 1 , NP
            AA = AA/10.0
 50      CONTINUE
      END IF
      X = AA
      RETURN
C
C ERROR RETURN
C
 60   IERR = 1
      RETURN
      END
**==nfield.f  
C
C FUNCTION 'NFIELD'
C =================
C
C&&& ROTAVATA VERSION ALLOWS COMMA SEPARATORS
C
C RETURNS THE NUMBER OF CHARACTER FIELDS PRESENT IN A TEXT STRING
C (CHARACTER STRINGS ARE ASSUMED TO BE SEPARATED BY ONE OR MORE SPACES
C  OR BY COMMAS)
C
      FUNCTION NFIELD(ITEXT,IMAX)
C
C PARAMETERS
C
C       ITEXT (I) CHARACTER ARRAY CONTAINING THE TEXT STRING
C        IMAX (I) THE NO. OF CHARACTERS IN THE TEXT STRING
C
C
C SPECIFICATION STATEMENTS
C
      CHARACTER*1 ITEXT(IMAX),ISP,ICOMMA
      DATA ISP/' '/,ICOMMA/','/
C
C CODE
C
      NFIELD = 0
      I = 0
 20   I = I + 1
      IF (I.GT.IMAX) RETURN
      IF (ITEXT(I).EQ.ISP .OR. ITEXT(I).EQ.ICOMMA) GO TO 20
      NFIELD = NFIELD + 1
 30   I = I + 1
      IF (I.GT.IMAX) RETURN
      IF (ITEXT(I).NE.ISP .AND. ITEXT(I).NE.ICOMMA) GO TO 30
      GO TO 20
      END
