      logical function equal(Q1,Q2,SQ1,SQ2,TOL)
c
c Function to return equality of two quantities Q1 and Q2
c   with standard deviations SQ1 and SQ2 with tolerance TOL
c
      implicit none
      real*4 q1,q2,sq1,sq2,tol,diff,sdiff
c
c Calculate difference bewteen Q1 and Q2 and standard deviation of
c   the difference
c
      diff=q1-q2
      sdiff = sqrt(sq1**2+sq2**2)
      if( abs(diff) .le. tol*sdiff )then
c Q1 and Q2 are close to equal
          equal = .TRUE.
      else
c Q1 and Q2 are not equal
          equal = .FALSE.
      endif
      return
      end
      options/extend
      subroutine bravais
c
c Routine to find centered BRAVAIS lattice from reduced cell.
c   The method of Niggli (1928) as reproduced by Buerger (1957)
c   is employed.
c
      implicit none
      include 'single.inc'
      real*4 s(3,2),s11,s22,s33,s23,s31,s12,sigs(3,2),ss11,ss22,ss33
     1,ss23,ss31,ss12,sigcell,rad,tol,soff,abssoff,sigsum,trns(3,3)
     2,ubt(3,3),g(3,3),det,sprod
      integer*4 i,j
      common /lsqcom/sigcell(6)
      logical equal
      equivalence (s(1,1),s11),(s(2,1),s22),(s(3,1),s33),(s(1,2),s23)
     1,(s(2,2),s31),(s(3,2),s12),(sigs(1,1),ss11),(sigs(2,1),ss22)
     2,(sigs(3,1),ss33),(sigs(1,2),ss23),(sigs(2,2),ss31),(sigs(3,2),ss12)
      data rad/57.295/
c
c Do least-squares fitting of lattice constants and get sigmas
c
      call lsq(1)
c Initialize transformation matrix
      do i=1,3
          do j=1,3
              trns(j,i)=0.0
          enddo
      enddo
c
c Calculate G and Fill in S-matrix and Sigmas from G and SIGCELL
c
      call MatTsp(ub,ubt)
      call MatMul(ubt,ub,g)
      call minv(g,g,det)
      do i=1,3
          s(i,1) = g(i,i)
          sigs(i,1)=2.0*cell(i)*sigcell(i)
	  if(sigs(i,1).lt.1.0e-5)sigs(i,1)=1.0e-6*s(i,1)
      enddo
      s23=g(2,3)
      s31=g(3,1)
      s12=g(1,2)
      ss23=abs(s23*sqrt(sigcell(2)**2/s22+sigcell(3)**2/s33+sigcell(4)**2
     1/sin(cell(4)/rad)**2))
      if(ss23.lt.1.0e-5)ss23=sqrt(ss22*ss33)
      ss31=abs(s31*sqrt(sigcell(3)**2/s33+sigcell(1)**2/s11+sigcell(5)**2
     1/sin(cell(5)/rad)**2))
      if(ss31.lt.1.0e-5)ss31=sqrt(ss33*ss11)
      ss12=abs(s12*sqrt(sigcell(1)**2/s11+sigcell(2)**2/s22+sigcell(6)**2
     1/sin(cell(6)/rad)**2))
      if(ss12.lt.1.0e-5)ss12=sqrt(ss11*ss22)
      type 10
10    format(/'$Bravais Lattice Search - Type Tolerance in ESD (DEF=3)?')
      read(5,'(f10.0)',end=9999)tol
      if(tol.le.0.0)tol=3.0
c Get sum of off-diagonal elements, sum of absolute values, and sigma
      soff=s23+s31+s12
      abssoff=abs(s12)+abs(s31)+abs(s23)
      sigsum=sqrt(ss23**2+ss31**2+ss12**2)
      if(equal(abssoff,0.0,sigsum,0.0,tol))then
c All off-diagonal elements are zero - return immediately
          type 20
20        format(/' Initial Cell is Orthogonal - No Transformation')
          return
      endif
      if(equal(s11,s22,ss11,ss22,tol))then
c s11 = s22
          if(equal(s22,s33,ss22,ss33,tol))then
c s11 = s22 = s33
              if(soff.gt.tol*sigsum)then
