

      subroutine pots_junc_t( ifield, utest, vtest, ifldpnt, 
     &                        isourc, inode, cvpot, cspot, imgprob)

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

      integer ifield, ifldpnt, isourc, inode, imgprob
      real utest, vtest
      complex cvpot(3), cspot, cftmp(3)
      complex cvtemp(8), cstemp(8), cfac1
      complex cj, cg
      real ru(3), rv(3), ruu(3), rvv(3), ruv(3)
      real get_r, get_rs, eps, ftemp(3), rhohat(3,8)
      real ak, pi
      integer k, ipnt, np_pot
      integer numimgs, imgnum
      real rsint, terml, xsi_f, eta_f, g, rho, rho_s
      real fac2, u_n, v_n, u, v, ri, rjac, rsint2, us, vs
      real xsi, eta, rs, r(8), w, termlr
      real posfield(3), posvects(3,8)
      data eps/1.e-3/
      integer iget_num_int_points
      real get_rho, get_rho_s, get_l, vecmag, fget_weight, min
      integer iget_images

      common /kkk/ ak, pi


c the source patch is triangular, the field patch may be either a wire
c  or a triangle.

c initialize...
      cj = cmplx(0.,1.)
      do k = 1,3
        cvpot(k) = 0.
      enddo
      cspot = 0.
      
      do k=1,8
        cvtemp(k) = 0.
        cstemp(k) = 0.
      end do

      call rvec(ifield,utest,vtest,posfield)
 
c   if source=field, the field patch is then known to be a triangle.
c   need to remove the singularity.
      if(isourc .eq. ifield) then  
        call sigint(isourc, utest, vtest, rsint)
        termlr = get_l(isourc,utest,vtest,inode)

        call uv_trans(isourc,utest,vtest,xsi_f,eta_f)
        if(inode .eq. 4) then
          g = 1. - (1.-xsi_f)*(1.-xsi_f)
        elseif(inode .eq. 5) then
          g = 1. - (1.-eta_f)*(1.-eta_f)
        elseif(inode .eq. 6) then
          g = 1. - (1.-(1.-xsi_f-eta_f))*(1.-(1.-xsi_f-eta_f))
        else
          write(6,*) '*** internal error ***'
          write(6,*) 'wrong basis reference in pots_junc_t'
          stop
        endif
        rho = get_rho(isourc,inode,utest,vtest)

        fac2 = g/rho
      endif

      if(inode .gt. 3) then
        u_n = tzpt(2*(inode-3),isourc)
        v_n = tzpt(2*(inode-3)-1,isourc)
        ri = get_r(isourc,u_n,v_n,ifield,utest,vtest)
        cfac1 = cexp(-cj*ak*ri)/ri

        u = datctr(1,isourc)
        v = datctr(2,isourc)
        call tanvecs(isourc,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,1) = ftemp(k) / vecmag(ftemp)
        end do
        numimgs = iget_images(imgprob,rhohat(1,1),2)

        call tanvecs(isourc,u_n,v_n,ru,rv,ruu,rvv,ruv)
        call crossp(ru,rv,ftemp)
        rjac =  vecmag(ftemp)      !jacobian at r1 

        call sigint(isourc, u_n, v_n, rsint2)

      else
          write(6,*) '*** internal error ***'
          write(6,*) 'wrong basis reference in pots_junc_t'
         stop
      endif

