Subject: Contouring layers.

From: [email-masked] (Dave Kennison)
Date: Wed, 3 Jun 92 09:56:58 MDT

I am attaching to this message two demo/test programs with the latest in
layered contours.  The first demonstrates how to do them using AREAS and
the second how to do them using cell arrays.  I have changed the common
block in CPMPXY so that you can put the contour plot in any parallelogram
in 3-space and then project that parallelogram from an "eye" position 
into a projection plane perpendicular to the line of sight to a point
you're looking at.  

The AREAS version seems to have problems with colors "leaking".  Once 
the graphics conference is over, I will have to track down why this is
so.  I suspect that it has to do with inaccuracy in the 3-D to 2-D 
projection equations.  I will probably have to use double precision in
that part of the code, the way I now do in MAPTRN.  

THE AREAS VERSION:

      PROGRAM LAYERS
C
C Declare a common block for communication with the routines CPMPXY
C and TDINVS.
C
	COMMON /TDCOM1/ XA11,YA11,ZA11,XDM1,YDM1,ZDM1,XD1N,YD1N,ZD1N
C
C Declare required data arrays and workspace arrays.
C
	DIMENSION ZDAT(23,14,5),RWRK(1000),IWRK(1000),IAMA(12000)
        DIMENSION IASF(13)
        DIMENSION XCRA(1000),YCRA(1000)
        DIMENSION IAIA(10),IGIA(10)
C
C Declare arrays to hold the list of indices and the list of labels
C required by the label-bar routine.
C
        DIMENSION LIND(14)
        CHARACTER*10 LLBS(15)
C
C Declare the routine which will color the areas.
C
        EXTERNAL COLRAM
C
C Define an array for GKS aspect source flags.
C
        DATA IASF / 13*1 /
C
C Define the list of indices required by the label-bar routine.
C
        DATA LIND / 2,3,4,5,6,7,8,9,10,11,12,13,14,15 /
C
C Open GKS.
C
        CALL OPNGKS
C
C Turn off the clipping indicator.
C
        CALL GSCLIP (0)
C
C Set all aspect source flags to "individual".
C
        CALL GSASF (IASF)
C
C Force solid fill.
C
        CALL GSFAIS (1)
C
C Define color indices.
C
        CALL DFCLRS
C
C Do an appropriate call to SET.
C
	CALL SET    (.10,.90,.15,.95,-.9,.9,-.9,.9,1)
C
C Initialize the 3-D to 2-D projection package.  The first three
C arguments specify the eye position, and the last three specify
C the position of the point looked at.  The line from the eye to
C the point looked at is the "line of sight".  The plane passing
C through the point looked at and perpendicular to the line of
C sight is the projection plane.
C
	CALL TDINIT (4.,3.,1.75,.5,.5,.5)
C
C Tell CONPACK we've done the call to SET.
C
	CALL CPSETI ('SET - DO-SET-CALL FLAG',0)
C
C Tell CONPACK exactly what contour levels to use.  We want 13 levels,
C splitting the range into 14 equal bands, one for each of the 14
C colors available.
C
	CALL CPSETI ('CLS - CONTOUR LEVEL SELECTOR',1)
	CALL CPSETR ('CIS - CONTOUR INTERVAL SELECTOR',1.)
	CALL CPSETR ('CMN - CONTOUR MINIMUM',1.)
	CALL CPSETR ('CMX - CONTOUR MAXIMUM',13.)
C
C Tell CONPACK to generate X and Y coordinates in the range from 0 to 1.
C
	CALL CPSETR ('XC1 - X COORDINATE AT I = 1',0.)
	CALL CPSETR ('XCM - X COORDINATE AT I = M',1.)
	CALL CPSETR ('YC1 - Y COORDINATE AT J = 1',0.)
	CALL CPSETR ('YCN - Y COORDINATE AT J = N',1.)
C
C Tell CONPACK to use mapping 3, which, in the version of CPMPXY
C supplied, is the mapping from a particular layer of the unit cube
C to the plotter frame.
C
	CALL CPSETI ('MAP - MAPPING FUNCTION',3)
C
C Tell CONPACK to draw the boundary of the data grid.
C
	CALL CPSETI ('PAI - PARAMETER ARRAY INDEX',-1)
	CALL CPSETI ('CLU - CONTOUR LEVEL USE FLAG',1)
C
C Tell CONPACK not to generate any labels on contour lines.
C
	CALL CPSETI ('LLP - LINE LABEL POSITIONING',0)
C
C Spin up the random number generator.
C
	DO 100 I=1,123
	  TEMP=FRAN()
  100   CONTINUE
C
C Generate data at the bottom and top levels; we'll interpolate at
C the intermediate levels.
C
	CALL GENDAT (ZDAT(1,1,1),23,23,14,20,20,0.,14.)
	CALL GENDAT (ZDAT(1,1,5),23,23,14,20,20,0.,14.)
C
C Interpolate data at the intermediate levels.
C
	DO 103 ILEV=2,4
	  VLEV=.25*REAL(ILEV-1)
	  DO 102 I=1,23
	    DO 101 J=1,14
	      ZDAT(I,J,ILEV)=(1.-VLEV)*ZDAT(I,J,1)+VLEV*ZDAT(I,J,5)
  101       CONTINUE
  102     CONTINUE
  103   CONTINUE
C
C Define the corners of the parallelogram containing the contours.
C
	XA11=0.
	YA11=0.
C
	XDM1=1.
	YDM1=0.
	ZDM1=0.
C
	XD1N=0.
	YD1N=1.
	ZD1N=0.
C
C Do five plots, each at a different level in the unit cube.
C
	DO 104 ILEV=1,5
C
C Set the vertical level of the contour plot.
C
	  ZA11=.25*REAL(ILEV-1)
