      program locscl
c
c Robert H. Blessing
c Hauptman-Woodward Institute
c 73 High St.
c Buffalo, NY 14203, USA
c Tel. 716-856-9600, ext. 335
c blessing@hwi.buffalo.edu
c
c Dec. 1995,..., Feb. 1999
c
c Scaling by least-squares fitted local scale factors
c
c   q = <F1(hkl)/F2(hkl)>
c
c to improve the accuracy of structure factor magnitude differences
c
c   Delta(F) = F1 - q*F2
c
c for SIR F1(native)/F2(derivative) or SAS F1(+h+k+l)/F2(-h-k-l) pairs.
c
c The local scale factors are evaluated in approximately equidimensional
c local blocks of
c
c   n = {[2*m(h) + 1] X [2*m(k) + 1] X [2*m(l) + 1]} - 1
c
c reciprocal lattice points surrounding, but excluding, each point hkl.
c
c B.W. Matthews and E.W. Czerwinski (1975).  Local Scaling:  a Method to
c Reduce Systematic Errors in Isomorphous Replacement and Anomalous
c Scattering Measurements.  Acta Cryst. A31, 480-497.
c
c W.A. Hendrickson (1988).  Structural Analysis by the Method of
c Multiwavelength Anomalous Diffraction.  In "Crystallographic
c Computing 4", edited by N.W. Isaacs and M.R. Taylor, pp. 97-108.
c Oxford:  Oxford Univ. Press.  And references cited therein.
c
      character atime*8,adate*9,title*72,type*3,ptgp*5,
     & file1*40,file2*40,form1*11,form2*11
      parameter (imax=999999)
      dimension data(10,0:imax)
      dimension cell(6),ginv(3,3)
      data ptgp,itype1,itype2,cutoff,cell,mh,mk,ml,i012,zlimit
     & /'     ',14*0/
c      atime='hh:mm:ss'
c      adate='dd-mmm-yy'
      call time (atime)
      call date (adate)
c
c Read control data file "locscl.dat".
c
      io1=10
      io2=20
      io3=30
      ilp=60
      open (unit=io1,file='locscl.dat',status='old')
      read (io1,'(a)') title
      read (io1,'(a)') type
      if (type.eq.'sir') type='SIR'
      if (type.eq.'sas') type='SAS'
      read (io1,'(a)') file1
      if (type.eq.'SIR') read (io1,'(a)') file2
      read (io1,*,end=1) itype1
      if (type.eq.'SIR') read (io1,*,end=1) itype2
      read (io1,*) cutoff
      read (io1,'(a)') ptgp
      read (io1,*) (cell(i),i=1,6)
      read (io1,*,end=1) mh,mk,ml
      read (io1,*,end=1) i012
      read (io1,*,end=1) zlimit
 1    close (unit=io1,status='keep')
c
c Maximum local block dimensions are
c (2*mh + 1)*(2*mk + 1)*(2*ml + 1) = 25*25*25.
c
      if (max(mh,mk,ml).gt.12) then
        write (6,'(/1x,''Choose smaller mh, mk, and/or ml .le. 12.'')')
        stop
      end if
c-----------------------------------------------------------------------
c
c Control data:
c ------------
c
c title   = Job title/crystal name
c
c type    = SAS or SIR = Type of data set pair
c
c  file1  = Name of SAS reflections file OR SIR Native reflections file
c (file2  = Name of SIR Derivative reflections file)
c
c  itype1 = File type (see below) of file1
c (itype2 = File type of file2)
c
c Input reflection data files:
c ---------------------------
c
c   For SAS data, one file:
c       ------------------
c   file1 containing both F(+h+k+l) and F(-h-k-l) in separate reflection
c   records, in any order.
c
c   For SIR data, two files:
c       -------------------
c   file1 for the native data, and
c   file2 for the derivative data.  In both files the reflection records
c   may be in any order, and the two files need not be of the same type.
c
c
c   The files are read in subroutine "read1" according to the values of
c   itype1 and itype2:
c   ------     ------
c
c   itype = 0  ASCII, formatted     h,k,l,Fsq,sigFsq,F,sigF,E,sigE
c           1    "        "         h,k,l,Fsq,sigFsq,F,sigF
c           2    "        "         h,k,l,Fsq,sigFsq       
c           3    "        "         h,k,l,           F,sigF,E,sigE
c           4    "        "         h,k,l,           F,sigF       
c           5  Binary, unformatted  h,k,l,Fsq,sigFsq,F,sigF,E,sigE
c           6    "          "       h,k,l,Fsq,sigFsq,F,sigF
c           7    "          "       h,k,l,Fsq,sigFsq
c           8    "          "       h,k,l,           F,sigF,E,sigE
c           9    "          "       h,k,l,           F,sigF
c
c
c cutoff = minimum F/sigma(F) threshold for input data to be included in
c          the scaling calculations and the locally scaled output files.
c
c          Default is cutoff = 3.  To apply cutoff = 0, enter 
c          cutoff .lt. 0.
c
c ptgp   = Point group symbol
c
c          Must be one of the following case-sensitive character values.
c
c          Tricl. Monocl. Ortho. Tetrag. Trig. Rhomb. Cub.
c          ------ ------- ------ ------- ----- ------ ----
c          1      2       222    4       3     R3     23
c          -1     m       mm2    -4      -3    R-3    m3
c                 2/m     mmm    4/m   
c                                        312   R32    432
c                                422     321   R3m    -43m
c                                4mm     31m   R-3m   m3m
c                                -42m    3m1
c                                -4m2    -31m
c                                4/mmm   -3m1
c 
c                                        Hexag.
c                                        ------
c                                        6
c                                        -6
c                                        6/m
c
c                                        622
c                                        6mm
c                                        -6m2
c                                        -62m
c                                        6/mmm
c
c cell     = lattice parameters (a, b, c, alpha, beta, gamma)
c
c mh,mk,ml = semidimensions  m(h), m(k), m(l)  if fixed semidimensions
c            are desired for the moving local block of
c
c            n = [2*m(h) + 1]*[2*m(k) + 1]*[2*m(l) + 1] - 1
c
c            reciprocal lattice points surrounding, but excluding, each
c            point hkl.
c
c            The default (mh = mk = ml = 0) is to automatically vary and
c            statistically optimize the semidimensions.
c
c            To bypass the local scale factor fitting and apply constant
c            local scale factors of unity, supply  mh = mk = ml = -1.
c
c            If fixed semidimensions are chosen, choose
c            mh, mk, ml .le. 12.
c
c i012     = 0  Keep constant the  mean of F1 and F2.
c            1   "      "      "  value of F1.
c            2   "      "      "    "    " F2.
c
c            The default  i012 = 0  is recommended.
c
c zlimit   = limiting value of
c
c            abs(Delta)          abs(Delta)
c            ---------- = -----------------------
c              sigma      1.25*median[abs(Delta)]
c
c                                     abs[F1 - F2 - median(F1 - F2)]
c                       = ------------------------------------------- ,
c                         1.25*median{abs[F1 - F2 - median(F1 - F2)]}
c
c            for reflection pairs to be accepted for the output file.
c            The variable zlimit is used to exclude data with unreliably
c            large values of  abs(F1 - F2)  in the tails of the data-set
c            Delta(F) distribution.
c
c            The zlimit test assumes that the data-set
c            z = Delta/sigma distribution should approximate a
c            zero-mean, unit-variance normal distribution, for which
c            values of  z .lt. -zlimit  or  z .gt. +zlimit  are
c            extremely improbable.
c
c            With the default  zlimit = 0 ,  the zlimit exclusion is not
c            applied.  A recommended trial value is  zlimit = 6.
c
c-----------------------------------------------------------------------
      open (unit=ilp,file='locscl.lp',status='new')
      write (ilp,'(/1x,''Program locscl.  '',a,'', '',a//1x,
     &''Least-squares local scale factors''/1x,
     &''  q = <F1(hkl)/F2(hkl)>''/1x,
     &''to improve the accuracy of structure factor magnitude differe'',
     &''nces''/1x,
     &''  Delta(F) = F1 - q*F2''/1x,
     &''for SIR [F1(native)/F2(derivative)] or SAS [F1(+h+k+l)/F2(-h-'',
     &''k-l)] pairs.''//1x,
     &''The local scale factors are evaluated in approximately equidi'',
     &''mensional''/1x,''local blocks of''/1x,
     &''  n = {[2*delta(h) + 1] X [2*delta(k) + 1] X [2*delta(l) + 1]'',
     &''} - 1''/1x,
     &''reciprocal lattice points surrounding, but excluding, each po'',
     &''int hkl.'')') atime,adate
      write (ilp,'(/1x,''Job title:''/1x,a)') title
      write (ilp,'(/1x,''Lattice parameters:''/1x,
     &''         a         b         c     alpha      beta     gamma''/
     &1x,3f10.4,3f10.3)') cell
      call metric (cell,ginv)
      if (cutoff.eq.0) cutoff=3
      if (cutoff.lt.0) cutoff=0
c
c For SIR data, open native and derivative reflection data files.
c
      if (type.eq.'SIR') then
        write (ilp,'(/1x,
     &  ''SIR File-1 [F1(Native)]:''/1x,a//1x,
     &  ''SIR File-2 [F2(Derivative)]:''/1x,a//1x,
     &  ''File-1 type:  itype1 = '',i2/1x,
     &  ''File-2 type:  itype2 = '',i2)') file1,file2,itype1,itype2
        form1='formatted'
        form2='formatted'
        if (itype1.ge.5) form1='unformatted'
        if (itype2.ge.5) form2='unformatted'
        open (unit=io1,file=file1,form=form1,status='old')
        open (unit=io2,file=file2,form=form2,status='old')
c
c For SAS data, prepare separate +h+k+l and -h-k-l reflection data files
c from one input file.
c
      else if (type.eq.'SAS') then
        write (ilp,'(/1x,
     &  ''SAS Input File [F1(+h+k+l) and F2(-h-k-l)]:''/1x,a//1x,
     &  ''File type:  itype1 = '',i2)') file1,itype1
        form1='formatted'
        if (itype1.ge.5) form1='unformatted'
        open (unit=io3,file=file1,status='old',form=form1)
        open (unit=io2,file='data.neg',status='new',form='unformatted')
        open (unit=io1,file='data.pos',status='new',form='unformatted')
        call bijvoet (io1,io2,io3,ilp,itype1,ptgp)
        close (unit=io3,status='keep')
        rewind io2
        rewind io1
        itype1=8
        itype2=8
      end if
c
c Store reflection data pairs in memory.
c
      call datain (io1,io2,itype1,itype2,cutoff,ihmin,ikmin,ilmin,nk,nl,
     & ginv,imax,data,jmin,jmax,type,ptgp,ilp)
      if (type.eq.'SIR') then
        close (unit=io1,status='keep')
        close (unit=io2,status='keep')
      else if (type.eq.'SAS') then
        close (unit=io1,status='delete')
        close (unit=io2,status='delete')
      end if
c
c Global pre-scaling of SIR pairs
c
      if (type.eq.'SIR') call pscale (imax,data,jmin,jmax,cutoff,ginv,
     & ptgp,ihmin,ikmin,ilmin,nk,nl,ilp)
c
c Local scaling
c
      write (ilp,'(//1x,''Local Scaling''/1x,''-------------'')')
      if (mh.lt.0.or.mk.lt.0.or.ml.lt.0) then
c
c Apply fixed local scale factors  q = 1.
c
        write (ilp,'(/1x,''All local scale factors fixed at  q = 1.'')')
        go to 15
      else if (mh.eq.0.and.mk.eq.0.and.ml.eq.0) then
c
c Optimize the local block semidimensions mh, mk, and ml to maximize
c   delta = sqrt(rmsd**2 - esd**2) - esd ,
c where
c   rmsd  = <(q - <q>)**2>**(1/2) ,
c   esd   = <sigma(q)**2>**(1/2) ,
c and
c   sigma(q) = sqrt{sum w*[(F1/F2) - q]**2/[(n - 1)*sum w]} .
c
c Both rmsd and esd tend to decrease with increasing local block size.
c Rmsd**2 estimates the total, systematic-plus-random, variance.
c Esd is the rms error of fit, and esd**2 estimates the random variance.
c Thus delta estimates the difference between the systematic and random
c local variation.
c
        write (ilp,'(1x,
     &''Optimize the local block semidimensions mh, mk, and ml to max'',
     &''imize''/1x,
     &''  delta = sqrt(rmsd**2 - esd**2) - esd ,''/1x,''where''/1x,
     &''  rmsd  = <(q - <q>)**2>**(1/2) ,''/1x,
     &''  esd   = <sigma(q)**2>**(1/2) ,''/1x,''and''/1x,
     &''  sigma(q) = sqrt{sum w*[(F1/F2) - q]**2/[(n - 1)*sum w]} .'')')
        write (ilp,'(/1x,     
     &''Both rmsd and esd tend to decrease with increasing local bloc'',
     &''k size.''/1x,
     &''Rmsd**2 estimates the total, systematic-plus-random, variance.''
     &/1x,
     &''Esd is the rms error of fit, and esd**2 estimates the random '',
     &''variance.''/1x,
     &''Thus delta estimates the difference between the systematic an'',
     &''d random''/1x,''local variation.'')')
        write (ilp,'(/1x,
     &'' m  mh mk ml  nh nk nl     n    <q>  rmsd   esd delta delta(m'',
     &'')-''/1x,
     &''                                        %     %     %  delta('',
     &''m-1)''/1x,
     &'' -  -- -- --  -- -- --     -    ---  ----   --- ----- -------'',
     &''----'')')
        a=cell(1)
        b=cell(2)
        c=cell(3)
        amax=max(a,b,c)
        a=a/amax
        b=b/amax
        c=c/amax
        delta2=0
        delta1=0
        delta=0
        m=0
        do i=1,12
          delta2=delta1
          delta1=delta
          m=m+1
          mh=max(nint(m*a),1)
          mk=max(nint(m*b),1)
          ml=max(nint(m*c),1)
          call qfit (imax,data,jmin,jmax,ihmin,ikmin,ilmin,nk,nl,
     &    mh,mk,ml,nrej,ndata,qmin,sigmin,qmax,sigmax,qmean,rmsd,esd,
     &    wqmean,wrmsd,wesd)
          delta=(sqrt(max((wrmsd**2-wesd**2),0.0))-wesd)
          write (ilp,'(1x,i2,1x,3i3,1x,3i3,i6,f7.4,3f6.2,e12.3)')
     &    m,mh,mk,ml,
     &    (2*mh+1),(2*mk+1),(2*ml+1),((2*mh+1)*(2*mk+1)*(2*ml+1))-1,
     &    wqmean,100*wrmsd/wqmean,100*wesd/wqmean,100*delta/wqmean,
     &    delta-delta1
          if (m.gt.2.and.delta1.gt.max(delta2,delta)) then
            m=m-1
            mh=max(nint(m*a),1)
            mk=max(nint(m*b),1)
            ml=max(nint(m*c),1)
            go to 5
          end if
        end do
        go to 10
      else if (mh.gt.0.or.mk.gt.0.or.ml.gt.0) then