c perform gaussian integration...
      np_pot = iget_num_int_points(isourc,2) !2=for potentials
      do ipnt = 1,np_pot
        call get_uv(isourc,ipnt,2,us,vs)

        call rvec(isourc,us,vs,posvects(1,1))
        numimgs = iget_images(imgprob,posvects(1,1),1)

        call uv_trans(isourc,us,vs,xsi,eta)

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

        rho = get_rho(isourc,inode,us,vs)
        rho_s = get_rho_s(isourc,inode,us,vs)
        w = fget_weight(isourc,ipnt,2)
        terml = get_l(isourc,us,vs,inode)
      
        call get_r_mult(posfield,posvects(1,1),numimgs,r)

        do imgnum=1,numimgs

          if(imgnum .eq. 1) then
            if (isourc .eq. ifield) then
              rs = get_rs(isourc,us,vs,utest,vtest)
              if( ak*min(r(imgnum),rs) .lt. eps ) then
                cvtemp(imgnum) = cvtemp(imgnum) + 
     &                           w*(-g*cj*ak/rho - cfac1/rho_s)
                cstemp(imgnum) = cstemp(imgnum) + w*(-terml*cj*ak)
              else
                cg = cexp(-cj*ak*r(imgnum))/r(imgnum)
                cvtemp(imgnum) = cvtemp(imgnum) + 
     &                           w*(g*cg/rho - cfac1/rho_s - fac2/rs) 
                cstemp(imgnum) = cstemp(imgnum) + 
     &                           w*(terml*cg - termlr/rs)
              endif
            else
              cg = cexp(-cj*ak*r(imgnum))/r(imgnum)
              cvtemp(imgnum) = cvtemp(imgnum) + 
     &                         w*(g*cg/rho - cfac1/rho_s)
              cstemp(imgnum) = cstemp(imgnum) + w*(terml*cg)
            endif
          else
            cg = cexp(-cj*ak*r(imgnum))/r(imgnum)
            cvtemp(imgnum) = cvtemp(imgnum) + 
     &                       w*(g*cg/rho - cfac1/rho_s)
            cstemp(imgnum) = cstemp(imgnum) + w*(terml*cg)
          endif
        end do
      end do
      
