      subroutine gentesters
c********************************************************************
c purpose: to compute the testing functions.
c     ctsta(k,i,p) is the integral of the basis function for the i-th
c                  edge in the p-th patch & k=x,y,z components.
c     ctstp(i,p)   is the integral of the divergence of the basis 
c                  function of the i-th edge in the p-th patch.
c
c  note: these terms do contain the approiate scale factors 
c        (ie, including aflag)
c********************************************************************
c
      include 'ctparm.inc'
      include 'ctgeom.inc'
      include 'ctgeom2.inc' 
      include 'gauss.inc'
      include 'ctgeomjunc.inc'
      include 'symm.inc'

      real ak, pi
      complex ctsta, ctstp
      real ru_wire(3),rv_wire(3)
      real wf(3), wdf, vtemp(3),ftemp(3)
      real ru(3),rv(3),ruu(3),rvv(3),ruv(3),rhohat(3)
      dimension ta(3)
      integer ijk, k, np, numbas, ik, jl, i
      real ta, tp, u_n, v_n, u, v, rjac, rsint, xsi, eta, g
      real rho, rho_s, w, terml
      complex cvtemp, cstemp
      integer n(3)
      data n/3,1,2/
      integer iget_num_int_points, iget_num_basis
      integer iget_matrix_position
      real fget_weighted_divbasis, vecmag, get_rho, get_rho_s
      real fget_weight, get_l

      common/tst/ ctsta(3,6,mxfaces),ctstp(6,mxfaces)
      common/kkk/ ak,pi


c.. loop over source=field triangle

      do 900 ijk = 1,nfaces
       if(ntype(ijk) .eq. ep_pattypwire) then
         call tanvecs(ijk, .5, 0., ru_wire, rv_wire, ruu,rvv,ruv)
         call vecscal(.5,ru_wire,vtemp)
         do k=1,3
           ctsta(k,1,ijk) = vtemp(k)
           ctsta(k,2,ijk) = vtemp(k)
           ctsta(k,3,ijk) = -vtemp(k)
           ctsta(k,4,ijk) = vtemp(k)
         end do
         ctstp(1,ijk) = -1.
         ctstp(2,ijk) = 1.
         ctstp(3,ijk) = 1.
         ctstp(4,ijk) = 1.
       else
        np = iget_num_int_points(ijk,2) !use potential integration points
        numbas = iget_num_basis(ijk)

        do 460 ik=1,numbas
          jl=iget_matrix_position(ijk,ik)
          if( (ik .gt. 3  .and.  jl .eq. 0) ) go to 460

          if(ik .le. 3) then
            do k = 1,3
              ta(k) = 0.
            enddo
            tp = 0.

c perform gaussian integration...
            do i = 1,np
c compute the basis function & its divergence 
              call get_weighted_basis(ijk,ik,i,2,wf)
              wdf = fget_weighted_divbasis(ijk,ik,i,2)

              do k = 1,3
                ta(k) = ta(k) + wf(k)
              enddo
              tp = tp + wdf
            enddo

c note: the 1/area term in f & div(f) cancels with the area
c       factor needed in the change of coord in the gauss integration
 
              do k = 1,3
                ctsta(k,n(ik),ijk) = ta(k)
              enddo
              ctstp(n(ik),ijk) = tp
            else
              u_n = tzpt(2*(ik-3),ijk)
              v_n = tzpt(2*(ik-3)-1,ijk)
              u = datctr(1,ijk)
              v = datctr(2,ijk)
              call tanvecs(ijk,u,v,ru,rv,ruu,rvv,ruv)
              do k=1,3
                ftemp(k) = (u-u_n)*ru(k) + (v-v_n)*rv(k)
              end do
              do k=1,3
                rhohat(k) = ftemp(k)/vecmag(ftemp)
              end do

              call tanvecs(ijk,u_n,v_n,ru,rv,ruu,rvv,ruv)
              call crossp(ru,rv,ftemp)
              rjac = vecmag(ftemp)

              call sigint(ijk,u_n,v_n,rsint)

              cvtemp = 0.
              cstemp = 0.

c perform gaussian integration...
              do i = 1,np
                call get_uv(ijk,i,2,u,v)
                call uv_trans(ijk,u,v,xsi,eta)

                if(ik .eq. 4) then
                  g = 1. - (1.-xsi)*(1.-xsi)
                elseif(ik .eq. 5) then
                  g = 1. - (1.-eta)*(1.-eta)
                elseif(ik .eq. 6) then
                  g = 1. - (xsi+eta)*(xsi+eta)
                else
                  write(6,*) '*** internal error ***'
                  write(6,*) 'wrong basis index in gentesters'
                  stop
                endif

                rho = get_rho(ijk,ik,u,v)
                rho_s = get_rho_s(ijk,ik,u,v)
                w = fget_weight(ijk,i,2)
                terml = get_l(ijk,u,v,ik)

                cvtemp = cvtemp + w*(g/rho - 1./rho_s)
                cstemp = cstemp + w*terml
              enddo

              do k=1,3
                ctsta(k,ik,ijk) = rhohat(k) * rjac/alpha(ik-3,ijk) *
     &                         (cvtemp*aa(ijk) + rsint)
              end do 
              ctstp(ik,ijk) = cstemp * rjac/alpha(ik-3,ijk)*aa(ijk)
            endif
 460      continue
        endif
 900  continue


      return
      end      

