      subroutine ymatrx(eta_alpha,cz,imgprob)
c******************************************************************
c  purpose: to fill the nunknsxnunkns square impedance matrix cz 
c           for the mfie operator.
c           the self term is the j/2 + curvature term. 
c  notes:
c        1. the indices of the z-matrix cz(i,j) are for the i-th
c	 row & j-th column
c
c	 2. this version breaks up the non-singular self-patch
c	 integral into 3 integrals according to the location of 
c        the field point inside the patch
c
c   uses version of basis function that satisfies f.n = 1/|t|
c   (no t-vector terms in f)
c
c******************************************************************
      implicit complex (c)
      parameter (np=7)		!max no of int pts/triangle
      include 'ctparm.inc'
      include 'ctgeom.inc' 
      include 'ctgeom2.inc' 
      dimension cz(mxedgs,mxedgs)
      common/kkk/ ak,pi
      integer n(3)
      dimension rus(3,np),rvs(3,np),us(np),vs(np)
      dimension ruf(3,np),rvf(3,np),uf(np),vf(np)
      dimension dum(3),xnf(3,np),rjac(np)
      dimension r(3),rp(3,8),delr(3,np,np,8)
      complex chint(3),dgterm(np,np,8),crdotp,ctemp(3)
      dimension f_fld(3),f_src(3),fxn(3),fxdelr(3)
      dimension f_tmp(3,8),f_source(3,np,8)
      dimension si1(np), si2(np), si3(np)
      complex cint(3,np)
      data n/3,1,2/,eps/1.e-3/,cj/(0.,1.)/

      i4pts = 1
      i7pts = 2

c.. loop over source triangle

      do 900 ijk = 1,nfaces
        call tridata(ijk,i2,i3,i4,n1,n2,n3)
c..     for each numerical integration point, compute quantities
c       which are only dependent on the source patch.
        nps = iget_num_int_points(ijk,i7pts)	!7 points
        do j = 1,nps
          call get_uv(ijk,j,i7pts,us(j),vs(j))
          call tanvec_uandv(ijk,us(j),vs(j),rus(1,j),rvs(1,j))
        enddo
 
c-- loop over field triangles...

       do 800 ij = 1,nfaces
         call tridata(ij,j2,j3,j4,nj1,nj2,nj3)
         nint_tst = iget_num_int_points(ij,i4pts)	!4 points

c... compute & store terms which only depend on source & field patch
c    combination (different for self vs non-self terms)
         if(ij.eq.ijk) then
c1 only do 1 point testing here
           do i = 1,1  		  	      !(field loop)
             call ac_trans(ij,0.33333,0.33333,uf(i),vf(i))
             call tanvec_uandv(ij,uf(i),vf(i),ruf(1,i),rvf(1,i))
             call outnormal(ij, ruf(1,i), rvf(1,i), xnf(1,i) )
             call crossp( ruf(1,i), rvf(1,i), dum )	!sign of dum
c							!doesn't matter
             rjac(i) = sqrt(dum(1)**2 + dum(2)**2 + dum(3)**2)
           enddo

         else		!ijk not = ij

           do i = 1,nint_tst		!(field loop)
             call get_uv(ij,i,i4pts,uf(i),vf(i))
             call tanvec_uandv(ij,uf(i),vf(i),ruf(1,i),rvf(1,i))
             call outnormal(ij, ruf(1,i), rvf(1,i), xnf(1,i) )
             call rvec(ij,uf(i),vf(i),r)
             do j = 1,nps	!(source loop)
c..           get vector to source point & its images
               call rvec(ijk,us(j),vs(j),rp(1,1))
               numimgs = iget_images(imgprob,rp(1,1),1)
               do img = 1,numimgs
                 bigr = 0.
                 do k = 1,3
                   delr(k,j,i,img) = r(k)-rp(k,img)
                   bigr = bigr + delr(k,j,i,img)**2
                 enddo
                 bigr = sqrt(bigr)
                 carg = cj*ak*bigr
                 dgterm(j,i,img) = (1.+carg)*cexp(-carg)/bigr**3
               enddo	!image  loop
             enddo 	!source loop
           enddo      	!field  loop
        
         endif