C
C Initialize the drawing of the contour plot.
C
	  CALL CPRECT (ZDAT(1,1,ILEV),23,23,14,RWRK,1000,IWRK,1000)
C
C Initialize the area map and put the contour lines into it.
C
	  CALL ARINAM (IAMA,12000)
	  CALL CPCLAM (ZDAT(1,1,ILEV),RWRK,IWRK,IAMA)
C
C Color the map.
C
	  CALL ARSCAM (IAMA,XCRA,YCRA,1000,IAIA,IGIA,10,COLRAM)
C
C Put black contour lines over the colored map.
C
	  CALL SFLUSH
	  CALL GSPLCI (0)
	  CALL CPCLDR (ZDAT(1,1,ILEV),RWRK,IWRK)
	  CALL SFLUSH
	  CALL GSPLCI (1)
C
  104   CONTINUE
C
C Draw a label bar for the plot, relating colors to values.
C
        CALL CPGETR ('ZMN',ZMIN)
        CALL CPGETR ('ZMX',ZMAX)
C
	DO 105 I=1,15
          CALL CPSETR ('ZDV - Z DATA VALUE',
     +                 ZMIN+REAL(I-1)*(ZMAX-ZMIN)/14.)
          CALL CPGETC ('ZDV - Z DATA VALUE',LLBS(I))
  105   CONTINUE
C
        CALL LBSETI ('CBL - COLOR OF BOX LINES',0)
	CALL LBLBAR (0,.05,.95,.05,.15,14,1.,.5,LIND,0,LLBS,15,1)
C
C Put a label at the top of the plot.
C
	CALL SFLUSH
	CALL GSLWSC (2.)
	CALL PLCHHQ (CFUX(.5),CFUY(.95),'LAYERED CONTOURS',.025,0.,0.)
C
C Advance the frame.
C
        CALL FRAME
C
C Close GKS.
C
        CALL CLSGKS
C
C Done.
C
        STOP
C
      END

      SUBROUTINE COLRAM (XCRA,YCRA,NCRA,IAIA,IGIA,NAIA)
C
        DIMENSION XCRA(*),YCRA(*),IAIA(*),IGIA(*)
C
C The arrays XCRA and YCRA, for indices 1 to NCRA, contain the X and Y
C coordinates of points defining a polygon.  The area identifiers in
C the array IAIA, each with an associated group identifier in the array
C IGIA, tell us whether the polygon is to be color-filled or not.
C
C
C Assume the polygon will be filled until we find otherwise.
C
        IFLL=1
C
C If any of the area identifiers is negative, don't fill the polygon.
C
        DO 101 I=1,NAIA
          IF (IAIA(I).LT.0) IFLL=0
  101   CONTINUE
C
C Otherwise, fill the polygon in the color implied by its area
C identifier relative to edge group 3 (the contour-line group).
C
        IF (IFLL.NE.0) THEN
          IFLL=0
          DO 102 I=1,NAIA
            IF (IGIA(I).EQ.3) IFLL=IAIA(I)
  102     CONTINUE
          IF (IFLL.GT.0.AND.IFLL.LT.15) THEN
            CALL GSFACI (IFLL+1)
            CALL GFA (NCRA-1,XCRA,YCRA)
          END IF
        END IF
C
C Done.
C
        RETURN
C
      END

      SUBROUTINE GENDAT (DATA,IDIM,M,N,MLOW,MHGH,DLOW,DHGH)
C
C This is a routine to generate test data for two-dimensional graphics
C routines.  Given an array "DATA", dimensioned "IDIM x 1", it fills
C the sub-array ((DATA(I,J),I=1,M),J=1,N) with a two-dimensional field
C of data having approximately "MLOW" lows and "MHGH" highs, a minimum
C value of exactly "DLOW" and a maximum value of exactly "DHGH".
C
C "MLOW" and "MHGH" are each forced to be greater than or equal to 1
C and less than or equal to 25.
C
C The function used is a sum of exponentials.
C
        DIMENSION DATA(IDIM,1),CCNT(3,50)
C
        FOVM=9./FLOAT(M)
        FOVN=9./FLOAT(N)
C
        NLOW=MAX0(1,MIN0(25,MLOW))
        NHGH=MAX0(1,MIN0(25,MHGH))
        NCNT=NLOW+NHGH
C
        DO 101 K=1,NCNT
          CCNT(1,K)=1.+(FLOAT(M)-1.)*FRAN()
          CCNT(2,K)=1.+(FLOAT(N)-1.)*FRAN()
          IF (K.LE.NLOW) THEN
            CCNT(3,K)=-1.
          ELSE
            CCNT(3,K)=+1.
          END IF
  101   CONTINUE
C
        DMIN=+1.E36
        DMAX=-1.E36
        DO 104 J=1,N
          DO 103 I=1,M
            DATA(I,J)=.5*(DLOW+DHGH)
            DO 102 K=1,NCNT
              TEMP=-((FOVM*(FLOAT(I)-CCNT(1,K)))**2+
     +               (FOVN*(FLOAT(J)-CCNT(2,K)))**2)
              IF (TEMP.GE.-20.) DATA(I,J)=DATA(I,J)+
     +            .5*(DHGH-DLOW)*CCNT(3,K)*EXP(TEMP)
  102       CONTINUE
            DMIN=AMIN1(DMIN,DATA(I,J))
            DMAX=AMAX1(DMAX,DATA(I,J))
  103     CONTINUE
  104   CONTINUE
C
        DO 106 J=1,N
          DO 105 I=1,M
            DATA(I,J)=(DATA(I,J)-DMIN)/(DMAX-DMIN)*(DHGH-DLOW)+DLOW
  105     CONTINUE
  106   CONTINUE
C
        RETURN
