program denzox
c
c Robert H. Blessing
c Hauptman-Woodward Institute
c 73 High St.
c Buffalo, NY 14203-1196, USA
c Tel: 716-856-9600, ext. 335
c e-mail:  blessing@hwi.buffalo.edu
c
c From a concatenated set of "denzo.x" files either:
c (0) Get fully recorded reflections; or
c (1) Scale both fully and partially recorded reflections and sum the
c     scaled partials.
c
      character atime*8,adate*9,title*80,file1*80,file2*80,filek*80
      real*8 fsqfull,sigfull,fsqsump,sigsump
c
c Maximum number of frames
c
      parameter (mframe=2000)
      dimension phii(mframe),scalek(mframe),sigmak(mframe)
      data scalek,sigmak /mframe*1, mframe*0/
c
c Maximum number of partial reflections
c
      parameter (mdata=2000000)
      integer data
      dimension data(mdata),index(mdata)
c
c Number of frames spanned by partial refnlections,
c n = 1, 2,..., 10 or more
c
      dimension nn(10)
      data nn /10*0/
      dimension ub(3,3),ginv(3,3),h(3),xyz0(3),r(3,3),xyz(3),u0(3),u1(3)
      data deg,rad /57.29578, 0.01745329/
      call time (atime)
      call date (adate)
      io1=10
      io2=20
      io3=30
      io5=50
      io6=60
c
c Read and echo the "denzox.dat" control data input file.
c
      open (io5,file='denzox.dat',status='old')
      open (io6,file='denzox.lp',status='new')
      data title/' '/, itype/1/, filek/' '/, file1/'frames.x'/,
     & file2/'data.denzox'/
      read (io5,'(a)') title
      write (io6,'(/1x,
     &''Program denzox.  '',a,'', '',a//1x,a//1x,
     &''From a concatenated set of "denzo.x" files either:''/1x,
     &''(itype = 0)  Get fully recorded reflections; or''/1x,
     &''(itype = 1)  Scale both fully and partially recorded reflecti'',
     &''ons and''/1x,
     &''             sum the scaled partials.'')') atime,adate,title
      read (io5,*) itype
      write (io6,'(/1x,''itype = '',i1)') itype
      if (itype.ne.0.and.itype.ne.1) itype=0
      if (itype.eq.1) then
        read (io5,'(a)') filek
        open (io1,file=filek,status='old',err=15)
        write (io6,'(/1x,
     &''Interframe scale factors file:''/1x,a//1x,
     &''iframe   scalek(i) sigmak(i)''/1x,
     &''------   --------- ---------'')') filek
        do i=1,mframe
          read (io1,*,err=20,end=25)       j,scalek(j),sigmak(j)
          write (io6,'(1x,i6,e12.4,e9.2)') j,scalek(j),sigmak(j)
        end do
 25     close (io1,status='keep')
        go to 11
 20     close (io1,status='keep')
        stop 'Stop.  Error while reading scale factors file.'
 15     close (io1,status='keep')
        write (io6,'(/1x,
     &''Scale factors equal to unity assumed for all frames.'')')
 11     continue
      end if
      read (io5,'(a)') file1
      read (io5,'(a)') file2
      write (io6,'(/1x,
     &''Input file of concatenated "Denzo.x" files:''/1x,a)') file1
      write (io6,'(/1x,
     &''Output file (ih, ik, il, Fsq, sigmaFsq, iframe):''/1x,a)') file2
      open (io1,file=file1,status='old')
      open (io2,file=file2,status='new')
      open (io3,status='scratch',form='unformatted',access='direct',
     & recl=6)
      iframe=1
      nrefln=0
      noverf=0
      negsig=0
      nfull=0
      npart=0
      nsump=0
      fsqfull=0
      sigfull=0
      fsqsump=0
      sigsump=0
      imin=+999
      imax=-999
      jmin=+999
      jmax=-999
      kmin=+999
      kmax=-999
      lmin=+999
      lmax=-999
      iend=0
 1    call readx (io6,io1,iend,j,k,l,fsq,sigfsq,iframe,ipart,phii,
     & nrefln,noverf,negsig,ub,ginv,flambda)
      if (iend.ne.0) go to 9
c
c Scale fully and partially recorded reflections.
c
      sigfsq=sqrt((fsq*sigmak(iframe))**2+(scalek(iframe)*sigfsq)**2)
      fsq=scalek(iframe)*fsq
      if (ipart.eq.0) then
        nfull=nfull+1
        fsqfull=fsqfull+abs(fsq)
        sigfull=sigfull+sigfsq
        write (io2,1000) j,k,l,fsq,sigfsq,iframe
      else
        npart=npart+1
        if (itype.eq.1) then
          imin=min(imin,iframe)
          imax=max(imax,iframe)
          jmin=min(jmin,j)
          jmax=max(jmax,j)
          kmin=min(kmin,k)
          kmax=max(kmax,k)
          lmin=min(lmin,l)
          lmax=max(lmax,l)
          write (io3,rec=npart) j,k,l,fsq,sigfsq,iframe
        end if
      end if
 1000 format (1x,3i5,2e15.7,i5)
      go to 1
 9    continue
      close (io1,status='keep')
      nframe=iframe-1
      write (io6,'(/1x,
     &''Nframes = '',i8/1x,
     &''Nreflns = '',i8/1x,
     &''Noverfl = '',i8/1x,
     &''NegSig  = '',i8/1x,
     &''Nparts  = '',i8/1x,
     &''Nfulls  = '',i8)') nframe,nrefln,noverf,negsig,npart,nfull
      if (itype.eq.0.or.npart.lt.2) go to 5
      if (npart.gt.mdata) then
        write (io6,'(/1x,
     &''Nparts = '',i8,'' exceeds Mparts = '',i8/1x,
     &''Only Mparts partial reflections can be sorted.'')') npart,mdata
        npart=mdata
      end if
c
c Sort partials on packed h, k, l, iframe.
c
      ni=imax-imin+1
      nl=lmax-lmin+1
      nk=kmax-kmin+1
      do i=1,npart
        read (io3,rec=i) j,k,l,fsq,sigfsq,iframe
        data(i)=(j-jmin)*nk*nl*ni+(k-kmin)*nl*ni+(l-lmin)*ni
     &            +(iframe-imin)
      end do
      call sort (npart,data,index)
