c$Id: cubiclin.f,v 1.3 2002/09/09 23:32:01 sanjay Exp $
        subroutine cubiclin(isw,eps,theta,ud,sig,kk)
c      * * F E A P * * A Finite Element Analysis Program

c....  Copyright (c) 1999: Sanjay Govindjee

c-----[--.----+----.----+----.-----------------------------------------]
c     Purpose: User Constitutive Model Routine for a single
c              crystal cubic material small strain formulation

c     Input:
c          isw = task switch (3,6,4,6) = (stress+tangent,stress,stress,stress)
c          eps = small strain
c          theta = vol-strain 
c          ud  = user properties [ alpha, beta, gamma, a(3), b(3), c(3) ]

c     Output:
c          sig = Cauchy stresses
c          kk  = spatial tangent

c-----[--.----+----.----+----.-----------------------------------------]



	implicit none

        integer isw, i, j, k, l, m, n
        integer vtot(2,6), ttov(3,3)

        real*8 ud(*)
	real*8 eps(6),sig(6),kk(6,6),theta

        real*8 alpha, beta, gamma, tmpx
        real*8 aaxis(3), baxis(3), caxis(3), iden(6)

        save

        data vtot /1,1, 2,2, 3,3, 1,2, 2,3, 3,1/
        data ttov /1,4,6,4,2,5,6,5,3/
        data iden /1.d0,1.d0,1.d0,0.d0,0.d0,0.d0/

c       Basic Tolerance



        if(isw.eq.3 .or. isw.eq.6 .or. isw.eq.4 .or. 
     &          isw.eq.8 ) then

c       Set up the moduli

        alpha = ud(1)
        beta  = ud(2)
        gamma = ud(3)
        tmpx = alpha - beta - 2.d0*gamma

        aaxis(1) = ud(4)
        aaxis(2) = ud(5)
        aaxis(3) = ud(6)
        baxis(1) = ud(7)
        baxis(2) = ud(8)
        baxis(3) = ud(9)
        caxis(1) = ud(10)
        caxis(2) = ud(11)
        caxis(3) = ud(12)



c       Moduli computation

        call pzero(kk,36)

        do m = 1,6
         i = vtot(1,m)
         j = vtot(2,m)
         do n = 1,6
           l = vtot(1,n)
           k = vtot(2,n)
           kk(m,n) = kk(m,n) + beta*iden(m)*iden(n)
     &             + gamma*(iden(ttov(i,k))*iden(ttov(j,l))
     &             +        iden(ttov(i,l))*iden(ttov(j,k)))
     &             + tmpx*(aaxis(i)*aaxis(j)*aaxis(k)*aaxis(l)  
     &             +       baxis(i)*baxis(j)*baxis(k)*baxis(l)  
     &             +       caxis(i)*caxis(j)*caxis(k)*caxis(l))
         end do
        end do

        
        do i = 1,6
          sig(i) = 0.d0
          do j = 1,6
           sig(i) = sig(i) + kk(i,j)*eps(j)
          end do
        end do

        else
         write(*,*) 'ISW ',isw,' not implemented'
        endif


	end ! subroutine