C
      END

      FUNCTION FRAN ()
        DIMENSION RSEQ (100)
        SAVE ISEQ
        DATA RSEQ / .749,.973,.666,.804,.081,.483,.919,.903,.951,.960 ,
     +              .039,.269,.270,.756,.222,.478,.621,.063,.550,.798 ,
     +              .027,.569,.149,.697,.451,.738,.508,.041,.266,.249 ,
     +              .019,.191,.266,.625,.492,.940,.508,.406,.972,.311 ,
     +              .757,.378,.299,.536,.619,.844,.342,.295,.447,.499 ,
     +              .688,.193,.225,.520,.954,.749,.997,.693,.217,.273 ,
     +              .961,.948,.902,.104,.495,.257,.524,.100,.492,.347 ,
     +              .981,.019,.225,.806,.678,.710,.235,.600,.994,.758 ,
     +              .682,.373,.009,.469,.203,.730,.588,.603,.213,.495 ,
     +              .884,.032,.185,.127,.010,.180,.689,.354,.372,.429 /
        DATA ISEQ / 0 /
        ISEQ=MOD(ISEQ,100)+1
        FRAN=RSEQ(ISEQ)
        RETURN
      END

      SUBROUTINE DFCLRS
C
C Define a set of RGB color triples for colors 1 through 15.
C
        DIMENSION RGBV(3,15)
C
C Define the RGB color triples needed below.
C
        DATA RGBV / 1.00 , 1.00 , 1.00 ,
     +              0.70 , 0.70 , 0.70 ,
     +              0.75 , 0.50 , 1.00 ,
     +              0.50 , 0.00 , 1.00 ,
     +              0.00 , 0.00 , 1.00 ,
     +              0.00 , 0.50 , 1.00 ,
     +              0.00 , 1.00 , 1.00 ,
     +              0.00 , 1.00 , 0.60 ,
     +              0.00 , 1.00 , 0.00 ,
     +              0.70 , 1.00 , 0.00 ,
     +              1.00 , 1.00 , 0.00 ,
     +              1.00 , 0.75 , 0.00 ,
     +              1.00 , 0.38 , 0.38 ,
     +              1.00 , 0.00 , 0.38 ,
     +              1.00 , 0.00 , 0.00 /
C
C Define 16 different color indices, for indices 0 through 15.  The
C color corresponding to index 0 is black and the color corresponding
C to index 1 is white.
C
        CALL GSCR (1,0,0.,0.,0.)
C
        DO 101 I=1,15
          CALL GSCR (1,I,RGBV(1,I),RGBV(2,I),RGBV(3,I))
  101   CONTINUE
C
C Done.
C
        RETURN
C
      END

      SUBROUTINE CPMPXY (IMAP,XINP,YINP,XOTP,YOTP)
C
C By default, CONPACK draws contour lines using X coordinates in the
C range from XAT1 to XATM and Y coordinates in the range from YAT1 to
C YATN.  Setting the flag 'MAP' non-zero causes those coordinates to be
C transformed by calling the subroutine CPMPXY.  By default, 'MAP' = 1
C selects the EZMAP transformations and 'MAP' = 2 selects the polar
C coordinate transformations.  The user of CONPACK may replace this
C subroutine as desired to transform the final coordinates and thus to
C transform the objects drawn.
C
C NOTE:  As of 4/25/91, this routine calls a new EZMAP routine, called
C MAPTRA, instead of MAPTRN.  The new routine returns 1.E12 for points
C which project outside the EZMAP perimeter.
C
C NOTE:  As of 1/14/92, this routine is being changed so that, when IMAP
C is negated, the inverse mapping is requested:  (XINP,YINP) is a point
C in the current user coordinate system; (XOTP,YOTP) is returned and is
C the point which would be carried into (XINP,YINP) by the mapping
C numbered ABS(IMAP).
C
C An additional convention has been adopted which will allow CONPACK to
C find out whether a given inverse transformation is available.  A call
C of the form
C
C       CALL CPMPXY (0,REAL(IMAP),RFLG,DUM1,DUM2)
C
C will return information in RFLG about the mapping numbered IMAP, as
C follows:
C
C   RFLG       forward mapping defined      inverse mapping defined
C   ----       -----------------------      -----------------------
C    0.                 no                            no
C    1.                yes                            no
C    2.                 no                           yes
C    3.                yes                           yes
C
C Versions of CPMPXY that have not been updated to include these new
C features should continue to work for a period of time, but ought to
C be updated eventually.
C
C Note:  As of 06/02/92, this routine is being updated to handle a
C perspective transformation.  The corner of the data array with index
C (1,1) is mapped into (XA11,YA11,ZA11), the corner with index (M,1) is
C mapped into (XA11+XDM1,YA11+YDM1,ZA11+ZDM1), and the corner with index
C (1,N) is mapped into (XA11+XD1N,YA11+YD1N,ZA11+ZD1N).  Thus, contours
C are mapped onto a parallelogram in 3-space; from there, they are
C mapped into a projection plane, using the perspective transformation
C defined by the current state of the 3-D C transformation routine
C TDPROJ (which the user has initialized by calling TDINIT).  These
C values are passed by the user in the following common block.
C
	COMMON /TDCOM1/ XA11,YA11,ZA11,XDM1,YDM1,ZDM1,XD1N,YD1N,ZD1N
C
C Sort out the different cases.  When IMAP = 0, CPMPXY is being asked
C what capabilities it has.  The IF tests distinguishing various values
C of XINP from one another makes no difference here, since inverses are
C available for all mappings done by this version; it is included as a
C model for those who pick up the routine and modify it (possibly by
C the addition of a mapping for which the inverse is not available).
C
	IF (IMAP.EQ.0) THEN
C
	  IF (INT(XINP).EQ.1) THEN
	    YINP=3.
	  ELSE IF (INT(XINP).EQ.2) THEN
	    YINP=3.
	  ELSE IF (INT(XINP).EQ.3) THEN
	    YINP=3.
	  ELSE
	    YINP=3.
	  END IF
