      subroutine zmatrx(c1,c2,cz,imgprob)

      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)
      integer imgprob
      real ak, pi
      integer i, j, k, isourc, ifield, isrcnode, ifldnode
      integer numsrcbasis, numfldbasis, isrcmatrixpos, ifldmatrixpos
      integer ipnt, il, jl, np_test
      real uf, vf, ftestdiv, rifv, rifs
      real ftemp(3), ru(3), rv(3), ruu(3), rvv(3), ruv(3)
      real rleft(3), rright(3)
      real vtemp(3),vtester(3)
      complex cvpot(3), cspot, cvtmp, cstmp
      complex cj, ctsta, ctstp
      integer n(3)
      data n/3,1,2/,cj/(0.,1.)/
      integer iget_num_basis, iget_matrix_position
      integer iget_num_int_points
      real fget_weighted_divbasis, vecmag

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

 
c..   loop over source triangle
      do 900 isourc = 1,nfaces
        numsrcbasis = iget_num_basis(isourc)

c--    loop over field triangles...
        do 800 ifield = 1,nfaces
          numfldbasis = iget_num_basis(ifield)

          call patch_interactions(isourc,ifield,imgprob)

c--      loop over source nodes (nodes w/ functions on source patch)
          do 460 isrcnode=1,numsrcbasis
            isrcmatrixpos = iget_matrix_position(isourc,isrcnode)
            if(isrcmatrixpos .eq. 0) goto 460

            if(ntype(ifield) .ne. ep_pattypwire) then
c  source patch is not a wire
              do 470 ifldnode=numfldbasis/2+1, numfldbasis
c  junction testing functions on triangular patch...
c  use moments approximation
                ifldmatrixpos = iget_matrix_position(ifield,ifldnode)
                if(ifldmatrixpos .eq. 0) goto 470

                ipnt = 0
                uf = datctr(1,ifield)
                vf = datctr(2,ifield)
                call pots(ifield,uf,vf,ipnt,isourc,isrcnode,
     &                     cvpot,cspot,imgprob)

                cvtmp = cmplx(0.,0.)
                cstmp = cmplx(0.,0.)
                do k=1,3
                  cvtmp = cvtmp + ctsta(k,ifldnode,ifield)*cvpot(k)
                end do
                cstmp = ctstp(ifldnode,ifield)*cspot

                jl = isrcmatrixpos
                il = ifldmatrixpos

                cz(il,jl) = cz(il,jl) + ( c1*cvtmp - c2*cstmp )

470           continue

              np_test = iget_num_int_points(ifield,1)  !1=for testing
              do ipnt=1,np_test
                call get_uv(ifield,ipnt,1,uf,vf)

                call pots(ifield,uf,vf,ipnt,isourc,isrcnode,
     &                     cvpot,cspot,imgprob)

                do 450 ifldnode=1,numfldbasis/2
                  ifldmatrixpos = iget_matrix_position(ifield,ifldnode)
                  if(ifldmatrixpos .eq. 0) goto 450

                  call get_weighted_basis(ifield,ifldnode,ipnt,1,
     &                     vtester)
                  ftestdiv = fget_weighted_divbasis(ifield,ifldnode,
     &                     ipnt,1)

                  cvtmp = cmplx(0.,0.)
                  cstmp = cmplx(0.,0.)

                  do k=1,3
                    cvtmp = cvtmp + cvpot(k)*vtester(k)
                  end do
                  cstmp = cstmp + cspot*ftestdiv

c
c the z-matrix is stored on a packed form, which assumes that it is
c symmetric. thus, cz(i,j) and cz(j,i) are averaged and stored in the
c upper triangular part of the matrix.
                  jl = isrcmatrixpos
                  il = ifldmatrixpos

                  cz(il,jl) = cz(il,jl) + ( c1*cvtmp - c2*cstmp )

  450           continue
              end do      
            else
              call vecscal(1., datnod(1,nconn_wire(1,ifield)), rleft)
              call vecscal(1., datnod(1,nconn_wire(2,ifield)), rright)
              call vecdif(rright, rleft, ru)
              rifv = vecmag(ru)/2.
              rifs = .5

              np_test = iget_num_int_points(ifield,1)  !1=for testing
              do ipnt=1,np_test
                call get_uv(ifield,ipnt,1,uf,vf)

                call pots(ifield,uf,vf,ipnt,isourc,isrcnode,
     &                     cvpot,cspot,imgprob)

                do 480 ifldnode=1,numfldbasis
                  ifldmatrixpos = iget_matrix_position(ifield,ifldnode)
                  if(ifldmatrixpos .eq. 0) goto 480

                  call get_weighted_basis(ifield,ifldnode,ipnt,1,
     &                     vtester)
                  ftestdiv = fget_weighted_divbasis(ifield,ifldnode,
     &                     ipnt,1)

                  cvtmp = cmplx(0.,0.)
                  cstmp = cmplx(0.,0.)

                  do k=1,3
                    cvtmp = cvtmp + cvpot(k)*vtester(k)
                  end do
                  cstmp = cstmp + cspot*ftestdiv

c
c the z-matrix is stored on a packed form, which assumes that it is
c symmetric. thus, cz(i,j) and cz(j,i) are averaged and stored in the
c upper triangular part of the matrix.
                  jl = isrcmatrixpos
                  il = ifldmatrixpos

                  cz(il,jl) = cz(il,jl) + ( c1*cvtmp*rifv - 
     &                                      c2*cstmp*rifs )

  480           continue
              end do      
            endif
  460     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
      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