c
c Use fixed, input local block semidimensions.
c
        write (ilp,'(1x,
     &  ''Fixed, input local block semidimensions:''//1x,
     &  ''   mh mk ml  nh nk nl     n''/1x,
     &  ''   -- -- --  -- -- --     -''/1x,2x,3i3,1x,3i3,i6)')
     &  mh,mk,ml,
     &  (2*mh+1),(2*mk+1),(2*ml+1),
     &  ((2*mh+1)*(2*mk+1)*(2*ml+1))-1
        go to 5
      end if
 5    continue
c
c Use either input or optimized local block dimensions in a final call
c to subroutine qfit.
c
      call qfit (imax,data,jmin,jmax,ihmin,ikmin,ilmin,nk,nl,mh,mk,ml,
     & nrej,ndata,qmin,sigmin,qmax,sigmax,qmean,rmsd,esd,wqmean,wrmsd,
     & wesd)
 10   continue
      write (ilp,'(/1x,
     &''Local Scaling Results''/1x,''---------------------''/1x,
     &''Local blocks of up to''/1x,
     &''  n = {[2*m(h) + 1] X [2*m(k) + 1] X [2*m(l) + 1]} - 1''/1x,
     &''reciprocal lattice points surrounding, but excluding, each po'',
     &''int hkl''/1x,
     &''  m(h) = '',i3,''    n(h) = '',i3/1x,
     &''  m(k) = '',i3,''    n(k) = '',i3/1x,
     &''  m(l) = '',i3,''    n(l) = '',i3,''    n = '',i6)')
     & mh,(2*mh+1),
     & mk,(2*mk+1),
     & ml,(2*ml+1),((2*mh+1)*(2*mk+1)*(2*ml+1))-1
      write (ilp,'(/1x,
     &''Nrej = '',i8,''  reflection pairs rejected due to local scale'',
     &'' factor''/1x,
     &''       '',8x,''  error-of-fit uncertainty''//1x,
     &''N    = '',i8,''  locally scaled F1(hkl) and F2(hkl) reflectio'',
     &''n pairs'')') nrej,ndata
      write (ilp,'(/1x,
     &''qmin = '',f8.4,''  sigma(qmin) = '',f8.4/1x,
     &''qmax = '',f8.4,''  sigma(qmax) = '',f8.4//1x,
     &''                          Unit-Weighted  Reciprocal-Variance-'',
     &''Weighted''/1x,
     &''                          -------------  --------------------'',
     &''--------''/1x,
     &''mean = <q>                   = '',2f8.4/1x,
     &''rmsd = <(q - <q>)**2>**(1/2) = '',2f8.4/1x,
     &''esd  = <sigma(q)**2>**(1/2)  = '',2f8.4//1x,
     &''delta = 100%*[sqrt(wrmsd**2 - wesd**2) - wesd]/wavg = '',f6.2,
     &''%'')')
     & qmin,sigmin,qmax,sigmax,qmean,wqmean,rmsd,wrmsd,esd,wesd,
     & 100*(sqrt(wrmsd**2-wesd**2)-wesd)/wqmean
 15   continue
c
c Apply local scale factors.
c
      call qscale (imax,data,jmin,jmax,ihmin,ikmin,ilmin,nk,nl,io1,io2,
     & i012,type,ptgp,ndata,nc,ilp)
c
c Compile locally scaled Delta(F) distribution statistics and histogram,
c and write output reflection files.
c
      call output (io1,io2,ilp,type,ptgp,ndata,nc,zlimit,atime,adate,
     & title)
c
c Compile subset agreement statistics.
c
      call agree (io1,ilp,ginv)
c
c Describe the output files.
c
      write (ilp,'(//1x,''Output Files''/1x,''------------'')')
      if (i012.eq.0) write (ilp,'(1x,
     &''Output data files have F1 and locally scaled F2 re-scaled to '',
     &''give''/1x,''the same mean as the unscaled F1 and F2.'')')
      if (i012.eq.1) write (ilp,'(1x,''Output data files have unchang'',
     &''ed F1 and locally scaled F2.'')')
      if (i012.eq.2) write (ilp,'(1x,''Output data files have locally'',
     &'' scaled F1 and unchanged F2.'')')
      write (ilp,'(/1x,
     &''The "data.locscl" output file is sorted on the Laue-group uni'',
     &''que indices.''/1x,
     &''    -------------''/1x,
     &''The file contains header records that label the data records:''/
     &1x,
     &''  h,k,l,F1,sigma(F1),F2,sigma(F2),E1,sigma(E1),E2,sigma(E2),''/
     &1x,''where''/1x,
     &''  1 designates the Native     or the +h+k+l data set, and''/1x,
     &''  2     "       "  Derivative  "  "  -h-k-l   "   " .'')')
      write (ilp,'(/1x,
     &''The "sort.locscl" output file is sorted on''/1x,
     &''    -------------''//1x,
     &''              abs[Delta(E)]              abs(E1 - E2)''/1x,
     &''  abs[z(E)] = ------------- = -------------------------------'',
     &''--''/1x,
     &''              sigma(Delta)    sqrt[sigma(E1)**2 + sigma(E2)**'',
     &''2]''//1x,
     &''or, in the absence of E values, the analogous z(F).''/1x,
     &''The file contains header records that label the data records:''/
     &1x,
     &''  h,k,l,F1,sigma(F1),F2,sigma(F2),E1,sigma(E1),E2,sigma(E2),z'',
     &''(E or F).'')')
      if (type.eq.'SAS') write (ilp,'(/1x,
     &''The "hphn.locscl" output file has adjacent separate records f'',
     &''or the SAS''/1x,
     &''    -------------''/1x,''Bijvoet pairs.''/1x,
     &''The file contains header records that label the data records:''/
     &1x,
     &''  +h,+k,+l,F1,sigma(F1),E1,sigma(E1);''/1x,
     &''  -h,-k,-l,F2,sigma(F2),E2,sigma(E2).'')')
      write (ilp,'(/)')
      stop 'Program locscl finis'
      end
c-----------------------------------------------------------------------
      subroutine agree (io1,ilp,ginv)
      dimension ginv(3,3)
      dimension d(0:20),n(20),sumdf(20),sumf1(20),sumf2(20),sumde(20),
     & sumae(20),zmin(10)
      real*8 sumdf,sumf1,sumf2,sumde,sumae
      write (ilp,'(//1x,
     &''Agreement Statistics''/1x,''--------------------''/1x,
     &''Q(F) = <F1>/<F2>''/1x,
     &''R(F) = <abs(F1 - F2)>/<0.5*(F1 + F2)>''/1x,
     &''R(E) = <abs(E1 - E2)>/<0.5*(E1 + E2)>''//1x,
     &''Statistics in equal-volume resolution shells with''/1x,
     &''dmax .gt. d .ge. dmin:''//1x,
     &''   dmax   dmin      N   Q(F)   R(F)   R(E)''/1x,
     &''   ----   ----      -   ----   ----   ----'')')
      open (unit=io1,file='data.locscl',status='old')
c
c Read dummy header records.
c
      do i=1,4
        read (io1,'(a1)') dummy
      end do
      dlimit=1e9
      ndata=0
      do i=1,1e6
        read (io1,*,end=1) ih,ik,il,f1,s1,f2,s2,e1,t1,e2,t2
        dlimit=min(dlimit,1/(2*sinthl(ih,ik,il,ginv)))
        ndata=ndata+1
      end do
 1    continue
      do i=20,1,-1
        d(i)=1/((float(i)/20)**(1/3.0)*(1/dlimit))
      end do
      do i=1,20
        n(i)=0
        sumdf(i)=0
        sumf1(i)=0
        sumf2(i)=0
        sumde(i)=0
        sumae(i)=0
      end do
      rewind io1
      do i=1,4
        read (io1,'(a1)') dummy
      end do
      do j=1,ndata
        read (io1,*) ih,ik,il,f1,s1,f2,s2,e1,t1,e2,t2
        di=1/(2*sinthl(ih,ik,il,ginv))
        do i=1,20
          if (di.ge.d(i)) then
            n(i)=n(i)+1
            sumdf(i)=sumdf(i)+abs(f1-f2)
            sumf1(i)=sumf1(i)+f1
            sumf2(i)=sumf2(i)+f2
            sumde(i)=sumde(i)+abs(e1-e2)
            sumae(i)=sumae(i)+e1+e2
            go to 2
          end if
        end do
 2      continue
      end do
      d(0)=999.99
      do i=1,20
        qf=0
        rf=0
        re=0
        if (sumf2(i).ne.0) qf=sumf1(i)/sumf2(i)
        if (sumf1(i)+sumf2(i).ne.0)
     &    rf=sumdf(i)/(0.5*(sumf1(i)+sumf2(i)))
        if (sumae(i).ne.0) re=sumde(i)/(0.5*sumae(i))
        write (ilp,'(1x,2f7.2,i7,3f7.3)') d(i-1),d(i),n(i),qf,rf,re
      end do
      write (ilp,'(/1x,
     &''Statistics for subsets with Y .ge. Ymin, where''/1x,
     &''Y(E) = min{[E1/sigma(E1)], [E2/sigma(E2)]},''/1x,
     &''or, in the absence of E-values, the analogous Y(F):''//1x,
     &''   Ymin      N   Q(F)   R(F)   R(E)''/1x,
     &''   ----      -   ----   ----   ----'')')
      data zmin /0, 1, 2, 3, 4, 6, 8, 10, 30, 100/
      do i=1,10
        n(i)=0
        sumdf(i)=0
        sumf1(i)=0
        sumf2(i)=0
        sumde(i)=0
        sumae(i)=0
      end do
      rewind io1
      do i=1,4
        read (io1,'(a1)') dummy
      end do
      read (io1,*) ih,ik,il,f1,s1,f2,s2,e1,t1,e2,t2
      if (e1.gt.0.and.e2.gt.0) then
        ie=1
      else
        ie=0
      end if
      rewind io1
      do i=1,4
        read (io1,'(a1)') dummy
      end do
      do j=1,ndata
        read (io1,*) ih,ik,il,f1,s1,f2,s2,e1,t1,e2,t2
        if (ie.eq.1) then
          zi=min(e1/t1,e2/t2)
        else
          zi=min(f1/s1,f2/s2)
        end if
        do i=1,10
          if (zi.ge.zmin(i)) then
            n(i)=n(i)+1
            sumdf(i)=sumdf(i)+abs(f1-f2)
            sumf1(i)=sumf1(i)+f1
            sumf2(i)=sumf2(i)+f2
            sumde(i)=sumde(i)+abs(e1-e2)
            sumae(i)=sumae(i)+e1+e2
          end if
        end do
      end do
      do i=1,10
        qf=0
        rf=0
        re=0
        if (sumf2(i).ne.0) qf=sumf1(i)/sumf2(i)
        if (sumf1(i)+sumf2(i).ne.0)
     &    rf=sumdf(i)/(0.5*(sumf1(i)+sumf2(i)))
        if (sumae(i).ne.0) re=sumde(i)/(0.5*sumae(i))
        write (ilp,'(1x,2i7,3f7.3)') int(zmin(i)),n(i),qf,rf,re
      end do
      write (ilp,'(/1x,
     &''Statistics for subsets with abs(Z) .ge. Zmin, where''/1x,
     &''Z(E) = (E1 - E2)/sqrt[sigma(E1)**2 + sigma(E2)**2],''/1x,
     &''or, in the absence of E-values, the analogous Z(F):''//1x,
     &''   Zmin      N   Q(F)   R(F)   R(E)''/1x,
     &''   ----      -   ----   ----   ----'')')
      do i=1,10
        n(i)=0
        sumdf(i)=0
        sumf1(i)=0
        sumf2(i)=0
        sumde(i)=0
        sumae(i)=0
      end do
      rewind io1
      do i=1,4
        read (io1,'(a1)') dummy
      end do
      do j=1,ndata
        read (io1,*) ih,ik,il,f1,s1,f2,s2,e1,t1,e2,t2
        if (ie.eq.1) then
          zi=(e1-e2)/sqrt(t1**2+t2**2)
        else
          zi=(f1-f2)/sqrt(s1**2+s2**2)
        end if
        do i=1,10
          if (abs(zi).ge.zmin(i)) then
            n(i)=n(i)+1
            sumdf(i)=sumdf(i)+abs(f1-f2)
            sumf1(i)=sumf1(i)+f1
            sumf2(i)=sumf2(i)+f2
            sumde(i)=sumde(i)+abs(e1-e2)
            sumae(i)=sumae(i)+e1+e2
          end if
        end do
      end do
      do i=1,10
        qf=0
        rf=0
        re=0
        if (sumf2(i).ne.0) qf=sumf1(i)/sumf2(i)
        if (sumf1(i)+sumf2(i).ne.0)
     &    rf=sumdf(i)/(0.5*(sumf1(i)+sumf2(i)))
        if (sumae(i).ne.0) re=sumde(i)/(0.5*sumae(i))
        write (ilp,'(1x,2i7,3f7.3)') int(zmin(i)),n(i),qf,rf,re
      end do
      close (unit=io1,status='keep')
      return
      end
c-----------------------------------------------------------------------
      subroutine bijvoet (io1,io2,io3,ilp,itype,ptgp)
c
c Prepare separate files of +h+k+l reflection data and Bijvoet-related
c -h-k-l anti-reflection data from a file of (acentric) point-group
c unique data.
c
      character ptgp*5,laue*5
      integer h,hp,hn
      call decode (ptgp,laue,iptgp,ilaue)
      write (ilp,'(/1x,''Point group:  '',a/1x,''Laue group:   '',a)')
     & ptgp,laue
      nc=0
      n1=0
      n2=0
      iend=0
      do i=1,1e9
        call read1 (itype,io3,iend,h,k,l,f,sigf,e,sige)
        if (iend.ne.0) go to 9
c
c Transform to point-group unique indicies.
c
        call equiv (ptgp,h,k,l)
c
c Acentric (ic = 0) or centric (ic = 1) reflection?
c
        call csym (ptgp,h,k,l,ic)
        if (ic.eq.1) then
c
c Write centric reflections to both the +h+k+l reflection file (1) and
c the -h-k-l anti-reflection file (2).
c
          nc=nc+1
          write (io1) h,k,l,f,sigf,e,sige
          n1=n1+1
          write (io2) h,k,l,f,sigf,e,sige
          n2=n2+1
        else
c
c Test acentric reflections by comparing point-group and Laue-group
c unique indicies.
c
          lh=h
          lk=k
          ll=l
          call equiv (laue,lh,lk,ll)
          hp=+lh
          kp=+lk
          lp=+ll
          hn=-lh
          kn=-lk
          ln=-ll
          call equiv (ptgp,hp,kp,lp)
          call equiv (ptgp,hn,kn,ln)
          if (h.eq.hp.and.k.eq.kp.and.l.eq.lp) then
c
c Write to +h+k+l reflection file.
c
            write (io1) hp,kp,lp,f,sigf,e,sige
            n1=n1+1
          else if (h.eq.hn.and.k.eq.kn.and.l.eq.ln) then