c add in singular integral term...
      if(isourc.eq.ifield) then
        do imgnum=1,numimgs
          do k=1,3
            cvpot(k) = rhohat(k,imgnum)*rjac/alpha(inode-3,isourc) * 
     &                 (cvtemp(imgnum)*aa(isourc) + cfac1*rsint2 + 
     &                             fac2*rsint)
          end do
          cspot = rjac/alpha(inode-3,isourc) * 
     &                 (cstemp(imgnum)*aa(isourc) + termlr*rsint)
        end do
      else
        do imgnum=1,numimgs
          do k=1,3
            cvpot(k) = rhohat(k,imgnum)*rjac/alpha(inode-3,isourc) * 
     &                 (cvtemp(imgnum)*aa(isourc) + cfac1*rsint2)
          end do
          cspot = rjac/alpha(inode-3,isourc) * 
     &            cstemp(imgnum)*aa(isourc)
        end do
      endif


      return
      end





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

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

      integer ipatch, ibasref, k
      real u, v, u_n, v_n
      real ru(3),rv(3),ruu(3),rvv(3),ruv(3)
      real vec(3)
      real dotp

      call tanvecs(ipatch,u,v,ru,rv,ruu,rvv,ruv)
      if(ibasref .gt. 3) then
        u_n = tzpt(2*(ibasref-3),ipatch)
        v_n = tzpt(2*(ibasref-3)-1,ipatch)
      else
        u_n = tzpt(2*ibasref,ipatch)
        v_n = tzpt(2*ibasref-1,ipatch)
      endif

      do k=1,3
        vec(k) = (u-u_n)*ru(k) + (v-v_n)*rv(k)
      end do

      get_rho = sqrt(dotp(vec,vec))

      return
      end



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

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

      integer ipatch, ibasref
      real u, v, u_n, v_n, hu, hv
      real ru(3),rv(3),ruu(3),rvv(3),ruv(3)
      real dotp

      if(ibasref .gt. 3) then
        u_n = tzpt(2*(ibasref-3),ipatch)
        v_n = tzpt(2*(ibasref-3)-1,ipatch)
      else
        u_n = tzpt(2*ibasref,ipatch)
        v_n = tzpt(2*ibasref-1,ipatch)
      endif

      call tanvecs(ipatch,u_n,v_n,ru,rv,ruu,rvv,ruv)
      hu = sqrt(dotp(ru,ru))
      hv = sqrt(dotp(rv,rv))

      get_rho_s = sqrt( ((u-u_n)*hu)**2  +
     &                  (u-u_n)*(v-v_n)*dotp(ru,rv)  +
     &                  ((v-v_n)*hv)**2 )

      return
      end




      subroutine uv_trans(ipatch,u,v,xsi,eta)

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

      integer ipatch
      real u, v, xsi, eta
      real u1, v1, u2, v2, u3, v3, twoa
      real utest, vtest

      u1 = tzpt(2,ipatch)
      v1 = tzpt(1,ipatch)

      u2 = tzpt(4,ipatch)
      v2 = tzpt(3,ipatch)

      u3 = tzpt(6,ipatch)
      v3 = tzpt(5,ipatch)

      twoa = (u1-u3)*(v2-v3) - (u2-u3)*(v1-v3)

      xsi = 1./twoa * ( (u-u3)*(v2-v3) - (v-v3)*(u2-u3) )
      eta = 1./twoa * (-(u-u3)*(v1-v3) + (v-v3)*(u1-u3) )

      return
      end




      real function get_l(isourc,u,v,inode)
      
      include 'ctparm.inc'
      include 'ctgeom2.inc'

      integer isourc, inode
      real u, v, u1, v1, u2, v2, u3, v3, twoa
      real g, ac, dacdu, dacdv, rho, gamma_u, gamma_v
      real xsi, eta, u_n, v_n, drhodu, drhodv
      real get_rho

      u1 = tzpt(2,isourc)
      v1 = tzpt(1,isourc)
      u2 = tzpt(4,isourc)
      v2 = tzpt(3,isourc)
      u3 = tzpt(6,isourc)
      v3 = tzpt(5,isourc)
      twoa = (u1-u3)*(v2-v3) - (u2-u3)*(v1-v3)

      call uv_trans(isourc,u,v,xsi,eta)

      if(inode .eq. 4) then
        g = 1. - (1.-xsi)*(1.-xsi)
        ac = xsi
        dacdu =  1./twoa * (v2-v3)
        dacdv = -1./twoa * (u2-u3)
      elseif(inode .eq. 5) then
        g = 1. -(1.-eta)*(1.-eta)
        ac = eta
        dacdu = -1./twoa * (v1-v3)
        dacdv =  1./twoa * (u1-u3)
      elseif(inode .eq. 6) then
        g = 1. - (1.-(1.-xsi-eta))*(1.-(1.-xsi-eta))
        ac = 1. - xsi - eta
        dacdu = -1./twoa * (v2-v1)
        dacdv =  1./twoa * (u2-u1)
      endif


      if(inode .gt. 3) then
        u_n = tzpt(2*(inode-3),isourc)
        v_n = tzpt(2*(inode-3)-1,isourc)
      else
        write(6,*) '*** internal error ***'
        write(6,*) 'wrong basis index in get_l'
        stop
      endif

      rho = get_rho(isourc,inode,u,v)
      call rhoderiv(isourc,inode,u,v,drhodu,drhodv)

      gamma_u = (u-u_n)*( 2.*dacdu*(1.-ac)/(ac*(2.-ac)) +
     &                    1./(u-u_n) - 2.*drhodu/rho )

      gamma_v = (v-v_n)*( 2.*dacdv*(1.-ac)/(ac*(2.-ac)) +
     &                    1./(v-v_n) - 2.*drhodv/rho )

      get_l = g/(rho*rho) * (gamma_u + gamma_v)


      return
      end




      subroutine rhoderiv(isourc,inode,u,v,drhodu,drhodv)

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

      integer isourc, inode
      real u, v, drhodu, drhodv, u_n, v_n
      real deltu, deltv, rho
      real ru(3),rv(3),ruu(3),rvv(3),ruv(3)
      real get_rho, dotp


      u_n = tzpt(2*(inode-3),isourc)
      v_n = tzpt(2*(inode-3)-1,isourc)

      deltu = (u-u_n)
      deltv = (v-v_n)
      rho = get_rho(isourc,inode,u,v)

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

      drhodu = 1./rho * ( deltu*dotp(ru,ru) +
     &                    deltu*deltu*dotp(ruu,ru) +
     &                    deltv*deltv*dotp(ruv,rv) +
     &                    deltv*dotp(ru,rv) +
     &                    deltu*deltv*dotp(ruu,rv) +
     &                    deltu*deltv*dotp(ru,ruv) )

      drhodv = 1./rho * ( deltv*dotp(rv,rv) +
     &                    deltv*deltv*dotp(rvv,rv) +
     &                    deltu*deltu*dotp(ruv,ru) +
     &                    deltu*dotp(rv,ru) +
     &                    deltu*deltv*dotp(rvv,ru) +
     &                    deltu*deltv*dotp(rv,ruv) )

      return
      end