C
C Handle the EZMAP case ...
C
	ELSE IF (ABS(IMAP).EQ.1) THEN
C
	  IF (IMAP.GT.0) THEN
	    CALL MAPTRA (YINP,XINP,XOTP,YOTP)
	  ELSE
	    CALL MAPTRI (XINP,YINP,YOTP,XOTP)
	  END IF
C
C ... the polar coordinate case ...
C
	ELSE IF (ABS(IMAP).EQ.2) THEN
C
	  IF (IMAP.GT.0) THEN
	    XOTP=XINP*COS(.017453292519943*YINP)
	    YOTP=XINP*SIN(.017453292519943*YINP)
	  ELSE
	    XOTP=SQRT(XINP*XINP+YINP*YINP)
	    YOTP=57.2957795130823*ATAN2(YINP,XINP)
	  END IF
C
C ... the layered case ...
C
	ELSE IF (ABS(IMAP).EQ.3) THEN
C
	  IF (IMAP.GT.0) THEN
	    XTMP=XA11+XINP*XDM1+YINP*XD1N
	    YTMP=YA11+XINP*YDM1+YINP*YD1N
	    ZTMP=ZA11+XINP*ZDM1+YINP*ZD1N
	    CALL TDPROJ (XTMP,YTMP,ZTMP,XOTP,YOTP)
	  ELSE
	    CALL TDINVS (XINP,YINP,XOTP,YOTP)
	  END IF
C
C ... and the identity case.
C
	ELSE
C
	  XOTP=XINP
	  YOTP=YINP
C
	END IF
C
C Done.
C
	RETURN
C
      END

      SUBROUTINE TDINIT (XEYE,YEYE,ZEYE,XOBJ,YOBJ,ZOBJ)
C
C This routine initializes the projection package.
C
	COMMON /TDCOM2/ XE,YE,ZE,XO,YO,ZO
	COMMON /TDCOM2/ A1,B1,C1,D1,A2,B2,C2,D2,E2,A3,B3,C3,D3,E3
C
C Copy the coordinates of the eye position into labelled common.
C
	XE=XEYE
	YE=YEYE
	ZE=ZEYE
C
C Copy the coordinates of the point looked at into labelled common.
C
	XO=XOBJ
	YO=YOBJ
	ZO=ZOBJ
C
C Compute the coefficients in the equation defining the projection
C plane.
C
	A1=XO-XE
	B1=YO-YE
	C1=ZO-ZE
	D1=XO*XE+YO*YE+ZO*ZE-XO*XO-YO*YO-ZO*ZO
	T1=-(A1*XE+B1*YE+C1*ZE+D1)
	A1=A1/T1
	B1=B1/T1
	C1=C1/T1
	D1=D1/T1
C
C Compute the coefficients in the equation defining a plane passing
C through the eye and the point looked at and perpendicular to the XY
C plane.
C
	A2=YO-YE
	B2=XE-XO
	C2=0.
	D2=XO*YE-XE*YO
	T2=SQRT(A2*A2+B2*B2+C2*C2)
	A2=A2/T2
	B2=B2/T2
	C2=C2/T2
	D2=D2/T2
	E2=A2*XE+B2*YE+C2*ZE+D2
C
C Compute the coefficients in the equation defining a plane passing
C through the eye and the point looked at and perpendicular to each
C of the first two planes.
C
	A3=(XO-XE)*(ZO-ZE)
	B3=(YO-YE)*(ZO-ZE)
	C3=-((XO-XE)*(XO-XE)+(YO-YE)*(YO-YE))
	D3=(XE*XE+YE*YE)*ZO+(XO*XO+YO*YO)*ZE-(XE*XO+YE*YO)*(ZE+ZO)
	T3=-SQRT(A3*A3+B3*B3+C3*C3)
	A3=A3/T3
	B3=B3/T3
	C3=C3/T3
	D3=D3/T3
	E3=A3*XE+B3*YE+C3*ZE+D3
C
C Done.
C
	RETURN
C
      END

      SUBROUTINE TDPROJ (XI3D,YI3D,ZI3D,XI2D,YI2D)
C
C This routine finds the projection of the point (XI3D,YI3D,ZI3D) in
C the projection plane and returns the distances of the projection
C from each of two planes perpendicular to that, which are the desired
C X and Y coordinates in the projection plane.
C
	COMMON /TDCOM2/ XE,YE,ZE,XO,YO,ZO
	COMMON /TDCOM2/ A1,B1,C1,D1,A2,B2,C2,D2,E2,A3,B3,C3,D3,E3
C
C What's involved is just algebra.
C
	DX=XI3D-XE
	DY=YI3D-YE
	DZ=ZI3D-ZE
C
	TT=1./(A1*DX+B1*DY+C1*DZ)
C
	XP=XE+DX*TT
	YP=YE+DY*TT
	ZP=ZE+DZ*TT
C
	XI2D=A2*XP+B2*YP+C2*ZP+D2
	YI2D=A3*XP+B3*YP+C3*ZP+D3
C
C Done.
C
	RETURN
C
      END

      SUBROUTINE TDINVS (XI2D,YI2D,XDAT,YDAT)
C
C This routine, given the projected coordinates of a point (XI2D,YI2D),
C returns the coordinates of the point (XDAT,YDAT), in the rectangle on
C which the data are defined, that projected into it.
C
C The variables in the following common block define the parallelogram
C into which the data rectangle maps.
C
	COMMON /TDCOM1/ XA11,YA11,ZA11,XDM1,YDM1,ZDM1,XD1N,YD1N,ZD1N
C
C The variables in the following common block define the mapping from
C 3-space to 2-space.
C
	COMMON /TDCOM2/ XE,YE,ZE,XO,YO,ZO
	COMMON /TDCOM2/ A1,B1,C1,D1,A2,B2,C2,D2,E2,A3,B3,C3,D3,E3
