      subroutine zmatrx(c1,c2,cz,imgprob)
      implicit none

      include 'ctparm.inc'
      include 'ctgeom.inc' 
      include 'ctgeom2.inc' 
      include 'gauss.inc'
      include 'ctgeomwire.inc'
      include 'ctgeomjunc.inc'

      complex c1, c2
      complex cz(mxedgs,mxedgs)
      complex cvpot(3), cspot, cvtemp
      complex ctsta, ctstp
      integer imgprob
      real vtester(3), ftestdiv, rleft(3), rright(3), ru(3)
      real ak, pi,cj
      integer iedge2node(3)
      data iedge2node/2,3,1/
      data cj/(0.0,1.0)/
      logical*1 ibasistype(mxedgs)
      integer inumpat(mxedgs),ipatlist(mxedgs,10)
c     integer irelbasis(mxfaces,mxedgs)
      integer isrcnode,ifldnode,i,j,k,m,n,isp,ifp,ifpnt
      integer iget_num_basis,iget_matrix_position,iget_num_int_points
      integer isrcpatch,ifldpatch,ep_pattypewire
      real utest,vtest,rifv,rifs,vecmag,ftestdiv
      real fget_weighted_divbasis

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

      do i=1,nfaces
        m = iget_num_basis(i)
        n = m/2
        do j=1,n
           k = iget_matrix_position(i,j)
           if(k .gt. 0) then
              ibasistype(k) = .false.
              inumpat(k) = inumpat(k) + 1
              ipatlist(k,inumpat(k)) = i           
           endif
        end do
        do j=n+1,m
           k = iget_matrix_position(i,j)
           if(k .gt. 0) then
              ibasistype(k) = .true.
              inumpat(k) = inumpat(k) + 1
              ipatlist(k,inumpat(k)) = i           
           endif
        end do
      end do
      
c      do i=1,nunkns
c        do j=1,nfaces
c          do k=1,iget_num_basis(j)
c            if(iget_matrix_position(j,k) .eq. i) then
c              irelbasis(j,i) = k
c            endif
c          end do
c        end do
c      end do


      do 900 i=1,nunkns
         do 800 j=1,nunkns
            do 700 isp=1,inumpat(i)
               isrcpatch = ipatlist(i,isp)
               do k=1,iget_num_basis(isrcpatch)
                  if(iget_matrix_position(isrcpatch,k).eq.i) then
                     isrcnode = k
                     goto 111
                  endif
               end do
 111           continue
               
               do 600 ifp=1,inumpat(j)
                  ifldpatch = ipatlist(j,ifp)
                  do k=1,iget_num_basis(ifldpatch)
                     if(iget_matrix_position(ifldpatch,k) .eq. j) then
                        ifldnode = k
                        goto 222
                     endif
                  end do
 222              continue
                  
                  if(ntype(ifldpatch) .ne. EP_PATTYPEWIRE  .and.
     &                 ibasistype(j)) then
c     /* field patch is a triangle and the testing    *
c     * function is a junction , use moments approx. */
                     ifpnt = 0
                     utest = datctr(1,ifldpatch)
                     vtest = datctr(2,ifldpatch)
                     call pots(ifldpatch,utest,vtest,ifpnt,
     &                    isrcpatch,isrcnode,
     &                    cvpot,cspot,imgprob)
                     
cg    ifldnode = irelbasis(ifldpatch,j)
                     cvtemp = ctsta(1,ifldnode,ifldpatch)*cvpot(1) +
     &                    ctsta(2,ifldnode,ifldpatch)*cvpot(2) +
     &                    ctsta(3,ifldnode,ifldpatch)*cvpot(3)
                     
                     cz(j,i) = cz(j,i) + ( c1*cvtemp - 
     &                    c2*ctstp(ifldnode,ifldpatch)*cspot )
                     
                  elseif(ntype(ifldpatch) .eq. EP_PATTYPWIRE) then