c
c Write to -h-k-l anti-reflection file, but with indices +h+k+l.
c
            write (io2) hp,kp,lp,f,sigf,e,sige
            n2=n2+1
          end if
        end if
      end do
 9    continue
      write (ilp,'(/1x,
     &''n1 = '',i6,'' F1(+h+k+l) SAS data''/1x,
     &''n2 = '',i6,'' F2(-h-k-l) SAS data'')') n1,n2
      if (nc.gt.0) write (ilp,'(/1x,
     &''nc = '',i6,'' centric reflections stored in''/1x,
     &''     '',6x,'' both the F1(+h+k+l) and F2(-h-k-l) subsets.''/1x,
     &''     '',6x,'' ----                ---'')') nc
      return
      end
c-----------------------------------------------------------------------
      subroutine csym (g,h,k,l,c)
c
c Table of centrosymmetric rows and zones in non-centrosymmetric
c reciprocal lattice point groups
c
c c = 0 for acentric distributions
c c = 1 for centric distributions
c
c Tricl. Monocl. Ortho.  Tetrag.    Trig.     Rhomb.    Cub.
c
c (1) 1  (3) 2   (6) 222  (9) 4     (17)  3   (33) R3   (38) 23
c (2) -1 (4) m   (7) mm2 (10) -4    (18)  -3  (34) R-3  (39) m3
c        (5) 2/m (8) mmm (11) 4/m
c                                   (19) 312  (35) R32  (40) 432
c                        (12) 422   (20) 321  (36) R3m  (41) -43m
c                        (13) 4mm   (21) 31m  (37) R-3m (42) m3m
c                        (14) -42m  (22) 3m1
c                        (15) -4m2  (23) -31m
c                        (16) 4/mmm (24) -3m1
c
c                                   Hexag.
c
c                                   (25) 6
c                                   (26) -6
c                                   (27) 6/m
c
c                                   (28) 622
c                                   (29) 6mm
c                                   (30) -6m2
c                                   (31) -62m
c                                   (32) 6/mmm
c
c A total of 42 cases is tabulated because there are ten cases of
c alternative axes for seven of the 32 crystallographic point groups:
c
c -42m  -4m2
c 3     R3
c -3    R-3
c 312   321   R32
c 31m   3m1   R3m
c -31m  -3m1  R-3m
c -62m  -6m2
c
      character g*5
      integer h,c
      c=0
c
c Triclinic
c
      if      (g.eq.'1    ') then
        c=0
      else if (g.eq.'-1   ') then
        c=1
c
c Monoclinic (b-axis unique)
c
      else if (g.eq.'2    ') then
        if (k.eq.0) c=1
      else if (g.eq.'m    ') then
        if (h.eq.0.and.l.eq.0) c=1
      else if (g.eq.'2/m  ') then
        c=1
c
c Orthorhombic
c
      else if (g.eq.'222  ') then
        if (h.eq.0.or.k.eq.0.or.l.eq.0) c=1
      else if (g.eq.'mm2  ') then
        if (l.eq.0) c=1
      else if (g.eq.'mmm  ') then
        c=1
c
c Tetragonal
c
      else if (g.eq.'4    ') then
        if (l.eq.0) c=1
      else if (g.eq.'-4   ') then
        if (l.eq.0) c=1
        if (h.eq.0.and.k.eq.0) c=1
      else if (g.eq.'4/m  ') then
        c=1
      else if (g.eq.'422  ') then
        if (l.eq.0.or.k.eq.0.or.h.eq.0) c=1
        if (abs(k).eq.abs(h)) c=1
      else if (g.eq.'4mm  ') then
        if (l.eq.0) c=1
      else if (g.eq.'-42m ') then
        if (l.eq.0.or.k.eq.0.or.h.eq.0) c=1
      else if (g.eq.'-4m2 ') then
        if (l.eq.0) c=1
        if (abs(k).eq.abs(h)) c=1
      else if (g.eq.'4/mmm') then
        c=1
c
c Trigonal, hexagonal axes
c
      else if (g.eq.'3    ') then
        c=0
      else if (g.eq.'-3   ') then
        c=1
      else if (g.eq.'312  ') then
        if (k.eq.h.or.k.eq.-2*h) c=1
      else if (g.eq.'321  ') then
        if (k.eq.0.or.h.eq.0) c=1
      else if (g.eq.'31m  ') then
        if ((k.eq.0.and.l.eq.0).or.(h.eq.0.and.l.eq.0)) c=1
      else if (g.eq.'3m1  ') then
        c=0
      else if (g.eq.'-31m ') then
        c=1
      else if (g.eq.'-3m1 ') then
        c=1
c
c Hexagonal
c
      else if (g.eq.'6    ') then
        if (l.eq.0) c=1
      else if (g.eq.'-6   ') then
        if (h.eq.0.and.k.eq.0) c=1
      else if (g.eq.'6/m  ') then
        c=1
      else if (g.eq.'622  ') then
        if (h.eq.0.or.k.eq.0.or.l.eq.0) c=1
        if (k.eq.h.or.k.eq.-2*h) c=1
      else if (g.eq.'6mm  ') then
        if (l.eq.0) c=1
      else if (g.eq.'-6m2 ') then
        if (k.eq.h.or.k.eq.-2*h) c=1
      else if (g.eq.'-62m ') then
        if (h.eq.0.or.k.eq.0) c=1
      else if (g.eq.'6/mmm') then
        c=1
c
c Trigonal, rhombohedral axes
c
      else if (g.eq.'R3   ') then
        c=0
      else if (g.eq.'R-3  ') then
        c=1
      else if (g.eq.'R32  ') then
        if (h.eq.k.or.k.eq.l.or.l.eq.h) c=1
      else if (g.eq.'R3m  ') then
        if (l.eq.0.and.k.eq.-h) c=1
        if (k.eq.0.and.l.eq.-h) c=1
        if (h.eq.0.and.l.eq.-k) c=1
      else if (g.eq.'R-3m ') then
        c=1
c
c Cubic
c
      else if (g.eq.'23   ') then
        if (h.eq.0.or.k.eq.0.or.l.eq.0) c=1
      else if (g.eq.'m3   ') then
        c=1
      else if (g.eq.'432  ') then
        if (h.eq.0.or.k.eq.0.or.l.eq.0) c=1
        if (abs(h).eq.abs(k).or.abs(k).eq.abs(l).or.abs(l).eq.abs(h))
     &  c=1
      else if (g.eq.'-43m ') then
        if (h.eq.0.or.k.eq.0.or.l.eq.0) c=1
      else if (g.eq.'m3m  ') then
        c=1
      else
        stop ' Invalid point group symbol in subroutine csym.'
      end if
      return
      end
c-----------------------------------------------------------------------
      subroutine datain (iunit1,iunit2,itype1,itype2,cutoff,
     & hmin,kmin,lmin,nk,nl,ginv,imax,data,jmin,jmax,type,ptgp,ilp)
c
c Store reflection data in memory.
c
      integer h,hmin,hmax,hkl
      character type*3,ptgp*5
      dimension data(10,0:imax),ginv(3,3)
      real*8 sumdf,sumf1,sumf2,sumde,sume,sum1sq,sum2sq
c
c Count data and find index limits.
c
      write (ilp,'(//1x,
     &''Data Selection for Local Scaling''/1x,
     &''--------------------------------''/1x,
     &''Minimum F/sigma(F) cutoff:''/1x,
     &''cutoff = '',f5.2)') cutoff
      hmin=+1e6
      hmax=-1e6
      kmin=+1e6
      kmax=-1e6
      lmin=+1e6
      lmax=-1e6
      m1=0
      n1=0
      iend=0
      do i=1,1e6
        call read1 (itype1,iunit1,iend,h,k,l,f,sigmaf,e,sigmae)
        if (iend.ne.0) go to 1
        m1=m1+1
        if (sigmaf.gt.0.and.f.ge.cutoff*sigmaf) then
          n1=n1+1
          hmin=min(hmin,h)
          hmax=max(hmax,h)
          kmin=min(kmin,k)
          kmax=max(kmax,k)
          lmin=min(lmin,l)
          lmax=max(lmax,l)
        end if
      end do
 1    continue
      m2=0
      n2=0
      iend=0
      do i=1,1e6
        call read1 (itype2,iunit2,iend,h,k,l,f,sigmaf,e,sigmae)
        if (iend.ne.0) go to 2
        m2=m2+1
        if (sigmaf.gt.0.and.f.ge.cutoff*sigmaf) then
          n2=n2+1
          hmin=min(hmin,h)
          hmax=max(hmax,h)
          kmin=min(kmin,k)
          kmax=max(kmax,k)
          lmin=min(lmin,l)
          lmax=max(lmax,l)
        end if
      end do
 2    continue
      nl=lmax-lmin+1
      nk=kmax-kmin+1
      write (ilp,'(/1x,
     &''m1 = '',i6,'' total F1 data''/1x,
     &''m2 = '',i6,'' total F2 data''//1x,
     &''n1 = '',i6,'' F1 data with F1 .ge. cutoff*sigma(F1)''/1x,
     &''n2 = '',i6,'' F2 data with F2 .ge. cutoff*sigma(F2)'')')
     & m1,m2,n1,n2
c
c Zero the data array.
c
      do i=0,imax
      do j=1,10
        data(j,i)=0
      end do
      end do
c
c Store data by packed index.
c
      jmin=imax
      jmax=0
      rewind iunit1
      do i=1,m1
        call read1 (itype1,iunit1,iend,h,k,l,f,sigmaf,e,sigmae)
        if (sigmaf.gt.0.and.f.ge.cutoff*sigmaf) then
          hkl=(h-hmin)*nk*nl+(k-kmin)*nl+(l-lmin)
          if (hkl.gt.imax) then
            write (ilp,1000) imax
          else
            jmin=min(jmin,hkl)
            jmax=max(jmax,hkl)
            data(1,hkl)=f
            data(2,hkl)=sigmaf
            data(5,hkl)=e
            data(6,hkl)=max(sigmae,e*sigmaf/f)
          end if
        end if
      end do
      rewind iunit2
      do i=1,m2
        call read1 (itype2,iunit2,iend,h,k,l,f,sigmaf,e,sigmae)
        if (sigmaf.gt.0.and.f.ge.cutoff*sigmaf) then
          hkl=(h-hmin)*nk*nl+(k-kmin)*nl+(l-lmin)
          if (hkl.gt.imax) then
            write (ilp,1000) imax
          else
            jmin=min(jmin,hkl)
            jmax=max(jmax,hkl)
            data(3,hkl)=f
            data(4,hkl)=sigmaf
            data(7,hkl)=e
            data(8,hkl)=max(sigmae,e*sigmaf/f)
          end if
        end if
      end do
 1000 format (/1x,
     &'Warning:  Packed index exceeds imax = ',i6/1x,
     &'Not all reflections can be stored in the data array.')
c
c Eliminate single, unpaired reflections from the data array.
c
      j1=0
      j2=0
      do j=jmin,jmax
        f1=data(1,j)
        f2=data(3,j)
        if (f1.le.0.xor.f2.le.0) then
          if (f1.gt.0) j1=j1+1
          if (f2.gt.0) j2=j2+1
          do i=1,8
            data(i,j)=0
          end do
        end if
      end do
      write (ilp,'(/1x,
     &''j1 = '',i6,'' single, unpaired F1 data eliminated from the da'',
     &''ta set''/1x,
     &''j2 = '',i6,'' single, unpaired F2 data eliminated from the da'',
     &''ta set'')') j1,j2
c
c Compile data statistics before scaling.
c
      n=0
      sumdf=0
      sumf1=0
      sumf2=0
      sumde=0
      sume=0
      sum1sq=0
      sum2sq=0
      smin=9
      smax=0
      nc=0
      do j=jmin,jmax
        f1=data(1,j)
        f2=data(3,j)
        if (f1.gt.0.and.f2.gt.0) then
c
c Accumulate agreement statistics.
c
          n=n+1
          sumdf=sumdf+abs(f1-f2)
          sumf1=sumf1+f1
          sumf2=sumf2+f2
          e1=data(5,j)
          e2=data(7,j)
          if (e1.gt.0.and.e2.gt.0) then
            sumde=sumde+abs(e1-e2)
            sume=sume+e1+e2
            sum1sq=sum1sq+e1**2
            sum2sq=sum2sq+e2**2
          end if
c
c Unpack indices, hkl = (h - hmin)*nk*nl + (k - kmin)*nl + (l - lmin).
c
          hkl=j
          h=hkl/(nk*nl)
          k=(hkl-h*nk*nl)/nl
          l=hkl-h*nk*nl-k*nl
          h=h+hmin
          k=k+kmin
          l=l+lmin
c
c Calculate sin(theta)/lambda.
c
          s=sinthl(h,k,l,ginv)
          smin=min(smin,s)
          smax=max(smax,s)
c
c Count centric reflections in SAS data sets.
c
          if (type.eq.'SAS') then
            call csym (ptgp,h,k,l,ic)
            if (ic.eq.1) nc=nc+1
          end if
c
c Store an initial local scale factor of unity.
c
          data(9,j)=1
        end if
      end do
      q=sumf1/sumf2
      rf=sumdf/(0.5*(sumf1+sumf2))
      if (sume.gt.0) then
        re=sumde/(0.5*sume)
        e1sq=sum1sq/n
        e2sq=sum2sq/n
      end if
      write (ilp,'(/1x,
     &''N  = '',i6,'' reflection pairs accepted for scaling'')') n
      if (nc.gt.0) write (ilp,'(/1x,
     &''nc = '',i6,'' centric reflections stored in''/1x,
     &''     '',6x,'' both the F1(+h+k+l) and F2(-h-k-l) subsets.''/1x,
     &''     '',6x,'' ----                ---'')') nc
      write (ilp,'(/1x,
     &''s = sin(theta)/lambda                  d = lambda/[2*sin(thet'',
     &''a)]''/1x,
     &''smin = '',f7.4,'' reciprocal Angstroms    dmax = '',f5.1,
     &''  Angstroms''/1x,
     &''smax = '',f7.4,''                         dmin = '',f6.2)')
     & smin,1/(2*smin),smax,1/(2*smax)
      write (ilp,'(/1x,
     &''Before scaling:''/1x,
     &''Q(F) = <F1>/<F2>                      = '',e10.3/1x,
     &''R(F) = <abs(F1 - F2)>/<0.5*(F1 + F2)> = '',e10.3/1x,
     &''R(E) = <abs(E1 - E2)>/<0.5*(E1 + E2)> = '',e10.3/1x,
     &''<E1**2> = '',f6.3/1x,''<E2**2> = '',f6.3)')
     & q,rf,re,e1sq,e2sq
      if (cutoff.gt.0) write (ilp,'(/1x,
     &''Values of <E**2> .gt. 1 are possible since minimum F/sigma(F)'',
     &'' = '',f5.2)') cutoff
      return
      end
c-----------------------------------------------------------------------
      subroutine decode (ptgp,laue,iptgp,ilaue)
