C EXAMPLE 11.8:  Plot image and transform of 8 by 8 unit source
C in 64 by 64 (otherwise zero) array.
C
      PARAMETER (N=64) 
      DOUBLE PRECISION  A(N,N),A2(N,N),ERR,ERMAX,DATA
      COMPLEX*16  IMAGE(N,N),IMAGE2(N,N),W(3*N+8)
      LOGICAL   FORWD
C
C  NOTE:  DOUBLE PRECISION COMPLEX IS NOT FORTRAN STANDARD.
C         MANY COMPILER VENDORS MAY FAIL TO IMPLEMENT IT OR MAY
C         IMPLEMENT IT IN A DIFFERENT MANNER.  IF YOUR COMPILER
C         DOES NOT ACCEPT THE COMPLEX*16 STATEMENT, WE 
C         RECOMMEND THAT YOU CHECK WITH YOUR VENDOR'S
C         INSTRUCTIONS AND MODIFY THE ABOVE ACCORDINGLY.
C
      WRITE (*,*) 'COMPUTING...'
      LDA = N
C Set up data. IMAGE is original. IMAGE2 is scaled by (-1)**(I+J)
C to place Fourier coefficients in the correct place for viewing.
      DO 1 I = 1,N
        DO 1 J = 1,N
          A(I,J) = 0.0D0
          IF (I.GE.(N/2-4) .AND. I.LE.(N/2+4) .AND. J.GE.(N/2-4)
     *         .AND. J.LE.(N/2+4)) A(I,J) = 1.0D0
          IMAGE(I,J) = DCMPLX(A(I,J),0.0D0)
          IMAGE2(I,J) = IMAGE(I,J)*(-1)**(I-1+J-1)
    1 CONTINUE
C
C Plot image  with axonometric plotting routine.
C
CCCCCC      CALL AXNPLT (A,LDA,N)
C
C Compute Fourier transform of IMAGE and IMAGE2
C
      FORWD = .TRUE.
      CALL DCFT2D (N,IMAGE,LDA,W,FORWD)
      CALL DCFT2D (N,IMAGE2,LDA,W,FORWD)
C
C Compute magnitude of components of transforms.
C Actual transforms are unscaled and need to be divided by N*N
C to be correct.
C
      DO 2 I = 1,N
        DO 2 J = 1,N
           A(I,J) = ABS(IMAGE(I,J))
           A2(I,J) = ABS(IMAGE2(I,J))
    2 CONTINUE
C
C PLOT MODULOUS OF TRANSFORM
C
CCCCC      CALL AXNPLT (A,LDA,N)
CCCCC      CALL AXNPLT (A2,LDA,N)
C
C Compute inverse transform of image and see if it agrees with original.
C
      FORWD = .FALSE.
      CALL DCFT2D (N,IMAGE,LDA,W,FORWD)
      ERMAX = 0.0D0
      DO 3 I = 1,N
        DO 3 J = 1,N
           DATA = 0.0D0
           IF (I.GE.(N/2-4) .AND. I.LE.(N/2+4) .AND. J.GE.(N/2-4)
     *         .AND. J.LE.(N/2+4)) DATA = 1.0D0
           ERR = ABS( DATA-ABS(IMAGE(I,J))/(N*N) )
           IF (ERR .GT. ERMAX) ERMAX = ERR
    3 CONTINUE
      WRITE (*,*) 'DCFT2D RESULTS (EX 11.8: PLOTS HAVE BEEN SKIPPED)'
      WRITE (*,'(2X,A,1X,D18.10)') 'MAXIMUM ERROR IS ',ERMAX
C
      WRITE (*,*)
      WRITE (*,*) 'REFERENCE RESULTS FROM IBM PC/AT'
      WRITE (*,*) ' MAXIMUM ERROR IS    0.1885421720E-15'
C
      END
C
C      SUBROUTINE AXNPLT(A,LDA,N)
CC
CC       AXONOMETRIC PLOT OF 2 D ARRAY A
CC
CC       N: SIZE OF A, N.LE.350
CC       LDA: DIMENSION OF FIRST COMPONENT OF A IN CALLING PROGRAM
CC
CC       ASSUMPTIONS:    0.0 .LE. A(I,J) .LE 1.0
CC
C      INTEGER LDA,J,K,ISAMAX,N,IRET,I
C      DBLE A(LDA,*),WX(350),WY(350),X(350),Y(350),YMX(350),X0,DEL,YMAX
C      DBLE YM,YX
CC
CC       FIND THE LARGEST Y SO CAN SCALE
C      DO 100 J=1,N
C        K=ISAMAX(N,A(1,J),1)
C        YMX(J)=A(K,J)
C  100 CONTINUE
C      K=ISAMAX(N,YMX,1)
C      YMAX=YMX(K)
C      X0=0.1
C      DEL=.4
C      CALL AGRAF0(IRET,0)
C      CALL MINMAX(IRET,0.0, 1.0, 0.0, 1.0, YM, YX)
C      CALL HOWPLT(IRET,0,1,15)
C      DO 2 I=1,N
C        DO 1 J=1,N
C           X(J)=(J-1.)/(N-1.)*DEL+X0+(I-1.)/(N-1.)*DEL
C           Y(J)=(A(I,J)/YMAX)*DEL+X0+(I-1.)/(N-1.)*DEL
C    1   CONTINUE
C        IF(I.EQ.1)THEN
C           CALL NOWAIT
C           CALL PLOT1(IRET,N,X,Y,WX,WY)
C        ELSE
C           IF(I.EQ.N)CALL WAIT
C           CALL HOWPLT(IRET,0,1,15)
C           CALL PLOT(IRET,N,X,Y,WX,WY)
C        ENDIF
C    2 CONTINUE
C      RETURN
C      END