      module elem20node
!       These constants are read and set in setup(), which should be
!       called just once at the beginning of the FEA analysis program
!      
      implicit none
      save
      integer, parameter :: NCOEFS = 115, NDETS = 1140, NTERMS = 73968
      integer, dimension(NCOEFS) :: dqjac, drjac, dtjac
      integer, dimension(NDETS) :: ix, iy, iz
      integer, dimension(NTERMS) :: jx, jy, jz, jc
      double precision :: det3ijk(18,19,20), cmult(NTERMS)
      end module
      
      program tst
      use elem20node
      implicit none
      integer i,j,npts
      double precision coord(3,20),coef(NCOEFS)
      call setup()
      npts = 1
      do j = 1, npts
         call random_number(coord)
         call Jacobian_deg(coord,coef)
         print 10,(i,coef(i),i=1,NCOEFS)
      end do
   10 format(i3,2x,ES12.4)
      end program

      subroutine setup()
      use elem20node
      implicit none
      integer i, nn, nd
      open(10,file='djac.txt',status='old')
      read(10,*)(dqjac(i),drjac(i),dtjac(i),i=1,NCOEFS)
      close(10)
      
      open(10,file='idx.txt',status='old')
      read(10,*)(ix(i),iy(i),iz(i),i=1,NDETS)
      close(10)
      
      open(10,file='det3s.txt',status='old')
      do i=1,NTERMS
         read(10,*)jx(i),jy(i),jz(i),jc(i),nn,nd
         cmult(i)=dble(nn)/nd
      end do
      close(10)
      return
      end subroutine setup
      
      subroutine Jacobian_deg(coord,coef)
      use elem20node
      implicit none
      integer i
      double precision coord(3,20), coef(115)
      coef(1:115) = 0d0
      do i=1,NDETS
         det3ijk(ix(i),iy(i),iz(i))=det3(ix(i),iy(i),iz(i))
      end do
      do i=1,NTERMS
         coef(jc(i))=coef(jc(i))+det3ijk(jx(i),jy(i),jz(i))*cmult(i)
      end do
      
      return

      contains

      double precision function det3(i,j,k)
      implicit none
      integer i,j,k
      det3 = coord(1,i)*(coord(2,j)*coord(3,k)-coord(2,k)*coord(3,j)) -
     1       coord(2,i)*(coord(1,j)*coord(3,k)-coord(1,k)*coord(3,j)) +
     2       coord(3,i)*(coord(1,j)*coord(2,k)-coord(1,k)*coord(2,j))
      end function det3
      end subroutine Jacobian_deg