C
C What's involved is just algebra.
C
	F2=E2-XI2D
	F3=E3-YI2D
C
	A12=A1*F2+A2
	B12=B1*F2+B2
	C12=C1*F2+C2
C
	A13=A1*F3+A3
	B13=B1*F3+B3
	C13=C1*F3+C3
C
	P1=A12*XDM1+B12*YDM1+C12*ZDM1
	Q1=A12*XD1N+B12*YD1N+C12*ZD1N
	R1=A12*(XA11-XE)+B12*(YA11-YE)+C12*(ZA11-ZE)
C
	P2=A13*XDM1+B13*YDM1+C13*ZDM1
	Q2=A13*XD1N+B13*YD1N+C13*ZD1N
	R2=A13*(XA11-XE)+B13*(YA11-YE)+C13*(ZA11-ZE)
C
	XDAT=(Q1*R2-Q2*R1)/(P1*Q2-P2*Q1)
	YDAT=(P2*R1-P1*R2)/(P1*Q2-P2*Q1)
C
C Done.
C
	RETURN
C
      END

THE CELL-ARRAY VERSION:

      PROGRAM LAYERS
C
C Declare a common block for communication with the routines CPMPXY
C and TDINVS.
C
	COMMON /TDCOM1/ XA11,YA11,ZA11,XDM1,YDM1,ZDM1,XD1N,YD1N,ZD1N
C
C Declare required local data arrays and workspace arrays.
C
	DIMENSION ZDAT(23,14,5),RWRK(1000),IWRK(1000),ICRA(512,512)
        DIMENSION IASF(13)
C
C Declare arrays to hold the list of indices and the list of labels
C required by the label-bar routine.
C
        DIMENSION LIND(14)
        CHARACTER*10 LLBS(15)
C
C Define an array for GKS aspect source flags.
C
        DATA IASF / 13*1 /
C
C Define the list of indices required by the label-bar routine.
C
        DATA LIND / 2,3,4,5,6,7,8,9,10,11,12,13,14,15 /
C
C Open GKS.
C
        CALL OPNGKS
C
C Turn off the clipping indicator.
C
        CALL GSCLIP (0)
C
C Set all aspect source flags to "individual".
C
        CALL GSASF (IASF)
C
C Force solid fill.
C
        CALL GSFAIS (1)
C
C Define color indices.
C
        CALL DFCLRS
C
C Do an appropriate call to SET.
C
	CALL SET    (.10,.90,.15,.95,-.9,.9,-.9,.9,1)
C
C Initialize the 3-D to 2-D projection package.  The first three
C arguments specify the eye position, and the last three specify
C the position of the point looked at.  The line from the eye to
C the point looked at is the "line of sight".  The plane passing
C through the point looked at and perpendicular to the line of
C sight is the projection plane.
C
	CALL TDINIT (4.,3.,1.75,.5,.5,.5)
C
C Tell CONPACK we've done the call to SET.
C
	CALL CPSETI ('SET - DO-SET-CALL FLAG',0)
C
C Tell CONPACK exactly what contour levels to use.  We want 13 levels,
C splitting the range into 14 equal bands, one for each of the 14
C colors available.
C
	CALL CPSETI ('CLS - CONTOUR LEVEL SELECTOR',1)
	CALL CPSETR ('CIS - CONTOUR INTERVAL SELECTOR',1.)
	CALL CPSETR ('CMN - CONTOUR MINIMUM',1.)
	CALL CPSETR ('CMX - CONTOUR MAXIMUM',13.)
C
C Tell CONPACK to generate X and Y coordinates in the range from 0 to 1.
C
	CALL CPSETR ('XC1 - X COORDINATE AT I = 1',0.)
	CALL CPSETR ('XCM - X COORDINATE AT I = M',1.)
	CALL CPSETR ('YC1 - Y COORDINATE AT J = 1',0.)
	CALL CPSETR ('YCN - Y COORDINATE AT J = N',1.)
C
C Tell CONPACK to use mapping 3, which, in the version of CPMPXY
C supplied, is the mapping from a particular layer of the unit cube
C to the plotter frame.
C
	CALL CPSETI ('MAP - MAPPING FUNCTION',3)
C
C Tell CONPACK what out-of-range value to expect from CPMPXY.
C
        CALL CPSETR ('ORV - OUT-OF-RANGE VALUE',1.E12)
C
C Tell CONPACK how to map area identifiers into color indices when
C CPCICA is called.
C
	CALL CPSETI ('CAF - CELL ARRAY FLAG',1)
C
C Tell CONPACK what area identifier to use for areas outside the
C grid.
C
        CALL CPSETI ('PAI - PARAMETER ARRAY INDEX',-1)
	CALL CPSETI ('AIA - AREA IDENTIFIER OUTSIDE THE GRID',-1)
C
C Spin up the random number generator.
C
	DO 100 I=1,123
	  TEMP=FRAN()
  100   CONTINUE
C
C Generate data at the bottom and top levels; we'll interpolate at
C the intermediate levels.
C
	CALL GENDAT (ZDAT(1,1,1),23,23,14,20,20,0.,14.)
	CALL GENDAT (ZDAT(1,1,5),23,23,14,20,20,0.,14.)
C
C Interpolate data at the intermediate levels.
C
	DO 103 ILEV=2,4
	  VLEV=.25*REAL(ILEV-1)
	  DO 102 I=1,23
	    DO 101 J=1,14
	      ZDAT(I,J,ILEV)=(1.-VLEV)*ZDAT(I,J,1)+VLEV*ZDAT(I,J,5)
  101       CONTINUE
  102     CONTINUE
  103   CONTINUE
C
C Define the corners of the parallelogram containing the contours.
C
	XA11=0.
	YA11=0.
	XDM1=1.
	YDM1=0.
	XD1N=0.
	YD1N=1.
	ZDM1=0.
	ZD1N=0.