c
c Point group and Laue group index numbers
c
c There are ten cases of alternative axes for seven of the 32
c crystallographic point groups, thus a total of 42 cases is tabulated.
c
c Similarly, there are three cases of alternative axes for two of the 11
c Laue groups, so 14 cases are tabulated.
c
      character ptgp*5,laue*5,symbol(42)*5
      data symbol
     &/'1    ','-1   ',
     & '2    ','m    ','2/m  ',
     & '222  ','mm2  ','mmm  ',
     & '4    ','-4   ','4/m  ',
     & '422  ','4mm  ','-42m ','-4m2 ','4/mmm',
     & '3    ','-3   ',
     & '312  ','321  ','31m  ','3m1  ','-31m ','-3m1 ',
     & '6    ','-6   ','6/m  ',
     & '622  ','6mm  ','-6m2 ','-62m ','6/mmm',
     & 'R3   ','R-3  ',
     & 'R32  ','R3m  ','R-3m ',
     & '23   ','m3   ',
     & '432  ','-43m ','m3m  '/
c
c Point group index numbers:
c
c Triclinic     1,  2;
c Monoclinic    3,  4,  5;
c Orthorhombic  6,  7,  8;
c Tetragonal    9, 10, 11,
c              12, 13, 14, 15, 16;
c Trigonal     17, 18,
c              19, 20, 21, 22, 23, 24;
c Hexagonal    25, 26, 27,
c              28, 29, 30, 31, 32;
c Rhombohedral 33, 34,
c              35, 36, 37;
c Cubic        38, 39,
c              40, 41, 42.
c 
      do i=1,42
        if (ptgp.eq.symbol(i)) go to 1
      end do
      stop ' Invalid input point group symbol.'
 1    iptgp=i
c
c Laue group
c
      if (i.ge. 1) ilaue=2
      if (i.ge. 3) ilaue=5
      if (i.ge. 6) ilaue=8
      if (i.ge. 9) ilaue=11
      if (i.ge.12) ilaue=16
      if (i.ge.17) ilaue=18
      if (i.eq.19.or.i.eq.21.or.i.eq.23) ilaue=23
      if (i.eq.20.or.i.eq.22.or.i.eq.24) ilaue=24
      if (i.ge.25) ilaue=27
      if (i.ge.28) ilaue=32
      if (i.ge.33) ilaue=34
      if (i.ge.35) ilaue=37
      if (i.ge.38) ilaue=39
      if (i.ge.40) ilaue=42
      laue=symbol(ilaue)
      return
      end
c-----------------------------------------------------------------------
      subroutine equiv (g,h,k,l)
c
c Equivalent, unique reflection indices for the crystallographic point
c groups
c
c Tricl. Monocl. Ortho.  Tetrag.    Trig.     Rhomb.    Cub.
c
c (1) 1  (3) 2   (6) 222  (9) 4     (17)  3   (33) R3   (38) 23
c (2) -1 (4) m   (7) mm2 (10) -4    (18)  -3  (34) R-3  (39) m3
c        (5) 2/m (8) mmm (11) 4/m
c                                   (19) 312  (35) R32  (40) 432
c                        (12) 422   (20) 321  (36) R3m  (41) -43m
c                        (13) 4mm   (21) 31m  (37) R-3m (42) m3m
c                        (14) -42m  (22) 3m1
c                        (15) -4m2  (23) -31m
c                        (16) 4/mmm (24) -3m1
c
c                                   Hexag.
c
c                                   (25) 6
c                                   (26) -6
c                                   (27) 6/m
c
c                                   (28) 622
c                                   (29) 6mm
c                                   (30) -6m2
c                                   (31) -62m
c                                   (32) 6/mmm
c
c See International Tables for Crystallography (1983).  Volume A, pp.
c 750-770.
c
c A total of 42 cases is tabulated because there are ten cases of
c alternative axes for seven of the 32 crystallographic point groups:
c
c -42m  -4m2
c 3     R3
c -3    R-3
c 312   321   R32
c 31m   3m1   R3m
c -31m  -3m1  R-3m
c -62m  -6m2
c
      character g*5
      integer h,x,y,z
c
c Triclinic
c
      if      (g.eq.'1    ') then
        return
      else if (g.eq.'-1   ') then
        if (l.gt.0) return
        h=-h
        k=-k
        l=-l
        if (l.gt.0) return
        if (k.gt.0) return
        h=-h
        k=-k
        if (k.gt.0) return
        if (h.ge.0) return
        h=-h
        return
c
c Monoclinic (b-axis unique)
c
      else if (g.eq.'2    ') then
 2      if (l.gt.0) return
        h=-h
        l=-l
        if (l.gt.0) return
        if (h.ge.0) return
        h=-h
        return
      else if (g.eq.'m    ') then
        k=abs(k)
        return
      else if (g.eq.'2/m  ') then
        k=abs(k)
        go to 2
c
c Orthorhombic
c
      else if (g.eq.'222  ') then
        if (h.eq.0.or.k.eq.0.or.l.eq.0) then
          h=abs(h)
          k=abs(k)
          l=abs(l)
          return
        end if
        if (h.ge.0.and.k.ge.0) return
        x=h
        y=k
        z=l
        h=-x
        k=-y
        l=z
        if (h.ge.0.and.k.ge.0) return
        h=-x
        k=y
        l=-z
        if (h.ge.0.and.k.ge.0) return
        h=x
        k=-y
        l=-z
        return
      else if (g.eq.'mm2  ') then
        h=abs(h)
        k=abs(k)
        return
      else if (g.eq.'mmm  ') then
        h=abs(h)
        k=abs(k)
        l=abs(l)
        return
c
c Tetragonal
c
      else if (g.eq.'4    ') then
 4      if (h.eq.0.and.k.eq.0) return
        if (h.gt.0.and.k.ge.0) return
        x=h
        y=k
        h=-y
        k=x
        if (h.gt.0.and.k.ge.0) return
        h=-x
        k=-y
        if (h.gt.0.and.k.ge.0) return
        h=y
        k=-x
        return
      else if (g.eq.'-4   ') then
        if (h.eq.0.and.k.eq.0) then
          l=abs(l)
          return
        end if
        if (h.gt.0.and.k.ge.0) return
        x=h
        y=k
        h=-x
        k=-y
        if (h.gt.0.and.k.ge.0) return
        h=y
        k=-x
        l=-l
        if (h.gt.0.and.k.ge.0) return
        h=-y
        k=x
        return
      else if (g.eq.'4/m  ') then
        l=abs(l)
        go to 4
      else if (g.eq.'422  ') then
        do n=0,1
          if (h.eq.0.or.k.eq.0.or.abs(h).eq.abs(k)) l=abs(l)
          if (h.ge.k.and.k.ge.0) return
          x=h
          y=k
          h=-x
          k=-y
          if (h.ge.k.and.k.ge.0) return
          h=-y
          k=x
          if (h.ge.k.and.k.ge.0) return
          h=y
          k=-x
          if (h.ge.k.and.k.ge.0) return
          h=-x
          k=y
          l=-l
        end do
      else if (g.eq.'4mm  ') then
 41     h=abs(h)
        k=abs(k)
        if (h.ge.k) return
        x=h
        h=k
        k=x
        return
      else if (g.eq.'-42m ') then
        do n=0,1
          if (h.eq.0.or.k.eq.0) l=abs(l)
          if (h.ge.k.and.k.ge.0) return
          x=h
          y=k
          h=-x
          k=-y
          if (h.ge.k.and.k.ge.0) return
          h=y
          k=-x
          l=-l
          if (h.ge.k.and.k.ge.0) return
          h=-y
          k=x
          if (h.ge.k.and.k.ge.0) return
          h=-x
          k=y
        end do
      else if (g.eq.'-4m2 ') then
        h=abs(h)
        k=abs(k)
        if (h.eq.k) l=abs(l)
        if (h.ge.k) return
        x=h
        h=k
        k=x
        l=-l
        return
      else if (g.eq.'4/mmm') then
        l=abs(l)
        go to 41
c
c Trigonal - hexagonal axes
c
      else if (g.eq.'3    ') then
 3      if (h.eq.0.and.k.eq.0) return
        if (h.lt.0.and.k.ge.0) return
        x=h
        y=k
        i=-h-k
        h=i
        k=x
        if (h.lt.0.and.k.ge.0) return
        h=y
        k=i
        return
      else if (g.eq.'-3   ') then
        if (h.eq.0.and.k.eq.0) then
          l=abs(l)
          return
        end if
        do n=0,1
          if (h.lt.0.and.k.ge.0.and.(l.gt.0.or.(l.eq.0.and.h.lt.-k)))
     &    return
          x=h
          y=k
          i=-h-k
          h=i
          k=x
          if (h.lt.0.and.k.ge.0.and.(l.gt.0.or.(l.eq.0.and.h.lt.-k)))
     &    return
          h=y
          k=i
          if (h.lt.0.and.k.ge.0.and.(l.gt.0.or.(l.eq.0.and.h.lt.-k)))
     &    return
          h=-x
          k=-y
          l=-l
        end do
      else if (g.eq.'312  ') then
        do n=0,1
          if (h.eq.0.or.k.eq.0.or.-h.eq.k) l=abs(l)
          if (h.ge.0.and.k.ge.0) return
          x=h
          y=k
          i=-h-k
          h=i
          k=x
          if (h.ge.0.and.k.ge.0) return
          h=y
          k=i
          if (h.ge.0.and.k.ge.0) return
          h=-y
          k=-x
          l=-l
        end do
      else if (g.eq.'321  ') then
        do n=0,1
          if (h.eq.k.or.-2*h.eq.k.or.-h.eq.2*k) l=abs(l)
          if (h.ge.0.and.h.ge.k.and.h.ge.-2*k) return
          x=h
          y=k
          i=-h-k
          h=i
          k=x
          if (h.ge.0.and.h.ge.k.and.h.ge.-2*k) return
          h=y
          k=i
          if (h.ge.0.and.h.ge.k.and.h.ge.-2*k) return
          h=y
          k=x
          l=-l
        end do
      else if (g.eq.'31m  ') then
        do n=0,1
          if (h.ge.0.and.h.ge.k.and.h.ge.-2*k) return
          x=h
          y=k
          i=-h-k
          h=i
          k=x
          if (h.ge.0.and.h.ge.k.and.h.ge.-2*k) return
          h=y
          k=i
          if (h.ge.0.and.h.ge.k.and.h.ge.-2*k) return
          h=y
          k=x
        end do
      else if (g.eq.'3m1  ') then
        do n=0,1
          if (h.ge.0.and.k.ge.0) return
          x=h
          y=k
          i=-h-k
          h=i
          k=x
          if (h.ge.0.and.k.ge.0) return
          h=y
          k=i
          if (h.ge.0.and.k.ge.0) return
          h=-y
          k=-x
        end do
      else if (g.eq.'-31m ') then
        do n=0,1
          if (h.eq.0.or.k.eq.0.or.-h.eq.k) l=abs(l)
          if (h.ge.0.and.h.ge.k.and.h.ge.-2*k.and.(l.gt.0.or.(l.eq.0
     &    .and.k.ge.0))) return
          x=h
          y=k
          i=-h-k
          h=i
          k=x
          if (h.ge.0.and.h.ge.k.and.h.ge.-2*k.and.(l.gt.0.or.(l.eq.0
     &    .and.k.ge.0))) return
          h=y
          k=i
          if (h.ge.0.and.h.ge.k.and.h.ge.-2*k.and.(l.gt.0.or.(l.eq.0
     &    .and.k.ge.0))) return
          h=y
          k=x
          if (h.ge.0.and.h.ge.k.and.h.ge.-2*k.and.(l.gt.0.or.(l.eq.0
     &    .and.k.ge.0))) return
          h=x
          k=i
          if (h.ge.0.and.h.ge.k.and.h.ge.-2*k.and.(l.gt.0.or.(l.eq.0
     &    .and.k.ge.0))) return
          h=i
          k=y
          if (h.ge.0.and.h.ge.k.and.h.ge.-2*k.and.(l.gt.0.or.(l.eq.0
     &    .and.k.ge.0))) return
          h=-x
          k=-y
          l=-l
        end do
      else if (g.eq.'-3m1 ') then
        do n=0,1
          if (h.eq.k.or.-2*h.eq.k.or.-h.eq.2*k) l=abs(l)
          if (h.ge.0.and.k.ge.0.and.(l.gt.0.or.(l.eq.0.and.h.ge.k)))
     &    return
          x=h
          y=k
          i=-h-k
          h=i
          k=x
          if (h.ge.0.and.k.ge.0.and.(l.gt.0.or.(l.eq.0.and.h.ge.k)))
     &    return
          h=y
          k=i
          if (h.ge.0.and.k.ge.0.and.(l.gt.0.or.(l.eq.0.and.h.ge.k)))
     &    return
          h=y
          k=x
          l=-l
          if (h.ge.0.and.k.ge.0.and.(l.gt.0.or.(l.eq.0.and.h.ge.k)))
     &    return
          h=x
          k=i
          if (h.ge.0.and.k.ge.0.and.(l.gt.0.or.(l.eq.0.and.h.ge.k)))
     &    return
          h=i
          k=y
          if (h.ge.0.and.k.ge.0.and.(l.gt.0.or.(l.eq.0.and.h.ge.k)))
     &    return
          h=-x
          k=-y
        end do
c
c Hexagonal
c
      else if (g.eq.'6    ') then
 6      if (h.eq.0.and.k.eq.0) return
        if (h.gt.0.and.k.ge.0) return
        x=h
        y=k
        i=-h-k
        h=i
        k=x
        if (h.gt.0.and.k.ge.0) return
        h=y
        k=i
        if (h.gt.0.and.k.ge.0) return
        h=-x
        k=-y
        if (h.gt.0.and.k.ge.0) return
        h=-i
        k=-x
        if (h.gt.0.and.k.ge.0) return
        h=-y
        k=-i
        return
      else if (g.eq.'-6   ') then
        l=abs(l)
        go to 3
      else if (g.eq.'6/m  ') then
        l=abs(l)
        go to 6
      else if (g.eq.'622  ') then
        do n=0,1
          if (h.eq.0.or.k.eq.0.or.abs(h).eq.abs(k).or.abs(2*h).eq.abs(k)
     &    .or.abs(h).eq.abs(2*k)) l=abs(l)
          if (h.ge.k.and.k.ge.0) return
          x=h
          y=k
          i=-h-k
          h=i
          k=x
          if (h.ge.k.and.k.ge.0) return
          h=y
          k=i
          if (h.ge.k.and.k.ge.0) return
          h=-x
          k=-y
          if (h.ge.k.and.k.ge.0) return
          h=-i
          k=-x
          if (h.ge.k.and.k.ge.0) return
          h=-y
          k=-i
          if (h.ge.k.and.k.ge.0) return
          h=y
          k=x
          l=-l
        end do
      else if (g.eq.'6mm  ') then
 61     do n=0,1
          if (h.ge.k.and.k.ge.0) return
          x=h
          y=k
          i=-h-k
          h=i
          k=x
          if (h.ge.k.and.k.ge.0) return
          h=y
          k=i
          if (h.ge.k.and.k.ge.0) return
          h=-x
          k=-y
          if (h.ge.k.and.k.ge.0) return
          h=-i
          k=-x
          if (h.ge.k.and.k.ge.0) return
          h=-y
          k=-i
          if (h.ge.k.and.k.ge.0) return
          h=y
          k=x
        end do
      else if (g.eq.'-6m2 ') then
        l=abs(l)
        if (h.eq.0.and.k.eq.0) return
        do n=0,1
          if (h.gt.0.and.k.ge.0) return
          x=h
          y=k
          i=-h-k
          h=i
          k=x
          if (h.gt.0.and.k.ge.0) return
          h=y
          k=i
          if (h.gt.0.and.k.ge.0) return
          h=-y
          k=-x
        end do
      else if (g.eq.'-62m ') then
        l=abs(l)
        if (h.eq.0.and.k.eq.0) return
        do n=0,1
          if (h.gt.0.and.h.ge.k.and.h.gt.-2*k) return
          x=h
          y=k
          i=-h-k
          h=i
          k=x
          if (h.gt.0.and.h.ge.k.and.h.gt.-2*k) return
          h=y
          k=i
          if (h.gt.0.and.h.ge.k.and.h.gt.-2*k) return
          h=y
          k=x
        end do
      else if (g.eq.'6/mmm') then
        l=abs(l)
        go to 61
