

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

      include 'ctparm.inc'
      include 'ctgeom.inc'
      include 'speed.inc'
      include 'symm.inc'

      integer ifield, ifldpnt, isourc, inode, imgprob
      real utest, vtest
      complex cvpot(3), cspot, cftmp(3)
      complex cj, cg1
      real rfield(3),basisvects(3,8),basistest(3)
      real divbasis(8)
      real ru(3), rv(3), ruu(3), rvv(3), ruv(3)
      real get_rs, eps, ftemp(3)
      real rsint, us, vs, rs, ak, pi, w
      integer k, np_pot, ipnt, numimgs, imgnum, i
      real divbasistest
      real fget_weighted_divbasis, fget_divbasis, fget_weight
      integer iget_num_int_points, iget_div_images
      integer iget_images

      data eps/1.e-3/
      common /kkk/ ak, pi


      cj = cmplx(0.,1.)
      do k = 1,3
         cvpot(k) = 0.
      enddo
      cspot = 0.

      if(isourc .eq. ifield) then  
        call sigint(isourc, utest, vtest, rsint)
      endif

      np_pot = iget_num_int_points(isourc,2)
      do ipnt = 1,np_pot

c weighted basis function vector at source point
        call get_weighted_basis(isourc,inode,ipnt,2,basisvects(1,1)) 
c divergence of basis function at source point
        divbasis(1) = fget_weighted_divbasis(isourc,inode,ipnt,2)

        numimgs = 1
c images of weighted basis function vector and divergence
        if(nplanes .gt. 0) then
          numimgs = iget_images(imgprob, basisvects(1,1), 2)
          numimgs = iget_div_images(imgprob, divbasis)
        endif


        do imgnum=1, numimgs
          if(imgnum .eq. 1) then

            if (isourc .eq. ifield) then
              call get_uv(isourc,ipnt,2,us,vs)
              rs = get_rs(ifield, us, vs, utest, vtest)
              call get_basis(isourc,inode,utest,vtest,basistest)
              divbasistest = fget_divbasis(isourc,inode,utest,vtest)

              if( ak*min(rimgs(imgnum,ipnt,ifldpnt),rs) .gt. 
     &                              eps ) then
                cg1 = cg(imgnum,ipnt,ifldpnt)
                w = fget_weight(isourc,ipnt,2)
                do k = 1,3
                  cvpot(k) = cvpot(k) + basisvects(k,imgnum)*cg1 -
     &                       w*basistest(k)/rs
                enddo
                cspot = cspot + divbasis(imgnum)*cg1 - w*divbasistest/rs
              else
                do k = 1,3
                  cvpot(k) = cvpot(k) - basisvects(k,imgnum)*cj*ak
                enddo
                cspot = cspot - divbasis(imgnum)*cj*ak
              endif
            else
              cg1 = cg(imgnum,ipnt,ifldpnt)
              do k = 1,3
                cvpot(k) = cvpot(k) + basisvects(k,imgnum)*cg1
              enddo
              cspot = cspot + divbasis(imgnum)*cg1
            endif
          else
            cg1 = cg(imgnum,ipnt,ifldpnt)
            do k = 1,3
              cvpot(k) = cvpot(k) + basisvects(k,imgnum)*cg1
            enddo
            cspot = cspot + divbasis(imgnum)*cg1
          endif
        end do
      end do
      
      if(isourc.eq.ifield) then
        do k = 1,3
          cvpot(k) = cvpot(k) + basistest(k)*rsint/aa(isourc)
        enddo
        cspot = cspot + divbasistest*rsint/aa(isourc)
      endif

      return
      end





      integer function iget_images(imgprob, vects, itype)

      include 'ctparm.inc'
      include 'symm.inc'

      integer imgprob, itype, i, numimgs, itemp, j
      real vects(3,8), result(3)

      do i=2,8
        vects(1,i) = 0.
        vects(2,i) = 0.
        vects(3,i) = 0.
      end do
      numimgs = 1

      do i=1,nplanes
        itemp = numimgs
        do j=1,itemp
          call reflect_vector(isymplane(i), vects(1,j), result)

          if(itype .eq. 2 .and. iplanetype(imgprob,i) .eq. 1) then
c pec plane and not doing positions
            call vecscal(-1., result, vects(1,numimgs+1))
          else
            call vecscal( 1., result, vects(1,numimgs+1))
          endif

          numimgs = numimgs + 1
        end do
      end do

      iget_images = numimgs

      return
      end


      

            


      subroutine reflect_vector(iplane, vectinit, vectrefl)

      integer iplane
      real vectinit(3), vectrefl(3)
      integer i
      real normvec(3), tempvec(3), d
      real dotp

      do i=1,3
        normvec(i) = 0.
      end do

      if(iplane .eq. 1) then               ! xy plane
      	normvec(3) = 1.
      elseif(iplane .eq. 2) then           ! xz plane
        normvec(2) = 1.
      elseif(iplane .eq. 3) then           ! yz plane
        normvec(1) = 1.
      else
        write(6,*) '*** internal error ***'
        write(6,*) 'unknown symmetry plane in reflect_vector'
        stop
      endif

      d = dotp(normvec,vectinit)
      call vecscal(-2.*d, normvec, tempvec)
      call veclin(1., vectinit, 1., tempvec, vectrefl)

      return
      end





      subroutine get_r_mult(posfield, posvects, numimg, rimgs)

      integer numimg, i
      real posfield(3), posvects(3,8), rimgs(8)

      do i=1,numimg
        rimgs(i) = sqrt( (posfield(1)-posvects(1,i))**2 +
     &                   (posfield(2)-posvects(2,i))**2 +
     &                   (posfield(3)-posvects(3,i))**2 )
      end do

      return
      end
                         




      subroutine add_basis(imgprob)

      include 'ctparm.inc'
      include 'ctgeom.inc'
      include 'symm.inc'
      
      integer imgprob, numchanged, i
      integer ichanged(mxedgs)
      integer iplanenum, iedge


      if(numchanged .gt. 0) then
        do i=1,numchanged
          ibkp(ichanged(i)) = 0
        end do
        nunkns = nunkns - numchanged
        numchanged = 0
      endif

      do iplanenum=1,nplanes
        if(iplanetype(imgprob,iplanenum) .eq. 1) then  ! pec 
          do iedge=1,nedges
            if(ionpln(iedge) .eq. iplanenum) then
              nunkns = nunkns + 1
              numchanged = numchanged + 1
              ichanged(numchanged) = iedge
              ibkp(iedge) = nunkns
            endif
          end do
        endif
      end do

      return
      end



      integer function iget_div_images(imgprob, divbasis)

      include 'ctparm.inc'
      include 'ctgeom.inc'
      include 'symm.inc'

      integer imgprob, i, j, numimgs, itemp
      real divbasis(8)

      
      do i=2,8
        divbasis(i) = 0.
      end do
      numimgs = 1

      do i=1,nplanes
        itemp = numimgs
        do j=1,itemp
          if(iplanetype(imgprob,i) .eq. 1) then         !pec
            divbasis(numimgs+1) = -divbasis(j)
          elseif(iplanetype(imgprob,i) .eq. 2) then     !pmc
            divbasis(numimgs+1) = divbasis(j)
          else
            write(6,*) '*** internal error ***'
 	    write(6,*) 'unknown symmetry plane type',
     %                 ' in iget_div_image'
            stop
          endif

          numimgs = numimgs + 1
        end do
      end do

      iget_div_images = numimgs

      return
      end
