subroutine orient(mu,nu,kcu,gg)
c
c CONTACT: ap@marc.de
c      
c* * * * * *
c
c     user subroutine for input of orientation for anisotropic option.
c
c     mu           element number
c     nu           integration point number
c     kcu          layer number
c     gg           orientation matrix
c
c* * * * * *
c
c PURPOSE
c  enter preferred directions for 3D (anisotropic) problems:
c  alpha     = angle over which the cylindrical system is
c              rotated. the rotation is about the radial
c              direction.
c  when alpha is zero, the preferred system is cylindrical.
c  x1(ijk)   = coordinates of 1st point of cylinder axis
c  x2(ijk)   = coordinates of 2nd point of cylinder axis
c  xip(ijk)  = integration point coordinates
c  gg(1,ijk) = 1st preferred direction is radial
c  gg(2,ijk) = 2nd preferred direction is the circumferential
c              rotated over angle alpha about radial direction.
c  gg(3,ijk) = 3rd preferred direction is the axial
c              rotated over angle alpha about radial direction.
c
c USAGE
c  Define the cylinder axis through the "ORIENTATION" model
c  definition option by flagging UORIENT. The fields for the
c  three components of user vector 1 define the 1st point on
c  the cylinder axis. The fields for the three components of
c  user vector 2 define the 2nd point of the cylinder axis.
c  Define the angle alpha through the "ORIENTATION" option
c  (the field for the orientation angle).
c
c* * * * * *
c
      implicit real*8 (a-h,o-z)                                               dp
      dimension gg(3,3),m(2),xip(12),x1(3),x2(3),hh(3,3),e1(3),e2(3)
c
      parameter(pi=3.14159265358979d0)
c
c some common blocks referencing marc internal variables.
      include '../common/space'
      include '../common/heat'
      include '../common/lass'
      include '../common/dimen'
      include '../common/arrays'
      include '../common/array4'
c
c define a data base off set
      lofr = (n-1)*nelstr
c
c x1, x2 define 1st and 2nd point of desired cylinder axis
c
c get 1st point on the cylinder axis
      x1(1) = vars(iuvec1+lofr)
      x1(2) = vars(iuvec1+lofr+1)
      x1(3) = vars(iuvec1+lofr+2)
c
c get 2nd point on the cylinder axis
      x2(1) = vars(iuvec2+lofr)
      x2(2) = vars(iuvec2+lofr+1)
      x2(3) = vars(iuvec2+lofr+2)
c
c use el. integr. point coordinates to define the transformations
      la1 = icrxpt + (n-1)*nelstr + (nn-1)*ncrd
      do ijk=1,ncrd
       xip(ijk) = vars(la1+ijk-1)
      enddo
c
c 3rd pref. direction is along cylinder axis
      anorm = 0d0
      do ijk=1,ncrd
       gg(3,ijk) = x2(ijk) - x1(ijk)
       anorm = anorm + gg(3,ijk)**2
      enddo
      anorm = dsqrt(anorm)
      do ijk=1,ncrd
       gg(3,ijk) = gg(3,ijk)/anorm
      enddo
c
c 2nd pref. direction is circumferential
      gg(2,1) = gg(3,2)*(xip(3)-x1(3)) - gg(3,3)*(xip(2)-x1(2))
      gg(2,2) = gg(3,3)*(xip(1)-x1(1)) - gg(3,1)*(xip(3)-x1(3))
      gg(2,3) = gg(3,1)*(xip(2)-x1(2)) - gg(3,2)*(xip(1)-x1(1))
      anorm = 0d0
      do ijk=1,ncrd
       anorm = anorm + gg(2,ijk)**2
      enddo
      anorm = dsqrt(anorm)
      do ijk=1,ncrd
       gg(2,ijk) = gg(2,ijk)/anorm
      enddo
c
c 1st pref. direction is radial (outward)
      gg(1,1) = gg(2,2)*gg(3,3) - gg(2,3)*gg(3,2)
      gg(1,2) = gg(2,3)*gg(3,1) - gg(2,1)*gg(3,3)
      gg(1,3) = gg(2,1)*gg(3,2) - gg(2,2)*gg(3,1)
c
c get rotation angle from input
      alpha = vars(iangle+lofr)
      alpha = alpha/180d0*pi
      co    = dcos(alpha)
      si    = dsin(alpha)
c
c rotate about radial axis over angle alpha
      do i1=1,3
       e1(i1) = gg(2,i1)
       e2(i1) = gg(3,i1)
      enddo
      do i1=1,3
       gg(2,i1) =  co*e1(i1) + si*e2(i1)
       gg(3,i1) = -si*e1(i1) + co*e2(i1)
      enddo
c
      return
      end