c
c Trigonal - rhombohedral axes
c
      else if (g.eq.'R3   ') then
 39     stop ' Re-index on hexagonal axes.'
      else if (g.eq.'R-3  ') then
        go to 39
      else if (g.eq.'R32  ') then
        go to 39
      else if (g.eq.'R3m  ') then
        go to 39
      else if (g.eq.'R-3m ') then
        go to 39
c
c Cubic
c
      else if (g.eq.'23   ') then
        if (h.eq.k.and.k.eq.abs(l)) return
        do n=0,2
          if (h.gt.abs(l).and.k.ge.abs(l)) return
          x=h
          y=k
          z=l
          h=-x
          k=-y
          l=z
          if (h.gt.abs(l).and.k.ge.abs(l)) return
          h=-x
          k=y
          l=-z
          if (h.gt.abs(l).and.k.ge.abs(l)) return
          h=x
          k=-y
          l=-z
          if (h.gt.abs(l).and.k.ge.abs(l)) return
          h=z
          k=x
          l=y
        end do
      else if (g.eq.'m3   ') then
        h=abs(h)
        k=abs(k)
        l=abs(l)
        if (h.eq.k.and.k.eq.l) return
        do n=0,2
          if (h.gt.l.and.k.ge.l) return
          x=h
          y=k
          z=l
          h=z
          k=x
          l=y
          if (h.gt.l.and.k.ge.l) return
          h=y
          k=z
          l=x
        end do
      else if (g.eq.'432  ') then
        if (abs(h).eq.abs(k).and.abs(k).eq.abs(l)) then
          h=abs(h)
          k=abs(k)
          l=abs(l)
          return
        end if
        do n=0,2
          if (h.gt.l.and.k.ge.l.and.l.ge.0) return
          x=h
          y=k
          z=l
          h=z
          k=x
          l=y
          if (h.gt.l.and.k.ge.l.and.l.ge.0) return
          h=y
          k=z
          l=x
          if (h.gt.l.and.k.ge.l.and.l.ge.0) return
          h=-x
          k=-y
          l=z
          if (h.gt.l.and.k.ge.l.and.l.ge.0) return
          h=z
          k=-x
          l=-y
          if (h.gt.l.and.k.ge.l.and.l.ge.0) return
          h=-y
          k=z
          l=-x
          if (h.gt.l.and.k.ge.l.and.l.ge.0) return
          h=-x
          k=y
          l=-z
          if (h.gt.l.and.k.ge.l.and.l.ge.0) return
          h=-z
          k=-x
          l=y
          if (h.gt.l.and.k.ge.l.and.l.ge.0) return
          h=y
          k=-z
          l=-x
          if (h.gt.l.and.k.ge.l.and.l.ge.0) return
          h=x
          k=-y
          l=-z
          if (h.gt.l.and.k.ge.l.and.l.ge.0) return
          h=-z
          k=x
          l=-y
          if (h.gt.l.and.k.ge.l.and.l.ge.0) return
          h=-y
          k=-z
          l=x
          if (h.gt.l.and.k.ge.l.and.l.ge.0) return
          h=y
          k=x
          l=-z
        end do
      else if (g.eq.'-43m ') then
        do n=0,3
          if (h.ge.k.and.k.ge.abs(l)) return
          x=h
          y=k
          z=l
          h=z
          k=x
          l=y
          if (h.ge.k.and.k.ge.abs(l)) return
          h=y
          k=z
          l=x
          if (h.ge.k.and.k.ge.abs(l)) return
          h=-x
          k=-y
          l=z
          if (h.ge.k.and.k.ge.abs(l)) return
          h=z
          k=-x
          l=-y
          if (h.ge.k.and.k.ge.abs(l)) return
          h=-y
          k=z
          l=-x
          if (h.ge.k.and.k.ge.abs(l)) return
          h=-x
          k=y
          l=-z
          if (h.ge.k.and.k.ge.abs(l)) return
          h=-z
          k=-x
          l=y
          if (h.ge.k.and.k.ge.abs(l)) return
          h=y
          k=-z
          l=-x
          if (h.ge.k.and.k.ge.abs(l)) return
          h=x
          k=-y
          l=-z
          if (h.ge.k.and.k.ge.abs(l)) return
          h=-z
          k=x
          l=-y
          if (h.ge.k.and.k.ge.abs(l)) return
          h=-y
          k=-z
          l=x
          if (h.ge.k.and.k.ge.abs(l)) return
          h=y
          k=x
          l=z
        end do
      else if (g.eq.'m3m  ') then
        h=abs(h)
        k=abs(k)
        l=abs(l)
        if (h.ge.k.and.k.ge.l) return
        x=h
        y=k
        z=l
        h=z
        k=x
        l=y
        if (h.ge.k.and.k.ge.l) return
        h=y
        k=z
        l=x
        if (h.ge.k.and.k.ge.l) return
        h=z
        k=y
        l=x
        if (h.ge.k.and.k.ge.l) return
        h=x
        k=z
        l=y
        if (h.ge.k.and.k.ge.l) return
        h=y
        k=x
        l=z
        return
      else
        stop ' Invalid point group symbol in subroutine equiv.'
      end if
      end
c-----------------------------------------------------------------------
      subroutine esym (g,h,k,l,e,m)
c
c Reciprocal lattice point degeneracies, e = epsilon(hkl), and
c multiplicities, m = m(hkl), for the crystallographic point groups
c
c As tabulated by:
c
c Hitoshi Iwasaki and Tetsuzo Ito (1977).  Values of epsilon for
c obtaining normalized structure factors.  Acta Cryst. A33, 227-229.
c
c Tricl. Monocl. Ortho.  Tetrag.    Trig.     Rhomb.    Cub.
c
c (1) 1  (3) 2   (6) 222 (9)  4     (17) 3    (33) r3   (38) 23
c (2) -1 (4) m   (7) mm2 (10) -4    (18) -3   (34) r-3  (39) m3
c        (5) 2/m (8) mmm (11) 4/m   
c                                   (19) 312  (35) r32  (40) 432
c                        (12) 422   (20) 321  (36) r3m  (41) -43m
c                        (13) 4mm   (21) 31m  (37) r-3m (42) m3m
c                        (14) -42m  (22) 3m1
c                        (15) -4m2  (23) -31m
c                        (16) 4/mmm (24) -3m1
c 
c                                    hexag.
c
c                                   (25) 6
c                                   (26) -6
c                                   (27) 6/m
c
c                                   (28) 622
c                                   (29) 6mm
c                                   (30) -6m2
c                                   (31) -62m
c                                   (32) 6/mmm
c
c A total of 42 cases is tabulated because there are ten cases of
c alternative axes for seven of the 32 crystallographic point groups:
c
c -42m  -4m2
c 3     r3
c -3    r-3
c 312   321   r32
c 31m   3m1   r3m
c -31m  -3m1  r-3m
c -62m  -6m2
c
      character g*5
      integer h
      e=1
      if      (g.eq.'1    ') then
        m=1
      else if (g.eq.'-1   ') then      
        m=2
      else if (g.eq.'2    ') then
        m=2
        if (h.eq.0.and.l.eq.0) e=2
      else if (g.eq.'m    ') then
        m=2
        if (k.eq.0) e=2
      else if (g.eq.'2/m  ') then
        m=4
        if (k.eq.0) e=2
        if (h.eq.0.and.l.eq.0) e=2
      else if (g.eq.'222  ') then
        m=4
        if (k.eq.0.and.l.eq.0) e=2
        if (h.eq.0.and.l.eq.0) e=2
        if (h.eq.0.and.k.eq.0) e=2
      else if (g.eq.'mm2  ') then
        m=4
        if (h.eq.0) e=2
        if (k.eq.0) e=2
        if (h.eq.0.and.k.eq.0) e=4
      else if (g.eq.'mmm  ') then
        m=8
        if (h.eq.0) e=2
        if (k.eq.0) e=2
        if (l.eq.0) e=2
        if (k.eq.0.and.l.eq.0) e=4
        if (h.eq.0.and.l.eq.0) e=4
        if (h.eq.0.and.k.eq.0) e=4
      else if (g.eq.'4    ') then
        m=4
        if (h.eq.0.and.k.eq.0) e=4
      else if (g.eq.'-4   ') then
        m=4
        if (h.eq.0.and.k.eq.0) e=2
      else if (g.eq.'4/m  ') then
        m=8
        if (l.eq.0) e=2
        if (h.eq.0.and.k.eq.0) e=4
      else if (g.eq.'422  ') then
        m=8
        if (l.eq.0.and.k.eq.h.or.k.eq.-h) e=2
        if (l.eq.0.and.k.eq.0.or.h.eq.0) e=2
        if (h.eq.0.and.k.eq.0) e=4
      else if (g.eq.'4mm  ') then
        m=8
        if (k.eq.0.or.h.eq.0) e=2
        if (k.eq.h.or.k.eq.-h) e=2
        if (h.eq.0.and.k.eq.0) e=8
      else if (g.eq.'-42m ') then
        m=8
        if (k.eq.h.or.k.eq.-h) e=2
        if (l.eq.0.and.(k.eq.0.or.h.eq.0)) e=2
        if (h.eq.0.and.k.eq.0) e=4
      else if (g.eq.'-4m2 ') then
        m=8
        if (k.eq.0.or.h.eq.0) e=2
        if (l.eq.0.and.(k.eq.h.or.k.eq.-h)) e=2
        if (h.eq.0.and.k.eq.0) e=4
      else if (g.eq.'4/mmm') then
        m=16
        if (k.eq.0.or.h.eq.0.or.l.eq.0) e=2
        if (k.eq.h.or.k.eq.-h) e=2
        if (l.eq.0.and.(k.eq.h.or.k.eq.-h)) e=4
        if (l.eq.0.and.(k.eq.0.or.h.eq.0)) e=4
        if (h.eq.0.and.k.eq.0) e=8
      else if (g.eq.'3    ') then
        m=3
        if (h.eq.0.and.k.eq.0) e=3
      else if (g.eq.'-3   ') then
        m=6
        if (h.eq.0.and.k.eq.0) e=3
      else if (g.eq.'312  ') then
        m=6
        if (l.eq.0.and.k.eq.0.or.h.eq.0.or.k.eq.-h) e=2
        if (h.eq.0.and.k.eq.0) e=3
      else if (g.eq.'321  ') then
        m=6
        if (l.eq.0.and.(k.eq.h.or.k.eq.-2*h.or.h.eq.-2*k)) e=2
        if (h.eq.0.and.k.eq.0) e=3
      else if (g.eq.'31m  ') then
        m=6
        if (k.eq.h.or.k.eq.-2*h.or.h.eq.-2*k) e=2
        if (h.eq.0.and.k.eq.0) e=6
      else if (g.eq.'3m1  ') then
        m=6
        if (k.eq.0.or.h.eq.0.or.k.eq.-h) e=2
        if (h.eq.0.and.k.eq.0) e=6
      else if (g.eq.'-31m ') then
        m=12
        if (k.eq.h.or.k.eq.-2*h.or.h.eq.-2*k) e=2
        if (l.eq.0.and.(k.eq.0.or.h.eq.0.or.k.eq.-h)) e=2
        if (h.eq.0.and.k.eq.0) e=6
      else if (g.eq.'-3m1 ') then
        m=12
        if (l.eq.0.and.(k.eq.h.or.k.eq.-2*h.or.h.eq.-2*k)) e=2
        if (k.eq.0.or.h.eq.0.or.k.eq.-h) e=2
        if (h.eq.0.and.k.eq.0) e=6
      else if (g.eq.'6    ') then
        m=6
        if (h.eq.0.and.k.eq.0) e=6
      else if (g.eq.'-6   ') then
        m=6
        if (l.eq.0) e=2
        if (h.eq.0.and.k.eq.0) e=3
      else if (g.eq.'6/m  ') then
        m=12
        if (l.eq.0) e=2
        if (h.eq.0.and.k.eq.0) e=6
      else if (g.eq.'622  ') then
        m=12
        if (l.eq.0.and.(k.eq.h.or.k.eq.-2*h.or.h.eq.-2*k)) e=2
        if (l.eq.0.and.(k.eq.0.or.h.eq.0.or.k.eq.-h)) e=2
        if (h.eq.0.and.k.eq.0) e=6
      else if (g.eq.'6mm  ') then
        m=12
        if (k.eq.h.or.k.eq.-2*h.or.h.eq.-2*k) e=2
        if (k.eq.0.or.h.eq.0.or.k.eq.-h) e=2
        if (h.eq.0.and.k.eq.0) e=12
      else if (g.eq.'-6m2 ') then
        m=12
        if (l.eq.0) e=2
        if (k.eq.0.or.h.eq.0.or.k.eq.-h) e=2
        if (l.eq.0.and.(k.eq.0.or.h.eq.0.or.k.eq.-h)) e=4
        if (h.eq.0.and.k.eq.0) e=6
      else if (g.eq.'-62m ') then
        m=12
        if (l.eq.0) e=2
        if (k.eq.h.or.k.eq.-2*h.or.h.eq.-2*k) e=2
        if (l.eq.0.and.(k.eq.h.or.k.eq.-2*h.or.h.eq.-2*k)) e=4
        if (h.eq.0.and.k.eq.0) e=6
      else if (g.eq.'6/mmm') then
        m=24
        if (l.eq.0) e=2
        if (k.eq.h.or.k.eq.-2*h.or.h.eq.-2*k) e=2
        if (l.eq.0.and.(k.eq.h.or.k.eq.-2*h.or.h.eq.-2*k)) e=4
        if (k.eq.0.or.h.eq.0.or.k.eq.-h) e=2
        if (l.eq.0.and.(k.eq.0.or.h.eq.0.or.k.eq.-h)) e=4
        if (h.eq.0.and.k.eq.0) e=12
      else if (g.eq.'r3   ') then
        m=3
        if (k.eq.h.and.l.eq.h) e=3
      else if (g.eq.'r-3  ') then
        m=6
        if (k.eq.h.and.l.eq.h) e=3
      else if (g.eq.'r32  ') then
        m=6
        if (l.eq.0.and.k.eq.-h) e=2
        if (k.eq.0.and.l.eq.-h) e=2
        if (h.eq.0.and.k.eq.-l) e=2
        if (k.eq.h.and.l.eq.h) e=3
      else if (g.eq.'r3m  ') then
        m=6
        if (k.eq.h.or.l.eq.h.or.l.eq.k) e=2
        if (abs(k).eq.abs(h).and.abs(l).eq.abs(h)) e=2
        if (k.eq.h.and.l.eq.h) e=6
      else if (g.eq.'r-3m ') then
        m=12
        if (k.eq.h.or.l.eq.h.or.k.eq.l) e=2
        if (abs(k).eq.abs(h).and.abs(l).eq.abs(h)) e=2
        if (k.eq.-h.and.l.eq.0) e=2
        if (l.eq.-h.and.k.eq.0) e=2
        if (l.eq.-k.and.h.eq.0) e=2
        if (k.eq.h.and.l.eq.h) e=6
      else if (g.eq.'23   ') then
        m=12
        if (abs(k).eq.abs(h).and.abs(l).eq.abs(h)) e=3
        if (k.eq.0.and.l.eq.0) e=2
        if (h.eq.0.and.l.eq.0) e=2
        if (h.eq.0.and.k.eq.0) e=2
      else if (g.eq.'m3   ') then
        m=24
        if (abs(k).eq.abs(h).and.abs(l).eq.abs(h)) e=3
        if (h.eq.0.or.k.eq.0.or.l.eq.0) e=2
        if (k.eq.0.and.l.eq.0) e=4
        if (h.eq.0.and.l.eq.0) e=4
        if (h.eq.0.and.k.eq.0) e=4
      else if (g.eq.'432  ') then
        m=24
        if (abs(k).eq.abs(h).and.abs(l).eq.abs(h)) e=3
        if (l.eq.0.and.abs(k).eq.abs(h)) e=2
        if (k.eq.0.and.abs(l).eq.abs(h)) e=2
        if (h.eq.0.and.abs(k).eq.abs(l)) e=2
        if (k.eq.0.and.l.eq.0) e=4
        if (h.eq.0.and.l.eq.0) e=4
        if (h.eq.0.and.k.eq.0) e=4
      else if (g.eq.'-43m ') then
        m=24
        if (abs(k).eq.abs(h).or.abs(l).eq.abs(h).or.abs(k).eq.abs(l))
     &  e=2
        if (abs(k).eq.abs(h).and.abs(l).eq.abs(h)) e=6
        if (k.eq.0.and.l.eq.0) e=4
        if (h.eq.0.and.l.eq.0) e=4
        if (h.eq.0.and.k.eq.0) e=4
      else if (g.eq.'m3m  ') then
        m=48
        if (abs(k).eq.abs(h).or.abs(l).eq.abs(h).or.abs(k).eq.abs(l))
     &  e=2
        if (abs(k).eq.abs(h).and.abs(l).eq.abs(h)) e=6
        if (h.eq.0.or.k.eq.0.or.l.eq.0) e=2
        if (l.eq.0.and.abs(k).eq.abs(h)) e=4
        if (k.eq.0.and.abs(l).eq.abs(h)) e=4
        if (h.eq.0.and.abs(k).eq.abs(l)) e=4
        if (k.eq.0.and.l.eq.0) e=8
        if (h.eq.0.and.l.eq.0) e=8
        if (h.eq.0.and.k.eq.0) e=8
      end if
      m=m/e
      return
      end
