      subroutine calclattice(u,v,bmat, m_rlv, h,k,l,
     &			zeta,eta,hklcen,uint,vint)
C      CALL rlattice(uuv,vuv,bmatuv,m,hh,kk,ll,zeta,eta)
      
C     Calculate reciprocal lattice after u, v, and bmat have been
C     transformed to u-v space. (i.e. u || x, u x v || z) in terms
C     of h(i,j), k(i,j), l(i,j), zeta(i,j), and eta(i,j).

C     Also want the reciprocal lattice centred about xcen and ycen.
C     Can do this by changing the offset n.

      implicit none
      include 'arraysize.inc'
c...	passed variables
      real u(3),v(3),bmat(3,3)
      integer m_rlv
      integer uint(3),vint(3)
      real zeta(m_rlv_max,m_rlv_max)
      real eta(m_rlv_max,m_rlv_max)
      integer h(m_rlv_max,m_rlv_max)
      integer k(m_rlv_max,m_rlv_max)
      integer l(m_rlv_max,m_rlv_max)
      integer hklcen(3)

c...	local variables
      integer n,i,j
      integer uirr(3),virr(3)
      character*10 inttitle,redtitle
      integer len1,len2

c	November 2001
      real	temp1, temp2, temp3
      
C     The reciprocal lattice in this scattering plane is
C     all combinations of reciprocal lattice vectors such
C     that (h*A + k*B + l*C).(u x v)=0
C     or i*bmatuv(*,1).uxv(*)+j*bmatuv(*,2).uxv(*)+k*bmatuv(*,3).uxv(*)=0
C     A similar result should come from linear,integer, addition of u
C     and v, and contain the reciprocal lattice vectors which both
C     satisfy u.(u x v)=0 and v.(u x v)=0.       
C
C     if u and v are not reduced to the lowest common denominator then
C     those points will be simply skipped      
C     is there any good way to reduce to the shortest vectors that
C     still define this plane and won't miss any points, like an
C     irreducible set ?
      
      call irreduc(uint,vint,uirr,virr)
      if ((uint(1).ne.uirr(1)).or.(uint(2).ne.uirr(2)).or.
     &  (uint(3).ne.uirr(3))) then
        call title(uint,inttitle,len1)
        call title(uirr,redtitle,len2)
        write(*,'(3a)') inttitle(1:len1),' reduces to ',redtitle(1:len2)
      endif
      if ((vint(1).ne.virr(1)).or.(vint(2).ne.virr(2)).or.
     &  (vint(3).ne.virr(3))) then
        call title(vint,inttitle,len1)
        call title(virr,redtitle,len2)
        write(*,'(3a)') inttitle(1:len1),' reduces to ',redtitle(1:len2)
      endif
      
      n=int(m_rlv/2)+1
        do i=1,m_rlv
          do j=1,m_rlv
            h(i,j)=hklcen(1)+(i-n)*uirr(1)+(j-n)*virr(1)
            k(i,j)=hklcen(2)+(i-n)*uirr(2)+(j-n)*virr(2)
            l(i,j)=hklcen(3)+(i-n)*uirr(3)+(j-n)*virr(3)
	    temp1 = float(h(i,j))
	    temp2 = float(k(i,j))
	    temp3 = float(l(i,j))
            call qspace(bmat, temp1, temp2, temp3, 
     &        zeta(i,j),eta(i,j))
c01Nov            call qspace(bmat,float(h(i,j)),float(k(i,j)),
c01Nov     &        float(l(i,j)),zeta(i,j),eta(i,j))
cDEbug	    write (*,*) 'calclattice:', i, j, h(i,j), k(i,j), l(i,j)
          enddo
        enddo


      return
      end
C***************************************************************

      subroutine qspace(bmat,h,k,l,x,y)

C     Gives x,y,z components of reciprocal lattice points for
C     Miller indicies h,k,l and reciprocal lattice vectors
C     bmat(*,1), bmat(*,2), and bmat(*,3).
      
      implicit none
c...	passed variables
      real bmat(3,3),x,y
      real h,k,l
c...	local variables
      real z
      
      x=h*bmat(1,1)+k*bmat(1,2)+l*bmat(1,3)
      y=h*bmat(2,1)+k*bmat(2,2)+l*bmat(2,3)
      z=h*bmat(3,1)+k*bmat(3,2)+l*bmat(3,3)

