      subroutine transform(ui,vi,bmat,u3,v3,b1mat)

      implicit none
c...	passed variables
      integer ui(3),vi(3)
      real bmat(3,3),b1mat(3,3)
      real u3(3),v3(3)
      
c...	local variables      
      integer i
      real tmat(3,3), b2mat(3,3)
      real u(3),v(3),u1(3),v1(3),u2(3),v2(3)
      real theta,pi
      real temp
      
      pi=4.0*atan(1.0)

C     The desired transformation is to a coordinate system where
C     x or zeta is perpendicular to uu, vv is a postive rotation
C     of A3 away from uu, and z is hence parallel to uu x vv.
C
C     Currently b[1] is parallel to x and b[2] lies in the x-y
C     plane because of the generating routine

C     Generate u and v in the current coordinate system.

      do i=1,3
        u(i)=float(ui(1))*bmat(i,1)+float(ui(2))*bmat(i,2)+
     &    float(ui(3))*bmat(i,3)
        v(i)=float(vi(1))*bmat(i,1)+float(vi(2))*bmat(i,2)+
     &    float(vi(3))*bmat(i,3)
      enddo
      
C     What is one of the transformation matrices that will align
C     u and x?  Will take 2 steps.

C      write(*,18) 'u is now ',u(1),' i+',u(2),' j+',u(3),' k'
C      write(*,18) 'v is now ',v(1),' i+',v(2),' j+',v(3),' k'
      
C     Rotation about z-axis with quadrant conventions.

      call calctheta(u(1),u(2),theta)

C      write(*,*) theta*180/pi
      
C     -theta appears in the transformation matix because we
C     want the transformation rotate the coordinate axies
C     so that u will lie in the x'-z' plane

      temp = -theta
      call rotate(temp,3,tmat)
c01Nov      call rotate(-theta,3,tmat)
      
      call mmmult(tmat,bmat,3,3,3,b1mat)
      call mvmult(tmat,u,3,3,u1)
      call mvmult(tmat,v,3,3,v1)
      
C      write(*,18) 'u is now ',u1(1),' i+',u1(2),' j+',u1(3),' k'
C      write(*,18) 'v is now ',v1(1),' i+',v1(2),' j+',v1(3),' k'


C     Rotation about y1-axis
C     This one is a bit different since with the y1-axis coming
C     out of the page x1-z1 are left-handed.  Best plan is to
C     choose theta such that it is correct for the +z1 axis
C     and then rotate pi/2 further to +x1.

      call calctheta(u1(3),u1(1),theta)

C      write(*,*) theta*180.0/pi
      temp = -theta+pi/2.0
      call rotate(temp, 2,tmat)
c01Nov      call rotate(-theta+pi/2.0,2,tmat)
      
      call mmmult(tmat,b1mat,3,3,3,b2mat)
      call mvmult(tmat,u1,3,3,u2)
      call mvmult(tmat,v1,3,3,v2)
      
C      write(*,18) 'u is now ',u2(1),' i+',u2(2),' j+',u2(3),' k'
C      write(*,18) 'v is now ',v2(1),' i+',v2(2),' j+',v2(3),' k'
 18   format(a9,f6.3,a3,f6.3,a3,f6.3,a2)

C     final rotation about x''-axis
      
      call calctheta(v2(2),v2(3),theta)

      temp = -theta
      call rotate(temp,1,tmat)
c01Nov      call rotate(-theta,1,tmat)

      call mmmult(tmat,b2mat,3,3,3,b1mat)
      call mvmult(tmat,u2,3,3,u3)
      call mvmult(tmat,v2,3,3,v3)

C      write(*,18) 'u is now ',u3(1),' i+',u3(2),' j+',u3(3),' k'
C      write(*,18) 'v is now ',v3(1),' i+',v3(2),' j+',v3(3),' k'

      return
      end

C*************************************************************

      subroutine rotate(theta,i,r)

C     Generates the passive transformation matrix for a 
C     counterclockwise rotation of the axes by theta radians 
C     about the i-th axis.

      implicit none
c...	passed variables
      integer i
      real theta, r(3,3)

c...	local variables
      integer j,k
      
C     j and k are the cyclic indicies from i
C     look for the code in generate.f
      
      call ijk(i,j,k)
      r(i,i)=1.0
      r(i,j)=0.0
      r(i,k)=0.0
      r(j,i)=0.0
      r(j,j)=cos(theta)
      r(j,k)=-sin(theta)
      r(k,i)=0.0
      r(k,j)=sin(theta)
      r(k,k)=cos(theta)

      return
      end

C****************************************************************

      subroutine calctheta(x,y,theta)

C     Find the angle theta from the +x-axis to vector
C     xi + yj over a range of -pi to pi (an extended
C     inverse trig function)

      implicit none
c...	passed variables
      real x,y,theta

      if ((x.eq.0.0).and.(y.eq.0.0)) then
        theta=0.0
        return
      endif

C     use arccosine to avoid any problems with x=0

C     Quadrants I and II are fine as is and theta values
C     will range from 0 to pi
      
      if (y.ge.0.0) then
        theta= acos(x/sqrt(x**2.0+y**2.0))
      else

C     Quadrants III and IV will need a negative sign since transform
C     would go to -x otherwise.

        theta=-acos(x/sqrt(x**2.0+y**2.0))

      endif

      return
      end