c-----------------------------------------------------------------------
      subroutine fsqtof (fsq,sigfsq,f,sigf)
      if (fsq.lt.0) then
        f=0
      else
        f=sqrt(fsq)
      end if
      if (fsq.le.sigfsq) then
        sigf=0.5*sqrt(sigfsq)
      else
        sigf=0.5*sigfsq/f
      end if
      return
      end
c-----------------------------------------------------------------------
      subroutine metric (a,gb)
c
c Calculate reciprocal lattice parameters and direct and reciprocal
c lattice metric tensors from direct lattice parameters.
c
      dimension a(6),b(6),ga(3,3),gb(3,3)
      pi=acos(-1.0)
      deg=180/pi
      rad=pi/180
      ca=cos(a(4)*rad)
      cb=cos(a(5)*rad)
      cg=cos(a(6)*rad)
      sa=sin(a(4)*rad)
      sb=sin(a(5)*rad)
      sg=sin(a(6)*rad)
      v=a(1)*a(2)*a(3)*sqrt(1-ca**2-cb**2-cg**2+2*ca*cb*cg)
      b(1)=a(2)*a(3)*sa/v
      b(2)=a(1)*a(3)*sb/v
      b(3)=a(1)*a(2)*sg/v
      b(4)=(cb*cg-ca)/(sb*sg)
      b(5)=(ca*cg-cb)/(sa*sg)
      b(6)=(ca*cb-cg)/(sa*sb)
      ga(1,1)=a(1)*a(1)
      ga(1,2)=a(1)*a(2)*cg
      ga(1,3)=a(1)*a(3)*cb
      ga(2,1)=ga(1,2)
      ga(2,2)=a(2)*a(2)
      ga(2,3)=a(2)*a(3)*ca
      ga(3,1)=ga(1,3)
      ga(3,2)=ga(2,3)
      ga(3,3)=a(3)*a(3)
      gb(1,1)=b(1)*b(1)
      gb(1,2)=b(1)*b(2)*b(6)
      gb(1,3)=b(1)*b(3)*b(5)
      gb(2,1)=gb(1,2)
      gb(2,2)=b(2)*b(2)
      gb(2,3)=b(2)*b(3)*b(4)
      gb(3,1)=gb(1,3)
      gb(3,2)=gb(2,3)
      gb(3,3)=b(3)*b(3)
      b(4)=acos(b(4))*deg
      b(5)=acos(b(5))*deg
      b(6)=acos(b(6))*deg
      return
      end
c-----------------------------------------------------------------------
      subroutine minv3 (a,b,d)
c
c Inverse of a 3 x 3 matrix
c
      dimension a(3,3),b(3,3)
c
c Determinant
c
      d=a(1,1)*a(2,2)*a(3,3)+a(1,2)*a(2,3)*a(3,1)+a(1,3)*a(2,1)*a(3,2)
     &  -a(3,1)*a(2,2)*a(1,3)-a(3,2)*a(2,3)*a(1,1)-a(3,3)*a(2,1)*a(1,2)
      if (d.eq.0) return
c
c Matrix of minors
c
      b(1,1)=a(2,2)*a(3,3)-a(3,2)*a(2,3)
      b(1,2)=a(2,1)*a(3,3)-a(3,1)*a(2,3)
      b(1,3)=a(2,1)*a(3,2)-a(3,1)*a(2,2)
      b(2,1)=a(1,2)*a(3,3)-a(3,2)*a(1,3)
      b(2,2)=a(1,1)*a(3,3)-a(3,1)*a(1,3)
      b(2,3)=a(1,1)*a(3,2)-a(3,1)*a(1,2)
      b(3,1)=a(1,2)*a(2,3)-a(2,2)*a(1,3)
      b(3,2)=a(1,1)*a(2,3)-a(2,1)*a(1,3)
      b(3,3)=a(1,1)*a(2,2)-a(2,1)*a(1,2)
c
c Matrix of co-factors (signed minors)
c
      do i=1,3
      do j=1,3
        b(i,j)=(-1)**(i+j)*b(i,j)
      end do
      end do
c
c Adjoint matrix (transposed matrix of co-factors)
c
      do i=1,3-1
      do j=i+1,3
        t=b(i,j)
        b(i,j)=b(j,i)
        b(j,i)=t
      end do
      end do
c
c Inverse matrix
c
      do i=1,3
      do j=1,3
        b(i,j)=b(i,j)/d
      end do
      end do
      return
      end
c-----------------------------------------------------------------------
      subroutine mv (a,b,c)
c
c Square matrix-column vector multiplication
c
c a(i,j)*b(j) = c(i)
c
      dimension a(3,3),b(3),c(3)
      do i=1,3
        c(i)=0
      do j=1,3
        c(i)=c(i)+a(i,j)*b(j)
      end do
      end do
      return
      end
c-----------------------------------------------------------------------
      subroutine output (ifile,jfile,ilp,type,ptgp,ndata,nc,zlimit,
     & atime,adate,title)
c
c Compile Delta(F) distribution statistics and histogram,
c optionally trim data from the distribution tails, and
c write the hkl-sorted and Delta/sigma(Delta)-sorted output files.
c
      parameter (imax=99999,jmax=25)
      character atime*8,adate*9,title*72,type*3,ptgp*5
      dimension z(imax),index(imax),xh(jmax),yh(jmax)
c
c Find  median(F1 - F2).
c
      rewind ifile
 1000 format (1x,3i5,2(e12.5,e9.2),2(2f7.3))
      do i=1,ndata
        read (ifile,1000) ih,ik,il,f1,s1,f2
        z(i)=f1-f2
      end do
      call sort (ndata,z,index)
      zmedian=z(index(nint(0.5*ndata)))
c
c Compile  Delta(F) = F1 - F2 - median(F1 - F2)  histogram.
c
      zmean=0
      do i=1,ndata
        d=z(i)-zmedian
        zmean=zmean+d
        z(i)=d
      end do
      zmin=z(index(1))
      zmax=z(index(ndata))
      zmean=zmean/ndata
      do j=1,jmax
        xh(j)=0
        yh(j)=0
      end do
      j=0
      zj=zmin
      d=(zmax-zmin)/jmax
      do i=1,ndata
        zi=z(index(i))
        do while (zi.ge.zj)
          j=j+1
          zj=zj+d
        end do
        if (j.le.jmax) then
          xh(j)=xh(j)+zi
          yh(j)=yh(j)+1
        end if
      end do
      xh(jmax)=xh(jmax)+zi
      yh(jmax)=yh(jmax)+1
c
c In SAS cases, reduce the count in the histogram bin that contains the
c centric-pair zero-differences.
c
      if (type.eq.'SAS') then
        j=0
        zj=zmin
        do while (zj.lt.0)
          j=j+1
          zj=zj+d
        end do
        yh(j)=yh(j)-nc
      end if
      yhmax=0
      do j=1,jmax
        if (yh(j).gt.0) then
          xh(j)=xh(j)/yh(j)
          yhmax=max(yhmax,yh(j))
        end if
      end do
c
c Print statistics summary.
c        
      write (ilp,'(//1x,
     &''Delta(F) = F1 - F2 - median(F1 - F2) Distribution Statistics''/
     &1x,
     &''------------------------------------------------------------''/
     &1x,
     &''Before trimming distribution tails from the data set:''/1x,
     &''------''/1x,
     &''  N                    = '',i10,''  Delta(F) values''/1x,
     &''  median (F1 - F2)     = '',e10.3/1x,
     &''  minimum Delta(F)     = '',e10.3/1x,
     &''  maximum              = '',e10.3/1x,
     &''  mean                 = '',e10.3)')
     & ndata,zmedian,zmin,zmax,zmean
c
c Store  median(F1 - F2).
c
      f0=zmedian
c
c Find the median and mean of the nonzero abs[Delta(F)] values.
c
      nz=0
      zmean=0
      do i=1,ndata
        d=abs(z(i))
        if (d.eq.0) nz=nz+1
        zmean=zmean+d
        z(i)=d
      end do
      call sort (ndata,z,index)
      zmedian=z(index(nint(0.5*(ndata+nz))))
      zmean=zmean/(ndata-nz)
      write (ilp,'(1x,
     &''  median abs[Delta(F)] = '',e10.3/1x,
     &''  mean                 = '',e10.3)')
     & zmedian,zmean
c
c Print Delta(F) histogram.
c
      write (ilp,'(/1x,
     &''Delta(F) distribution before trimming distribution tails''/1x,
     &''--------------------------------------------------------''/1x,
     &''     n  <Delta(F)>''/1x,''     -  ----------'')')
      do j=1,jmax
        n=nint(50*yh(j)/yhmax)
        write (ilp,'(1x,i6,e12.3,1x,50a1)') nint(yh(j)),xh(j),
     &  ('X',i=1,n)
      end do
c
c Write the hkl-sorted, locally scaled output file "data.locscl" and
c optionally reject data pairs with abnormal Delta(F) values in the
c extreme tails of the Delta(F) distribution.
c
      sigma=1.25*zmedian
      dlimit=zlimit*sigma
      nrej=0
      nrejn=0
      nrejp=0
      rewind ifile
      open (unit=jfile,file='data.locscl',status='new')
      write (jfile,'(1x,
     &''Program locscl.  '',a,'', '',a,''.  '',a//1x,
     &''    h    k    l          F1  sigmaF1          F2  sigmaF2'',
     &''    E1 sigmaE1    E2 sigmaE2''/1x,
     &''    -    -    -          --  -------          --  -------'',
     &''    -- -------    -- -------'')') atime,adate,title
      j=0
      zmean=0
      do i=1,ndata
        read (ifile,1000) ih,ik,il,f1,s1,f2,s2,e1,t1,e2,t2
        d=f1-f2-f0
        if (dlimit.gt.0.and.abs(d).gt.dlimit) then
          nrej=nrej+1
          if (d.lt.0) nrejn=nrejn+1
          if (d.gt.0) nrejp=nrejp+1
          if (nrej.eq.1) write (ilp,'(/1x,
     &''Reflection pairs trimmed from the distribution tails:''/1x,
     &''-----------------------------------------------------''/1x,
     &''    h    k    l          F1  sigmaF1          F2  sigmaF2'',
     &''    E1 sigmaE1    E2 sigmaE2 (F1-F2)/sqrt(sigF1**2+sigF2**2)''/
     &1x,
     &''    -    -    -          --  -------          --  -------'',
     &''    -- -------    -- ------- ------------------------------'')')
          if (nrej.le.100) write (ilp,
     &    '(1x,3i5,2(e12.5,e9.2),2(2f7.3),e10.2)')
     &    ih,ik,il,f1,s1,f2,s2,e1,t1,e2,t2,(f1-f2)/sqrt(s1**2+s2**2)
        else
          write (jfile,1000) ih,ik,il,f1,s1,f2,s2,e1,t1,e2,t2
          j=j+1
          z(j)=d
          zmean=zmean+d
        end if
      end do
      close (unit=ifile,status='delete')
      endfile jfile
      i=ifile
      ifile=jfile
      jfile=i