C
C Initialize the contents of the cell array.
C
	DO 105 I=1,512
	  DO 104 J=1,512
	    ICRA(I,J)=0
  104     CONTINUE
  105   CONTINUE
C
C Do five plots, each at a different level in the unit cube.
C
	DO 106 ILEV=1,5
C
C Set the vertical level of the contour plot.
C
	  ZA11=.25*REAL(ILEV-1)
C
C Initialize the drawing of the contour plot.
C
	  CALL CPRECT (ZDAT(1,1,ILEV),23,23,14,RWRK,1000,IWRK,1000)
C
C Add color indices to the cell array.
C
	  CALL CPCICA (ZDAT(1,1,ILEV),RWRK,IWRK,ICRA,512,512,512,
     +                                             .1,.15,.9,.95)
C
  106   CONTINUE
C
C Plot the cell array.
C
	CALL GCA (CFUX(.1),CFUY(.15),CFUX(.9),CFUY(.95),512,512,1,1,
     +                                                 512,512,ICRA)
C
C Draw a label bar for the plot, relating colors to values.
C
        CALL CPGETR ('ZMN',ZMIN)
        CALL CPGETR ('ZMX',ZMAX)
C
	DO 107 I=1,15
          CALL CPSETR ('ZDV - Z DATA VALUE',
     +                 ZMIN+REAL(I-1)*(ZMAX-ZMIN)/14.)
          CALL CPGETC ('ZDV - Z DATA VALUE',LLBS(I))
  107   CONTINUE
C
        CALL LBSETI ('CBL - COLOR OF BOX LINES',0)
	CALL LBLBAR (0,.05,.95,.05,.15,14,1.,.5,LIND,0,LLBS,15,1)
C
C Put a label at the top of the plot.
C
	CALL SFLUSH
	CALL GSLWSC (2.)
	CALL PLCHHQ (CFUX(.5),CFUY(.95),'LAYERED CONTOURS',.025,0.,0.)
C
C Advance the frame.
C
        CALL FRAME
C
C Close GKS.
C
        CALL CLSGKS
C
C Done.
C
        STOP
C
      END

      SUBROUTINE GENDAT (DATA,IDIM,M,N,MLOW,MHGH,DLOW,DHGH)
C
C This is a routine to generate test data for two-dimensional graphics
C routines.  Given an array "DATA", dimensioned "IDIM x 1", it fills
C the sub-array ((DATA(I,J),I=1,M),J=1,N) with a two-dimensional field
C of data having approximately "MLOW" lows and "MHGH" highs, a minimum
C value of exactly "DLOW" and a maximum value of exactly "DHGH".
C
C "MLOW" and "MHGH" are each forced to be greater than or equal to 1
C and less than or equal to 25.
C
C The function used is a sum of exponentials.
C
        DIMENSION DATA(IDIM,1),CCNT(3,50)
C
        FOVM=9./FLOAT(M)
        FOVN=9./FLOAT(N)
C
        NLOW=MAX0(1,MIN0(25,MLOW))
        NHGH=MAX0(1,MIN0(25,MHGH))
        NCNT=NLOW+NHGH
C
        DO 101 K=1,NCNT
          CCNT(1,K)=1.+(FLOAT(M)-1.)*FRAN()
          CCNT(2,K)=1.+(FLOAT(N)-1.)*FRAN()
          IF (K.LE.NLOW) THEN
            CCNT(3,K)=-1.
          ELSE
            CCNT(3,K)=+1.
          END IF
  101   CONTINUE
C
        DMIN=+1.E36
        DMAX=-1.E36
        DO 104 J=1,N
          DO 103 I=1,M
            DATA(I,J)=.5*(DLOW+DHGH)
            DO 102 K=1,NCNT
              TEMP=-((FOVM*(FLOAT(I)-CCNT(1,K)))**2+
     +               (FOVN*(FLOAT(J)-CCNT(2,K)))**2)
              IF (TEMP.GE.-20.) DATA(I,J)=DATA(I,J)+
     +            .5*(DHGH-DLOW)*CCNT(3,K)*EXP(TEMP)
  102       CONTINUE
            DMIN=AMIN1(DMIN,DATA(I,J))
            DMAX=AMAX1(DMAX,DATA(I,J))
  103     CONTINUE
  104   CONTINUE
C
        DO 106 J=1,N
          DO 105 I=1,M
            DATA(I,J)=(DATA(I,J)-DMIN)/(DMAX-DMIN)*(DHGH-DLOW)+DLOW
  105     CONTINUE
  106   CONTINUE
C
        RETURN
C
      END

      FUNCTION FRAN ()
        DIMENSION RSEQ (100)
        SAVE ISEQ
        DATA RSEQ / .749,.973,.666,.804,.081,.483,.919,.903,.951,.960 ,
     +              .039,.269,.270,.756,.222,.478,.621,.063,.550,.798 ,
     +              .027,.569,.149,.697,.451,.738,.508,.041,.266,.249 ,
     +              .019,.191,.266,.625,.492,.940,.508,.406,.972,.311 ,
     +              .757,.378,.299,.536,.619,.844,.342,.295,.447,.499 ,
     +              .688,.193,.225,.520,.954,.749,.997,.693,.217,.273 ,
     +              .961,.948,.902,.104,.495,.257,.524,.100,.492,.347 ,
     +              .981,.019,.225,.806,.678,.710,.235,.600,.994,.758 ,
     +              .682,.373,.009,.469,.203,.730,.588,.603,.213,.495 ,
     +              .884,.032,.185,.127,.010,.180,.689,.354,.372,.429 /
        DATA ISEQ / 0 /
        ISEQ=MOD(ISEQ,100)+1
        FRAN=RSEQ(ISEQ)
        RETURN
      END

      SUBROUTINE DFCLRS
