       PROGRAM c3mmrcs
c*********************************************************************
c
c   purpose:  to compute the rcs of a target using the method of 
c             moments.
c             uses curved triangles with multiple parameterizations.
c
c  modified to compute rcs's for multiple polarizations 5/91
c  includes timing routines. 
c  single precision version.
c*******************************************************************
c
c    			   warning
c
c    this  code  contains  technical  data  whose  export  is
c    restricted  by the  arms export control act (22 usc 2751
c    et seq)  or the  export administration act  of  1979, as
c    amended  ( 50 usc 2401 et seq ).   violations  of  these
c    export  laws are  subject to  severe  civil and criminal
c    penalties.   distributions  limited  to   us  government
c    agencies and  their contractors. other requests for this 
c    code shall be referred to syracuse research corporation,
c    merrill lane, syracuse new york 13210-4080.
c
c*******************************************************************
c
c      implicit complex (c)

      include 'ctparm.inc'
      include 'dspparam.inc'
      include 'ctgeom.inc'	
      include 'ctgeom2.inc'	
      include 'mmprms2.inc'
      include 'mmloops.inc'
      include 'mmiounit.inc'
      include 'mmiounitd.inc'
      include 'ctgeomwire.inc'
      include 'ctgeomjunc.inc'
      include 'symm.inc'
      include 'symmd.inc'
      include 'wirechar.inc'

      real pi
      integer k_tim, nrcs, npoloop, imgprob, kt, kpt, ntpol
      integer kr, kpr, nrpol
      integer i, j, iunit
      real alamda, rcond, aztx, thetai, phii, azrx, thetar, phir
      real hmag, ak
      integer ips(mxedgs)
      real tpr(3,3)
      real*4 rcs(1000)
      complex htheta,hphi,ertheta,erphi
      complex cz(mxedgs,mxedgs)
      complex ci(mxedgs),cvtemp(mxedgs)
      complex cvsct(2,2,1000), chinc(2,2,1000), cs(1000)
      complex chtheta(1000),chphi(1000)
      complex c1, c2, ceth, ceph
      character*80 str
      logical zfill,yfill,wires,loads
      integer index, iwire
      complex input_voltage(mxfaces), input_current(mxfaces)
      complex input_impedance(mxfaces)
      real powerin
      character*1 operator
c      real timediff,timearray(2),dtime

      common /param/ thetai,phii
      common /kkk/ ak,pi
      common /field/ htheta,hphi,alamda
      common /geomopt/ wires,loads
      common /ieoper/ operator
 
      k_tim = 0		!for timing
      k_all = 0         !for timing
      nrcs = 0          ! final # of rcs values computed

      call initprm

      if( npols .eq. 1 ) then
        npoloop=1
      else
        npoloop=2		!go through tx & rx = v,h
      endif

      call ontimer(k_all)

      call rdcfil

      call initcon(nfaces,alamda,c1,c2)

      call optcheck(zfill,yfill,eta_alpha)


      call define_junction()

      call speedup()

      call gentesters

      do k = 1, 1000
         do j = 1,2
            do i = 1,2
               cvsct(i,j,k) = cmplx(0.,0.)
            end do
         end do
      end do

      do imgprob=1, 2**nplanes
        call add_basis(imgprob)

        write(6,*) 'image problem #',imgprob
        write(6,*) 'number of unknowns = ',nunkns

        do j = 1,nunkns		!always initialize [z] = 0
          do i = 1,nunkns
            cz(i,j) = cmplx(0.,0.)
          enddo
        enddo

        call ontimer(k_tim)

        if(zfill) then
           call ontimer(k_tim)
c           timediff = dtime(timearray)
           call zmatrx(c1,c2,cz,imgprob)
c           timediff = dtime(timearray)
c           print *,'Zmatrix fill(in sec.) =',timediff
           call shotimer(k_tim,str)
           write(6,*) ' time to fill [z] from zmatrx is:'
           write(6,'(a80)') str
        endif
        if(yfill) then
           call ontimer(k_tim)
           call ymatrx(eta_alpha,cz,imgprob)
           call shotimer(k_tim,str)
           write(6,*) ' time to fill [z] from ymatrx is:'
           write(6,'(a80)') str
        endif
        if(loads) then
           call ontimer(k_tim)
           call sload(cz)
           call shotimer(k_tim,str)
           write(6,*) ' time to fill [z] from sload is:'
           write(6,'(a80)') str
        endif

          call ontimer(k_tim)
        call cgeco(cz,mxedgs,nunkns,ips,rcond,cvtemp)
          call shotimer(k_tim,str)
          write(6,*) ' time to factor [z] from cgeco is:'
          write(6,'(a80)') str
          write(6,*) ' condition number (rcond) = ', rcond

        call tcomp(pitch,roll,tpr)
        aztx = azts
        nrcs = 0

          call ontimer(k_tim)

        if(operator .eq. 'a') then