c
c Sum scaled partials.
c
      read (io3,rec=index(1)) ih,ik,il,fsq,sigfsq,iframe
      n=1
      p=fsq
      q=sigfsq**2
      w=max(fsq,sigfsq)
      x=w
      y=w*iframe
      do i=2,npart
        read (io3,rec=index(i)) jh,jk,jl,fsq,sigfsq,jframe
c
c Summable partials must have identical indices and be on consecutive
c frames.
c
        if (jh.eq.ih.and.jk.eq.ik.and.jl.eq.il.and.jframe.eq.iframe+1)
     &  then
c
c Depending on crystal mosaicity, oscillation range, and hkl coordinates,
c a reflection might span more than two consecutive frames.
c
          n=n+1
          p=p+fsq
          q=q+sigfsq**2
          w=max(fsq,sigfsq)
          x=x+w
          y=y+w*jframe
        else
          n=min(n,10)
          nn(n)=nn(n)+1
          if (n.ge.2) then
            nsump=nsump+1
            q=sqrt(q)
            fsqsump=fsqsump+abs(p)
            sigsump=sigsump+q
c
c Assign summed partials the intensity-weighted average frame number.
c
            write (io2,1000) ih,ik,il,p,q,nint(y/x)
          end if
          n=1
          p=fsq
          q=sigfsq**2
          w=max(fsq,sigfsq)
          x=w
          y=w*jframe
          ih=jh
          ik=jk
          il=jl
        end if
        iframe=jframe
      end do
      write (io6,'(/1x,
     &''Nsump   = '',i8,'' summed partial reflns.''/1x,
     &''Ntotal  = '',i8,'' full + summed partial reflns.'')')
     & nsump,(nfull+nsump)
      write (io6,'(/1x,
     &''Number of  Number of''/1x,
     &'' Frames     Summed''/1x,
     &''Spanned    Partials''/1x,
     &''---------  ---------'')')
      write (io6,'(1x,'' 1'',8x,i10,'' rejected single partials'')')
     & nn(1)
      do n=2,9
        write (io6,'(1x,i2,8x,i10)') n,nn(n)
      end do
      write (io6,'(1x,''10 or more'',i10)') nn(10)
      write (io6,'(/1x,
     &''<sigma(Fsq)>/<Fsq> = '',e10.3,'' for Nfull = '',i8,
     &'' full reflns.''/1x,
     &''<sigma(Fsq)>/<Fsq> = '',e10.3,'' for Nsump = '',i8,
     &'' summed partial reflns.'')')
     & sigfull/fsqfull,nfull,sigsump/fsqsump,nsump
 5    continue
      stop
      end
c-----------------------------------------------------------------------
      subroutine metric (a,b,gb,v)
c
c Reciprocal lattice parameters and direct and reciprocal lattice metric
c 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 mm (a,b,c)
c
c Square matrix-matrix multiplication
c
c a(i,k)*b(k,j) = c(i,j)
c
      dimension a(3,3),b(3,3),c(3,3)
      do i=1,3
      do j=1,3
        c(i,j)=0
      do k=1,3
        c(i,j)=c(i,j)+a(i,k)*b(k,j)
      end do
      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 readx (io6,ifile,iend,ih,ik,il,fsq,sigfsq,iframe,ipart,
     & phii,nrefln,noverf,negsig,ub,ginv,flambda)
c
c Read and interpret a concatenated set of denzo.x files.
c
      save
      dimension phii(1),w(3,3),wa(3,3),rphinv(3,3),rphiu(3,3)
      dimension a(3,3),suma(3,3),asq(3,3),u(3,3),sumu(3,3),usq(3,3),
     & at(3,3),ata(3,3),ut(3,3),utu(3,3),row(3),col(3),
     & b(3,3),ub(3,3),uh(3,3),
     & spin(3),sumspin(3),spinsq(3),vert(3),sumvert(3),vertsq(3),
     & sumcell(6),cellsq(6),rot(3),sumrot(3),rotsq(3)
      dimension esda(3,3),esdu(3,3),esdspn(3),esdvrt(3),esdcll(6),
     & esdrot(3)
      real*8 suma,asq,sumu,usq,sumspin,spinsq,sumvert,vertsq,
     & sumcell,cellsq,sumrot,rotsq,sumpol,polsq,sumwl,wlsq,sumeta,etasq,
     & sumdist,distsq
      dimension cell(6),rlcell(6),ginv(3,3)
      character title*80,record*80,word*9,flag*2,c1*1,c2*1,c3*1,c4*1,
     & c5*1,c6*1
      dimension ihkl(3,13)
      data ihkl /1,0,0, 0,1,0, 0,0,1, 1,1,0, -1,1,0, 1,0,1, -1,0,1,
     & 0,1,1, 0,-1,1, 1,1,1, -1,1,1, 1,-1,1, -1,-1,1/
      data rad /0.01745329/
      if (nrefln.gt.0) go to 5
c
c On first entry, initialize frame variables to be accumulated.
c
      n=0
      do i=1,3
      do j=1,3
        suma(i,j)=0
        sumu(i,j)=0
        asq(i,j)=0
        usq(i,j)=0
        rphinv(i,j)=0
      end do
      end do
      rphinv(3,3)=1
      sumpol=0
      polsq=0
      sumwl=0
      wlsq=0
      sumeta=0
      etasq=0
      do i=1,3
        sumspin(i)=0
        sumvert(i)=0
        sumcell(i)=0
        sumcell(i+3)=0
        sumrot(i)=0
        spinsq(i)=0
        vertsq(i)=0
        cellsq(i)=0
        cellsq(i+3)=0
        rotsq(i)=0
      end do
      sumdist=0
      distsq=0
      read (ifile,'(a)',end=9) title
      write (io6,'(/1x,''First frame title:''/1x,a)') title
      write (io6,'(/1x,
     &''Frame oscillation limits and reflection populations:''//1x,
     &''iframe    phi1    phi2 abs(1-2) avg(phi) nfull   npart''/1x,
     &''------    ----    ---- -------- -------- -----   -----'')')
      rewind ifile
 1    continue
c
c Start reading a new frame.
c
      nfull=0
      npart=0