c unsymmetrical scalers are positive
                  if(.not.equal(s23,s31,ss23,ss31,tol))then
c symmetrical scalers are accidentally equal - check condition IIIc
                      go to 3000
                  else
                      if(.not.equal(s23,s12,ss23,ss12,tol))then
c symmetrical scalers are accidentally equal - This is condition IIb
                          go to 2000
                      else
                          if(equal(s12,0.5*s11,ss12,0.5*ss11,tol))then
c F-centered cubic - Transform with matrix 44C
                              type *,'F-centered Cubic (44C)'
                              trns(1,1)=1.0
                              trns(1,2)=-1.0
                              trns(1,3)=1.0
                              trns(2,1)=1.0
                              trns(2,2)=1.0
                              trns(2,3)=-1.0
                              trns(3,1)=-1.0
                              trns(3,2)=1.0
                              trns(3,3)=1.0
                              call trans(trns,1)
                              return
                          else
c R-centered - transform with 49C
                              type *,'R-centered (49C)'
                              trns(1,1)=1.0
                              trns(1,2)=-1.0
                              trns(2,1)=-1.0
                              trns(2,3)=1.0
                              trns(3,1)=-1.0
                              trns(3,2)=-1.0
                              trns(3,3)=-1.0
                              call trans(trns,1)
                              return
                          endif
                      endif
                  endif
              else
c Unsymmetrical scalers are negative
                  if(.not.equal(s11,abs(soff),ss11,sigsum,tol))then
                      if(equal(abs(s23),abs(s31),ss23,ss31,tol).and.
     1                       equal(abs(s23),abs(s12),ss23,ss12,tol))then
c R-centered - Transform with matrix 49D
                          type *,'R-centered (49D)'
                          trns(1,1)=1.0
                          trns(1,2)=-1.0
                          trns(2,1)=-1.0
                          trns(2,3)=1.0
                          trns(3,1)=-1.0
                          trns(3,2)=-1.0
                          trns(3,3)=-1.0
                          call trans(trns,1)
                          return
                      else
                          go to 2000
                      endif
                  else
c Sum of Unsymmetrical scaler = s11 - Count Number that are equal
		      i=0	
                      if(equal(abs(s23),abs(s31),ss23,ss31,tol))i=i+1
                      if(equal(abs(s23),abs(s12),ss23,ss12,tol))i=i+1
                      if(equal(abs(s12),abs(s31),ss12,ss31,tol))i=i+1
                      if(i.eq.3)then
c I-centered cubic - transform with matrix 44B
                          type *,'I-centered Cubic (44B)'
                          trns(1,1)=1.0
                          trns(1,3)=1.0
                          trns(2,1)=1.0
                          trns(2,2)=1.0
                          trns(3,2)=1.0
                          trns(3,3)=1.0
                          call trans(trns,1)
                          return
                      endif
		      if(i.eq.1)then
c I-centered tetragonal - transform with matrix 45D
                          type *,'I-centered Tetragonal (45D-P)'
                          trns(1,2)=1.0
                          trns(1,3)=1.0
                          trns(2,1)=1.0
                          trns(2,2)=1.0
                          trns(3,1)=-1.0
                          trns(3,3)=-1.0
                          call trans(trns,1)
                          return
                      else
c I-centered orthorhombic - Transform with 52A
                          type *,'I-centered Orthorhombic (52A)'
                          trns(1,1)=1.0
                          trns(1,3)=1.0
                          trns(2,1)=1.0
                          trns(2,2)=1.0
                          trns(3,2)=1.0
                          trns(3,3)=1.0
                          call trans(trns,1)
                          return
                      endif
                  endif
              endif
          else
c  S11 = S22 .ne. S33
2000          if(soff.gt.tol*sigsum)then
c unsymmetrical scalers are positive
                  if(equal(s23,s31,ss23,ss31,tol))then
                      if(equal(s23,s12,ss23,ss12,tol)
     1			.and.equal(s12,0.5*s11,ss12,0.5*ss11,tol))then
                          type *,'R-centered (49B)'
                          trns(1,1)=1.0
                          trns(2,1)=-1.0
                          trns(2,2)=1.0
                          trns(3,1)=-1.0
                          trns(3,2)=-1.0
                          trns(3,3)=3.0
                          call trans(trns,1)
                          return
                      else
                          type *,'C-centered Monoclinic (55A-P)'
                          trns(1,1)=-1.0
                          trns(1,2)=-1.0
                          trns(2,1)=1.0
                          trns(2,2)=-1.0
                          trns(3,3)=1.0
                          call trans(trns,1)
                          return
                      endif
                  else
                      go to 4000
                  endif
              else
