      subroutine radint(ci,theta,phi,imgprob,ceth,cephi)


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

      integer imgprob
      real theta, phi
      complex ci(1), htheta, hphi
      complex ctsta, ctstp
      real r(3,8)
      real khat(3), escat(3)
      complex cescat(3), carg, ctemp(8), ceth, cephi
      real rmoment(3,8), ak, pi
      integer ijk, numimgs, numimgs2, iii, k, i1, iref, icmp
      real argmnt, alamda
      integer n(3), numbas
      data n/3,1,2/
      integer iget_images, iget_matrix_position
      real dotp

      common/kkk/ak,pi
      common/field/htheta,hphi,alamda
      common/tst/ ctsta(3,6,mxfaces),ctstp(6,mxfaces)

      cescat(1) = 0.
      cescat(2) = 0.
      cescat(3) = 0.

      khat(1) = sin(theta)*cos(phi)
      khat(2) = sin(theta)*sin(phi)   
      khat(3) = cos(theta)

      do 499 ijk=1,nfaces
        call rvec_indexed( ijk, 0, 1, r(1,1))

        numimgs = iget_images(imgprob, r(1,1), 1)

        do iii=1,numimgs
          argmnt = dotp(r(1,iii), khat)
          carg = cmplx(0.0,ak*argmnt)
          ctemp(iii) = cexp(carg)
        end do

        numbas = iget_num_basis(ijk)
        do 460 k=1,numbas
          i1 = iget_matrix_position(ijk,k)
          if (i1 .eq. 0) go to 460

          if(ntype(ijk) .ne. 4) then
            if(k .gt. 3) then
              iref = k
            else
              iref = n(k)
            endif
          else
            iref = k
          endif


          do icmp=1,3
            rmoment(icmp,1) = real(ctsta(icmp,iref,ijk))
          end do

          numimgs2 = iget_images(imgprob,rmoment(1,1),2)

          do iii=1,numimgs2
            do icmp = 1,3
              cescat(icmp) = cescat(icmp) + 
     &                       ctemp(iii)*ci(i1)*rmoment(icmp,iii)
            end do
          end do

 460    continue
 499  continue


      ceth  =  cescat(1)*cos(theta)*cos(phi) + 
     &         cescat(2)*cos(theta)*sin(phi) - 
     &         cescat(3)*sin(theta)

      cephi = -cescat(1)*sin(phi) + cescat(2)*cos(phi)


      return
      end