c     /* testing on a wire "patch" */
                     call vecscal(1., datnod(1,
     &                    nconn_wire(1,ifldpatch)),rleft) 
                     call vecscal(1., 
     &                    datnod(1,nconn_wire(2,ifldpatch)),rright) 
                     call vecdif(rright, rleft, ru)
                     rifv = vecmag(ru)/2.
                     rifs = .5
                     
                     do 490 ifpnt=1,iget_num_int_points(ifldpatch,1)
                        call get_uv(ifldpatch,ifpnt,1,utest,vtest)
                        call pots(ifldpatch,utest,vtest,ifpnt,
     &                       isrcpatch,isrcnode,
     &                       cvpot,cspot,imgprob)
                        call get_weighted_basis(ifldpatch,
     &                       ifldnode,ifpnt,1,
     &                       vtester)
                        ftestdiv = fget_weighted_divbasis(ifldpatch,
     &                       ifldnode,ifpnt,1)
                        
                        cvtemp = vtester(1)*cvpot(1) +
     &                       vtester(2)*cvpot(2) +
     &                       vtester(3)*cvpot(3)
                        
                        cz(j,i) = cz(j,i) + ( c1*cvtemp*rifv - 
     &                       c2*cspot*ftestdiv*rifs )
                        
 490                 continue      
                  else
                     do 500 ifpnt=1,iget_num_int_points(ifldpatch,1)
                        call get_uv(ifldpatch,ifpnt,1,utest,vtest)
                        call pots(ifldpatch,utest,vtest,ifpnt,
     &                       isrcpatch,isrcnode,
     &                       cvpot,cspot,imgprob)
                        call get_weighted_basis(ifldpatch,
     &                       ifldnode, ifpnt,1,vtester)
                        ftestdiv = fget_weighted_divbasis(ifldpatch,
     &                       ifldnode, ifpnt,1)
                        cvtemp = cvpot(1)*vtester(1) + 
     &                       cvpot(2)*vtester(2) +
     &                          cvpot(3)*vtester(3)
                        
                        cz(j,i) = cz(j,i)+c1*cvtemp-
     &                       c2*cspot*ftestdiv
 500                 continue
                  endif
 600           continue
 700        continue
 800     continue
 900  continue
         
      return
      end


      integer function iget_matrix_position( ipatch, ibasisref )

      include 'ctparm.inc'
      include 'ctgeom.inc'
      include 'ctgeom2.inc'
      include 'ctgeomwire.inc'
      include 'ctgeomjunc.inc'

      integer ipatch, ibasisref, ipos, iedge
      integer inode2edge(3)
      data inode2edge /3,1,2/

      if(ntype(ipatch) .ne. ep_pattypwire) then
        if(ibasisref .gt. 3) then
          ipos = ijuncmap(ibasisref-3,ipatch)
        else
          iedge = inode2edge(ibasisref)
          ipos = ibkp( nbound(iedge,ipatch) )
        endif
      else
        if(ibasisref .gt. 2) then
          ipos = ijuncmap(ibasisref-2,ipatch)
        else
          ipos = ibkp_wire( nconn_wire(ibasisref,ipatch) )
        endif
      endif

      iget_matrix_position = ipos

      return
      end



      integer function iget_num_int_points(ipatch,itype)

      include 'ctparm.inc'
      include 'gauss.inc'

      integer ipatch, itype

      iget_num_int_points = numintpts(ipatch,itype)

      return
      end


      

      subroutine get_uv(ipatch,int_pt_ref,itype,u,v)

      include 'ctparm.inc'
      include 'ctgeom.inc'
      include 'ctgeom2.inc'
      include 'gauss.inc'

      integer ipatch, int_pt_ref, itype
      real u, v
      real xsi, eta

      if(int_pt_ref .eq. 0) then
c centroid of patch...
c datctr hols u,v coordinates of centroid
        u = datctr(1,ipatch)
        v = datctr(2,ipatch)
      else
        xsi = evalptsx(ipatch,itype,int_pt_ref)
        eta = evalptse(ipatch,itype,int_pt_ref)

        if(ntype(ipatch) .ne. ep_pattypwire) then
          call ac_trans(ipatch,xsi,eta,u,v)
        else
          u = (xsi+1.)/2.
          v = 0.
        endif
      endif

      return
      end





      subroutine pots(ifield, uf, vf, ifldpnt, isourc, ibasisref,
     &              cvpot, cspot, imgprob)

      include 'ctparm.inc'
      include 'ctgeom2.inc'

      integer ifield, ifldpnt, isourc, ibasisref, imgprob
      real uf, vf
      complex cvpot(3), cspot