c Unsymmetrical scalers are negative
                  if(.not.equal(abs(soff),s11,sigsum,ss11,tol))then
                      if(.not.equal(abs(s23),abs(s31),ss23,ss31,tol))then
                          go to 4000
                      else
                          if(.not.equal(abs(s23),0.0,ss23,0.0,tol))then
                              type *,'C-centered Monoclinic (55A)'
                              trns(1,1)=1.0
                              trns(1,2)=1.0
                              trns(2,1)=-1.0
                              trns(2,2)=1.0
                              trns(3,3)=1.0
                              call trans(trns,1)
                              return
                          else
                              if(equal(abs(s12),0.5*s11,ss12,0.5*ss11,tol))then
				  type *,'P-centered Hexagonal (48A)'
                                  return
                              else
                                  type *,'C-centered Orthorhombic (50D)'
                                  trns(1,1)=-1.0
                                  trns(1,2)=1.0
                                  trns(2,3)=1.0
                                  trns(3,1)=1.0
                                  trns(3,2)=1.0
                                  call trans(trns,1)
                                  return
                              endif
                          endif
                      endif
                  else
                      if(equal(abs(s12),0.0,ss12,0.0,tol))then
                          type *,'I-centered Tetragonal (45C)'
                          trns(1,1)=1.0
                          trns(2,2)=1.0
                          trns(3,1)=1.0
                          trns(3,2)=1.0
                          trns(3,3)=2.0
                          call trans(trns,1)
                          return
                      else
                          if(equal(abs(s23),abs(s31),ss23,ss31,tol))then
                              type *,'F-centered Orthorhombic (51A)'
                              trns(1,1)=1.0
                              trns(1,2)=-1.0
                              trns(2,1)=1.0
                              trns(2,2)=1.0
                              trns(2,3)=2.0
                              trns(3,1)=-1.0
                              trns(3,2)=-1.0
                              call trans(trns,1)
                              return
                          else
                              type *,'I-centered Monoclinic (57B)'
                              trns(1,2)=-1.0
                              trns(1,3)=-1.0
                              trns(2,1)=1.0
                              trns(2,2)=1.0
                              trns(3,1)=1.0
                              trns(3,3)=1.0
                              call trans(trns,1)
                              return
                          endif
                      endif
                  endif
              endif
          endif
      else
c S11 .ne. S22
3000      if(equal(s22,s33,ss22,ss33,tol))then
c ******************
c S11 .ne. S22 = S33
c ******************
              if(soff.gt.tol*sigsum)then
c Unsymmetric scalers are positive
                  if(.not.equal(s31,s12,ss31,ss12,tol))go to 4000
                  if(equal(s31,0.5*s11,ss31,0.5*ss11,tol))then
                      if(equal(s23,0.25*s11,ss23,0.25*ss11,tol))then
                          type *,'I-centered Tetragonal (45E)'
                          trns(1,2)=-1.0
                          trns(1,3)=1.0
                          trns(2,1)=1.0
                          trns(2,2)=-1.0
                          trns(2,3)=-1.0
                          trns(3,1)=1.0
                          call trans(trns,1)
                          return
                      else
                          type *,'I-centered Orthorhombic (52B)'
                          trns(1,2)=-1.0
                          trns(1,3)=1.0
                          trns(2,1)=-1.0
                          trns(2,2)=1.0
                          trns(2,3)=1.0
                          trns(3,1)=-1.0
                          call trans(trns,1)
                          return
                      endif
                  else
                      type *,'C-centered Monoclinic (55B-P)'
                      trns(1,2)=-1.0
                      trns(1,3)=-1.0
                      trns(2,2)=1.0
                      trns(2,3)=-1.0
                      trns(3,1)=1.0
                      call trans(trns,1)
                      return
                  endif
              else
