      subroutine expw(cv,imgprob)

      include 'ctparm.inc'
      include 'dspparam.inc'
      include 'ctgeom.inc'
      include 'ctgeom2.inc'
      include 'ctgeomwire.inc'
      include 'ctgeomjunc.inc'
      include 'symm.inc'
      include 'wirechar.inc'

      integer imgprob
      real ak, pi, theta, phi
      complex ctsta, ctstp
      complex cv(1)
      complex htheta,hphi,etheta,ephi,ex,ey,ez,edott
      complex ceinc(3,8)
      real khat(3,8), ereal(3,8), eimag(3,8)
      integer i, j, ij, jl, ir, iii, j2, j3, j4, nj1, nj2, nj3
      integer numimgs, node, iref
      real r, argmnt, alamda
      complex cimp, cetemp
      integer n(3), k
      dimension r(3)
      data n/3,1,2/
      integer iget_images, iget_matrix_position
      real dotp
      complex cdotp
      character*1 operator
      integer iwire, index

      common /kkk/ ak,pi
      common /param/ theta,phi
      common /field/ htheta,hphi,alamda
      common /tst/ ctsta(3,6,mxfaces), ctstp(6,mxfaces)
      common /ieoper/ operator

       do 10 i=1,nunkns
   10  cv(i)=0.


      if(operator .eq. 'a') then
c       calculating antenna power pattern.  excitation is given by
c       feed points on wires.
        do iwire=1,numwires
          ifeednode = inumfeedref(iwire)
          index = ibkp_wire(ifeednode)
          if(index .eq. 0) then
c           feed node doesn't lie on 'wire'.  check to see if
c           it is on a junction node.
            do i=1,numjunctions
              if(nodejunction(i) .eq. ifeednode) then
c               the feed node lies on this junction node
                index = ibkp_junc(i)
              endif
            end do
          endif

          if(index .gt. 0) then
            cv(index) = cmplx( wamp(iwire)*cos(wphase(iwire)),
     &                         wamp(iwire)*sin(wphase(iwire)) )
          endif
        end do 
      else
c       calculating scattering (rcs).  ignore feed point information
c       and assume plane wave excitation.

c       free space impedance
        cimp=120.*pi

c       calculated incident field in theta and phi
        etheta =  hphi*cimp
        ephi   = -htheta*cimp

c       compute the amplitude of the incident electric field in
c       cartesian coordinates 
        ceinc(1,1) = etheta*cos(theta)*cos(phi)-ephi*sin(phi)
        ceinc(2,1) = etheta*cos(theta)*sin(phi)+ephi*cos(phi)
        ceinc(3,1) =-etheta*sin(theta)

c       reflect real and imaginary parts of incident wave separately
c       and recombine to get reflected complex incident fields
        do k = 1,3
          ereal(k,1) = 1./(2**nplanes)*(real(ceinc(k,1)))
        enddo
        numimgs = iget_images(imgprob, ereal(1,1), 2)

        do k = 1,3
          eimag(k,1) = 1./(2**nplanes)*(aimag(ceinc(k,1)))
        enddo
        numimgs = iget_images(imgprob, eimag(1,1), 2)

        do i=1,numimgs
          do j=1,3
            ceinc(j,i) = cmplx( ereal(j,i), eimag(j,i) )
          end do
        end do

c       incident directions need to be reflected as position vectors
        khat(1,1) = sin(theta)*cos(phi)
        khat(2,1) = sin(theta)*sin(phi)
        khat(3,1) = cos(theta)
        numimgs = iget_images(imgprob, khat(1,1), 1)

c       compute excitation on wires and surfaces
        do i=1,numimgs
          do 499 ij=1,nfaces
            if( ntype(ij) .eq. ep_pattypwire) then
              call rvec(ij,.5,0.,r)
              do 460 ir=1,4
                jl = iget_matrix_position(ij,ir)
                if(jl.eq.0) go to 460
                edott = cdotp(ceinc(1,i),ctsta(1,ir,ij))
                argmnt = dotp(r, khat(1,i))
                cetemp=edott*cexp(cmplx(0.0,-ak*argmnt))
                cv(jl)=cv(jl)+cetemp
460	      continue
            else

              call tridata(ij,j2,j3,j4,nj1,nj2,nj3)
c             excitation on patches is computed via the moments.  need
c             to get the centroid of the patch
              call rvec_indexed(ij,0,1,r)
              do 450 ir=1,6
                jl = iget_matrix_position(ij,ir)
                if(jl.eq.0) go to 450
                if(ir .gt. 3) then
                  iref = ir
                else
                  iref = n(ir)
                endif

                edott = cdotp(ceinc(1,i),ctsta(1,iref,ij))
                argmnt = dotp(r, khat(1,i))
                cetemp=edott*cexp(cmplx(0.0,-ak*argmnt))
                cv(jl)=cv(jl)+cetemp

  450         continue
            endif      !wire or triangle
  499     continue     ! face loop
        end do         !image loop
      endif            !antenna or scattering

      return
      end





      complex function cdotp(cvec, cvec2)

      complex cvec(3), cvec2(3)
      complex ctemp
      integer i

      ctemp = cmplx(0.,0.)

      do i=1,3
        ctemp = ctemp + cvec(i) * cvec2(i)
      end do

      cdotp = ctemp

      return
      end
