**==aa0001.f  
C
C PROGRAM 'LAUENORM'   
C ==================
C
C JOHN W. CAMPBELL, APRIL 1986
C
C UPDATED TO DECONVOLUTE MULTIPLES USING WAVELENGTH NORMALISATION CURVE (1992)
C
C UPDATED TO PROCESS MORE THAN ONE FILM PACK, SEPTEMBER 1986
C OPTION FOR 'A' FILM DATA FROM GE1 FILES AS INPUT, MARCH 1990
C
C FORTRAN 77
C
C READ ONE OR MORE LAUE DATA FILES, SPLIT THE DATA INTO WAVELENGTH BINS,
C SCALE THESE TOGETHER USING REFLECTION OVERLAPS, CURVE FIT
C THE RESULTING SCALE FACTORS AND SCALE THE FILM PACKS. THE WAVELENGTH
C BIN SCALING AND FILM PACK SCALING PROCESS IS ITERATED. FINALLY
C THE LAUE DATA ARE NORMALISED AND OUTPUT.
C
C
C
C FORMAT OF REFLECTION ARRAY 'IA'   (6 WORDS PER REFLECTION)
C
C    WORD 1      2**16(H+128)+2**8(K+128)+L+128
C    WORD 2      10000.0*LAMBDA
C    WORD 3      10000.0*THETA           (THETA IN RADIANS)
C    WORD 4      AI                      (INTENSITY) (AS A REAL VALUE)
C    WORD 5      ISIG + 2**16*KSYM       (SIG(I) AND SYMMETRY NO.)
C    WORD 6      2**16*IB+2**6*IPACK+
C                2*IFILM+IREJ            (IB=LAMBDA BIN NUMBER,IPACK=
C                                         FILM PACK NUMBER, IFILM IS
C                                         THE NUMBER OF THE FILM IN
C                                         THE PACK, IREJ IS
C                                         A FLAG =0 TO INCLUDE DATA
C                                         IN SCALE CALCULATION, =1 TO
C                                         EXCLUDE DATA FROM SCALE
C                                         FACTOR CALCULATION)
C
C
C WHEN THE LAUE DATA ARE INITIALLY READ IN THEY ARE STORED IN
C IN THE ARRAY IA AS DESCRIBED ABOVE EXCEPT THAT WORD 5 CONTAINS
C ISIG ONLY AND WORD 6 CONTAINS IPACK & FILM NUMBER ONLY.
C
C ARRAY TO HOLD MULTIPLES DATA
C
C 6 WORDS/REFLECTION ('MUL' (AMUL))
C
C    WORD 1      2**16(HN128)+2**8(KN+128)+LN+128  (PACKED NODAL INDICES)
C    WORD 2      NODAL LAMBDA (REAL)
C    WORD 3      2**6*IPACK+2*IFILM  (cf SINGLES ARRAY)
C    WORD 4      INTENSITY  (REAL)       
C    WORD 5      ISIG + 2**16*KSYM       (SIG(I) AND SYMMETRY NO.)
C    WORD 6      2**20*MULT+2**10*MINHRM+MAXHRM  (MULTIPLICITY, MINIMUM AND
C                                                 MAXIMUM HARMONIC NUMBERS)
C
C WHEN THE LAUE DATA ARE INITIALLY READ IN THEY ARE STORED IN
C IN THE ARRAY MUL AS DESCRIBED ABOVE EXCEPT THAT WORD 5 CONTAINS
C ISIG ONLY, ALSO WORD 6 CONTAINS ONLY MULT AND MINHRM
C 
C
C PARAMETER STATEMENTS
C
C          MXL = MAXIMUM NUMBER OF LAUE RANGES.
C          MXB = MAXIMUM NUMBER OF LAMBDA BINS ALLOWED AND MAXIMUM
C                NUMBER OF FILM PACKS ALLOWED.
C          MXR = MAXIMUM NUMBER OF REFLECTIONS.
C          MXMUL = MAXIMUM NUMBER OF MULTIPLES WHICH MAY BE STORED
C
      PARAMETER(MXL=8,MXB=100,MXR=300000,MXMUL=100000)
C
C
C SPECIFICATION STATEMENTS
C ------------------------
C
      LOGICAL BAT(MXB),END,BREF,OVLP(MXB,MXB),ERR,CFLAG,AFSCAL
      LOGICAL CBFLAG,CLFLAG
      LOGICAL ABSNT
      LOGICAL IFMULT
      DIMENSION NREFB(MXB),YMIN(MXB),YMAX(MXB),IOVLAP(MXB,MXB)
      DIMENSION NREFP(MXB),FFBAR(MXB),NBAR(MXB),QSCALE(MXB)
      DIMENSION IA(6*MXR)
      DIMENSION MUL(6*MXMUL),AMUL(6*MXMUL)
      DIMENSION ICENT(7),IHKL(3),ABSL(6)
      DIMENSION QVAL(11)
      INTEGER ISORT(5)
      CHARACTER*10 SPGNAM
      CHARACTER*10 PGNAME
      CHARACTER*6 PGN
      CHARACTER*1 LTYPE
      CHARACTER*70 TITL
      CHARACTER*80 HSTRNG(2)
      CHARACTER*11 OVHDR(2)
      CHARACTER*80 ITEXT
      CHARACTER*30 LABLS(7)
      CHARACTER*1 CTYPS(7)
      CHARACTER*30 LABROT(7)
      CHARACTER*1 CROT(7)
      CHARACTER*20 HEADR
      CHARACTER*8 LAMCOD(4)
      CHARACTER*9 OPT(2)
      CHARACTER*5 SCODE(5)
      CHARACTER*7 LFIN
      CHARACTER*15 INTT(2)
      DIMENSION CELL(6)
      REAL RSYMM(4,4,96)
      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),CFLAG,
     *              CBFLAG(MXB),CLFLAG(MXL)
      COMMON /INCUR/QMINL(MXL),QMAXL(MXL),ICOEFF(MXL),Q(MXL,11),
     *              LAMTYP,LRANGE
      COMMON /BINS/NRANGE,NBTOT,NVAL(MXB),SCALE(MXB),BTEMP(MXB),
     *ALAM(MXB),CSCALE(MXB)
      COMMON /PACKS/NPACK,NPVAL(MXB),PACKSC(MXB),PACKTF(MXB),
     *CPSCAL(MXB),CPTFAC(MXB)
      COMMON /MULDAT/ NMULT
      COMMON /SYSA/ ISYS,KSYS(3)
      COMMON /SYMM1/PGNAME
      EQUIVALENCE (IAA,AAA)
      EQUIVALENCE (MUL(1),AMUL(1))
      DATA LABLS/'H','K','L','F','SIGF','DANO','SIGDANO'/
      DATA CTYPS/'H','H','H','F','Q','D','Q'/
      DATA LABROT/'H','K','L','M/ISYM','BATCH','I','SIGI'/
      DATA CROT/'H','H','H','Y','B','J','Q'/
      DATA BAT/100*.FALSE./,NREFB/100*0/
      DATA HEADR/'ITERATION NUMBER   $'/
      DATA LFIN/'LAUE000'/
      DATA OVHDR/'LAMBDA BINS','FILM  PACKS'/
      DATA LAMCOD/'INPUT','UNITY','PEALECT','YORK'/
      DATA OPT/'NORMALISE','SCALE'/
      DATA SCODE/'UNITY','IAVE','SYMM','SYMMB','INPUT'/
      DATA NTWO16/65536/
      DATA NTWO8/256/
      DATA NTWO6/64/
      DATA NTWO20/1048576/
      DATA NTWO10/1024/
      DATA INTT/'BOX','PROFILE FITTED'/
      DATA HSTRNG/'Created by LAUENORM',' '/
      EXTERNAL IBPACK
      INTEGER IBPACK
C
C INITIALISATIONS
C ---------------
C
      ALMIN = 100000.0
      ALMAX = 0.0
      DLAM = 0.0
      SCALIN = 1.0
      SCALOP = 150.0
      NABS = 0
      NREF = 0
      NSCAL = 0
      NRIN = 0
      NREAD = 0
      NRIDX = 0
      NREJ = 0
      NREJL = 0
      NNEG = 0
      NSREJ = 0
      NMINP = 0
      NMULT = 0
      AFSCAL = .TRUE.
      CFLAG = .FALSE.
      IFAIL = 0
      CALL CCPOPN(5,'DATA',5,1,LDUM,IFAIL)
      CALL CCPOPN(6,'PRINTER',6,1,LDUM,IFAIL)
      CALL CCPVRS(6,'LAUENORM','07/11/94')
      CALL MTZINI
C
C READ DATA CONTROL CARDS AND OPEN OUTPUT MTZ FILE
C ------------------------------------------------
C
C DATA CARD 1   TITLE FOR OUTPUT FILE
C
      READ (5,6010) TITL
      WRITE (6,6220) TITL
C
C DATA CARD 2   CELL DIMENSIONS
C
      CALL FFREAD(5,ITEXT,80,6,.TRUE.)
      DO 20 I = 1 , 6
         CALL FVALUE(CELL(I))
 20   CONTINUE
      CALL SETRSL(CELL(1),CELL(2),CELL(3),CELL(4),CELL(5),CELL(6))
C
C DATA CARD 3   PROGRAM OPTION, NUMBER OF FILM PACKS ['GE', INTYP, ISPAT]
C
      CALL FFREAD(5,ITEXT,80,6,.TRUE.)
      IOPT = 1
      CALL TEXLIM(M1,M2)
      IF (M1.GT.0) THEN
         IF (ITEXT(M1:M2).EQ.OPT(1)) THEN
            IOPT = 1
         ELSE IF (ITEXT(M1:M2).EQ.OPT(2)) THEN
            IOPT = 2
         ELSE
            WRITE (6,6750) ITEXT(M1:M2)
            STOP
         END IF
      END IF
      WRITE (6,6530) OPT(IOPT)
      CALL IVALUE(NPACK)
      IF (NPACK.LE.0) NPACK = 1
      IF (NPACK.GT.MXB) THEN
         WRITE (6,6700) MXB
         STOP
      ELSE
         WRITE (6,6310) NPACK
         CALL TEXLIM(M1,M2)
         IF (M1.GT.0) THEN
            IF (ITEXT(M1:M1).NE.'G' .AND. ITEXT(M1:M1).NE.'g') THEN
               WRITE (6,6690) ITEXT(M1:M2)
               STOP
            ELSE
               AFSCAL = .FALSE.
               CALL IVALUE(INTYP)
               IF (INTYP.LE.0 .OR. INTYP.GT.2) INTYP = 2
               WRITE (6,6440) INTT(INTYP)
               CALL IVALUE(ISPAT)
               IF (ISPAT.NE.1) ISPAT = 0
               IF (ISPAT.EQ.0) WRITE (6,6590)
               IF (ISPAT.EQ.1) WRITE (6,6600)
            END IF
         END IF
C
C DATA CARD 4   INITITIAL LAMBDA CURVE TYPE
C
         CALL FFREAD(5,ITEXT,80,6,.TRUE.)
         LAMTYP = 2
         CALL TEXLIM(M1,M2)
         IF (M1.GT.0) THEN
            DO 30 L = 1 , 4
               LAMTYP = L
               IF (ITEXT(M1:M2).EQ.LAMCOD(L)) GO TO 40
 30         CONTINUE
            WRITE (6,6710) ITEXT(M1:M2)
            STOP
         END IF
      END IF
 40   WRITE (6,6120) LAMCOD(LAMTYP)
      ITORD = 0
      IF (LAMTYP.NE.2) ITORD = 1
      IF (ITORD.EQ.0) THEN
         WRITE (6,6250)
      ELSE
         WRITE (6,6260)
      END IF
C
C DATA CARD(S) 5  INPUT SCALING OPTIONS
C
      CALL FFREAD(5,ITEXT,80,6,.TRUE.)
      ISCTYP = 1
      CALL TEXLIM(M1,M2)
      IF (M1.GT.0) THEN
         DO 50 I = 1 , 5
            IF (ITEXT(M1:M2).EQ.SCODE(I)) GO TO 60
 50      CONTINUE
         WRITE (6,6760) ITEXT(M1:M2)
         STOP
 60      ISCTYP = I
      END IF
      DO 70 N = 1 , NPACK
         PACKSC(N) = 1.0
         PACKTF(N) = 0.0
         CPSCAL(N) = 1.0
         CPTFAC(N) = 0.0
 70   CONTINUE
      WRITE (6,6540) SCODE(ISCTYP)
      IF (ISCTYP.EQ.2) THEN
         CALL FVALUE(SCLAM1)
         CALL FVALUE(SCLAM2)
         IF (SCLAM2.GT.0.0 .AND. SCLAM2.LE.SCLAM1) THEN
            WRITE (6,6720) SCLAM1 , SCLAM2
            STOP
         ELSE
            IF (SCLAM2.EQ.0.0) THEN
               SCLAM1 = 0.95
               SCLAM2 = 1.05
            END IF
            WRITE (6,6130) SCLAM1 , SCLAM2
         END IF
      END IF
      IF (ISCTYP.EQ.3 .OR. ISCTYP.EQ.4) THEN
         CALL FVALUE(DLAM)
         IF (DLAM.LE.0.0) DLAM = 0.1
         WRITE (6,6550) DLAM
         NABS = 0
         DO 80 I = 1 , 6
            CALL FVALUE(X)
            IF (X.GT.0.0) THEN
               NABS = NABS + 1
               ABSL(NABS) = X
            END IF
 80      CONTINUE
         IF (NABS.GT.0) WRITE (6,6560) (ABSL(N),N=1,NABS)
      END IF
      IF (ISCTYP.EQ.5) THEN
         CALL GTVALS(5,ITEXT,80,NPACK,PACKSC,6,.TRUE.)
         DO 90 N = 1 , NPACK
            CPSCAL(N) = PACKSC(N)
 90      CONTINUE
         CALL GTVALS(5,ITEXT,80,NPACK,PACKTF,6,.TRUE.)
         DO 100 N = 1 , NPACK
            CPTFAC(N) = PACKTF(N)
 100     CONTINUE
         WRITE (6,6170) (N,PACKSC(N),PACKTF(N),N=1,NPACK)
      END IF
C
C DATA CARD 6   SPACE GROUP NUMBER,
C               NO. RANGES, NUMBER OF ITERATIONS,
C               NUMBER OF CYCLES REFINEMENT OF LAMBDA BIN SCALE FACTORS
C               PER ITERATION, NO. REFERENCE BIN, NUMBER OF CYCLES OF
C               REFINEMENT OF FILM PACK SCALE FACTORS PER ITERATION, NO.
C               OF REFERENCE FILM PACK.
C
      CALL FFREAD(5,ITEXT,80,6,.TRUE.)
      CALL IVALUE(NSPGRP)
      IF (NSPGRP.LE.0) NSPGRP = 1
      CALL IVALUE(NRANGE)
      IF (NRANGE.LE.0) NRANGE = 1
      IF (NRANGE.GT.MXL) THEN
         WRITE (6,6630)
         STOP
      ELSE
         CALL IVALUE(NITER)
         IF (NITER.LT.1) NITER = 1
         CALL IVALUE(NBCYC)
         IF (NBCYC.LT.0) NBCYC = 0
         CALL IVALUE(NBFIX)
         IF (NBFIX.LE.0) NBFIX = 1
         CALL IVALUE(NPCYC)
         IF (NPCYC.LT.0) NPCYC = 0
         CALL IVALUE(NPFIX)
         IF (NPFIX.LE.0) NPFIX = 1
         WRITE (6,6030) NSPGRP , NRANGE
         IF (IOPT.EQ.1) THEN
            IF (NPACK.GT.1) THEN
               WRITE (6,6150) NITER , NBCYC , NBFIX , NPCYC , NPFIX
            ELSE
               WRITE (6,6150) NITER , NBCYC , NBFIX
            END IF
         END IF
C
C DATA CARD 7   WEIGHTS FLAG, CURVE FIT SUPPRESSION COUNT,B-REFINE
C               SUPPRESSION FLAG, SIGREJ, IRFAC, IR2, IR3, SIGSCA
C
         CALL FFREAD(5,ITEXT,80,6,.TRUE.)
         CALL IVALUE(IWSCAL)
         CALL IVALUE(NCURV)
         CALL IVALUE(IBREF)
         CALL FVALUE(SIGREJ)
         CALL IVALUE(IRFAC)
         CALL IVALUE(IR2)
         CALL IVALUE(IR3)
         CALL FVALUE(SIGSCA)
         IF (IWSCAL.EQ.0) THEN
            WRITE (6,6290)
         ELSE
            WRITE (6,6300)
         END IF
         WRITE (6,6160) NCURV
         IF (NPACK.GT.1 .AND. IBREF.GT.0) WRITE (6,6110) IBREF
         IF (SIGREJ.GT.0.0) WRITE (6,6380) SIGREJ
         IF (SIGSCA.GT.0) WRITE (6,6430) SIGSCA
C
C DATA CARDS 8  BMINL BMAXL BINS NFIT IFEX  [IORD Q(0) .. Q(IORD)]
C
C  * NRANGE
C
         DO 140 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 650
C
C SET INITIAL LAMBDA CURVE SCALING COEFFICIENTS
C
            IF (LAMTYP.EQ.1) THEN
               CALL TEXLIM(M1,M2)
               IF (M1.EQ.0) GO TO 640
               CALL DECNUM(ITEXT,M1,M2,1,XX,IERR)
               IF (IERR.GT.0) GO TO 640
               IORD = NINT(XX)
               ICOEFF(NN) = IORD + 1
               NFIT(NN) = IORD
               IF (IORD.LT.0 .OR. IORD.GT.10) GO TO 640
               DO 120 I = 1 , ICOEFF(NN)
                  CALL TEXLIM(M1,M2)
                  IF (M1.EQ.0) THEN
                     CALL GTVALS(5,ITEXT,80,IORD-I+2,QVAL,6,.TRUE.)
                     K = 0
                     DO 110 J = I , ICOEFF(NN)
                        K = K + 1
                        Q(NN,J) = QVAL(K)
 110                 CONTINUE
                     GO TO 130
                  END IF
                  CALL DECNUM(ITEXT,M1,M2,2,Q(NN,I),IERR)
                  IF (IERR.GT.0) GO TO 640
 120           CONTINUE
 130           QMINL(NN) = BMINL(NN)
               QMAXL(NN) = BMAXL(NN)
               LRANGE = NRANGE
            ELSE
               CALL SETCUR
            END IF
            CLFLAG(NN) = .TRUE.
            IF (NFIT(NN).LE.0) THEN
               NFIT(NN) = 0
               CLFLAG(NN) = .FALSE.
               IF (IFEX(NN).GT.0) GO TO 650
            END IF
            IF (NFIT(NN).GT.10) NFIT(NN) = 10
            WRITE (6,6060) NN
            WRITE (6,6070) BMINL(NN) , BMAXL(NN)
            IF (NFIT(NN).GT.0) THEN
               WRITE (6,6080) BINS(NN) , NFIT(NN)
            ELSE
               WRITE (6,6580) BINS(NN)
            END IF
            IF (MOD(IFEX(NN),2).NE.0) WRITE (6,6400)
            IF (IFEX(NN).GT.1) WRITE (6,6410)
            IF (LAMTYP.EQ.1) WRITE (6,6330) IORD , 
     +                                 (Q(NN,I),I=1,ICOEFF(NN)+1)
 140     CONTINUE
         IF (NRANGE.EQ.1 .AND. BMINL(1).EQ.0.0 .AND. BMAXL(1).EQ.0.0)
     +          THEN
            WRITE (6,6280)
            GO TO 160
         END IF
         DO 150 NN = 1 , NRANGE
            IF (BMINL(NN).LT.0.0) GO TO 590
            IF (BMAXL(NN).LT.0.0) GO TO 590
            IF (BMAXL(NN).LE.BMINL(NN)) GO TO 590
            IF (NN.NE.NRANGE) THEN
               IF (BMAXL(NN).GE.BMINL(NN+1)) GO TO 590
            END IF
 150     CONTINUE
      END IF
C
C READ SYMMETRY OPERATORS
C
 160  CALL MSYMLB (10,NSPGRP,SPGNAM,PGN,NSYMP,NSYMM,RSYMM)
      LTYPE = SPGNAM(1:1)
      PGNAME = PGN(3:)
      CALL SYMMAT
      KSYS(1) = 0
      KSYS(2) = 0
      KSYS(3) = 0
      ISYS = 0
      IF (LTYPE.EQ.'P') THEN
         ISYS = 0
      ELSE IF (LTYPE.EQ.'A') THEN
         ISYS = 2
         KSYS(2) = 1
         KSYS(3) = 1
      ELSE IF (LTYPE.EQ.'B') THEN
         ISYS = 2
         KSYS(1) = 1
         KSYS(3) = 1
      ELSE IF (LTYPE.EQ.'C') THEN
         ISYS = 2
         KSYS(1) = 1
         KSYS(2) = 1
      ELSE IF (LTYPE.EQ.'I') THEN
         ISYS = 2
         KSYS(1) = 1
         KSYS(2) = 1
         KSYS(3) = 1
      ELSE IF (LTYPE.EQ.'F') THEN
         ISYS = 4
      ELSE IF (LTYPE.EQ.'R') THEN
         ISYS = 3
         KSYS(1) = -1
         KSYS(2) = 1
         KSYS(3) = 1
      END IF