C
C Define a set of RGB color triples for colors 1 through 15.
C
        DIMENSION RGBV(3,15)
C
C Define the RGB color triples needed below.
C
        DATA RGBV / 1.00 , 1.00 , 1.00 ,
     +              0.70 , 0.70 , 0.70 ,
     +              0.75 , 0.50 , 1.00 ,
     +              0.50 , 0.00 , 1.00 ,
     +              0.00 , 0.00 , 1.00 ,
     +              0.00 , 0.50 , 1.00 ,
     +              0.00 , 1.00 , 1.00 ,
     +              0.00 , 1.00 , 0.60 ,
     +              0.00 , 1.00 , 0.00 ,
     +              0.70 , 1.00 , 0.00 ,
     +              1.00 , 1.00 , 0.00 ,
     +              1.00 , 0.75 , 0.00 ,
     +              1.00 , 0.38 , 0.38 ,
     +              1.00 , 0.00 , 0.38 ,
     +              1.00 , 0.00 , 0.00 /
C
C Define 16 different color indices, for indices 0 through 15.  The
C color corresponding to index 0 is black and the color corresponding
C to index 1 is white.
C
        CALL GSCR (1,0,0.,0.,0.)
C
        DO 101 I=1,15
          CALL GSCR (1,I,RGBV(1,I),RGBV(2,I),RGBV(3,I))
  101   CONTINUE
C
C Done.
C
        RETURN
C
      END

      SUBROUTINE CPMPXY (IMAP,XINP,YINP,XOTP,YOTP)
C
C By default, CONPACK draws contour lines using X coordinates in the
C range from XAT1 to XATM and Y coordinates in the range from YAT1 to
C YATN.  Setting the flag 'MAP' non-zero causes those coordinates to be
C transformed by calling the subroutine CPMPXY.  By default, 'MAP' = 1
C selects the EZMAP transformations and 'MAP' = 2 selects the polar
C coordinate transformations.  The user of CONPACK may replace this
C subroutine as desired to transform the final coordinates and thus to
C transform the objects drawn.
C
C NOTE:  As of 4/25/91, this routine calls a new EZMAP routine, called
C MAPTRA, instead of MAPTRN.  The new routine returns 1.E12 for points
C which project outside the EZMAP perimeter.
C
C NOTE:  As of 1/14/92, this routine is being changed so that, when IMAP
C is negated, the inverse mapping is requested:  (XINP,YINP) is a point
C in the current user coordinate system; (XOTP,YOTP) is returned and is
C the point which would be carried into (XINP,YINP) by the mapping
C numbered ABS(IMAP).
C
C An additional convention has been adopted which will allow CONPACK to
C find out whether a given inverse transformation is available.  A call
C of the form
C
C       CALL CPMPXY (0,REAL(IMAP),RFLG,DUM1,DUM2)
C
C will return information in RFLG about the mapping numbered IMAP, as
C follows:
C
C   RFLG       forward mapping defined      inverse mapping defined
C   ----       -----------------------      -----------------------
C    0.                 no                            no
C    1.                yes                            no
C    2.                 no                           yes
C    3.                yes                           yes
C
C Versions of CPMPXY that have not been updated to include these new
C features should continue to work for a period of time, but ought to
C be updated eventually.
C
C Note:  As of 06/02/92, this routine is being updated to handle a
C perspective transformation.  The corner of the data array with index
C (1,1) is mapped into (XA11,YA11,ZA11), the corner with index (M,1) is
C mapped into (XA11+XDM1,YA11+YDM1,ZA11+ZDM1), and the corner with index
C (1,N) is mapped into (XA11+XD1N,YA11+YD1N,ZA11+ZD1N).  Thus, contours
C are mapped onto a parallelogram in 3-space; from there, they are
C mapped into a projection plane, using the perspective transformation
C defined by the current state of the 3-D C transformation routine
C TDPROJ (which the user has initialized by calling TDINIT).  These
C values are passed by the user in the following common block.
C
	COMMON /TDCOM1/ XA11,YA11,ZA11,XDM1,YDM1,ZDM1,XD1N,YD1N,ZD1N
C
C Sort out the different cases.  When IMAP = 0, CPMPXY is being asked
C what capabilities it has.  The IF tests distinguishing various values
C of XINP from one another makes no difference here, since inverses are
C available for all mappings done by this version; it is included as a
C model for those who pick up the routine and modify it (possibly by
C the addition of a mapping for which the inverse is not available).
C
	IF (IMAP.EQ.0) THEN
C
	  IF (INT(XINP).EQ.1) THEN
	    YINP=3.
	  ELSE IF (INT(XINP).EQ.2) THEN
	    YINP=3.
	  ELSE IF (INT(XINP).EQ.3) THEN
	    YINP=3.
	  ELSE
	    YINP=3.
	  END IF
C
C Handle the EZMAP case ...
C
	ELSE IF (ABS(IMAP).EQ.1) THEN
C
	  IF (IMAP.GT.0) THEN
	    CALL MAPTRA (YINP,XINP,XOTP,YOTP)
	  ELSE
	    CALL MAPTRI (XINP,YINP,YOTP,XOTP)
	  END IF
C
C ... the polar coordinate case ...
C
	ELSE IF (ABS(IMAP).EQ.2) THEN
C
	  IF (IMAP.GT.0) THEN
	    XOTP=XINP*COS(.017453292519943*YINP)
	    YOTP=XINP*SIN(.017453292519943*YINP)
	  ELSE
	    XOTP=SQRT(XINP*XINP+YINP*YINP)
	    YOTP=57.2957795130823*ATAN2(YINP,XINP)
	  END IF
C
C ... the layered case ...
C
	ELSE IF (ABS(IMAP).EQ.3) THEN