c--  loop over nodes of source patch
       do 460 ik=1,3
         i1 = nbound(n(ik),ijk)
         jl=ibkp(i1)
         if(jl.eq.0) go to 460
         u_n = tzpt(2*ik,ijk)
         v_n = tzpt(2*ik-1,ijk)

c--  if ij=ijk (self-term), then precompute the terms which do
c     not depend on the testing basis function
         if(ij.eq.ijk) then
         do i = 1,1
c   compute the surface integral by breaking it into 3 & just doing it
           call sing_intgl3(ijk,xnf(1,i),u_n,v_n,
     %                      uf(i),vf(i),aa(ijk),cint(1,i),imgprob)
         enddo	!(testing point loop)

        else	!ij not = ijk

c.. when the source & field patches are different, compute & store
c   the source basis function and its images
        do j = 1,nps
           do k = 1,3
              f_tmp(k,1) = (us(j)-u_n)*rus(k,j) + (vs(j)-v_n)*rvs(k,j)
           enddo
           numimgs = iget_images(imgprob,f_tmp(1,1),2)
           do i = 1,8
              do k = 1,3
                 f_source(k,j,i) = f_tmp(k,i)
              enddo
           enddo
        enddo

        endif

c--  loop over nodes of testing patch
       do 450 ir=1,3
         if(ir.eq.1) j1=j4
         if(ir.eq.2) j1=j2
         if(ir.eq.3) j1=j3
         il=ibkp(j1)
         if(il.eq..0) go to 450
         u_m = tzpt(2*ir,ij)
         v_m = tzpt(2*ir-1,ij)
c
        if ( ij .eq. ijk )  then		!self term 

          chsct = 0.
          do i = 1,1				!(field loop)
c---   form testing (field) and source basis functions at the field pt
            do k = 1,3
              f_fld(k) = (uf(i)-u_m)*ruf(k,i) + (vf(i)-v_m)*rvf(k,i)
              f_src(k) = (uf(i)-u_n)*ruf(k,i) + (vf(i)-v_n)*rvf(k,i)
            enddo
            ctmp1 = dotp(f_fld,f_src)/(rjac(i)*8.*aa(ijk))
            ctmp2 = c8rdotp(cint(1,i),f_fld)/(16.*pi)
            chsct = chsct +  1. * ( ctmp1 - ctmp2 )	 !(wgt=1)

          enddo
          cz(il,jl) = cz(il,jl) + eta_alpha * chsct * 
     %                            aflag(n(ik),ijk) * aflag(n(ir),ij)

        else

          chsct = 0.

          do i = 1,nint_tst	   !(testing done by num. integration)
c... form basis function on field patch 
            do k = 1,3
              f_fld(k) = (uf(i)-u_m)*ruf(k,i) + (vf(i)-v_m)*rvf(k,i)
              chint(k) = 0.
            enddo
            call crossp(f_fld,xnf(1,i),fxn)

            do j = 1,nps		!(integrate on source patch)
              w = fget_weight(ijk,j,i7pts)
              do img = 1,numimgs
                call crossp(f_source(1,j,img),delr(1,j,i,img),fxdelr)
                do k = 1,3
                   chint(k) = chint(k) + w * fxdelr(k)*dgterm(j,i,img)
                enddo
              enddo
            enddo

            wgt_tst = fget_weight(ij,i,i4pts)
            chsct = chsct + wgt_tst * c8rdotp(chint,fxn)

          enddo
          cz(il,jl) = cz(il,jl) - eta_alpha * chsct * 
     %                aflag(n(ik),ijk) * aflag(n(ir),ij) /(16.*pi)
        endif 
  450   continue
  460  continue
  800 continue
  900 continue
      return
      end

       subroutine outnormal(ipatch,ru,rv,xn)
c-------------------------------------------------------------------
c purpose: to compute the outward normal of the surface given the
c          2 tangent vectors ru & rv and the sign of the normal 
c          where  xn = normsign * ( ru x rv ) / | ru x rv |
c--------------------------------------------------------------------
       include 'ctparm.inc'
       include 'dspparam.inc'
       include 'ctgeom2.inc'
       common /nsign/ normsign(mxobj+mxwire)
       dimension ru(3),rv(3),xn(3)

       call crossp(ru,rv,xn)
       call normze(xn)
       do k = 1,3
         xn(k) = normsign(iparms(ipatch)) * xn(k)
       enddo

       return
       end


      subroutine sing_intgl3(ijk,xn,u_n,v_n,u0,v0,area,ci,imgprob)