c         antenna radiation does not need a loop over transmitter
          do j = 1,nunkns		!always initialize [v] = 0
            ci(j) = cmplx(0.,0.)
          enddo

          call expw(ci,imgprob)
          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
              input_voltage(iwire) = input_voltage(iwire) + ci(index)
            endif
          end do
            

          call cgesl(cz,mxedgs,nunkns,ips,ci,0)
          if(savej) call outci(ci)
          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
              input_current(iwire) = input_current(iwire) + ci(index)
            endif
          end do



c         receiver azimuth
          azrx = aztx
   
          do kr = 1,nstr           ! do over receivers for this x-mitter
            do kpr = 1,npoloop	 ! do over rx polarizations
              if(npols.eq.1) then
                nrpol = ntypol(2,1)
              else
                nrpol = kpr
              endif
c             get receiver polarization vector (ethetar,ephir)
              call parmconv(2,imode,azrx,tilt,nrpol,tpr,
     &                      thetar,phir,erphi,ertheta)

c             compute and store the incident h-field and scattered 
c             voltage (scattered field dotted with receiver polarization)
              call radint_nomom(ci,thetar,phir,imgprob,ceth,ceph)

c             update the number of "rcs's" computed
              if( kpr.eq.1 ) then
                nrcs = nrcs + 1
              endif
c             for antenna problems, the "scattered voltage" is stored 
c             in the first column.  note: the incident h field is stored
c             in the first row.
              chinc(1,1,nrcs) = htheta
              chinc(2,1,nrcs) = hphi
              cvsct(1,kpr,nrcs) = cvsct(1,kpr,nrcs) +
     &                                ceth*ertheta + ceph*erphi

            end do		       !rx pol loop
            azrx = azrx + azrdt      !increment receiver azimuth
          end do !receiver loop
        else
           kk = 1
          do 200 kt = 1,nstt
            do 195 kpt = 1,npoloop
              if( (kpt.eq.2) .and. (nstt.eq.1) ) nrcs=0
              if(npols.eq.1) then
                ntpol = ntypol(1,1)
              else
                ntpol = kpt
              endif
              call parmconv(1,imode,aztx,tilt,ntpol,tpr,
     &                      thetai,phii,hphi,htheta)

              do j = 1,nunkns		!always initialize [v] = 0
                ci(j) = cmplx(0.,0.)
              enddo

              if(zfill) call expw(ci,imgprob)
              if(yfill) call hinc(eta_alpha,ci,imgprob)

              goto 222
              print *,'nstt =',nstt,' npoloop= ',npoloop
              write(6,*) (j,kk,real(ci(j)),aimag(ci(j)), 
     &             j = 1,nunkns)
 1333         format(2(1x,'[',I3,',',I3,']','(',F9.4,',',F9.4,')'))

              kk = kk + 1
              print *,'-----------------'

 222          continue

              call cgesl(cz,mxedgs,nunkns,ips,ci,0)
c              if(savej) call outci(ci)

c-- determine the receiver azimuth depending on the mode
              if (imode .eq. 0) then          ! monostatic
                azrx = aztx
              elseif (imode .eq. 1) then      ! fixed bistatic angle
                azrx = aztx + biang
              elseif (imode .eq. 2) then      ! fixed transmitter
                azrx = azrs
              else                            ! fixed equiv. line of sight
                azrx = aztx + biang
              endif

              do 100 kr = 1,nstr           ! do over receivers for this x-mitter
                do 95 kpr = 1,npoloop	 ! do over rx polarizations
                  if(npols.eq.1) then
                    nrpol = ntypol(2,1)
                  else
                    nrpol = kpr
                  endif
                  call parmconv(2,imode,azrx,tilt,nrpol,tpr,
     &                          thetar,phir,erphi,ertheta)

c-- compute and store the incident h-field and scattered voltage
                  call radint(ci,thetar,phir,imgprob,ceth,ceph)