C
	  IF (IMAP.GT.0) THEN
	    XTMP=XA11+XINP*XDM1+YINP*XD1N
	    YTMP=YA11+XINP*YDM1+YINP*YD1N
	    ZTMP=ZA11+XINP*ZDM1+YINP*ZD1N
	    CALL TDPROJ (XTMP,YTMP,ZTMP,XOTP,YOTP)
	  ELSE
	    CALL TDINVS (XINP,YINP,XOTP,YOTP)
	  END IF
C
C ... and the identity case.
C
	ELSE
C
	  XOTP=XINP
	  YOTP=YINP
C
	END IF
C
C Done.
C
	RETURN
C
      END

      SUBROUTINE TDINIT (XEYE,YEYE,ZEYE,XOBJ,YOBJ,ZOBJ)
C
C This routine initializes the projection package.
C
	COMMON /TDCOM2/ XE,YE,ZE,XO,YO,ZO
	COMMON /TDCOM2/ A1,B1,C1,D1,A2,B2,C2,D2,E2,A3,B3,C3,D3,E3
C
C Copy the coordinates of the eye position into labelled common.
C
	XE=XEYE
	YE=YEYE
	ZE=ZEYE
C
C Copy the coordinates of the point looked at into labelled common.
C
	XO=XOBJ
	YO=YOBJ
	ZO=ZOBJ
C
C Compute the coefficients in the equation defining the projection
C plane.
C
	A1=XO-XE
	B1=YO-YE
	C1=ZO-ZE
	D1=XO*XE+YO*YE+ZO*ZE-XO*XO-YO*YO-ZO*ZO
	T1=-(A1*XE+B1*YE+C1*ZE+D1)
	A1=A1/T1
	B1=B1/T1
	C1=C1/T1
	D1=D1/T1
C
C Compute the coefficients in the equation defining a plane passing
C through the eye and the point looked at and perpendicular to the XY
C plane.
C
	A2=YO-YE
	B2=XE-XO
	C2=0.
	D2=XO*YE-XE*YO
	T2=SQRT(A2*A2+B2*B2+C2*C2)
	A2=A2/T2
	B2=B2/T2
	C2=C2/T2
	D2=D2/T2
	E2=A2*XE+B2*YE+C2*ZE+D2
C
C Compute the coefficients in the equation defining a plane passing
C through the eye and the point looked at and perpendicular to each
C of the first two planes.
C
	A3=(XO-XE)*(ZO-ZE)
	B3=(YO-YE)*(ZO-ZE)
	C3=-((XO-XE)*(XO-XE)+(YO-YE)*(YO-YE))
	D3=(XE*XE+YE*YE)*ZO+(XO*XO+YO*YO)*ZE-(XE*XO+YE*YO)*(ZE+ZO)
	T3=-SQRT(A3*A3+B3*B3+C3*C3)
	A3=A3/T3
	B3=B3/T3
	C3=C3/T3
	D3=D3/T3
	E3=A3*XE+B3*YE+C3*ZE+D3
C
C Done.
C
	RETURN
C
      END

      SUBROUTINE TDPROJ (XI3D,YI3D,ZI3D,XI2D,YI2D)
C
C This routine finds the projection of the point (XI3D,YI3D,ZI3D) in
C the projection plane and returns the distances of the projection
C from each of two planes perpendicular to that, which are the desired
C X and Y coordinates in the projection plane.
C
	COMMON /TDCOM2/ XE,YE,ZE,XO,YO,ZO
	COMMON /TDCOM2/ A1,B1,C1,D1,A2,B2,C2,D2,E2,A3,B3,C3,D3,E3
C
C What's involved is just algebra.
C
	DX=XI3D-XE
	DY=YI3D-YE
	DZ=ZI3D-ZE
C
	TT=1./(A1*DX+B1*DY+C1*DZ)
C
	XP=XE+DX*TT
	YP=YE+DY*TT
	ZP=ZE+DZ*TT
C
	XI2D=A2*XP+B2*YP+C2*ZP+D2
	YI2D=A3*XP+B3*YP+C3*ZP+D3
C
C Done.
C
	RETURN
C
      END

      SUBROUTINE TDINVS (XI2D,YI2D,XDAT,YDAT)
C
C This routine, given the projected coordinates of a point (XI2D,YI2D),
C returns the coordinates of the point (XDAT,YDAT), in the rectangle on
C which the data are defined, that projected into it.
C
C The variables in the following common block define the parallelogram
C into which the data rectangle maps.
C
	COMMON /TDCOM1/ XA11,YA11,ZA11,XDM1,YDM1,ZDM1,XD1N,YD1N,ZD1N
C
C The variables in the following common block define the mapping from
C 3-space to 2-space.
C
	COMMON /TDCOM2/ XE,YE,ZE,XO,YO,ZO
	COMMON /TDCOM2/ A1,B1,C1,D1,A2,B2,C2,D2,E2,A3,B3,C3,D3,E3
C
C What's involved is just algebra.
C
	F2=E2-XI2D
	F3=E3-YI2D
C
	A12=A1*F2+A2
	B12=B1*F2+B2
	C12=C1*F2+C2
C
	A13=A1*F3+A3
	B13=B1*F3+B3
	C13=C1*F3+C3
C
	P1=A12*XDM1+B12*YDM1+C12*ZDM1
	Q1=A12*XD1N+B12*YD1N+C12*ZD1N
	R1=A12*(XA11-XE)+B12*(YA11-YE)+C12*(ZA11-ZE)
C
	P2=A13*XDM1+B13*YDM1+C13*ZDM1
	Q2=A13*XD1N+B13*YD1N+C13*ZD1N
	R2=A13*(XA11-XE)+B13*(YA11-YE)+C13*(ZA11-ZE)
C
	XDAT=(Q1*R2-Q2*R1)/(P1*Q2-P2*Q1)
	YDAT=(P2*R1-P1*R2)/(P1*Q2-P2*Q1)
C
C Done.
C
	RETURN
C
      END