c*********************************************************************
      subroutine initcon(nfaces,alamda,c1,c2)
c
c purpose: to initialize constants used in the mom solution
c
c      implicit complex (c)

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

      real alamda, pi, ak
      complex c1, c2

      common /kkk/ ak,pi
      
c-- set up some constants
      pi = 2. * acos(0.)

      ak = 2.*pi/alamda
      c1 = ak * 30. * (0.0,1.0)
      c2 = 30. * (0.0,1.0)/ak

c-- for every element (triangle or wire) in the model, set up the
c   numerical integration values
      do k = 1,nfaces
        if(ntype(k).eq.ep_pattypwire) then
            call set_int_wire(k)
        else
            call set_int_tria(k)
        endif
      enddo

      return
      end


      subroutine set_int_tria(ipat)

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

      integer ipat, itype, num
      real xsi(mxintpt),eta(mxintpt),w(mxintpt)
      real a1, b1, a2, b2

c**************************************************************
c testing integration points
      itype = 1

c 4 point testing
      num = 4
      xsi(1) = 1./3.
      xsi(2) = .6
      xsi(3) = .2
      xsi(4) = .2
      eta(1) = 1./3.
      eta(2) = .2
      eta(3) = .6
      eta(4) = .2
      w(1) = -27./48.
      w(2) = 25./48.
      w(3) = w(2)
      w(4) = w(2)

      call set_num_int_points(ipat,itype,num,xsi,eta,w)

c******************************************************************
c    potential integration points
      itype = 2
      num = 7

      a1 = .0597158717
      b1 = .4701420641
      a2 = .7974269253
      b2 = .1012865073

      xsi(1) = 1./3.
      xsi(2) = a1
      xsi(3) = b1
      xsi(4) = b1
      xsi(5) = a2
      xsi(6) = b2
      xsi(7) = b2
      eta(1) = 1./3.
      eta(2) = b1
      eta(3) = a1
      eta(4) = b1
      eta(5) = b2
      eta(6) = a2
      eta(7) = b2
      w(1) = .2250000000
      w(2) = .1323941527
      w(3) = w(2)
      w(4) = w(2)
      w(5) = .1259391805
      w(6) = w(5)
      w(7) = w(6)

      call set_num_int_points(ipat,itype,num,xsi,eta,w)
c      call set_num_int_points(ipat,1,num,xsi,eta,w)  !7 point testing too

      return
      end





      subroutine set_int_wire(ipat)

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

      integer ipat, itype, num
      real xsi(mxintpt),eta(mxintpt),w(mxintpt)

c*******************************************************************
c 7 point potentials
c 
      itype = 2
      num = 7

      xsi(1) = -.9491079123
      xsi(2) = -.7415311855
      xsi(3) = -.4058451513
      xsi(4) =  0.00
      xsi(5) = -xsi(3)
      xsi(6) = -xsi(2)
      xsi(7) = -xsi(1)
      eta(1) = 0.00
      eta(2) = 0.00
      eta(3) = 0.00
      eta(4) = 0.00
      eta(5) = 0.00
      eta(6) = 0.00
      eta(7) = 0.00
      w(1) = .1294849661
      w(2) = .2797053914
      w(3) = .3818300505
      w(4) = .4179591836
      w(5) = w(3)
      w(6) = w(2)
      w(7) = w(1)

      call set_num_int_points(ipat,itype,num,xsi,eta,w)
c 7 point testing as well
      itype = 1
      call set_num_int_points(ipat,itype,num,xsi,eta,w)


      return
      end



      subroutine set_num_int_points(ipatch,itype,
     &                         num,eval_x,eval_e,w)

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

      integer ipatch, itype, num, i
      real eval_x(num),eval_e(num), w(num)

      numintpts(ipatch,itype) = num
      do i=1,num
        evalptsx(ipatch,itype,i) = eval_x(i)
        evalptse(ipatch,itype,i) = eval_e(i)
        weights(ipatch,itype,i) = w(i)
      end do

      return
      end
        



      subroutine define_junction

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

      integer ipatch, i, j, k, node, n
      real u_n, v_n, eta, xsi, u, v, v1mag, v2mag
      real rleft(3), rright(3), r_n(3)
      real vec1(3), vec2(3), argument, angle, totalangle
      real vecmag, dotp
      logical check_wire_nodes

      do ipatch=1,nfaces
        if(ntype(ipatch) .ne. ep_pattypwire) then
          do i=1,3
            u_n = tzpt(2*i,ipatch)
            v_n = tzpt(2*i-1,ipatch)