cg      real ru(3)

      if(ntype(isourc) .ne. ep_pattypwire) then
        if(ibasisref .gt. 3) then
          call pots_junc_t(ifield,uf,vf,ifldpnt,
     &                     isourc,ibasisref,
     &                     cvpot,cspot,imgprob)
        else
          call pots_tria(ifield,uf,vf,ifldpnt,
     &                   isourc,ibasisref,
     &                   cvpot,cspot,imgprob)
        endif
      else
        call pots_wire(ifield,uf,vf,isourc,ibasisref,
     &                 cvpot,cspot,imgprob)
      endif

      return
      end




      subroutine get_basis(ipatch,ibasref,u,v,basis)

      include 'ctparm.inc'
      include 'ctgeom.inc'
      include 'ctgeom2.inc'

      integer ipatch, ibasref
      real u, v, basis(3)
      integer k
      real u_n, v_n, fact
      real ru(3), rv(3), ruu(3), rvv(3), ruv(3)
      integer n(3)
      data n/3,1,2/
      real vecmag

      call tanvecs(ipatch,u,v,ru,rv,ruu,rvv,ruv)

      if(ntype(ipatch) .ne. ep_pattypwire) then
        if(ibasref .gt. 3) then
          write(6,*) '*** internal error ***'
          write(6,*) 'wrong basis function reference in get_basis'
          write(6,*) 'cannot return junction function value'
          stop
        endif
        u_n = tzpt(2*ibasref,ipatch)
        v_n = tzpt(2*ibasref-1,ipatch)
        do k=1,3
          basis(k) = (u-u_n)*ru(k) + (v-v_n)*rv(k)
          basis(k) = basis(k) * aflag(n(ibasref),ipatch)/2.
        end do
      else
        if(ibasref .eq. 2  .or.  ibasref .eq. 4) then
          do k=1,3
            basis(k) = u*ru(k) / vecmag(ru)
          end do
        elseif(ibasref .eq. 1  .or.  ibasref .eq. 3) then
          if(ibasref .eq. 3) fact=-1.
          if(ibasref .eq. 1) fact=1.
          do k=1,3
            basis(k) = fact * (1.-u)*ru(k) / vecmag(ru)
          end do
        else
          write(6,*) '*** internal error ***'
          write(6,*) 'unknown basis reference ',ibasref
          write(6,*) 'in get_basis'
          stop
        endif
      endif

      return
      end



      real function fget_divbasis(ipatch,ibasref,u,v)

      include 'ctparm.inc'
      include 'ctgeom.inc'
      include 'ctgeom2.inc'

      integer ipatch, ibasref
      real u, v
      real ru(3), rv(3), ruu(3), rvv(3), ruv(3)
      integer n(3)
      data n/3,1,2/

      call tanvecs(ipatch,u,v,ru,rv,ruu,rvv,ruv)

      if(ntype(ipatch) .ne. ep_pattypwire) then
        if(ibasref .gt. 3) then
          write(6,*) '*** internal error ***'
          write(6,*) 'wrong basis reference in fget_divbasis'
          write(6,*) 'cannot return divergence of junction function'
          stop
        endif
        fget_divbasis = 1. * aflag(n(ibasref),ipatch)
      else
        if(ibasref .eq. 1) then
          fget_divbasis = -1.
        elseif(ibasref .eq. 2) then
          fget_divbasis = 1.
        elseif(ibasref .eq. 3) then
          fget_divbasis = 1.
        elseif(ibasref .eq. 4) then
          fget_divbasis = 1.
        else
          write(6,*) '*** internal error ***'
          write(6,*) 'unknown basis reference ',ibasref
          write(6,*) 'in fget_divbasis'
          stop
        endif
      endif

      return
      end




      subroutine get_weighted_basis(ipatch,ibasref,intref,itype,
     &                             basis)

      include 'ctparm.inc'
      include 'speed.inc'

      integer ipatch, ibasref, intref, itype, k
      real basis(3), w
      real fget_weight

      w = fget_weight(ipatch,intref,itype)

      do k=1,3
        basis(k) = w * bas_f(k,ipatch,intref,itype,ibasref)
      end do

      return
      end


      real function fget_weighted_divbasis(ipatch,ibasref,intref,itype)

      include 'ctparm.inc'
      include 'speed.inc'

      integer ipatch, ibasref, intref, itype
      real w
      real fget_weight

      w = fget_weight(ipatch,intref,itype)

      fget_weighted_divbasis = w * divbas_f(ipatch,intref,itype,ibasref)

      return
      end


      real function fget_weight(ipatch,intref,itype)

      include 'ctparm.inc'
      include 'gauss.inc'

      integer ipatch, intref, itype

      fget_weight = weights(ipatch,itype,intref)

      return
      end





      integer function iget_num_basis(ipatch)

      include 'ctparm.inc'
      include 'ctgeom2.inc'

      integer ipatch, iret

      if(ntype(ipatch) .eq. ep_pattypwire) then
        iret = 4
      else
        iret = 6
      endif

      iget_num_basis = iret

      return
      end


      subroutine rvec_indexed(ipat,int_pt_ref,itype,r)

      include 'ctparm.inc'
      include 'speed.inc'

      integer ipat, int_pt_ref,itype
      real r(3)

      r(1) = r_vect(1,ipat,int_pt_ref,itype)
      r(2) = r_vect(2,ipat,int_pt_ref,itype)
      r(3) = r_vect(3,ipat,int_pt_ref,itype)


      return
      end
      
            