C     If everything is correct then h,k,l and bmat should be such
C     that there is no z-component to the scattering.

      if (abs(z).gt.1.1e-5) then
        write (*,*)
     &  'Reciprocal lattice point out-of scattering plane'
        write(*,*) x,y,z,h,k,l
      endif
      
      return
      end
C
C*************************************************************      
C
      subroutine irreduc(u,v,uirr,virr)

      implicit none

c...	passed variables
      integer u(3),v(3),uirr(3),virr(3)

c...	local variables
      integer diff(3)
      integer ucoeff,vcoeff
      integer j, k
      integer n, diffsq, vsq
      
C  generate irreducible pair of vectors uirr,virr as close
C  to right angles as possible that are still in the
C  same plane as u and v

C  first check for any common denominators using division by primes

      do j=1,3
        uirr(j)=u(j)
        virr(j)=v(j)
      enddo
      
C      write(*,*) (uirr(i),i=1,3)
      call lcd(uirr,n)
C      write(*,*) (virr(i),i=1,3)      
      call lcd(virr,n)

C what is the possibility that an integer sum of the peaks will yield
C zero for one of the indices?

 11   continue
      do j=1,3
        call mult(uirr(j),virr(j),ucoeff,vcoeff)
C        write(*,*) 'Coefficients are ',ucoeff,vcoeff
        if ((ucoeff.eq.0).or.(vcoeff.eq.0)) goto 25
        do k=1,3
          diff(k)=ucoeff*uirr(k)+vcoeff*virr(k)
        enddo

C        write(*,*) 'Linear combination of u and v ',(diff(k),k=1,3)
        call lcd(diff,n)
C        if (n.ne.1) then
        diffsq=0
        vsq=0
        do k=1,3
          diffsq=diffsq+diff(k)**2
          vsq=vsq+virr(k)**2
        enddo
        if (diffsq.lt.vsq) then
          do k=1,3
            virr(k)=diff(k)
          enddo
          goto 11
        endif
  25  enddo
      
      return
      end
C
C*****************************************************
C      
      subroutine lcd(uirr,fac)
      implicit none

c...	constants
      integer nprime
      parameter (nprime=26)

c...	passed variables
      integer fac,uirr(3)

c...	local variables
      integer i,j
      integer prime(nprime)
      logical comm

      data (prime(j),j=1,nprime)/2,3,5,7,11,13,17,19,23,29,31,37,
     &  41,43,47,53,59,61,67,71,73,79,83,89,97,101/
      
      fac=1

      if ((uirr(1).eq.0).and.(uirr(2).eq.0).and.(uirr(3).eq.0))
     &  goto 16


  10  continue
      do i=1,nprime
        comm=.true.
        do j=1,3
           comm=(prime(i)*(uirr(j)/prime(i)).eq.uirr(j)).and.comm
        enddo
        if (comm) then
C          write(*,*) (uirr(j),j=1,3),' is being divided by ',prime(i)
          do j=1,3
            uirr(j)=uirr(j)/prime(i)
          enddo
          fac=fac*prime(i)
          go to 10
        endif
      enddo

C 16   if (fac.eq.1) write(*,*) 'No common factors for ',(uirr(j),j=1,3)
C include a quick check of negatives
 16   comm=.true.
      do j=1,3
        comm=(-uirr(j).eq.abs(uirr(j))).and.comm
      enddo
      if (comm) then
        do j=1,3
          uirr(j)=-uirr(j)
        enddo
      endif
      
      return
      end

C*******************************************************************
      
      subroutine mult(x,y,a,b)
      implicit none

c...	passed variables
      integer x,y,a,b
c...	local variables
      integer n(3),fac

C find the solution to ax+by=0
C with no common factor between a and b

      if ((x.eq.0).and.(y.eq.0)) then
        a=0                   ! both a and b indeterminant so this
        b=0                   ! really a flag not to do any reduction
        goto 66
      else if (x.eq.0) then
        b=0                   ! a indeterminant
        a=1
        goto 66
      else if (y.eq.0) then
        a=0                   ! b indeterminant
        b=1
        goto 66
      endif

C so a and b are determined up to some multiple and the solution
C a=y, b=-x is fine, dump it into lcd to check for any common
C factors     
      n(1)=y
      n(2)=-x
      n(3)=0
      call lcd(n,fac)
      a=n(1)
      b=n(2)

 66   return
      end