c
c Read title line.
c
      read (ifile,'(a)',end=9) record
      if (record.ne.title) then
        title=record
        write (io6,'(/1x,''New frame title:''/1x,a)') title
      end if
c
c Read three lines of orientation matrix rows.
c
      do i=1,3
        read (ifile,'(a)',end=9) record
c
c Check for decimal points in (3f15.8,3f10.6) fields.
c
        read (record,'(3(6x,a1,8x),3(3x,a1,6x))') c1,c2,c3,c4,c5,c6
        if (c1.ne.'.'.or.c2.ne.'.'.or.c3.ne.'.'
     &  .or.c4.ne.'.'.or.c5.ne.'.'.or.c6.ne.'.') go to 1
        read (record,'(3f15,3f10)') (a(i,j),j=1,3),(u(i,j),j=1,3)     
      end do
c
c Read one line of oscillation frame parameters.
c
      read (ifile,'(a)',end=9) record
 5    continue
c
c Read a reflection record.
c
      read (ifile,'(a)',end=9) record
c
c Reflection records for each frame are terminated by a (3i4) record
c (i,j,k) with i = 999 and j, k = integration box dimensions.
c
      read (record,'(3i4,a2)',err=1) i,j,k,flag
      if (i.eq.999.or.(flag.ne.' 0'.and.flag.ne.' 1')) go to 6
c
c For reflection records, check for leading blanks in (i4) h, k, and l
c fields.
c
      read (record,'(3(a1,3x))') c1,c2,c3
      if (c1.ne.' '.or.c2.ne.' '.or.c3.ne.' ') go to 5
c
c Check for full/part flag = 0/1.
c
      read (record,'(3(4x),a2)') flag
      if (flag.ne.' 0'.and.flag.ne.' 1') go to 5
c
c Count reflection records.
c
      nrefln=nrefln+1
c
c Check for "overflow" records by checking for decimal points in the
c (f8.1) fsq and (f6.1) sigfsq fields.
c
      read (record,'(3(4x),2x,(6x,a1,1x),8x,7x,(4x,a1,1x))') c1,c2
      if (c1.ne.'.'.or.c2.ne.'.') noverf=noverf+1
c
c Read a full or partial reflection record.
c
      read (record,'(3i4,i2,f8.0,8x,7x,f6.0)',err=1) ih,ik,il,ipart,fsq,
     & sigfsq
c
c Reject any measurements with nonpositive sigma(fsq).
c
      if (sigfsq.le.0) then
        negsig=negsig+1
        go to 5
      end if
c
c Count fully and partially recorded reflections.
c
      if (ipart.eq.0) then
        nfull=nfull+1
      else
        npart=npart+1
      end if
c
c Return data for one reflection to the calling program.
c
      return
 6    continue
c
c After all of a frame's reflection records have been read, accumulate
c the frame's orientation data.
c
      n=n+1
      phii(n)=0.5*(phi1+phi2)
c
c The Denzo A-matrix is about constant from frame to frame, but the
c Denzo U-matrix changes with phi-rotation.
c
      cosphi=cos(phii(n)*rad)
      sinphi=sin(phii(n)*rad)
      rphinv(1,1)= cosphi
      rphinv(1,2)=+sinphi
      rphinv(2,1)=-sinphi
      rphinv(2,2)= cosphi
      do i=1,3
      do j=1,3
        rphiu(i,j)=u(i,j)
      end do
      end do
      call mm (rphinv,rphiu,u)
      do i=1,3
      do j=1,3
        suma(i,j)=suma(i,j)+a(i,j)
        sumu(i,j)=sumu(i,j)+u(i,j)
        asq(i,j)=asq(i,j)+a(i,j)**2
        usq(i,j)=usq(i,j)+u(i,j)**2
      end do
      end do
c
c Read to end of the frame, and accumulate the other frame variables.
c
      iend=0
      do while (iend.eq.0)
        read (ifile,'(a)',end=9) record
        read (record,'(a)') word
c
c The last record in each "frame.x" file gives the "crossfire"
c parameters.
c
        if (word.eq.'crossfire') iend=1
c
c Read monochromator polarization ratio.
c
        if (word.eq.'monochrom') then
          read (record,'(13x,f8)') pol
          sumpol=sumpol+pol
          polsq=polsq+pol**2
        end if
c
c Record oscillation range phi1 and phi2 and reflection population for
c each frame.
c
        if (word.eq.'oscillati') then
          read (record,'(17x,f8,4x,f8)') phi1,phi2
          write (io6,'(1x,i6,4f8.2,2i8)')
     &    iframe,phi1,phi2,abs(phi1-phi2),0.5*(phi1+phi2),nfull,npart
        end if
c
c Read wavelength.
c       
        if (word.eq.'wavelengt') then
          read (record,'(10x,f8)') wl
          sumwl=sumwl+wl
          wlsq=wlsq+wl**2
        end if
c
c Read mosaicity.
c       
        if (word.eq.'mosaicity') then
          read (record,'(9x,f8)') eta
          sumeta=sumeta+eta
          etasq=etasq+eta**2
        end if
c
c Read spindle axis [hkl] and "vertical" axis [hkl].
c
        if (word.eq.'spindle a') then
          read (record,'(12x,3f4,14x,3f4)') (spin(i),i=1,3),(vert(i),
     &    i=1,3)
          do i=1,3
            sumspin(i)=sumspin(i)+spin(i)
            sumvert(i)=sumvert(i)+vert(i)
            spinsq(i)=spinsq(i)+spin(i)**2
            vertsq(i)=vertsq(i)+vert(i)**2
          end do
        end if
c
c Read unit cell dimensions.
c
        if (word.eq.'unit cell') then
          read (record,'(9x,6f9)') (cell(i),i=1,6)
          do i=1,6
            sumcell(i)=sumcell(i)+cell(i)
            cellsq(i)=cellsq(i)+cell(i)**2
          end do
        end if
c
c Read crystal rotation angles.
c
        if (word.eq.'crystal r') then
          read (record,'(7x,3(5x,f9))') (rot(i),i=1,3)
          do i=1,3
            sumrot(i)=sumrot(i)+rot(i)
            rotsq(i)=rotsq(i)+rot(i)**2
          end do
        end if