c
c Recompile Delta(F) distribution statistics and histogram.
c
      if (nrej.eq.0) then
        write (ilp,'(/1x,
     &''No data trimmed from distribution tails with  zlimit = '',e10.3)
     &') zlimit
      else
        ndata=ndata-nrej
        call sort (ndata,z,index)
        zmin=z(index(1))
        zmax=z(index(ndata))
        zmedian=z(index(nint(0.5*ndata)))
        zmean=zmean/ndata
        do j=1,jmax
          xh(j)=0
          yh(j)=0
        end do
        j=0
        zj=zmin
        d=(zmax-zmin)/jmax
        do i=1,ndata
          zi=z(index(i))
          do while (zi.ge.zj)
            j=j+1
            zj=zj+d
          end do
          if (j.le.jmax) then
            xh(j)=xh(j)+zi
            yh(j)=yh(j)+1
          end if
        end do
        xh(jmax)=xh(jmax)+zi
        yh(jmax)=yh(jmax)+1
        if (type.eq.'SAS') then
          j=0
          zj=zmin
          do while (zj.lt.0)
            j=j+1
            zj=zj+d
          end do
          yh(j)=yh(j)-nc
        end if
        yhmax=0
        do j=1,jmax
          if (yh(j).gt.0) then
            xh(j)=xh(j)/yh(j)
            yhmax=max(yhmax,yh(j))
          end if
        end do
        write (ilp,'(/1x,
     &''Delta(F) Distribution Statistics''/1x,
     &''--------------------------------''/1x,
     &''After trimming distribution tails:''/1x,
     &''-----''//1x,
     &''               maximum abs[Delta(F)]''/1x
     &''  zlimit = ------------------------- = '',e10.3/1x
     &''           1.25*median abs[Delta(F)]''//1x,
     &''  Nrej                 = '',i10,''  rejected data pairs''/1x,
     &''  Nrejn                = '',i10,''  negative Delta(F)''/1x,
     &''  Nrejp                = '',i10,''  positive''//1x,
     &''  N                    = '',i10,''  remaining data pairs''/1x,
     &''  minimum Delta(F)     = '',e10.3/1x,
     &''  maximum              = '',e10.3/1x,
     &''  mean                 = '',e10.3)')
     & zlimit,nrej,nrejn,nrejp,ndata,zmin,zmax,zmean
        zmean=0
        do i=1,ndata
          d=abs(z(i))
          zmean=zmean+d
          z(i)=d
        end do
        call sort (ndata,z,index)
        zmedian=z(index(nint(0.5*(ndata+nz))))
        zmean=zmean/(ndata-nz)
        write (ilp,'(1x,
     &  ''  median abs[Delta(F)] = '',e10.3/1x,
     &  ''  mean                 = '',e10.3)')
     &  zmedian,zmean
        write (ilp,'(/1x,
     &  ''Delta(F) distribution after trimming distribution tails''/1x,
     &  ''-------------------------------------------------------''/1x,
     &  ''     n  <Delta(F)>''/1x,''     -  ----------'')')
        do j=1,jmax
          n=nint(50*yh(j)/yhmax)
          write (ilp,'(1x,i6,e12.3,1x,50a1)') nint(yh(j)),xh(j),
     &    ('X',i=1,n)
        end do
      end if
c
c Write an output file sorted on
c
c             abs[Delta(E)]              abs(E1 - E2)
c abs[z(E)] = ------------- = --------------------------------- 
c             sigma(Delta)    sqrt[sigma(E1)**2 + sigma(E2)**2]
c
c or, in the absence of E values, the analogous z(F).
c
      rewind ifile
      do i=1,4
        read (ifile,'(a1)') dummy
      end do
      read (ifile,*) ih,ik,il,f1,s1,f2,s2,e1,t1,e2,t2
      if (e1.gt.0.and.e2.gt.0) then
        ie=1
      else
        ie=0
      end if
      rewind ifile
      do i=1,4
        read (ifile,'(a1)') dummy
      end do
      open (unit=jfile,status='scratch',form='unformatted',
     & access='direct',recl=11)
      do i=1,ndata
        read (ifile,*) ih,ik,il,f1,s1,f2,s2,e1,t1,e2,t2
        if (ie.eq.1) then
          z(i)=abs(e1-e2)/sqrt(t1**2+t2**2)
        else
          z(i)=abs(f1-f2)/sqrt(s1**2+s2**2)
        end if
        write (jfile,rec=i) ih,ik,il,f1,s1,f2,s2,e1,t1,e2,t2
      end do
      call sort (ndata,z,index)
      close (unit=ifile,status='keep')
      open (unit=ifile,file='sort.locscl',status='new')
      write (ifile,'(1x,
     &''Program locscl.  '',a,'', '',a,''.  '',a//1x,
     &''    h    k    l          F1  sigmaF1          F2  sigmaF2'',
     &''    E1 sigmaE1    E2 sigmaE2 Delta/sigma(Delta)''/1x,
     &''    -    -    -          --  -------          --  -------'',
     &''    -- -------    -- ------- ------------------'')')
     &atime,adate,title
      do i=ndata,1,-1
        j=index(i)
        read  (jfile,rec=j) ih,ik,il,f1,s1,f2,s2,e1,t1,e2,t2
        if (ie.eq.1) then
          z(j)=sign(z(j),(e1-e2))
        else
          z(j)=sign(z(j),(f1-f2))
        end if
        write (ifile,'(1x,3i5,2(e12.5,e9.2),2(2f7.3),f8.2)')
     &  ih,ik,il,f1,s1,f2,s2,e1,t1,e2,t2,z(j)
      end do
      close (unit=ifile,status='keep')
      close (jfile,status='delete')
c
c Write an SAS output file with separate records for +h+k+l and -h-k-l.
c
      if (type.eq.'SAS') then
        open (unit=ifile,file='data.locscl',status='old')
        open (unit=jfile,file='hphn.locscl',status='new')
        write (jfile,'(1x,
     &  ''Program locscl.  '',a,'', '',a,''.  '',a//1x,
     &  ''    h    k    l           F   sigmaF      E sigmaE''/1x,
     &  ''    -    -    -           -   ------      - ------'')')
     &  atime,adate,title
        do i=1,4
          read (ifile,'(a1)') dummy
        end do
        do i=1,ndata
          read (ifile,*) ih,ik,il,f1,s1,f2,s2,e1,t1,e2,t2
          write (jfile,'(1x,3i5,e12.5,e9.2,2f7.3)')
     &     +ih,+ik,+il,f1,s1,e1,t1
          call csym (ptgp,ih,ik,il,ic)
          if (ic.eq.0) write (jfile,'(1x,3i5,e12.5,e9.2,2f7.3)')
     &     -ih,-ik,-il,f2,s2,e2,t2
        end do
        close (jfile,status='keep')
      end if
      close (ifile,status='keep')
      return
      end
c-----------------------------------------------------------------------
      subroutine pscale (imax,data,minhkl,maxhkl,cutoff,ginv,ptgp,
     & ihmin,ikmin,ilmin,nk,nl,ilp)
c
c Global pre-scaling via a logarithmically linearized least-squares fit
c of empirical scaling parameters p0, p1, and p2 to minimize
c   chi**2 = sum w*[(F1/F2) - p]**2,
c where
c   p = p0*exp(p1*s**2 + p2*s**4)    and    s = sin(theta)/lambda.
c
c p1 approximates <Biso2> - <Biso1>, and
c p2 approximates <(Biso1 - <Biso1>)**2> - <(Biso2 - <Biso2>)**2>.
c
      character ptgp*5
      dimension data(10,0:imax),ginv(3,3),a(3,3),ainv(3,3),b(3),c(3)
      real*8 sumw,sumx,sumy,sumx2,sumx3,sumx4,sumy2,sumxy,sumx2y,
     & sump,sumv,sumpsq
      n=0
      sumw=0
      sumx=0
      sumy=0
      sumx2=0
      sumx3=0
      sumx4=0
      sumy2=0
      sumxy=0
      sumx2y=0
      do ihkl=minhkl,maxhkl
c
c Get F and sigma(F) values.
c
        f1=data(1,ihkl)
        s1=data(2,ihkl)
        f2=data(3,ihkl)
        s2=data(4,ihkl)
c
c Are both reflections significant?
c
        if (s1.gt.0.and.s2.gt.0.and.f1.gt.cutoff*s1.and.f2.gt.cutoff*s2)
     &  then
c
c Unpack indices,
c
c ihkl = (ih - ihmin)*nk*nl + (ik - ikmin)*nl + (il - ilmin).
c
          ih=ihkl/(nk*nl)
          ik=(ihkl-ih*nk*nl)/nl
          il=ihkl-ih*nk*nl-ik*nl
          ih=ih+ihmin
          ik=ik+ikmin
          il=il+ilmin
c
c Calculate sin(theta)/lambda.
c
          s=sinthl(ih,ik,il,ginv)
c
c Get point group multiplicity.
c
          call esym (ptgp,ih,ik,il,epsilon,m)
c
c Acentric (ic = 0) or centric (ic = 1) reflection?
c
          call csym (ptgp,ih,ik,il,ic)
c
c Double the multiplicity for acentric reflections.
c
          if (ic.eq.0) m=2*m
c
c Store data for logarithmic fit.
c                                    
          x=s**2
          y=log(f1/f2)
          w=m/((s1/f1)**2+(s2/f2)**2)
c
c Accumulate least-squares sums.
c
          n=n+m
          sumw=sumw+w
          sumx=sumx+w*x
          sumy=sumy+w*y
          sumx2=sumx2+w*x**2
          sumx3=sumx3+w*x**3
          sumx4=sumx4+w*x**4
          sumy2=sumy2+w*y**2
          sumxy=sumxy+w*x*y
          sumx2y=sumx2y+w*x**2*y
        end if
      end do
c
c Store and then invert the normal matrix
c
      a(1,1)=sumw
      a(1,2)=sumx
      a(1,3)=sumx2
      a(2,1)=a(1,2)
      a(2,2)=sumx2
      a(2,3)=sumx3
      a(3,1)=a(1,3)
      a(3,2)=a(2,3)
      a(3,3)=sumx4
      call minv3 (a,ainv,det)
c
c Store the normal vector and then evaluate the parameters and their
c esd's.
c
      b(1)=sumy
      b(2)=sumxy
      b(3)=sumx2y
      call mv (ainv,b,c)
      p0=c(1)
      p1=c(2)
      p2=c(3)
c
c chisq = sum w*(y - p)**2
c       = sum w*y**2 + sum w*p**2 - 2*sum w*p*y
c
c p     = p0 + p1*x + p2*x**2
c
c p**2  = p0**2 + (p1*x + p2*x**2)**2 + 2*p0*(p1*x + p2*x**2)
c       = p0**2 + p1**2*x**2 + p2**2*x**4 + 2*p1*p2*x**3 
c                                     + 2*p0*p1*x + 2*p0*p2*x**2
c
      chisq=sumy2+p0**2*sumw+p1**2*sumx2+p2**2*sumx4
     &              +2*(p0*p1*sumx+p0*p2*sumx2+p1*p2*sumx3)
     &                -2*(p0*sumy+p1*sumxy+p2*sumx2y)
      zsq=chisq/(n-3)
      v0=zsq*ainv(1,1)
      v1=zsq*ainv(2,2)
      v2=zsq*ainv(3,3)
      cov01=zsq*ainv(1,2)
      cov02=zsq*ainv(1,3)
      cov12=zsq*ainv(2,3)
c
c If the p2 coefficient of the term quadratic in s**2 is not
c statistically significant, resort to a scaling function linear in s**2.
c
      if (abs(p2).lt.sqrt(v2)) then
        det=sumw*sumx2-sumx*sumx
        p0=(sumx2*sumy-sumx*sumxy)/det
        p1=(-sumx*sumy+sumw*sumxy)/det
        chisq=sumy2+p0**2*sumw+p1**2*sumx2
     &                +2*p0*p1*sumx
     &                  -2*(p0*sumy+p1*sumxy)
        zsq=chisq/(n-2)
        v0=zsq*sumx2/det
        v1=zsq*sumw/det
        cov01=zsq*(-sumx/det)
        p2=0
        v2=0
        cov02=0
        cov12=0
c
c If the p1 coefficient of the term linear in s**2 is not statistically
c significant, resort to a constant scaling function.
c
        if (abs(p1).lt.sqrt(v1)) then
          p0=sumy/sumw
          p1=0
          chisq=sumy2+p0**2*sumw-2*p0*sumy
          zsq=chisq/(n-1)
          v0=zsq/sumw
          v1=0
          cov01=0
        end if
      end if
c
c Scale the F2(derivative) values, keeping constant the F1(native)
c values.
c
      sumx=0
      sumy=0
      pmin=+1e9
      pmax=-1e9
      sumn=0
      sump=0
      sumv=0
      sumpsq=0
      do ihkl=minhkl,maxhkl
        f1=data(1,ihkl)
        s1=data(2,ihkl)
        f2=data(3,ihkl)
        s2=data(4,ihkl)
        if (s1.gt.0.and.s2.gt.0.and.f1.ge.cutoff*s1.and.f2.ge.cutoff*s2)
     &  then
          ih=ihkl/(nk*nl)
          ik=(ihkl-ih*nk*nl)/nl
          il=ihkl-ih*nk*nl-ik*nl
          ih=ih+ihmin
          ik=ik+ikmin
          il=il+ilmin
          s=sinthl(ih,ik,il,ginv)
          p=p0+p1*s**2+p2*s**4
          v=v0+v1*s**4+v2*s**8+2*(cov01*s**2+cov02*s**4+cov12*s**6)
c
c Convert from logarithmic fitting scale to antilogarithmic data scale.
c
          p=exp(p)
          v=p**2*v
c
c Apply scaling function, and compile scaling statistics.
c
          s2=sqrt((p*s2)**2+f2**2*v)
          f2=p*f2
          data(3,ihkl)=f2
          data(4,ihkl)=s2
          sumx=sumx+f1
          sumy=sumy+f2
          if (p.lt.pmin) then
            pmin=p
            sigmin=sqrt(v)
          end if
          if (p.gt.pmax) then
            pmax=p
            sigmax=sqrt(v)
          end if
          sumn=sumn+1
          sump=sump+p
          sumv=sumv+v
          sumpsq=sumpsq+p**2
        end if
      end do
c
c Re-scale the re-scaled data to adjust from the reciprocal-variance-
c weighted, logarithmic fitting scale to the unit-weighted data scale.
c
      q=sumx/sumy
      pmin=pmin*q
      pmax=pmax*q
      sigmin=sigmin*q
      sigmax=sigmax*q
      sump=sump*q
      sumv=sumv*q**2
      sumpsq=sumpsq*q**2
      sumx=0
      sumy=0
      sumw=0
      do ihkl=minhkl,maxhkl
        f1=data(1,ihkl)
        s1=data(2,ihkl)
        f2=data(3,ihkl)
        s2=data(4,ihkl)
        if (s1.gt.0.and.s2.gt.0.and.f1.ge.cutoff*s1.and.f2.ge.cutoff*s2)
     &  then
          f2=f2*q
          s2=s2*q
          data(3,ihkl)=f2
          data(4,ihkl)=s2
          sumw=sumw+abs(f1-f2)
          sumx=sumx+f1
          sumy=sumy+f2
        end if
      end do          
      write (ilp,'(//1x,
     &''Global Pre-scaling''/1x,''------------------''/1x,
     &''Empirical scaling parameters p0, p1, and p2 obtained by a''/1x,
     &''logarithmically linearized least-squares fit to minimize''/1x,
     &''  chi**2 = sum w*[(F1/F2) - p]**2 ,''/1x,''where''/1x,
     &''  p = p0*exp(p1*s**2 + p2*s**4)    and    s = sin(theta)/lamb'',
     &''da .''//1x,
     &''p1 approximates <Biso(2)> - <Biso(1)>, and''/1x,
     &''p2 approximates <(Biso(1) - <Biso(1)>)**2> - <(Biso(2) - <Bis'',
     &''o(2)>)**2>.''//1x,
     &''p0 = '',e10.3/1x,
     &''p1 = '',e10.3/1x,
     &''p2 = '',e10.3,''    sqrt(abs(p2)) = '',e10.3)')
     & q*exp(p0),p1,p2,sqrt(abs(p2))
      write (ilp,'(/1x,
     &''Pre-scaling Statistics:''/1x,
     &''-----------------------''/1x,
     &''pmin = '',e10.3,''  sigma(pmin) = ''e10.3/1x,
     &''pmax = '',e10.3,''  sigma(pmax) = ''e10.3//1x,
     &''<p>                   = '',e10.3/1x,
     &''<sigma(p)**2>**(1/2)  = '',e10.3/1x,
     &''<(p - <p>)**2>**(1/2) = '',e10.3//1x,
     &''Q = <F1>/<F2*p>                        = '',e10.3/1x,
     &''R = <abs(F1 - F2*p)>/<0.5*(F1 + F2*p)> = '',e10.3)')
     & pmin,sigmin,pmax,sigmax,
     & sump/sumn,sqrt(sumv/sumn),sqrt((sumpsq/sumn)-(sump/sumn)**2),
     & sumx/sumy,sumw/(0.5*(sumx+sumy))
      return
      end
