program diffE
c
c Robert H. Blessing
c Hauptman-Woodward Institute
c 73 High Street
c Buffalo, New York 14203, USA
c Telephone:  (716) 856-9600, extension 335
c Electronic Mail:  blessing@hwi.buffalo.edu
c
c July 1996,..., February 1999
c
c Evaluate SIR or SAS difference E values.
c
      character atime*8,adate*9,atitle*80,btime*8,bdate*9,btitle*60,
     & file*60,locscl*60,type*3,evaldat(2)*60
      character syst*12,latt*1,ptgp*5,laue*5,xray*2,el*2
      dimension cell(6),astar(6),ginv(3,3),h(3)
      dimension el(40),xi(40),zi(40),a1(40),b1(40),a2(40),b2(40),a3(40),
     & b3(40),a4(40),b4(40),c0(40),fp(40),fpp(40)
      double precision sumf,sumfsq,sumfpp,sumfppsq
      double precision fder,fsqder,fnat,fsqnat
      dimension data(100000),index(100000)
      data data,index /200000*0/
      data nlist /100/
      call time (atime)
      call date (adate)
c      atime='hh:mm:ss'
c      adate='dd-mmm-yy'
      io1=10
      io2=20
      ilp=60
c
c Read control data file "diffe.dat".
c
      open (unit=io1,file='diffe.dat',status='old')
      read (io1,'(a)') atitle
      read (io1,'(a)') locscl
      read (io1,'(a)') type
      if (type.eq.'sir') type='SIR'
      if (type.eq.'sas') type='SAS'
      read (io1,'(a)') evaldat(1)
      if (type.eq.'SIR') read (io1,'(a)') evaldat(2)
      data tmax,xmin,ymin,zmin,zmax,smin,smax /7*0/
      read (io1,*,end=5) tmax
      read (io1,*,end=5) xmin
      read (io1,*,end=5) ymin
      read (io1,*,end=5) zmin
      read (io1,*,end=5) zmax
      read (io1,*,end=5) smin,smax
 5    close (unit=io1,status='keep')
