Program CT
*------------------------------------------------------------------*
*                                                                  *
*  A First Course in Computational Physics,   EXERCISE 6.18        *
*                                                                  *
*  YOUR NAME and identifying information go here.                  *
*                                                                  *
*  The program CT performs image reconstruction from               *
*  given "x-ray" projections.                                      *
*                                                                  *
*                                                                  *
*                     last revised:  January 1, 1993    pld        *
*                 current revision:  today's date       abc        *
*                                                  (your initials) *
*                                                                  *
*------------------------------------------------------------------*
*
       double precision image(51,51), xi(32), projection(32)
       double precision pi, phi, k
       complex*16 temporary(32)
       integer log_2_N, N, halfN, num_phi, i, j, count, 
     +       index, pixels 
       parameter ( pi = 3.14159 26535 89793 d0, log_2_N = 5, 
     +            num_phi=12, pixels = 51 ) 
*
* Step 0: Initialize everything, zero the image
*
       N = 2**log_2_N
       halfN = N/2

       DO i = 1, pixels
          DO j = 1, pixels
             image(i,j) = 0.d0
          END DO
       END DO
       ...
*
*  If your computer is equipped with VGA, you might try 
*  using a grey scale instead of the default color scheme.
*
*  With the VGA you have, typically, 16 colors available to 
*  you to use.  These colors are referred to by their INDEX,
*  as we've previously seen.  But we can change the color 
*  associated with a particular index, using the graphics
*  command
*
*          call color_map( index, red, green, blue )
*
*  where the arguments are all integers.  This "maps" the color 
*  having particular RED, GREEN, and BLUE components to INDEX.
*  The components range in value from 0 to 63, so that there is a 
*  total of 64*64*64=262,144 possible colors from which to choose.
*
*  Greys are obtained by having equal components of red, green, and
*  blue.  The following redefines the palette of available colors
*  to a graduated scale of greys:
*
*       call color_map(  0,  0,  0,  0 ) 
*       call color_map(  1,  7,  7,  7 )
*       call color_map(  2, 11, 11, 11 )
*       call color_map(  3, 15, 15, 15 )
*       call color_map(  4, 19, 19, 19 )
*       call color_map(  5, 23, 23, 23 )
*       call color_map(  6, 27, 27, 27 )
*       call color_map(  7, 31, 31, 31 )
*       call color_map(  8, 35, 35, 35 )
*       call color_map(  9, 39, 39, 39 )
*       call color_map( 10, 43, 43, 43 )
*       call color_map( 11, 47, 47, 47 )
*       call color_map( 12, 51, 51, 51 )
*       call color_map( 13, 55, 55, 55 )
*       call color_map( 14, 59, 59, 59 )
*       call color_map( 15, 63, 63, 63 )
*
*
* Loop over Steps 1-5:
*
       DO  count = 1, num_phi
          phi =  ...
*
* STEP 1: Get the total projection at phi
*
          delta_xi = 0.1d0    
          DO  i = 1, N
             xi(i) = -1.6d0 + delta_xi * dble(i-1)
             projection(i) = ct_projection(xi(i),phi)
          END DO
*
* STEP 2: Fourier transform
*
          DO i = 1, halfN
             temporary(  i      ) = projection(i + halfN)
             temporary(i + halfN) = projection(  i      )
          END DO

          call fft( temporary, log_2_N, 0)
*
* STEP 3: Multiply by the absolute value of k
*         
          delta_k = 2.d0 * pi / (N * delta_x)
          DO i = 1, halfN
             k                =  dble(i-1) * delta_k
             temporary(i)     = abs(k) * temporary(i) 
             k                = -dble( i ) * delta_k
             temporary(N+1-i) = abs(k) * temporary(N+1-i)
          END DO

          call fft( temporary, log_2_N, 1)

          DO i = 1, halfN
             projection( i     ) = REAL( temporary(i+halfN) )
             projection(i+halfN) = REAL( temporary( i     ) )
          END DO

*
* STEP 4: Update IMAGE
*
          DO i = 1, pixels
             x = -1.d0 + dble(i-1)*2.d0 / dble(pixels-1)
             DO j = 1, pixels
                y = -1.d0 + dble(j-1)*2.d0 / dble(pixels-1)
                xi_value = x * cos(phi) + y * sin(phi)
* 
*       Find "index" such that xi(index) < xi < xi(index+1)
*
                index =  (xi_value + 1.6d0)/delta_xi + 1.d0
*
*       Interpolate
* 
                image(i,j) = image(i,j) 
     +       +  ((xi_value - xi(index))*projection(index+1) +
     +           (xi(index+1)-xi_value)*projection( index ) )
     +           /delta_xi
             END DO
          END DO
*
* STEP 5: Display IMAGE
*
          DO i = 1, pixels
             DO j = 1, pixels 
              < index = ???, depending upon image(i,j) >
                call color( index )
* 
* Remember, Subroutine PIXEL counts from top to bottom!
*
                call pixel( i+100, 200-j )
             END DO
          END DO
*
* STEP 6: Repeat steps 1-5 
*
       END DO
       end