#define ITERS 100     /* number of data sets        */
#define NCHAN 4      /* number of antenna channels */
#define NRANGE 10    /* number of range elements   */
#define NDOP 512     /* number of doppler elements */
#define MDOP 9       /* log2 of doppler elements   */
#define NTARGETS 4   /* number of targets          */
#define K0 1         /* Doppler window is K0..K1   */
#define K1 40

#define CHANPROCS NCHAN    /* processors for channel dimension */

#define R1 1
#define C1 1
#define R2 2
#define C2 5
#define RANGEPROCS R2*C2  /* processors for range dimension */

c 
c radar.fx - radar benchmark program (1D distribution)
c 
c Dave O'Hallaron, Carnegie Mellon University, July 1993
c
c Adapted from G. Shaw et al, "Multiprocessors for Radar Signal Processing",
c Lincoln Laboratories, MIT, Lexington, MA, Tech Report 961, 17 Nov 1992.
c
      program radar

c     major distributed arrays 
      real inpr(2,NDOP,NRANGE,NCHAN)  ! input sensor array (before c.t.)
      integer det(NDOP,NRANGE,NCHAN)  ! output detect array

c     precomputed FFT data structures 
      real w(NDOP)             ! hamming window 
      integer brt(NDOP)        ! bit reverse table
      real twiddles(2,NDOP/2)  ! coefficients

c     target ranges and frequencies
      integer range0(NTARGETS)    
      integer doppler0(NTARGETS)  

c     misc items
      integer det2 (K1-K0+1,NRANGE) 
      integer k, iter
      real sigma0, thresh, const2, const1
      real pi, a0
      real cos, real

c     timing variables 
      integer fxcellid
      integer fclocks, tclocks, rclocks, dclocks
      integer fusecs, tusecs, rusecs, dusecs, usecs
      integer fpercent, tpercent, rpercent, dpercent

      common /timing/ fclocks, tclocks, rclocks, dclocks

c     data layout statements
      template tr(NCHAN)
      template td(NRANGE)
      align inpr(l,k,j,i) with tr(i)
      align det(k,j,i) with td(j)   

c     tasking change
c     distribute tr(BLOCK(CHANPROCS))
      distribute tr(BLOCK(R1*C1))

      distribute td(BLOCK(R2*C2))

c     initialize constants
      pi = 3.141592653589793
      sigma0 = 1.13 
      thresh = 2.421
      const2 = thresh/(NRANGE*(K1-K0+1)-1+thresh)
      const1 = 2*pi/NDOP
      
c     target frequencies
      doppler0(1) = 6
      doppler0(2) = 16
      doppler0(3) = 26
      doppler0(4) = 36

c     target ranges
      range0(1) = 3
      range0(2) = 5
      range0(3) = 7
      range0(4) = 9

c     onetime setup for FFT
      call gen_bit_reverse_table(brt,NDOP)
      call gen_W_table(twiddles,NDOP,MDOP,NDOP/2)

c     precompute the hamming window
      a0 = 0.5*(NDOP-1)
      do k = 1,NDOP
         w(k) = 0.54 + 0.46*cos(const1*(k-1-a0))
      enddo

c     create the input data set (one channel per node)
      call dgen(inpr, range0, doppler0, sigma0)

c
c     The radar benchmark
c
      tclocks = 0 ! transpose time
      fclocks = 0 ! fft/mag time
      rclocks = 0 ! reduce time
      dclocks = 0 ! detect time
      
c      call timer_start()
      enter tasking      
      do iter = 1,ITERS
         call sample(inpr)
	 output(inpr)
	 processor(R1,C1)
	 origin(0,0)
	 
         call compute(inpr, det, brt, twiddles, w, const2, det2)
	 input(inpr)
	 processor(R2,C2)
	 origin (0,1),(2,1),(4,1),(6,1)

c         call dummy()
c	 processor(7,1)
c	 origin(1,0)