c-----------------------------------------------------------------------
c
c Control data:
c -------------
c
c atitle (a) = Job title
c
c locscl (a) = Name of the locally scaled reflections file "data.locscl"
c              or "sort.locscl" from the program "locscl", an ascii,
c              formatted file read by program "diffe" under free format
c              with records:
c
c              h,k,l,F1,sigma(F1),F2,sigma(F2),E1,sigma(E1),E2,sigma(E2)
c
c type (a)   = SIR or SAS
c
c              For an SIR case,
c              the E1 data should be those for the native crystal, and
c              the E2 data should be those for the derivative crystal.
c
c              For an SAS case,
c              the E1 data should correspond to +h+k+l, and
c              the E2 data should correspond to -h-k-l.
c
c evaldat(1) (a) = Name of the "eval.dat" file from program "levy" (or
c                  program "rogers") for native SIR data or SAS data.
c                                        ------ ---
c evaldat(2) (a) = Name of the "eval.dat" file for derivative SIR data.
c                                                  ---------- ---
c                  For an SAS case, omit the evaldat(2) record.
c
c                  From the "eval.dat" file(s) the program gets the
c                  crystal data, the unit cell chemical composition, and
c                  the atomic scattering factor coefficients.
c 
c tmax (*) = Maximum value of
c
c            abs(Delta)          abs(Delta)
c            ---------- = -----------------------
c              sigma      1.25*median[abs(Delta)]
c
c                                     abs[E1 - E2 - median(E1 - E2)]
c                       = ------------------------------------------- ,
c                         1.25*median{abs[E1 - E2 - median(E1 - E2)]}
c
c            for reflection pairs to be accepted for renormalization.
c            The variable tmax is used to exclude data with unreliably
c            large values of abs(E1 - E2) in the tails of the data-set
c            (E1 - E2) distribution.
c
c            The  tmax  test assumes that the data-set  t = Delta/sigma
c            distribution should approximate a zero-mean, unit-variance
c            normal distribution, for which  values of t .lt. -tmax  or
c            t .gt. +tmax  are extremely improbable.
c
c            Default value:            tmax = 99.9
c            Recommended trial value:  tmax = 6
c
c xmin (*) = Minimum value for
c
c            min[E1/sigma(E1), E2/sigma(E2)]
c
c            for reflection pairs to be included in the program "diffE"
c            processing.  Default value:  xmin = 0.  Recommended trial
c            value:  xmin = 3.
c
c ymin (*) = Minimum value for
c
c            abs(E1 - E2)/sqrt[sigma(E1)**2 + sigma(E2)**2]
c
c            for reflection pairs to be included in the program "diffE"
c            processing.  Default value:  ymin = 0.  Recommended trial
c            value:  ymin = 1.
c
c zmin (*) = Minimum value for
c
c            diffE/sigma(diffE)
c
c            for reflection pairs to be included in the diffE-sorted
c            output reflections file "sort.diffe".  Default value:
c            zmin = 0.  Recommended  trial value:  zmin = 3.
c
c zmax (*) = Maximum value for
c
c            (diffE - diffEmax)/sigma(diffE)
c
c            for reflection pairs to be included in the diffE-sorted
c            output reflections file "sort.diffe".  Default value:
c            zmax = 0.  Largest recommended trial value:  zmax = 3.
c
c            For SIR cases
c
c                                   abs{sum[f(Der)] - sum[f(Nat)}
c        diffEmax = -------------------------------------------------- ,
c                   sqrt(epsilon*abs{sum[f**2(Der)] - sum[f**2(Nat)]})
c
c            and for SAS cases
c
c                            sum(fpp)
c        diffEmax = ------------------------- .
c                   sqrt[epsilon*sum(fpp**2)]
c
c smin, smax (*) = Minimum and maximum limits for s = sin(theta)/lambda
c                  for reflection pairs to be included in the program
c                  "diffE" processing.
c
c                  For default values smin = smax = 0, the exclusion
c                  test is not applied.
c
c                  Exclusion based on sin(theta)/lambda limits is
c                  applied only if smin .gt. 0 .or. smax .gt. 0.
c
c The user might need to try several runs of the program, experimenting
c with the input cutoff data selection variables tmax, xmin, ymin, and
c smax, in order to arrive an acceptable distribution of diffE values,
c as judged by the plots of the diffE distributions and the lists of
c distribution statistics recorded in the "diffe.lp" output file.
c
c In addition to using the output cutoff variable zmin to reject data
c that have large estimated standard uncertainties, the user should
c inspect the (descending-order) diffE-ranked, formatted, ascii output
c reflection file "sort.diffe", and delete any reflection pairs at the
c top of the list that have diffE larger than about three or four.  Such
c large diffE values are very improbable and are almost certainly due to
c measurement errors.
c
c-----------------------------------------------------------------------
      open (unit=ilp,file='diffe.lp',status='unknown')
      write (ilp,600) atime,adate,atitle
 600  format ('1'/1x,'Program diffE.  ',a,', ',a,'.  ',a)
      if (type.eq.'SIR') then
        write (ilp,921)
 921    format (/1x,'SIR difference E values,'//1x,
     &'        abs[sqrt{sum[f**2(Der)]}*E(Der) - sqrt{sum[f**2(Nat)]',
     &'}*E(Nat)]'/1x,
     &'diffE = -----------------------------------------------------',
     &'--------- ,'/1x,
     &'              q*sqrt[abs{sum[f**2(Der)] - sum[f**2(Nat)]}]'//
     &1x,
     &'where q is a renormalization scaling function such that <diff',
     &'E**2> = 1.')
        write (ilp,922) locscl,evaldat(2),evaldat(1)
 922    format (//1x, 'Input files:'/1x, '------------'/1x,
     &'    Locally scaled SIR reflection pairs file:  ',a/1x,
     &'                   ---'/1x,
     &'    Derivative "eval.dat" file:                ',a/1x,
     &'    ----------'/1x,
     &'    Native "eval.dat" file:                    ',a/1x,
     &'    ------')
      end if
      if (type.eq.'SAS') then
        write (ilp,923)
 923    format (/1x, 'SAS difference E values,' //1x,
     &'        sqrt{sum[(f0 + fp)**2 + fpp**2])}*abs[E(+h) - E(-h)]'/1x,
     &'diffE = ---------------------------------------------------- ,'
     &/1x, '                       2*q*sqrt[sum(fpp**2)]'//1x,
     &'where q is a renormalization scaling function such that <diff',
     &'E**2> = 1.')
        write (ilp,924) locscl,evaldat(1)
 924    format (//1x, 'Input files:'/1x, '----- -----'/1x,
     &'    Locally scaled SAS reflection pairs file:  ',a/1x,
     &'                   ---'/1x,
     &'    "eval.dat" file:                           ',a)
      end if
c
c Read and echo crystal and scattering factor data from "eval.dat"
c file(s) from program "levy" or program "rogers".
c
      if (type.eq.'SIR') n=2
      if (type.eq.'SAS') n=1
      i=0
      do j=n,1,-1
c
c In an SIR case, read first the "eval.dat" file for the derivative
c crystal (j = 2), and then the file for the native crystal (j = 1).
c
        open (unit=io1,file=evaldat(j),status='old')
        read (io1,'(1x,a)') btime
        read (io1,'(1x,a)') bdate
        read (io1,'(1x,a)') title
        read (io1,*)        itype
        read (io1,'(1x,a)') file
        read (io1,'(1x,a)') latt
        read (io1,'(1x,a)') ptgp
        read (io1,*)        cell
        read (io1,'(1x,a)') xray
        read (io1,*)        nel
        do k=1,nel
          i=i+1
          read (io1,'(1x,a2,f10.2,f6.0,4(1x,f10.6))')
     &     el(i),xi(i),zi(i),a1(i),b1(i),a2(i),b2(i)
          read (io1,'(1x,5(1x,f10.6),2x,2(1x,f10.6))')
     &     a3(i),b3(i),a4(i),b4(i),c0(i),fp(i),fpp(i)
        end do
        if (j.eq.2) nder=nel
        if (j.eq.1) nnat=nel
        call decode (ptgp,laue,syst,latt,mlatt,iptgp)
        if (type.eq.'SIR') then
          if (j.eq.2) write (ilp,913)
  913 format (/1x,'SIR Derivative crystal and scattering ',
     &'factor data:'/1x,'--------------')
          if (j.eq.1) write (ilp,914)
  914 format (/1x,'SIR Native crystal and scattering ',
     &'factor data:'/1x,'---------')
        end if
        write (ilp,898) btime,bdate,title
        write (ilp,897) file
        write (ilp,896) syst,latt,laue,ptgp
        if (ptgp.eq.laue) then
          write (ilp,8961)
        else
          write (ilp,8960)
        end if
        write (ilp,893) cell
        if (j.eq.2) k1=1
        if (j.eq.1) k1=1+nder
        k2=k1+nel-1
        write (ilp,892) (el(k),xi(k),k=k1,k2)
        write (ilp,891) xray,(el(k),fp(k),fpp(k),k=k1,k2)
        close (unit=io1,status='keep')
      end do
      call metric (cell,astar,ginv)
 898  format ('0Crystal and scattering factors data from "eval.dat" fi',
     &'le from program rogers or program levy job:'/'0',4x,a,1x,a,'.  ',
     &a)
 897  format ('0',4x,'Input reflection file:   ',a)
 896  format ('0',a/'0',a,'-lattice'/'0Laue group ',a/'0Point group ',a)
 8960 format ('0Noncentrosymmetric structure')
 8961 format ('0Centrosymmetric structure')
 893  format ('0Lattice Parameters'/'          a         b         c',
     &'     alpha      beta     gamma'/' ',3f10.4,3f10.3)
 892  format ('0Unit Cell Contents:'//(' ',a2,f10.2))
 891  format (/1x,'Radiation:  ',a2//1x,'Anomalous Dispersion Correcti',
     &'on Terms:'/1x,'--------------------------------------'/1x,
     &'el        fp       fpp'/1x,'--        --       ---'/
     &(1x,a2,2f10.3))
c
c Select input data and calculate initial difference E values.
c
      if (tmax.le.0) tmax=99.9
      if (xmin.lt.0) xmin=0
      if (ymin.lt.0) ymin=0
      if (smin.le.0) then
        smin=0
        dmax=999.9
      else
        dmax=1/(2*smin)
      end if
      if (smax.le.0) then
        smax=9.9999
        dmin=0
      else
        dmin=1/(2*smax)
      end if
      if (zmin.lt.0) zmin=0
      if (zmax.lt.0) zmax=0
      write (ilp,600) atime,adate,atitle
      write (ilp,925) tmax,xmin,ymin,smin,dmax,smax,dmin
 925  format (//1x, 'Data selection:'/1x, '---------------'/1x,
     &'Accept for input to difference renormalization only data pair',
     &'s with:'//1x,
     &'              abs[E1 - E2 - median(E1 - E2)]'/1x,
     &'  -------------------------------------------    .le. tmax,  ',
     &'  tmax = ',f4.1/1x,
     &'  1.25*median{abs[E1 - E2 - median(E1 - E2)]}'//1x,
     &'  min{[E1/sigma(E1)], [E2/sigma(E2)]}            .ge. xmin,  ',
     &'  xmin = ',f4.1//1x,
     &'  abs(E2 - E1)/sqrt[sigma(E2)**2 + sigma(E1)**2] .ge. ymin,  ',
     &'  ymin = ',f4.1//1x,
     &'  smin .le. sin(theta)/lambda .le. smax,                     ',
     &'  smin = ',f7.4,' reciprocal Angstroms      dmax = ',f5.1,
     &'  Angstroms'/1x,
     &'                                                             ',
     &'  smax = ',f7.4,'                           dmin = ',f6.2)
      write (ilp,926) zmin,zmax
 926  format (/1x,
     &'Accept for output to phasing trials only data with:'//1x,
     &'  diffE/sigma(diffE)              .ge. zmin,                 ',
     &'  zmin = ',f4.1//1x,
     &'  (diffE - diffEmax)/sigma(diffE) .le. zmax,                 ',
     &'  zmax = ',f4.1)
      if (type.eq.'SIR') then
        write (ilp,927)
 927    format (/1x,
     &'                             abs{sum[f(Der)] - sum[f(Nat)}'/1x,
     &'  diffEmax = ------------------------------------------------',
     &'--'/1x,
     &'             sqrt(epsilon*abs{sum[f**2(Der)] - sum[f**2(Nat)]',
     &'})')
      end if
      if (type.eq.'SAS') then
        write (ilp,928)
 928    format (/1x, '                   sum(fpp)'/1x,
     &'  diffEmax = -------------------------'/1x,
     &'             sqrt[epsilon*sum(fpp**2)]')
      end if
      open (unit=io1,file=locscl,status='old')
      open (unit=io2,status='scratch',form='unformatted')
      call datain (io1,io2,ilp,type,iptgp,tmax)
      n0=0
      n1=0
      n2=0
      n3=0
      na=0
      nc=0
      ndata=0
      rewind io1
 1    read (io1,end=9) j,k,l,f1,sigf1,f2,sigf2,e1,sige1,e2,sige2
c
c Apply input data selection tests.
c
      if (min(sige1,sige2).le.0) go to 1
      n0=n0+1
      if (min(e1/sige1,e2/sige2).lt.xmin) go to 1
      n1=n1+1
      if (abs(e1-e2)/sqrt(sige1**2+sige2**2).lt.ymin) go to 1
      n2=n2+1
      s=sinthl(j,k,l,ginv)
      if ((smin.gt.0.and.s.lt.smin).or.(smax.gt.0.and.s.gt.smax))
     & go to 1
      n3=n3+1
c
c Calculate initial difference E values.
c
      if (type.eq.'SIR') then
c
c For SIR data, calculate sum(f) and sum(f**2) for both the derivative
c (F2) and the native (F1) structures.
c
        nel=nder
        call fcalc (nel,xi,a1,b1,a2,b2,a3,b3,a4,b4,c0,fp,fpp,s,
     &   sumf,sumfsq,sumfpp,sumfppsq)
        fder=sumf
        fsqder=sumfsq
        i=nel+1
        nel=nnat
        call fcalc (nel,xi(i),a1(i),b1(i),a2(i),b2(i),a3(i),b3(i),
     &  a4(i),b4(i),c0(i),fp(i),fpp(i),s,sumf,sumfsq,sumfpp,sumfppsq)
        fnat=sumf
        fsqnat=sumfsq
        ddiffe=abs(sqrt(fsqder)*e2-sqrt(fsqnat)*e1)/
     &          sqrt(abs(fsqder-fsqnat))
        sigdif=sqrt(fsqder*sige2**2+fsqnat*sige1**2)/
     &           sqrt(abs(fsqder-fsqnat))
        call esym (iptgp,j,k,l,epsilon,m)
        epsilon=mlatt*epsilon
        difmax=abs(fder-fnat)/sqrt(epsilon*abs(fsqder-fsqnat))
      else
c
c For SAS data, calculate sum[(f0 + fp)**2 + fpp**2], sum(fpp), and
c sum(fpp**2).
c
        call fcalc (nel,xi,a1,b1,a2,b2,a3,b3,a4,b4,c0,fp,fpp,s,
     &   sumf,sumfsq,sumfpp,sumfppsq)
        ddiffe=sqrt(sumfsq)*abs(e1-e2)/(2*sqrt(sumfppsq))
        sigdif=sqrt(sumfsq*(sige1**2+sige2**2)/(4*sumfppsq))
        call esym (iptgp,j,k,l,epsilon,m)
        epsilon=mlatt*epsilon
        difmax=sumfpp/sqrt(epsilon*sumfppsq)
      end if
c
c Acentric or centric reflection pair?
c
      call csym (iptgp,j,k,l,ic)
      if (ic.eq.0) then
c
c Acentric pair
c
        na=na+1
        m=2*m
        var=1
      else
c
c Centric pair
c
        nc=nc+1
        var=2
      end if
      write (io2) j,k,l,f1,sigf1,f2,sigf2,e1,sige1,e2,sige2,
     & ddiffe,sigdif,m,var,epsilon,difmax
      ndata=ndata+1
      go to 1
 9    continue
      endfile io2
      close (unit=io1,status='delete')
      write (ilp,929) n0,tmax,n1,xmin,n2,ymin,n3,smin,dmax,smax,dmin
 929  format (//1x,
     &'n0 = ',i6,'  reflection pairs accepted with  abs(Delta)/sig',
     &'ma                               .le. tmax,  where  tmax = ',
     &f4.1/1x,
     &'n1 = ',i6,'  of the n0  pairs accepted with  min{[E1/sigma(',
     &'E1)], [E2/sigma(E2)]}            .ge. xmin,  where  xmin = ',
     &f4.1/1x,
     &'n2 = ',i6,'  of the n1  pairs accepted with  abs(E1 - E2)/s',
     &'qrt[sigma(E1)**2 + sigma(E2)**2] .ge. ymin,  where  ymin = ',
     &f4.1/1x,
     &'n3 = ',i6,'  of the n2  pairs accepted with  smin .le. sin(',
     &'theta)/lambda .le. smax,  where'/1x,
     &'             smin = ',f7.4, ' reciprocal Angstroms      dmax',
     &' = ',f5.1,'  Angstroms'/1x,
     &'             smax = ',f7.4,'                           dmin',
     &' = ',f6.2)
      write (ilp,930) na,nc
 930  format (//1x,
     &'na = ',i6,'  acentric reflection pairs accepted'/1x,
     &'nc = ',i6,'  centric reflection pairs accepted')
c
c Re-scale to renormalize to <diffE**2> = 1.
c
      open (unit=io1,status='scratch',form='unformatted')
      call renorm (ndata,ginv,io2,io1,ilp,smin,smax)
      i=io1
      io1=io2
      io2=i
      close (unit=io2,status='delete')
c
c Compile renormalized diffE distribution statistics.
c
      write (ilp,600) atime,adate,atitle
      write (ilp,931) ndata,smin,1/(2*smin),smax,1/(2*smax)
 931  format (/1x,
     &'E-statistics for N accepted input data with smin .le. sin(the',
     &'ta)/lambda .le. smax, where'/1x,
     &'N    = ',i7,' reflections'/1x,
     &'smin = ',f7.4,' reciprocal Angstroms      dmax = ', f5.1,
     &'  Angstroms'/1x,
     &'smax = ',f7.4,'                           dmin = ', f6.2)
      call estats (io1,ilp,smin,smax,iptgp,mlatt,ginv)
c
c Apply output data selection tests.
c
      nrej1=0
      nrej2=0
      smin=9
      smax=0
      rewind io1
      open (unit=io2,status='scratch',form='unformatted')
      do i=1,ndata
        read (io1) j,k,l,f1,sigf1,f2,sigf2,e1,sige1,e2,sige2,
     &  ddiffe,sigdif,m,var,epsilon,difmax
        irej=0
        if (ddiffe/sigdif.lt.zmin) then
          irej=1
          nrej1=nrej1+1
        end if
        if ((ddiffe-difmax)/sigdif.gt.zmax) then
          irej=1
          nrej2=nrej2+1
        end if
        if (irej.eq.0) then
          write (io2) j,k,l,f1,sigf1,f2,sigf2,e1,sige1,e2,sige2,
     &    ddiffe,sigdif,m,var,epsilon,difmax
          s=sinthl(j,k,l,ginv)
          smin=min(smin,s)
          smax=max(smax,s)
        end if
      end do
      ndata=ndata-(nrej1+nrej2)
      close (unit=io1,status='delete')
      endfile io2
      i=io1
      io1=io2
      io2=i
      write (ilp,600) atime,adate,atitle
      write (ilp,932) nrej1,zmin,nrej2,zmax
 932  format (/1x,
     &'Nreject1 = ',i6,' re-scaled, renormalized data rejected bec',
     &'ause  diffE/sigma(diffE)              .lt. zmin = ',f4.1/1x,
     &'Nreject2 = ',i6,' re-scaled, renormalized data rejected bec',
     &'ause  (diffE - diffEmax)/sigma(diffE) .gt. zmax = ',f4.1)
      write (ilp,933) ndata,smin,1/(2*smin),smax,1/(2*smax)
 933  format (/1x,
     &'Naccept  = ',i6,' re-scaled, renormalized data accepted'/1x,
     &'s = sin(theta)/lambda                  d=1/(2*s)'/1x,
     &'smin = ',f7.4,' reciprocal Angstroms    dmax = ',f5.1,
     &'  Angstroms'/1x,
     &'smax = ',f7.4,'                         dmin = ',f6.2)
c
c Sort on diffE.
c
      rewind io1
      open (unit=io2,status='scratch',form='unformatted',
     & access='direct',recl=13)
      do i=1,ndata
        read (io1) j,k,l,f1,sigf1,f2,sigf2,e1,sige1,e2,sige2,ddiffe,
     &  sigdif
        write (io2,rec=i) j,k,l,f1,sigf1,f2,sigf2,e1,sige1,e2,sige2,
     &  ddiffe,sigdif
        data(i)=ddiffe
      end do
      call sort (ndata,data,index)
      close (unit=io1,status='delete')
c
c Record the accepted diffE values in the diffE-sorted output file
c "sort.diffe", and list the largest accepted values.
c
      open (unit=io1,file='sort.diffe',status='new')
      write (io1,934) atime,adate,atitle
 934  format (1x,'Program diffE.  ',a,', ',a,'.  ',a)
      write (io1,935)
 935  format (/1x,
     &'   h   k   l        F1  sig(F1)        F2  sig(F2)',
     &'  E1 sig(E1)  E2 sig(E2)  diffE sig(diffE)'/1x,
     &'   -   -   -        --  -------        --  -------',
     &'  -- -------  -- -------  ----- ----------')
      write (ilp,936)
 936  format (/1x, 'Largest accepted diffE values:')
      write (ilp,937)
 937  format (/1x,
     &'   h   k   l sin(th)/l     F1  sig(F1)        F2  sig(F2)',
     &'  E1 sig(E1)  E2 sig(E2) diffE sig(diffE) rank'/1x,
     &'   -   -   - ---------     --  -------        --  -------',
     &'  -- -------  -- ------- ----- ---------- ----')
      nlist=min(nlist,ndata)
      n=0
      do i=ndata,1,-1
        read  (io2,rec=index(i))
     &       j,k,l,f1,sigf1,f2,sigf2,e1,sige1,e2,sige2,ddiffe,sigdif
        write (io1,938)
     &       j,k,l,f1,sigf1,f2,sigf2,e1,sige1,e2,sige2,ddiffe,sigdif
 938    format (1x,3i4,2(e10.3,e9.2),4f6.3,2x,2f8.3)
        if (n.le.nlist.and.ddiffe.ge.1) then
          s=sinthl(j,k,l,ginv)
          write  (ilp,939) j,k,l,s,f1,sigf1,f2,sigf2,e1,sige1,
     &                     e2,sige2,ddiffe,sigdif,ndata-i+1
 939      format (1x,3i4,f7.4,2(e10.3,e9.2),4f6.3,2x,2f6.3,i6)
          n=n+1
        end if
      end do
      close (unit=io2,status='delete')
      close (unit=io1,status='keep')
      write (ilp,940) ndata,zmin,zmax
 940  format (/1x, 'The output reflection file is an ASCII file ',
     &'named "sort.diffe".'//1x,
     &'The file contains N = ',i6,' accepted reflection pairs with'/
     &1x,
     &'  diffE/sigma(diffE)              .ge. zmin = ',e10.3/1x,
     &'and'/1x,
     &'  (diffE - diffEmax)/sigma(diffE) .le. zmax = ',e10.3//1x,
     &'The file is sorted in order of decreasing diffE, and it has h',
     &'eader records labeling the reflection pair records:'/1x,
     &'h,k,l, F1,sig(F1), F2, sig(F2), E1, sig(E1), E2, sig(E2), dif',
     &'fE, sig(diffE)')
      stop 'Program diffE finis'
      end
c-----------------------------------------------------------------------
      subroutine csym (ptgp,h,k,l,c)
c
c Table of centrosymmetric rows and zones in non-centrosymmetric
c reciprocal lattice point 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 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
c c = 0 for acentric distributions
c c = 1 for centric distributions
c
      integer ptgp,h,c
      c=0
      go to (1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,
     & 23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42) ptgp
 1    return
 2    c=1
      return
 3    if (k.eq.0) c=1
      return
 4    if (h.eq.0.and.l.eq.0) c=1
      return
 5    c=1
      return
 6    if (h.eq.0.or.k.eq.0.or.l.eq.0) c=1
      return
 7    if (l.eq.0) c=1
      return
 8    c=1
      return
 9    if (l.eq.0) c=1
      return
 10   if (l.eq.0) c=1
      if (h.eq.0.and.k.eq.0) c=1
      return
 11   c=1
      return
 12   if (l.eq.0.or.k.eq.0.or.h.eq.0) c=1
      if (abs(k).eq.abs(h)) c=1
      return
 13   if (l.eq.0) c=1
      return
 14   if (l.eq.0.or.k.eq.0.or.h.eq.0) c=1
      return
 15   if (l.eq.0) c=1
      if (abs(k).eq.abs(h)) c=1
      return
 16   c=1
      return
 17   return
 18   c=1
      return
 19   if (k.eq.h.or.k.eq.-2*h) c=1
      return
 20   if (k.eq.0.or.h.eq.0) c=1
      return
 21   if ((k.eq.0.and.l.eq.0).or.(h.eq.0.and.l.eq.0)) c=1
      return
 22   return
 23   c=1
      return
 24   c=1
      return
 25   if (l.eq.0) c=1
      return
 26   if (h.eq.0.and.k.eq.0) c=1
      return
 27   c=1
      return
 28   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
      return
 29   if (l.eq.0) c=1
      return
 30   if (k.eq.h.or.k.eq.-2*h) c=1
      return
 31   if (h.eq.0.or.k.eq.0) c=1
      return
 32   c=1
      return
 33   return
 34   c=1
      return
 35   if (h.eq.k.or.k.eq.l.or.l.eq.h) c=1
      return
 36   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
      return
 37   c=1
      return
 38   if (h.eq.0.or.k.eq.0.or.l.eq.0) c=1
      return
 39   c=1
      return
 40   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
      return
 41   if (h.eq.0.or.k.eq.0.or.l.eq.0) c=1
      return
 42   c=1
      return
      end
c-----------------------------------------------------------------------
      subroutine datain (ifile,jfile,ilp,type,iptgp,tmax)
c
c Compile Delta(E) distribution statistics and histogram,
c optionally trim data from the distribution tails, and
c write a working data file.
c
      character type*3
      parameter (imax=50000,jmax=25)
      dimension x(imax),index(imax),xh(jmax),yh(jmax)
c
c Read dummy header records.
c
      rewind ifile
      do i=1,4
        read (ifile,'(a1)') dummy
      end do
c
c Count data and find  median(E1 - E2).
c
      ndata=0
      nc=0
      do i=1,imax
        read (ifile,*,end=5) j,k,l,f1,sigf1,f2,sigf2,e1,sige1,e2,sige2
        ndata=ndata+1
        if (type.eq.'SAS') then
          call csym (iptgp,j,k,l,ic)
          if (ic.eq.1) nc=nc+1
        end if
        x(i)=e1-e2
      end do
 5    continue
      call sort (ndata,x,index)
      xmedian=x(index(nint(0.5*ndata)))
c
c Compile  Delta(E) = E1 - E2 - median(E1 - E2)  histogram.
c
      xmean=0
      do i=1,ndata
        d=x(i)-xmedian
        xmean=xmean+d
        x(i)=d
      end do
      xmin=x(index(1))
      xmax=x(index(ndata))
      xmean=xmean/ndata
      do j=1,jmax
        xh(j)=0
        yh(j)=0
      end do
      j=0
      xj=xmin
      d=(xmax-xmin)/jmax
      do i=1,ndata
        xi=x(index(i))
        do while (xi.ge.xj)
          j=j+1
          xj=xj+d
        end do
        if (j.le.jmax) then
          xh(j)=xh(j)+xi
          yh(j)=yh(j)+1
        end if
      end do
      xh(jmax)=xh(jmax)+xi
      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
        xj=xmin
        do while (xj.lt.0)
          j=j+1
          xj=xj+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,941) ndata,xmedian,xmin,xmax,xmean
 941  format (//1x,
     &'Delta(E) = E1 - E2 - median(E1 - E2) Distribution Statistics'/
     &1x,
     &'------------------------------------------------------------'/
     &1x,
     &'Before trimming distribution tails from the data set:' / 1x,
     &'------'/1x,
     &'  N                    = ', i6,'  Delta(E) values' / 1x,
     &'  median (E1 - E2)     = ', f6.3 / 1x,
     &'  minimum Delta(E)     = ', f6.3 / 1x,
     &'  maximum              = ', f6.3 / 1x,
     &'  mean                 = ', f6.3)
c
c Store median(E1 - E2).
c
      dmedian=xmedian
c
c Find the median and mean of the nonzero abs[Delta(E)] values.
c
      n0=0
      xmean=0
      do i=1,ndata
        d=abs(x(i))
        if (d.eq.0) n0=n0+1
        xmean=xmean+d
        x(i)=d
      end do
      call sort (ndata,x,index)
      xmedian=x(index(nint(0.5*(ndata+n0))))
      xmean=xmean/(ndata-n0)
      write (ilp,942) xmedian,xmean 
 942  format (1x,
     &'  median abs[Delta(E)] = ', f6.3 / 1x,
     &'  mean                 = ', f6.3)
c
c Print Delta(E) histogram.
c
      write (ilp,943)
 943  format (/1x,
     &'Delta(E) distribution before trimming distribution tails' / 1x,
     &'--------------------------------------------------------' / 1x,
     &'     n  <Delta(E)>' / 1x, '     -  ----------')
      do j=1,jmax
        n=nint(50*yh(j)/yhmax)
        write (ilp,944) nint(yh(j)),xh(j),('X',i=1,n)
 944    format (1x,i6,e12.3,1x,50a1)
      end do
c
c Write a working data file after optional rejection of data pairs with
c abnormal Delta(E) values in the extreme tails of the distribution.
c
      sigma=1.25*xmedian
      dlimit=tmax*sigma
      nrej=0
      nrejn=0
      nrejp=0
      n=0
      xmean=0
      rewind ifile
      do i=1,4
        read (ifile,'(a1)') dummy
      end do
      do i=1,ndata
        read (ifile,*) ih,ik,il,f1,sigf1,f2,sigf2,e1,sige1,e2,sige2
c
c Re-center the (E1 - E2) values, aportioning a larger fraction of the
c origin sfift, median(E1 - E2), to the more uncertain of E1 and E2.
c
        e1=e1-dmedian*sige1/(sige1+sige2)
        e2=e2+dmedian*sige2/(sige1+sige2)
c
c Re-average centric SAS pairs.
c
        if (type.eq.'SAS') then
          call csym (iptgp,ih,ik,il,ic)
          if (ic.eq.1) then
            e1=0.5*(e1+e2)
            e2=e1
          end if
        end if
        d=e1-e2
        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,945)
 945                     format (/1x,
     &'Reflection pairs trimmed from the tails of the distribution:'/
     &1x,
     &'------------------------------------------------------------'/
     &1x,
     &'   h   k   l    E1 sigmaE1    E2 sigmaE2 (E1-E2)/sqrt(sigE1**',
     &'2+sigE2**2)'/1x,
     &'   -   -   -    -- -------    -- ------- --------------------',
     &'-----------')
          if (nrej.le.1000) write (ilp,946) ih,ik,il,e1,sige1,
     &                        e2,sige2,(e1-e2)/sqrt(sige1**2+sige2**2)
 946                        format (1x,3i4,2(2x,2f6.3),e10.2)
        else
          write (jfile) ih,ik,il,f1,sigf1,f2,sigf2,e1,sige1,e2,sige2
          n=n+1
          x(n)=d
          xmean=xmean+d
        end if
      end do
      close (unit=ifile,status='keep')
      endfile jfile
      i=ifile
      ifile=jfile
      jfile=i