C
C DATA CARD 9   OUTPUT FILE FLAG, OUTPUT SCALE, CENTRIC FLAGS
C
      CALL FFREAD(5,ITEXT,80,6,.TRUE.)
      CALL IVALUE(IOUT)
      IF (IOUT.EQ.0) IOUT = 1
      IF (IOUT.EQ.-1) WRITE (6,6340)
      IF (IOUT.EQ.-2) WRITE (6,6420)
      IF (IOUT.EQ.1) WRITE (6,6350)
      IF (IOUT.EQ.2) WRITE (6,6360)
      IF (IOUT.GT.2) WRITE (6,6370) IOUT
      IF (IOUT.EQ.2) THEN
         DO 170 I = 1 , 7
            CALL IVALUE(ICENT(I))
 170     CONTINUE
         WRITE (6,6390) ICENT
      END IF
C
C DATA CARD 10 IFMUL, SCMAX, IWMUL, NEGTYP, MDIAG
C
      CALL FFREAD(5,ITEXT,80,6,.TRUE.)
      CALL IVALUE(IFMUL)
      IFMULT = .FALSE.
      IF (IFMUL.GT.0) IFMULT = .TRUE.
      CALL FVALUE(SCMAX)
      CALL IVALUE(IWMUL)
      CALL IVALUE(NEGTYP)
      CALL IVALUE(MDIAG)
      IF (MDIAG.GT.0) THEN
         IFAIL = 0
         CALL CCPOPN(7,'MULTDIAG',4,1,LDUM,IFAIL)
      END IF
C
C OPEN OUTPUT FILE(S)
C -------------------
C
      IF (IOUT.LT.0) THEN
         IFAIL = 0
         CALL CCPOPN(2,'SHELX',4,1,LDUM,IFAIL)
      END IF
      IF (IOUT.EQ.1 .OR. IOUT.EQ.2) THEN
         ISORT(1) = 1
         ISORT(2) = 2
         ISORT(3) = 3
         ISORT(4) = 0
         ISORT(5) = 0
         CALL LWOPEN (1,'HKLOUT')
         CALL LWTITL (1,TITL,0)
         IF (IOUT.EQ.1) THEN
            CALL LWCLAB (1,LABLS,5,CTYPS,0)
         ELSE
            CALL LWCLAB (1,LABLS,7,CTYPS,0)
         ENDIF
         CALL LWHIST (1,HSTRNG,1)
         CALL LWCELL (1,CELL)
         CALL LWSYMM (1,NSYMM,NSYMP,RSYMM,LTYPE,NSPGRP,SPGNAM,PGNAME)
         CALL LWSORT (1,ISORT)
      END IF
      IF (IOUT.GT.2) THEN
         ISORT(1) = 1
         ISORT(2) = 2
         ISORT(3) = 3
         ISORT(4) = 4
         ISORT(5) = 5
         CALL LWOPEN (1,'HKLOUT')
         CALL LWTITL (1,TITL,0)
         CALL LWCLAB (1,LABROT,7,CROT,0)
         CALL LWHIST (1,HSTRNG,1)
         CALL LWCELL (1,CELL)
         CALL LWSYMM (1,NSYMM,NSYMP,RSYMM,LTYPE,NSPGRP,SPGNAM,PGNAME)
         CALL LWSORT (1,ISORT)
         CALL LWBTIT (1,IOUT,TITL)
      END IF
C
C OPEN MULTIPLES OUTPUT FILE IF REQUIRED
C
      IF (IFMULT) THEN
         CALL LWOPEN (2,'HKLMULT')
         CALL LWTITL (2,TITL,0)
         CALL LWCLAB (2,LABLS,5,CTYPS,0)
         CALL LWHIST (2,HSTRNG,1)
         CALL LWCELL (2,CELL)
         CALL LWSYMM (2,NSYMM,NSYMP,RSYMM,LTYPE,NSPGRP,SPGNAM,PGNAME)
      END IF
C
C READ DATA FOR ALL LAUE FILM PACKS. STORE TEMPORARILY IN ARRAY 'IA'
C ------------------------------------------------------------------
C
C FIND MINIMUM AND MAXIMUM LAMBDA VALUES IN THE DATA
C --------------------------------------------------
C
      SIGMAX = 0.0
      DO 200 IP = 1 , NPACK
         FFBAR(IP) = 0.0
         NBAR(IP) = 0
         LUN = 20 + IP
         WRITE (LFIN(5:7),'(I3.3)') IP
         IFAIL = 0
         IF (AFSCAL) THEN
            CALL CCPOPN(LUN,LFIN,5,1,LDUM,IFAIL)
 180        READ (LUN,6020,END=190) IH , IK , IL , ALBDA , THETA , II , 
     +                              ISIG
            NREAD = NREAD + 1
            IF (SIGREJ.GT.0.0) THEN
               IF (II.LE.0) THEN
                  NNEG = NNEG + 1
                  GO TO 180
               END IF
               IF (II.LT.SIGREJ*ISIG) THEN
                  NSREJ = NSREJ + 1
                  GO TO 180
4              END IF
            ELSE
               IF (II.LE.0) NNEG = NNEG + 1
            END IF
            NRIN = NRIN + 1
            IF (NRIN.GT.MXR) GO TO 620
            JJ = 6*NRIN - 5
            IA(JJ) = NTWO16*(IH+128) + NTWO8*(IK+128) + IL + 128
            IA(JJ+1) = 10000.0*ALBDA
            IA(JJ+2) = 10000.0*THETA
            AAA = II
            IA(JJ+3) = IAA
            IA(JJ+4) = ISIG
            IA(JJ+5) = NTWO6*IP
            ALMAX = AMAX1(ALMAX,ALBDA)
            ALMIN = AMIN1(ALMIN,ALBDA)
            SIGMAX = AMAX1(FLOAT(ISIG),SIGMAX)
            GO TO 180
 190        CLOSE (UNIT=LUN)
         ELSE
            CALL CCPOPN(LUN,LFIN,5,4,19,IFAIL)
            CALL RDGE(LUN,IP,INTYP,SIGREJ,ISPAT,IA,NRIN,NREAD,NSREJ,
     +                NNEG,ALMIN,ALMAX,SIGMAX,MUL,NMINP)
            CLOSE (UNIT=LUN)
         END IF
 200  CONTINUE
      IF (SIGMAX.GT.32760.0) THEN
         SCALIN = 32760.0/SIGMAX
         WRITE (6,6510) SCALIN
      END IF
C
C CALCULATE BIN REQUIREMENTS
C
      NBTOT = 0
      IF (NRANGE.EQ.1 .AND. BMINL(1).EQ.0.0 .AND. BMAXL(1).EQ.0.0) THEN
         BMINL(1) = ALMIN
         BMAXL(1) = ALMAX
         IF (LAMTYP.EQ.1) THEN
            QMINL(1) = ALMIN
            QMAXL(1) = ALMAX
         END IF
      END IF
      DO 220 NN = 1 , NRANGE
         IF (BINS(NN).LE.0.0) GO TO 610
         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
         DO 210 INEW = 1 , NBIN(NN)
            IF (CLFLAG(NN)) THEN
               CBFLAG(NBTOT+INEW) = .TRUE.
            ELSE
               CBFLAG(NBTOT+INEW) = .FALSE.
            END IF
 210     CONTINUE
         NBTOT = NBTOT + NBIN(NN)
         IF (NBTOT.GT.MXB) GO TO 600
         WRITE (6,6090) NBIN(NN) , NN
 220  CONTINUE
      NPP = 0
      DO 230 NN = 1 , NRANGE
         NP(NN) = NPP
         NPP = NPP + NBIN(NN)
 230  CONTINUE
      WRITE (6,6040) NREAD
      WRITE (6,6050) ALMIN , ALMAX , NBTOT
      IF (NBTOT.GT.MXB) GO TO 600
      IF (NBFIX.GT.NBTOT) THEN
         WRITE (6,6680)
         STOP
      ELSE
         REWIND 1
         DO 240 N = 1 , NBTOT
            NVAL(N) = 0
            CSCALE(N) = 1.0
 240     CONTINUE
C
C 'RE-READ' LAUE DATA FROM TEMP STORAGE IN 'IA'  AND STORE IN ARRAY 'IA'
C ----------------------------------------------------------------------
C
         DO 270 NNI = 1 , NRIN
            IREJ = 0
            KK = 6*NNI - 5
            IPACK = IBPACK(2,IA(KK+5),ERR)
            IFILM = IBPACK(3,IA(KK+5),ERR)
            ALBDA = IA(KK+1)/10000.0
            IAA = IA(KK+3)
            AI = SCALIN*AAA
            ISIG = NINT(SCALIN*IA(KK+4))
            IF (ISIG.LT.1) ISIG = 1
            IF (ISIG.GT.32767) THEN
               NREJ = NREJ + 1
               GO TO 270
            END IF
            IF (SIGSCA.GT.0.0 .AND. AI.LT.SIGSCA*ISIG) IREJ = 1
            IF (AI.LE.0.0) IREJ = 1
            CALL GETHKL(IA(KK),IHKL)
            IH = IHKL(1)
            IK = IHKL(2)
            IL = IHKL(3)
            IF (PGNAME(1:1).EQ.'3' .OR. PGNAME(1:1).EQ.'6' .OR. 
     +          PGNAME(1:1).EQ.'4') THEN
               CALL TURN3(IH,IK,IL,KSYM)
            ELSE
               CALL TURNIP(IH,IK,IL,KSYM)
            END IF
            DO 250 NN = 1 , NRANGE
               IF (ALBDA.GE.BMINL(NN) .AND. ALBDA.LE.BMAXL(NN))
     +             GO TO 260
 250        CONTINUE
            NREJL = NREJL + 1
            GO TO 270
 260        IB = NBIN(NN)*(ALBDA-BMINL(NN))/(BMAXL(NN)-BMINL(NN)) + 1
            IF (IB.LT.1) IB = 1
            IF (IB.GT.NBIN(NN)) IB = NBIN(NN)
            NB = IB
            IB = IB + NP(NN)
            IF (.NOT.BAT(IB)) THEN
               YMIN(IB) = (NB-1)*(BMAXL(NN)-BMINL(NN))/NBIN(NN)
     +                    + BMINL(NN)
               YMAX(IB) = NB*(BMAXL(NN)-BMINL(NN))/NBIN(NN) + BMINL(NN)
               ALAM(IB) = 0.5*(YMIN(IB)+YMAX(IB))
               CSCALE(IB) = CURVIN(IB,ALAM(IB))
               SCALE(IB) = CSCALE(IB)
               BTEMP(IB) = 0.0
               IF (YMAX(IB).GT.BMAXL(NN)) YMAX(IB) = BMAXL(NN)
               BAT(IB) = .TRUE.
            END IF
            NREF = NREF + 1
            IF (IH.LT.-127 .OR. IH.GT.128) THEN
               NREF = NREF - 1
               NRIDX = NRIDX + 1
               GO TO 270
            END IF
            IF (IK.LT.-127 .OR. IK.GT.127) THEN
               NREF = NREF - 1
               NRIDX = NRIDX + 1
               GO TO 270
            END IF
            IF (IL.LT.-127 .OR. IL.GT.127) THEN
               NREF = NREF - 1
               NRIDX = NRIDX + 1
               GO TO 270
            END IF
            IF (IREJ.EQ.0) NSCAL = NSCAL + 1
            JJ = 6*NREF - 5
            IA(JJ) = NTWO16*(IH+128) + NTWO8*(IK+128) + IL + 128
            IA(JJ+1) = IA(KK+1)
            IA(JJ+2) = IA(KK+2)
            AAA = AI
            IA(JJ+3) = IAA
            IA(JJ+4) = NTWO16*KSYM + ISIG
            IA(JJ+5) = NTWO16*IB + NTWO6*IPACK + 2*IFILM + IREJ
            IF (IREJ.EQ.0) NVAL(IB) = NVAL(IB) + 1
            IF (IREJ.EQ.0) NPVAL(IPACK) = NPVAL(IPACK) + 1
            IF (ISCTYP.EQ.2) THEN
               IF (ALBDA.GE.SCLAM1 .AND. ALBDA.LE.SCLAM2) THEN
                  NBAR(IPACK) = NBAR(IPACK) + 1
                  FFBAR(IPACK) = FFBAR(IPACK) + AAA
               END IF
            END IF
            NREFB(IB) = NREFB(IB) + 1
            NREFP(IPACK) = NREFP(IPACK) + 1
 270     CONTINUE
C
C 'RE-READ' MULTIPLES LAUE DATA FROM TEMP STORAGE IN 'MUL'  AND STORE IN
C ----------------------------------------------------------------------
C ARRAY 'MUL'
C -----------
C
         IF (IFMULT.AND.NMINP.GT.0) THEN
            DO 275 NNI = 1 , NMINP
               KK = 6*NNI - 5
               IPACK = IBPACK(2,MUL(KK+4),ERR)
               IFILM = IBPACK(3,MUL(KK+4),ERR)
               IAA = MUL(KK+3)
               AI = SCALIN*AAA
               ISIG = NINT(SCALIN*MUL(KK+4))
               IF (ISIG.LT.1) ISIG = 1
               CALL GETHKL(MUL(KK),IHKL)
               IH = IHKL(1)
               IK = IHKL(2)
               IL = IHKL(3)
               IF (PGNAME(1:1).EQ.'3' .OR. PGNAME(1:1).EQ.'6' .OR. 
     +             PGNAME(1:1).EQ.'4') THEN
                   CALL TURN3(IH,IK,IL,KSYM)
               ELSE
                  CALL TURNIP(IH,IK,IL,KSYM)
               END IF
               NMULT = NMULT + 1
               JJ = 6*NMULT - 5
               MUL(JJ) = NTWO16*(IH+128) + NTWO8*(IK+128) + IL + 128
               MUL(JJ+1) = MUL(KK+1)
               MUL(JJ+2) = MUL(KK+2)
               AAA = AI
               MUL(JJ+3) = IAA
               MUL(JJ+4) = NTWO16*KSYM + ISIG
               MUL(JJ+5) = MUL(KK+5)
               CALL GETMUL (MUL(JJ+5),MULT,MINHRM,MAXHRM)
               NMUL = 1
               IHRM = MINHRM
 274           IHRM = IHRM + 1
               JJH = IHRM*IH
               JJK = IHRM*IK
               JJL = IHRM*IL
               IF (ISYS.NE.0) THEN
                  CALL SYSAB(JJH,JJK,JJL,ABSNT)
                  IF (ABSNT) GO TO 274
               END IF
               MAXHRM = IHRM
               NMUL = NMUL + 1
               IF (NMUL.LT.MULT) GO TO 274
               MUL(JJ+5) = NTWO20*MULT+NTWO10*MINHRM+MAXHRM
 275        CONTINUE
         END IF
C
C END OF DATA REACHED
C -------------------
C
         IF (SIGREJ.GT.0.0) THEN
            WRITE (6,6100) NREF , NSCAL , NREJ , NREJL , NRIDX , NNEG , 
     +                     SIGREJ , NSREJ
         ELSE
            WRITE (6,6102) NREF , NSCAL , NREJ , NREJL , NRIDX 
            WRITE (6,6103) NNEG
         END IF
         IF (IFMULT) WRITE (6,6105) NMULT
C
C PRINT BINS SUMMARY
C ------------------
C
         WRITE (6,6180) NBTOT
         CSS = CSCALE(NBFIX)
         DO 280 N = 1 , NBTOT
            CSCALE(N) = CSCALE(N)/CSS
            WRITE (6,6190) N , NREFB(N) , NVAL(N) , CSCALE(N) , YMIN(N)
     +                     , YMAX(N)
 280     CONTINUE
C
C PRINT FILM PACKS SUMMARY
C ------------------------
C
         IF (NPACK.GT.1) THEN
            WRITE (6,6230) NPACK
            DO 290 N = 1 , NPACK
               WRITE (6,6240) N , NREFP(N) , NPVAL(N)
 290        CONTINUE
         END IF
C
C FIND OVERLAPS
C -------------
C
         IF (NREF.EQ.0) STOP
         IF (NREF.GT.MXR) GO TO 620
         CALL REFSRT(IA,NREF,6,6*NREF)
         IF (ISCTYP.EQ.3 .OR. ISCTYP.EQ.4) CALL LAMSRT(IA,NREF)
         DO 390 IOV = 1 , 2
            NMAX = NBTOT
            IF (IOV.EQ.2) NMAX = NPACK
            DO 310 I = 1 , NMAX
               DO 300 J = 1 , NMAX
                  IOVLAP(I,J) = 0
 300           CONTINUE
 310        CONTINUE
            END = .FALSE.
            M1 = 1
            ILAST = IA(1)
            K = 0
 320        K = K + 1
            IF (K.GT.NREF) THEN
               END = .TRUE.
               GO TO 330
            END IF
            JJ = 6*K - 5
            IF (IA(JJ).EQ.ILAST) GO TO 320
 330        M2 = K - 1
            IF (M2.GT.M1) THEN
               DO 350 I = 1 , NMAX
                  DO 340 J = 1 , NMAX
                     OVLP(I,J) = .FALSE.
 340              CONTINUE
 350           CONTINUE
               DO 370 M = M1 , M2
                  I1 = IBPACK(IOV,IA(6*M),ERR)
                  IF (.NOT.(ERR)) THEN
                     MM1 = M + 1
                     IF (MM1.LE.M2) THEN
                        DO 360 MM = MM1 , M2
                           I2 = IBPACK(IOV,IA(6*MM),ERR)
                           IF (.NOT.(ERR)) THEN
                              IF (.NOT.OVLP(I1,I2)) THEN
                                 IOVLAP(I1,I2) = IOVLAP(I1,I2) + 1
                                 IF (I1.NE.I2) IOVLAP(I2,I1)
     +                               = IOVLAP(I1,I2)
                                 OVLP(I1,I2) = .TRUE.
                                 OVLP(I2,I1) = .TRUE.
                              END IF
                           END IF
 360                    CONTINUE
                     END IF
                  END IF
 370           CONTINUE
            END IF
            IF (END) THEN
C
C PRINT OVERLAPS TABLE
C --------------------
C
               WRITE (6,6450) OVHDR(IOV)
               WRITE (6,6460) (I,I=1,NMAX)
               DO 380 N = 1 , NMAX
                  WRITE (6,6470) N , (IOVLAP(N,K),K=1,NMAX)
 380           CONTINUE
            ELSE
               M1 = K
               ILAST = IA(JJ)
               GO TO 320
            END IF
 390     CONTINUE
C
C PRINT INITIAL R-FACTOR ANALYSES
C -------------------------------
C
         WRITE (6,6480)
         CALL RFACTS(IA,NREF,NBTOT,NPACK,.FALSE.,IR2,IR3)
C
C SORT MULTIPLES DATA ARRAY
C -------------------------
C
         IF (IFMULT) THEN
            CALL REFSRT(MUL,NMULT,6,6*NMULT)
            IF (MDIAG.GT.0) CALL MULDIA (MUL,NMULT)
         END IF
C
C WAVELENGTH BIN AND FILM PACK SCALING
C ====================================
C
C
C CALCULATE APPROXIMATE FILM PACK SCALE FACTORS BASED ON AVERAGE I
C ----------------------------------------------------------------
C VALUES IF REQUIRED
C ------------------
C
         IF (NPACK.GT.1 .AND. ISCTYP.EQ.2) THEN
            DO 400 I = 1 , NPACK
               IF (NBAR(I).EQ.0) GO TO 630
               FFBAR(I) = FFBAR(I)/NBAR(I)
 400        CONTINUE
            FFREF = FFBAR(NPFIX)
            DO 410 I = 1 , NPACK
               CPSCAL(I) = FFREF/FFBAR(I)
               CPTFAC(I) = 0.0
 410        CONTINUE
            WRITE (6,6140) (I,NBAR(I),CPSCAL(I),I=1,NPACK)
         END IF
C
C CALCULATE INITIAL SCALE FACTORS BETWEEN FILM PACKS BASED ON
C -----------------------------------------------------------
C EQUIVALENT REFLECTIONS MEASURED AT WAVELENGTHS WITHIN DLAM OF
C -------------------------------------------------------------
C EACH OTHER IF REQUIRED
C ----------------------
C
         IF (NPACK.GT.1 .AND. (ISCTYP.EQ.3 .OR. ISCTYP.EQ.4)) THEN
            WRITE (6,6570)
            BREF = .FALSE.
            IF (ISCTYP.EQ.4) BREF = .TRUE.
            CALL BSCALE(PACKSC,PACKTF,IOVLAP,NPACK,IA,NREF,NPCYC,BREF,
     +                  IWSCAL,NPFIX,0,2,DLAM,NABS,ABSL)
            DO 420 I = 1 , NPACK
               CPSCAL(I) = PACKSC(I)
               CPTFAC(I) = PACKTF(I)
 420        CONTINUE
         END IF
C
C ITERATIONS LOOP
C ===============
C
         BREF = .FALSE.
         DO 540 NITR = 1 , NITER
            IF (IBREF.GT.0 .AND. NITR.GE.IBREF) BREF = .TRUE.
            WRITE (HEADR(18:19),6320) NITR
            CALL HEADXX(6,HEADR,'$','1',2)
C
C REFINE INTER FILM PACK SCALE FACTORS HERE IF ITORD<>0
C -----------------------------------------------------
C
            IF (ITORD.NE.0 .AND. NPACK.GT.1) THEN
               DO 430 I = 1 , NPACK
                  PACKSC(I) = CPSCAL(I)
                  PACKTF(I) = CPTFAC(I)
 430           CONTINUE
               WRITE (6,6200)
               CALL BSCALE(PACKSC,PACKTF,IOVLAP,NPACK,IA,NREF,NPCYC,
     +                     BREF,IWSCAL,NPFIX,NITR,2,0.0,NABS,ABSL)
               DO 440 I = 1 , NPACK
                  CPSCAL(I) = PACKSC(I)
                  CPTFAC(I) = PACKTF(I)
 440           CONTINUE
               IF (IRFAC.GT.0) CALL RFACTS(IA,NREF,NBTOT,NPACK,.TRUE.,
     +             IR2,IR3)
            END IF