c         call dummy()
c	 processor(8,2)
c	 origin(0,6)
      enddo
      exit tasking

c
c     print results for the first channel
c

c     moved inside compute
c      det2 = det(K0:K1,:,1)
      if (fxcellid .eq. 0) then
         tusecs = (tclocks/20)/ITERS
         fusecs = (fclocks/20)/ITERS
         rusecs = (rclocks/20)/ITERS
         dusecs = (dclocks/20)/ITERS
         usecs = tusecs+fusecs+rusecs+dusecs

         tpercent = (real(tusecs)/real(usecs))*100.0
         fpercent = (real(fusecs)/real(usecs))*100.0
         rpercent = (real(rusecs)/real(usecs))*100.0
         dpercent = (real(dusecs)/real(usecs))*100.0

         print *, "1D Distribution"
         print *, NCHAN, " channels",  NRANGE, " range cells"
         print *, "tpose     :", tusecs, tpercent
         print *, "fft/mag   :", fusecs, fpercent
         print *, "reduce    :", rusecs, rpercent
         print *, "detect    :", dusecs, dpercent
         print *, "total(us) :", usecs
         print *, ""

         print *, det2(5,2), det2(6,2), det2(7,2), 
     $        det2(8,2), det2(9,2)
         print *, det2(5,3), det2(6,3), det2(7,3), 
     $        det2(8,3), det2(9,3)
         print *, det2(15,5), det2(16,5), det2(17,5), 
     $        det2(18,5), det2(19,5)
         print *, det2(25,7), det2(26,7), det2(27,7), 
     $        det2(28,7), det2(29,7)
         print *, det2(35,9), det2(36,9), det2(37,9), 
     $        det2(38,9), det2(39,9)
         print *, det2(35,10), det2(36,10), det2(37,10),
     $        det2(38,10), det2(39,10)
      endif

      end


c      subroutine dummy()
c      end      
c
c     sample - sample the input sensor data
c
      subroutine sample(inpr)
      real inpr(2,NDOP,NRANGE,NCHAN)

      template tr(NCHAN)
      align inpr(l,k,j,i) with tr(i)
      distribute tr(BLOCK(R1*C1))

      nocheck inpr
      return
      end

c
c     compute - process the input samples
c

c
c     for tasking, we eliminate the inpr array, since corner turn
c     is implicit

      subroutine compute(inpd, det, brt, twiddles, w, const2,det2)
ct      real inpr(2,NDOP,NRANGE,NCHAN)  ! sensor inputs (before corner turn)
      integer det(NDOP,NRANGE,NCHAN)  ! detect array

ct   added
      integer det2 (K1-K0+1,NRANGE)

      integer brt(NDOP)               ! bit reverse table
      real twiddles(2,NDOP/2)         ! coefficients
      real w(NDOP)                    ! hamming window 
      real const2                     ! precomputed threshhold constant
      
c     intermediate distributed arrays
      real inpd(2,NDOP,NRANGE,NCHAN)  ! sensor inputs (after corner turn)
      real mag(NDOP,NRANGE,NCHAN)     ! magnitude array

c     data layout statements
      template tr(NCHAN)
      template td(NRANGE)
ct      align inpr(l,k,j,i) with tr(i)
      align inpd(l,k,j,i) with td(j)
      align mag(k,j,i) with td(j)   
      align det(k,j,i) with td(j)   
ct      distribute tr(BLOCK(CHANPROCS))
      distribute td(BLOCK(R2*C2))
      nocheck inpd, det

c     intermediate replicated arrays
      real x(2,NDOP)              ! windowed FFT input 
      real y(2,NDOP)              ! scaled  FFT output 
      real localthresh(NCHAN)     ! threshholds 

c     misc
      integer i, j, k
      real sqrt, sum
      real tmpmag(K1-K0+1), magsum

      integer fclocks, tclocks, rclocks, dclocks
      common /timing/ fclocks, tclocks, rclocks, dclocks

