program DRSRFT
c>> 1996-06-19 DRSRFT Krogh  Minor change for C conversion.
c>> 1994-10-19 DRSRFT Krogh  Changes to use M77CON
c>> 1994-08-09 DRSRFT WVS Remove '0' from format
c>> 1993-02-04 DRSRFT CLL
c>> 1989-05-07 DRSRFT FTK, CLL
c>> 1989-05-04 DRSRFT FTK, CLL
c     Demo driver for SRFT -- Multi-dimensional real Fourier transform
c     ------------------------------------------------------------------
c--S replaces "?": DR?RFT, ?RFT
c     ------------------------------------------------------------------
      integer J, J1, J2, K, L, M(2), MS, N, N2, N4, ND
      real             A(16, 16), ONE, PI, S(3), SIG, SIGD, TEMP, ZERO
      parameter (PI = 3.1415926535897932384E0)
      parameter (ONE = 1.E0)
      parameter (ZERO = 0.E0)
      data M / 4, 4 /, ND / 2 /
c     ------------------------------------------------------------------
c                         Start of code -- Construct A
      N = 2 ** M(1)
      N2 = N /2
      N4 = N2 / 2
      SIGD = PI / N2
      do 20 J1 = 1, N
         do 10 J2 = 1, N
            A(J1, J2) = ZERO
            if (abs(J1-N2-1) + abs(J2-N2-1) .le. N4) A(J1, J2) = ONE
   10    continue
   20 continue
c     ------------------------------------------------------------------
c                Compute Fourier transform and apply sigma factors
      MS = 0
      call SRFT (A, 'A', M, ND, MS, S)
      do 50 J1 = 1, N, 2
         A(J1, N2+1) = ZERO
         A(J1+1, N2+1) = ZERO
         do 40 J2 = 1, N2
            SIG = ONE
            if (J1 .EQ. 1) then
c                            No change in SIG due to J
               if (J2 .NE. 1) then
                  J = J2 - 1
                  K = 1
               else
                  A(2, 1) = ZERO
                  go to 40
               end if
            else
               J = J1 / 2
               K = 0
            end if
c                         Get nontrivial sigma factors * SIG
   30       continue
            if (J .LT. N4) then
               TEMP = S(J)
            else if (J .EQ. N4) then
               TEMP = ONE
            else
               TEMP = S(N2-J)
            end if
            SIG = SIG * TEMP / (SIGD * real(J))
            if (K .EQ. 0) then
               if (J2 .NE. 1) then
                  J = J2 - 1
                  K = 1
                  go to 30
               end if
            else
c                                       Apply sigma factors
               if (J1 .EQ. 1) then
                  A(1, N-J2+2) = ZERO
                  A(2, N-J2+2) = ZERO
               else
                  A(J1, N-J2+2) = SIG * A(J1, N-J2+2)
                  A(J1+1, N-J2+2) = SIG * A(J1+1, N-J2+2)
               end if
            end if
            A(J1, J2) = SIG * A(J1, J2)
            A(J1+1, J2) = SIG * A(J1+1, J2)
   40    continue
   50 continue
      call SRFT (A, 'S', M, ND, MS, S)
      print '(/'' Smoothed A'')'
      do 60 L = 1, 9
         print'(9f8.4)', (A(L,N), N = 1, 9)
   60 continue
      stop
      end