c
c Recompile Delta(E) distribution statistics and histogram.
c
      if (nrej.eq.0) then
        write  (ilp,947) tmax
 947    format (/1x,
     &   'No data trimmed from distribution tails with  tmax = ',e10.3)
      else
        ndata=ndata-nrej
        call sort (ndata,x,index)
        xmin=x(index(1))
        xmax=x(index(ndata))
        xmedian=x(index(nint(0.5*ndata)))
        xmean=xmean/ndata
        do j=1,jmax
          xh(j)=0
          yh(j)=0
        end do
        j=0
        xj=xmin
        d=(xmax-xmin)/jmax
        do i=1,ndata
          xi=x(index(i))
          do while (xi.ge.xj)
            j=j+1
            xj=xj+d
          end do
          if (j.le.jmax) then
            xh(j)=xh(j)+xi
            yh(j)=yh(j)+1
          end if
        end do
        xh(jmax)=xh(jmax)+xi
        yh(jmax)=yh(jmax)+1
        if (type.eq.'SAS') then
          j=0
          xj=xmin
          do while (xj.lt.0)
            j=j+1
            xj=xj+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,948) tmax,nrej,nrejn,nrejp,ndata,xmin,xmax,xmean
 948    format (/1x,
     &'Delta(E) Distribution Statistics'/1x,
     &'--------------------------------'/1x,
     &'After trimming distribution tails:'/1x,
     &'-----'//1x,
     &'               maximum abs[Delta(E)]'/1x,
     &'  tmax  =  ------------------------- = ',e10.3/1x,
     &'           1.25*median abs[Delta(E)]'//1x,
     &'  Nrej                 = ', i6, '  rejected data pairs'/1x,
     &'  Nrejn                = ', i6, '  negative Delta(E)'/1x,
     &'  Nrejp                = ', i6, '  positive'//1x,
     &'  N                    = ', i6, '  remaining data pairs'/1x,
     &'  minimum Delta(E)     = ', f6.3 / 1x,
     &'  maximum              = ', f6.3 / 1x,
     &'  mean                 = ', f6.3)
        xmean=0
        do i=1,ndata
          d=abs(x(i))
          xmean=xmean+d
          x(i)=d
        end do
        call sort (ndata,x,index)
        xmedian=x(index(nint(0.5*(ndata+n0))))
        xmean=xmean/(ndata-n0)
        write (ilp,949) xmedian,xmean
 949    format (1x,
     &  '  median abs[Delta(E)] = ', f6.3 / 1x,
     &  '  mean                 = ', f6.3)
        write (ilp,950)
 950    format (/1x,
     &  'Delta(E) distribution after trimming distribution tails'/1x,
     &  '-------------------------------------------------------'/1x,
     &  '     n  <Delta(E)>'/1x,'     -  ----------')
        do j=1,jmax
          n=nint(50*yh(j)/yhmax)
          write (ilp,951) nint(yh(j)),xh(j),('X',i=1,n)
 951      format (1x,i6,e12.3,1x,50a1)
        end do
      end if
      return
      end