c-----------------------------------------------------------------------
      subroutine qfit (imax,data,minhkl,maxhkl,ihmin,ikmin,ilmin,
     & nk,nl,mh,mk,ml,nrej,ndata,qmin,sigmin,qmax,sigmax,qmean,rmsd,esd,
     & wqmean,wrmsd,wesd)
c
c Fit least-squares local scale factors in blocks of
c   n = [2*m(h) + 1]*[2*m(k) + 1]*[2*m(l) + 1] - 1
c reciprocal lattice points surrounding, but excluding, each point hkl.
c
      dimension data(10,0:imax)
      parameter (jmax=25)
      dimension t(4,-jmax:+jmax,-jmax:+jmax,-jmax:+jmax)
      real*8 sumw,sumx,sumxsq,sumq,sumqsq,sumv,wsum,wsumq,wsumqsq
c
c Loop through the data set.
c
      nrej=0
      ndata=0
      qmin=+1e9
      qmax=-1e9
      sumq=0
      sumqsq=0
      sumv=0
      wsum=0
      wsumq=0
      wsumqsq=0
      do ihkl=minhkl,maxhkl
c
c Get F and sigma(F) values.
c
        f1=data(1,ihkl)
        s1=data(2,ihkl)
        f2=data(3,ihkl)
        s2=data(4,ihkl)
        if (f1.gt.0.and.s1.gt.0.and.f2.gt.0.and.s2.gt.0) then
c
c Fetch local block
c
c   -mh .le. jh .le. +mh
c   -mk .le. jk .le. +mk
c   -ml .le. jj .le. +ml
c
c surrounding ih,ik,il stored by packed indices:
c
c   ihkl = (ih - ihmin)*nk*nl + (ik - ikmin)*nl + (il - ilmin)
c
c   jhkl = (ih - ihmin + jh)*nk*nl
c                             + (ik - ikmin + jk)*nl
c                                               + (il - ilmin + jl)
c
c   jhkl = ihkl + jh*nk*nl + jk*nl + jl
c
          do jl=-ml,+ml
          do jk=-mk,+mk
          do jh=-mh,+mh
            do j=1,4
              t(j,jh,jk,jl)=0
            end do
            jhkl=ihkl+jh*nk*nl+jk*nl+jl
            if (jhkl.ne.ihkl.and.
     &        minhkl.le.jhkl.and.jhkl.le.maxhkl) then
              do j=1,4
                t(j,jh,jk,jl)=data(j,jhkl)
              end do
            end if
          end do
          end do
          end do
c
c Accumulate local least-squares sums for reflections surrounding, but
c excluding, the given reflection.
c
          n=0
          sumw=0
          sumx=0
          sumxsq=0
          do jl=-ml,+ml
          do jk=-mk,+mk
          do jh=-mh,+mh
            if (jh.ne.0.or.jk.ne.0.or.jl.ne.0) then
              f1=t(1,jh,jk,jl)
              s1=t(2,jh,jk,jl)
              f2=t(3,jh,jk,jl)
              s2=t(4,jh,jk,jl)
              if (f1.gt.0.and.s1.gt.0.and.f2.gt.0.and.s2.gt.0) then
                x=f1/f2
                w=1/(x**2*((s1/f1)**2+(s2/f2)**2))
                n=n+1
                sumw=sumw+w
                sumx=sumx+w*x
                sumxsq=sumxsq+w*x**2
              end if
            end if
          end do
          end do
          end do
c
c Test for ill-determined local scaling data.
c
          if (n.lt.2.or.sumw.le.0.or.sumx.le.0) then
            q=0
            chisq=0
          else
c
c Evaluate the local scale factor, q = <F1/F2>, and the local mean-
c square error of fit.
c
            q=sumx/sumw
            chisq=sumxsq-2*q*sumx+q**2*sumw
          end if
c
c Test for ill-determined error of fit.
c
          if (chisq.le.0) then
            q=0
            sigmaq=0
            nrej=nrej+1
          else
c
c Evaluate the local error-of-fit uncertainty, sigma(q).
c
            sigmaq=sqrt(chisq/((n-1)*sumw))
c
c Accumulate the q distribution statistics.
c
            ndata=ndata+1
            sumq=sumq+q
            sumqsq=sumqsq+q**2
            sumv=sumv+sigmaq**2
            w=1/sigmaq**2
            wsum=wsum+w
            wsumq=wsumq+w*q
            wsumqsq=wsumqsq+w*q**2
            if (q.lt.qmin) then
              qmin=q
              sigmin=sigmaq
            end if
            if (q.gt.qmax) then
              qmax=q
              sigmax=sigmaq
            end if
          end if
c
c Store the local scale factor.
c
          data( 9,ihkl)=q
          data(10,ihkl)=sigmaq
        end if
      end do  
c
c Evaluate the q-distribution statistics.
c
      qmean=sumq/ndata
      rmsd=sqrt((sumqsq/ndata)-qmean**2)
      esd=sqrt(sumv/ndata)
      wqmean=wsumq/wsum
      wrmsd=sqrt((wsumqsq/wsum)-wqmean**2)
      wesd=sqrt(ndata/wsum)
      return
      end
c-----------------------------------------------------------------------      
      subroutine qscale (imax,data,minhkl,maxhkl,ihmin,ikmin,ilmin,
     & nk,nl,ifile,jfile,i012,type,ptgp,ndata,nc,ilp)
c
c Apply the local scale factors.
c
      dimension data(10,0:imax)
      character type*3,ptgp*5
      real*8 sumdf,sumde,sumf1,sumf2,sume,sum1sq,sum2sq,zmean
c
c Loop through the data set.
c
      sumdf=0
      sumde=0
      sumf1=0
      sumf2=0
      sume=0
      sum1sq=0
      sum2sq=0
      ndata=0
      nc=0
      open (unit=ifile,file='temp1.locscl',status='new')
      do ihkl=minhkl,maxhkl
c
c Get F and sigma(F), E and sigma(E), and q and sigma(q) values.
c
        f1=data(1,ihkl)
        s1=data(2,ihkl)
        f2=data(3,ihkl)
        s2=data(4,ihkl)
        e1=data(5,ihkl)
        t1=data(6,ihkl)
        e2=data(7,ihkl)
        t2=data(8,ihkl)
             q=data( 9,ihkl)
        sigmaq=data(10,ihkl)
c
c Are both reflections present?
c
        if (f1.gt.0.and.s1.gt.0.and.f2.gt.0.and.s2.gt.0.and.q.gt.0) then
c
c Unpack indices, hkl = (h - hmin)*nk*nl + (k - kmin)*nl + (l - lmin).
c
          ih=ihkl/(nk*nl)
          ik=(ihkl-ih*nk*nl)/nl
          il=ihkl-ih*nk*nl-ik*nl
          ih=ih+ihmin
          ik=ik+ikmin
          il=il+ilmin
c
c Evaluate and write out the locally scaled data.
c
          if (i012.eq.0) then
c
c Keep constant the mean of F1 and F2.
c
            p1=f1+f2
            p2=f1+f2*q
            v1=s1**2+s2**2
            v2=s1**2+(s2*q)**2+(f2*sigmaq)**2
            p=p1/p2
            v=(v1+v2*(p1/p2)**2)/p2**2
c
c The factor p = (F1 + F2)/(F1 + q*F2) preserves the mean of F1 and F2,
c i.e., F1 + F2 = p*(F1 + q*F2).
c
            s1=sqrt((p*s1)**2+f1**2*v)
            s2=sqrt((p*q*s2)**2+(f2*q)**2*v+(f2*p*sigmaq)**2)
            f1=f1*p
            f2=f2*p*q
            if (e1.gt.0.and.e2.gt.0) then
              t1=sqrt((p*t1)**2+e1**2*v)
              t2=sqrt((p*q*t2)**2+(e2*q)**2*v+(e2*p*sigmaq)**2)
              e1=e1*p
              e2=e2*p*q
            end if
          else if (i012.eq.1) then
c
c Keep constant the value of F1.
c
            s2=sqrt((q*s2)**2+(f2*sigmaq)**2)
            f2=q*f2
            if (e2.gt.0) then
              t2=sqrt((q*t2)**2+(e2*sigmaq)**2)
              e2=q*e2
            end if
          else
c
c Keep constant the value of F2.
c
            s1=sqrt((s1**2+(sigmaq*f1/q)**2)/q**2)
            f1=f1/q
            if (e1.gt.0) then
              t1=sqrt((t1**2+(sigmaq*e1/q)**2)/q**2)
              e1=e1/q
            end if
          end if
c
c For SAS data, average centric reflection pairs.
c
          if (type.eq.'SAS') then
            call csym (ptgp,ih,ik,il,ic)
            if (ic.eq.1) then
              nc=nc+1
              f1=0.5*(f1+f2)
              e1=0.5*(e1+e2)
              s1=0.5*sqrt(s1**2+s2**2)
              t1=0.5*sqrt(t1**2+t2**2)
              f2=f1
              e2=e1
              s2=s1
              t2=t1
            end if
          end if
c
c Accumulate agreement statistics.
c
          sumdf=sumdf+abs(f1-f2)
          sumde=sumde+abs(e1-e2)
          sumf1=sumf1+f1
          sumf2=sumf2+f2
          sume=sume+e1+e2
          sum1sq=sum1sq+e1**2
          sum2sq=sum2sq+e2**2
c
c Write a locally scaled output file.
c
          ndata=ndata+1
          write (ifile,1000) ih,ik,il,f1,s1,f2,s2,e1,t1,e2,t2
        end if
      end do
 1000 format (1x,3i5,2(e12.5,e9.2),2(2f7.3))
      endfile ifile
c
c Evaluate and print the agreement statistics.
c
      rf=sumdf/(0.5*(sumf1+sumf2))
      if (sume.gt.0) then
        re=sumde/(0.5*sume)
        e1sq=sum1sq/ndata
        e2sq=sum2sq/ndata
      end if
      write (ilp,'(/1x,''After scaling:''/1x,
     &''Q(F) = <F1>/<F2>                      = '',f7.4/1x,
     &''R(F) = <abs(F1 - F2)>/<0.5*(F1 + F2)> = '',f7.4/1x,
     &''R(E) = <abs(E1 - E2)>/<0.5*(E1 + E2)> = '',f7.4/1x,
     &''<E1**2> = '',f6.3/1x,''<E2**2> = '',f6.3)')
     & sumf1/sumf2,rf,re,e1sq,e2sq
      return
      end
c-----------------------------------------------------------------------
      subroutine read1 (itype,iunit,iend,h,k,l,f,sigf,e,sige)
c
c itype = 0  ASCII, formatted     h,k,l,Fsq,sigFsq,F,sigF,E,sigE
c         1    "        "         h,k,l,Fsq,sigFsq,F,sigF
c         2    "        "         h,k,l,Fsq,sigFsq       
c         3    "        "         h,k,l,           F,sigF,E,sigE
c         4    "        "         h,k,l,           F,sigF       
c         5  Binary, unformatted  h,k,l,Fsq,sigFsq,F,sigF,E,sigE
c         6    "          "       h,k,l,Fsq,sigFsq,F,sigF
c         7    "          "       h,k,l,Fsq,sigFsq
c         8    "          "       h,k,l,           F,sigF,E,sigE
c         9    "          "       h,k,l,           F,sigF
c
      integer h
      e=0
      sige=0
 1    if (itype.eq.0) read (iunit,*,end=9) h,k,l,y,sigy,f,sigf,e,sige
      if (itype.eq.1) read (iunit,*,end=9) h,k,l,y,sigy,f,sigf
      if (itype.eq.2) read (iunit,*,end=9) h,k,l,y,sigy
      if (itype.eq.3) read (iunit,*,end=9) h,k,l,       f,sigf,e,sige
      if (itype.eq.4) read (iunit,*,end=9) h,k,l,       f,sigf
      if (itype.eq.5) read (iunit,  end=9) h,k,l,y,sigy,f,sigf,e,sige
      if (itype.eq.6) read (iunit,  end=9) h,k,l,y,sigy,f,sigf
      if (itype.eq.7) read (iunit,  end=9) h,k,l,y,sigy
      if (itype.eq.8) read (iunit,  end=9) h,k,l,       f,sigf,e,sige
      if (itype.eq.9) read (iunit,  end=9) h,k,l,       f,sigf
      if (itype.eq.2.or.itype.eq.7) then
        if (sigy.le.0) go to 1
        call fsqtof (y,sigy,f,sigf)
      end if
      if (sigf.le.0) go to 1
      return
 9    iend=1
      return
      end
c-----------------------------------------------------------------------
      function sinthl(ih,ik,il,ginv)
      dimension h(3),ginv(3,3)
      h(1)=ih
      h(2)=ik
      h(3)=il
      q=0
      do i=1,3
      do j=1,3
        q=q+h(j)*ginv(i,j)*h(i)
      end do
      end do
      sinthl=0.5*sqrt(q)
      return
      end
c-----------------------------------------------------------------------
      subroutine sort (n,data,index)
c
c Indexes the array data(n) and returns the array index(n) sorted such
c that the values data(index(i)) are in ascending order for i = 1, 2,
c ..., n.  The input variables n and data(n) are not changed.
c
c Employs the "heapsort" algorithm.
c
c Fortran code quoted from William H. Press, Brian P. Flannery, Saul A.
c Teukolsky, and William T. Vetterling (1986).  Numerical Recipies:  The
c                                               --------- --------   ---
c Art of Scientific Computing, pp. 229-233.  Cambridge, England:
c --- -- ---------- ---------
c Cambridge University Press.
c
      dimension data(n),index(n)
      do i=1,n
        index(i)=i
      end do
      i=n/2+1
      j=n
 1    continue
      if (i.gt.1) then
        i=i-1
        inext=index(i)
        t=data(inext)
      else
        inext=index(j)
        t=data(inext)
        index(j)=index(1)
        j=j-1
        if (j.eq.1) then
          index(1)=inext
          return
        end if
      end if
      k=i
      l=i+i
      do while (l.le.j)
        if (l.lt.j) then
          if (data(index(l)).lt.data(index(l+1))) l=l+1
        end if
        if (t.lt.data(index(l))) then
          index(k)=index(l)
          k=l
          l=l+l
        else
          l=j+1
        end if
      end do
      index(k)=inext
      go to 1
      end
c-----------------------------------------------------------------------