c     corner turn (transpose) the input data 
c      call timer_start()
ct      idxperm(inpd,1,4,3,2) = idxperm(inpr,1,4,3,2)
c      tclocks = tclocks + timer_stop()

c     windowing and Doppler processing
      do i = 1, NCHAN
c         call timer_start()
         pdo j = 1,NRANGE
         pin inpd(:,:,j,:)
         pout mag(:,j,:)
            x(1,:) = inpd(1,:,j,i)*w
            x(2,:) = inpd(2,:,j,i)*w
            call fft(x,brt,twiddles,NDOP,MDOP,NDOP/2)
            y = x/NDOP
            tmpmag = y(1,K0:K1)*y(1,K0:K1) + 
     $           y(2,K0:K1)*y(2,K0:K1)
            mag(K0:K1,j,i) = sqrt(tmpmag)
            
         endpdo
c         fclocks = fclocks + timer_stop()
         
c     pdo-merge does not work with tasking
c
c     Doppler subselection and reduction
c         call timer_start()
c         pdo j = 1,NRANGE
c         pin mag(:,j,:)
c         pmvars magsum
c         pinit
c            magsum = 0.0
c         pbody
c            do k = K0,K1
c               magsum = magsum + mag(k,j,i)
c            enddo
c         pmerge
c            magsum = left(magsum) + right(magsum)
c         endpdo
c         rclocks = rclocks + timer_stop()

      magsum = sum(mag(K0:K1,:,i))
c      rclocks = rclocks + timer_stop()
         
c     constant false alarm rate (CFAR) detection 
c         call timer_start()
         localthresh(i) = const2*magsum
         pdo j = 1,NRANGE
         pin mag(:,j,:)
         pout det(:,j,:)
            do k = K0,K1
               if (mag(k,j,i) .gt. localthresh(i)) then
                  det(k,j,i) = 1
               else
                  det(k,j,i) = 0
               endif
            enddo
         endpdo
c         dclocks = dclocks + timer_stop()
      enddo 
      det2 = det(K0:K1,:,1)
c      call timer_start()
      return
      end

c 
c     dgen - generate the input data set 
c
      subroutine dgen(inpr, range0, doppler0, sigma0)

      real inpr(2,NDOP,NRANGE,NCHAN)
      integer range0(NTARGETS), doppler0(NTARGETS)
      real sigma0
      integer i, j, k, m, range, doppler
      real arg1, arg2, pi, gauss1, cos, sin
      
      template tr(NCHAN)
      align inpr(l,k,j,i) with tr(i)   
c      distribute tr(BLOCK(CHANPROCS))
      distribute tr(BLOCK(1))
      nocheck inpr

      pi = 3.141592653589793
      arg1 = 2.0*pi/NDOP

c     add noise to the jth range vector
      pdo i = 1,NCHAN
      pout inpr(:,:,:,i)
         do j = 1,NRANGE
            do k = 1,NDOP
               inpr(1,k,j,i) = gauss1(0.0,sigma0)
               inpr(2,k,j,i) = gauss1(0.0,sigma0)
            enddo
         enddo
      endpdo

c     add a signal to selected range vectors
      do m = 1, NTARGETS
         range = range0(m)
         doppler = doppler0(m)
         pdo i = 1, NCHAN
         pin inpr(:,:,:,i)
         pout inpr(:,:,:,i)
            do k = 1,NDOP
               arg2 = arg1*(k-1)*doppler
               inpr(1,k,range,i) = inpr(1,k,range,i) + cos(arg2)
               inpr(2,k,range,i) = inpr(2,k,range,i) + sin(arg2)
            enddo
         endpdo
      enddo

      end

c
c     gauss1 - crude approximation of normal distribution formed
c              from sum of 12 uniforms.
c
      real function gauss1(mean,sigma)
      real mean, sigma

      real x,random1,real
      integer i
      real x

      x = 0.0
      do i=1,12
         x = x + real(random1())
      enddo
      gauss1 = mean + sigma*(x - 6.0)
      end