c-----------------------------------------------------------------------
      subroutine decode (ptgp,laue,syst,latt,mlatt,iptgp)
c
c Crystal system, lattice type and multiplicity, and Laue group and
c point group symmetries
c
c There are ten cases of alternative axes for seven of the 32
c crystallographic point groups (see subroutine esym), thus a total of
c 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 syst*12,latt*1,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 'Error in point group symbol.'
 1    iptgp=i
c
c Laue group
c
      if (i.ge. 1) laue='-1   '
      if (i.ge. 3) laue='2/M  '
      if (i.ge. 6) laue='MMM  '
      if (i.ge. 9) laue='4/M  '
      if (i.ge.12) laue='4/MMM'
      if (i.ge.17) laue='-3   '
      if (i.ge.19) laue='-31M '
      if (i.eq.20.or.i.eq.22.or.i.eq.24) laue='-3M1 '
      if (i.ge.25) laue='6/M  '
      if (i.ge.28) laue='6/MMM'
      if (i.ge.33) laue='R-3  '
      if (i.ge.35) laue='R-3M '
      if (i.ge.38) laue='M3   '
      if (i.ge.40) laue='M3M  '
c
c crystal system
c
      if (i.ge. 1) syst='Triclinic'
      if (i.ge. 3) syst='Monoclinic'
      if (i.ge. 6) syst='Orthorhombic'
      if (i.ge. 9) syst='Tetragonal'
      if (i.ge.17) syst='Trigonal'
      if (i.ge.25) syst='Hexagonal'
      if (i.ge.33) syst='Rhombohedral'
      if (i.ge.38) syst='Cubic'
c
c lattice centering multiplicity
c
      if (latt.eq.'P') mlatt=1
      if (latt.eq.'A') mlatt=2
      if (latt.eq.'B') mlatt=2
      if (latt.eq.'C') mlatt=2
      if (latt.eq.'F') mlatt=4
      if (latt.eq.'I') mlatt=2
      if (latt.eq.'R') mlatt=1
      if (latt.eq.'R'.and.syst.ne.'Rhombohedral') mlatt=3
      return
      end
c-----------------------------------------------------------------------
      function delta(i,j)
c
c Kronecker delta
c
      delta=0
      if (i.eq.j) delta=1
      return
      end
c-----------------------------------------------------------------------
      subroutine estats (ifile,ilp,smin,smax,iptgp,mlatt,ginv)
