      PROGRAM PGEX16
C-----------------------------------------------------------------------
C Test program for PGPLOT: test of contouring routines PGCONT, PGCONS,
C and PGCONB.
C-----------------------------------------------------------------------
      REAL BLANK
      PARAMETER (BLANK=-1.2E20)
      INTEGER I,J
      REAL F(40,40),FMIN,FMAX,ALEV,TR(6),X,Y,R
      DATA TR/0.,1.,0.,0.,0.,1./
C
C Open device for graphics.
C
      CALL PGBEGIN(0,'?',1,1)
C
C Compute a suitable function.
C
      FMIN = F(1,1)
      FMAX = F(1,1)
      DO 20 I=1,40
          DO 10 J=1,40
              F(I,J) = COS(0.3*SQRT(I*2.)-0.4*J/3.)*COS(0.4*I/3)+
     1                     (I-J)/40.0
              FMIN = MIN(F(I,J),FMIN)
              FMAX = MAX(F(I,J),FMAX)
   10     CONTINUE
   20 CONTINUE
C
C (1) Test PGCONT.
C Clear the screen. Set up window and viewport.
C
      CALL PGPAGE
      CALL PGVPORT(0.05,0.95,0.05,0.95)
      CALL PGWINDOW(1.0,40.0,1.0,40.0)
      CALL PGBOX('bcts',0.0,0,'bcts',0.0,0)
      CALL PGMTEXT('t',1.0,0.0,0.0,'PGEX16: Contouring using PGCONT')
C
C Draw the map.  PGCONT is called once for each contour, using
C different line attributes to distinguish contour levels.
C
      DO 30 I=1,21
          ALEV = FMIN + (I-1)*(FMAX-FMIN)/20.0
          IF (MOD(I,5).EQ.0) THEN
              CALL PGSLW(3)
          ELSE
              CALL PGSLW(1)
          END IF
          IF (I.LT.10) THEN
              CALL PGSCI(2)
              CALL PGSLS(2)
          ELSE
              CALL PGSCI(3)
              CALL PGSLS(1)
          END IF
          CALL PGCONT(F,40,40,1,40,1,40,ALEV,-1,TR)
   30 CONTINUE
      CALL PGSLW(1)
      CALL PGSLS(1)
      CALL PGSCI(1)
C
C (2) Test PGCONS.
C Clear the screen. Set up window and viewport.
C
      CALL PGPAGE
      CALL PGBOX('bcts',0.0,0,'bcts',0.0,0)
      CALL PGMTEXT('t',1.0,0.0,0.0,'PGEX16: Contouring using PGCONS')
C
C Draw the map.  PGCONS is called once for each contour, using
C different line attributes to distinguish contour levels.
C
      DO 40 I=1,21
          ALEV = FMIN + (I-1)*(FMAX-FMIN)/20.0
          IF (MOD(I,5).EQ.0) THEN
              CALL PGSLW(3)
          ELSE
              CALL PGSLW(1)
          END IF
          IF (I.LT.10) THEN
              CALL PGSCI(2)
              CALL PGSLS(2)
          ELSE
              CALL PGSCI(3)
              CALL PGSLS(1)
          END IF
          CALL PGCONS(F,40,40,1,40,1,40,ALEV,-1,TR)
   40 CONTINUE
      CALL PGSLW(1)
      CALL PGSLS(1)
      CALL PGSCI(1)
C
C (3) Test PGCONB. This is the same as the test of PGCONS, except
C first we blank the data outside an annulus.
C
      DO 60 I=1,40
          DO 50 J=1,40
              R = SQRT((I-20.5)**2 + (J-20.5)**2)
              IF (R.GT.20.0 .OR. R.LT.3.0) F(I,J) = BLANK
   50     CONTINUE
   60 CONTINUE
C
      CALL PGPAGE
      CALL PGBOX('bcts',0.0,0,'bcts',0.0,0)
      CALL PGMTEXT('t',1.0,0.0,0.0,'PGEX16: Contouring using PGCONB')
      DO 80 I=1,21
          ALEV = FMIN + (I-1)*(FMAX-FMIN)/20.0
          IF (MOD(I,5).EQ.0) THEN
              CALL PGSLW(3)
          ELSE
              CALL PGSLW(1)
          END IF
          IF (I.LT.10) THEN
              CALL PGSCI(2)
              CALL PGSLS(2)
          ELSE
              CALL PGSCI(3)
              CALL PGSLS(1)
          END IF
          CALL PGCONB(F,40,40,1,40,1,40,ALEV,-1,TR,BLANK)
   80 CONTINUE
C
C Mark the blanked points for easy identification.
C
      CALL PGSCI(1)
      DO 100 I=1,40
          DO 90 J=1,40
              IF (F(I,J).EQ.BLANK) THEN
                  X = TR(1) + REAL(I)*TR(2) + REAL(J)*TR(3)
                  Y = TR(4) + REAL(I)*TR(5) + REAL(J)*TR(6)
                  CALL PGPOINT(1, X, Y, -1)
              END IF
   90     CONTINUE
  100 CONTINUE
C
C Close the device and exit.
C
      CALL PGEND
      END