c
c Read crystal-to-detector distance.
c
        if (word.eq.'distance ') then
          read (record,'(8x,f8)') dist
          sumdist=sumdist+dist
          distsq=distsq+dist**2
        end if
      end do
      iframe=iframe+1
      iend=0
      go to 1
 9    iend=1
c
c After all frames have been read, average the data collection
c variables.
c
      write (io6,'(/1x,
     &''Means and root-mean-square deviations from the means''/1x,
     &''====================================================''/1x,
     &''for the data collection variables''/1x,
     &''=================================''//1x,
     &''n = '',i6,'' frames'')') n
      f=sqrt(float(n)/(n-1))
      pol=sumpol/n
      polsq=polsq/n
      esdpol=f*sqrt(abs(polsq-pol**2))
      write (io6,'(/1x,
     &''polarization ratio and rmsd, r = [I(sigma) - I(pi)]/[I(sigma)'',
     &'' + I(pi)]''/1x,
     &''-------------------------------------------------------------'',
     &''--------''/1x,f7.4/1x,f7.4)') pol,esdpol
      wl=sumwl/n
      wlsq=wlsq/n
      esdwl=f*sqrt(abs(wlsq-wl**2))
      write (io6,'(/1x,
     &''wavelength and rmsd''/1x,
     &''-------------------''/1x,f8.5/1x,f8.5)') wl,esdwl
      flambda=wl
      eta=sumeta/n
      etasq=etasq/n
      esdeta=f*sqrt(abs(etasq-eta**2))
      write (io6,'(/1x,
     &''mosaicity and rmsd''/1x,
     &''------------------''/1x,f7.4/1x,f7.4)') eta,esdeta
      dist=sumdist/n
      distsq=distsq/n
      esddis=f*sqrt(abs(distsq-dist**2))
      write (io6,'(/1x,
     &''crystal-to-detector distance and rmsd''/1x,
     &''-------------------------------------''/1x,f9.3/1x,f9.3)') dist,
     & esddis
      do i=1,3
        spin(i)=sumspin(i)/n
        vert(i)=sumvert(i)/n
        spinsq(i)=spinsq(i)/n
        vertsq(i)=vertsq(i)/n
      end do
      do i=1,3
        esdspn(i)=f*sqrt(abs(spinsq(i)-spin(i)**2))
        esdvrt(i)=f*sqrt(abs(vertsq(i)-vert(i)**2))
      end do
      write (io6,'(/1x,
     &''spindle axis:  [hkl] = ['',3i2,'' ] rmsds = '',3e9.2/1x,
     &''------------'')') (nint(spin(i)),i=1,3),esdspn
      write (io6,'(/1x,
     &''vertical axis: [hkl] = ['',3i2,'' ] rmsds = '',3e9.2/1x,
     &''-------------'')') (nint(vert(i)),i=1,3),esdvrt
      do i=1,3
        rot(i)=sumrot(i)/n
        rotsq(i)=rotsq(i)/n
      end do
      do i=1,3
        esdrot(i)=f*sqrt(abs(rotsq(i)-rot(i)**2))
      end do
      write (io6,'(/1x,
     &''crystal rotations and rmsds''/1x,
     &''---------------------------''/1x,
     &''    phi-x    phi-y    phi-z''/1x,3f9.3/1x,3f9.3)') rot,esdrot
      do i=1,6
        cell(i)=sumcell(i)/n
        cellsq(i)=cellsq(i)/n
      end do
      do i=1,6
        esdcll(i)=f*sqrt(abs(cellsq(i)-cell(i)**2))
      end do
      write (io6,'(/1x,
     &''averaged lattice parameters and rmsds''/1x,
     &''-------------------------------------''/1x,
     &''        a        b        c    alpha     beta    gamma''/1x,
     & 6f9.3/1x,6f9.3)') cell,esdcll
      call metric (cell,rlcell,ginv,vcell)
      write (io6,'(/1x,
     &''  Vcell = '',e12.4/1x,
     &''1/Vcell = '',e12.4)') vcell,1/vcell
      write (io6,'(/1x,
     &''inverse metric matrix''/1x,
     &''---------------------''/(1x,3f12.8))') ((ginv(i,j),j=1,3),i=1,3)