c 
c random1 - uniformly distributed [0,1) real
c
      real function random1() 
      integer a, m, m1, b, a0, a1, b0, b1
      integer mod
      real real

      save a
      data a/11111111/

      m = 100000000
      m1 = 10000
      b = 31415821

      a1 = a/m1
      a0 = mod(a,m1)
      b1 = b/m1
      b0 = mod(b,m1)
      a0 = mod((a0*b1 + a1*b0),m1)*m1 + a0*b0
      a = mod((a0 + 1),m)
      random1 = real(a)/m
      end

c
c     gen_bit_reverse_table
c           Initialize bit reverse table (it actually differed by 1).
c       
c     Postcondition: br(i) = bit-reverse(i-1) + 1       
c
      subroutine gen_bit_reverse_table(brt,n)
      integer n

      integer brt(n)
      integer i, j, k, nv2

      nv2 = n/2
      j = 1
      brt(1) = j
      do i = 2, n
         k = nv2
 6       continue
         if (k .lt. j) then
            j = j - k
            k = k/2
            goto 6
         endif
         j = j + k
         brt(i) = j
      enddo
      return
      end

c
c     gen_W_table: generate powers of W.
c
c           Postcondition: W(i) = W**(i-1)            
c
      subroutine gen_W_table(W,n,m,ndv2)
      integer n,m,ndv2

      real W(2,ndv2)
      integer i
      real pi
      real wr,wi
      real ptr, pti
      real cos, sin
      
      pi = 3.141592653589793

      wr = cos(pi/ndv2)
      wi = -sin(pi/ndv2)
      W(1,1) = 1.0
      W(2,1) = 0.0
      ptr = 1.0
      pti = 0.0
      do i = 2,ndv2
         W(1,i) = ptr*wr - pti*wi
         W(2,i) = ptr*wi + pti*wr
         ptr = W(1,i)
         pti = W(2,i)
      enddo
      return
      end

      subroutine fft(a, brt, W, n, m, ndv2)
      integer n,m,ndv2

      real a(2,n)
      integer brt(n)
      real W(2,ndv2)

      integer i,j
      integer powerOfW
      integer sPowerOfW
      integer ijDiff
      integer stage
      integer stride
      integer first
      real pwr,pwi,tr,ti
      real ir, ii, jr, ji
      
      nocheck
c
c     bit reverse step
c
      do i = 1,n
         j = brt(i)
         if (i .lt. j) then
            tr = a(1,j)
            ti = a(2,j)
            a(1,j) = a(1,i)
            a(2,j) = a(2,i)
            a(1,i) = tr
            a(2,i) = ti
         endif
      enddo
c     
c     butterfly computations
c     
      ijDiff = 1
      stride = 2
      sPowerOfW = ndv2
      do stage = 1,m
c     Invariant: stride = 2 ** stage            
c     Invariant: ijDiff = 2 ** (stage - 1)
         
         first = 1
         do powerOfW = 1, ndv2, sPowerOfW
            pwr = W(1,powerOfW)
            pwi = W(2,powerOfW)
            
c     Invariant: pwr + sqrt(-1)*pwi = W**(powerOfW - 1)

            do i = first, n, stride
               j = i + ijDiff
               jr = a(1,j)
               ji = a(2,j)
               ir = a(1,i)
               ii = a(2,i)
               tr = jr*pwr - ji*pwi
               ti = jr*pwi + ji*pwr
               a(1,j) = ir - tr
               a(2,j) = ii - ti
               a(1,i) = ir + tr
               a(2,i) = ii + ti
            enddo

            first = first + 1
         enddo
            
         ijDiff = stride
         stride = stride * 2
         sPowerOfW = sPowerOfW / 2
      enddo
      return
      end