C
C REFINE INTER LAMBDA BIN SCALE FACTORS (OR USE LAMBDA CURVES FOR 'SCALE')
C ------------------------------------------------------------------------
C
            DO 450 I = 1 , NBTOT
               SCALE(I) = CSCALE(I)
 450        CONTINUE
            IF (IOPT.EQ.1) THEN
               WRITE (6,6210)
               CALL BSCALE(SCALE,BTEMP,IOVLAP,NBTOT,IA,NREF,NBCYC,
     +                     .FALSE.,IWSCAL,NBFIX,NITR,1,0.0,NABS,ABSL)
               IF (IRFAC.GT.0) THEN
                  DO 460 I = 1 , NBTOT
                     CSCALE(I) = SCALE(I)
 460              CONTINUE
                  CALL RFACTS(IA,NREF,NBTOT,NPACK,.TRUE.,IR2,IR3)
               END IF
               IF (NITR.GE.NCURV) THEN
                  CALL CURFIT
                  CFLAG = .TRUE.
                  CALL TABLE1
                  DO 470 I = 1 , NBTOT
                     IF (CBFLAG(I)) SCALE(I) = CURVAL(ALAM(I))
 470              CONTINUE
               END IF
            ELSE IF (IOPT.EQ.2) THEN
               WRITE (6,6270)
               CFLAG = .TRUE.
               DO 490 NN = 1 , NRANGE
                  NCOEFF(NN) = NFIT(NN) + 1
                  DO 480 I = 1 , NCOEFF(NN)
                     C(NN,I) = Q(NN,I)
 480              CONTINUE
 490           CONTINUE
               DO 500 I = 1 , NBTOT
                  IF (CBFLAG(I)) SCALE(I) = CURVAL(ALAM(I))
 500           CONTINUE
            END IF
            DO 510 I = 1 , NBTOT
               CSCALE(I) = SCALE(I)
 510        CONTINUE
            IF (ITORD.EQ.0) THEN
               IF ((NITR.GE.NCURV .AND. IRFAC.GT.0) .OR. NPACK.EQ.1)
     +             CALL RFACTS(IA,NREF,NBTOT,NPACK,.TRUE.,IR2,IR3)
            ELSE
               IF (NITR.GE.NCURV .OR. IRFAC.EQ.0 .OR. IOPT.EQ.2)
     +             CALL RFACTS(IA,NREF,NBTOT,NPACK,.TRUE.,IR2,IR3)
            END IF
C
C REFINE INTER FILM PACK SCALE FACTORS HERE IF ITORD=0
C ----------------------------------------------------
C
            IF (NPACK.NE.1) THEN
               IF (ITORD.EQ.0) THEN
                  DO 520 I = 1 , NPACK
                     PACKSC(I) = CPSCAL(I)
                     PACKTF(I) = CPTFAC(I)
 520              CONTINUE
                  WRITE (6,6200)
                  CALL BSCALE(PACKSC,PACKTF,IOVLAP,NPACK,IA,NREF,NPCYC,
     +                        BREF,IWSCAL,NPFIX,NITR,2,0.0,NABS,ABSL)
                  DO 530 I = 1 , NPACK
                     CPSCAL(I) = PACKSC(I)
                     CPTFAC(I) = PACKTF(I)
 530              CONTINUE
                  CALL RFACTS(IA,NREF,NBTOT,NPACK,.TRUE.,IR2,IR3)
               END IF
            END IF
 540     CONTINUE
C
C PRINT TABLE OF CUMULATIVE INVERSE LAMBDA BIN SCALE FACTORS
C ----------------------------------------------------------
C
         QSMAX = -100000.0
         DO 550 I = 1 , NBTOT
            QSCALE(I) = 1.0/CSCALE(I)
            IF (QSCALE(I).GT.QSMAX) QSMAX = QSCALE(I)
 550     CONTINUE
         CALL HEADXX(6,'NORMALISED INVERSE LAMBDA CURVE$','$','1',1)
         WRITE (6,6490)
         DO 560 I = 1 , NBTOT
            QSCALE(I) = QSCALE(I)/QSMAX
            WRITE (6,6500) I , ALAM(I) , QSCALE(I)
 560     CONTINUE
         WRITE (6,'('' '')')
C
C SCALE DATA AND WRITE OUTPUT FILE
C --------------------------------
C
         DO 570 NR = 1 , NREF
            JJ = 6*NR - 5
            ALBDA = IA(JJ+1)/10000.0
            THETA = IA(JJ+2)/10000.0
            IBIN = IBPACK(1,IA(JJ+5),ERR)
            IPACK = IBPACK(2,IA(JJ+5),ERR)
            SC = SCFAC(IBIN,IPACK,ALBDA,THETA)
            IAA = IA(JJ+3)
            AAA = AAA*SC
            IA(JJ+3) = IAA
            KSYM = IA(JJ+4)/NTWO16
            ISIG = IA(JJ+4) - NTWO16*KSYM
            ISIG = NINT(SC*ISIG)
            IF (ISIG.EQ.0) ISIG = 1
            IF (ISIG.GT.32767) ISIG = 32767
            IA(JJ+4) = NTWO16*KSYM + ISIG
 570     CONTINUE
C
C WRITE OUTPUT FILE
C -----------------
C
         CALL OPFILE(IOUT,IA,NREF,ICENT,SCALOP)
C
C DECONVOLUTE HARMONICS
C ---------------------
C
         IF (IFMULT) CALL DECONV(MUL,NMULT,SCMAX,IWMUL,NEGTYP,MDIAG)
         STOP
      END IF
 590  WRITE (6,6640) (BMINL(NN),BMAXL(NN),NN=1,NRANGE)
      STOP
 600  WRITE (6,6650) MXB
      STOP
 610  WRITE (6,6660) BINS(NN) , NN
      STOP
 620  WRITE (6,6670) MXR
      STOP
 630  WRITE (6,6730) I
      STOP
 640  WRITE (6,6740)
      STOP
 650  WRITE (6,6770) NN
      STOP