c Unsymmetrical scalers are negative
                  if(.not.equal(abs(s31),abs(s12),ss31,ss12,tol))go to 4000
                  if(equal(abs(s12),0.0,ss12,0.0,tol))then
                      if(equal(abs(s23),0.5*s22,ss23,0.5*ss22,tol))then
                          type *,'P-centered Hexagonal (48B)'
                          trns(1,2)=1.0
                          trns(2,3)=1.0
                          trns(3,1)=1.0
                          call trans(trns,1)
                          return
                      else
                          type *,'C-centered Orthorhombic (50E)'
                          trns(1,2)=1.0
                          trns(1,3)=1.0
                          trns(2,2)=-1.0
                          trns(2,3)=1.0
                          trns(3,1)=1.0
			  call trans(trns,1)
			  return
                      endif
                  else
                      if(.not.equal(abs(s12),s11/3.0,ss12,ss11/3.0,tol))then
                          type *,'C-centered Monoclinic (55B)'
                          trns(1,2)=1.0
                          trns(1,3)=1.0
                          trns(2,2)=-1.0
                          trns(2,3)=1.0
                          trns(3,1)=1.0
                          call trans(trns,1)
                          return
                      else
			  sprod=sqrt(ss11*ss11/36.+0.25*ss22*ss22)
                          if(equal(abs(s23),0.5*s22-s11/6.,ss23,sprod,tol))then
                              type *,'R-centered (49E)'
                              trns(1,1)=1.0
                              trns(1,2)=2.0
                              trns(1,3)=1.0
                              trns(2,2)=-1.0
                              trns(2,3)=1.0
                              trns(3,1)=1.0
                              call trans(trns,1)
                              return
                          else
                              go to 4000
                          endif
                      endif
                  endif
              endif
          else
c S11 .ne. S22 .ne. S33
4000          if(soff.gt.tol*sigsum)then
c Unsymmetrical scalers are positive
                  if(equal(s12,0.5*s11,ss12,0.5*ss11,tol))then
                      if(equal(s31,0.5*s11,ss31,0.5*ss11,tol))then
                          if(equal(s23,0.25*s11,ss23,0.25*ss11,tol))then
                              type *,'F-centered Orthorhombic (51B)'
                              trns(1,1)=-1.0
                              trns(1,2)=2.0
                              trns(2,1)=-1.0
                              trns(2,3)=2.0
                              trns(3,1)=1.0
                              call trans(trns,1)
                              return
                          else
                              type *,'I-centered Monoclinic (57C)'
			      trns(1,2)=1.0
			      trns(1,3)=-1.0
			      trns(2,1)=-1.0
			      trns(3,1)=-1.0
			      trns(3,2)=1.0
			      trns(3,3)=1.0
			      call trans(trns,1)
                              return
                          endif
                      else
                          if(equal(s23,0.5*s31,ss23,0.5*ss31,tol))then
                              type *,'C-centered Monoclinic (56C)'
                              trns(1,1)=-1.0
                              trns(2,1)=1.0
                              trns(2,2)=-2.0
                              trns(3,3)=1.0
                              call trans(trns,1)
                              return
                          else
                              go to 5000
                          endif
                      endif
                  else
                      if(equal(s31,0.5*s11,ss31,0.5*ss11,tol))then
                          if(equal(s23,0.5*s12,ss23,0.5*ss12,tol))then
                              type *,'C-centered Monoclinic (56C)'
                              trns(1,1)=1.0
                              trns(2,1)=-1.0
                              trns(2,2)=2.0
                              trns(3,3)=1.0
                              call trans(trns,1)
                              return
                          else
                              go to 5000
                          endif
                      else
                          if(equal(s31,0.5*s12,ss31,0.5*ss12,tol).and.
     1                          equal(s23,0.5*s22,ss31,0.5*ss22,tol))then
                              type *,'C-centered Monoclinic (56B)'
                              trns(1,2)=-1.0
                              trns(2,2)=1.0
                              trns(2,3)=-2.0
                              trns(3,1)=1.0
                              call trans(trns,1)
                              return
                          else
                              go to 5000
                          endif
                      endif
                  endif
              else
