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  preferred system is spherical for 3D problems
c  x1(ijk)   = coordinates of center point of sphere
c  x2(ijk)   = coordinates of north pole of sphere
c  xip(ijk)  = integration point coordinates
c  gg(1,ijk) = 1st preferred direction is radial
c  gg(2,ijk) = 2nd preferred direction is circumferential
c  gg(3,ijk) = 3rd preferred direction is normal to 1 and 2
c
c* * * * * *
c
      implicit real*8 (a-h,o-z)                                               dp
      dimension gg(3,3),m(2),xip(12),x1(3),x2(3)
c
      include '../common/space'
      include '../common/heat'
      include '../common/lass'
      include '../common/dimen'
      include '../common/arrays'
      include '../common/array4'
c
c define center point of sphere
      x1(1) = 0d0
      x1(2) = 0d0
      x1(3) = 0d0
c
c define the north pole of sphere
      x2(1) = 0d0
      x2(2) = 0d0
      x2(3) = 1d0
c
c define some epsilon value
      epsilon = 1d-6
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 1st pref. direction is radially from center to integr. point
      anorm = 0d0
      do ijk=1,ncrd
       gg(1,ijk) = xip(ijk) - x1(ijk)
       anorm = anorm + gg(1,ijk)**2
      enddo
      anorm = dsqrt(anorm)
      if (anorm.lt.epsilon) then
c if integr. point and center point coincide don't transform
       gg(1,1) = 1d0
       gg(1,2) = 0d0
       gg(1,3) = 0d0
       gg(2,1) = 0d0
       gg(2,2) = 1d0
       gg(2,3) = 0d0
       gg(3,1) = 0d0
       gg(3,2) = 0d0
       gg(3,3) = 1d0
       return
      endif
      do ijk=1,ncrd
       gg(1,ijk) = gg(1,ijk)/anorm
      enddo
c
c sphere axis is directed from center to north pole and stored in x1
      anorm = 0d0
      do ijk=1,ncrd
       x1(ijk) = x2(ijk) - x1(ijk)
       anorm = anorm + x1(ijk)**2
      enddo
      anorm = dsqrt(anorm)
      do ijk=1,ncrd
       x1(ijk) = x1(ijk)/anorm
      enddo
c
      anorm = 0d0
      do ijk=1,ncrd
       anorm = anorm + x1(ijk)*gg(1,ijk)
      enddo
      anorm = dabs(anorm)
      if (anorm.gt.1d0-epsilon) then
c
c if 1st dir. and polar axis dir. coincide then choose some 2nd dir.
       if (dabs(gg(1,1)).lt.epsilon) then
        gg(2,1) =  0d0
        gg(2,2) = -gg(1,3)
        gg(2,3) =  gg(1,2)
       else
        gg(2,1) =  gg(1,2)
        gg(2,2) = -gg(1,1)
        gg(2,3) =  0d0
       endif
      else
c
c 2nd pref. direction is circumferential around sphere axis
       gg(2,1) = x1(2)*gg(1,3) - x1(3)*gg(1,2)
       gg(2,2) = x1(3)*gg(1,1) - x1(1)*gg(1,3)
       gg(2,3) = x1(1)*gg(1,2) - x1(2)*gg(1,1)
      endif
      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)/anrom
      enddo
c
c 3rd pref. direction is right handed normal to the other two
      gg(3,1) = gg(1,2)*gg(2,3) - gg(1,3)*gg(2,2)
      gg(3,2) = gg(1,3)*gg(2,1) - gg(1,1)*gg(2,3)
      gg(3,3) = gg(1,1)*gg(2,2) - gg(1,2)*gg(2,1)
c
      return
      end