      subroutine calcki(bmat,h,ei,ef,q,ki,fail,phi,psi,pn)
      implicit none
      include 'spur_warn.inc'
      include 'arraysize.inc'
      include 'spur_scan.inc'
            
c...	common blocks
      real branchdeg
      common/specparm/ branchdeg
      
c...	functions
      real dot
      
c...	passed variables
      real bmat(3,3),h(3),q(3),ei,ef
      real phi,psi
      real ki(2)
      logical fail
      integer pn

c...	local variables      
      real magki,magkf,magq
      real alpha,beta,pi
      real kf(2)
      real cosphi
      logical cross
      real branch
      integer i, j
      
      
      
      pi=4.0*atan(1.0)
      fail=.false.
      cross=.false.
C     How are hc, kc, lc resolved in terms of zeta and eta?
C     Q= h*A + k*B + l*C, in terms of (hkl) and recip vec A,B,C
C     Q= h(1)*bmatuv(*,1)+h(2)*bmatuv(*,2)+h(3)+bmatuv(*,3)

      call mvmult(bmat,h,3,3,q)
      magq=sqrt(dot(q,q))

cDebug      write (*,*) 'calcki: ei, ef =', ei, ef
cDebug      write (*,*) 'calcki: q=', (q(i), i=1,3)
cDebug      write (*,*) 'calcki: magq=', magq

c...	Looks like some memory is getting stomped on.
c...	If change a4 to 115., then get q(3) ne 0 errors

cDebug      if (pn .eq. 1) then
cDebug	do i = 1, 3
cDebug	  write (*,*) 'calcki: bmat(', i, ', j)=',  (bmat(i,j), j=1,3)
cDebug	enddo
cDebug	write (*,*) 'calcki: H,K,L=',  (h(i), i=1,3)
cDebug        write (*,*) 'calcki: q=', (q(i), i=1,3), ' magq = ', magq
cDebug      endif

      if (abs(q(3)).gt.1e-5) then
         write(*,'(a,i2)')'calcki: Q out of scattering plane! Point#',pn
c...		 			 123456789_123456789_
	 call set_warning(wbit_qout, 1, 'Q out of plane!     ')	 
         fail=.true.
         if (pn .eq. 1) then
	   do i = 1, 3
	    write (*,*) 'calcki: bmat(', i, ', j)=',  (bmat(i,j), j=1,3)
	   enddo
	   write (*,*) 'calcki: H,K,L=',  (h(i), i=1,3)
           write (*,*) 'calcki: q=', (q(i), i=1,3), ' magq = ', magq
         endif
         return
      endif

C      write(*,'(a,f5.1,1x,f5.1)') 'Incident and final energies: ',ei,ef

C     Calculate the lengths of the scattering vectors.
C     ei and ef are already positive      

      magki=sqrt(ei/k2_to_E)
      magkf=sqrt(ef/k2_to_E)
c      magki=sqrt(ei/2.07213)
c      magkf=sqrt(ef/2.07213)

C      write(*,'(a,f5.2,1x,f5.2,1x,f5.2)')
C     &  'Incident, final, and Q wavevectors: ',magki,magkf,magq


C     Scattering angle phi can be calculated simply from the
C     cosine law applied to the three scattering vectors.
      
C  check to make sure the triangle can be closed

      cosphi=(magki**2+magkf**2-magq**2)/(2.0*magki*magkf)
      if(abs(cosphi).gt.1.0) then
         write(*,'(2a,i2)') 'Unable to close scattering triangle! ',
     &    'Point #',pn
c...		 			 123456789_123456789_
	 call set_warning(wbit_qout, 1, 'Can''t close triangl ')	 
         fail=.true.
         a4_ok(pn) = .false.
         return
      else
        a4_ok(pn) = .true.
      endif

      phi=acos(cosphi)
C      write(*,*) 'cos(phi), phi are ',cosphi,phi*180/pi,cos(phi)

C     Once phi is known then the sine law can be used to
C     calculate the angle between from Q to ki which I have
C     called alpha.
C     Alpha will be positive in the sense of the positive
C     angle of the spectrometer but in the u-v RH coordinate
C     system it will be a negative (CCW) angle.
      
C changed this because range of asin does not allow for obtuse angles
C      alpha=asin(magkf*sin(phi)/magq)
      alpha=acos((magki**2+magq**2-magkf**2)/(2.0*magki*magq))
C this will give the other orientation
C      alpha=-acos((magki**2+magq**2-magkf**2)/(2.0*magki*magq))
C      write(*,'(a,f7.3)') 'Alpha-angle is ',alpha*180/pi

C     The next unknown is the angle between Q and u which
C     is obtained from the expression cos beta = a.b/|a||b|.
C     Since we now have the reciprocal lattice vectors in
C     a space where u // x -> cos beta = Qzeta/magQ.
C     This will return values between 0 and pi (quadrants
C     1 and 2) but will either need a negative sign or a
C     shift by pi in quadrants 3 and 4.

      call calctheta(q(1),q(2),beta)

C      write(*,'(a,f8.3)') 'Beta-angle is ',beta*180/pi
      
C     to calculate ki in the RH u-v coordinate system it
C     is a vector with length ki and angle beta-alpha wrt u

C      write(*,'(a,f8.3)') 'Beta-alpha angle is ',(beta-alpha)*180/pi

      ki(1)=magki*cos(beta-alpha)
      ki(2)=magki*sin(beta-alpha)
      kf(1)=ki(1)-q(1)
      kf(2)=ki(2)-q(2)

      if (abs(kf(1)**2.0+kf(2)**2.0-magkf**2).gt.4.0e-4) then
         write(*,'(a)')'Problem with kf calculation.'
         write(*,'(a,i2)') 'Result is likely incorrect at point# ',pn
         write(*,*) 'calc, actual ',kf(1)**2.0+kf(2)**2.0,magkf**2
         write(*,*) 'diff ',kf(1)**2.0+kf(2)**2.0-magkf**2
c...		 			 123456789_123456789_
	 call set_warning(wbit_qout, 1, 'Problem with kf calc')	 
         fail=.true.
C         return
      endif
      
C     A3 on the other hand is really pi/2 plus this angle
C     modulo "branch cut" conventions to prevent A3 wrapping.
C     The current range of psi is perhaps ill-defined since
C     it depends on the quadrant of Q and then switches dis-
C     continuously by -2 pi. 

      psi=beta-alpha+pi/2

C     branch set to 360 will have no effect in this config

      branch=pi/180.0*(branchdeg)
c01Dec      branch=pi/180*(branchdeg)
      if ((psi.le.branch).and.(.not.cross)) then
         write(*,'(a,f6.2,a,i2)') 'A3 has crossed a branch cut at',
     &        branchdeg,' at point #',pn
         cross=.true.
         psi=2*pi-psi
      endif

C      write(*,'(a,f7.3,a,f8.3,a1)') 'A4 (phi) is ',phi*180/pi,
C     &  ' and A3 (psi) is ',psi*180/pi,'.'
      
      return
      end