c---------------------------------------------------------------------
c purpose: to compute the 'magnetic vector potential integral' term
c          (nxa) when the source & field patches are the same. this is
c          done by sub-dividing the triangle into 3 smaller patches
c          using the field point as a 4th node. then 7 point gauss
c	   quad integration is used in each sub-triangle.
c          for the image terms (when symmetry planes are used), the
c          image terms will not be singular, so only 7 pt integration 
c          in the image triangles is used.
c          all results are summed into ci.
c---------------------------------------------------------------------
       implicit complex (c)
       common /kkk/ ak,pi
       complex ci(3)
       dimension xn(3), ru(3),rv(3),r(3),rp(3,8)
       dimension dr(3)
       dimension f_src(3,8)
       dimension fxdr(3),term(3)
       data cj/(0.,1.)/

       i7pts = 2

       call rvec(ijk,u0,v0,r)

       do k = 1,3
         ci(k) = 0.
       enddo
       xndfv = dotp(xn,fv)

       nps = iget_num_int_points(ijk,i7pts)		!7 points

c--  now divide the triangle into 3 sub-patches & integrate...
       do ll = 1,3		!3 sub triangles
c--  use the field point as the interior node to divide the patch with
         call areacomp(ijk,ll,u0,v0,area3)
         do j = 1,nps		!(integrate on portion of source patch)
            call ac_trans_3tri(ijk,ll,j,u0,v0,u,v)
            call tanvec_uandv(ijk,u,v,ru,rv)
            call rvec(ijk,u,v,rp(1,1))
            call vecdif(r,rp(1,1),dr)
            bigr = sqrt(dotp(dr,dr))
            carg = cj*ak*bigr
            cgarg = (1.+carg)*cexp(-carg)/bigr**3
c         (note: the complete, exact gns fn = cgarg * dr)

c--  form source basis fn at (u,v)
           do k = 1,3
             f_src(k,1) = (u-u_n)*ru(k) + (v-v_n)*rv(k)
           enddo

c--  form the triple cross product
          call crossp(f_src(1,1),dr,fxdr)
          call crossp(xn,fxdr,term)

          w = fget_weight(ijk,j,i7pts)
          do k = 1,3
           ci(k) = ci(k) + w * ( term(k) * cgarg )*area3
          enddo
        enddo		!integration point loop
       enddo		!3 sub_triangle loop


c-- the image terms will now be included...

      do j = 1,nps
         w = fget_weight(ijk,j,i7pts)
         call get_uv(ijk,j,i7pts,u,v)

c        form the position vector to the source point & its images
         call rvec(ijk,u,v,rp(1,1))
         numimgs = iget_images(imgprob,rp(1,1),1)
         
c        form source basis fn at (u,v) along with its images
         call tanvec_uandv(ijk,u,v,ru,rv)
         do k = 1,3
            f_src(k,1) = (u-u_n)*ru(k) + (v-v_n)*rv(k)
         enddo
         numimgs = iget_images(imgprob,f_src(1,1),2)
  
         do img = 2,numimgs		!don't include original here
            call vecdif(r,rp(1,img),dr)
            bigr = sqrt(dotp(dr,dr))
            carg = cj*ak*bigr
            cgarg = (1.+carg)*cexp(-carg)/bigr**3
c           form the triple cross product
            call crossp(f_src(1,img),dr,fxdr)
            call crossp(xn,fxdr,term)
  
            do k = 1,3
             ci(k) = ci(k) + w * ( term(k) * cgarg )*area
            enddo
         enddo		!image loop
       enddo		!source point integration loop


c--  now divide by the patch area (which is part of the basis fn)
c    since this term is not cancled by the area term in the gaussian
c    quadrature integrals.
       do k = 1,3
         ci(k) = ci(k) / area
       enddo

       return
       end


      subroutine areacomp(ipat,ktri,uin,vin,area)
