

      subroutine sload(cz)
c******************************************************************
c  purpose: to fill the nunknsxnunkns square impedance matrix cz 
c           with the material load term
c
c notes: 1. the term <fm,fn> will not have any contributions from
c        the images of fn.
c        2. the inner product has been computed anayltically for
c        wires. the wire portion of a junction basis function can
c        be loaded.
c        3. the inner product is done using 4 integration points for
c	 patch basis functions. 
c        4. loads on the junction basis function patches are not
c        currently allowed.
c******************************************************************
      implicit complex (c)
      parameter (np=7)
      include 'ctparm.inc'
      include 'ctgeom.inc' 
      include 'ctgeom2.inc' 
      include 'matl.inc'
      include 'dspparam.inc'
      include 'wirechar.inc'
      include 'ctgeomwire.inc'

      dimension cz(mxedgs,mxedgs)
      real    ru(3), rv(3), u(np), v(np)
      real    dum(3), rjac(np)
      real    f_fld(3), f_src(3)
      logical wire_type, patch_type
      complex patch_load
      integer iwirenum

      i4pts = 1		!use 4 point integration for triangles

c.. add in effects of distributed loads on surfaces (and wires)
c.. loop over triangles & wires
      do 900 ijk = 1,nfaces
        patch_load = cmatlst(ixmat(ijk))
        if( cabs(patch_load) .eq. 0. ) go to 900	       !no load 

        numbasis = iget_num_basis(ijk) 
        if( ntype(ijk) .eq. ep_pattypwire ) then
           wire_type = .true. 
           patch_type = .false.
        else
           wire_type = .false. 
           patch_type = .true.
        endif

c.. for triangles, precompute the (u,v) points & the jacobians
        if(patch_type) then
          nint_tst = iget_num_int_points(ijk,i4pts)
          do i = 1,nint_tst
             call get_uv(ijk,i,i4pts,u(i),v(i))
             call tanvec_uandv(ijk,u(i),v(i),ru,rv)
             call crossp(ru,rv,dum)
             rjac(i) = sqrt(dotp(dum,dum))
          enddo
        endif

c--  loop over nodes of source patch
        do 460 isnod = 1,numbasis
          jl = iget_matrix_position(ijk,isnod)
          if(jl.eq.0) go to 460
          if( patch_type .and. isnod.gt.3) go to 460   !no junction patches

c--  loop over nodes of testing patch
          do 450 ifnod = 1,numbasis
            il = iget_matrix_position(ijk,ifnod)
            if(il.eq..0) go to 450
            if( patch_type .and. isnod.gt.3) go to 460  !no junction patches
c
            if(wire_type) then
              call basis_innerprod_on_wire(ijk,isnod,ifnod,term)
            else
              term = 0.
              do i = 1,nint_tst	
                wgt_tst = fget_weight(ijk,i,i4pts)
c---   form inner product of source & testing basis fns
                call get_basis(ijk,isnod,u(i),v(i),f_src)
                call get_basis(ijk,ifnod,u(i),v(i),f_fld)
c---   note: here, f_src & f_fld are missing 1/(jacobian*area) term.
c      also, one 1/jac*area gets cancelled in the integration.
                tmp1 = dotp(f_fld,f_src) / ( aa(ijk) * rjac(i) )
                term = term +  wgt_tst * tmp1
              enddo
            endif

            cz(il,jl) = cz(il,jl) + patch_load * term

  450     continue  ! field node

  460   continue    ! source node

  900 continue   ! patch loop


c     add in effects of lumped loads on wires.  note: lumped loads
c     only occur at wire feed points, therefore there can be only
c     one lumped load per wire.
      do iwirenum=1,numwires
        index = inumfeedref(iwirenum)
        if(index .gt. 0) then
          cz(index,index) = cz(index,index) + wlump(iwirenum)
        endif
      end do


      return
      end




      subroutine basis_innerprod_on_wire(iwire,isnod,ifnod,answer)
c----------------------------------------------------------------------
c purpose: to compute the inner product of two basis functions on a
c          wire (iwire), and the 2 basis fns are associated with
c          node indices isnod & ifnod. the result is returned in answer.
c notes: 1. i*nod = (1,2) for the wire basis fns & (3,4) for the 
c           junction basis function.
c----------------------------------------------------------------------

      include 'ctparm.inc'
      include 'ctgeom.inc'
      include 'ctgeomwire.inc'
      real rlen(3)

c the integral int(u^2) = int((1-u)^2) from [0,1] = 1/3
c the integral int(u*(1-u))            from [0,1] = 1/6
c depending on isnod & ifnod, get the correct integral value...
      ai1 = 1./3.
      ai2 = 1./6.
      if(ifnod.eq.isnod) then
        aint = ai1
      elseif( abs(isnod-ifnod).eq.2 ) then	! 1 & 3  or  2 & 4
        aint = ai1
      else
        aint = ai2
      endif
 
      iflg = 1
      if(isnod.eq.3 .xor. ifnod.eq.3) iflg = -1

c... get the length of the line segment
      call vecdif(datnod(1,nconn_wire(2,iwire)),
     %            datnod(1,nconn_wire(1,iwire)), rlen)
      alength = vecmag(rlen)

      answer = alength * aint * iflg

      return
      end