C
C FORMAT STATEMENTS
C -----------------
C
 6010 FORMAT (A)
 6020 FORMAT (3I5,2F10.5,2I10)
 6030 FORMAT (/,' SPACE GROUP NUMBER =',I4,//,
     +        ' NUMBER OF LAMBDA RANGES = ',I1)
 6040 FORMAT (//,' NUMBER OF REFLECTIONS READ FROM INPUT LAUE FILES =',
     +        I6)
 6050 FORMAT (//,' FROM INPUT LAUE FILE LAMBDA RANGE =',F6.3,' TO',F6.3,
     +        //,' NUMBER OF BATCHES =',I4)
 6060 FORMAT (////,' WAVELENGTH RANGE',I2)
 6070 FORMAT (//,' MINIMUM LAMBDA =',F5.3,'A  MAXIMUM LAMBDA =',F5.3,
     +        'A')
 6080 FORMAT (//,' BINS PARAMETER =',F7.3,///,
     +        ' ORDER FOR CURVE FITTING =',I3)
 6090 FORMAT (//,' NUMBER OF LAMBDA BINS TO BE USED FOR NORMALISATION ='
     +        ,I4,' FOR RANGE',I2,/)
 6100 FORMAT (//,' NUMBER OF REFLECTIONS ACCEPTED             =',I6,/,
     +        ' NUMBER OF REFLECTIONS FOR SCALING          =',I6,/,
     +        ' NUMBER REJECTED (SIGMA > 32767 )           =',I6,/,
     +        ' NUMBER REJECTED (OUTSIDE OF LAMBDA RANGES) =',I6,/,
     +        ' NUMBER REJECTED (OUTSIDE INDEX LIMITS)     =',I6,/,
     +        ' NUMBER REJECTED (NON-POSITIVE INTENSITIES) =',I6,/,
     +        ' NUMBER REJECTED (LESS THAN',F6.2,' SIGMA)    =',I6)
 6102 FORMAT (//,' NUMBER OF REFLECTIONS ACCEPTED             =',I6,/,
     +        ' NUMBER OF REFLECTIONS FOR SCALING          =',I6,/,
     +        ' NUMBER REJECTED (SIGMA > 32767 )           =',I6,/,
     +        ' NUMBER REJECTED (OUTSIDE OF LAMBDA RANGES) =',I6,/,
     +        ' NUMBER REJECTED (OUTSIDE INDEX LIMITS)     =',I6)
 6103 FORMAT (//' NUMBER OF NON-POSITIVE INTENSITIES INCLUDED =',I6)
 6105 FORMAT (//,' NUMBER OF MULTIPLES STORED = ',I6)
 6110 FORMAT (/,' REFINE FILM PACK TEMPERATURE FACTORS FROM ITERATION',
     +        ' NUMBER',I3)
 6120 FORMAT (//,' INITIAL LAMBDA CURVE TYPE: ',A)
 6130 FORMAT (//,' CALCULATE INITIAL FILM PACK SCALE FACTORS FROM',
     +        ' AVERAGE INTENSITIES IN THE LAMBDA RANGE',F7.3,' TO',
     +        F7.3)
 6140 FORMAT (//,' INITIAL FILM PACK SCALING FACTORS FROM AVERAGE',
     +        ' INTENSITIES',//,'   PACK   NO. OF    SCALE',/,
     +        '    NO.   REFLNS    FACTOR',/,(/,I7,I9,F10.3))
 6150 FORMAT (/,' NUMBER OF ITERATIONS OF SCALING PROCEDURE    =',I3,/,
     +        /,' FOR LAMBDA BIN SCALING:',/,
     +        ' NUMBER OF CYCLES OF REFINEMENT PER ITERATION =',I3,/,
     +        ' NUMBER OF REFERENCE BIN FOR SCALING          =',I3,:,/,
     +        /,' FOR FILM PACK SCALING:',/,
     +        ' NUMBER OF CYCLES OF REFINEMENT PER ITERATION =',I3,/,
     +        ' NUMBER OF REFERENCE PACK FOR SCALING         =',I3,/)
 6160 FORMAT (/,' FIT CURVES TO THE WAVELENGTH BIN SCALE FACTORS',
     +        ' FROM ITERATION NUMBER',I3)
 6170 FORMAT (//,' INPUT FILM PACK SCALE AND TEMPERATURE FACTORS',//,
     +        '  PACK    SCALE   TEMPERATURE',/,
     +        '   NO.    FACTOR     FACTOR',/,(/,' ',I5,F9.3,F12.3))
 6180 FORMAT (//,' TOTAL NUMBER OF BINS FOR SCALING =',I4,//,
     +        ' SUMMARY OF WAVELENGTH BINS DATA:',//,
     +        ' IBIN  NO.REFLNS NO._FOR_SCALING   SCALE  RANGE',/)
 6190 FORMAT (' ',I4,I11,I16,F8.3,'  WAVELENGTH RANGE',F7.3,' TO',F7.3)
 6200 FORMAT (//,' REFINE FILM PACK SCALE FACTORS')
 6210 FORMAT (//,' REFINE LAMBDA BIN SCALE FACTORS')
 6220 FORMAT (//,' TITLE: ',A)
 6230 FORMAT (//,' TOTAL NUMBER OF FILM PACKS FOR SCALING =',I4,//,
     +        ' SUMMARY OF FILM PACKS DATA:',//,
     +        ' IPACK  NO.REFLNS NO._FOR_SCALING ',/)
 6240 FORMAT (' ',I5,I11,I16)
 6250 FORMAT (//,' DURING ITERATIONS REFINE NORMALISATION CURVE THEN',
     +        ' FILM PACK SCALE FACTORS')
 6260 FORMAT (//,' DURING ITERATIONS REFINE FILM PACK SCALE FACTORS',
     +        ' THEN NORMALISATION CURVE')
 6270 FORMAT (//,
     +        ' APPLY PACK SCALE FACTORS AND INPUT NORMALISATION CURVES'
     +        )
 6280 FORMAT (//,' LAMBDA LIMITS WILL BE TAKEN FROM INPUT LAUE FILE')
 6290 FORMAT (/,' USE SIGMAS IN WEIGHTING REFLECTIONS FOR SCALE FACTOR',
     +        ' REFINEMENTS')
 6300 FORMAT (/,' USE UNIT WEIGHTS IN WEIGHTING REFLECTIONS FOR SCALE',
     +        ' FACTOR REFINEMENTS')
 6310 FORMAT (/,' NUMBER OF FILM PACKS =',I4)
 6320 FORMAT (I2)
 6330 FORMAT (//,' POLYNOMIAL COEFFICIENTS   Q 0  TO  Q',I2,//,
     +        (' ',6E20.8))
 6340 FORMAT (/,' OUTPUT RELECTION FILE WILL CONTAIN MERGED REFLECTION',
     +        ' DATA (F AND SIGF) IN SHELX FORMAT')
 6350 FORMAT (/,' OUTPUT REFLECTION FILE WILL CONTAIN MERGED',
     +        ' REFLECTION DATA WITHOUT ANOMALOUS DATA')
 6360 FORMAT (/,' OUTPUT REFLECTION FILE WILL CONTAIN MERGED',
     +        ' REFLECTION DATA INCLUDING ANOMALOUS DATA')
 6370 FORMAT (/,' OUTPUT REFLECTION FILE WILL CONTAIN UNMERGED',
     +        ' REFLECTION DATA WITH A BATCH NUMBER OF',I6)
 6380 FORMAT (/,' REJECT REFLECTIONS WITH INTENSITIES <',F6.3,' SIGMA')
 6390 FORMAT (/,' CENTRIC ZONES FLAGS:',7I2)
 6400 FORMAT (/,' ADD EXTRA EXTRAPOLATED POINT AT  LOW WAVELENGTH ',
     +        'END OF RANGE BEFORE CURVE FITTING')
 6410 FORMAT (/,' ADD EXTRA EXTRAPOLATED POINT AT HIGH WAVELENGTH ',
     +        'END OF RANGE BEFORE CURVE FITTING')
 6420 FORMAT (/,' OUTPUT REFLECTION FILE WILL CONTAIN UNMERGED',
     +        ' REFLECTION DATA (F, SIGF, LAMBDA) IN SHELX FORMAT')
 6430 FORMAT (/,' EXCLUDE REFLECTIONS WITH I <',F5.1,'*SIG(I) FROM',
     +      ' CALCULATING THE LAMBDA CURVE AND INTER-PACK SCALE FACTORS'
     +      )
 6440 FORMAT (/,' ''A'' FILM DATA INPUT FROM GE1 FILE(S) - INTEGRATION',
     +        ' TYPE: ',A)
 6450 FORMAT (///,' ',A,' OVERLAP ANALYSIS',//)
 6460 FORMAT (10X,20I6)
 6470 FORMAT (/,' ',I5,4X,20I6,(/,10X,20I6))
 6480 FORMAT (///,' R-FACTOR ANALYSES BEFORE SCALING')
 6490 FORMAT (//,' BIN   LAMBDA   SCALE',/)
 6500 FORMAT (' ',I3,F9.3,F8.4)
 6510 FORMAT (//,' *WARNING* PROGRAM RESCALED INPUT LAUE VALUES BY',
     +        F8.3)
 6520 FORMAT (//,' *WARNING* PROGRAM RESCALED OUTPUT INTENSITIES BY',
     +        F8.3)
 6530 FORMAT (///,' PROGRAM OPTION SELECTED: ',A)
 6540 FORMAT (//,' INPUT SCALING OPTION: ',A)
 6550 FORMAT (//,' DELTA LAMBDA FOR ACCEPTING REFLECTIONS FOR INITIAL',
     +        ' SCALING =',F6.2)
 6560 FORMAT (//,' ABSORPTION EDGES:',6F8.3)
 6570 FORMAT ('1',///,' CALCULATION OF INITIAL FILM PACK SCALING FROM',
     +        ' EQUIVALENT REFLECTIONS MEASURED WITHIN DELTA LAMBDA')
 6580 FORMAT (//,' BINS PARAMETER =',F7.3,///,
     +        ' NO CURVE FITTING FOR THIS RANGE')
 6590 FORMAT (/,' SPATIAL OVERLAP DATA NOT INCLUDED')
 6600 FORMAT (/,' SPATIAL OVERLAP DATA INCLUDED')
C
 6630 FORMAT (/,' **MAXIMUM OF ',I2,' RANGES ALLOWED**')
 6640 FORMAT (/,' **LAMBDA RANGE LIMITS MUST BE IN INCREASING POSITIVE',
     +        ' VALUES OF LAMBDA**:',12F7.3)
 6650 FORMAT (/,' **MAX OF ',I3,' BINS ALLOWED**')
 6660 FORMAT (/,' **BINS PARAMETER MUST BE GREATER THAN 0, BINS =',
     +        F10.5,' FOR RANGE',I2,' **')
 6670 FORMAT (/,' **MAXIMUM OF',I6,' REFLECTIONS EXCEEDED**')
 6680 FORMAT (/,' **BIN TO BE ASSIGNED A SCALE FACTOR OF 1.0 IS',
     +        ' ABSENT**')
 6690 FORMAT (/,' **UNKNOWN OPTION: ',A,' **')
 6700 FORMAT (/,' **MAXIMUM OF ',I3,' FILM PACKS ALLOWED**')
 6710 FORMAT (/,' **LAMBDA CURVE TYPE ',A,' INVALID**')
 6720 FORMAT (/,' **SCLAM1 > SCLAM2**:',2F6.2)
 6730 FORMAT (/,' **NO REFLECTION DATA FOR SCALING BETWEEN SCLAM1',
     +        ' AND SCLAM2 FOR FILM PACK',I3,' **')
 6740 FORMAT (/,' **INVALID CURVE FITTING PARAMETER SPECIFICATION**')
 6750 FORMAT (/,' **PROGRAM OPTION ',A,' INVALID**')
 6760 FORMAT (/,' **INPUT SCALING TYPE ',A,' INVALID**')
 6770 FORMAT (/,' **INVALID VALUE OF IFEX FOR RANGE',I2,' **')
      END
**==setcur.f  
C
C SUBROUTINE 'SETCUR'
C ===================
C
C SET INITIAL NORMALISATION CURVE DEPENDING ON OPTION SELECTED
C
      SUBROUTINE SETCUR
C
C PARAMETERS
C
C NONE
C
      PARAMETER(MXL=8)
C
C SPECIFICATION STATEMENTS
C ------------------------
C
      DIMENSION KCOEFF(3,2),RMIN(3,2),RMAX(3,2),R(11,3,2)
      COMMON /INCUR/QMINL(MXL),QMAXL(MXL),ICOEFF(MXL),Q(MXL,11),
     *              LAMTYP,LRANGE
      DATA KCOEFF
     P           /  4,  5,  6,
     Y              5,  5,  5/
      DATA RMIN
     P           / 0.450, 0.490, 0.920,
     Y             0.300, 0.490, 0.920/
      DATA RMAX
     P           / 0.490, 0.920, 2.600,
     Y             0.490, 0.920, 1.500/
      DATA R
     P           /  -803.23, 3562.90,-3843.70,  -10.17,   0.000, 6*0,
     P               386.64,-1623.42, 2568.06,-1802.63,  472.85, 6*0,
     P              -120.16,  489.63, -760.96,  573.90, -211.65,
     P                30.69,5*0,
     Y               0.1978E05,-0.1797E06, 0.6304E06,-0.1004E07,
     Y               0.6082E06,6*0,
     Y               0.4578E03,-0.1604E04, 0.2032E04,-0.1049E04,
     Y               0.1640E03,6*0,
     Y               0.3358E03,-0.9594E03, 0.1034E04,-0.4967E03,
     Y               0.8977E02,6*0/
C
C SET CURVE PARAMETERS FOR INITIAL CURVE
C --------------------------------------
C
      IF (LAMTYP.GT.2) THEN
C
C PEA LECTIN OR YORK CURVE
C ------------------------
C
         LRANGE = 3
         II = LAMTYP - 2
         DO 30 NN = 1 , 3
            ICOEFF(NN) = KCOEFF(NN,II)
            QMINL(NN) = RMIN(NN,II)
            QMAXL(NN) = RMAX(NN,II)
            DO 20 I = 1 , 11
               Q(NN,I) = R(I,NN,II)
 20         CONTINUE
 30      CONTINUE
         RETURN
      ELSE
C
C UNITY
C -----
C
         LRANGE = 1
         ICOEFF(1) = 1
         Q(1,1) = 1.0
         DO 40 I = 2 , 11
            Q(1,I) = 0.0
 40      CONTINUE
         QMINL(1) = 0.0
         QMAXL(1) = 10000.0
         RETURN
      END IF
      END
**==curvin.f  
C
C FUNCTION 'CURVIN'
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. THIS FUNCTION IS FOR THE INITIAL SCALING CURVES
C OF LAUENORM
C
      FUNCTION CURVIN(IB,X)
C
C PARAMETERS
C
C          IB (I)   THE LAMBDA BIN NUMBER
C           X (I)   THE LAMBDA VALUE
C
C RETURNS   CURVIN, THE SCALING FUNCTION VALUE.
C
      PARAMETER(MXL=8)
C
C SPECIFICATION STATEMENTS
C ------------------------
C
      COMMON /INCUR/QMINL(MXL),QMAXL(MXL),ICOEFF(MXL),Q(MXL,11),
     *              LAMTYP,LRANGE
C
C CALCULATE REQUIRED SCALING VALUE
C --------------------------------
C
      NN = 1
      IF (X.GE.QMINL(1)) THEN
         DO 20 N = 1 , LRANGE
            NN = N
            IF (X.GE.QMINL(NN) .AND. X.LE.QMAXL(NN)) GO TO 30
 20      CONTINUE
      END IF
      WRITE (6,6010) IB
 30   CURVIN = Q(NN,ICOEFF(NN))
      DO 40 J = ICOEFF(NN) - 1 , 1 , -1
         CURVIN = CURVIN*X + Q(NN,J)
 40   CONTINUE
      RETURN
C
C FORMAT STATEMENTS
C -----------------
C
 6010 FORMAT (/,' *WARNING* AVERAGE LAMBDA VALUE FOR BIN',I3,
     +        ' IS OUTSIDE RANGE OF REFERENCE CURVES')
      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
**==lamsrt.f  
C
C SUBROUTINE 'LAMSRT'
C ===================
C
C SORT EQUIVALENTS INTO LAMBDA ORDER WITHIN LAUENORM 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 = 6*K - 5
      IF (IA(JJ).EQ.ILAST) GO TO 20
 30   M2 = K - 1
      MM = 6*M1 - 11
      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 + 6*L
               L2 = L1 + 6
               IF (IA(L1+1).GT.IA(L2+1)) THEN
                  SWAP = .TRUE.
                  DO 40 II = 0 , 5
                     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
**==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  (ASSUMES 'NWRD' WORDS/REFLECTION)
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)   NUMBER OF WORDS/REFLECTION IN ARRAY
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-1) + 1
      J2 = NWRD*(I2-1) + 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
**==bscale.f  
C
C SUBROUTINE 'BSCALE'
C ===================
C
C SCALE A SET OF BATCHES OF DATA TOGETHER USING OVERLAPPING
C REFLECTIONS. THE CODE IS BASED ON CODE TAKEN FROM THE PROGRAM
C ROTAVATA
C
      SUBROUTINE BSCALE(SCALE,BTEMP,IOVLAP,NBTOT,IA,NREF,NCYC,BREF,IWT
     *,JUNITY,NITR,ITYP,DLAM,NABS,ABSL)
C
C PARAMETERS
C
C       SCALE (I/O) ARRAY HOLDING SCALE FACTORS FOR THE BATCHES
C                   (DIMENSIONED TO AT LEAST NBTOT)
C       BTEMP (I/O) ARRAY HOLDING TEMPERATURE FACTORS FOR THE BATCHES
C                   (DIMENSIONED TO AT LEAST NBTOT)
C      IOVLAP (O)   ARRAY RETURNING NUMBER OF OVERLAPS BETWEEN ALL
C                   PAIRS OF BINS AS USED IN THE SCALING PROCESS
C                   (DIMENSIONED NBTOT*NBTOT)
C       NBTOT (I)   THE NUMBER OF BINS TO BE SCALED
C          IA (I)   ARRAY HOLDING LAUE REFLECTION DATA (SEE LAUENORM
C                   MAIN PROGRAM FOR FORMAT)
C        NREF (I)   THE NUMBER OF REFLECTIONS IN ARRAY 'IA'
C        NCYC (I)   THE NUMBER OF REFINEMENT CYCLES
C        BREF (I)   LOGICAL FLAG =.TRUE.  REFINE TEMPERATURE FACTORS
C                                =.FALSE. DO NOT
C         IWT (I)   FLAG = 0  WEIGHTING OF 1/SD*SD
C                        = 1  CONSTANT WEIGHTS OF 0.001
C      JUNITY (I)   NUMBER OF THE BIN WHOSE SCALE FACTOR IS TO BE SET
C                   TO 1.0
C        NITR (I)   ITERATION NUMBER (FOR OUTPUT ANNOTATION)
C        ITYP (I)   FLAG =1 LAMBDA BIN SCALING
C                        =2 FILM PACK SCALING
C        DLAM (I)   IF DLAM > 0.0 THEN TREAT REFLECTION EQUIVALENTS,
C                   MEASURED AT WAVELENGTHS MORE THAN DLAM DIFFERENT
C                   OR AT OPPOSITE SIDES OF AN ABSORPTION EDGE, AS
C                   IF THEY WERE SEPARATE REFLECTIONS. FOR THIS OPTION
C                   THE EQUIVALENTS FOR A REFLECTION MUST BE SORTED
C                   BY WAVELENGTH.
C        NABS (I)   NUMBER OF ABSORPTION EDGES
C        ABSL (I)   ARRAY CONTAINING NABS ABSORPTION EDGE WAVELENGTHS
C
C PARAMETER STATEMENTS
C
C           MXB = MAXIMUM NUMBER OF LAMBDA BINS ALLOWED
C           MXD = MAXIMUM DIMENSION FOR MATRIX ARRAYS (2*MXB)
C           MXM = MAXIMUM NUMBER OF MEASUREMENTS/REFLECTION
C           MXR = MAXIMUM NUMBER OF REFLECTIONS
C
      PARAMETER (MXB=100,MXD=2*MXB,MXM=40,MXR=300000)
C
C SPECIFICATION STATEMENTS
C ------------------------
C
      LOGICAL ISTOR,IEND,BREF,ERR
      DIMENSION SCALE(NBTOT),BTEMP(NBTOT),IOVLAP(NBTOT,NBTOT)
      DIMENSION SCX(MXB),BTX(MXB)
      DIMENSION IA(6*NREF)
      DIMENSION INTER(MXB),ISER(MXB),SDK(MXB),SDB(MXB)
      DIMENSION COEFC(MXD),VAL(MXD),COEFS(MXD,MXD),VEC(MXD,MXD)
      DIMENSION WORK(8*MXD),DP(MXD)
      DIMENSION FF(MXM),WSD(MXM),EX(MXM),GG(MXM),NO(MXM)
      DIMENSION ABSL(*)
      EQUIVALENCE (IAA,AAA)
      DATA NTWO16/65536/
C
C INITIALISATIONS
C ---------------
C
      IF (NCYC.LE.0) RETURN
      DO 20 I = 1 , NBTOT
         SCX(I) = 1.0
         BTX(I) = 0.0
 20   CONTINUE
      NBLOCK = MXD
      NPARMS = 2*NBTOT
      IF (.NOT.BREF) NPARMS = NBTOT
C
C LOOP THROUGH REFINEMENT CYCLES
C ------------------------------
C
      DO 350 ICYC = 1 , NCYC
C
C SCALE DOWN BTX BY 50.0
C
         DO 30 K = 1 , NBTOT
            BTX(K) = 0.02*BTX(K)
 30      CONTINUE
C
C INITIALIZE FOR EVERY BLOCK
C
         DO 40 N = 1 , NBTOT
            INTER(N) = 0
 40      CONTINUE
         DO 60 N = 1 , NBLOCK
            COEFC(N) = 0.
            VAL(N) = 0.
            DO 50 J = 1 , NBLOCK
               COEFS(J,N) = 0.
               VEC(J,N) = 0.
 50         CONTINUE
 60      CONTINUE
         DO 80 I = 1 , NBTOT
            DO 70 J = 1 , NBTOT
               IOVLAP(J,I) = 0
 70         CONTINUE
 80      CONTINUE
         ISTOR = .FALSE.
         IEND = .FALSE.
         NR = 0
         ND2 = 0
         WD2 = 0
         KOUNTA = 0
         FBAR = 0.0
C
C ACCESS NEXT REFLECTION
C
 90      NR = NR + 1
         IF (NR.GT.NREF) THEN
C
C END OF REFLECTIONS, SET EOF FLAG AND PROCESS LAST REFLECTION
C
            IEND = .TRUE.
            IF (.NOT.ISTOR) GO TO 200
            IF (MSER.LE.1 .OR. INEG.NE.0) GO TO 200
            GO TO 160
         ELSE
            JJ = 6*NR - 5
            KHKL = IA(JJ)
            ALBDA = IA(JJ+1)/10000.0
            THETA = IA(JJ+2)/10000.0
            IB = IBPACK(1,IA(JJ+5),ERR)
            IF (ERR) GO TO 90
            IPACK = IBPACK(2,IA(JJ+5),ERR)
            SC = SCFAC(IB,IPACK,ALBDA,THETA)
            IAA = IA(JJ+3)
            RMODF = AAA*SC
            KSYM = IA(JJ+4)/NTWO16
            SDIN = IA(JJ+4) - NTWO16*KSYM
            SDIN = SDIN*SC
            N1 = IBPACK(ITYP,IA(JJ+5),ERR)
C
C BEGIN TO PROCESS REFLECTION JUST ACCESSED.
C IF THIS IS A NEW HKL DEAL WITH PREVIOUS HKL FIRST
C AND THEN STORE NEW DATA FOR THIS REFLECTION.
C
            IF (ISTOR) THEN
               IF (KHKL.NE.IHKL) GO TO 150
C
C EXTRA CHECKS IF DLAM > 0.0
C
               IF (DLAM.GT.0.0) THEN
                  IF (ABS(ALBDA-ALAM).GT.DLAM) GO TO 150
                  IF (NABS.GT.0) THEN
                     DO 100 III = 1 , NABS
                        IF (ABSL(III).GE.ALAM .AND. ABSL(III).LE.ALBDA)
     +                      GO TO 150
 100                 CONTINUE
                  END IF
               END IF
C
C TEST FOR HKL CHANGE; THEN IF SAME REFLECTION
C TEST FOR NEW BATCH OR REPEATED BATCH FOR THIS SAME REFLECTION
C
               DO 110 I1 = 1 , MSER
                  IF (N1.EQ.ISER(I1)) GO TO 140
 110           CONTINUE
               GO TO 130
            END IF
         END IF
C
C INITIALISE A NEW HKL
C
 120     TH = IA(JJ+2)/10000.0
         ALAM = IA(JJ+1)/10000.0
         TH = SIN(TH)/ALAM
         TH = TH*TH
         TH100 = TH*100.
         IHKL = KHKL
         MSER = 0
         ISTOR = .TRUE.
         INEG = 0
         FMODA = 0.
         FFMODA = 0.
         SDSUM = 0.
         NOSUM = 0
         F = 0.
         W = 0.
C
C INITIALISE A NEW BATCH FOR THE SAME HKL
C
C NOTE THAT TERM INVOLVING TH IS HERE SCALED UP BY 50.0
C
 130     MSER = MSER + 1
         ISER(MSER) = N1
         GG(MSER) = EXP(TH100*BTX(N1))/SCX(N1)
C
C PARTIAL DERIVATIVE  DGG(MSER)/DSCALE(N1)
C
         EX(MSER) = -GG(MSER)/SCX(N1)
         FF(MSER) = 0.
         WSD(MSER) = 0.
         NO(MSER) = 0
         I1 = MSER
C
C STORE AND SUM MEASUREMENTS OF REFLECTION
C SUM FOR AVERAGE OF ALL REFLECTIONS FROM SAME BATCH
C
 140     WJ = 0.0001
         IF (IWT.EQ.0) WJ = 1./(SDIN*SDIN)
         FF(I1) = FF(I1) + RMODF*WJ
         WSD(I1) = WSD(I1) + WJ
         NO(I1) = NO(I1) + 1
C
C SUMS FOR MEAN I(H) OVER ALL BATCHES
C
         G = GG(I1)
         F = F + RMODF*WJ*G
         W = W + WJ*G*G
         RMODF = RMODF/G
         SDIN = SDIN/G
         FMODA = FMODA + RMODF
         FFMODA = FFMODA + RMODF*RMODF
         SDSUM = SDSUM + SDIN*SDIN
         NOSUM = NOSUM + 1
         KOUNTA = KOUNTA + 1
         FBAR = FBAR + RMODF
         IF (RMODF.LE.0) THEN
C
C FLAG NEGATIVE INTENSITY - TO EXCLUDE THIS HKL FROM REFINEMENT.
C
            INEG = 1
         END IF
         GO TO 90
C
C NEW REFLECTION FOUND IN ARRAY 'IA'
C SO DEAL WITH PREVIOUS REFLECTION FOR WHICH ALL DATA IS NOW COMPLETE
C
 150     ISTOR = .FALSE.
C
C EXCLUDE ALL OBSERVATIONS FOR ANY HKL HAVING A NEGATIVE OBSERVATION
C OR ANY HKL FOR WHICH ALL OBSERVATIONS COME FROM A SINGLE BATCH.
C
         IF (INEG.NE.0) GO TO 120
         IF (MSER.LE.1) GO TO 120
 160     FMODQ = FMODA/FLOAT(NOSUM)
C
C WEIGHTED MEAN INTENSITY AND VARIANCE FOR ALL OBSERVNS. OF THIS HKL
C
         W = 1.0/W
         F = F*W
C
C WEIGHTED UNSCALED MEAN INTENSITY AND WEIGHT FOR EACH BATCH SEPARATELY
C
         DO 170 I = 1 , MSER
            FF(I) = FF(I)/WSD(I)
 170     CONTINUE
         IF (INEG.EQ.0) THEN
            IF (MSER.NE.1) THEN
C
C OUTER LOOP ROUND BATCHES
C
               DH = TH100
               DO 190 I = 1 , MSER
                  L = ISER(I)
                  FL = FF(I)
                  GL = GG(I)
                  WL = WSD(I)
                  DI = FL - GL*F
                  DHOGLP = DH*SCX(L)
C
C SUMS FOR RESIDUAL ANALYSIS
C
                  WD = WL*DI*DI
                  INTER(L) = INTER(L) + 1
                  LP = L + NBTOT
                  WD2 = WD2 + WD
                  ND2 = ND2 + 1
                  WD = WL*W*EX(I)
                  S1 = WD*(FL-GL*F)
                  S2 = WD*FL*F
C
C SUM RHS ELEMENTS
C
                  SR = WL*F*DI*EX(I)
                  COEFC(L) = COEFC(L) + SR
                  IF (BREF) COEFC(LP) = COEFC(LP) - SR*DHOGLP
C
C EXTRA COMPONENT OF MATRIX DIAGONAL
C
                  S3 = F*EX(I)
                  SR = WL*S3*S3
                  COEFS(L,L) = COEFS(L,L) + SR
                  IF (BREF) THEN
                     SR = SR*DHOGLP
                     COEFS(LP,LP) = COEFS(LP,LP) + SR*(DHOGLP)
C
C EXTRA COMPONENT FOR CROSS TERMS
C
                     COEFS(L,LP) = COEFS(L,LP) - SR
                  END IF
C
C INNER LOOP FOR MATRIX ELEMENTS
C
                  DO 180 J = I , MSER
                     M = ISER(J)
                     DHOGMP = DH*SCX(M)
                     SR = WSD(J)*EX(J)*(S1*FF(J)-S2*GG(J))
                     MP = M + NBTOT
                     COEFS(M,L) = COEFS(M,L) + SR
                     IOVLAP(L,M) = IOVLAP(L,M) + 1
                     IF (L.NE.M) IOVLAP(M,L) = IOVLAP(L,M)
                     IF (BREF) THEN
                        IF (L.NE.M) COEFS(M,LP) = COEFS(M,LP)
     +                      - SR*DHOGLP
                        SR = SR*DHOGMP
                        COEFS(MP,L) = COEFS(MP,L) - SR
                        COEFS(MP,LP) = COEFS(MP,LP) + SR*DHOGLP
                     END IF
 180              CONTINUE
C
 190           CONTINUE
            END IF
         END IF
C
C
C CHECK END OF REFLECTION DATA
C
         IF (.NOT.(IEND)) GO TO 120
C
C END OF ALL REFLECTIONS, NOW SOLVE EQUATIONS
C
 200     IF (DLAM.GT.0.0 .AND. ICYC.EQ.1) THEN
            WRITE (6,6050) (NQ,NQ=1,NBTOT)
            DO 210 NQ = 1 , NBTOT
               WRITE (6,6060) NQ , (IOVLAP(NQ,KQ),KQ=1,NBTOT)
 210        CONTINUE
         END IF
C
C BEGIN FOX & HOLMES CALCULATIONS USING ARRAYS SETUP EARLIER IN CYCLE.
C
C CHECK FOR ANY NULL ROWS OR COLUMNS IN MATRIX
C
         IF (ICYC.LE.1) THEN
            DO 220 L = 1 , NBTOT
               IF (INTER(L).EQ.0) GO TO 370
 220        CONTINUE
         END IF
C
C CALCULATE SCALE FACTOR TO BRING MEAN DIAGONAL OF MATRIX TO 1.0
C
         S2 = 0.0
         DO 230 I = 1 , NBLOCK
            S2 = S2 + COEFS(I,I)
 230     CONTINUE
         S2 = FLOAT(NBLOCK)/S2
C
C SCALE AND SYMMETRISE MATRIX.
C
         COEFS(1,1) = COEFS(1,1)*S2
         DO 250 I = 2 , NBLOCK
            COEFS(I,I) = COEFS(I,I)*S2
            K = I - 1
            DO 240 J = 1 , K
               COEFS(J,I) = (COEFS(J,I)+COEFS(I,J))*S2
               COEFS(I,J) = COEFS(J,I)
 240        CONTINUE
 250     CONTINUE
C
C CALCULATE MATRIX OF EIGENVECTORS VEC, AND EIGENVALUES VAL
C
         CALL EA06C(COEFS,VAL,VEC,NPARMS,NBLOCK,NBLOCK,WORK)
C FIND TWO SMALLEST EIGENVALUES
C
         VMIN1 = 1000.
         VMIN2 = 1000.
         I1 = 1
         DO 260 I = 1 , NPARMS
            VALI = VAL(I)
            IF (VALI.LE.VMIN1) THEN
               I2 = I1
               VMIN2 = VMIN1
               VMIN1 = VALI
               I1 = I
            ELSE IF (VALI.LE.VMIN2) THEN
               VMIN2 = VALI
               I2 = I
            END IF
 260     CONTINUE
C
C SET EIGENVALUES TO ZERO
C
         VAL(I1) = 0.
         IF (.NOT.BREF) THEN
            J = 1
         ELSE
            VAL(I2) = 0.
C      WRITE(6,3003) J,I1,I2
            J = 2
         END IF
C      WRITE(6,3003) J,I1
C
C MULTIPLY RHS (COEFC) BY VEC TRANSPOSE AND DIVIDE BY VAL
C (I.E. MULTIPLY BY INVERSE DIAGONAL MATRIX)
C TO GET TRANSFORMED SHIFT VECTOR
C
         DO 280 I = 1 , NPARMS
            DP(I) = 0.0
            IF (VAL(I).NE.0.0) THEN
               DO 270 J = 1 , NPARMS
                  DP(I) = DP(I) + VEC(J,I)*COEFC(J)
 270           CONTINUE
               DP(I) = (DP(I)/VAL(I))*S2
            END IF
 280     CONTINUE
C
C MEAN RESIDUAL = (WD**2)/(M-N)    USED AS SCALE FACTOR FOR VARIANCES
C
         WD2 = WD2/FLOAT(ND2-NPARMS)
         S = S2*WD2
C
C TRANSFORM SHIFT VECTOR BACK TO TRUE FRAME,
C AND RE-USE ARRAY COEFS FOR VARIANCE/COVARIANCE MATRIX
C VARIANCE/COVARIANCE MATRIX = VEC*(D**-1)*VEC TRANSPOSE
C
C APART FROM ITS DIAGONAL ELEMENTS, THE REMAINDER OF SYMMETRICAL
C MATRIX COEFS IS ONLY NEEDED FOR PRINTING OUT
C ONLY HALF OF IT IS REALLY NEEDED EVEN THEN
C
         DO 310 I = 1 , NPARMS
            SJ = 0.0
            DO 300 J = 1 , NPARMS
               SJ = SJ + VEC(I,J)*DP(J)
               SK = 0.0
               IF (I.EQ.J) THEN
                  DO 290 K = 1 , NPARMS
                     IF (VAL(K).NE.0.0) THEN
                        SK = SK + VEC(I,K)*VEC(J,K)/VAL(K)
                     END IF
 290              CONTINUE
               END IF
               COEFS(J,I) = SK*S
 300        CONTINUE
            COEFC(I) = SJ
 310     CONTINUE
         DO 320 I = 1 , NPARMS
            DP(I) = SQRT(COEFS(I,I))
 320     CONTINUE
         II = JUNITY
         J = JUNITY
         S1 = SCX(II)/(SCX(II)+COEFC(J))
         IF (BREF) THEN
            S2 = COEFC(NBTOT+J)
         END IF
C
C APPLY SHIFTS & RESET SCX,BTX,SDK,SDB. STORE RATIO SHIFT/SD IN VAL
C
         DO 330 I = 1 , NBTOT
            N1 = I
            N2 = NBTOT + I
            S = (SCX(N1)+COEFC(I))*S1
            COEFC(I) = S - SCX(N1)
            SCX(N1) = S
            SDK(N1) = DP(I)
            S = ABS(COEFC(I))
C           SHIFTK = AMAX1(S,SHIFTK)
            VAL(I) = S/DP(I)
C           DELK = AMAX1(VAL(I),DELK)
            IF (BREF) THEN
               S3 = COEFC(N2) - S2
               COEFC(N2) = 50.0*S3
               BTX(N1) = 50.0*(BTX(N1)+S3)
               SDB(N1) = 50.0*DP(N2)
               S3 = ABS(S3)
               BSHIFT = AMAX1(50.0*S3,BSHIFT)
               VAL(N2) = S3/DP(N2)
               DELB = AMAX1(VAL(N2),DELB)
            ELSE
               BTX(N1) = 50.0*BTX(N1)
               SDB(N1) = 0.0
            END IF
 330     CONTINUE
C
C WRITE OUT K, B, SHIFTS AND ERRORS FOR REFINED BATCHS
C
         IF (BREF) THEN
            WRITE (6,6010) NITR , ICYC
         ELSE
            WRITE (6,6030) NITR , ICYC
         END IF
         S1 = 0.0
         S2 = 0.0
         S3 = 0.0
         S4 = 0.0
         S5 = 0.0
         S6 = 0.0
         S7 = 0.0
         S8 = 0.0
         M = (NBTOT+1)/2
         N = NBTOT/2
         DO 340 I = 1 , M
            JP = I
            KP = JP + N
            J = I + NBTOT
            II = I + N
            JJ = J + N
            S1 = COEFC(I)
            S3 = VAL(I)
            S5 = COEFC(II)
            S7 = VAL(II)
            IF (BREF) THEN
               S2 = COEFC(J)
               S4 = VAL(J)
               S6 = COEFC(JJ)
               S8 = VAL(JJ)
            END IF
            IF (BREF) THEN
               WRITE (6,6020) JP , SCX(JP)*SCALE(JP) , BTX(JP)
     +                        + BTEMP(JP) , S1 , S2 , SDK(JP) , SDB(JP)
     +                        , INTER(JP) , KP , SCX(KP)*SCALE(KP) , 
     +                        BTX(KP) + BTEMP(KP) , S5 , S6 , SDK(KP) , 
     +                        SDB(KP) , INTER(KP)
            ELSE
               WRITE (6,6040) JP , SCX(JP)*SCALE(JP) , S1 , SDK(JP) , 
     +                        INTER(JP) , KP , SCX(KP)*SCALE(KP) , S5 , 
     +                        SDK(KP) , INTER(KP)
            END IF
 340     CONTINUE
C
C END OF CYCLES LOOP
C ------------------
C
 350  CONTINUE
      DO 360 I = 1 , NBTOT
         SCALE(I) = SCALE(I)*SCX(I)
         BTEMP(I) = BTEMP(I) + BTX(I)
 360  CONTINUE
      RETURN
C
C ERROR CONDITIONS
C ----------------
C
 370  WRITE (6,6070)
      STOP
C
C FORMAT STATEMENTS
C -----------------
C
 6010 FORMAT (///' ITERATION NUMBER',I4,'  REFINEMENT CYCLE',I4,//,
     +        '  BATCH  KSCALE   BTEMP  KSHIFT  BSHIFT  ESDK  ESDB',
     +        '  NOBS   ',
     +        '  BATCH  KSCALE   BTEMP  KSHIFT  BSHIFT  ESDK  ESDB',
     +        '  NOBS')
 6020 FORMAT (I7,2(F8.4,F8.3),F6.3,F6.2,I6,3X,I7,2(F8.4,F8.3),F6.3,F6.2,
     +        I6)
 6030 FORMAT (///' ITERATION NUMBER',I4,'  REFINEMENT CYCLE',I4,//,
     +        '  BATCH  KSCALE  KSHIFT  ESDK  NOBS   ',
     +        '  BATCH  KSCALE  KSHIFT  ESDK  NOBS  ')
 6040 FORMAT (I7,2(F8.4),F6.3,I6,3X,I7,2(F8.4),F6.3,I6)
 6050 FORMAT (//,' NUMBERS OF OVERLAPPING REFLECTIONS',//,(10X,20I6))
 6060 FORMAT (/,' ',I5,4X,20I6,/,(10X,20I6))
C
 6070 FORMAT (//,' **MISSING DATA FOR SOME BATCHES**')
      END
**==rfacts.f  
C
C SUBROUTINE 'RFACTS'
C ===================
C
C CALCULATE MERGING R-FACTORS FOR THE PROGRAM 'LAUENORM'
C
      SUBROUTINE RFACTS(IA,NREF,NBTOT,NPACK,IFSCAL,IR2,IR3)
C
C PARAMETERS
C
C          IA (I)   ARRAY HOLDING REFLECTION DATA (SEE MAIN PROGRAM
C                   FOR DETAILS OF THE FORMAT)
C        NREF (I)   NUMBER OF LAUE REFLECTIONS PRESENT
C       NBTOT (I)   NUMBER OF LAMBDA BINS INTO WHICH THE DATA ARE SPLIT
C       NPACK (I)   NUMBER OF FILM PACKS
C      IFSCAL (I)   LOGICAL FLAG =.TRUE.   APPLY BIN SCALING FACTORS
C                                =.FALSE.  DO NOT APPLY SCALE FACTORS
C                   DIVIDED
C         IR2 (I)   FLAG =1 OUTPUT 2'ND TABLE OF R-FACTORS
C                        =0 DO NOT
C         IR3 (I)   FLAG =1 OUTPUT 3'RD TABLE OF R-FACTORS
C                        =0 DO NOT
C
C PARAMETER STATEMENTS
C
C       MXB = MAXIMUM NUMBER OF LAMBDA BINS ALLOWED
C       MXR = MAXIMUM NUMBER OF REFLECTIONS
C
      PARAMETER(MXB=100,MXR=300000)
C
C SPECIFICATION STATEMENTS
C ------------------------
C
      LOGICAL IFSCAL,ERR
      DIMENSION IA(6*NREF)
      DIMENSION RFAC1(MXB,MXB),RFAC2(MXB,MXB),RFAC3(MXB,MXB)
      DIMENSION RSUM1(MXB,MXB),RSUM2(MXB,MXB),RSUM3(MXB,MXB)
      DIMENSION NREF1(MXB,MXB),NREF2(MXB,MXB),NREF3(MXB,MXB)
      DIMENSION RFACP(MXB,MXB),RSUMP(MXB,MXB),NREFP(MXB,MXB)
      LOGICAL RF(3),END
      DIMENSION SIGWI(3),SIGW(3),RN(3),RD(3),NR(3),AVE(3),DEL(3)
      DIMENSION SIGI(3),NI(3),RRN(3),RRD(3),AV(3),DL(3)
      EQUIVALENCE (IAA,AAA)
      DATA NTWO16/65536/
C
C INITIALISATIONS
C ---------------
C
      DO 30 I = 1 , NBTOT
         DO 20 J = 1 , NBTOT
            RFAC1(I,J) = 0.0
            RSUM1(I,J) = 0.0
            NREF1(I,J) = 0
            RFAC2(I,J) = 0.0
            RSUM2(I,J) = 0.0
            NREF2(I,J) = 0
            RFAC3(I,J) = 0.0
            RSUM3(I,J) = 0.0
            NREF3(I,J) = 0
 20      CONTINUE
 30   CONTINUE
      DO 40 I = 1 , 3
         RN(I) = 0.0
         RRN(I) = 0.0
         RD(I) = 0.0
         RRD(I) = 0.0
         NR(I) = 0
 40   CONTINUE
      DO 60 I = 1 , NPACK
         DO 50 J = 1 , NPACK
            RFACP(I,J) = 0.0
            RSUMP(I,J) = 0.0
            NREFP(I,J) = 0
 50      CONTINUE
 60   CONTINUE
      END = .FALSE.
      M1 = 1
      ILAST = IA(1)
      K = 0
C
C STEP THROUGH REFLECTIONS
C ------------------------
C
 70   K = K + 1
      IF (K.GT.NREF) THEN
         END = .TRUE.
         GO TO 80
      END IF
      JJ = 6*K - 5
      IF (IA(JJ).EQ.ILAST) GO TO 70
 80   M2 = K - 1
      IF (M2.GT.M1) THEN
         DO 120 IM = M1 , M2
            IP = 6*IM - 5
            ALM = IA(IP+1)/10000.0
            THETA = IA(IP+2)/10000.0
            IBIN = IBPACK(1,IA(IP+5),ERR)
            IPACK = IBPACK(2,IA(IP+5),ERR)
            SC = 1.0
            IF (IFSCAL) SC = SCFAC(IBIN,IPACK,ALM,THETA)
            IAA = IA(IP+3)
            AI = AAA*SC
            KSYM = IA(IP+4)/NTWO16
            SI = (IA(IP+4)-NTWO16*KSYM)*SC
            DO 90 I = 1 , 3
               SIGWI(I) = 0.0
               SIGW(I) = 0.0
               SIGI(I) = 0.0
               NI(I) = 0
               RF(I) = .FALSE.
 90         CONTINUE
C
C FIND THE MEAN INTENSITIES FOR THE REFLECTION REQUIRED FOR CALCULATING
C THE R-FACTORS. THE MEANS ARE (1) THE MEAN OF ALL MEASUREMENTS,
C (2) THE MEAN OF ALL MEASUREMENTS OF THE SAME SIGN AS THE CURRENT
C REFLECTION (I.E. I+ OR I-) AND (3) THE MEAN OF MEASUREMENTS WITH THE
C SAME SIGN AND MEASURED WITHINN 0.1 ANGSTROM OF THE SAME WAVELENGTH.
C
            DO 100 JM = M1 , M2
               IF (IM.NE.JM) RF(1) = .TRUE.
               JP = 6*JM - 5
               AALM = IA(JP+1)/10000.0
               AATH = IA(JP+2)/10000.0
               JBIN = IBPACK(1,IA(JP+5),ERR)
               JPACK = IBPACK(2,IA(JP+5),ERR)
               SSC = 1.0
               IF (IFSCAL) SSC = SCFAC(JBIN,JPACK,AALM,AATH)
               IAA = IA(JP+3)
               AAI = AAA*SSC
               KKSYM = IA(JP+4)/NTWO16
               SSI = (IA(JP+4)-NTWO16*KKSYM)*SSC
               WWI = 1/(SSI*SSI)
               SIGWI(1) = SIGWI(1) + WWI*AAI
               SIGW(1) = SIGW(1) + WWI
               SIGI(1) = SIGI(1) + AAI
               NI(1) = NI(1) + 1
               IF (IM.LT.JM) THEN
                  RFAC1(IBIN,JBIN) = RFAC1(IBIN,JBIN) + ABS(AI-AAI)
                  RSUM1(IBIN,JBIN) = RSUM1(IBIN,JBIN) + AI + AAI
                  NREF1(IBIN,JBIN) = NREF1(IBIN,JBIN) + 1
                  RFAC1(JBIN,IBIN) = RFAC1(IBIN,JBIN)
                  RSUM1(JBIN,IBIN) = RSUM1(IBIN,JBIN)
                  NREF1(JBIN,IBIN) = NREF1(IBIN,JBIN)
                  RFACP(IPACK,JPACK) = RFACP(IPACK,JPACK) + ABS(AI-AAI)
                  RSUMP(IPACK,JPACK) = RSUMP(IPACK,JPACK) + AI + AAI
                  NREFP(IPACK,JPACK) = NREFP(IPACK,JPACK) + 1
                  RFACP(JPACK,IPACK) = RFACP(IPACK,JPACK)
                  RSUMP(JPACK,IPACK) = RSUMP(IPACK,JPACK)
                  NREFP(JPACK,IPACK) = NREFP(IPACK,JPACK)
               END IF
               IF (MOD(KSYM,2).EQ.MOD(KKSYM,2)) THEN
C
C SAME SIGN
C
                  SIGWI(2) = SIGWI(2) + WWI*AAI
                  SIGW(2) = SIGW(2) + WWI
                  SIGI(2) = SIGI(2) + AAI
                  NI(2) = NI(2) + 1
                  IF (IM.NE.JM) RF(2) = .TRUE.
                  IF (IM.LT.JM) THEN
                     RFAC2(IBIN,JBIN) = RFAC2(IBIN,JBIN) + ABS(AI-AAI)
                     RSUM2(IBIN,JBIN) = RSUM2(IBIN,JBIN) + AI + AAI
                     NREF2(IBIN,JBIN) = NREF2(IBIN,JBIN) + 1
                     RFAC2(JBIN,IBIN) = RFAC2(IBIN,JBIN)
                     RSUM2(JBIN,IBIN) = RSUM2(IBIN,JBIN)
                     NREF2(JBIN,IBIN) = NREF2(IBIN,JBIN)
                  END IF
                  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
                     SIGI(3) = SIGI(3) + AAI
                     NI(3) = NI(3) + 1
                     IF (IM.NE.JM) RF(3) = .TRUE.
                     IF (IM.LT.JM) THEN
                        RFAC3(IBIN,JBIN) = RFAC3(IBIN,JBIN)
     +                     + ABS(AI-AAI)
                        RSUM3(IBIN,JBIN) = RSUM3(IBIN,JBIN) + AI + AAI
                        NREF3(IBIN,JBIN) = NREF3(IBIN,JBIN) + 1
                        RFAC3(JBIN,IBIN) = RFAC3(IBIN,JBIN)
                        RSUM3(JBIN,IBIN) = RSUM3(IBIN,JBIN)
                        NREF3(JBIN,IBIN) = NREF3(IBIN,JBIN)
                     END IF
                  END IF
               END IF
 100        CONTINUE
C
C ACCUMULATE VALUES FOR R-FACTORS
C
            DO 110 I = 1 , 3
               AVE(I) = SIGWI(I)/SIGW(I)
               DEL(I) = ABS(AI-AVE(I))
               AV(I) = SIGI(I)/NI(I)
               DL(I) = ABS(AI-AV(I))
               IF (RF(I)) THEN
                  RN(I) = RN(I) + DEL(I)
                  RD(I) = RD(I) + AVE(I)
                  RRN(I) = RRN(I) + DL(I)
                  RRD(I) = RRD(I) + AV(I)
                  NR(I) = NR(I) + 1
               END IF
 110        CONTINUE
C
 120     CONTINUE
      END IF
      IF (END) THEN
C
C END OF REFLECTIONS REACHED, OUTPUT TABLES OF R-FACTORS
C ------------------------------------------------------
C
         DO 140 I = 1 , NBTOT
            DO 130 J = 1 , NBTOT
               IF (NREF1(I,J).GT.0) RFAC1(I,J) = RFAC1(I,J)/RSUM1(I,J)
               IF (NREF2(I,J).GT.0) RFAC2(I,J) = RFAC2(I,J)/RSUM2(I,J)
               IF (NREF3(I,J).GT.0) RFAC3(I,J) = RFAC3(I,J)/RSUM3(I,J)
 130        CONTINUE
 140     CONTINUE
         DO 160 I = 1 , NPACK
            DO 150 J = 1 , NPACK
               IF (NREFP(I,J).GT.0) RFACP(I,J) = RFACP(I,J)/RSUMP(I,J)
 150        CONTINUE
 160     CONTINUE
         DO 170 I = 1 , 3
            IF (NR(I).GT.0) THEN
               RN(I) = RN(I)/RD(I)
               RRN(I) = RRN(I)/RRD(I)
            END IF
 170     CONTINUE
         WRITE (6,6010)
         WRITE (6,6020)
         CALL ROUT(RFAC1,NREF1,NBTOT)
         IF (IR2.NE.0) THEN
            WRITE (6,6030)
            CALL ROUT(RFAC2,NREF2,NBTOT)
         END IF
         IF (IR3.NE.0) THEN
            WRITE (6,6040)
            CALL ROUT(RFAC3,NREF3,NBTOT)
         END IF
         WRITE (6,6050)
         WRITE (6,6060) RRN(1)
         WRITE (6,6090) NR(1)
         WRITE (6,6070) RRN(2)
         WRITE (6,6090) NR(2)
         WRITE (6,6080) RRN(3)
         WRITE (6,6090) NR(3)
         WRITE (6,6100) RN
         WRITE (6,6110)
         WRITE (6,6020)
         CALL ROUT(RFACP,NREFP,NPACK)
         RETURN
      ELSE
         M1 = K
         ILAST = IA(JJ)
         GO TO 70
      END IF
C
C FORMAT STATEMENTS
C -----------------
C
 6010 FORMAT (//,' LAUE INTER-BIN MERGING R-FACTORS = ',
     +        'SUM(ABS(I(I)-I(J)))/SUM((I(I)+I(J))')
 6020 FORMAT (/,' ALL PAIRS OF MEASUREMENTS')
 6030 FORMAT (/,' PAIRS OF MEASUREMENTS OF THE SAME SIGN (I+ OR I-)')
 6040 FORMAT (/,' PAIRS OF MEASUREMENTS OF THE SAME SIGN AND WITH',
     +        ' WAVELENGTHS WITHIN 0.1 ANGSTROM')
 6050 FORMAT (//,' LAUE MERGING R-FACTORS = SUM(ABS(I-IMEAN)/SUM(I)',//,
     +        '                          (SUM OVER ALL MEASUREMENTS)')
 6060 FORMAT (//,' RFACT1 =',F6.3,
     +        '     IMEAN FROM ALL MEASUREMENTS FOR',
     +        ' THE REFLECTION (I+ AND I-)')
 6070 FORMAT (//,' RFACT2 =',F6.3,'     IMEAN FROM ALL MEASUREMENTS OF',
     +        ' THE SAME SIGN (I+ OR I-)')
 6080 FORMAT (//,' RFACT3 =',F6.3,'     IMEAN FROM ALL MEASUREMENTS OF',
     +        ' THE SAME SIGN (I+ OR I-) AND WITH LAMBDA WITHIN 0.1A')
 6090 FORMAT (/,' N =',I6)
 6100 FORMAT (//,' (WITH WEIGHTED IMEAN RFACT1 =',F6.3,'   RFACT2 =',
     +        F6.3,'   RFACT3 =',F6.3,')')
 6110 FORMAT (//,' LAUE INTER-FILMPACK MERGING R-FACTORS = ',
     +        'SUM(ABS(I(I)-I(J)))/SUM((I(I)+I(J))')
      END
**==rout.f  
C
C SUBROUTINE 'ROUT'
C =================
C
C OUTPUT A TABLE OF INTER BIN/FILMPACK R-FACTORS
C
      SUBROUTINE ROUT(RFAC,NREF,NBTOT)
C
C PARAMETERS
C
C        RFAC (I)   ARRAY HOLDING R-FACTORS
C        NREF (I)   ARRAY HOLDING NUMBERS OF REFLECTIONS
C       NBTOT (I)   NUMBER OF BINS/FILMPACKS CALCULATED
C
C PARAMETER STATEMENTS
C
C               MXB = MAXIMUM NUMBER OF LAMBDA BINS/FILMPACKS ALLOWED
C
      PARAMETER(MXB=100)
C
C SPECIFICATION STATEMENTS
C ------------------------
C
      DIMENSION RFAC(MXB,MXB),NREF(MXB,MXB)
C
C OUTPUT TABLE
C ------------
C
      WRITE (6,6010) (N,N=1,NBTOT)
      NLINE = (NBTOT-1)/16 + 1
      DO 30 I = 1 , NBTOT
         K1 = -15
         DO 20 N = 1 , NLINE
            K1 = K1 + 16
            K2 = K1 + 15
            IF (K2.GT.NBTOT) K2 = NBTOT
            IF (N.EQ.1) THEN
               WRITE (6,6020) I , (RFAC(I,K),K=K1,K2)
            ELSE
               WRITE (6,6030) (RFAC(I,K),K=K1,K2)
            END IF
            WRITE (6,6040) (NREF(I,K),K=K1,K2)
 20      CONTINUE
 30   CONTINUE
      RETURN
C
C FORMAT STATEMENTS
C -----------------
C
 6010 FORMAT (/,10X,20I6)
 6020 FORMAT (/,' ',I5,4X,20F6.3)
 6030 FORMAT (10X,20F6.3)
 6040 FORMAT (10X,20I6)
      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(MXL=8,MXB=100)
C
C SPECIFICATION STATEMENTS
C ------------------------
C
      LOGICAL CFLAG,CBFLAG,CLFLAG
      CHARACTER*1 SEL
      DOUBLE PRECISION COEFFS(11),RMS(11),CW(100),CR(100),CLAM(100)
      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),CFLAG,
     *              CBFLAG(MXB),CLFLAG(MXL)
      COMMON /BINS/NRANGE,NBTOT,NVAL(MXB),SCALE(MXB),BTEMP(MXB),
     *ALAM(MXB),CSCALE(MXB)
C
C PREPARE POINTS FOR CURVE FITTING, ELIMINATING BINS WITHOUT DATA
C ---------------------------------------------------------------
C
      CALL HEADXX(6,'CURVE FITTING RESULTS$','$','0',2)
      DO 60 NN = 1 , NRANGE
         WRITE (6,6010) NN
         IF (.NOT.CLFLAG(NN)) THEN
            WRITE (6,6090)
            GO TO 60
         END IF
         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 = SCALE(NB1+2)
            AL1 = ALAM(NB1+1)
            SC1 = SCALE(NB1+1)
            AL2 = ALAM(NB1)
            SC2 = SCALE(NB1)
            NPOINT = 1
            CW(NPOINT) = 1.0
            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) = 1.0
            CR(NPOINT) = SCALE(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 = SCALE(NB2-2)
            AL1 = ALAM(NB2-1)
            SC1 = SCALE(NB2-1)
            AL2 = ALAM(NB2)
            SC2 = SCALE(NB2)
            NPOINT = NPOINT + 1
            CW(NPOINT) = 1.0
            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)
         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,6100)
      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,//,
     +        (' ',6D20.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,F11.3)
 6090 FORMAT (/,' NO CURVE FITTING DONE FOR THIS RANGE')
 6100 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,MXB=100)
C
C SPECIFICATION STATEMENTS
C ------------------------
C
      LOGICAL CFLAG,CBFLAG,CLFLAG
      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),CFLAG,
     *              CBFLAG(MXB),CLFLAG(MXL)
      COMMON /BINS/NRANGE,NBTOT,NVAL(MXB),SCALE(MXB),BTEMP(MXB),
     *ALAM(MXB),CSCALE(MXB)
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
      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
**==scfac.f  
C
C FUNCTION 'SCFAC'
C ================
C
C RETURN A REFLECTION SCALE FACTOR
C
      FUNCTION SCFAC(IB,IPACK,ALBDA,THETA)
C
C PARAMETERS
C
C        IB (I)   LAMBDA BIN NUMBER
C     IPACK (I)   FILM PACK NUMBER
C     ALBDA (I)   LAMBDA VALUE
C     THETA (I)   THETA VALUE
C
      PARAMETER(MXL=8,MXB=100)
C
C SPECIFICATION STATEMENTS
C ------------------------
C
      LOGICAL CFLAG,CBFLAG,CLFLAG
      COMMON /BINS/NRANGE,NBTOT,NVAL(MXB),SCALE(MXB),BTEMP(MXB),
     *ALAM(MXB),CSCALE(MXB)
      COMMON /PACKS/NPACK,NPVAL(MXB),PACKSC(MXB),PACKTF(MXB),
     *CPSCAL(MXB),CPTFAC(MXB)
      COMMON /CURVE/CMIN(MXL),CMAX(MXL),NCOEFF(MXL),C(MXL,11),CFLAG,
     *              CBFLAG(MXB),CLFLAG(MXL)
C
C FIND THE SCALE FACTOR
C ---------------------
C
      IF (CFLAG .AND. CBFLAG(IB)) THEN
         SCFAC = CURVAL(ALBDA)
      ELSE
         SCFAC = CSCALE(IB)
      END IF
      IF (CPTFAC(IPACK).NE.0.0) THEN
         STOL = SIN(THETA)/ALBDA
         STOL2 = STOL*STOL
         SCPAK = CPSCAL(IPACK)*EXP(-2.0*CPTFAC(IPACK)*STOL2)
      ELSE
         SCPAK = CPSCAL(IPACK)
      END IF
      SCFAC = SCFAC*SCPAK
      RETURN
      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
**==table1.f  
C
C SUBROUTINE 'TABLE1'
C ===================
C
C OUTPUT A TABLE AND GRAPH SHOWING THE RESULTS OF THE CURVE FITTING
C
      SUBROUTINE TABLE1
C
C PARAMETERS
C
C NONE
C
      PARAMETER(MXL=8,MXB=100)
C
C SPECIFICATION STATEMENTS
C ------------------------
C
      LOGICAL CFLAG,CBFLAG,CLFLAG
      CHARACTER*32 GRLIN,GR0,GR1
      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),CFLAG,
     *              CBFLAG(MXB),CLFLAG(MXL)
      COMMON /BINS/NRANGE,NBTOT,NVAL(MXB),SCALE(MXB),BTEMP(MXB),
     *ALAM(MXB),CSCALE(MXB)
      DATA GR0/'................................'/
      DATA GR1/'.                              .'/
C
C OUTPUT STATISTICS
C -----------------
C
      CALL HEADXX(6,'CURVE FITTING ANALYSIS$','$','1',2)
      WRITE (6,6010)
      RMAX = 1E-10
      RMIN = 1E+10
      DO 20 N = 1 , NBTOT
         IF (CBFLAG(N)) THEN
            CV = CURVAL(ALAM(N))
            IF (CV.GT.RMAX) RMAX = CV
            IF (SCALE(N).GT.RMAX) RMAX = SCALE(N)
            IF (CV.LT.RMIN) RMIN = CV
            IF (SCALE(N).LT.RMIN) RMIN = SCALE(N)
         END IF
 20   CONTINUE
      DO 30 NN = 1 , NRANGE
         IF (CLFLAG(NN)) THEN
            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))
         END IF
 30   CONTINUE
      DO 70 NN = 1 , NRANGE
         IF (CLFLAG(NN)) THEN
            WRITE (6,6070)
            WRITE (6,6050) GR0
            GRLIN = GR1
            CV = CURVAL(BMINL(NN))
            EXT1 = CV
            IP = 15
            IF (RMAX.GT.RMIN) IP = NINT(29*(CV-RMIN)/(RMAX-RMIN)) + 2
            GRLIN(IP:IP) = 'E'
            WRITE (6,6050) GRLIN
            DO 50 NA = 1 , NBIN(NN)
               N = NP(NN) + NA
               X = 0.0
               IF (ALAM(N).GE.CMIN(NN) .AND. ALAM(N).LE.CMAX(NN)) THEN
                  X = CURVAL(ALAM(N))
               END IF
               GRLIN = GR1
               IP1 = 15
               IP2 = 15
               IF (RMAX.GT.RMIN) THEN
                  IP1 = NINT(29*(SCALE(N)-RMIN)/(RMAX-RMIN)) + 2
                  CV = CURVAL(ALAM(N))
                  IP2 = NINT(29*(CV-RMIN)/(RMAX-RMIN)) + 2
               END IF
               IF (IP1.EQ.IP2) THEN
                  GRLIN(IP1:IP1) = '*'
               ELSE
                  GRLIN(IP1:IP1) = 'X'
                  GRLIN(IP2:IP2) = '+'
               END IF
               WRITE (6,6020) N , ALAM(N) , NVAL(N) , SCALE(N) , X , 
     +                        GRLIN
               IF (NA.NE.NBIN(NN)) WRITE (6,6050) 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,6050) GRLIN
            WRITE (6,6050) GR0
            WRITE (6,6080) EXT1 , EXT2
         ELSE
            DO 60 NA = 1 , NBIN(NN)
               N = NP(NN) + NA
               WRITE (6,6040) N , ALAM(N) , NVAL(N) , SCALE(N)
 60         CONTINUE
         END IF
 70   CONTINUE
      WRITE (6,6030)
      RETURN
C
C FORMAT STATEMENTS
C -----------------
C
 6010 FORMAT (//,'  BIN  AVERAGE  NO. OF ',' SCALE      CURVE FITTED',/,
     +        '  NO.   LAMBDA  REFLNS','  FACTOR     SCALE FACTOR   ')
 6020 FORMAT (' ',I4,1X,F7.3,I7,F10.3,F14.3,12X,A)
 6030 FORMAT (//,' KEY:   E = EXTRAPOLATED POINT ON GRAPH',/,
     +        '        X = POINT USED FOR CURVE FITTING',/,
     +        '        + = CALCULATED POINT FROM FITTED CURVE',/,
     +        '        * = POINTS OVERLAP (OR X AND + OVERPRINTED)',/)
 6040 FORMAT (/,' ',I4,1X,F7.3,I7,F10.3,'         NONE')
 6050 FORMAT (' ',55X,A)
 6060 FORMAT ('+',55X,A)
 6070 FORMAT (//)
 6080 FORMAT (/,' EXTRAPOLATED SCALE AT FIRST POINT =',F10.3,/,
     +        ' EXTRAPOLATED SCALE AT FINAL POINT =',F10.3)
      END
**==opfile.f  
C
C SUBROUTINE 'OPFILE'
C ===================
C
C WRITE THE OUTPUT MTZ FILE (MERGED OR UNMERGED)
C
      SUBROUTINE OPFILE(IOUT,IA,NREF,ICENT,SCALOP)
C
C PARAMETERS
C
C       IOUT (I)  OUTPUT FILE TYPE = 1 MERGED, NO ANOMALOUS
C                                  = 2 MERGED WITH ANOMALOUS
C                                  > 2 UNMERGED MULTIRECORD OUTPUT
C                                    WITH BATCH SERIAL NUMBER 'IOUT'
C                                  =-1 SHELX OUTPUT FILE (MERGED, NO
C                                      ANOMALOUS)
C                                  =-2 SHELX OUTPUT (UNMERGED DATA WITH
C                                      LAMBDA VALUES)
C         IA (I)  ARRAY HOLDING REFLECTION DATA (SEE MAIN PROGRAM FOR
C                 DETAILS)
C       NREF (I)  NUMBER OF REFLECTIONS IN 'IA'
C      ICENT (I)  7 WORD ARRAY HOLDING CENTRIC ZONES FLAGS
C     SCALOP (I)  SCALE FACTOR FOR OUTPUT F'S FOR MERGED DATA
C
C SPECIFICATION STATEMENTS
C ------------------------
C
      LOGICAL CENTRC,CSYMM,END,SWAP
      DIMENSION IA(6*NREF)
      DIMENSION IDATA(8),ADATA(8),ICENT(7),HKL(3)
      DIMENSION SIGWI(3),SIGW(3),KN(3),ISYM(100),IPOINT(100)
      EQUIVALENCE (IAA,AAA)
      DATA NTWO16/65536/
C
C INVERT POINT GROUP MATRICES
C----------------------------
C
      CALL SYMINV
C
C WRITE MERGED OUTPUT DATA FILE
C -----------------------------
C
      IF (IOUT.GT.2) THEN
C
C WRITE UNMERGED DATA OUTPUT FILE (MTZ)
C -------------------------------------
C
C
C INITIALISATIONS
C
         END = .FALSE.
         NLARGE = 0
         NROUT = 0
         M1 = 1
         ILAST = IA(1)
         ADATA(8) = 0.0
         ADATA(5) = IOUT
         K = 0
         GO TO 80
      ELSE IF (IOUT.EQ.-2) THEN
C
C WRITE UNMERGED DATA OUTPUT 'SHELX' TYPE FILE WITH LAMBDAS
C ---------------------------------------------------------
C
C
C INITIALISATIONS
C
         END = .FALSE.
         NROUT = 0
         M1 = 1
         ILAST = IA(1)
         K = 0
         GO TO 160
      ELSE
C
C INITIALISATIONS
C
         CSYMM = .TRUE.
         END = .FALSE.
         NROUT = 0
         M1 = 1
         ILAST = IA(1)
         K = 0
      END IF
C
C STEP THROUGH THE REFLECTIONS
C
 20   K = K + 1
      IF (K.GT.NREF) THEN
         END = .TRUE.
         GO TO 30
      END IF
      JJ = 6*K - 5
      IF (IA(JJ).EQ.ILAST) GO TO 20
 30   M2 = K - 1
      DO 40 I = 1 , 3
         SIGWI(I) = 0.0
         SIGW(I) = 0.0
         KN(I) = 0
 40   CONTINUE
      MM = 6*M1 - 5
      CALL GETHKL(IA(MM),IDATA)
      ADATA(1) = IDATA(1)
      ADATA(2) = IDATA(2)
      ADATA(3) = IDATA(3)
      ALM = IA(MM+1)/10000.0
      TH = IA(MM+2)/10000.0
      STH = SIN(TH)/ALM
      IF (IOUT.EQ.2) CSYMM = CENTRC(IDATA,ICENT)
      ADATA(6) = 0.0
      ADATA(7) = 0.0
      DO 50 M = M1 , M2
         JP = 6*M - 5
         IAA = IA(JP+3)
         AI = AAA
         KSYM = IA(JP+4)/NTWO16
         SIGI = IA(JP+4) - NTWO16*KSYM
         W = 1.0/(SIGI*SIGI)
         SIGWI(1) = SIGWI(1) + W*AI
         SIGW(1) = SIGW(1) + W
         KN(1) = KN(1) + 1
         IF (.NOT.CSYMM) THEN
            KS = 3 - MOD(KSYM,2)
            SIGWI(KS) = SIGWI(KS) + W*AI
            SIGW(KS) = SIGW(KS) + W
            KN(KS) = KN(KS) + 1
         END IF
 50   CONTINUE
      FF = SIGWI(1)/SIGW(1)
      IF (FF.LT.0.0) GO TO 70
      SS = SQRT(1.0/SIGW(1))
      F = SQRT(FF)
      SIGF = SQRT(FF+SS) - SQRT(FF)
      IF (IOUT.EQ.-1) THEN
         WRITE (2,6020) IDATA(1) , IDATA(2) , IDATA(3) , F , SIGF
         NROUT = NROUT + 1
         GO TO 70
      END IF
      ADATA(4) = SCALOP*F
      ADATA(5) = SCALOP*SIGF
      JF = KN(2)*KN(3)
      IF (JF.GT.0) THEN
         FF1 = SIGWI(2)/SIGW(2)
         IF (FF1.LT.0.0) GO TO 70
         SS1 = SQRT(1.0/SIGW(2))
         FF2 = SIGWI(3)/SIGW(3)
         IF (FF2.LT.0.0) GO TO 70
         SS2 = SQRT(1.0/SIGW(3))
         SDF1 = SQRT(FF1+SS1) - SQRT(FF1)
         SDF2 = SQRT(FF2+SS2) - SQRT(FF2)
         DF = SQRT(FF1) - SQRT(FF2)
         SDF = SQRT(SDF1*SDF1+SDF2*SDF2)
         ADATA(6) = SCALOP*DF
         ADATA(7) = SCALOP*SDF
      END IF
      CALL LWREFL (1,ADATA)
      NROUT = NROUT + 1
 70   IF (END) THEN
         IF (IOUT.GT.0) CALL LWCLOS (1,1)
         WRITE (6,6010) NROUT
         RETURN
      ELSE
         M1 = K
         ILAST = IA(JJ)
         GO TO 20
      END IF
C
C STEP THROUGH THE REFLECTIONS
C
 80   K = K + 1
      IF (K.GT.NREF) THEN
         END = .TRUE.
         GO TO 90
      END IF
      JJ = 6*K - 5
      IF (IA(JJ).EQ.ILAST) GO TO 80
 90   M2 = K - 1
      IF (M2-M1.GT.99) GO TO 230
      MM = 6*M1 - 5
      CALL GETHKL(IA(MM),IDATA)
      ADATA(1) = IDATA(1)
      ADATA(2) = IDATA(2)
      ADATA(3) = IDATA(3)
      NVAL = 0
      DO 100 M = M1 , M2
         NVAL = NVAL + 1
         IP = 6*M - 5
         KSYM = IA(IP+4)/NTWO16
         ISYM(NVAL) = KSYM
         IPOINT(NVAL) = NVAL
 100  CONTINUE
C
C SORT EACH REFLECTION ON M/ISYM (=KSYM)
C
      IF (NVAL.GE.2) THEN
         NMAX = NVAL - 1
         DO 120 N = 1 , NMAX
            LMAX = NVAL - N
            SWAP = .FALSE.
            DO 110 L = 1 , LMAX
               IF (ISYM(L).GT.ISYM(L+1)) THEN
                  SWAP = .TRUE.
                  II = ISYM(L)
                  ISYM(L) = ISYM(L+1)
                  ISYM(L+1) = II
                  II = IPOINT(L)
                  IPOINT(L) = IPOINT(L+1)
                  IPOINT(L+1) = II
               END IF
 110        CONTINUE
            IF (.NOT.SWAP) GO TO 130
 120     CONTINUE
      END IF
 130  DO 150 N = 1 , NVAL
         II = M1 - 1 + IPOINT(N)
         IP = 6*II - 5
         IAA = IA(IP+3)
         ADATA(6) = AAA
         KSYM = IA(IP+4)/NTWO16
         ADATA(4) = KSYM
         ADATA(7) = IA(IP+4) - NTWO16*KSYM
         CALL LWREFL (1,ADATA)
         NROUT = NROUT + 1
 150  CONTINUE
      IF (END) THEN
         IF (IOUT.GT.0) CALL LWCLOS (1,1)
         WRITE (6,6010) NROUT
         RETURN
      ELSE
         M1 = K
         ILAST = IA(JJ)
         GO TO 80
      END IF
C
C STEP THROUGH THE REFLECTIONS
C
 160  K = K + 1
      IF (K.GT.NREF) THEN
         END = .TRUE.
         GO TO 170
      END IF
      JJ = 6*K - 5
      IF (IA(JJ).EQ.ILAST) GO TO 160
 170  M2 = K - 1
      IF (M2-M1.LE.99) THEN
         MM = 6*M1 - 5
         CALL GETHKL(IA(MM),IDATA)
         ADATA(1) = IDATA(1)
         ADATA(2) = IDATA(2)
         ADATA(3) = IDATA(3)
         NVAL = 0
         DO 180 M = M1 , M2
            NVAL = NVAL + 1
            IP = 6*M - 5
            KSYM = IA(IP+4)/NTWO16
            ISYM(NVAL) = KSYM
            IPOINT(NVAL) = NVAL
 180     CONTINUE
C
C SORT EACH REFLECTION ON M/ISYM (=KSYM)
C
         IF (NVAL.GE.2) THEN
            NMAX = NVAL - 1
            DO 200 N = 1 , NMAX
               LMAX = NVAL - N
               SWAP = .FALSE.
               DO 190 L = 1 , LMAX
                  IF (ISYM(L).GT.ISYM(L+1)) THEN
                     SWAP = .TRUE.
                     II = ISYM(L)
                     ISYM(L) = ISYM(L+1)
                     ISYM(L+1) = II
                     II = IPOINT(L)
                     IPOINT(L) = IPOINT(L+1)
                     IPOINT(L+1) = II
                  END IF
 190           CONTINUE
               IF (.NOT.SWAP) GO TO 210
 200        CONTINUE
         END IF
 210     DO 220 N = 1 , NVAL
            II = M1 - 1 + IPOINT(N)
            IP = 6*II - 5
            IAA = IA(IP+3)
            FF = AAA
            IF (FF.LT.0.0) GO TO 220
            KSYM = IA(IP+4)/NTWO16
            SS = IA(IP+4) - NTWO16*KSYM
            F = SQRT(FF)
            SIGF = SQRT(FF+SS) - SQRT(FF)
            ALAM = IA(IP+1)/10000.0
            CALL DECSYM(IDATA,HKL,KSYM)
            WRITE (2,6020) NINT(HKL(1)), NINT(HKL(2)), NINT(HKL(3)),
     +                     F , SIGF , ALAM
            NROUT = NROUT + 1
 220     CONTINUE
         IF (END) THEN
            WRITE (6,6010) NROUT
            RETURN
         ELSE
            M1 = K
            ILAST = IA(JJ)
            GO TO 160
         END IF
      END IF
C
C ERROR CONDITIONS
C ================
C
 230  WRITE (6,6030)
      STOP
C
C FORMAT STATEMENTS
C -----------------
C
 6010 FORMAT (//,' NUMBER OF REFLECTIONS WRITTEN TO OUTPUT FILE  =',I7)
 6020 FORMAT (3I5,3F8.2)
 6030 FORMAT (/,' **MORE THAN 100 MEASUREMENTS FOUND FOR A REFLECTION**'
     +        )
      END
**==gethkl.f  
C
C SUBROUTINE 'GETHKL'
C ===================
C
C THIS SUBROUTINE UNPACKS THE INDICES FOR A LAUE REFLECTION
C
      SUBROUTINE GETHKL(IDX,IDXL)
C
C PARAMETERS
C
C      IDX (I)   WORD CONTAINING PACKED INDICES
C     IDXL (O)   3 WORD ARRAY CONTAINING H, K, L.
C
C SPECIFICATION STATEMENTS
C ------------------------
C
      DIMENSION IDXL(3)
      DATA NTWO16/65536/,NTWO8/256/
C
C UNPACK INDICES
C --------------
C
      I = IDX
      IDXL(1) = I/NTWO16
      I = I - NTWO16*IDXL(1)
      IDXL(1) = IDXL(1) - 128
      IDXL(2) = I/NTWO8
      IDXL(3) = I - NTWO8*IDXL(2) - 128
      IDXL(2) = IDXL(2) - 128
      RETURN
      END
**==getmul.f  
C
C SUBROUTINE 'GETMUL'
C ===================
C
C THIS SUBROUTINE UNPACKS THE MULTIPLICITY & HARMONIC LIMITS FOR
C A LAUE MULTIPLE REFLECTION
C
      SUBROUTINE GETMUL(IDX,MULT,MINHRM,MAXHRM)
C
C PARAMETERS
C
C      IDX (I)   WORD CONTAINING PACKED DATA 2**20*MULT+2**10*MINHRM+MAXHRM
C     MULT (O)   RETURNS MULTIPLICITY
C   MINHRM (O)   RETURNS MINIMUM HARMONIC
C   MAXHRM (O)   RETURNS MAXIMUM HARMONIC
C
C SPECIFICATION STATEMENTS
C ------------------------
C
      DATA NTWO20/1048576/,NTWO10/1024/
C
C UNPACK FLAGS
C ------------
C
      I = IDX
      MULT = I/NTWO20
      I = I - NTWO20*MULT
      MINHRM = I/NTWO10
      MAXHRM = I - NTWO10*MINHRM
      RETURN
      END
**==ibpack.f  
C
C FUNCTION 'IBPACK'
C =================
C
C RETURN THE LAMBDA BIN OR FILM PACK NUMBER OR FILM NUMBER DECODED FROM THE
C PACKED WORD IN THE REFLECTION ARRAY 'IA'
C
      FUNCTION IBPACK(ITYP,IAWORD,ERR)
C
C PARAMETERS
C
C        ITYP (I)   FLAG =1, RETURN LAMBDA BIN NUMBER
C                        =2, RETURN FILM PACK NUMBER
C                        =3, RETURN FILM NUMBER WITHIN PACK
C      IAWORD (I)   WORD FROM ARRAY CONTAINING PACKED BIN NUMBER,
C                   PACK NUMBER, FILM NUMBER AND SCALING REJECT FLAG.
C         ERR (O)   LOGICAL FLAG =.TRUE.  REJECT SCALING FLAG SET
C                                =.FALSE. REJECT FLAG CLEAR
C
C SPECIFICATION STATEMENTS
C ------------------------
C
      LOGICAL ERR
      DATA NTWO16/65536/
      DATA NTWO6/64/
C
C GET BIN/PACK NUMBER
C -------------------
C
      ERR = .FALSE.
      IF (MOD(IAWORD,2).NE.0) ERR = .TRUE.
      II = IAWORD/NTWO16
      IF (ITYP.EQ.1) THEN
         IBPACK = II
         RETURN
      END IF
      JJ = IAWORD - NTWO16*II
      II = JJ/NTWO6
      IF (ITYP.EQ.2) THEN
         IBPACK = II
         RETURN
      END IF
      JJ = JJ - NTWO6*II
      IBPACK = JJ/2
      RETURN
      END
**==ea06c.f  
      SUBROUTINE EA06C(A,VALUE,VECTOR,M,IA,IV,W)
C       18/03/70 LAST LIBRARY UPDATE
      REAL A(IA,M),VALUE(M),VECTOR(IV,M),W(*)
      M1 = M + 1
      W(1) = A(1,1)
      IF (M.LT.2) THEN
      ELSE IF (M.EQ.2) THEN
         W(2) = A(2,2)
         W(4) = A(2,1)
      ELSE
         CALL MC04B(A,W,W(M1),M,IA,W(M+M1))
      END IF
      CALL EA08C(W,W(M1),VALUE,VECTOR,M,IV,W(M+M1))
      IF (M.LE.2) RETURN
      DO 50 L = 1 , M
         DO 40 II = 3 , M
            I = M - II + 1
            IF (W(M1+I).NE.0.0) THEN
               PP = 0.
               I1 = I + 1
               DO 20 K = I1 , M
                  PP = PP + A(I,K)*VECTOR(K,L)
 20            CONTINUE
               PP = PP/(A(I,I+1)*W(M1+I))
               DO 30 K = I1 , M
                  VECTOR(K,L) = VECTOR(K,L) + PP*A(I,K)
 30            CONTINUE
            END IF
 40      CONTINUE
 50   CONTINUE
      RETURN
      END
**==mc04b.f  
      SUBROUTINE MC04B(A,ALPHA,BETA,M,IA,Q)
      DIMENSION A(IA,*),ALPHA(*),BETA(*),Q(*)
      ALPHA(1) = A(1,1)
      DO 30 J = 2 , M
         J1 = J - 1
         DO 20 I = 1 , J1
            A(I,J) = A(J,I)
 20      CONTINUE
         ALPHA(J) = A(J,J)
 30   CONTINUE
      M1 = M - 1
      M2 = M - 2
      DO 120 I = 1 , M2
         PP = 0.0
         I1 = I + 1
         DO 40 J = I1 , M
            PP = PP + A(I,J)**2
 40      CONTINUE
         PP1 = SQRT(PP)
         IF (A(I,I+1).LT.0) THEN
            BETA(I+1) = PP1
         ELSE
            BETA(I+1) = -PP1
         END IF
         IF (PP.GT.0) THEN
            H = PP - BETA(I+1)*A(I,I+1)
            A(I,I+1) = A(I,I+1) - BETA(I+1)
            DO 70 KI = I1 , M
               QJ = 0.0
               DO 50 KJ = I1 , KI
                  QJ = QJ + A(KJ,KI)*A(I,KJ)
 50            CONTINUE
               IF (KI.LT.M) THEN
                  I2 = KI + 1
                  DO 60 KJ = I2 , M
                     QJ = QJ + A(KI,KJ)*A(I,KJ)
 60               CONTINUE
               END IF
               Q(KI) = QJ/H
 70         CONTINUE
            BIGK = 0.0
            DO 80 KJ = I1 , M
               BIGK = BIGK + A(I,KJ)*Q(KJ)
 80         CONTINUE
            BIGK = BIGK/(2.0*H)
            DO 90 KJ = I1 , M
               Q(KJ) = Q(KJ) - BIGK*A(I,KJ)
 90         CONTINUE
            DO 110 KI = I1 , M
               DO 100 KJ = KI , M
                  A(KI,KJ) = A(KI,KJ) - Q(KI)*A(I,KJ) - Q(KJ)*A(I,KI)
 100           CONTINUE
 110        CONTINUE
         END IF
 120  CONTINUE
      DO 130 I = 2 , M
         H = ALPHA(I)
         ALPHA(I) = A(I,I)
         A(I,I) = H
 130  CONTINUE
      BETA(M) = A(M-1,M)
      RETURN
      END
**==ea08c.f  
      SUBROUTINE EA08C(A,B,VALUE,VEC,M,IV,W)
C       19/03/70 LAST LIBRARY UPDATE
      REAL A(M),B(M),VALUE(M),VEC(*),W(*)
      DATA EPS/  9.6E-7/,A34/0./
C     THIS USES QR ITERATION TO FIND THE EIGENVALUES AND EIGENVECTORS
C  OF THE SYMMETRIC TRIDIAGONAL MATRIX WHOSE DIAGONAL ELEMENTS ARE
C  A(I),I=1,M AND OFF-DIAGONAL ELEMENTS ARE B(I),I=2,M.  THE ARRAY
C  W IS USED FOR WORKSPACE AND MUST HAVE DIMENSION AT LEAST 2*M.
C  WE TREAT VEC AS IF IT HAD DIMENSIONS (IV,M).
      CALL EA09C(A,B,W(M+1),M,W)
C     SET VEC TO THE IDENTITY MATRIX.
      DO 30 I = 1 , M
         VALUE(I) = A(I)
         W(I) = B(I)
         K = (I-1)*IV + 1
         L = K + M - 1
         DO 20 J = K , L
            VEC(J) = 0.
 20      CONTINUE
         VEC(K+I-1) = 1.
 30   CONTINUE
      ITER = 0
      IF (M.EQ.1) RETURN
      N2 = M
C     EACH QR ITERATION IS PERFORMED OF ROWS AND COLUMNS N1 TO N2
 40   DO 50 II = 2 , N2
         N1 = 2 + N2 - II
         IF (ABS(W(N1)).LE.(ABS(VALUE(N1-1))+ABS(VALUE(N1)))*EPS)
     +       GO TO 60
 50   CONTINUE
      N1 = 1
 60   IF (N2.NE.N1) THEN
         ROOT = W(M+N2)
         ITER = ITER + 1
         N2M1 = N2 - 1
         A22 = VALUE(N1)
         A12 = A22 - ROOT
         A23 = W(N1+1)
         A13 = A23
         DO 80 I = N1 , N2M1
            A33 = VALUE(I+1)
            IF (I.NE.N2M1) A34 = W(I+2)
            S = SIGN(SQRT(A12*A12+A13*A13),A12)
            SI = A13/S
            CO = A12/S
            JK = I*IV + 1
            J1 = JK - IV
            J2 = J1 + MIN0(M,I+ITER) - 1
            DO 70 JI = J1 , J2
               V1 = VEC(JI)
               V2 = VEC(JK)
               VEC(JI) = V1*CO + V2*SI
               VEC(JK) = V2*CO - V1*SI
               JK = JK + 1
 70         CONTINUE
            IF (I.NE.N1) W(I) = S
            A11 = CO*A22 + SI*A23
            A12 = CO*A23 + SI*A33
            A13 = SI*A34
            A21 = CO*A23 - SI*A22
            A22 = CO*A33 - SI*A23
            A23 = CO*A34
            VALUE(I) = A11*CO + A12*SI
            A12 = -A11*SI + A12*CO
            W(I+1) = A12
            A22 = A22*CO - A21*SI
 80      CONTINUE
         VALUE(N2) = A22
         GO TO 40
      ELSE
         N2 = N2 - 1
         IF (N2.GT.1) GO TO 40
C     RAYLEIGH QUOTIENT
         DO 100 J = 1 , M
            K = (J-1)*IV
            XX = VEC(K+1)**2
            XAX = XX*A(1)
            DO 90 I = 2 , M
               XX = XX + VEC(K+I)**2
               XAX = XAX + VEC(K+I)*(2.*B(I)*VEC(K+I-1)+A(I)*VEC(K+I))
 90         CONTINUE
            VALUE(J) = XAX/XX
 100     CONTINUE
         RETURN
      END IF
      END
**==ea09c.f  
      SUBROUTINE EA09C(A,B,VALUE,M,OFF)
C       18/03/70 LAST LIBRARY UPDATE
      REAL A(M),B(M),VALUE(M),OFF(M)
      DATA A34/0./,EPS/9.6E-7/
      VALUE(1) = A(1)
      IF (M.EQ.1) RETURN
      DO 20 I = 2 , M
         VALUE(I) = A(I)
         OFF(I) = B(I)
 20   CONTINUE
C     EACH QR ITERATION IS PERFORMED OF ROWS AND COLUMNS N1 TO N2
      N2 = M
 30   IF (N2.LE.1) RETURN
      DO 40 II = 2 , N2
         N1 = 2 + N2 - II
         IF (ABS(OFF(N1)).LE.(ABS(VALUE(N1-1))+ABS(VALUE(N1)))*EPS)
     +       GO TO 50
 40   CONTINUE
      N1 = 1
 50   IF (N2.NE.N1) THEN
C     ROOT  IS THE EIGENVALUE OF THE BOTTOM 2*2 MATRIX THAT IS NEAREST
C     TO THE LAST MATRIX ELEMENT AND IS USED TO ACCELERATE THE
C     CONVERGENCE
         BB = (VALUE(N2)-VALUE(N2-1))*0.5
         CC = OFF(N2)*OFF(N2)
         SBB = 1.
         IF (BB.LT.0.) SBB = -1.
         ROOT = VALUE(N2) + CC/(BB+SBB*SQRT(BB*BB+CC))
         N2M1 = N2 - 1
         A22 = VALUE(N1)
         A12 = A22 - ROOT
         A23 = OFF(N1+1)
         A13 = A23
         DO 60 I = N1 , N2M1
            A33 = VALUE(I+1)
            IF (I.NE.N2M1) A34 = OFF(I+2)
            S = SQRT(A12*A12+A13*A13)
            SI = A13/S
            CO = A12/S
            IF (I.NE.N1) OFF(I) = S
            A11 = CO*A22 + SI*A23
            A12 = CO*A23 + SI*A33
            A13 = SI*A34
            A21 = CO*A23 - SI*A22
            A22 = CO*A33 - SI*A23
            A23 = CO*A34
            VALUE(I) = A11*CO + A12*SI
            A12 = -A11*SI + A12*CO
            OFF(I+1) = A12
            A22 = A22*CO - A21*SI
 60      CONTINUE
         VALUE(N2) = A22
      ELSE
         N2 = N2 - 1
      END IF
      GO TO 30
      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
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
**==headxx.f  
C
C SUBROUTINE 'HEADXX'
C ===================
C
C THIS SUBROUTINE WAS IS USED TO OUTPUT A HEADING
C (CUT DOWN FROM ORIGINAL VERSION WHICH ALLOWED UNDERLINING & BOLD PRINT)
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 TEXT
C                          =2, DOUBLE SPACED TEXT
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
**==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
**==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
**==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
**==rdge.f  
C
C SUBROUTINE 'RDGE'
C =================
C
C READ LAUE DATA FROM A GENERATE FILE FOR A NUMBER OF FILMS
C
      SUBROUTINE RDGE(LUN,IP,INTYP,SIGREJ,ISPAT,IA,NRIN,NREAD,
     *                NSREJ,NNEG,ALMIN,ALMAX,SIGMAX,MUL,NMINP)
C
C PARAMETERS
C
C       LUN (I)   LOGICAL UNIT NUMBER  FOR READING 'GE1' FILE
C        IP (I)   PACK NUMBER (1 TO NPACK)
C     INTYP (I)   INTEGRATION TYPE 1=BOX, 2=PROFILE FITTED
C    SIGREJ (I)   OMIT REFLECTIONS WITH I<SIGREJ*SIG(I)
C     ISPAT (I)   =1 INCLUDE SPATIAL OVERLAP DATA, =0 DO NOT
C        IA (I/O) ARRAY TO STORE REFLECTION DATA
C      NRIN (I/O) NUMBER OF REFLECTIONS STORED AFTER READ
C     NREAD (I/O) NUMBER OF REFLECTIONS READ (MEASURED VALUES)
C     NSREJ (I/O) COUNT OF I<SIGREJ*SIG(I) REJECTS
C     NNEG  (I/O) COUNT OF REFLECTIONS WITH I<=0
C     ALMIN (I/O) MINIMUM LAMBDA VALUE IN DATA
C     ALMAX (I/O) MAXIMUM LAMBDA VALUE IN DATA
C    SIGMAX (I/O) MAXIMUM SIGMA VALUE
C       MUL (I/O) ARRAY TO STORE MULTIPLESV REFLECTION DATA
C     NMINP (I/O) NUMBER OF MULTIPLES READ  
C
C PARAMETER STATEMENTS
C
C       MXR = MAXIMUM NUMBER OF REFLECTIONS
C       MXMUL = MAXIMUM NUMBER OF MULTIPLES WHICH MAY BE STORED
C
      PARAMETER(MXR=300000,MXMUL=100000)
C
C
C SPECIFICATION STATEMENTS
C ------------------------
C
      DIMENSION IA(6*MXR)
      DIMENSION MUL(6*MXMUL)
      DIMENSION CELL(6)
      LOGICAL REJ
      EQUIVALENCE (AAA,IAA)
      EQUIVALENCE (BBB,IBB)
      DATA NTWO16/65536/
      DATA NTWO8/256/
      DATA NTWO6/64/
      DATA NTWO20/1048576/
      DATA NTWO10/1024/
C
C CODE
C ----
C
      AKP = 0.105
      ATAU = 1.0
      SCIN = 100.0
      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) RETURN
      CALL GETINT(LUN,IREC,IFILM,INTYP,ISPAT,IH,IK,IL,ALBDA,II,ISIG,X,Y,
     +            REJ,MULT,MINHRM)
      IF (REJ) GO TO 20
      NREAD = NREAD + 1
      STHL2 = STHLSQ(IH,IK,IL)
      STH = ALBDA*SQRT(STHL2)
      THETA = ASIN(STH)
      IF (SIGREJ.GT.0.0) THEN
         IF (II.LT.SIGREJ*ISIG) THEN
            NSREJ = NSREJ + 1
            GO TO 20
         END IF
         IF (II.LE.0) THEN
            NNEG = NNEG + 1
            GO TO 20
         END IF
      ELSE
         IF (II.LE.0) NNEG = NNEG + 1
      END IF
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
      IF (MULT.EQ.1) THEN
         NRIN = NRIN + 1
         IF (NRIN.GT.MXR) THEN
C
C ERROR CONDITION
C
            WRITE (6,6010) MXR
            STOP
         ELSE
            JJ = 6*NRIN - 5
            IA(JJ) = NTWO16*(IH+128) + NTWO8*(IK+128) + IL + 128
            IA(JJ+1) = 10000.0*ALBDA
            IA(JJ+2) = 10000.0*THETA
            IA(JJ+3) = IAA
            IA(JJ+4) = ISIG
            IA(JJ+5) = NTWO6*IP + IFILM
            ALMAX = AMAX1(ALMAX,ALBDA)
            ALMIN = AMIN1(ALMIN,ALBDA)
            SIGMAX = AMAX1(FLOAT(ISIG),SIGMAX)
            GO TO 20
         END IF
C
C STORE MULTIPLES
C
      ELSE
         NMINP = NMINP + 1
         IF (NMINP.GT.MXMUL) THEN
            WRITE (6,6020) MXMUL
            STOP
         ELSE
            JJ = 6*NMINP - 5
            MUL(JJ) = NTWO16*(IH+128) + NTWO8*(IK+128) + IL + 128
            BBB = ALBDA 
            MUL(JJ+1) = IBB
            MUL(JJ+2) = NTWO6*IP + IFILM
            MUL(JJ+3) = IAA
            MUL(JJ+4) = ISIG
            MUL(JJ+5) = NTWO20*MULT + NTWO10*MINHRM
            SIGMAX = AMAX1(FLOAT(ISIG),SIGMAX)
            GO TO 20
         END IF
      END IF
C
C FORMAT STATEMENTS
C
 6010 FORMAT (/,' **MAXIMUM REFLECTION CAPACITY OF',I6,' EXCEEDED**')
 6020 FORMAT (/,' **MAXIMUM MULTIPLES CAPACITY OF',I6,' EXCEEDED**')
      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,MULT,MINMUL)
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  (NODAL FOR A  MULTIPLE)
C          IK (O)   "
C          IL (O)   "
C       ALBDA (O)   WAVELENGTH  (NODAL FOR A MULTIPLE)
C          II (O)   INTENSITY
C        ISIG (O)   SIGMA(INTENSITY)
C         REJ (O)   REJECT FLAG = .TRUE. REJECT, = .FALSE. DO NOT
C        MULT (O)   MULTIPLICITY
C      MINMUL (O)   MINIMUM HARMONIC FOR A MULTIPLE
C
C
C SPECIFICATION STATEMENTS
C ------------------------
C
      LOGICAL REJ, MULTIP
      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 SPATIAL OVERLAPS
C
      CALL CCPBYI(MULTF,IMULT,1)
      CALL CCPBYI(MHKL,IMULT,-4)
      MULT = NBITST(MULTF(1),0,5)
      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.
      IF (MULT.GT.1) THEN
         CALL NODHKL(IH,IK,IL,NH,NK,NL)
         MINMUL = 0
         IF (NH.NE.0) THEN
            MINMUL = IH/NH
         ELSE IF (NK.NE.0) THEN
            MINMUL = IK/NK
         ELSE IF (NL.NE.0) THEN
            MINMUL = IL/NL
         ENDIF
         IH = NH
         IK = NK
         IL = NL
         ALBDA = MINMUL*ALBDA
      END IF
      RETURN
      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
C LAST UPDATE:
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C******************************
C**       NODHKL             **
C******************************
C
C
      SUBROUTINE NODHKL (H, K, L, HN, KN, LN)
C
C
C Purpose: To calculate the nodal indices of a given reflection
C
C Author:  Based on code from GENLAUE in NODAL_INDEX routine by Clifton et al.
C
C Arguments:
C
      INTEGER   H
      INTEGER   K
      INTEGER   L
      INTEGER   HN
      INTEGER   KN
      INTEGER   LN
C
C H    (R)   h index
C K    (R)   k index
C L    (R)   l index
C HN   (W)   Nodal h index
C KN   (W)   Nodal k index
C LN   (W)   Nodal l index
C
C====== FUNCTION TYPES:
C
      INTEGER HCF
C
C====== LOCALS:
C
      INTEGER   HCF1, HCF2
      INTEGER   H1, K1, L1
C
C------------------------------------------------------------------------------
C
      H1 = H
      K1 = K
      L1 = L
      HCF1 = HCF(H1,K1)
      HCF2 = HCF(HCF1,L1)
      HN = ISIGN(H1/HCF2,H1)
      KN = ISIGN(K1/HCF2,K1)
      LN = ISIGN(L1/HCF2,L1)
C
      RETURN
      END
C LAST UPDATE:
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C*************************
C**       HCF           **
C*************************
C
C
      INTEGER    FUNCTION  HCF (II,JJ)
C
C
C Purpose: To calculate h.c.f. of 2 integers
C
C Author:  Based on code from GENLAUE by Clifton et al.
C
C Arguments:
C
      INTEGER    II
      INTEGER    JJ
C 
C II    (R)   First integer
C JJ    (R)   Second integer
C
C====== LOCALS:
C
      INTEGER    I, J, ISAV, K
C
C------------------------------------------------------------------------------
C
      I = II
      J = JJ
      IF (I.EQ.J) THEN
         HCF = I
         RETURN
      END IF

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

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

 20   K = MOD(I,J)

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

      I = J
      J = K
      GO TO 20

      END
C
C SUBROUTINE 'MULDIA'
C ===================
C
C THIS SUBROUTINE IS USED TO PRINT DETAILS OF THE MULTIPLES DATA FOR
C DIAGNOSTIC PURPOSES.
C
      SUBROUTINE MULDIA (MUL,NMULT)
C
C PARAMETERS
C
C        MUL (I)   INTEGER ARRAY CONTAINING THE MULTIPLES REFLECTION DATA.
C      NMULT (I)   NUMBER OF MULTIPLES
C
C SPECIFICATION STATEMENTS
C ------------------------
C
      DIMENSION MUL(*)
      DIMENSION IHKL(3)
      EQUIVALENCE (IAA,AAA)
      DATA NTWO16/65536/
C
C CODE
C ----
C
      IF (NMULT.LE.0) RETURN
      WRITE (7,6010)
      DO 100 I = 1, NMULT
         JJ = 6*I - 5
         CALL GETHKL (MUL(JJ),IHKL)
         IAA = MUL(JJ+1)
         CALL GETMUL(MUL(JJ+5),MULT,MINHRM,MAXHRM)
         KSYM = MUL(JJ+4)/NTWO16
         ALBDA = AAA
         WRITE(7,6020)IHKL(1),IHKL(2),IHKL(3),KSYM,ALBDA,MULT,MINHRM,MAXHRM
100   CONTINUE
      WRITE (7,'('' '')')
      RETURN
C
C FORMAT STATEMENTS
C -----------------
C
 6010 FORMAT(//,'MULTIPLES DATA',//,
     + '    NH   NK   NL  KSYM    NOD-LAM  MULT MINHRM MAXHRM',/)
 6020 FORMAT(' ',3I5,I6,F11.3,I6,2I7)
      END
C
C SUBROUTINE 'DECONV'
C ===================
C
C THIS SUBROUTINE IS USED TO DECONVOLUTE THE MULTIPLES DATA 
C
      SUBROUTINE DECONV (MUL,NMULT,SCMAX,IWMUL,NEGTYP,MDIAG)
C
C PARAMETERS
C
C        MUL (I)   INTEGER ARRAY CONTAINING THE MULTIPLES REFLECTION DATA.
C      NMULT (I)   NUMBER OF MULTIPLES
C      SCMAX (I)   MAXIMUM LAMBDA CURVE SCALING FACTOR FOR INCLUSION OF
C                  TERMS (0.0 = NO LIMIT)
C      IWMUL (I)   = 0, USE UNIT WEIGHTS IN EQUATIONS
C                  = 1, USE WEIGHTS OF 1/SIGMA(I)**2
C     NEGTYP (I)   = 0, ELIMINATE RESULTS FROM EQUATIONS WITH NEGATIVE RESULTS
C                  = 1, INCLUDE RECALCULATED POSITIVE COMPONENTS
C                  = 2, ONLY OUTPUT THESE
C      MDIAG (I)   = 0, DO NOT WRITE MULTIPLES DIAGNOSTICS FILE
C                  > 0, WRITE MULTIPLES DIAGNOSTICS FILE
C
C     PARAMETER MXEQN   MAXIMUM NO. OF EQUATIONS ALLOWED FOR A REFLECTION
C     PARAMETER MXTRM   MAXIMUM NO> OF TERMA ALLOWED IN A EQUATION
C
      PARAMETER (MXB=100,MXEQN=200,MXTRM=100)
C
C SPECIFICATION STATEMENTS
C ------------------------
C
      DIMENSION MUL(*)
      LOGICAL EQN(MXEQN)
      LOGICAL ALLOK, NEGVAL, TERMS, ERR
      INTEGER IHARM(MXTRM)
      COMMON /DECON/COEF(MXEQN,MXTRM),AFF(MXEQN)
      REAL    AXX(MXTRM)
      REAL    ASIGXX(MXTRM)
      REAL    ALM(MXEQN,MXTRM)
      REAL    DSPAC(MXEQN,MXTRM)
      REAL    F(MXTRM)
      REAL    ASIGFF(MXEQN)
      REAL    WT(MXEQN)
      REAL    ADATA(10)
      LOGICAL END, ABSNT, DONE, ENTRY
      INTEGER NDHKL(3)
      INTEGER IPACK
      EQUIVALENCE (IAA,AAA)
      COMMON /PACKS/NPACK,NPVAL(MXB),PACKSC(MXB),PACKTF(MXB),
     *CPSCAL(MXB),CPTFAC(MXB)
      COMMON /SYSA/ ISYS,KSYS(3)
      DATA NTWO16/65536/
      EXTERNAL IBPACK
      INTEGER IBPACK
C
C INITIALISE
C
      NUNIQ = 0
      NOK = 0
      NFEW = 0
      NSNMTX = 0
      NNEG = 0
      NPOS = 0
      NBADSC = 0
      END = .FALSE.
      IF (NMULT.LE.0) RETURN
      M1 = 1
      ILAST = MUL(1)
      K = 0
C
C STEP THROUGH THE REFLECTIONS
C
 20   K = K + 1
      IF (K.GT.NMULT) THEN
         END = .TRUE.
         GO TO 30
      END IF
      JJ = 6*K - 5
      IF (MUL(JJ).EQ.ILAST) GO TO 20
 30   M2 = K - 1
      MM = 6*M1 - 5
      CALL GETHKL(MUL(MM),NDHKL)
      NUNIQ = NUNIQ + 1
      NEQN = M2-M1+1
      IF (NEQN.GT.MXEQN) THEN
         WRITE(6,6100) MXEQN
         STOP
      END IF
      DO 35 I = 1, NEQN
         EQN(I) = .TRUE.
 35   CONTINUE
      NTERMS = 0
      MAXHM = 0
      STHL2 = STHLSQ(NDHKL(1),NDHKL(2),NDHKL(3))
      DO 36 M = M1, M2
         JP = 6*M - 5
         CALL GETMUL (MUL(JP+5),MULT,MINHRM,MAXHRM)
         IF (MAXHRM.GT.MAXHM) MAXHM = MAXHRM
         IAA = MUL(JP+3)
         AFF(M-M1+1) = AAA
         ISIG = MUL(JP+4)/NTWO16
         ASIGFF(M-M1+1) = ISIG
         IF (IWMUL.GT.0) THEN
            WT(M-M1+1) = 1.0/(FLOAT(ISIG)*FLOAT(ISIG))
         ELSE
            WT(M-M1+1) = 1.0
         END IF
 36   CONTINUE
      DO 60 IHRM = 1, MAXHM
         ENTRY = .FALSE.
         IH = IHRM*NDHKL(1)
         IK = IHRM*NDHKL(2)
         IL = IHRM*NDHKL(3)
         IF (ISYS.NE.0) THEN
            CALL SYSAB(IH,IK,IL,ABSNT)
            IF (ABSNT) GO TO 60
         END IF
         STHL2 = STHLSQ(IH,IK,IL)
         DO 50 M = M1 , M2
            IEQN = M - M1 + 1
            COEF(IEQN,NTERMS+1) = 0.0
            ALM(IEQN,NTERMS+1) = 0.0
            DSPAC(IEQN,NTERMS+1) = 0.0
            JP = 6*M - 5
            CALL GETMUL (MUL(JP+5),MULT,MINHRM,MAXHRM)
            IF (IHRM.LT.MINHRM.OR.IHRM.GT.MAXHRM) GO TO 50
            IAA = MUL(JP+1)
            ALNOD = AAA
            ALBDA = ALNOD/IHRM
            SCAL = CURVAL(ALBDA)
            IF (SCAL.GT.0.0.AND.(SCMAX.LE.0.0.OR.SCAL.LE.SCMAX))THEN
               IPACK = IBPACK(2,MUL(JP+2),ERR)
               IF (CPTFAC(IPACK).NE.0.0) THEN
                  SCPAK = CPSCAL(IPACK)*EXP(-2.0*CPTFAC(IPACK)*STHL2)
                     ELSE
                  SCPAK = CPSCAL(IPACK)
               END IF
               COEF(IEQN,NTERMS+1) = 1.0/(SCAL*SCPAK)
               ALM(IEQN,NTERMS+1)=ALBDA
               DSPAC(IEQN,NTERMS+1)=SQRT(1.0/(4.0*STHL2))
               ENTRY = .TRUE.
            ELSE
               COEF(IEQN,NTERMS+1) = -1.0
               EQN(IEQN) = .FALSE.
            END IF
 50      CONTINUE
         IF (ENTRY) THEN
            NTERMS = NTERMS + 1
            IHARM(NTERMS) = IHRM
            IF (NTERMS.GE.MXTRM) THEN
               WRITE (6,6110) MXTRM
               STOP
            END IF
         END IF 
 60   CONTINUE
C
C ELIMINATE EQUATIONS WITH MISSING SCALE FACTORS
C
      NELIM = 0
      III = 0
      DO 64 I = 1, NEQN
         IF (EQN(I)) THEN
            III = III + 1
         ELSE
            NELIM = NELIM + 1
            GO TO 64
         ENDIF
         IF (III.NE.I) THEN
            AFF(III) = AFF(I)
            ASIGFF(III) = ASIGFF(I)
            WT(III) = WT(I)
            DO 63 J = 1, NTERMS
               COEF(III,J) = COEF(I,J)
               ALM(III,J) = ALM(I,J)
               DSPAC(III,J) = DSPAC(I,J)
 63         CONTINUE
         END IF
 64   CONTINUE
      NEQN = NEQN - NELIM
      NBADSC = NBADSC + NELIM
C
C====== Eliminate any terms for which all coefficients are now all 0.0
C
      IF (NEQN.GT.0) THEN
         III = 0
         DO 75 I = 1, NTERMS
            TERMS = .FALSE.
            DO 72 J = 1, NEQN
               IF (COEF(J,I).NE.0.0) TERMS = .TRUE.
72          CONTINUE
            IF (TERMS) THEN
               III = III + 1
               IF (III.NE.I) THEN
                  DO 74 J=1,NEQN
                     COEF(J,III) = COEF(J,I)
                     ALM(J,III) = ALM(J,I)
                     DSPAC(J,III) = DSPAC(J,I)
                     IHARM(III) = IHARM(I)
74                CONTINUE
               END IF
            END IF
75       CONTINUE
      END IF
      NTERMS = III
      IF (NEQN.GE.NTERMS.AND.NTERMS.GT.0) THEN
         NOK =  NOK + 1
      ELSE
         NFEW = NFEW + 1
      END IF
      IF (MDIAG.GT.0.AND.NEQN.GE.NTERMS.AND.NTERMS.GT.0) THEN
         WRITE (7,6010) NDHKL(1), NDHKL(2), NDHKL(3), NTERMS,
     +   (IHARM(I),I=1,NTERMS)
         DO 65 I=1,NEQN
            WRITE (7,6020) AFF(I), ASIGFF(I),
     +      (COEF(I,J),ALM(I,J),DSPAC(I,J),J=1,NTERMS)
 65      CONTINUE
      END IF
      NEGVAL = .FALSE.
100   CALL EQNSOL (NTERMS,NEQN,COEF,MXEQN,AFF,WT,AXX,ASIGXX,IERR)
      IF (IERR.EQ.0) THEN
         ALLOK = .TRUE.
         DO 66 I=1,NTERMS
            IF (AXX(I).GT.0.0) THEN
               F(I)=SQRT(AXX(I))
            ELSE
               F(I)= 0.0
               ALLOK = .FALSE.
            END IF
 66      CONTINUE
         IF (MDIAG.GT.0) THEN
            WRITE(7,6040) (AXX(I),I=1,NTERMS)
            WRITE(7,6050) (ASIGXX(I),I=1,NTERMS)
         END IF
         IF (ALLOK) THEN
            DO 68 I = 1, NTERMS
               ADATA(1) = IHARM(I)*NDHKL(1)
               ADATA(2) = IHARM(I)*NDHKL(2)
               ADATA(3) = IHARM(I)*NDHKL(3)
               ADATA(4) = F(I)
               ADATA(5) = ASIGXX(I)/(2.0*F(I))
               IF (NEGVAL) THEN
                  IF (NEGTYP.GT.0) CALL LWREFL(2,ADATA)
               ELSE
                  IF (NEGTYP.LT.2) CALL LWREFL(2,ADATA)
               END IF
 68         CONTINUE
            IF (NEGVAL) THEN
               NNEG = NNEG + 1
            ELSE
               NPOS = NPOS + 1
            END IF
         ELSE
            IF (MDIAG.GT.0) THEN
               WRITE(7,6030)
            END IF
            NEGVAL = .TRUE.
            III = 0
            NTNEG = 0
            DO 250 I = 1, NTERMS
               IF (AXX(I).GE.0.0) THEN
                  III = III + 1
               ELSE
                  NTNEG = NTNEG + 1
                  GO TO 250
               END IF
               IF (III.NE.I) THEN
                  DO 240 J = 1, NEQN
                     COEF(J,III) = COEF (J,I)
                     ALM(J,III) = ALM(J,I)
                     DSPAC(J,III) = DSPAC(J,I)
                     IHARM(III) = IHARM(I)
 240              CONTINUE
               END IF
 250        CONTINUE
            NTERMS = NTERMS - NTNEG
            GO TO 100
         END IF
      ELSE IF (IERR.EQ.1) THEN
         NSNMTX = NSNMTX + 1
      END IF
      IF (.NOT.END) THEN
         M1 = K
         ILAST = MUL(JJ)
         GO TO 20
      END IF
      CALL LWCLOS(2,1)
C
C OUTPUT COUNTS
C
      WRITE(6,6180)
      WRITE(6,6190) SCMAX, NBADSC
      WRITE(6,6200)
      WRITE(6,6210) NUNIQ
      WRITE(6,6220) NOK
      WRITE(6,6230) NFEW
      WRITE(6,6300)
      WRITE(6,6310) NPOS
      WRITE(6,6320) NNEG
      WRITE(6,6330) NSNMTX
      RETURN
C
C FORMAT STATEMENTS
C -----------------
C
 6010 FORMAT(/,' MULTIPLE',3I5,' NTERMS =',I3,' HARMONICS: ',(10I5))
 6020 FORMAT(' I =',F10.2,' SIG(I) =',F6.2,' COEFFS(LAM[D]):',
     +(10(F7.3,'(',F6.3,'[',F6.3,'])')))
 6030 FORMAT(' ELIMINATE NEGATIVE TERMS AND RE-SOLVE')
 6040 FORMAT(' INTENSITIES: ',(10F10.1))
 6050 FORMAT(' SIGMA(I)''S : ',(10F10.1))
 6100 FORMAT(//,' **MAXIMUM NUMBER OF EQUATIONS FOR MULTIPLES',
     + 'DECONVOLUTION OF',I4,' EXCEEDED**')
 6110 FORMAT(//,' **MAXIMUM NUMBER OF TERMS FOR MULTIPLES',
     + 'DECONVOLUTION OF',I4,' EXCEEDED**')
 6180 FORMAT(///,' M U L T I P L E S    D E C O N V O L U T I O N',//)
 6190 FORMAT(' NUMBER OF EQUATIONS ELIMINATED BECAUSE LAMBDA',
     + ' NORMALISATION FACTOR OUT OF LAMBDA RANGE, < 0.0 or >',
     + F6.1,' =',I5)
 6200 FORMAT(//,' SETS OF MULTIPLES FOR ATTEMPTED DECONVOLUTION',/)
 6210 FORMAT(' NUMBER OF UNIQUE SETS  OF  MULTIPLES  FOUND:',I5)
 6220 FORMAT(' NUMBER OF UNIQUES WITH SUFFICIENT EQUATIONS:',I5)
 6230 FORMAT(' NUMBER OF UNIQUES WITH  TOO  FEW  EQUATIONS:',I5)
 6300 FORMAT(//,' SOLUTION OF THE HARMONIC DECONVOLUTION EQUATIONS',/)
 6310 FORMAT(' NUMBER OF SETS SOLVED WITH  POSITIVE  TERMS  ONLY:',I5)
 6320 FORMAT(' NUMBER OF SETS SOLVED AFTER ELIMINATING -VE TERMS:',I5)
 6330 FORMAT(' NUMBER OF UNSOLVED SETS  WITH  SINGULAR  MATRICES:',I5)
      END
C
C SUBROUTINE 'EQNSOL'
C ===================
C
C THIS SUBROUTINE IS USED TO SOLVE EQUATIONS (BOTH SIDES ALSO MULTIPLIED
C BY WEIGHTING FACTORS)
C
C C(1,1)*X(1) + C(1,2)*X(2) + ... C(1,N)*X(N) = G(1)
C C(2,1)*X(1) + C(2,2)*X(2) + ... C(2,N)*X(N) = G(2)
C      .            .                  .         .
C      .            .                  .         .
C C(M,1)*X(1) + C(M,2)*X(2) + ... C(M,N)*X(N) = G(M)
C
      SUBROUTINE EQNSOL (N,M,C,IDIM,G,WT,X,SIGX,IERR)
C
C PARAMETERS 
C
C        N (I)   NUMBER OF UNKNOWNS
C        M (I)   NUMBER OF EQUATIONS
C        C (I)   EQUATIONS COEFFICIENTS MATRIX
C     IDIM (I)   FIRST DIMENSION OF 'C' ARRAY
C        G (I)   OBSERVED VALUES
C       WT (I)   WEIGHTS OF OBSERVEED VALUES
C        X (O)   RETURNS X VALUES
C     SIGX (O)   RETURNS SIGMA(X) VALUES
C     IERR (O)   ERROR FLAG 0 = OK
C                           1 = SINGULAR MATRIX
C                           -100 = FEWER EQUATIONS THAN TERMS
C                           -200 = NO UNKNOWNS
C
C
C     PARAMETER MXEQN   MAXIMUM NO. OF EQUATIONS ALLOWED FOR A REFLECTION
C     PARAMETER MXTRM   MAXIMUM NO> OF TERMA ALLOWED IN A EQUATION
C
      PARAMETER (MXEQN=200,MXTRM=100)
C
C
C SPECIFICATION STATEMENTS
C ------------------------
C
      REAL C(IDIM,*),G(*),WT(*),X(*),SIGX(*)
      DOUBLE PRECISION A(MXEQN,MXTRM),AT(MXTRM,MXEQN)
      DOUBLE PRECISION AN(MXTRM,MXTRM),AI(MXTRM,MXTRM),RHS(MXTRM)
      DOUBLE PRECISION DET
C
C INITIALISE
C
      IERR = 0
      IF (N.LE.0) THEN
         IERR = -200
         RETURN
      END IF
      IF (M.LT.N) THEN
         IERR = -100
         RETURN
      END IF
C
C SET UP 'A' MATRIX
C
      DO 20 I = 1, M
         DO 10 J = 1, N
            A(I,J) = C(I,J)
 10      CONTINUE
 20   CONTINUE
C
C GET TRANSPOSE OF 'A' IN 'AT'
C
      DO 40 I = 1, N
         DO 30 J = 1, M
            AT(I,J) = A(J,I)
30       CONTINUE
40    CONTINUE
C
C GET WEIGHTED 'A'
C
      DO 45 I = 1, M
         DO 42 J = 1, N
            A(I,J) = WT(I)*A(I,J)
 42      CONTINUE
 45   CONTINUE
C
C SET NORMAL MATRIX
C
      DO 70 I = 1,N
         DO 60 J = 1,N
            AN(I,J) = 0.0
            DO 50 K = 1, M
               AN(I,J)=AN(I,J)+AT(I,K)*A(K,J)
 50         CONTINUE   
 60      CONTINUE
 70   CONTINUE
C
C GET RHS ELEMENTS
C
      DO 90 I =1,N
         RHS(I) = 0.0
         DO 80 K = 1, M
            RHS(I) = RHS(I) + AT(I,K)*WT(K)*G(K)
 80      CONTINUE
 90   CONTINUE
C
C INVERT NORMAL MATRIX
C
      CALL INVMAT(AN,AI,N,MXTRM,6,DET,IERR)
      IF (IERR.NE.0) RETURN
C
C GET SOLUTIONS
C
      DO 120 I = 1, N
         X(I) = 0.0
         DO 110 K = 1, N
            X(I) = X(I) + AI(I,K)*RHS(K)
 110     CONTINUE
 120  CONTINUE
C
C GET SIGMAS
C
      DO 130 I = 1, N
         SIGX(I) = AI(I,I)
         IF (SIGX(I).GE.0) THEN
            SIGX(I) = SQRT(SIGX(I))
         ELSE
            SIGX(I) = 0.0
         END IF
 130  CONTINUE
      RETURN
      END
C
C SUBROUTINE 'LSQSOL'
C ===================
C
      SUBROUTINE LSQSOL (N,M,C,IDIM,G,X,IERR)
C
C PARAMETERS 
C
C        N (I)   NUMBER OF UNKNOWNS
C        M (I)   NUMBER OF EQUATIONS
C        C (I)   EQUATIONS COEFFICIENTS MATRIX
C     IDIM (I)   FIRST DIMENSION OF 'C' ARRAY
C        G (I)   OBSERVED VALUES
C        X (O)   RETURNS X VALUES
C     IERR (O)   ERROR FLAG
C
C
C     PARAMETER MXEQN   MAXIMUM NO. OF EQUATIONS ALLOWED FOR A REFLECTION
C     PARAMETER MXTRM   MAXIMUM NO> OF TERMA ALLOWED IN A EQUATION
C
      PARAMETER (MXEQN=200,MXTRM=100)
C
C
C SPECIFICATION STATEMENTS
C ------------------------
C
      REAL C(IDIM,*),G(*),X(*)
      DOUBLE PRECISION FVECC(MXEQN),W((MXEQN+2*MXTRM+3)*MXTRM+2*MXEQN)
      DOUBLE PRECISION XX(MXTRM),SHFLIM,H,DAMP,FMIN,ACC
      EXTERNAL CALCI
C
C INITIALISE
C
      SHFLIM = 0.01
      H = 1.0
      DAMP = 0.5
      ACC = 0.001
      IERR = 0
      IF (M.LT.N) THEN
         IERR = -100
         RETURN
      END IF
      GAVE = 0.0
      DO 10 I = 1, M
         GAVE = GAVE + G(I)
10    CONTINUE
      GAVE = GAVE/M
      DO 20 I = 1, N
         XX(I) = GAVE/2.0
20    CONTINUE
C      CALL LSQMIN(M,N,XX,SHFLIM,50,H,DAMP,CALCI,FMIN,FVECC,IFLAG,W)
      CALL SSQMIN(M,N,XX,ACC,500,CALCI,FMIN,FVECC,IFLAG,W)
      FM = FMIN
      WRITE(6,6000)FM,IFLAG
      DO 30 I = 1, N
         X(I) = XX(I)*XX(I)
30    CONTINUE
      RETURN
6000  FORMAT(' LSQMIN RETURNS FMIN =',F7.3,'  IFLAG =',I8)
      END
C
C SUBROUTINE 'CALCI'
C ==================
C
C CALCULATE FUNCTIONS FOR LEAST SQUARES REFINEMENT OF MULTIPLES COMPONENTS
C
      SUBROUTINE CALCI(M,N,FVECC,X)
C
      PARAMETER (MXEQN=200,MXTRM=100)
C
C SPECIFICATION STATEMENTS
C ------------------------
C
      DOUBLE PRECISION FVECC(M),X(N)
      COMMON /DECON/COEF(MXEQN,MXTRM),AFF(MXEQN)
C
C CODE
C ----
C
      DO 100 I = 1, M
         FVECC(I) = 0.0
         DO 50 J = 1, N
            FVECC(I) = FVECC(I) + COEF(I,J)*X(J)*X(J)
50       CONTINUE
         FVECC(I) = FVECC(I) - AFF(I)
100   CONTINUE
      RETURN
      END