c-- update the number of "rcs's" computed based on whether nstt=1 or not
                  if( nstt.eq.1 .and. kpr.eq.1 ) then
                    nrcs = nrcs + 1
                  elseif ( nstt.ge.1 .and. kpt.eq.1 .and. 
     &                        kpr.eq.1 ) then
                    nrcs = nrcs + 1
                  endif
                  chinc(1,kpt,nrcs) = htheta
                  chinc(2,kpt,nrcs) = hphi
                  cvsct(kpt,kpr,nrcs) = cvsct(kpt,kpr,nrcs) +
     &                                  ceth*ertheta + ceph*erphi

  95            continue		       !rx pol loop
                azrx = azrx + azrdt      !increment receiver azimuth
 100          continue
 195        continue	  	       !tx pol loop

            if (imode.eq.3) biang = biang+angdt       !fixed equiv. los
            aztx = aztx + aztdt
 200      continue
          endif		!antenna or scatterer

          call shotimer(k_tim,str)
          write(6,*) ' time to solve for all [i] & e-scat is:'
          write(6,'(a80)') str

      end do                   !image problems

c         do k = 1, nrcs
c            do j = 1,2
c               do i = 1,2
c                  print *,'[',i,',',j,',',k,']',chinc(i,j,k)
c               end do
c            end do
c         end do

      if(operator .eq. 'a') then
c       for antenna problems we need to compute the input impedance
c       and total power input.
        powerin = 0.
        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
            powerin = powerin + abs( real( input_voltage(iwire)*
     &                conjg(input_current(iwire)) ) )
            input_impedance(iwire) = input_voltage(iwire) / 
     *                               input_current(iwire)
          endif
        end do
      endif


c     compute & output the rcs or power gain for each tx/rx pol 
      if(npols.eq.1) then
c       v and h-inc were stored in 'first' terms of cvsct & chinc
        do i = 1,nrcs
          hmag = sqrt(cabs(chinc(1,1,i))**2 + cabs(chinc(2,1,i))**2)
          if(operator .ne. 'a') then
            rcs(i) = pi * (cabs(cvsct(1,1,i)) /
     &                   ((alamda**2)*hmag))**2         !db per lambda sq.
          else
c           the radiated power for an antenna problem is given by
c                         1/(eta*pin) * |escat * er|**2 
c           where eta is the free space impedance = 377
c                 escat is the scattered e field
c                 er is the unit normal in the direction of the
c                      receiver polarization
c                 pin is the total input power
            rcs(i) = 377.*ak*ak/(4.*pi) * cabs(cvsct(1,1,i))**2/
     &                                       powerin
          endif

          if (rcs(i) .le. 0.) then
            rcs(i) = - 400.
          else
           rcs(i) = 10. * log10(rcs(i))
          endif
        enddo
        iunit = 91
        write(iunit,1000) (rcs(i),i=1,nrcs)       
      else
        do j = 1,npols
           print *,'This branch $$$$$$'
c--       get v & h-inc based on stored scattering matrix
          call getv(cvsct,ntypol(1,j),ntypol(2,j),nrcs,cs)
          call gethinc(chinc,ntypol(1,j),nrcs,chtheta,chphi)
          do i = 1,nrcs
            if(operator .ne. 'a') then
              hmag = sqrt(cabs(chtheta(i))**2 + cabs(chphi(i))**2)
              rcs(i) = pi * (cabs(cs(i)) /
     &                   ((alamda**2)*hmag))**2         !db per lambda sq.
            else
c           the radiated power for an antenna problem is given by
c                         1/(eta*pin) * |escat * er|**2 
c           where eta is the free space impedance = 377
c                 escat is the scattered e field
c                 er is the unit normal in the direction of the
c                      receiver polarization
c                 pin is the total input power
              rcs(i) = 4.*pi/377. * cabs(cs(i)) * cabs(cs(i)) /
     &                              powerin
            endif

            if (rcs(i) .le. 0.) then
              rcs(i) = - 400.
            else
              rcs(i) = 10. * log10(rcs(i))
            endif
          enddo
          iunit = 91 + j - 1
          write(iunit,1000) (rcs(i),i=1,nrcs)
        enddo
      endif
 1000 format(10f7.2)
c
 999  continue
      call shotimer(k_all,str)
      write(6,*) ' total time is:'
      write(6,*) str


c     for antenna problem we need to write out the computed input impedances
      if(operator .eq. 'a') then
        write(6,*) 
        write(6,*)
        write(6,'(20x,a30)') 'antenna input impedances:'
        write(6,*)
        write(6,'(12x,a10,5x,a10,4x,a10)')'antenna','real','imag'
        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
            write(6,'(16x,i5,9x,f7.2,7x,f7.2)') iwire,
     &                       real(input_impedance(iwire)),
     &                       aimag(input_impedance(iwire))
          endif
        end do
      endif


c      call exit
      stop 'program ends !'
      end