c
c E-statistics                  Theoretical Moments
c                           Acentric Centric Hypercentric
c
c a1 = <E>                     0.866   0.798   0.718
c a2 = <E**2>                  1.000   1.000   1.000 
c a3 = <E**3>                  1.329   1.596   1.916 
c a4 = <E**4>                  2.000   3.000   4.500 
c a5 = <E**5>                  3.323   6.383  12.260 
c a6 = <E**6>                  6.000  15.000  37.500 
c b1 = <abs(E**2 - 1)>         0.736   0.968   1.145 
c b2 = <(E**2 - 1)**2>         1.000   2.000   3.500 
c b3 = <(E**2 - 1)**3>         2.000   8.000  26.000 
c c3 = <abs(E**2 - 1)**3>      2.415   8.691  26.903
c
c f1  = fraction E .ge. 1      0.368   0.320
c f2                    2      0.018   0.050
c f3                    3      0.0001  0.003
c
c f05 = fraction E .lt. 0.5    0.22    0.38
c f15 =            .gt. 1.5    0.105   0.13
c
      dimension f1(3),f2(3),f3(3),ni(7),a1(7),a2(7),a3(7),a4(7),a5(7),
     & a6(7),b1(7),b2(7),b3(7),c3(7),ginv(3,3),h(3),si(10),nj(10),
     & sj(10),e1(10),e2(10),f05(10),f15(10)
      data f1,f2,f3,ni,a1,a2,a3,a4,a5,a6,b1,b2,b3,c3,nj,sj,e1,e2,f05,f15
     & /146*0/
      dimension w(50),x(50),y(50),z(50)
