      subroutine hinc(eta_alpha,cv,imgprob)
c********************************************************************
c
c  purpose : to compute the excitation matrix for curved triangles
c            for the mfie. the testing is implemented with gauss int.
c            v = (eta/alpha) * <f,nxh>
c            note: the excitation matrix may be non-zero on input.	
c********************************************************************
c
      implicit complex(c)
c
      include 'ctparm.inc'
      include 'ctgeom.inc'
      include 'ctgeom2.inc'
      include 'symm.inc'

      complex cv(1)
      complex htheta,hphi,cxnxh(3)
      integer n(3)
      dimension r(3),ru(3),rv(3),xn(3),f(3)
      integer imgprob
      real ak, pi, theta, phi
      complex hcart(3,8)
      real khat(3,8), hreal(3,8), himag(3,8)
      real r, argmnt, alamda
      integer iget_images, iget_matrix_position
      data n/3,1,2/,cj/(0.,1.)/
c
      common/kkk/ ak,pi
      common/param/ theta,phi
      common/field/ htheta,hphi,alamda
c
c
      i4pts = 1

      ct1=cmplx(cos(theta),0.0)
      ct2=cmplx(sin(theta),0.0)
      cphi1=cmplx(cos(phi),0.0)
      cphi2=cmplx(sin(phi),0.0)
c
c compute the amplitude of incident magnetic field in
c cartesian coord.
c
      hcart(1,1) = htheta*ct1*cphi1-hphi*cphi2
      hcart(2,1) = htheta*ct1*cphi2+hphi*cphi1
      hcart(3,1) =-htheta*ct2
c
      do k = 1,3
        hreal(k,1) = 1./(2**nplanes)*(real(hcart(k,1)))
      enddo
      numimgs = iget_images(imgprob, hreal(1,1), 3)

      do k = 1,3
        himag(k,1) = 1./(2**nplanes)*(aimag(hcart(k,1)))
      enddo
      numimgs = iget_images(imgprob, himag(1,1), 3)

      do i=1,numimgs
        do j=1,3
          hcart(j,i) = cmplx( hreal(j,i), himag(j,i) )
        end do
      end do

      khat(1,1) = sin(theta)*cos(phi)
      khat(2,1) = sin(theta)*sin(phi)
      khat(3,1) = cos(theta)

      numimgs = iget_images(imgprob, khat(1,1), 1)

      do 500 i=1,numimgs
      do 499 ij=1,nfaces

        call tridata(ij,j2,j3,j4,nj1,nj2,nj3)
c
        do 450 ir=1,3
          if(ir.eq.1) j1=j4
          if(ir.eq.2) j1=j2
          if(ir.eq.3) j1=j3     
          jl=ibkp(j1)			!(still need ibkp for pole holes)
          if(jl.eq.0) go to 450

c  perform gauss integration 
          chi = 0.
          nint_tst = iget_num_int_points(ij,i4pts)		!4 points
          do j = 1,nint_tst
            call get_uv(ij,j,i4pts,u,v)
            call tanvec_uandv(ij,u,v,ru,rv)
            call outnormal(ij,ru,rv,xn)
            call rccrossp(xn,hcart(1,i),cxnxh)

c     set up the basis function
            u_n = tzpt(2*ir,ij)
            v_n = tzpt(2*ir-1,ij)
            do k = 1,3
              f(k) = (u-u_n)*ru(k)+(v-v_n)*rv(k)
            enddo

            call rvec(ij,u,v,r)
            argmnt = dotp(r,khat(1,i))

            chtmp = c8rdotp(cxnxh,f)*cexp(-cj*ak*argmnt)
            wgt_tst = fget_weight(ij,j,i4pts)
            chi = chi + chtmp * wgt_tst
          enddo            
          cv(jl)=cv(jl)+ eta_alpha * chi*aflag(n(ir),ij)/2.
  450   continue
  499 continue
  500 continue
      return
      end

      complex function c8rdotp(ca,b)
c--------------------------------------------------------
c purpose: to compute the dot product of a c*16 3-vector
c          and a r*4 3-vector.
c--------------------------------------------------------
      complex ca(3)
      real*4 b(3)
      c8rdotp=0.
      do k = 1,3
        c8rdotp = c8rdotp + ca(k)*b(k)
      enddo
      return
      end

      subroutine rccrossp(b,ca,cc)
c---------------------------------------------------------
c purpose: to compute cross product of a r*4 3-vector &
c          a c*16 3-vector: cc = b x ca
c---------------------------------------------------------
      complex ca(3),cc(3)
      real*4 b(3)

      cc(1) = b(2)*ca(3) - b(3)*ca(2)
      cc(2) = b(3)*ca(1) - b(1)*ca(3)
      cc(3) = b(1)*ca(2) - b(2)*ca(1)

      return
      end
