      PROGRAM PGDEM9
C-----------------------------------------------------------------------
C Test program for PGPLOT: test of imaging routine PGPIXL.
C-----------------------------------------------------------------------
      INTEGER PGOPEN
      INTEGER N, NCOL, NLEV
      PARAMETER (N=64, NCOL=32, NLEV=9)
      INTEGER I,J,CI1,CI2
      REAL F(N,N),FMIN,FMAX,R,G,B,CLEV(NLEV),TR(6)
      INTEGER IA(N,N)
C
C Compute a suitable function.
C
      DO 15 I=1,N
         DO 10 J=1,N
            F(I,J) = COS(0.6*SQRT(I*2.)-0.4*J/3.)*COS(0.4*I/3)+
     1           (I-J)/REAL(N)
 10      CONTINUE
 15   CONTINUE
C
C Find max and min
C
      FMIN = F(1,1)
      FMAX = F(1,1)
      DO 18 I=1,N
         DO 17 J=1,N
            FMIN = MIN(F(I,J),FMIN)
            FMAX = MAX(F(I,J),FMAX)
 17      CONTINUE
 18   CONTINUE
C
C Convert to color indices
C
      DO 25 I=1,N
         DO 24 J=1,N
            IA(I,J) = (F(I,J)-FMIN)/(FMAX-FMIN)*(NCOL-1)+16
 24      CONTINUE
 25   CONTINUE
C
C Open plot device and set up coordinate system. We will plot the
C image within a unit square.
C
      IF (PGOPEN('?') .LT. 1) STOP
      CALL PGQCOL(CI1, CI2)
      IF (CI2.LT. 15+NCOL) THEN
         WRITE (*,*) 'This program requires a device with at least',
     1        15+NCOL,' colors'
         STOP
      END IF
      CALL PGPAGE
      CALL PGSCR(0, 0.0, 0.3, 0.2)
      CALL PGSVP(0.05,0.95,0.05,0.95)
      CALL PGWNAD(0.0, 1.0, 0.0, 1.0)
C
C Set up a color palette using NCOL indices from 16 to 15+NCOL.
C
      DO 30 I=1,NCOL
         R = REAL(I-1)/REAL(NCOL-1)*0.8 + 0.2
         G = MAX(0.0, 2.0*REAL(I-1-NCOL/2)/REAL(NCOL-1))
         B = 0.2 + 0.4*REAL(NCOL-I)/REAL(NCOL)
         CALL PGSCR(I+15, R, G, B)
 30   CONTINUE
C
C Use PGPIXL to plot the image.
C
      CALL PGPIXL(IA,N,N, 1, N, 1, N, 0.0, 1.0, 0.0, 1.0)
C
C Annotation.
C
      CALL PGSCI(1)
      CALL PGMTXT('t',1.0,0.0,0.0,'Test of PGPIXL')
      CALL PGBOX('bcnts',0.0,0,'bcnts',0.0,0)
C
C Overlay a contour map.
C
      TR(1) = -1.0/REAL(N-1)
      TR(2) = 1.0/REAL(N-1)
      TR(3) = 0.0
      TR(4) = -1.0/REAL(N-1)
      TR(5) = 0.0
      TR(6) = 1.0/REAL(N-1)
      DO 40 I=1,NLEV
         CLEV(I) = FMIN + (FMAX-FMIN)*REAL(I)/REAL(NLEV)
 40   CONTINUE
      CALL PGCONT(F, N, N, 1, N, 1, N, CLEV, NLEV, TR)
C
C Done.
C
      CALL PGCLOS
      END