c
c Moments of the E-distribution
c
      rewind ifile
      nrec=0
      na=0
      nc=0
 1    read (ifile,end=9) j,k,l,(dummy,i=1,8),e,sige
      nrec=nrec+1
      s=sinthl(j,k,l,ginv)
      if (smin.le.s.and.s.le.smax) then
        call esym (iptgp,j,k,l,epsilon,m)
        call csym (iptgp,j,k,l,ic)
        if (ic.eq.0) m=2*m
        do 10 i=1,7
          if (i.eq.2.and.ic.ne.0) go to 10
          if (i.eq.3.and.ic.ne.1) go to 10
          if (i.eq.4.and.(j.eq.0.or.k.eq.0.or.l.eq.0)) go to 10
          if (i.eq.5.and.j.ne.0) go to 10
          if (i.eq.6.and.k.ne.0) go to 10
          if (i.eq.7.and.l.ne.0) go to 10
          ni(i)=ni(i)+m
          a1(i)=a1(i)+m*e
          a2(i)=a2(i)+m*e**2
          a3(i)=a3(i)+m*e**3
          a4(i)=a4(i)+m*e**4
          a5(i)=a5(i)+m*e**5
          a6(i)=a6(i)+m*e**6
          b1(i)=b1(i)+m*abs(e**2-1)
          b2(i)=b2(i)+m*(e**2-1)**2
          b3(i)=b3(i)+m*(e**2-1)**3
          c3(i)=c3(i)+m*abs(e**2-1)**3
          if (i.le.3) then
            if (e.ge.1) f1(i)=f1(i)+m
            if (e.ge.2) f2(i)=f2(i)+m
            if (e.ge.3) f3(i)=f3(i)+m
            if (e.lt.0.5) f05(i)=f05(i)+m
            if (e.gt.1.5) f15(i)=f15(i)+m
          end if
 10     continue
        if (ic.eq.0) na=na+1
        if (ic.eq.1) nc=nc+1
      end if
      go to 1
 9    continue
      do i=1,7
        if (ni(i).gt.0) then
          a1(i)=a1(i)/ni(i)
          a2(i)=a2(i)/ni(i)
          a3(i)=a3(i)/ni(i)
          a4(i)=a4(i)/ni(i)
          a5(i)=a5(i)/ni(i)
          a6(i)=a6(i)/ni(i)
          b1(i)=b1(i)/ni(i)
          b2(i)=b2(i)/ni(i)
          b3(i)=b3(i)/ni(i)
          c3(i)=c3(i)/ni(i)
        end if
      end do
      do i=1,3
        if (ni(i).gt.0) then
          f1(i)=f1(i)/ni(i)
          f2(i)=f2(i)/ni(i)
          f3(i)=f3(i)/ni(i)
          f05(i)=f05(i)/ni(i)
          f15(i)=f15(i)/ni(i)
        end if
      end do
      write (ilp,900)
      write (ilp,901) a1, 0.866,  0.798,  0.718
      write (ilp,902) a2, 1.000,  1.000,  1.000
      write (ilp,903) a3, 1.329,  1.596,  1.916
      write (ilp,904) a4, 2.000,  3.000,  4.500
      write (ilp,905) a5, 3.323,  6.383, 12.260
      write (ilp,906) a6, 6.000, 15.000, 37.500
      write (ilp,907) b1, 0.736,  0.968,  1.145
      write (ilp,908) b2, 1.000,  2.000,  3.500
      write (ilp,909) b3, 2.000,  8.000, 26.000
      write (ilp,910) c3, 2.415,  8.691, 26.903
      write (ilp,911) ni
      write (ilp,912) f1, 0.368,  0.320,
     &                f2, 0.018,  0.050,
     &                f3, 0.0001, 0.003,
     &                (f05(i),i=1,3), 0.22,  0.38,
     &                (f15(1),i=1,3), 0.105, 0.13
 900  format ('0                                           Experimenta',
     &'l Data Averages                            Theoretical Moments'/
     &'                       ------------------------------------------
     &--------------------------  -------------------------------'/' Ave
     &rage               All Data  Acentric   Centric       hkl       0k
     &l       h0l       hk0  Acentric   Centric Hypercentric'/
     &'                       --------  --------   -------       ---',
     &'       ---       ---       ---  --------   ------- ------------')
 901  format (' <E>                 ',10f10.3)
 902  format (' <E**2>              ',10f10.3)
 903  format (' <E**3>              ',10f10.3)
 904  format (' <E**4>              ',10f10.3)
 905  format (' <E**5>              ',10f10.3)
 906  format (' <E**6>              ',10f10.3)
 907  format (' <abs(E**2 - 1)>     ',10f10.3)
 908  format (' <(E**2 - 1)**2>     ',10f10.3)
 909  format (' <(E**2 - 1)**3>     ',10f10.3)
 910  format (' <abs(E**2 -1)**3>   ',10f10.3)
 911  format ('                       ----------------------------------
     &----------------------------------  ------------------------------
     &-'/' Mult.-Wtd. No.      ',7i10)
 912  format ('0                                        Experimental',
     &'              Theoretical'/'                                 ----
     &------------------------  ------------------'/
     &'                                 All Data  Acentric   Centric  Ac
     &entric   Centric'/'                                 --------  ----
     &----   -------  --------   -------'/
     &' Fraction of data with E .ge. 1',5f10.3/
     &'                              2',5f10.3/
     &'                              3',f11.4,f10.4,f9.3,f11.4,f9.3/
     &'                                 ----------------------------  --
     &----------------'/
     &' Fraction of data with E .lt. 0.5',f7.2,1x,4(f9.2,1x)/
     &'                       E .gt. 1.5',f7.2,1x,2(f9.2,1x),f10.3,f9.2/
     &'                                 ----------------------------  --
     &----------------')
c
c E-statistics in equal-volume shells of sin(theta)/lambda
c
      n=10
      do i=1,n
        si(i)=(float(i)/n)**(1/3.0)*smax
      end do
      rewind ifile
      do 2 irec=1,nrec
        read (ifile) j,k,l,(dummy,i=1,8),e,sige
        s=sinthl(j,k,l,ginv)
        if (smin.le.s.and.s.le.smax) then
          call esym (iptgp,j,k,l,epsilon,m)
          call csym (iptgp,j,k,l,ic)
          if (ic.eq.0) m=2*m
          do i=1,n
            if (s.le.si(i)) then
              nj(i)=nj(i)+m
              sj(i)=sj(i)+m*s
              e1(i)=e1(i)+m*e
              e2(i)=e2(i)+m*e**2
              if (e.lt.0.5) f05(i)=f05(i)+m
              if (e.gt.1.5) f15(i)=f15(i)+m
              go to 2
            end if
          end do
        end if
 2    continue
      do i=1,n
        if (nj(i).gt.0) then
          sj(i)=sj(i)/nj(i)
          e1(i)=e1(i)/nj(i)
          e2(i)=e2(i)/nj(i)
          f05(i)=f05(i)/nj(i)
          f15(i)=f15(i)/nj(i)
        end if
      end do
      write (ilp,920) (sj(i),e2(i),e1(i),f05(i),f15(i),nj(i),i=1,n)
 920  format ('0E-statistics in equal-volume shells of s = sin(theta)/la
     &mbda'/'0       <s>    <E**2>       <E> f(E < 0.5) f(E > 1.5) Mult.
     &-Wtd. No.'/'        ---    ------       --- ---------- ---------',
     &'- --------------'/(' ',3f10.3,2f10.2,i10))
c
c E-statistics in index parity groups
c
      do i=1,8
        nj(i)=0
        e1(i)=0
        e2(i)=0
        f05(i)=0
        f15(i)=0
      end do
      rewind ifile
      do irec=1,nrec
        read (ifile) j,k,l,(dummy,i=1,8),e,sige
        s=sinthl(j,k,l,ginv)
        if (smin.le.s.and.s.le.smax) then
          call esym (iptgp,j,k,l,epsilon,m)
          call csym (iptgp,j,k,l,ic)
          if (ic.eq.0) m=2*m
          call parity (i,j,k,l)
          nj(i)=nj(i)+m
          e1(i)=e1(i)+m*e
          e2(i)=e2(i)+m*e**2
          if (e.lt.0.5) f05(i)=f05(i)+m
          if (e.gt.1.5) f15(i)=f15(i)+m
        end if
      end do
      do i=1,8
        if (nj(i).gt.0) then
          e1(i)=e1(i)/nj(i)
          e2(i)=e2(i)/nj(i)
          f05(i)=f05(i)/nj(i)
          f15(i)=f15(i)/nj(i)
        end if
      end do
      write (ilp,931) (e2(i),e1(i),f05(i),f15(i),nj(i),i=1,8)
 931  format ('0E-statistics for index parity groups'/
     &'0    Parity    <E**2>       <E> f(E < 0.5) f(E > 1.5) Mult.-Wtd',
     &'. No.'/'     ------    ------       --- ---------- ---------- ---
     &-----------'/
     &'        eee',2f10.3,2f10.2,i10/
     &'        eeo',2f10.3,2f10.2,i10/
     &'        eoe',2f10.3,2f10.2,i10/
     &'        oee',2f10.3,2f10.2,i10/
     &'        eoo',2f10.3,2f10.2,i10/
     &'        oeo',2f10.3,2f10.2,i10/
     &'        ooe',2f10.3,2f10.2,i10/
     &'        ooo',2f10.3,2f10.2,i10)
c
c Plot E-distributions, P(E) vs. E, where 0 < P < 1 and 0 < E < 5.0.
c
      if (na.gt.0) then
c
c Acentric distribution
c
        n=min(50,nint(0.1*na))
        do i=1,n
          x(i)=0
          y(i)=0
        end do
        dx=5.0/n
        rewind ifile
        do 20 irec=1,nrec
          read (ifile) j,k,l,(dummy,i=1,8),e,sige
          s=sinthl(j,k,l,ginv)
          if (s.lt.smin.or.s.gt.smax) go to 20
          call csym (iptgp,j,k,l,ic)
          if (ic.eq.1) go to 20
          call esym (iptgp,j,k,l,epsilon,m)
          m=2*m
          do i=1,n
            if (e.le.i*dx) then
              x(i)=x(i)+m*e
              y(i)=y(i)+m
              go to 20
            end if
          end do
 20     continue
c
c Evaluate bin averages, and eliminate any empty bins.
c
        m=n
        n=0
        do i=1,m
          if (y(i).gt.0) then
            n=n+1
            x(n)=x(i)/y(i)
            y(n)=y(i)
          end if
        end do
c
c Ensure increasing x values.
c
        m=n
        n=1
        do i=2,m
          if (x(i).gt.x(i-1)) then
            n=n+1
            x(n)=x(i)
            y(n)=y(i)
          end if
        end do
c
c Normalize bin areas.
c
        a=y(1)*(x(2)-x(1))+y(n)*(x(n)-x(n-1))
        do i=2,n-1
          a=a+0.5*y(i)*(x(i+1)-x(i-1))
        end do
        do i=1,n
          y(i)=y(i)/a
        end do
        write (ilp,600)
        call plota (n,x,y,ilp)
      end if
 600  format ('1'//1x,10x,'Program diffE.  Acentric E-distribution'/)
      if (nc.gt.0) then
c
c Centric distribution
c
        n=min(50,nint(0.1*nc))
        do i=1,n
          x(i)=0
          y(i)=0
        end do
        dx=5.0/n
        rewind ifile
        do 21 irec=1,nrec
          read (ifile) j,k,l,(dummy,i=1,8),e,sige
          s=sinthl(j,k,l,ginv)
          if (s.lt.smin.or.s.gt.smax) go to 21
          call esym (iptgp,j,k,l,epsilon,m)
          call csym (iptgp,j,k,l,ic)
          if (ic.eq.0) go to 21
          do i=1,n
            if (e.le.i*dx) then
              x(i)=x(i)+m*e
              y(i)=y(i)+m
              go to 21
            end if
          end do
 21     continue
c
c Evaluate bin averages, and eliminate any empty bins.
c
        m=n
        n=0
        do i=1,m
          if (y(i).gt.0) then
            n=n+1
            x(n)=x(i)/y(i)
            y(n)=y(i)
          end if
        end do
c
c Ensure increasing x values.
c
        m=n
        n=1
        do i=2,m
          if (x(i).gt.x(i-1)) then
            n=n+1
            x(n)=x(i)
            y(n)=y(i)
          end if
        end do
c
c Normalize bin areas.
c
        a=y(1)*(x(2)-x(1))+y(n)*(x(n)-x(n-1))
        do i=2,n-1
          a=a+0.5*y(i)*(x(i+1)-x(i-1))
        end do
        do i=1,n
          y(i)=y(i)/a
        end do
        write (ilp,610)
        call plotc (n,x,y,ilp)
      end if
 610  format (1h1//1x,10x,'Program diffE.  Centric E-distribution'/)
c
c Plot <E**2> versus <sin(theta)/lambda>.
c
      n=min(50,nint(float(nrec)/50))
      do i=1,n
        w(i)=(float(i)/n)**(1/3.0)*smax
      end do
      do i=1,n
        x(i)=0
        y(i)=0
        z(i)=0
      end do
      rewind ifile
      do 3 irec=1,nrec
        read (ifile) j,k,l,(dummy,i=1,8),e,sige
        s=sinthl(j,k,l,ginv)
        if (s.lt.smin.or.s.gt.smax) go to 3
        call esym (iptgp,j,k,l,epsilon,m)
        call csym (iptgp,j,k,l,ic)
        if (ic.eq.0) m=2*m
        do i=1,n
          if (s.le.w(i)) then
            x(i)=x(i)+m*s
            y(i)=y(i)+m*e**2
            z(i)=z(i)+m
            go to 3
          end if
        end do
 3    continue
c
c Evaluate bin averages, and eliminate any empty bins.
c
      m=n
      n=0
      do i=1,m
        if (z(i).gt.0) then
          n=n+1
          x(n)=x(i)/z(i)
          y(n)=y(i)/z(i)
        end if
      end do
c
c Ensure increasing x values.
c
      m=n
      n=1
      do i=2,m
        if (x(i).gt.x(i-1)) then
          n=n+1
          x(n)=x(i)
          y(n)=y(i)
        end if
      end do        
      write (ilp,630)
      call plots (n,x,y,smax,ilp)
 630  format (1h1//1x,10x,'Program diffE.  <E**2> versus <sin(theta)/lamb
     &da> distribution'/)
      return
      end
c-----------------------------------------------------------------------
      subroutine esym (ptgp,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
      integer ptgp,h
      e=1
      go to (1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,
     & 23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42) ptgp
 1    m=1
      go to 99
 2    m=2
      go to 99
 3    m=2
      if (h.eq.0.and.l.eq.0) e=2
      go to 99
 4    m=2
      if (k.eq.0) e=2
      go to 99
 5    m=4
      if (k.eq.0) e=2
      if (h.eq.0.and.l.eq.0) e=2
      go to 99
 6    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
      go to 99
 7    m=4
      if (h.eq.0) e=2
      if (k.eq.0) e=2
      if (h.eq.0.and.k.eq.0) e=4
      go to 99
 8    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
      go to 99
 9    m=4
      if (h.eq.0.and.k.eq.0) e=4
      go to 99
 10   m=4
      if (h.eq.0.and.k.eq.0) e=2
      go to 99
 11   m=8
      if (l.eq.0) e=2
      if (h.eq.0.and.k.eq.0) e=4
      go to 99
 12   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
      go to 99
 13   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
      go to 99
 14   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
      go to 99
 15   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
      go to 99
 16   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
      go to 99
 17   m=3
      if (h.eq.0.and.k.eq.0) e=3
      go to 99
 18   m=6
      if (h.eq.0.and.k.eq.0) e=3
      go to 99
 19   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
      go to 99
 20   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
      go to 99
 21   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
      go to 99
 22   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
      go to 99
 23   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
      go to 99
 24   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
      go to 99
 25   m=6
      if (h.eq.0.and.k.eq.0) e=6
      go to 99
 26   m=6
      if (l.eq.0) e=2
      if (h.eq.0.and.k.eq.0) e=3
      go to 99
 27   m=12
      if (l.eq.0) e=2
      if (h.eq.0.and.k.eq.0) e=6
      go to 99
 28   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
      go to 99
 29   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
      go to 99
 30   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
      go to 99
 31   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
      go to 99
 32   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
      go to 99
 33   m=3
      if (k.eq.h.and.l.eq.h) e=3
      go to 99
 34   m=6
      if (k.eq.h.and.l.eq.h) e=3
      go to 99
 35   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
      go to 99
 36   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
      go to 99
 37   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
      go to 99
 38   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
      go to 99
 39   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
      go to 99
 40   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
      go to 99
 41   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
      go to 99
 42   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
 99   m=m/e
      return
      end
c-----------------------------------------------------------------------
      subroutine fcalc (n,x,a1,b1,a2,b2,a3,b3,a4,b4,c0,fp,fpp,s,
     & sumf,sumfsq,sumfpp,sumfppsq)
      dimension x(1),a1(1),b1(1),a2(1),b2(1),a3(1),b3(1),a4(1),b4(1),
     & c0(1),fp(1),fpp(1)
      double precision sumf,sumfsq,sumfpp,sumfppsq
      sumf=0
      sumfsq=0
      sumfpp=0
      sumfppsq=0
      do i=1,n
        xi=x(i)
        fi=a1(i)*exp(-b1(i)*s**2)
     &       +a2(i)*exp(-b2(i)*s**2)
     &          +a3(i)*exp(-b3(i)*s**2)
     &             +a4(i)*exp(-b4(i)*s**2)+c0(i)
        fpi=fp(i)
        fppi=fpp(i)
        fi=sqrt((fi+fpi)**2+fppi**2)
        sumf=sumf+xi*fi
        sumfsq=sumfsq+xi*fi**2
        sumfpp=sumfpp+xi*fppi
        sumfppsq=sumfppsq+xi*fppi**2
      end do
      return   
      end
c-----------------------------------------------------------------------
      subroutine fsqtof (fsq,sigfsq,f,sigf)
      if (fsq.le.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 ftofsq (f,sigf,fsq,sigfsq)
      if (f.lt.0) f=0
      if (sigf.lt.0) sigf=0
      fsq=f**2
      sigfsq=max(2*f*sigf,4*sigf**2)
      return
      end
c-----------------------------------------------------------------------
      subroutine metric (a,b,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)
      rad=pi/180
      deg=180/pi
      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 parity (i,j,k,l)
c
c 1 eee
c 2 eeo
c 3 eoe
c 4 oee
c 5 eoo
c 6 oeo
c 7 ooe
c 8 ooo
c
      p=mod(abs(j),2)
      q=mod(abs(k),2)
      r=mod(abs(l),2)
      if (p.eq.0.and.q.eq.0.and.r.eq.0) i=1
      if (p.eq.0.and.q.eq.0.and.r.eq.1) i=2
      if (p.eq.0.and.q.eq.1.and.r.eq.0) i=3
      if (p.eq.1.and.q.eq.0.and.r.eq.0) i=4
      if (p.eq.0.and.q.eq.1.and.r.eq.1) i=5
      if (p.eq.1.and.q.eq.0.and.r.eq.1) i=6
      if (p.eq.1.and.q.eq.1.and.r.eq.0) i=7
      if (p.eq.1.and.q.eq.1.and.r.eq.1) i=8
      return
      end
c-----------------------------------------------------------------------
      subroutine plota (n,x,y,ilp)
      parameter (nx=100,ny=50)
      character aa*1
      dimension aa(0:nx,0:ny),x(n),y(n)
c
c Blank out the plot array.
c
      do  ix=0,nx
      do  iy=0,ny
        aa(ix,iy)=' '
      end do
      end do
c
c Fill in the grid marks.
c
      do ix=0,nx
        aa(ix, 0)='.'
        aa(ix,ny)='.'
      end do
      do ix=0,nx,10
        aa(ix, 0)='+'
        aa(ix,ny)='+'
      end do
      do iy=0,ny
        aa( 0,iy)='.'
        aa(nx,iy)='.'
      end do
      do iy=0,ny,10
        aa( 0,iy)='+'
        aa(nx,iy)='+'
      end do
c
c Set the plot limits, 0 < E < 5.0 and 0 < p(E) < 1.
c
      xmin=0
      xmax=5
      ymin=0
      ymax=1
      xdel=xmax-xmin
      ydel=ymax-ymin
c
c Fill in the distribution function.
c
      dx=xdel/(2*nx)
      do i=0,(2*nx)
        xi=i*dx
        yi=2*xi*exp(-xi**2)
        ix=nint(nx*(xi-xmin)/xdel)
        iy=nint(ny*(ymax-yi)/ydel)
        ix=max(ix,0)
        ix=min(ix,nx)
        iy=max(iy,0)
        iy=min(iy,ny)
        aa(ix,iy)='*'
      end do
c
c Fill in the experimental data points.
c
      do i=1,n
        ix=nint(nx*(x(i)-xmin)/xdel)
        iy=nint(ny*(ymax-y(i))/ydel)
        ix=max(ix,0)
        ix=min(ix,nx)
        iy=max(iy,0)
        iy=min(iy,ny)
        aa(ix,iy)='O'
      end do
c
c Print the plot array.
c
      dy=(ymax-ymin)/ny
      do iy=0,ny
        if (mod(iy,10).eq.0) then
          write (ilp,100) ymax-iy*dy,(aa(ix,iy),ix=0,nx)
        else
          write (ilp,101)            (aa(ix,iy),ix=0,nx)
        end if
      end do
      dx=(xmax-xmin)/10
      write (ilp,102) (xmin+ix*dx,ix=0,10)
  100 format (' ',f10.1,101a1)
  101 format (' ',10x,  101a1)
  102 format (' ',1x,11f10.1)
      write (ilp,'(/1x,10x,''Acentric Distribution:  p(E) = 2*E*exp(-E**
     &2),  "*" Theoretical,  "o" Experimental'')')
      return
      end
c-----------------------------------------------------------------------
      subroutine plotc (n,x,y,ilp)
      parameter (nx=100,ny=50)
      character aa*1
      dimension aa(0:nx,0:ny),x(n),y(n)
c
c Blank out the plot array.
c
      do  ix=0,nx
      do  iy=0,ny
        aa(ix,iy)=' '
      end do
      end do
c
c Fill in the grid marks.
c
      do ix=0,nx
        aa(ix, 0)='.'
        aa(ix,ny)='.'
      end do
      do ix=0,nx,10
        aa(ix, 0)='+'
        aa(ix,ny)='+'
      end do
      do iy=0,ny
        aa( 0,iy)='.'
        aa(nx,iy)='.'
      end do
      do iy=0,ny,10
        aa( 0,iy)='+'
        aa(nx,iy)='+'
      end do
c
c Set the plot limits, 0 < E < 5.0 and 0 < p(E) < 1.
c
      xmin=0
      xmax=5
      ymin=0
      ymax=1
      xdel=xmax-xmin
      ydel=ymax-ymin
c
c Fill in the distribution function.
c
      dx=xdel/(2*nx)
      c=sqrt(2/3.141593)
      do i=0,(2*nx)
        xi=i*dx
        yi=c*exp(-0.5*xi**2)
        ix=nint(nx*(xi-xmin)/xdel)
        iy=nint(ny*(ymax-yi)/ydel)
        ix=max(ix,0)
        ix=min(ix,nx)
        iy=max(iy,0)
        iy=min(iy,ny)
        aa(ix,iy)='*'
      end do
c
c Fill in the experimental data points.
c
      do i=1,n
        ix=nint(nx*(x(i)-xmin)/xdel)
        iy=nint(ny*(ymax-y(i))/ydel)
        ix=max(ix,0)
        ix=min(ix,nx)
        iy=max(iy,0)
        iy=min(iy,ny)
        aa(ix,iy)='O'
      end do
c
c Print the plot array.
c
      dy=(ymax-ymin)/ny
      do iy=0,ny
        if (mod(iy,10).eq.0) then
          write (ilp,100) ymax-iy*dy,(aa(ix,iy),ix=0,nx)
        else
          write (ilp,101)            (aa(ix,iy),ix=0,nx)
        end if
      end do
      dx=(xmax-xmin)/10
      write (ilp,102) (xmin+ix*dx,ix=0,10)
  100 format (' ',f10.1,101a1)
  101 format (' ',10x,  101a1)
  102 format (' ',1x,11f10.1)
      write (ilp,'(/1x,10x,''Centric Distribution:  p(E) = sqrt(2/pi)*ex
     &p(-0.5*E**2),  "*" Theoretical,  "O" Experimental'')')
      return
      end
c-----------------------------------------------------------------------
      subroutine plots (n,x,y,smax,ilp)
      parameter (nx=100,ny=50)
      character aa*1
      dimension aa(0:nx,0:ny),x(n),y(n),ypp(50)
c
c Blank out the plot array.
c
      do  ix=0,nx
      do  iy=0,ny
        aa(ix,iy)=' '
      end do
      end do
c
c Fill in the grid marks.
c
      do ix=0,nx
        aa(ix, 0)='.'
        aa(ix,ny)='.'
      end do
      do ix=0,nx,10
        aa(ix, 0)='+'
        aa(ix,ny)='+'
      end do
      do iy=0,ny
        aa( 0,iy)='.'
        aa(nx,iy)='.'
      end do
      do iy=0,ny,10
        aa( 0,iy)='+'
        aa(nx,iy)='+'
      end do
c
c Set the plot limits, 0 < sin(theta)/lambda < smax and 0 < E(s) < 2.5.
c
      xmin=0
      xmax=smax
      ymin=0
      ymax=2.5
      xdel=xmax-xmin
      ydel=ymax-ymin
c
c Mark E**2 = 1.
c
      do ix=0,nx
        aa(ix,30)='.'
      end do
c
c Fill in the spline-interpolated data points.
c
      call spline (n,x,y,ypp)
      dx=(x(n)-x(1))/(2*nx)
      do i=0,(2*nx)
        xi=x(1)+i*dx
        call splint (n,x,y,ypp,xi,yi)
        ix=nint(nx*(xi-xmin)/xdel)
        iy=nint(ny*(ymax-yi)/ydel)
        ix=max(ix,0)
        ix=min(ix,nx)
        iy=max(iy,0)
        iy=min(iy,ny)
        aa(ix,iy)='*'
      end do
c
c Fill in the experimental data points.
c
      do i=1,n
        ix=nint(nx*(x(i)-xmin)/xdel)
        iy=nint(ny*(ymax-y(i))/ydel)
        ix=max(ix,0)
        ix=min(ix,nx)
        iy=max(iy,0)
        iy=min(iy,ny)
        aa(ix,iy)='O'
      end do
c
c Print the plot array.
c
      dy=(ymax-ymin)/ny
      do iy=0,ny
        if (mod(iy,10).eq.0) then
          write (ilp,100) ymax-iy*dy,(aa(ix,iy),ix=0,nx)
        else
          write (ilp,101)            (aa(ix,iy),ix=0,nx)
        end if
      end do
      dx=(xmax-xmin)/10
      write (ilp,102) (xmin+ix*dx,ix=0,10)
  100 format (' ',f10.1,101a1)
  101 format (' ',10x,  101a1)
  102 format (' ',1x,11f10.4)
      return
      end
c-----------------------------------------------------------------------
      subroutine renorm (ndata,ginv,ifile,jfile,ilp,smin,smax)
c
c Global re-scaling via a logarithmically linearized least-squares fit
c of the coefficients q0, q1, qnd q2 of an empirical renormalization
c function,
c   q = q0*exp(q1*s**2 + q2*s**4) ,
c in which
c   s = sin(theta)/lambda .
c The fit minimizes
c   chi**2 = sum w*[diffE**2 - q**2]**2
c so that
c   <(diffE/q)**2> = 1.
c
      dimension ginv(3,3),dummy(8)
      dimension a(3,3),ainv(3,3),b(3),c(3)
      real*8 sumdsq,summ,sumw,sumx,sumy,sumx2,sumx3,sumx4,sumy2,sumxy,
     & sumx2y,sumq,sumv,sumqsq
      smin=9
      smax=0
      sumdsq=0
      summ=0
      n=0
      sumw=0
      sumx=0
      sumy=0
      sumx2=0
      sumx3=0
      sumx4=0
      sumy2=0
      sumxy=0
      sumx2y=0
      rewind ifile
      do i=1,ndata
        read (ifile) j,k,l,dummy,ddiffe,sigdif,m,var
        s=sinthl(j,k,l,ginv)
        smin=min(smin,s)
        smax=max(smax,s)
        sumdsq=sumdsq+m*ddiffe**2
        summ=summ+m
        if (ddiffe.gt.0) then
          x=s**2
          y=ddiffe**2
          sigy=2*ddiffe*sigdif
c
c Logarithmically linearize the least-squares problem.
c
          w=m/((sigy/y)**2+var)
          y=log(y)
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
      write (ilp,960) sumdsq/summ,ndata,smin,1/(2*smin),smax,
     &                1/(2*smax)
 960  format (//1x,
     &'Before renormalization:'/1x,
     &'-----------------------'/1x,
     &'<diffE**2> = ',e10.3/1x,
     &'Ndata      = ',i10,' reflection pairs'/1x,
     &'s = sin(theta)/lambda                  d=1/(2*s)'/1x,
     &'smin = ',f7.4,' reciprocal Angstroms    dmax = ',f5.1,
     &'  Angstroms'/1x,
     &'smax = ',f7.4,'                         dmin = ',f6.2)
     &
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)
      q0=c(1)
      q1=c(2)
      q2=c(3)
c
c chisq = sum w*(y - q)**2
c       = sum w*y**2 + sum w*q**2 - 2*sum w*q*y
c
c q     = q0 + q1*x + q2*x**2
c
c q**2  = q0**2 + (q1*x + q2*x**2)**2 + 2*q0*(q1*x + q2*x**2)
c       = q0**2 + q1**2*x**2 + q2**2*x**4 + 2*q1*q2*x**3 
c                                     + 2*q0*q1*x + 2*q0*q2*x**2
c
      chisq=sumy2+q0**2*sumw+q1**2*sumx2+q2**2*sumx4
     &              +2*(q0*q1*sumx+q0*q2*sumx2+q1*q2*sumx3)
     &                -2*(q0*sumy+q1*sumxy+q2*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 q2 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(q2).lt.sqrt(v2)) then
        det=sumw*sumx2-sumx*sumx
        q0=(sumx2*sumy-sumx*sumxy)/det
        q1=(-sumx*sumy+sumw*sumxy)/det
        chisq=sumy2+q0**2*sumw+q1**2*sumx2
     &                -2*(q0*sumy+q1*sumxy-q0*q1*sumx)
        zsq=chisq/(n-2)
        v0=zsq*sumx2/det
        v1=zsq*sumw/det
        cov01=zsq*(-sumx/det)
        q2=0
        v2=0
        cov02=0
        cov12=0
c
c If the q1 coefficient of the term linear in s**2 is not statistically
c significant, resort to a constant scaling function.
c
        if (abs(q1).lt.sqrt(v1)) then
          q0=sumy/sumw
          q1=0
          chisq=sumy2+q0**2*sumw-2*q0*sumy
          zsq=chisq/(n-1)
          v0=zsq/sumw
          v1=0
          cov01=0
        end if
      end if
c
c Re-scale the data, and compile re-scaling statistics.
c
      sumdsq=0
      summ=0
      sumq=0
      sumv=0
      sumqsq=0
      qmin=+1e9
      qmax=-1e9
      rewind ifile
      rewind jfile
      do i=1,ndata
        read (ifile) j,k,l,dummy,ddiffe,sigdif,m,var,epsilon,difmax
        s=sinthl(j,k,l,ginv)
        q=q0+q1*s**2+q2*s**4
        v=v0+v1*s**4+v2*s**8+2*(cov01*s**2+cov02*s**4+cov12*s**6)
c
c Convert from log fitting scale to antilog data scale.
c
        q=exp(q)
        v=q**2*v
c
c Convert from diffE**2 to diffE scale.
c
        v=v/(4*q)
        q=sqrt(q)
        sigdif=sqrt((ddiffe/q)**2*((sigdif/ddiffe)**2+(v/q**2)))
        ddiffe=ddiffe/q
        sumdsq=sumdsq+m*ddiffe**2
        summ=summ+m
        sumq=sumq+m*q
        sumv=sumv+m*v
        sumqsq=sumqsq+m*q**2
        if (q.lt.qmin) then
          qmin=q
          sigmin=sqrt(v)
        end if
        if (q.gt.qmax) then
          qmax=q
          sigmax=sqrt(v)
        end if
        write (jfile) j,k,l,dummy,ddiffe,sigdif,m,var,epsilon,difmax
      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 with the unit-weighted average, <diffE**2> = 1.
c
      q=1/sqrt(sumdsq/summ)
      n1=0
      n2=0
      n3=0
      rewind ifile
      rewind jfile
      do i=1,ndata
        read (jfile) j,k,l,dummy,ddiffe,sigdif,m,var,epsilon,difmax
        ddiffe=q*ddiffe
        sigdif=q*sigdif
        write (ifile) j,k,l,dummy,ddiffe,sigdif,m,var,epsilon,difmax
        if (ddiffe.ge.  sigdif) n1=n1+1
        if (ddiffe.ge.2*sigdif) n2=n2+1
        if (ddiffe.ge.3*sigdif) n3=n3+1
      end do
c
c Convert the renormalization, re-scaling parameters and statistics from
c the logarithmic diffE**2 scale to the direct diffE scale.
c
      q0=sqrt(exp(q0))/q
      q1=q1/2
      q2=q2/2
      write (ilp,971) q0,q1,q2,sqrt(abs(q2))
 971  format (//1x,
     &'Renormalization Re-scaling:'/1x,
     &'---------------------------'/1x,
     &'Global re-scaling via a logarithmically linearized least-squa',
     &'res fit'/1x,
     &'of the coefficients q0, q1, qnd q2 of an empirical renormaliz',
     &'ation'/1x,
     &'function,'/1x,
     &'  q = q0*exp(q1*s**2 + q2*s**4) ,'/1x,
     &'in which'/1x,
     &'  s = sin(theta)/lambda .'//1x,
     &'The fit minimizes'/1x,
     &'  chi**2 = sum w*(diffE**2 - q**2)**2 ,'/1x,
     &'so that'/1x,
     &'  <(diffE/q)**2> = 1 .'//1x,
     &'q0 = ',e10.3/1x,
     &'q1 = ',e10.3/1x,
     &'q2 = ',e10.3,'    sqrt(abs(q2)) = ',e10.3)
      qmin=qmin/q
      qmax=qmax/q
      sigmin=sigmin/q
      sigmax=sigmax/q
      sumq=sumq/q
      sumv=sumv/q**2
      sumqsq=sumqsq/q**2
      write (ilp,972) qmin,sigmin,qmax,sigmax,sumq/summ,
     &           sqrt(sumv/summ),sqrt((sumqsq/summ)-(sumq/summ)**2)
 972  format (/1x,
     &'Renormalization Statistics:'/1x,
     &'---------------------------'/1x,
     &'qmin = ',e10.3,'  sigma(qmin) = ',e10.3/1x,
     &'qmax = ',e10.3,'  sigma(qmax) = ',e10.3//1x,
     &'<q>                   = ',e10.3/1x,
     &'<sigma(q)**2>**(1/2)  = ',e10.3/1x,
     &'<(q - <q>)**2>**(1/2) = ',e10.3)
      write (ilp,973) ndata,n1,n2,n3
 973  format (/1x,
     &'N0 = ',i6,' renormalized pairs'/1x,
     &'N1 = ',i6,' renormalized pairs with diffE .ge. 1*sigma(diff',
     &'E)'/1x,
     &'N2 = ', i6,'                                    2'/1x,
     &'N3 = ', i6,'                                    3')
      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 adapted 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
        indext=index(i)
        t=data(indext)
      else
        indext=index(j)
        t=data(indext)
        index(j)=index(1)
        j=j-1
        if (j.eq.1) then
          index(1)=indext
          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)=indext
      go to 1
      end
c-----------------------------------------------------------------------
      subroutine spline (n,x,y,ypp)
c
c Natural cubic spline fit to a tabulated function y(i) = y(x(i)), i =
c 1, 2,..., n, with x(1) .lt. x(2) .lt. ... .lt. x(n).  returns the
c tabulated second derivatives y'' = d2y/dx2 = ypp(i), i = 1, 2,..., n,
c with ypp(1) = ypp(n) = 0.
c
c Fortran code adapted from William H. Press, Brian P. Flannery, Saul A.
c Teukolsky, and William T. Vetterling (1986).  Numerical Recipies:  The
c Art of Scientific Computing, pp. 86-89.  Cambridge, England:
c Cambridge University Press.
c
c Must have nmax .ge. n.
c
      parameter (nmax=1000)
      dimension x(n),y(n),ypp(n),u(nmax)
      ypp(1)=0
      ypp(n)=0
      u(1)=0
      u(n)=0
      do i=2,n-1
        p=(x(i)-x(i-1))/(x(i+1)-x(i-1))
        q=p*ypp(i-1)+2
        ypp(i)=(p-1)/q
        u(i)=(6*((y(i+1)-y(i))/(x(i+1)-x(i))-(y(i)-y(i-1))/
     &   (x(i)-x(i-1)))/(x(i+1)-x(i-1))-p*u(i-1))/q
      end do
      do i=n-1,1,-1
        ypp(i)=ypp(i)*ypp(i+1)+u(i)
      end do
      return
      end
c-----------------------------------------------------------------------
      subroutine splint (n,x,y,ypp,xi,yi)
c
c Cubic spline interpolation of a tabulated function y(i) = y(x(i)), i =
c 1, 2,..., n, using the tabulated second derivatives ypp(i) from
c subroutine spline.
c
c Fortran code adapted from William H. Press, Brian P. Flannery, Saul A.
c Teukolsky, and William T. Vetterling (1986).  Numerical Recipies:  The
c Art of Scientific Computing, pp. 86-89.  Cambridge, England:
c Cambridge University Press.
c
      dimension x(n),y(n),ypp(n)
      if (xi.le.x(1)) then
        yi=y(1)
        return
      end if
      if (xi.ge.x(n)) then
        yi=y(n)
        return
      end if
      i1=1
      i2=n
      do while (i2-i1.gt.1)
        i=(i2+i1)/2
        if (x(i).gt.xi) then
          i2=i
        else
          i1=i
        end if
      end do
      h=x(i2)-x(i1)
      a=(x(i2)-xi)/h
      b=(xi-x(i1))/h
      yi=a*y(i1)+b*y(i2)+((a**3-a)*ypp(i1)+(b**3-b)*ypp(i2))*h**2/6
      return
      end
c-----------------------------------------------------------------------