c---------------------------------------------------------------------
c purpose: to compute the area for patch ipat and
c          sub-triangle ktri in (u,v) coordinates.
c          (uin,vin) is a point inside the triangle.
c the 3 sub triangles are:
c	(1) node 1 - node 2 - (u,v)_inside
c       (2) node 2 - node 3 - (u,v)_inside
c	(3) node 3 - node 1 - (u,v)_inside
c
c and the 'third' node in a triangle is the origin in area coordinates
c---------------------------------------------------------------------

      include 'ctparm.inc'
      include 'ctgeom2.inc'
      dimension r0(2), dr_xsi(2), dr_eta(2)

      r0(1) = uin
      r0(2) = vin
      if(ktri.eq.1) then
         dr_xsi(1) = tzpt(2,ipat) - r0(1)
         dr_eta(1) = tzpt(4,ipat) - r0(1)
         dr_xsi(2) = tzpt(1,ipat) - r0(2)
         dr_eta(2) = tzpt(3,ipat) - r0(2)
      elseif(ktri.eq.2) then
         dr_xsi(1) = tzpt(4,ipat) - r0(1)
         dr_eta(1) = tzpt(6,ipat) - r0(1)
         dr_xsi(2) = tzpt(3,ipat) - r0(2)
         dr_eta(2) = tzpt(5,ipat) - r0(2)
      elseif(ktri.eq.3) then
         dr_xsi(1) = tzpt(6,ipat) - r0(1)
         dr_eta(1) = tzpt(2,ipat) - r0(1)
         dr_xsi(2) = tzpt(5,ipat) - r0(2)
         dr_eta(2) = tzpt(1,ipat) - r0(2)
      endif

      area = abs( dr_xsi(1)*dr_eta(2) - dr_xsi(2)*dr_eta(1) )/2.

      return
      end

      subroutine ac_trans_3tri(ipat,ktri,j_point,uin,vin,u,v)
c---------------------------------------------------------------------
c purpose: to convert a point in the local area coordinates (xsi,eta)
c (which corresponds to point j_point in the 7-point integration scheme)
c for patch ipat and sub-triangle ktri into (u,v) coordinates.
c
c the 3 sub triangles are:
c	node 1 - node 2 - (u,v)_inside
c       node 2 - node 3 - (u,v)_inside
c	node 3 - node 1 - (u,v)_inside
c
c and the 'third' node in a triangle is the origin in area coordinates
c---------------------------------------------------------------------
      include 'ctparm.inc'
      include 'ctgeom2.inc'
      include 'gauss.inc'
      dimension r0(2), dr_xsi(2), dr_eta(2)
      i7pts = 2

      xsi = evalptsx( ipat, i7pts, j_point )
      eta = evalptse( ipat, i7pts, j_point )

      r0(1) = uin
      r0(2) = vin
      if(ktri.eq.1) then
         dr_xsi(1) = tzpt(2,ipat) - r0(1)
         dr_eta(1) = tzpt(4,ipat) - r0(1)
         dr_xsi(2) = tzpt(1,ipat) - r0(2)
         dr_eta(2) = tzpt(3,ipat) - r0(2)
      elseif(ktri.eq.2) then
         dr_xsi(1) = tzpt(4,ipat) - r0(1)
         dr_eta(1) = tzpt(6,ipat) - r0(1)
         dr_xsi(2) = tzpt(3,ipat) - r0(2)
         dr_eta(2) = tzpt(5,ipat) - r0(2)
      elseif(ktri.eq.3) then
         dr_xsi(1) = tzpt(6,ipat) - r0(1)
         dr_eta(1) = tzpt(2,ipat) - r0(1)
         dr_xsi(2) = tzpt(5,ipat) - r0(2)
         dr_eta(2) = tzpt(1,ipat) - r0(2)
      endif
      u = r0(1) + xsi * dr_xsi(1)  +  eta * dr_eta(1) 
      v = r0(2) + xsi * dr_xsi(2)  +  eta * dr_eta(2) 
c      u = tzpt(6,ipat) + xsi * ( tzpt(2,ipat) - tzpt(6,ipat) ) +
c     %                   eta * ( tzpt(4,ipat) - tzpt(6,ipat) )
c      v = tzpt(5,ipat) + xsi * ( tzpt(1,ipat) - tzpt(5,ipat) ) +
c     %                   eta * ( tzpt(3,ipat) - tzpt(5,ipat) )

      return
      end
