Subject: Layered contouring.

From: [email-masked] (Dave Kennison)
Date: Mon, 1 Jun 92 19:42:25 MDT

The code is not done yet, but I have a first cut at layered contours.  The
code for an example is attached.  The version of CPMPXY in this example 
will allow you to do color fill of the layers using CPCLAM and AREAS, but
it will not allow you to use the new routine CPCICA.  I am going to try
to get that done tonight or tomorrow, but I probably won't make it, and
it will then be Thursday before I get back to it.  

A lot more needs to be done with this.  It will probably become a part 
of the 3.2 version of CONPACK (in some fashion, if only as an example).

I had always thought it was trivial to make CONPACK do this.  It wasn't
bad, but it wasn't all that trivial, either.

EXAMPLE FOLLOWS:

      PROGRAM LAYERS
C
C Declare a common block for communication with the routine CPMPXY.
C
	COMMON /MYCOMN/ VLEV
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 XXINIT (4.,3.,2.,.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,111
	  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 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
	  VLEV=.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 Declare a common block in which the main program can communicate
C the vertical level at which the current layer sits.
C
	COMMON /MYCOMN/ VLEV
C
C Sort out the different cases.
C
	IF (IMAP.EQ.0) THEN
	  IF (INT(XINP).EQ.1.OR.INT(XINP).EQ.2) THEN
	    YINP=3.
	  ELSE IF (INT(XINP).EQ.3) THEN
	    YINP=1.
	  ELSE
	    YINP=3.
	  END IF
C
C Handle the EZMAP case ...
C
	ELSE IF (ABS(IMAP).EQ.1) THEN
	  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
	  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
	  CALL XXPROJ (XINP,YINP,VLEV,XOTP,YOTP)
C
C ... and the identity case.
C
	ELSE
	  XOTP=XINP
	  YOTP=YINP
	END IF

C
C Done.
C
	RETURN
C
      END

      SUBROUTINE XXINIT (XEYE,YEYE,ZEYE,XOBJ,YOBJ,ZOBJ)
C
C This routine initializes the projection package.
C
	COMMON /XXCOMN/ XE,YE,ZE,XO,YO,ZO
	COMMON /XXCOMN/ A1,B1,C1,D1,T1,A2,B2,C2,D2,T2,A3,B3,C3,D3,T3
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)
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)
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)
C
C Done.
C
	RETURN
      END

      SUBROUTINE XXPROJ (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 /XXCOMN/ XE,YE,ZE,XO,YO,ZO
	COMMON /XXCOMN/ A1,B1,C1,D1,T1,A2,B2,C2,D2,T2,A3,B3,C3,D3,T3
C
C What's involved is just algebra.
C
	DX=XI3D-XE
	DY=YI3D-YE
	DZ=ZI3D-ZE
C
	TT=T1/(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)/T2
	YI2D=(A3*XP+B3*YP+C3*ZP+D3)/T3
C
C Done.
C
	RETURN
C
      END