c            call rvec(ipatch,u_n,v_n,r_n)
            do k=1,3
              node = nbound(i+3,ipatch)
              r_n(k) = datnod(k,node)
            end do
            if(i .eq. 1) then
              eta = 0.
              xsi = .99
              call ac_trans(ipatch,xsi,eta,u,v)
              call rvec(ipatch,u,v,rleft)
              eta = .01
              call ac_trans(ipatch,xsi,eta,u,v)
              call rvec(ipatch,u,v,rright)
            elseif(i .eq. 2) then
              xsi = .01
              eta = .99
              call ac_trans(ipatch,xsi,eta,u,v)
              call rvec(ipatch,u,v,rleft)
              xsi = 0.
              call ac_trans(ipatch,xsi,eta,u,v)
              call rvec(ipatch,u,v,rright)
            elseif(i .eq. 3) then
              xsi = 0.
              eta = .01
              call ac_trans(ipatch,xsi,eta,u,v)
              call rvec(ipatch,u,v,rleft)
              eta = 0.
              xsi = .01
              call ac_trans(ipatch,xsi,eta,u,v)
              call rvec(ipatch,u,v,rright)
            endif
      
            call vecdif(r_n, rleft, vec1)
            call vecdif(r_n, rright, vec2)

            v1mag = vecmag(vec1)
            v2mag = vecmag(vec2)

            if(abs(argument) .gt. 1.e-7) then
 	      argument = dotp(vec1,vec2)/(v1mag*v2mag)
              if(abs(argument) .gt. 1.) 
     &                      argument = argument/abs(argument)
            endif
            angle =  abs( acos(argument) )

            alpha(i,ipatch) = angle
          end do
        else
c patch is a wire...make sure its nodes are not on a symm plane
          if( check_wire_nodes(ipatch,ibadnode) ) then
          write(6,*) '*** error ***'
          write(6,*) 'wires have been found in the data file with'
          write(6,*) 'one or more of their nodes on a symmetry plane'
          write(6,*)
          write(6,*) 'specifically:'
          write(6,*) '        wire segment number: ',ipatch
          write(6,*) '        node number:         ',ibadnode
          write(6,*) 'exiting...'
          stop
          endif
        endif
      end do


      do n=1,numjunctions
        node = nodejunction(n)
        totalangle = 0.
        do ipatch=1,nfaces
          if(ntype(ipatch) .eq. ep_pattypwire) then
            do i=1,2
              if(nconn_wire(i,ipatch) .eq. node) then
                ijuncmap(i,ipatch) = ibkp_junc(n)
              endif
            end do
          else
            do i=1,3
              if(nbound(i+3,ipatch) .eq. node) then
                ijuncmap(i,ipatch) = ibkp_junc(n)
                totpats(n) = totpats(n) + 1
                if( cabs(cmatlst(ixmat(ipatch))) .ne. 0. ) then
                  write(6,*) '*** error ***'
                  write(6,*) 'triangular patch found which has'
                  write(6,*) 'a load and also an attached wire.'
                  write(6,*) 
                  write(6,*) 'current versions of this code do not'
                  write(6,*) 'support the specification of materials'
                  write(6,*) 'around a junction point.'
                  write(6,*)
                  write(6,*) 'error caused by patch ',ipatch
                  write(6,*) 'which has a load impedance of',
     &                        cmatlst(ixmat(ipatch))
                  write(6,*) 'the node to which a wire is attached is'
                  write(6,*) 'number ',node
                  write(6,*) 'exiting...'
                  stop
                endif
              endif
            end do
          endif
        end do
      end do


      do n=1,numjunctions
        node = nodejunction(n)
        do ipatch=1,nfaces
          if(ntype(ipatch) .ne. ep_pattypwire) then
            do i=1,3
              if(nbound(i+3,ipatch) .eq. node) then
                alpha(i,ipatch) = alpha(i,ipatch)*totpats(n)
              endif
            end do
          endif
        end do
      end do


      return 
      end




      logical function check_wire_nodes(ipat,ibadnode)

      include 'ctparm.inc'
      include 'ctgeomwire.inc'

      integer ipat, ibadnode
      logical on_plane, result
      real eps
      logical check_node
      
      eps = 1.e-6
      result = .false.

      do i=1,2
        node_num = nconn_wire(i,ipat)
        if( check_node(node_num) ) then
          on_plane = .true.
          ibadnode = node_num
        endif
      end do

      check_wire_nodes = on_plane

      return
      end


      logical function check_node(node)

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

      integer node      
      logical on_plane
      real eps

      eps = 1.e-6
      on_plane = .false.

      do iplane = 1,nplanes
        if(isymplane(iplane) .eq. 1) then
c x-y plane... if z value less than eps node is on plane
          if( abs(datnod(3,node)) .lt. eps ) then
            on_plane = .true.
          endif

        elseif(isymplane(iplane) .eq. 2) then
c x-z plane... if y value less than eps node is on plane
 	  if( abs(datnod(2,node)) .lt. eps ) then
            on_plane = .true.
          endif

        elseif(isymplane(iplane) .eq. 3) then
c y-z plane... if x value less than eps node is on plane
          if( abs(datnod(1,node)) .lt. eps ) then
            on_plane = .true.
          endif
        else
c unknown plane type... stop
          write(6,*) '*** internal error***'
          write(6,*) ' unknown plane type',  isymplane(iplane)
          write(6,*) ' in check_node.'
          stop
        endif
      end do  ! planes

      check_node = on_plane

      return
      end
      