c Unsymmetrical Scalers < 0
		  sprod=0.5*sqrt(ss11**2+ss22**2)
                  if(equal(abs(soff),0.5*(s11+s22),sigsum
     1                  ,sprod,tol))then
                      if(equal(abs(s12),0.0,ss12,0.0,tol))then
                          type *,'I-centered Orthorhombic (52C)'
                          trns(1,1)=-1.0
                          trns(2,2)=1.0
                          trns(3,1)=-1.0
                          trns(3,2)=-1.0
                          trns(3,3)=-2.0
                          call trans(trns,1)
                          return
                      else
                          type *,'I-centered Monoclinic (57A)'
                          trns(1,1)=-1.0
                          trns(2,1)=1.0
                          trns(2,2)=1.0
                          trns(2,3)=2.0
                          trns(3,2)=1.0
                          call trans(trns,1)
                          return
                      endif
                  endif
                  if(equal(abs(s23),0.5*s22,ss23,ss22*0.5,tol))then
                      if(equal(s12,0.0,ss12,0.0,tol))then
                          if(equal(abs(s31),0.0,ss31,0.0,tol))then
                              type *,'C-centered Orthorhombic (50F)'
                              trns(1,2)=1.0
                              trns(2,2)=-1.0
                              trns(2,3)=-2.0
                              trns(3,1)=-1.0
                              call trans(trns,1)
                              return
                          else
                              type *,'C-centered Monoclinic (54B)'
                              trns(1,2)=1.0
                              trns(1,3)=2.0
                              trns(2,2)=1.0
                              trns(3,1)=-1.0
                              call trans(trns,1)
                              return
                          endif
                      else
                          go to 5000
                      endif
                  endif
                  if(equal(abs(s12),0.5*s11,ss12,0.5*ss11,tol))then
                      if(equal(abs(s31),0.0,ss31,0.0,tol))then
                          if(equal(abs(s23),0.0,ss23,0.0,tol))then
                              type *,'C-centered Orthorhombic (50B)'
                              trns(1,1)=-1.0
                              trns(2,1)=-1.0
                              trns(2,2)=-2.0
                              trns(3,3)=1.0
                              call trans(trns,1)
                              return
                          else
                              type *,'C-centered Monoclinic (54A)'
                              trns(1,1)=1.0
                              trns(1,2)=2.0
                              trns(2,1)=-1.0
                              trns(3,3)=1.0
                              call trans(trns,1)
                              return
                          endif
                      else
                          go to 5000
                      endif
                  endif
                  if(equal(abs(s31),0.5*s11,ss31,0.5*ss11,tol))then
                      if(equal(abs(s12),0.0,ss12,0.0,tol))then
                          if(equal(abs(s23),0.0,ss23,0.0,tol))then
                              type *,'C-centered Orthorhombic (50A)'
                              trns(1,1)=1.0
                              trns(2,1)=-1.0
                              trns(2,3)=-2.0
                              trns(3,2)=1.0
                              call trans(trns,1)
                              return
                          else
                              type *,'C-centered Monoclinic (54C)'
                              trns(1,1)=1.0
                              trns(1,3)=2.0
                              trns(2,1)=1.0
                              trns(3,2)=1.0
                              call trans(trns,1)
                              return
                          endif
                      else
                          go to 5000
                      endif
                  endif
                  if(equal(abs(s23),0.0,ss23,0.0,tol))then
                      if(equal(abs(s12),0.0,ss12,0.0,tol))then
                          go to 5000
                      endif
                      if(equal(abs(s31),0.0,ss31,0.0,tol))then
                          type *,'P-centered Monoclinic (53C)'
                          trns(1,1)=1.0
                          trns(2,3)=-1.0
                          trns(3,2)=1.0
                          call trans(trns,1)
                          return
                      endif
                  else
                      if(equal(abs(s31),0.0,ss31,0.0,tol)
     1                      .and.equal(abs(s12),0.0,ss12,0.0,tol))then
                          type *,'P-centered Monoclinic (53B)'
                          trns(1,2)=1.0
                          trns(2,1)=-1.0
                          trns(3,3)=1.0
                          call trans(trns,1)
                          return
                      endif
                  endif
              endif
          endif
      endif
5000  type *,'No Transformation Performed'
9999  return
      end