c
c Rossman and Blow (1962) reciprocal space orthogonalization matrix
c
c Rossman, Michael G., and Blow, D.M (1962).  The Detection of Sub-Units
c Within the Crystallographic Asymmetric Unit.  Acta Cryst. 15, 24-31.
c                                               ----------- ==
c
      omega=acos((cos(cell(5)*rad)-cos(cell(4)*rad)*cos(cell(6)*rad))/
     &              (sin(cell(4)*rad)*sin(cell(6)*rad)))
      a11= 1/(cell(1)*sin(cell(6)*rad)*sin(omega))
      a12= 0
      a13= 0
      a21= 1/(cell(2)*tan(cell(4)*rad)*tan(omega))
     &      -1/(cell(2)*tan(cell(6)*rad)*sin(omega))
      a22= 1/cell(2)
      a23=-1/(cell(2)*tan(cell(4)*rad))
      a31=-1/(cell(3)*sin(cell(4)*rad)*tan(omega))
      a32= 0
      a33= 1/(cell(3)*sin(cell(4)*rad))
      write (io6,'(/1x,
     &''Rossman and Blow (1962) reciprocal space orthogonalization ma'',
     &''trix alpha''/1x,
     &''-------------------------------------------------------------'',
     &''----------'',/1x,
     &''x-axis parallel to b cross c,''/1x,
     &''y-axis parallel to b, and''/1x,
     &''z-axis parallel to x cross y.''//1x,
     &''omega = acos{[cos(beta) - cos(alpha)*cos(gamma)]/[sin(alpha)*'',
     &''sin(gamma)]}''//1x,
     &''alpha(1,1) =  1/[a*sin(gamma)*sin(omega)]''/1x,
     &''alpha(1,2) =  0''/1x,
     &''alpha(1,3) =  0''/1x,
     &''alpha(2,1) =  1/[b*tan(alpha)*tan(omega)] - 1/[b*tan(gamma)*s'',
     &''in(omega)]''/1x,
     &''alpha(2,2) =  1/b''/1x,
     &''alpga(2,3) = -1/[b*tan(alpha)]''/1x,
     &''alpha(3,1) = -1/[c*sin(alpha)*tan(omega)]''/1x,
     &''alpha(3,2) =  0''/1x,
     &''alpha(3,3) =  1/[c*sin(alpha)]''//(1x,3f12.8))')
     & a11,a12,a13,a21,a22,a23,a31,a32,a33
c
c Busing and Levy (1967) reciprocal space orthogonalization matrix
c
c Busing, William R., and Levy, Henri, A. (1967).  Angle Calculations
c for 3- and 4-Circle X-Ray and Neutron Diffractometers.  Acta Cryst.
c                                                         -----------
c 22, 457-464.
c ==
c
      write (io6,'(/1x,
     &''Busing-Levy (1967) reciprocal space orthogonalization matrix B''
     &/1x,
     &''--------------------------------------------------------------''
     &/1x,
     &''x-axis parallel to a_star,''/1x,
     &''y-axis in the plane of a_star and b_star, and''/1x,
     &''z-axis perpendicular to that plane.''//1x,
     &''    [ a_star b_star*cos(gamma_star)  c_star*cos(beta_star)   '',
     &''         ]''/1x,
     &''B = [ 0      b_star*sin(gamma_star) -c_star*sin(beta_star)*co'',
     &''s(alpha) ]''/1x,
     &''    [ 0      0                       1/c                     '',
     &''         ]''//(1x,3f12.8))')
     & rlcell(1), rlcell(2)*cos(rlcell(6)*rad), rlcell(3)*cos(rlcell(5)*
     & rad),
     & 0.0, rlcell(2)*sin(rlcell(6)*rad), -rlcell(3)*sin(rlcell(5)*rad)*
     & cos(cell(4)*rad),
     & 0.0, 0.0, 1/cell(3)
c
c Denzo A-matrix
c
      do i=1,3
      do j=1,3
        a(i,j)=suma(i,j)/n
        asq(i,j)=asq(i,j)/n
      end do
      end do
      do i=1,3
      do j=1,3
        esda(i,j)=f*sqrt(abs(asq(i,j)-a(i,j)**2))
      end do
      end do
      write (io6,'(/1x,
     &''averaged Denzo A-matrix                 rmsds''/1x,
     &''-----------------------                 -----'')')
      do i=1,3
        write (io6,'(1x,3e12.4,5x,3e9.2)')
     &  (a(i,j),j=1,3),(esda(i,j),j=1,3)
      end do
      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)
      write (io6,'(/1x,''det(A) = '',e12.4)') d
      do i=1,3
      do j=1,3
        at(i,j)=a(j,i)
      end do
      end do
      call mm (at,a,ata)
      write (io6,'(/1x,
     &''(A-transpose)*A''/1x,
     &''---------------''/(1x,3f12.8))') ((ata(i,j),j=1,3),i=1,3)
c
c Denzo U-matrix
c
      do i=1,3
      do j=1,3
        u(i,j)=sumu(i,j)/n
        usq(i,j)=usq(i,j)/n
      end do
      end do
      do i=1,3
      do j=1,3
        esdu(i,j)=f*sqrt(abs(usq(i,j)-u(i,j)**2))
      end do
      end do
      write (io6,'(/1x,
     &''averaged Denzo (Phi-inverse)*U matrix   rmsds''/1x,
     &''-------------------------------------   -----'')')
      do i=1,3
        write (io6,'(1x,3e12.4,5x,3e9.2)')
     &  (u(i,j),j=1,3),(esdu(i,j),j=1,3)
      end do
      d=u(1,1)*u(2,2)*u(3,3)+u(1,2)*u(2,3)*u(3,1)+u(1,3)*u(2,1)*u(3,2)
     &  -u(3,1)*u(2,2)*u(1,3)-u(3,2)*u(2,3)*u(1,1)-u(3,3)*u(2,1)*u(1,2)
      write (io6,'(/1x,''det(U) = '',e12.4)') d
      do i=1,3
      do j=1,3
        ut(i,j)=u(j,i)
      end do
      end do
      call mm (ut,u,utu)
      write (io6,'(/1x,
     &''(U-transpose)*U''/1x,
     &''---------------''/(1x,3f12.8))') ((utu(i,j),j=1,3),i=1,3)
      do i=1,3
        col(i)=0
        row(i)=0
      do j=1,3
        col(i)=col(i)+u(j,i)**2
        row(i)=row(i)+u(i,j)**2
      end do
      end do
      write (io6,'(/1x,
     &''column sum[u(i,j)**2] = '',3e12.4/1x
     &''row    sum[u(i,j)**2] = '',3e12.4)')
     & (col(i),i=1,3),(row(i),i=1,3)
      call mm (ut,a,b)
      write (io6,'(/1x,
     &''(U-transpose)*A = (U-transpose)*U*A0 = A0''/1x,
     &''-----------------------------------------''/(1x,3f12.8))')
     & ((b(i,j),j=1,3),i=1,3)
c
c Transform the Denzo A-matrix from Denzo to Rossman XYZ-axes.
c
c [x]   [a11 a12 a13] [h]
c [y] = [a21 a22 a23] [k]
c [z]   [a31 a32 a33] [l]
c
c phi       ^             x(D) .                y(R) .
c spindle   .                  .                     .
c           .                  .                     .
c (right-   .                  .                     .   x(R)
c handed    .                  .                     .  .
c rotation) .                  .                     . .
c           .                  .                     ..
c           . . . . . . >      . . . . . . z(D)      . . . . . . z(R)
c                  X-rays     .
c                            .
c                           .
c                          y(D)
c
c [x(R)]   [ 0 -1  0 ] [x(D)]
c [y(R)] = [ 1  0  0 ] [y(D)]
c [z(R)]   [ 0  0  1 ] [z(D)]
c
c Rossman, Michael G. (1979).  Processing Oscillation Diffraction Data
c for Very Large Unit Cells with an Automatic Convolution Technique and
c Profile Fitting.  J. Appl. Cryst., 12, 225-238.
c                   ---------------  ==
c
      do i=1,3
      do j=1,3
        w(i,j)=0
      end do
      end do
      w(1,2)=-1
      w(2,1)=+1
      w(3,3)=+1
c
c Apparently, Denzo prints the A-matrix referred to Rossman XYZ axes,
c not to the Denzo spindle/beam XYZ axes described in the Denzo "HKL
c Manual."  So, make the transformation matrix w(3,3) a unit matrix.
c
      do i=1,3
      do j=1,3
        w(i,j)=0
      end do
      end do
      w(1,1)=1
      w(2,2)=1
      w(3,3)=1
      call mm (w,a,wa)
c
c Transform from Rossman to Wonacott XYZ-axes.
c
c phi       ^             y(R) .                z(W) .
c spindle   .                  .                     .
c           .                  .                     .
c (right-   .                  .   x(R)              .   y(W)
c handed    .                  .  .                  .  .
c rotation) .                  . .                   . .
c           .                  ..                    ..
c           . . . . . . >      . . . . . . z(R)      . . . . . . x(W)
c                  X-rays
c
c [x(W)]   [ 0  0  1 ] [x(R)]
c [y(W)] = [ 1  0  0 ] [y(R)]
c [z(W)]   [ 0  1  0 ] [z(R)]
c
c Wonacott, A.J. (1977).  Geometry of the Rotation Method. In The
c                                                             ---
c Rotation Method in Crystallography, pp. 75-103, edited by U.W. Arndt
c ----------------------------------
c and A.J. Wonacott.  Amsterdam:  North-Holland Pub. Co.
c
      do i=1,3
      do j=1,3
        w(i,j)=0
      end do
      end do
      w(1,3)=1
      w(2,1)=1
      w(3,2)=1
      call mm (w,wa,ub)
      write (io6,'(/1x,
     &''Busing-Levy orientation matrix UB referred to Wonacott XYZ ax'',
     &''es''/1x,
     &''-------------------------------------------------------------'',
     &''--''/1x,
     &''[x]   [ub(1,1) ub(1,2) ub(1,3)] [h]''/1x,
     &''[y] = [ub(2,1) ub(2,2) ub(2,3)] [k]''/1x,
     &''[z]   [ub(3,1) ub(3,2) ub(3,3)] [l]''//
     & (1x,3e12.4))') ((ub(i,j),j=1,3),i=1,3)
      d=ub(1,1)*ub(2,2)*ub(3,3)
     &    +ub(1,2)*ub(2,3)*ub(3,1)
     &      +ub(1,3)*ub(2,1)*ub(3,2)
     &        -ub(3,1)*ub(2,2)*ub(1,3)
     &          -ub(3,2)*ub(2,3)*ub(1,1)
     &            -ub(3,3)*ub(2,1)*ub(1,2)
      write (io6,'(/1x,''det(UB) = '',e12.4)') d
c
c Try transformation by rotations about each of the Cartesian axes
c (x, y, z), i.e., the permutations of paiwise (yz, xz, xy) reversals of
c the axial directions.
c
      do i=1,3
      do j=1,3
        w(i,j)=0
      end do
      end do
      w(1,1)=+1
      w(2,2)=-1
      w(3,3)=-1
      call mm (w,ub,ut)
      write (io6,'(/1x,''+X,-Y,-Z''//(1x,3e12.4))')
     & ((ut(i,j),j=1,3),i=1,3)
      w(1,1)=-1
      w(2,2)=+1
      w(3,3)=-1
      call mm (w,ub,ut)
      write (io6,'(/1x,''-X,+Y,-Z''//(1x,3e12.4))')
     & ((ut(i,j),j=1,3),i=1,3)
      w(1,1)=-1
      w(2,2)=-1
      w(3,3)=+1
      call mm (w,ub,ut)
      write (io6,'(/1x,''-X,-Y,+Z''//(1x,3e12.4))')
     & ((ut(i,j),j=1,3),i=1,3)
c
c Transform from Wonacott to Hamilton XYZ-axes, and transpose to obtain
c the Hamilton U-matrix.
c
c phi       ^             z(W) .               -z(H) .
c spindle   .                  .                     .
c           .                  .                     .
c (right-   .                  .   y(W)              .   y(H)
c handed    .                  .  .                  .  .
c rotation) .                  . .                   . .
c           .                  ..                    ..
c           . . . . . . >      . . . . . . x(W)      . . . . . . -x(H)
c                  X-rays
c
c [x(H)]   [-1  0  0 ] [x(W)]
c [y(H)] = [ 0  1  0 ] [y(W)]
c [z(H)]   [ 0  0 -1 ] [z(W)]
c
c Hamilton, Walter C. (1974).  Angle Settings for Four-Circle
c Diffractometers. In International Tables for X-Ray Crystallography,
c                     ----------------------------------------------
c Vol. IV, pp. 273-284, edited by James A. Ibers and Walter C. Hamilton.
c Birmingham, England:  The Kynoch Press.
c
      do i=1,3
      do j=1,3
        w(i,j)=0
      end do
      end do
      w(1,1)=-1
      w(2,2)= 1
      w(3,3)=-1
      call mm (w,ub,uh)
      do i=1,3-1
      do j=i+1,3
        x=uh(i,j)
        uh(i,j)=uh(j,i)
        uh(j,i)=x
      end do
      end do
      write (io6,'(/1x,
     &''Hamilton orientation matrix UH''/1x,
     &''------------------------------''/(1x,3e12.4))')
     & ((uh(i,j),j=1,3),i=1,3)
c
c List pricipal zone axis setting angles.
c
      write (io6,'(/1x,
     &''Hamilton setting angles for principal zone axes''//1x,
     &'' [ h   k   l ]         omega     chi     phi''/1x,
     &''   -   -   -           -----     ---     ---'')')
      do i=1,13
        ih=ihkl(1,i)
        ik=ihkl(2,i)
        il=ihkl(3,i)
        call setang (ih,ik,il,uh,ginv,flambda,twoth,omega,chi,phi)
        write (io6,'(1x,3i4,8x,3f8.2)') ih,ik,il,omega,chi,phi
      end do
c
c Find zone axis nearest the spindle axis.
c
      xmin=90
      do ih=-3,+3
      do ik=-3,+3
      do il= 0,+3
        if (ih.ne.0.or.ik.ne.0.or.il.ne.0) then
          call setang (ih,ik,il,uh,ginv,flambda,twoth,omega,chi,phi)
          x=90-abs(chi)
          if (x.lt.xmin) then
            xmin=x
            jh=ih
            jk=ik
            jl=il
            omegaj=omega
            chij=chi
            phij=phi
          end if
        end if
      end do
      end do
      end do
      write (io6,'(/1x,3i4,8x,3f8.2)') jh,jk,jl,omegaj,chij,phij
c
c Compute orientation matrix from lattice parameters and Denzo crystal
c rotations.
c
c      call ucalc (io6,cell,spin,vert,rot(1),rot(2),rot(3),u)
c
c 2 Jan. '97.  Comment out the call to ucalc, because the Denzo "HKL
c Manual" gives confusing and apparently incorrect descriptions of the 
c relationships among the laboratory xyz-axes; the instrumental spindle,
c "vertical", and X-ray axes; and the phi-x, phi-y, and phi-z crystal
c rotations.
c
      return
      end
c-----------------------------------------------------------------------
      subroutine setang (ih,ik,il,u,ginv,flambda,twoth,omega,chi,phi)
c
c Setting angles for Hamilton (1974) International Tables diffractometer
c axes in bisecting position, psi = omega = 0, at positive two-theta
c
      dimension h(3),u(3,3),x(3),ginv(3,3)
      data pi,deg /3.141593, 57.29578/
      h(1)=ih
      h(2)=ik
      h(3)=il
      call vm (h,u,x)
      dstar=2*sinthl(ih,ik,il,ginv)
      do i=1,3
        x(i)=x(i)/dstar
        if (x(i).lt.-1) x(i)=-1
        if (x(i).gt.+1) x(i)=+1
      end do
      twoth=2*asin(0.5*dstar*flambda)
      omega=0
      chi=asin(-x(3))
c
c Choose chi in the interval -pi/2 .le. chi .le. + pi/2.
c
      chi=mod(chi,2*pi)
      if (chi.lt.-pi) chi=chi+2*pi
      if (chi.gt.+pi) chi=chi-2*pi
      if (chi.lt.-pi/2) chi=-pi-chi
      if (chi.gt.+pi/2) chi=+pi-chi
      coschi=cos(chi)
      cosphi=x(2)/coschi
      sinphi=x(1)/coschi
      phi=atan2(sinphi,cosphi)
      twoth=twoth*deg
      omega=omega*deg
      phi=phi*deg
      chi=chi*deg
      twoth=a180(twoth)
      omega=a180(omega)
      phi=a180(phi)
      chi=a180(chi)
      return
      end
c-----------------------------------------------------------------------
      function a180(a)
      a=mod(a,360.0)
      if (a.le.-180) a=a+360
      if (a.gt.+180) a=a-360
      a180=a
      return
      end
c-----------------------------------------------------------------------
      function a360(a)
      a=mod(a,360.0)
      if (a.lt.0) a=a+360
      a360=a
      return
      end
c-----------------------------------------------------------------------
      function sinthl(ih,ik,il,ginv)
      dimension ginv(3,3),h(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
      integer data,t
      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-----------------------------------------------------------------------
      subroutine ucalc (io6,a,u1,u2,phix,phiy,phiz,u)
c
c Orientation matrix from lattice parameters and Denzo crystal rotation
c angles
c
c x = U*h
c
c     ( x )               ( h )
c x = ( y )    and    h = ( k )
c     ( z )               ( l )
c
c U = R(phiz)*R(phiy)*R(phix)*U0*A0
c
c      ( astar  bstar*cos(gammastar)  cstar*cos(betastar)            )
c A0 = ( 0      bstar*sin(gammastar) -cstar*sin(betastar)*cos(alpha) )
c      ( 0      0                     1/c                            )
c        
c      ( u1(1)  u1(2)  u1(3) )
c U0 = ( u2(1)  u2(2)  u2(3) )
c      ( u3(1)  u3(2)  u3(3) )
c
c where u1, u2, and u3 are unit vectors that define crystal-fixed
c Cartesian axes.
c
c In the Denzo convention:
c
c   u3 is parallel to the incident radiation beam.
c
c   u1 is parallel to the spindle axis, from the crystal to the
c   goniometer base, i.e., parallel to the detector frame x-axis.
c   u1 is designated by the reciprocal lattice axis nearest u1.
c
c   u2 is the so-called "vertical" axis parallel to the detector frame
c   y-axis; it is designated by the reciprocal lattice axis nearest the
c   cross product u3 X u1.
c
c For right-handed rotations about right-handed xyz-axes:
c                                             
c                                              +z
c                                         +z'  :        +y'
c                                          .   :       .
c           (  1      0          0     )    .  :     .
c R(phix) = (  0  cos(phix) -sin(phix) )     . :   .   phi-x
c           (  0  sin(phix)  cos(phix) )      .: .  \
c                                              :..........+y
c
c
c                                              +x
c                                         +x'  :        +z'
c                                          .   :       .
c           (  cos(phiy)  0  sin(phiy) )    .  :     .
c R(phiy) = (      0      1      0     )     . :   .   phi-y
c           ( -sin(phix)  0  cos(phix) )      .: .  \
c                                              :..........+z
c
c
c                                              +y
c                                         +y'  :        +x'
c                                          .   :       .
c           (  cos(phiz) -sin(phiz)  0 )    .  :     .
c R(phiz) = (  sin(phiz)  cos(phiz)  0 )     . :   .   phi-z
c           (      0          0      1 )      .: .  \
c                                              :..........+x
c
c
c The x,y,z-axes are orthonormal laboratory axes that are coincident
c with the u1,u2,u3-axes at phi(spindle) = 0.
c
c According to the goniometer head illustration in "The HKL Manual" by
c Gewirth, Otwinowski, and Minor (4th ed., Feb. 1995), the Denzo
c convention for positive senses for the x-, y-, and z-rotations is such
c that, if the x- and y-rotations are right-handed, the z-rotation is
c left-handed.
c
      dimension a(6),b(6),ginv(3,3),u1(3),u2(3),u3(3),v1(3),v2(3),v3(3),
     & a0(3,3),a1(3,3),u0(3,3),u(3,3),q(3,3),r(3,3)
      rad=0.017453293
      write (io6,'(/1x,
     &''Orientation matrix from lattice parameters and crystal rotati'',
     &''ons''/1x,
     &''============================================================='',
     &''==='')')
c
c Reciprocal space orthogonalization matrix
c
      call metric (a,b,ginv,vcell)
      a0(1,1)= b(1)
      a0(1,2)= b(2)*cos(b(6)*rad)
      a0(1,3)= b(3)*cos(b(5)*rad)
      a0(2,1)= 0
      a0(2,2)= b(2)*sin(b(6)*rad)
      a0(2,3)=-b(3)*sin(b(5)*rad)*cos(a(4)*rad)
      a0(3,1)= 0
      a0(3,2)= 0
      a0(3,3)= 1/a(3)
      write (io6,'(/1x,
     &''A0 Busing-Levy reciprocal space orthogonalization matrix''/1x,
     &''--------------------------------------------------------''/
     & (1x,3f12.6))') ((a0(i,j),j=1,3),i=1,3)
c
c Transposed inverse of reciprocal space orthogonalization matrix
c
      call minv3 (a0,a1,d)
      do i=1,3-1
      do j=i+1,3
        t=a1(i,j)
        a1(i,j)=a1(j,i)
        a1(j,i)=t
      end do
      end do
      write (io6,'(/1x,
     &''A0-inverse-transpose''/1x,
     &''--------------------''/
     & (1x,3f12.4))') ((a1(i,j),j=1,3),i=1,3)
c
c Orthogonal components of spindle axis [hkl]
c
      call mv (a0,u1,v1)
c
c Orthogonal components of "vertical" axis [hkl]
c
      call mv (a0,u2,v2)
c
c Unitary orientation matrix
c
      call unitv (v1)
      call unitv (v2)
      call vprod (v1,v2,v3)
      call unitv (v3)
      call vprod (v3,v1,v2)
      call unitv (v2)
      call vprod (v2,v3,v1)
      call unitv (v1)
      do j=1,3
        u0(1,j)=v1(j)
        u0(2,j)=v2(j)
        u0(3,j)=v3(j)
      end do
      write (io6,'(/1x,
     &''U0 unitary orientation matrix''/1x,
     &''-----------------------------''/
     & (1x,3f12.6))') ((u0(i,j),j=1,3),i=1,3)
c
c Misorientation corrections
c
      u(1,1)=1
      u(1,2)=0
      u(1,3)=0
      u(2,1)=0
      u(2,2)=1
      u(2,3)=0
      u(3,1)=0
      u(3,2)=0
      u(3,3)=1
      c=cos(phix*rad)
      s=sin(phix*rad)
c
c        ( 1  0  0 )
c r(x) = ( 0  c -s )    right-handed rotation about x.
c        ( 0 +s  c )
c
      r(1,1)= 1
      r(1,2)= 0
      r(1,3)= 0
      r(2,1)= 0
      r(2,2)= c
      r(2,3)=-s
      r(3,1)= 0
      r(3,2)=+s
      r(3,3)= c
      call mm (r,u,q)
      c=cos(phiy*rad)
      s=sin(phiy*rad)
c
c        ( c  0 +s )
c r(y) = ( 0  1  0 )    right-handed rotation about y.
c        (-s  0  c )
c
      r(1,1)= c
      r(1,2)= 0
      r(1,3)=+s
      r(2,1)= 0
      r(2,2)= 1
      r(2,3)= 0
      r(3,1)=-s
      r(3,2)= 0
      r(3,3)= c
      call mm (r,q,u)
      c=cos(phiz*rad)
      s=sin(phiz*rad)
c
c        ( c -s  0 )
c r(z) = (+s  c  0 )    right-handed rotation about z.
c        ( 0  0  1 )
c
      r(1,1)= c
      r(1,2)=-s
      r(1,3)= 0
      r(2,1)=+s
      r(2,2)= c
      r(2,3)= 0
      r(3,1)= 0
      r(3,2)= 0
      r(3,3)= 1
      call mm (r,u,q)
      write (io6,'(/1x,
     &''Rz*Ry*Rx misorientation rotation matrix''/1x,
     &''---------------------------------------''/
     & (1x,3f12.6))') ((q(i,j),j=1,3),i=1,3)
      call mm (u0,a0,r)
      call mm (q,r,u)
      write (io6,'(/1x,
     &''U = Rz*Ry*Rx*U0*A0 orientation matrix''/1x,
     &''-------------------------------------''/
     & (1x,3f12.6))') ((u(i,j),j=1,3),i=1,3)
      do i=1,3
      do j=1,3
        q(i,j)=u(j,i)
      end do
      end do
      call mm (q,u,r)
      write (io6,'(/1x,
     &''U-transpose*U''/1x,
     &''-------------''/ (1x,3f12.6))') ((r(i,j),j=1,3),i=1,3)
      d=u(1,1)*u(2,2)*u(3,3)+u(1,2)*u(2,3)*u(3,1)+u(1,3)*u(2,1)*u(3,2)
     &   -u(3,1)*u(2,2)*u(1,3)-u(3,2)*u(2,3)*u(1,1)-u(3,3)*u(2,1)*u(1,2)
      write (io6,'(/1x,''det(U) = '',e12.4)') d
      return
      end
c-----------------------------------------------------------------------
      subroutine unitv (x)
c
c Normalize to unit length a vector referred to Cartesian axes.
c
      dimension x(3)
      r=sqrt(x(1)**2+x(2)**2+x(3)**2)
      x(1)=x(1)/r
      x(2)=x(2)/r
      x(3)=x(3)/r
      return
      end
c-----------------------------------------------------------------------
      subroutine vm (a,b,c)
c
c Row vector-square matrix multiplication
c
c a(i)*b(i,j) = c(j)
c
      dimension a(3),b(3,3),c(3)
      do j=1,3
        c(j)=0
      do i=1,3
        c(j)=c(j)+a(i)*b(i,j)
      end do
      end do
      return
      end
c-----------------------------------------------------------------------
      subroutine vprod (x,y,z)
c
c Cross product of vectors referred to Cartesian axes
c
      dimension x(3),y(3),z(3)
      z(1)=x(2)*y(3)-x(3)*y(2)
      z(2)=x(3)*y(1)-x(1)*y(3)
      z(3)=x(1)*y(2)-x(2)*y(1)
      return
      end
c-----------------------------------------------------------------------