C
C $Id: TdpackDemo.f,v 1.7 2006/11/10 21:52:14 kennison Exp $
C                                                                      
C                Copyright (C)  2000
C        University Corporation for Atmospheric Research
C                All Rights Reserved
C
C This file is free software; you can redistribute it and/or modify
C it under the terms of the GNU General Public License as published
C by the Free Software Foundation; either version 2 of the License, or
C (at your option) any later version.
C
C This software is distributed in the hope that it will be useful, but
C WITHOUT ANY WARRANTY; without even the implied warranty of
C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
C General Public License for more details.
C
C You should have received a copy of the GNU General Public License
C along with this software; if not, write to the Free Software
C Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
C USA.
C
      PROGRAM TDDEMO
C
C Declare parameters specifying the maximum sizes of the data-defining
C arrays and the size of the triangle-list array.
C
        PARAMETER (IMAX=81,JMAX=81,KMAX=81,MTRI=80000)
C
C Declare local dimensioned variables to hold data defining a simple
C surface or an isosurface.
C
        DIMENSION U(IMAX),V(JMAX),W(KMAX),S(IMAX,JMAX),F(IMAX,JMAX,KMAX)
C
C Declare the local dimensioned variable in which the triangle list is
C to be constructed and a couple of temporary variables used in sorting
C the list.
C
        DIMENSION RTRI(10,MTRI),RTWK(MTRI,2),ITWK(MTRI)
C
C Declare some character variables into which to encode numeric and
C informational labels for the axes.
C
        CHARACTER*128 UNLB,VNLB,WNLB,UILB,VILB,WILB,PILB
C
C Declare some character variables in which formats for encoding the
C numeric labels are kept.
C
        CHARACTER*16 UFMT,VFMT,WFMT
C
C Define a character variable into which to read a file name.
C
        CHARACTER*128 FLNM
C
C Define a character variable into which to read commands.
C
        CHARACTER*1 COMD
C
C Declare a variable to hold the names of the various mark types.
C
        CHARACTER*15 IMKT(5)
C
C Define the conversion constant from degrees to radians.
C
        DATA DTOR / .017453292519943 /
C
C Set the default handedness of the coordinate system.  This must match
C the default value of the TDPACK internal parameter 'HND'.  Zero means
C a right-handed system, one a left-handed system.
C
        DATA IHND / 0 /
C
C Set the default value of the parameter that says whether or not the
C X workstation will be updated as changes are made.  (When making a
C movie, it's a good idea to set this to zero.)
C
        DATA IUXW / 1 /
C
C Set the default value of the surface selector.
C
        DATA ISAS / -1 /
C
C Set the default values of the dimensions of the data arrays.
C
        DATA IDIM,JDIM,KDIM / 21,21,21 /
C
C Set the default value of the flag that says whether the basic color
C scheme will be white on black (IBOW=0) or black on white (IBOW=1).
C
        DATA IBOW / 0 /
C
C Set the default values of parameters determining the eye position, the
C position of the point looked at, and the field of view.
C
        DATA ANG1,ANG2,REYE,RAPT,FOVP / -35.,25.,2.9,0.,20. /
C
C Set the default value of the TDPACK character multiplier.
C
        DATA RCS1 / 1.25 /
C
C Set the default value of the parameter that says to compute new center
C positions for the exponentials used to generate dummy surfaces and
C isosurfaces.  Basically, with ICEN = 0, the same surface will be
C generated by consecutive calls to GENDAT or to TDGNDT, even if the
C dimensions of the data array are changed, but, with ICEN = 1, a
C different surface will be generated each time.
C
        DATA ICEN / 0 /
C
C Set the default cut-off values for generated isosurfaces.
C
        DATA FIS1,FIS2 / 95.,95. /
C
C Set the default value of the multiplier for the random number term
C in the code used to compute the surface/isosurface data.
C
        DATA RNDM / 0. /
C
C Set the default value of the parameter that says what kinds of grids
C are to be drawn.
C
        DATA IGDF / 13 /
C
C Set the default value of parameters that specify the grid spacing.
C
        DATA UGIN,VGIN,WGIN / 1.,1.,1. /
C
C Set the default value of parameters that specify the label spacing.
C
        DATA ULIN,VLIN,WLIN / 5.,5.,5. /
C
C Set the default values of parameters specifying how the numeric labels
C are to be encoded.
C
        DATA UFMT / '(21(1X,F5.1))' /
        DATA VFMT / '(21(1X,F5.1))' /
        DATA WFMT / '(21(1X,F5.1))' /
C
C Set the default values of the U, V, and W axis labels and the plot
C label.
C
        DATA UILB / 'U Coordinate Values' /
        DATA VILB / 'V Coordinate Values' /
        DATA WILB / 'W Coordinate Values' /
C
        DATA PILB / ' ' /
C
C Set the default value of the parameter that controls fading of axis
C labels for rotation angles near multiples of 90 degrees.
C
        DATA IFDE / 0 /
C
C Set the default values of parameters used to position the plot label.
C
        DATA PILX,PILY / 0. , .54 /
C
C Set the default value of the flag that says whether or not to draw a
C 3D grid of marks.
C
        DATA IMRK / 0 /
C
C Set the default value of the mark size parameter.
C
        DATA SMRK / .025 /
C
C Set the default value of parameters that specify the mark spacing.
C
        DATA UMKS,VMKS,WMKS / 5.,5.,5. /
C
C Set the default value of the parameter that determines how the
C triangles are ordered.  This will be used as the final argument in
C calls to TDOTRI.
C
        DATA IORD / 0 /
C
C Define default values of the shading parameters.
C
        DATA SHDE,SHDR,ISHD,ANG3,ANG4
     +       /.15,  1.,   0,  0.,  0./
C
C Define default values of the rendering-style parameters.
C
        DATA IFCB,IFCT,ILCB,ILCT,IDTR,USIN,VSIN,WSIN
     +       /  4,   2, 116, 116,   0,  1.,  1.,  1./
C
C Set a couple of counters that keep track of how many NCGM and
C PostScript frames have been saved.
C
        DATA NFNG,NFPS / 0,0 /
C
C Set the default value of the stereo-view flag.  The user can set it
C negative to get stereo views on a single frame or positive to get
C stereo views on separate frames.  In either case, its absolute value
C is the angle, in degrees, between the lines of sight of the two views.
C
        DATA STVF / 0. /
C
C When STVF is set negative. WOSW is the width of each stereo window,
C as a fraction of the width of the plotter frame.  If this value is
C set greater than .5, the windows will overlap slightly.
C
        DATA WOSW / .5 /
C
C Define the names of the various mark types.
C
        DATA IMKT / 'a tetrahedron.','an octahedron.','a cube.',
     +                         'an icosahedron.', 'a "sphere".'/
C
C Open GKS.
C
        CALL GOPKS (6,0)
C
C Open and activate an X workstation.
C
        CALL GOPWK (1,0,8)
        CALL GACWK (1)
C
C Open and activate an NCGM workstation.
C
        CALL GOPWK (2,0,1)
        CALL GACWK (2)
C
C Open and activate a PostScript workstation.
C
        CALL GOPWK (3,0,20)
        CALL GACWK (3)
C
C Turn clipping off.
C
        CALL GSCLIP (0)
C
C Select a "fill area interior style" of "solid".
C
        CALL GSFAIS (1)
C
C Double the line width.
C
        CALL GSLWSC (2.)
C
C Define colors on all the workstations.
C
        DO 102 IWID=1,3
          CALL TDCLRS (IWID,IBOW,SHDE,SHDR,101,116,0)
  102   CONTINUE
C
C Save the values of SHDE and SHDR.
C
        SSHE=SHDE
        SSHR=SHDR
C
C Deactivate the NCGM and PostScript workstations.
C
        CALL GDAWK (2)
        CALL GDAWK (3)
C
C Select font number 25, turn on the outlining of filled fonts, set the
C line width to 1, and turn off the setting of the outline color.
C
        CALL PCSETI ('FN - FONT NUMBER',25)
        CALL PCSETI ('OF - OUTLINE FLAG',1)
        CALL PCSETR ('OL - OUTLINE LINE WIDTH',1.)
        CALL PCSETR ('OC - OUTLINE LINE COLOR',-1.)
C
C Compute values defining the surface to be depicted.
C
  103   PRINT * , 'Generating U data.'
C
        DO 104 I=1,IDIM
          U(I)=REAL(I-1)
  104   CONTINUE
C
        UMIN=U(   1)
        UMAX=U(IDIM)
        UMID=.5*(UMIN+UMAX)
C
        PRINT * , 'Generating V data.'
C
        DO 105 J=1,JDIM
          V(J)=REAL(J-1)
  105   CONTINUE
C
        VMIN=V(   1)
        VMAX=V(JDIM)
        VMID=.5*(VMIN+VMAX)
C
        PRINT * , 'Generating W data.'
C
        DO 106 K=1,KDIM
          W(K)=REAL(K-1)
  106   CONTINUE
C
        WMIN=W(   1)
        WMAX=W(KDIM)
        WMID=.5*(WMIN+WMAX)
C
C Generate data representing a simple surface, a sphere, a torus, two
C interlocking tori, three joined tori, or any of the above with some
C added random noise, depending on the value of ISAS.
C
        IF (ISAS.LT.0) THEN
          PRINT * , 'Generating a simple random surface.'
          CALL GENDAT (S,IMAX,IDIM,JDIM,23,23,WMIN,WMAX,ICEN)
          IF (RNDM.NE.0.) CALL ADDRAN (F,IMAX,IDIM,JDIM,WMIN,WMAX,RNDM)
        ELSE IF (ISAS.EQ.0) THEN
          PRINT * , 'Generating a random isosurface.'
          CALL TDGNDT (F,IMAX,JMAX,IDIM,JDIM,KDIM,23,23,0.,190.,ICEN)
          IF (RNDM.NE.0.) CALL TDADRN (F,IMAX,JMAX,IDIM,JDIM,KDIM,
     +                                               0.,190.,RNDM)
        ELSE IF (ISAS.GE.1.AND.ISAS.LE.5) THEN
          PRINT * , 'Generating a canned isosurface.'
          FMIN=+1.E36
          FMAX=-1.E36
          DO 109 I=1,IDIM
            UTMP=1.-2.*(REAL(IDIM-I)/REAL(IDIM-1))
            DO 108 J=1,JDIM
              VTMP=1.-2.*(REAL(JDIM-J)/REAL(JDIM-1))
              DO 107 K=1,KDIM
                WTMP=1.-2.*(REAL(KDIM-K)/REAL(KDIM-1))
                IF      (ISAS.EQ.1) THEN
                  F(I,J,K)=100.*SQRT(UTMP**2+VTMP**2+WTMP**2)
                ELSE IF (ISAS.EQ.2) THEN
                  F(I,J,K)=380.*
     +                     SQRT((SQRT(UTMP**2+VTMP**2)-.7)**2+WTMP**2)
                ELSE IF (ISAS.EQ.3) THEN
                  F(I,J,K)=380.*
     +                     MIN(SQRT((SQRT((1.35*UTMP+.35)**2+
     +                          (1.35*VTMP)**2)-.7)**2+(1.35*WTMP)**2),
     +                         SQRT((SQRT((1.35*UTMP-.35)**2+
     +                          (1.35*WTMP)**2)-.7)**2+(1.35*VTMP)**2))
                ELSE IF (ISAS.EQ.4) THEN
                  F(I,J,K)=760.*MIN(
     +                     SQRT((SQRT(UTMP**2+VTMP**2)-.7)**2+WTMP**2),
     +                     SQRT((SQRT(WTMP**2+UTMP**2)-.7)**2+VTMP**2),
     +                     SQRT((SQRT(VTMP**2+WTMP**2)-.7)**2+UTMP**2))
                ELSE IF (ISAS.EQ.5) THEN
                  F(I,J,K)=4.*FIS1*(UTMP**2+VTMP**2-.75*WTMP**2)
                END IF
                FMIN=MIN(FMIN,F(I,J,K))
                FMAX=MAX(FMAX,F(I,J,K))
  107         CONTINUE
  108       CONTINUE
  109     CONTINUE
          IF (RNDM.NE.0.) CALL TDADRN (F,IMAX,JMAX,IDIM,JDIM,KDIM,
     +                                             FMIN,FMAX,RNDM)
        END IF
C
C Get triangles representing the surface, to be rendered using rendering
C style 1.
C
  110   PRINT * , 'Generating triangles representing the surface.'
C
        NTRI=0
C
        IF (ISAS.LT.0.OR.(ISAS.GE.99.AND.ISAS.LE.105)) THEN
          CALL TDSTRI (U,IDIM,V,JDIM,S,IMAX,RTRI,MTRI,NTRI,1)
          IF (NTRI.LT.MTRI) THEN
            PRINT * , 'Percentage of triangle space used by TDSTRI:',
     +                                INT(100.*REAL(NTRI)/REAL(MTRI))
          ELSE
            PRINT * , 'Triangle space overflow on call to TDSTRI.'
          END IF
        ELSE
          IF (FIS1.EQ.FIS2) THEN
            CALL TDITRI (U,IDIM,V,JDIM,W,KDIM,F,IMAX,JMAX,FIS1,
     +                                        RTRI,MTRI,NTRI,1)
          ELSE
            CALL TDITRI (U,IDIM,V,JDIM,W,KDIM,F,IMAX,JMAX,FIS1,
     +                                        RTRI,MTRI,NTRI,2)
            CALL TDITRI (U,IDIM,V,JDIM,W,KDIM,F,IMAX,JMAX,FIS2,
     +                                        RTRI,MTRI,NTRI,3)
          END IF
          IF (NTRI.LT.MTRI) THEN
            PRINT * , 'Percentage of triangle space used by TDITRI:',
     +                                INT(100.*REAL(NTRI)/REAL(MTRI))
          ELSE
            PRINT * , 'Triangle space overflow on call to TDITRI.'
          END IF
        END IF
C
C If it's been requested, put a 3D grid of marks inside the viewing box.
C
        IF (IMRK.NE.0) THEN
C
          IF (UMKS.GT.0..AND.UMAX-UMIN.LE.20.*UMKS) THEN
            UEPS=1.E-3*(UMAX-UMIN)
            IST1=INT((UMIN-UEPS)/UMKS+.5+SIGN(.5,UMIN-UEPS))
            IST2=INT((UMAX+UEPS)/UMKS-.5+SIGN(.5,UMAX+UEPS))
          ELSE
            IST1=1
            IST2=0
          END IF
C
          IF (VMKS.GT.0..AND.VMAX-VMIN.LE.20.*VMKS) THEN
            VEPS=1.E-3*(VMAX-VMIN)
            JST1=INT((VMIN-VEPS)/VMKS+.5+SIGN(.5,VMIN-VEPS))
            JST2=INT((VMAX+VEPS)/VMKS-.5+SIGN(.5,VMAX+VEPS))
          ELSE
            JST1=1
            JST2=0
          END IF
C
          IF (WMKS.GT.0..AND.WMAX-WMIN.LE.20.*WMKS) THEN
            WEPS=1.E-3*(WMAX-WMIN)
            KST1=INT((WMIN-WEPS)/WMKS+.5+SIGN(.5,WMIN-WEPS))
            KST2=INT((WMAX+WEPS)/WMKS-.5+SIGN(.5,WMAX+WEPS))
          ELSE
            KST1=1
            KST2=0
          END IF
C
          IF (IST1.LE.IST2.AND.JST1.LE.JST2.AND.KST1.LE.KST2) THEN
            SIZE=SMRK*MIN(UMAX-UMIN,VMAX-VMIN,WMAX-WMIN)
            DO 903 I=IST1,IST2
              UVAL=UMKS*REAL(I)
              DO 902 J=JST1,JST2
                VVAL=VMKS*REAL(J)
                DO 901 K=KST1,KST2
                  WVAL=WMKS*REAL(K)
                  CALL TDMTRI (IMRK,UVAL,VVAL,WVAL,SIZE,
     +                         RTRI,MTRI,NTRI,4,
     +                         UMIN,VMIN,WMIN,UMAX,VMAX,WMAX)
  901           CONTINUE
  902         CONTINUE
  903       CONTINUE
          ELSE
            PRINT * , 'Omitting marks - either one of the specified'
            PRINT * , 'mark intervals is non-positive or the marks'
            PRINT * , 'would be too crowded.'
          END IF
C
        END IF
C
C Tell the user what the minima and maxima are.
C
  111   PRINT * , 'U minimum and maximum: ',UMIN,UMAX
        PRINT * , 'V minimum and maximum: ',VMIN,VMAX
        PRINT * , 'W minimum and maximum: ',WMIN,WMAX
C
C Define numeric and informational labels for the plot.
C
        PRINT * , 'Encoding labels for the plot.'
C
        IF (ULIN.GT.0.) THEN
          IF (UMAX-UMIN.LE.20.*ULIN) THEN
            UEPS=1.E-3*(UMAX-UMIN)
            IST1=INT((UMIN-UEPS)/ULIN+.5+SIGN(.5,UMIN-UEPS))
            IST2=INT((UMAX+UEPS)/ULIN-.5+SIGN(.5,UMAX+UEPS))
            WRITE (UNLB,UFMT) (ULIN*REAL(I),I=IST1,IST2)
          ELSE
            PRINT *, 'Omitting too-dense numeric labels in U direction.'
            UNLB=' '
          END IF
        END IF
C
        IF (VLIN.GT.0.) THEN
          IF (VMAX-VMIN.LE.20.*VLIN) THEN
            VEPS=1.E-3*(VMAX-VMIN)
            JST1=INT((VMIN-VEPS)/VLIN+.5+SIGN(.5,VMIN-VEPS))
            JST2=INT((VMAX+VEPS)/VLIN-.5+SIGN(.5,VMAX+VEPS))
            WRITE (VNLB,VFMT) (VLIN*REAL(J),J=JST1,JST2)
          ELSE
            PRINT *, 'Omitting too-dense numeric labels in V direction.'
            VNLB=' '
          END IF
        END IF
C
        IF (WLIN.GT.0.) THEN
          IF (WMAX-WMIN.LE.20.*WLIN) THEN
            WEPS=1.E-3*(WMAX-WMIN)
            KST1=INT((WMIN-WEPS)/WLIN+.5+SIGN(.5,WMIN-WEPS))
            KST2=INT((WMAX+WEPS)/WLIN-.5+SIGN(.5,WMAX+WEPS))
            WRITE (WNLB,WFMT) (WLIN*REAL(K),K=KST1,KST2)
          ELSE
            PRINT *, 'Omitting too-dense numeric labels in W direction.'
            WNLB=' '
          END IF
        END IF
C
C Set the overall character multiplier for TDPACK.
C
        CALL TDSETR ('CS1',RCS1)
C
C Define rendering styles 1, 2, and 3.
C
        PRINT * , 'Defining rendering styles.'
C
        IF (IFCB.LT.0.OR.IFCB.GT.7.OR.SHDE.EQ.0.) THEN
          IFC1=IFCB
          IFC2=IFCB
        ELSE
          IFC1=101+IFCB*16
          IFC2=116+IFCB*16
        END IF
C
        IF (IFCT.LT.0.OR.IFCT.GT.7.OR.SHDE.EQ.0.) THEN
          IFC3=IFCT
          IFC4=IFCT
        ELSE
          IFC3=101+IFCT*16
          IFC4=116+IFCT*16
        END IF
C
        IF (USIN.NE.0..AND.UMAX-UMIN.LE.50.*USIN) THEN
          USIM=1.
        ELSE
          IF (USIN.NE.0.)
     +    PRINT * , 'Omitting too-dense slice lines in U direction.'
          USIM=0.
        END IF
C
        IF (VSIN.NE.0..AND.VMAX-VMIN.LE.50.*VSIN) THEN
          VSIM=1.
        ELSE
          IF (VSIN.NE.0.)
     +    PRINT * , 'Omitting too-dense slice lines in V direction.'
          VSIM=0.
        END IF
C
        WSIM=1.
C
        IF (WSIN.NE.0..AND.WMAX-WMIN.LE.50.*WSIN) THEN
          WSIM=1.
        ELSE
          IF (WSIN.NE.0.)
     +    PRINT * , 'Omitting too-dense slice lines in W direction.'
          WSIM=0.
        END IF
C
        CALL TDSTRS (1,IFC1,IFC2,IFC3,IFC4,ILCB,ILCT,IDTR,
     +                      USIM*USIN,VSIM*VSIN,WSIM*WSIN)
        CALL TDSTRS (2, 109, 116,IFC3,IFC4,  -1,ILCT,IDTR,
     +                      USIM*USIN,VSIM*VSIN,WSIM*WSIN)
        CALL TDSTRS (3,IFC1,IFC2, 109, 116,ILCB,  -1,IDTR,
     +                      USIM*USIN,VSIM*VSIN,WSIM*WSIN)
        CALL TDSTRS (4, 109, 116, 213, 228,  -1,  -1,   0,
     +                             0.,       0.,       0.)
C
C Check for too-dense grid lines.
C
        PRINT * , 'Checking grid parameters.'
C
        IF (UGIN.GT.0.) THEN
          IF (UMAX-UMIN.LE.50.*UGIN) THEN
            UGIM=1.
          ELSE
            PRINT * , 'Omitting too-dense grid lines in U direction.'
            UGIM=0.
          END IF
        END IF
C
        IF (VGIN.GT.0.) THEN
          IF (VMAX-VMIN.LE.50.*VGIN) THEN
            VGIM=1.
          ELSE
            PRINT * , 'Omitting too-dense grid lines in V direction.'
            VGIM=0.
          END IF
        END IF
C
        IF (WGIN.GT.0.) THEN
          IF (WMAX-WMIN.LE.50.*WGIN) THEN
            WGIM=1.
          ELSE
            PRINT * , 'Omitting too-dense grid lines in W direction.'
            WGIM=0.
          END IF
        END IF
C
C Determine the positions of the eye and the point looked at.
C
        PRINT * , 'Determining eye position.'
C
        DOFB=SQRT((UMAX-UMIN)**2+(VMAX-VMIN)**2+(WMAX-WMIN)**2)
        DEYE=REYE*DOFB
        DAPT=RAPT*DOFB
C
        UEYE=UMID+DEYE*COS(DTOR*ANG1)*COS(DTOR*ANG2)
        VEYE=VMID+DEYE*SIN(DTOR*ANG1)*COS(DTOR*ANG2)
        WEYE=WMID+DEYE*SIN(DTOR*ANG2)
C
        UAPT=UMID-DAPT*COS(DTOR*ANG1)*COS(DTOR*ANG2)
        VAPT=VMID-DAPT*SIN(DTOR*ANG1)*COS(DTOR*ANG2)
        WAPT=WMID-DAPT*SIN(DTOR*ANG2)
C
C Compute some derivative quantities that are needed to put a label on
C the plot.
C
        UDNL=UMID-UEYE
        VDNL=VMID-VEYE
        WDNL=WMID-WEYE
C
        IF (IHND.EQ.0) THEN
          UDNX=VMID-VEYE
          VDNX=UEYE-UMID
          WDNX=0.
          UDNY=VDNX*WDNL-VDNL*WDNX
          VDNY=WDNX*UDNL-WDNL*UDNX
          WDNY=UDNX*VDNL-UDNL*VDNX
        ELSE
          UDNX=VEYE-VMID
          VDNX=UMID-UEYE
          WDNX=0.
          UDNY=VDNL*WDNX-VDNX*WDNL
          VDNY=WDNL*UDNX-WDNX*UDNL
          WDNY=UDNL*VDNX-UDNX*VDNL
        END IF
C
        DNOM=SQRT(UDNX*UDNX+VDNX*VDNX+WDNX*WDNX)
        UDCX=UDNX/DNOM
        VDCX=VDNX/DNOM
        WDCX=WDNX/DNOM
C
        DNOM=SQRT(UDNY*UDNY+VDNY*VDNY+WDNY*WDNY)
        UDCY=UDNY/DNOM
        VDCY=VDNY/DNOM
        WDCY=WDNY/DNOM
C
        CALL TDGETR ('CS1',CSM1)
        CSM2=CSM1*MIN(UMAX-UMIN,VMAX-VMIN,WMAX-WMIN)
C
C Initialize TDPACK.
C
        PRINT * , 'Initializing TDPACK.'
C
        CALL TDINIT (UEYE,VEYE,WEYE,UAPT,VAPT,WAPT,UAPT,VAPT,WAPT+1.,0.)
        CALL GETSET (XVPL,XVPR,YVPB,YVPT,XWDL,XWDR,YWDB,YWDT,LNLG)
C
C If the X workstation is not to be updated, jump to get user input.
C
        IF (IUXW.EQ.0) GO TO 112
C
C Clear the X workstation.
C
        CALL GCLRWK (1,1)
C
C Draw the plot label, if any.  The plot label is positioned in a
C rectangle that is perpendicular to the line of sight and passes
C through the midpoint of the data box; some algebra is required to
C come up with the direction cosines of the sides of that rectangle.
C We can then call TDPARA to define it as the current reference
C parallelogram and TDPLCH to draw a label relative to the rectangle.
C
        IF (PILB.NE.' ') THEN
          PRINT * , 'Drawing plot label.'
          CALL TDPARA (UMID-.5*DOFB*(UDCX+UDCY),
     +                 VMID-.5*DOFB*(VDCX+VDCY),
     +                 WMID-.5*DOFB*(WDCX+WDCY),
     +                 UDCX,VDCX,WDCX,UDCY,VDCY,WDCY)
          CALL TDPLCH ((.5+PILX)*DOFB,(.5+PILY)*DOFB,
     +                 PILB(1:LNBPCS(PILB)),.04*CSM2,0.,0.)
        END IF
C
C Put some simple labels on the axes of the plot.  If the last argument
C in the call to TDLBLS is a zero, all six outer edges of the projection
C of the box are labelled; if the argument is non-zero, then only three
C edges are labelled (which three depends on the sign of the argument),
C which gives one set of U axis labels, one set of V axis labels, and
C one set of W axis labels.  If the label fading flag is set and the
C viewing angle is close to a multiple of 90 degrees, the polyline and
C fill area color indices are set so as to do the fade (useful mostly
C when a movie is being made).
C
        IF (IGDF.NE.0) THEN
C
          PRINT * , 'Drawing axis labels.'
C
          IF (IFDE.NE.0) THEN
            ANGM=MOD(ANG1,90.)
            IF (ANGM.LT. 0.) ANGM=90.+ANGM
            IF (ANGM.GT.45.) ANGM=90.-ANGM
            MANG=101+INT(ANGM)
          ELSE
            MANG=116
          END IF
C
          IF (MANG.GT.101) THEN
            IF (MANG.GE.102.AND.MANG.LE.115) THEN
              CALL SFLUSH
              IF (IBOW.EQ.0) THEN
                CALL GSPLCI (217-MANG)
                CALL GSFACI (217-MANG)
              ELSE
                CALL GSPLCI (MANG)
                CALL GSFACI (MANG)
              END IF
            END IF
            CALL TDLBLS (UMIN,VMIN,WMIN,UMAX,VMAX,WMAX,
     +                   UNLB,VNLB,WNLB,UILB,VILB,WILB,1)
            IF (MANG.GE.102.AND.MANG.LE.115) THEN
              CALL SFLUSH
              CALL GSPLCI (1)
              CALL GSFACI (1)
            END IF
          END IF
C
C Draw the sides of the box that could be hidden.
C
          PRINT * , 'Drawing far sides of box.'
C
          CALL TDGRDS (UMIN,VMIN,WMIN,UMAX,VMAX,WMAX,
     +                 UGIM*UGIN,VGIM*VGIN,WGIM*WGIN,IGDF,1)
C
        END IF
C
C Draw all the triangles.
C
        PRINT * , 'Rendering surface.'
C
        CALL TDOTRI (RTRI,MTRI,NTRI,RTWK,ITWK,IORD)
        CALL TDDTRI (RTRI,MTRI,NTRI,ITWK)
C
C Draw the sides of the box that cannot be hidden.
C
        IF (IGDF.NE.0) THEN
C
          PRINT * , 'Drawing near sides of box.'
C
          CALL TDGRDS (UMIN,VMIN,WMIN,UMAX,VMAX,WMAX,
     +                 UGIM*UGIN,VGIM*VGIN,WGIM*WGIN,IGDF,0)
C
        END IF
C
C Flush buffers and update the X workstation.
C
        CALL SFLUSH
        CALL GUWK (1,0)
C
C Let the user enter a command.
C
  112   PRINT * , ' '
        PRINT * , 'Enter a command (H for help, Q for quit):'
C
        READ  '(A1)', COMD
C
        IF      (COMD.EQ.'D'.OR.COMD.EQ.'d') THEN
C
          IF (ISAS.GE. 99) ISAS=ISAS-100
          ISAS=MAX(-1,MIN(5,ISAS))
          PRINT * , ' '
          PRINT * , 'Change the type of data displayed (Y or N)?'
          READ  '(A1)', COMD
          IF (COMD.EQ.'Y'.OR.COMD.EQ.'y') THEN
            PRINT * , ' '
            PRINT * , 'The current surface type selector is ',ISAS
            PRINT * , ' '
            PRINT * , 'Enter -1 for a random surface, 0 for a random'
            PRINT * , 'isosurface, 1-5 for various canned isosurfaces,'
            PRINT * , '6 or greater to read a data file.'
            CALL TDRDIN (ITMP,ISAS)
            IF (ITMP.LE.5) THEN
              ISAS=MAX(-1,ITMP)
            ELSE
              PRINT * , ' '
              PRINT * , 'Enter the name of a data file:'
              READ '(A128)' , FLNM
              OPEN (13,FILE=FLNM(1:LNBPCS(FLNM)),STATUS='OLD',
     +                                         FORM='FORMATTED',ERR=113)
              PRINT * , 'File opened successfully.'
              READ (13,'(2I4)',ERR=113) IDIM,JDIM
              IF (IDIM.LT.1.OR.IDIM.GT.IMAX.OR.
     +            JDIM.LT.1.OR.JDIM.GT.JMAX) THEN
                PRINT * , 'Dimensions read are out of range!'
                GO TO 112
              END IF
              READ (13,'(6E12.0)',ERR=113,END=114)
     +                                     UMIN,UMAX,VMIN,VMAX,WMIN,WMAX
              UMID=.5*(UMIN+UMAX)
              VMID=.5*(VMIN+VMAX)
              WMID=.5*(WMIN+WMAX)
              READ (13,'(6E12.0)',ERR=113,END=114) (U(I),I=1,IDIM)
              READ (13,'(6E12.0)',ERR=113,END=114) (V(I),I=1,JDIM)
              READ (13,'(6E12.0)',ERR=113,END=114)
     +                                ((S(I,J),I=1,IDIM),J=1,JDIM)
              CLOSE (13)
              ISAS=ISAS+100
              GO TO 110
  113         PRINT * , ' '
              PRINT * , 'Error trying to read data file!'
              GO TO 112
  114         PRINT * , ' '
              PRINT * , 'Premature EOF on data file!'
              GO TO 112
            END IF
          END IF
          PRINT * , ' '
          PRINT * , 'Change the data dimensions (Y or N)?'
          READ  '(A1)', COMD
          IF (COMD.EQ.'Y'.OR.COMD.EQ.'y') THEN
            PRINT * , ' '
            PRINT * , 'Current dimensions are as follows:'
            PRINT * , '  1st data dimension: ',IDIM
            PRINT * , '  2nd data dimension: ',JDIM
            PRINT * , '  3rd data dimension: ',KDIM
            PRINT * , ' '
            PRINT * , 'Enter the new 1st data dimension:'
            CALL TDRDIN (IDIM,IDIM)
            IDIM=MAX(5,MIN(IMAX,IDIM))
            PRINT * , 'Enter the new 2nd data dimension:'
            CALL TDRDIN (JDIM,JDIM)
            JDIM=MAX(5,MIN(JMAX,JDIM))
            PRINT * , 'Enter the new 3rd data dimension:'
            CALL TDRDIN (KDIM,KDIM)
            JDIM=MAX(5,MIN(KMAX,KDIM))
          END IF
          IF (ISAS.LE.0) THEN
            PRINT * , ' '
            PRINT * , 'Use old centers for exponentials (Y or N)?'
            READ  '(A1)', COMD
            IF (COMD.EQ.'Y'.OR.COMD.EQ.'y') THEN
              ICEN=0
            ELSE
              ICEN=1
            END IF
          END IF
          GO TO 103
C
        ELSE IF (COMD.EQ.'G'.OR.COMD.EQ.'g') THEN
C
          PRINT * , ' '
          PRINT * , 'Current value of grid type selector: ',IGDF
          PRINT * , ' '
          PRINT * , 'Enter new value of grid type selector (a 2-digit'
          PRINT * , 'integer).  The first digit applies to the near'
          PRINT * , 'sides of the box and the second to the far sides'
          PRINT * , 'of the box.  Each digit is a 0 (nothing), a 1'
          PRINT * , '(perimeter only), a 2 (perimeter plus ticks),'
          PRINT * , 'or a 3 (perimeter plus grid lines):'
          CALL TDRDIN (IGDF,IGDF)
          IGDF=10*MAX(0,MIN(3,IGDF/10))+MAX(0,MIN(3,MOD(IGDF,10)))
          IF (IGDF.NE.0) THEN
            PRINT * , ' '
            PRINT * , 'Change spacing of grid lines (Y or N)?'
            READ  '(A1)', COMD
            IF (COMD.EQ.'Y'.OR.COMD.EQ.'y') THEN
              PRINT * , ' '
              PRINT * , 'Current spacing of grid lines:'
              PRINT * , '  Grid-line spacing in U: ',UGIN
              PRINT * , '  Grid-line spacing in V: ',VGIN
              PRINT * , '  Grid-line spacing in W: ',WGIN
              PRINT * , ' '
              PRINT * , 'Enter new value for grid-line spacing in U:'
              CALL TDRDRN (UGIN,UGIN)
              UGIN=MAX(0.,UGIN)
              PRINT * , 'Enter new value for grid-line spacing in V:'
              CALL TDRDRN (VGIN,VGIN)
              VGIN=MAX(0.,VGIN)
              PRINT * , 'Enter new value for grid-line spacing in W:'
              CALL TDRDRN (WGIN,WGIN)
              WGIN=MAX(0.,WGIN)
            END IF
            PRINT * , ' '
            PRINT * , 'Change spacing of numeric labels (Y or N)?'
            READ  '(A1)', COMD
            IF (COMD.EQ.'Y'.OR.COMD.EQ.'y') THEN
            PRINT * , ' '
              PRINT * , 'Current spacing of numeric labels:'
              PRINT * , '  Numeric-label spacing in U: ',ULIN
              PRINT * , '  Numeric-label spacing in V: ',VLIN
              PRINT * , '  Numeric-label spacing in W: ',WLIN
              PRINT * , ' '
              PRINT * , 'Enter new value for label spacing in U:'
              CALL TDRDRN (ULIN,ULIN)
              ULIN=MAX(0.,ULIN)
              PRINT * , 'Enter new value for label spacing in V:'
              CALL TDRDRN (VLIN,VLIN)
              VLIN=MAX(0.,VLIN)
              PRINT * , 'Enter new value for label spacing in W:'
              CALL TDRDRN (WLIN,WLIN)
              WLIN=MAX(0.,WLIN)
            END IF
            PRINT * , ' '
            PRINT * , 'Change formats for numeric labels (Y or N)?'
            READ  '(A1)', COMD
            IF (COMD.EQ.'Y'.OR.COMD.EQ.'y') THEN
            PRINT * , ' '
              PRINT * , 'Current formats for numeric labels:'
              PRINT * , '  Format for U labels: ',UFMT(1:LNBPCS(UFMT))
              PRINT * , '  Format for V labels: ',VFMT(1:LNBPCS(VFMT))
              PRINT * , '  Format for W labels: ',WFMT(1:LNBPCS(WFMT))
              PRINT * , ' '
              PRINT * , 'Formats must not be more than 16 characters'
              PRINT * , 'long and must not specify a line more than'
              PRINT * , '128 characters long.  Enter them carefully.'
              PRINT * , ' '
              PRINT * , 'Enter new U format:'
              READ '(A16)' , UFMT
              PRINT * , 'Enter new V format:'
              READ '(A16)' , VFMT
              PRINT * , 'Enter new W format:'
              READ '(A16)' , WFMT
            END IF
            PRINT * , ' '
            PRINT * , 'Change informational labels (Y or N)?'
            READ  '(A1)', COMD
            IF (COMD.EQ.'Y'.OR.COMD.EQ.'y') THEN
              PRINT * , ' '
              PRINT * , 'Current U axis label:',UILB(1:LNBPCS(UILB))
              PRINT * , 'Current V axis label:',VILB(1:LNBPCS(VILB))
              PRINT * , 'Current W axis label:',WILB(1:LNBPCS(WILB))
              PRINT * , 'Current plot label:  ',PILB(1:LNBPCS(PILB))
              PRINT * , ' '
              PRINT * , 'Enter new U axis label:'
              READ '(A128)' , UILB
              PRINT * , 'Enter new V axis label:'
              READ '(A128)' , VILB
              PRINT * , 'Enter new W axis label:'
              READ '(A128)' , WILB
              PRINT * , 'Enter new plot label:'
              READ '(A128)' , PILB
            END IF
            PRINT * , ' '
            PRINT * , 'Change informational-label positions (Y or N)?'
            READ  '(A1)', COMD
            IF (COMD.EQ.'Y'.OR.COMD.EQ.'y') THEN
              PRINT * , ' '
              PRINT * , 'Current plot label position (horizontal):',PILX
              PRINT * , 'Current plot label position (vertical):  ',PILY
              PRINT * , ' '
              PRINT * , 'Enter new plot label position (horizontal):'
              CALL TDRDRN (PILX,PILX)
              PILX=MAX(-1.,MIN(1.,PILX))
              PRINT * , 'Enter new plot label position (vertical):'
              CALL TDRDRN (PILY,PILY)
              PILY=MAX(-1.,MIN(1.,PILY))
            END IF
            PRINT * , ' '
            IF (IFDE.NE.0) THEN
              PRINT * , 'Axis labels are currently faded for azimuth ang
     +les near 90 degrees.'
            ELSE
              PRINT * , 'Axis labels are not currently faded for azimuth
     + angles near 90 degrees.'
            END IF
            PRINT * , ' '
            PRINT * , 'Fade labels for azimuth angles near 90 degrees (Y
     + or N)?'
            READ  '(A1)', COMD
            IF (COMD.EQ.'Y'.OR.COMD.EQ.'y') THEN
              IFDE=1
            ELSE
              IFDE=0
            END IF
          END IF
          GO TO 111
C
        ELSE IF (COMD.EQ.'H'.OR.COMD.EQ.'h') THEN
C
          PRINT * , ' '
          PRINT * , 'Possible commands are as follows:'
          PRINT * , '  D  =>  Data change'
          PRINT * , '  G  =>  Grid change'
          PRINT * , '  H  =>  Help'
          PRINT * , '  M  =>  Miscellaneous changes'
          PRINT * , '  Q  =>  Quit'
          PRINT * , '  R  =>  Refresh/rendering change'
          PRINT * , '  S  =>  Save plot (NCGM and/or PostScript)'
          PRINT * , '  V  =>  View change'
          GO TO 112
C
        ELSE IF (COMD.EQ.'M'.OR.COMD.EQ.'m') THEN
C
          PRINT * , ' '
          IF (IBOW.EQ.0) THEN
            PRINT * , 'The color scheme flag is now 0 (white on black).'
          ELSE
            PRINT * , 'The color scheme flag is now 1 (black on white).'
          END IF
          PRINT * , ' '
          PRINT * , 'Enter new color scheme flag (0 => WOB, 1 => BOW):'
          CALL TDRDIN (ITMP,IBOW)
          ITMP=MAX(0,MIN(1,ITMP))
          IF (ITMP.NE.IBOW) THEN
            IBOW=ITMP
            CALL GCLRWK (1,1)
            DO 115 IWID=1,3
              IF (IWID.NE.1) CALL GACWK (IWID)
              IF (IBOW.EQ.0) THEN
                CALL GSCR (IWID,0,0.,0.,0.) ! black
                CALL GSCR (IWID,1,1.,1.,1.) ! white
              ELSE
                CALL GSCR (IWID,0,1.,1.,1.) ! white
                CALL GSCR (IWID,1,0.,0.,0.) ! black
              END IF
              IF (IWID.NE.1) CALL GDAWK (IWID)
  115       CONTINUE
            SSHE=SHDE
            SSHR=SHDR
          END IF
C
          PRINT * , ' '
          IF (IUXW.EQ.0) THEN
            PRINT * , 'The X-workstation-update flag is now 0 (off).'
          ELSE
            PRINT * , 'The X-workstation-update flag is now 1 (on).'
          END IF
          PRINT * , ' '
          PRINT * , 'Enter new X-workstation-update flag (0 => off, 1 =>
     + on):'
          CALL TDRDIN (IUXW,IUXW)
          IUXW=MAX(0,MIN(1,IUXW))
C
          PRINT * , ' '
          IF (IHND.EQ.0) THEN
            PRINT * , 'The handedness flag is now 0 (right-handed).'
          ELSE
            PRINT * , 'The handedness flag is now 1 (left-handed).'
          END IF
          PRINT * , ' '
          PRINT * , 'Enter new handedness flag (0 => right, 1 => left):'
          CALL TDRDIN (IHND,IHND)
          IHND=MAX(0,MIN(1,IHND))
          CALL TDSETI ('HND',IHND)
C
          PRINT * , ' '
          PRINT * , 'Current value of stereo-view flag: ',STVF
          PRINT * , ' '
          PRINT * , 'Enter a new value of the stereo-view flag, which'
          PRINT * , 'is set non-zero to make the S command save stereo'
          PRINT * , 'views.  Use a negative value to get each stereo'
          PRINT * , 'pair on a single frame or a positive value to use'
          PRINT * , 'separate frames.  The absolute value is the angle,'
          PRINT * , 'in degrees, between the lines of sight of the two'
          PRINT * , 'views.  (A value of one or two degrees seems to'
          PRINT * , 'work best):'
          CALL TDRDRN (STVF,STVF)
          STVF=MAX(-10.,MIN(10.,STVF))
C
          PRINT * , ' '
          PRINT * , 'Current character size multiplier: ',RCS1
          PRINT * , ' '
          PRINT * , 'Enter new value for character size multiplier:'
          CALL TDRDRN (RCS1,RCS1)
          RCS1=MAX(.1,MIN(10.,RCS1))
C
          PRINT * , ' '
          PRINT * , 'Multiplier for random term in dummy data: ',RNDM
          PRINT * , ' '
          PRINT * , 'Enter new multiplier for random term:'
          CALL TDRDRN (RNDM,RNDM)
          RNDM=MAX(0.,RNDM)
C
          PRINT * , ' '
          PRINT * , 'Current triangle-ordering flag: ',IORD
          PRINT * , ' '
          PRINT * , 'Enter new triangle-ordering flag (0 => sort on  '
          PRINT * , 'distances to triangle centers, -1 => sort on    '
          PRINT * , 'distances to triangle farpoints, +1 => sort on  '
          PRINT * , 'distances to triangle farpoints and then run an '
          PRINT * , 'algorithm to fix ordering errors):'
          CALL TDRDIN (IORD,IORD)
          IORD=MAX(-1,MIN(+1,IORD))
C
          PRINT * , ' '
          PRINT * , 'Current flag for the grid of 3D marks: ',IMRK
          IF (IMRK.EQ.0) THEN
            PRINT * , 'No grid of 3D marks is to be generated.'
          ELSE
            PRINT * , 'Grid of 3D marks is to be generated, each ',
     +                      IMKT(IABS(IMRK))(1:LNBPCS(IMKT(IABS(IMRK))))
            IF (IMRK.LT.0) THEN
              PRINT * , 'Marks are not to be clipped by the data box.'
            ELSE
              PRINT * , 'Marks are to be clipped by the data box.'
            END IF
          END IF
          IMKO=IMRK
          SMKO=SMRK
          UMKO=UMKS
          VMKO=VMKS
          WMKO=WMKS
          PRINT * , ' '
          PRINT * , 'Enter new value of flag for the grid of 3D marks'
          PRINT * , '(0 => no marks, 1 => tetrahedron, 2 => octahedron,'
          PRINT * , '3 = > cube, 4 => icosahedron, 5 => sphere; use a'
          PRINT * , 'negated value to omit clipping by data box.'
          CALL TDRDIN (IMRK,IMRK)
          IMRK=MAX(-5,MIN(5,IMRK))
          IF (IMRK.NE.0) THEN
            PRINT * , ' '
            PRINT * , 'Current value of mark size parameter:   ',SMRK
            PRINT * , ' '
            PRINT * , 'Enter new mark size (as a fraction of the'
            PRINT * , 'smallest side of the data box):'
            CALL TDRDRN (SMRK,SMRK)
            SMRK=MAX(.001,MIN(1.,SMRK))
            PRINT * , ' '
            PRINT * , 'Current mark spacing parameters:',UMKS,VMKS,WMKS
            PRINT * , ' '
            PRINT * , 'Enter new mark spacing in U:'
            CALL TDRDRN (UMKS,UMKS)
            UMKS=MAX(0.,UMKS)
            PRINT * , 'Enter new mark spacing in V:'
            CALL TDRDRN (VMKS,VMKS)
            VMKS=MAX(0.,VMKS)
            PRINT * , 'Enter new mark spacing in W:'
            CALL TDRDRN (WMKS,WMKS)
            WMKS=MAX(0.,WMKS)
          END IF
C
          PRINT * , ' '
          PRINT * , 'Current isosurface cut-off value 1: ',FIS1
          PRINT * , 'Current isosurface cut-off value 2: ',FIS2
          PRINT * , ' '
          PRINT * , 'Isosurface cut-off values must be between 0'
          PRINT * , 'and 200 (the initial values were 95 and 95).'
          PRINT * , 'If the two values are equal, one isosurface'
          PRINT * , 'is specified; if the two values are different,'
          PRINT * , 'two isosurfaces are specified.'
          PRINT * , ' '
          PRINT * , 'Enter new isosurface cut-off value 1:'
          OFS1=FIS1
          CALL TDRDRN (FIS1,FIS1)
          FIS1=MAX(0.,MIN(200.,FIS1))
          PRINT * , 'Enter new isosurface cut-off value 2:'
          OFS2=FIS2
          CALL TDRDRN (FIS2,FIS2)
          FIS2=MAX(0.,MIN(200.,FIS2))
          IF (FIS1.GT.FIS2) THEN
            FIS0=FIS1
            FIS1=FIS2
            FIS2=FIS0
          END IF
C
          IF (IMRK.NE.IMKO.OR.SMRK.NE.SMKO.OR.UMKS.NE.UMKO.OR.
     +        UMKS.NE.UMKO.OR.VMKS.NE.VMKO.OR.WMKS.NE.WMKO.OR.
     +        FIS1.NE.OFS1.OR.FIS2.NE.OFS2) THEN
            GO TO 110
          ELSE
            GO TO 111
          END IF
C
        ELSE IF (COMD.EQ.'Q'.OR.COMD.EQ.'q') THEN
C
          CALL GDAWK (1)
          CALL GCLWK (1)
          CALL GCLWK (2)
          IF (NFPS.NE.0) THEN
            CALL GACWK (3)
            CALL GCLRWK (3,1)
            CALL GDAWK (3)
          END IF
          CALL GCLWK (3)
          CALL GCLKS
          STOP
C
        ELSE IF (COMD.EQ.'R'.OR.COMD.EQ.'r') THEN
C
          PRINT * , ' '
          PRINT * , 'Redrawing the picture.'
          PRINT * , ' '
          PRINT * , 'Change the shading parameters (Y or N)?'
          READ  '(A1)', COMD
          IF (COMD.EQ.'Y'.OR.COMD.EQ.'y') THEN
            PRINT * , ' '
            PRINT * , 'Current shading parameter values are as follows:'
            PRINT * , '  Shading parameter 1:   ',SHDE
            PRINT * , '  Shading parameter 2:   ',SHDR
            PRINT * , '  TDPACK parameter SHD:  ',ISHD
            PRINT * , '  Light source azimuth:  ',ANG3
            PRINT * , '  Light source elevation:',ANG4
            PRINT * , ' '
            PRINT * , 'Note that if both angles are zero, the light'
            PRINT * , 'source is at the position of the viewer.'
            PRINT * , ' '
            PRINT * , 'Enter new shading parameter 1 (0 => shading off,'
            PRINT * , 'near 0 => brighter colors, near 1 => pastels).'
            CALL TDRDRN (SHDE,SHDE)
            SHDE=MAX(0.,MIN(1.,SHDE))
            PRINT * , 'Enter new shading parameter 2 (near 0 => small'
            PRINT * , 'range of shades, near 1 => full range).'
            CALL TDRDRN (SHDR,SHDR)
            SHDR=MAX(0.,MIN(1.,SHDR))
            PRINT * , 'Enter new value of TDPACK internal parameter SHD'
            PRINT * , '(0 => shading over a 90-degree range of triangle'
            PRINT * , 'orientations, 1 => over a 180-degree range):'
            CALL TDRDIN (ISHD,ISHD)
            ISHD=MAX(0,MIN(1,ISHD))
            CALL TDSETI ('SHD',ISHD)
            PRINT * , 'Enter light source angle 1 (azimuth, in degrees):
     +'
            CALL TDRDRN (ANG3,ANG3)
            ANG3=SIGN(MOD(ABS(ANG3),360.),ANG3)
            PRINT * , 'Enter light source angle 2 (elevation, in degrees
     +):'
            CALL TDRDRN (ANG4,ANG4)
            ANG4=MAX(-89.999,MIN(+89.999,ANG4))
            IF (ANG3.NE.0..OR.ANG4.NE.0.) THEN
              UCLS=UMID+100.*DOFB*COS(DTOR*ANG3)*COS(DTOR*ANG4)
              VCLS=VMID+100.*DOFB*SIN(DTOR*ANG3)*COS(DTOR*ANG4)
              WCLS=WMID+100.*DOFB*SIN(DTOR*ANG4)
            ELSE
              UCLS=0.
              VCLS=0.
              WCLS=0.
            END IF
            CALL TDSETR ('LSU',UCLS)
            CALL TDSETR ('LSV',VCLS)
            CALL TDSETR ('LSW',WCLS)
            IF (SHDE.NE.0..AND.(SHDE.NE.SSHE.OR.SHDR.NE.SSHR)) THEN
              CALL GCLRWK (1,1)
              DO 117 IWID=1,3
                IF (IWID.NE.1) CALL GACWK (IWID)
                CALL TDCLRS (IWID,IBOW,SHDE,SHDR,101,116,0)
                IF (IWID.NE.1) CALL GDAWK (IWID)
  117         CONTINUE
              SSHE=SHDE
              SSHR=SHDR
            END IF
          END IF
          PRINT * , ' '
          PRINT * , 'Change the rendering parameters (Y or N)?'
          READ  '(A1)', COMD
          IF (COMD.EQ.'Y'.OR.COMD.EQ.'y') THEN
            PRINT * , ' '
            PRINT * , 'Usable color indices are as follows:'
            PRINT * , '  -1 => if line color, line not drawn; if fill co
     +lor, no fill done'
            PRINT * , '   0 => background or 101-116 => shades from all-
     +white to all-black'
            PRINT * , '   1 => foreground or 117-132 => shades of gray (
     +light to dark)'
            PRINT * , '   2 => red        or 133-148 => shades of red (l
     +ight to dark)'
            PRINT * ,  '   3 => green      or 149-164 => shades of green
     + (light to dark)'
            PRINT * , '   4 => blue       or 165-180 => shades of blue (
     +light to dark)'
            PRINT * , '   5 => cyan       or 181-196 => shades of cyan (
     +light to dark)'
            PRINT * , '   6 => magenta    or 197-212 => shades of magent
     +a (light to dark)'
            PRINT * , '   7 => yellow     or 213-228 => shades of yellow
     + (light to dark)'
            PRINT * , '  Colors with indices 117 through 228 are affecte
     +d by the values of'
            PRINT * , '  shading flags.  Using a fill color index 0 thro
     +ugh 7 when shading'
            PRINT * , '  is turned on implies using shades of that color
     + (for example, use'
            PRINT * , '  of the value 4 implies the use of color indices
     + 165 through 180).'
            PRINT * , ' '
            PRINT * , 'Current rendering values are as follows:'
            PRINT * , '  Bottom/top fill colors:      ',IFCB,IFCT
            PRINT * , '  Bottom/top line colors:      ',ILCB,ILCT
            PRINT * , '  Triangle-drawing flag:       ',IDTR
            PRINT * , '  U/V/W slice separations:     ',USIN,VSIN,WSIN
            PRINT * , ' '
            PRINT * , 'Enter new bottom fill color:'
            CALL TDRDIN (IFCB,IFCB)
            IFCB=MAX(-1,IFCB)
            PRINT * , 'Enter new top fill color:'
            CALL TDRDIN (IFCT,IFCT)
            IFCT=MAX(-1,IFCT)
            PRINT * , 'Enter new bottom line color:'
            CALL TDRDIN (ILCB,ILCB)
            ILCB=MAX(-1,ILCB)
            PRINT * , 'Enter new top line color:'
            CALL TDRDIN (ILCT,ILCT)
            ILCT=MAX(-1,ILCT)
            PRINT * , 'Enter new triangle-drawing flag (0 or 1):'
            CALL TDRDIN (IDTR,IDTR)
            IDTR=MAX(0,MIN(1,IDTR))
            PRINT * , 'Enter new U slice separation (0 => none):'
            CALL TDRDRN (USIN,USIN)
            USIN=MAX(0.,USIN)
            PRINT * , 'Enter new V slice separation (0 => none):'
            CALL TDRDRN (VSIN,VSIN)
            VSIN=MAX(0.,VSIN)
            PRINT * , 'Enter new W slice separation (0 => none):'
            CALL TDRDRN (WSIN,WSIN)
            WSIN=MAX(0.,WSIN)
          END IF
          GO TO 111
C
        ELSE IF (COMD.EQ.'S'.OR.COMD.EQ.'s') THEN
C
          PRINT * , ' '
          IF (STVF.EQ.0.) THEN
            PRINT * , 'Saving just the current image.  (An M command'
            PRINT * , 'can be used to turn on saving of stereo pairs.)'
          ELSE
            PRINT * , 'Saving a stereo pair of the current image.'
          END IF
          PRINT * , ' '
          PRINT * , 'Do you want to save to an NCAR CGM file (Y or N)?'
          READ  '(A1)', COMD
          IF (COMD.NE.'Y'.AND.COMD.NE.'y') THEN
            ISNG=0
          ELSE
            ISNG=1
            CALL GACWK (2)
            NFNG=NFNG+1
            IF (NFNG.GT.1) CALL GCLRWK (2,1)
          END IF
          PRINT * , ' '
          PRINT * , 'Do you want to save to a PostScript file (Y or N)?'
          READ  '(A1)', COMD
          IF (COMD.NE.'Y'.AND.COMD.NE.'y') THEN
            ISPS=0
          ELSE
            ISPS=1
            CALL GACWK (3)
            NFPS=NFPS+1
            IF (NFPS.GT.1) CALL GCLRWK (3,1)
          END IF
          IF (ISNG.NE.0.OR.ISPS.NE.0) THEN
            CALL GDAWK(1)
            OFFS=-(DEYE+DAPT)*TAN(DTOR*ABS(STVF)/2.)
  118       PRINT * , ' '
            IF      (OFFS.LT.0.) THEN
              PRINT * , 'Initializing TDPACK for left-eye view.'
            ELSE IF (OFFS.GT.0.) THEN
              PRINT * , 'Initializing TDPACK for right-eye view.'
            ELSE
              PRINT * , 'Initializing TDPACK for single view.'
            END IF
            CALL TDINIT (UEYE,VEYE,WEYE,UAPT,VAPT,WAPT,
     +                                  UAPT,VAPT,WAPT+1.,OFFS)
            CALL GETSET (XVPL,XVPR,YVPB,YVPT,XWDL,XWDR,YWDB,YWDT,LNLG)
            IF (OFFS.NE.0..AND.STVF.LT.0.) THEN
              IF (OFFS.LT.0.) THEN
                CALL SET  (1.-WOSW,1.,.5-.5*WOSW,.5+.5*WOSW,
     +                             XWDL,XWDR,YWDB,YWDT,LNLG)
              ELSE
                CALL SET  (  0., WOSW,.5-.5*WOSW,.5+.5*WOSW,
     +                             XWDL,XWDR,YWDB,YWDT,LNLG)
              END IF
            END IF
            IF (PILB.NE.' ') THEN
              PRINT * , 'Drawing plot label.'
              CALL TDPARA (UMID-.5*DOFB*(UDCX+UDCY),
     +                     VMID-.5*DOFB*(VDCX+VDCY),
     +                     WMID-.5*DOFB*(WDCX+WDCY),
     +                     UDCX,VDCX,WDCX,UDCY,VDCY,WDCY)
              CALL TDPLCH ((.5+PILX)*DOFB,(.5+PILY)*DOFB,
     +                     PILB(1:LNBPCS(PILB)),.04*CSM2,0.,0.)
            END IF
            IF (IGDF.NE.0) THEN
              PRINT * , 'Drawing axis labels.'
              IF (IFDE.NE.0) THEN
                ANGM=MOD(ANG1,90.)
                IF (ANGM.LT. 0.) ANGM=90.+ANGM
                IF (ANGM.GT.45.) ANGM=90.-ANGM
                MANG=101+INT(ANGM)
              ELSE
                MANG=116
              END IF
              IF (MANG.GT.101) THEN
                IF (MANG.GE.102.AND.MANG.LE.115) THEN
                  CALL SFLUSH
                  IF (IBOW.EQ.0) THEN
                    CALL GSPLCI (217-MANG)
                    CALL GSFACI (217-MANG)
                  ELSE
                    CALL GSPLCI (MANG)
                    CALL GSFACI (MANG)
                  END IF
                END IF
                CALL TDLBLS (UMIN,VMIN,WMIN,UMAX,VMAX,WMAX,
     +                       UNLB,VNLB,WNLB,UILB,VILB,WILB,1)
                IF (MANG.GE.102.AND.MANG.LE.115) THEN
                  CALL SFLUSH
                  CALL GSPLCI (1)
                  CALL GSFACI (1)
                END IF
              END IF
              PRINT * , 'Drawing far sides of box.'
              CALL TDGRDS (UMIN,VMIN,WMIN,UMAX,VMAX,WMAX,
     +                     UGIM*UGIN,VGIM*VGIN,WGIM*WGIN,IGDF,1)
            END IF
            PRINT * , 'Rendering surface.'
            CALL TDOTRI (RTRI,MTRI,NTRI,RTWK,ITWK,IORD)
            CALL TDDTRI (RTRI,MTRI,NTRI,ITWK)
            IF (IGDF.NE.0) THEN
              PRINT * , 'Drawing near sides of box.'
              CALL TDGRDS (UMIN,VMIN,WMIN,UMAX,VMAX,WMAX,
     +                     UGIM*UGIN,VGIM*VGIN,WGIM*WGIN,IGDF,0)
            END IF
            IF      (OFFS.LT.0.) THEN
              CALL PLCHHQ (.98*XWDL+.02*XWDR,.98*YWDB+.02*YWDT,'L',.006,
     +                                                            0.,0.)
            ELSE IF (OFFS.GT.0.) THEN
              CALL PLCHHQ (.02*XWDL+.98*XWDR,.98*YWDB+.02*YWDT,'R',.006,
     +                                                            0.,0.)
            END IF
            CALL SFLUSH
            IF (ISNG.NE.0) CALL GUWK (2,0)
            IF (ISPS.NE.0) CALL GUWK (3,0)
            IF (OFFS.LT.0.) THEN
              OFFS=-OFFS
              IF (STVF.GT.0.) THEN
                IF (ISNG.NE.0) CALL GCLRWK (2,1)
                IF (ISPS.NE.0) CALL GCLRWK (3,1)
              END IF
              GO TO 118
            END IF
            IF (ISNG.NE.0) CALL GDAWK (2)
            IF (ISPS.NE.0) CALL GDAWK (3)
            CALL GACWK(1)
          END IF
          GO TO 112
C
        ELSE IF (COMD.EQ.'V'.OR.COMD.EQ.'v') THEN
C
          PRINT * , ' '
          PRINT * , 'Changing view angle and distance.  Current'
          PRINT * , 'values are as follows:'
          PRINT * , '  Viewing angle 1 (azimuth):                 ',ANG1
          PRINT * , '  Viewing angle 2 (elevation):               ',ANG2
          PRINT * , '  Distance to eye (diagonal multiple):       ',REYE
          PRINT * , '  Distance to aim point (diagonal multiple): ',RAPT
          PRINT * , '  Field of view (in degrees):                ',FOVP
          PRINT * , ' '
          PRINT * , 'Enter viewing angle 1 (azimuth, in degrees):'
          CALL TDRDRN (ANG1,ANG1)
          ANG1=SIGN(MOD(ABS(ANG1),360.),ANG1)
          PRINT * , 'Enter viewing angle 2 (elevation, in degrees):'
          CALL TDRDRN (ANG2,ANG2)
          ANG2=MAX(-89.999,MIN(+89.999,ANG2))
          PRINT * , 'Enter distance to eye (diagonal multiple):'
          CALL TDRDRN (REYE,REYE)
          REYE=MAX(.51,MIN(51.,REYE))
          PRINT * , 'Enter distance to aim point (diagonal multiple):'
          CALL TDRDRN (RAPT,RAPT)
          RAPT=MAX(-.9*REYE,MIN(100.*REYE,RAPT))
          PRINT * , 'Enter desired field of view (in degrees):'
          CALL TDRDRN (FOVP,FOVP)
          FOVP=MAX(1.,MIN(179.,FOVP))
          CALL TDSETR ('FOV',FOVP)
          GO TO 111
C
        ELSE
C
          PRINT * , ' '
          PRINT * , 'Illegal command.  Try again.'
C
          GO TO 112
C
        END IF
C
      END


      SUBROUTINE GENDAT (DATA,IMAX,IDIM,JDIM,MLOW,MHGH,DLOW,DHGH,ICEN)
C
C This is a routine to generate test data for two-dimensional graphics
C routines.  Given an array "DATA", dimensioned "IMAX x 1", it fills
C the sub-array ((DATA(I,J),I=1,IDIM),J=1,JDIM) with a two-dimensional
C field of data having approximately "MLOW" lows and "MHGH" highs, a
C minimum value of exactly "DLOW" and a maximum value of exactly "DHGH".
C If ICEN is non-zero, new centers are computed for the exponentials,
C yielding a whole new field.
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(IMAX,*),CCNT(3,50)
C
        SAVE NLOW,NHGH,NCNT,CCNT
C
        DATA NCNT / 0 /
C
        IF (ICEN.NE.0.OR.NCNT.EQ.0) THEN
C
          NLOW=MAX0(1,MIN0(25,MLOW))
          NHGH=MAX0(1,MIN0(25,MHGH))
          NCNT=NLOW+NHGH
C
          DO 101 ICNT=1,NCNT
            CCNT(1,ICNT)=FRAN()
            CCNT(2,ICNT)=FRAN()
            IF (ICNT.LE.NLOW) THEN
              CCNT(3,ICNT)=-1.
            ELSE
              CCNT(3,ICNT)=+1.
            END IF
  101     CONTINUE
C
        END IF
C
        DMIN=+1.E36
        DMAX=-1.E36
C
        DO 104 I=1,IDIM
          XPOS=REAL(I-1)/REAL(IDIM-1)
          DO 103 J=1,JDIM
            YPOS=REAL(J-1)/REAL(JDIM-1)
            DATA(I,J)=0.
            DO 102 K=1,NCNT
              TEMP=-50.*((XPOS-CCNT(1,K))**2+(YPOS-CCNT(2,K))**2)
              IF (TEMP.GE.-20.) DATA(I,J)=DATA(I,J)+CCNT(3,K)*EXP(TEMP)
  102       CONTINUE
            DMIN=MIN(DMIN,DATA(I,J))
            DMAX=MAX(DMAX,DATA(I,J))
  103     CONTINUE
  104   CONTINUE
C
        DO 106 I=1,IDIM
          DO 105 J=1,JDIM
            DATA(I,J)=(DATA(I,J)-DMIN)/(DMAX-DMIN)*(DHGH-DLOW)+DLOW
  105     CONTINUE
  106   CONTINUE
C
        RETURN
C
      END


      SUBROUTINE ADDRAN (DATA,IMAX,IDIM,JDIM,DLOW,DHGH,RNDM)
C
        DIMENSION DATA(IMAX,*)
C
C This subroutine is called to add a random quantity to each element of
C the IDIM x JDIM array DATA, whose first FORTRAN dimension is IMAX.
C
C The magnitude of the random quantity is stated as a multiple (RNDM)
C of the range of the data.
C
        DMUL=RNDM*(DHGH-DLOW)
C
        DMIN=+1.E36
        DMAX=-1.E36
C
C Add random quantities to the data.
C
        DO 102 I=1,IDIM
          DO 101 J=1,JDIM
            DATA(I,J)=DATA(I,J)+DMUL*FRAN()
            DMIN=MIN(DMIN,DATA(I,J))
            DMAX=MAX(DMAX,DATA(I,J))
  101     CONTINUE
  102   CONTINUE
C
C Readjust the data to have the high and low values as before.
C
        DO 104 I=1,IDIM
          DO 103 J=1,JDIM
            DATA(I,J)=(DATA(I,J)-DMIN)/(DMAX-DMIN)*(DHGH-DLOW)+DLOW
  103     CONTINUE
  104   CONTINUE
C
C Done.
C
        RETURN
C
      END


      SUBROUTINE TDGNDT (DATA,IMAX,JMAX,IDIM,JDIM,KDIM,MLOW,MHGH,
     +                                            DLOW,DHGH,ICEN)
C
C This is a routine to generate test data for three-dimensional graphics
C routines.  Given an array "DATA", dimensioned "IMAX x JMAX x ...", it
C fills the sub-array (((DATA(I,J,K),I=1,MX),J=1,MY),K=1,MZ) with a
C three-dimensional field of data having approximately "MLOW" lows and
C "MHGH" highs, a minimum value of exactly "DLOW" and a maximum value
C of exactly "DHGH".  If ICEN is non-zero, new centers are computed for
C the exponentials, yielding a whole new field.
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(IMAX,JMAX,*),CCNT(4,50)
C
        SAVE NLOW,NHGH,NCNT,CCNT
C
        DATA NCNT / 0 /
C
        IF (ICEN.NE.0.OR.NCNT.EQ.0) THEN
C
          NLOW=MAX0(1,MIN0(25,MLOW))
          NHGH=MAX0(1,MIN0(25,MHGH))
          NCNT=NLOW+NHGH
C
          DO 101 ICNT=1,NCNT
            CCNT(1,ICNT)=FRAN()
            CCNT(2,ICNT)=FRAN()
            CCNT(3,ICNT)=FRAN()
            IF (ICNT.LE.NLOW) THEN
              CCNT(4,ICNT)=-1.
            ELSE
              CCNT(4,ICNT)=+1.
            END IF
  101     CONTINUE
C
        END IF
C
        DMIN=+1.E36
        DMAX=-1.E36
C
        DO 105 I=1,IDIM
          XPOS=REAL(I-1)/REAL(IDIM-1)
          DO 104 J=1,JDIM
            YPOS=REAL(J-1)/REAL(JDIM-1)
            DO 103 K=1,KDIM
              ZPOS=REAL(K-1)/REAL(KDIM-1)
              DATA(I,J,K)=0.
              DO 102 ICNT=1,NCNT
                TEMP=-50.*((XPOS-CCNT(1,ICNT))**2+
     +                     (YPOS-CCNT(2,ICNT))**2+
     +                     (ZPOS-CCNT(3,ICNT))**2)
                IF (TEMP.GE.-20.) DATA(I,J,K)=DATA(I,J,K)+
     +                                        CCNT(4,ICNT)*EXP(TEMP)
  102         CONTINUE
              DMIN=MIN(DMIN,DATA(I,J,K))
              DMAX=MAX(DMAX,DATA(I,J,K))
  103       CONTINUE
  104     CONTINUE
  105   CONTINUE
C
        DO 108 I=1,IDIM
          DO 107 J=1,JDIM
            DO 106 K=1,KDIM
              DATA(I,J,K)=DLOW+((DATA(I,J,K)-DMIN)/
     +                          (DMAX       -DMIN))*(DHGH-DLOW)
  106       CONTINUE
  107     CONTINUE
  108   CONTINUE
C
        RETURN
C
      END


      SUBROUTINE TDADRN (DATA,IMAX,JMAX,IDIM,JDIM,KDIM,DLOW,DHGH,RNDM)
C
        DIMENSION DATA(IMAX,JMAX,*)
C
C This subroutine is called to add a random quantity to each element of
C the IDIM x JDIM x KDIM array DATA, whose first two FORTRAN dimensions
C are IMAX and JMAX.
C
C The magnitude of the random quantity is stated as a multiple (RNDM)
C of the range of the data.
C
        DMUL=RNDM*(DHGH-DLOW)
C
        DMIN=+1.E36
        DMAX=-1.E36
C
C Add random quantities to the data.
C
        DO 103 I=1,IDIM
          DO 102 J=1,JDIM
            DO 101 K=1,KDIM
              DATA(I,J,K)=DATA(I,J,K)+DMUL*FRAN()
              DMIN=MIN(DMIN,DATA(I,J,K))
              DMAX=MAX(DMAX,DATA(I,J,K))
  101       CONTINUE
  102     CONTINUE
  103   CONTINUE
C
C Readjust the data to have the high and low values as before.
C
        DO 106 I=1,IDIM
          DO 105 J=1,JDIM
            DO 104 K=1,KDIM
              DATA(I,J,K)=(DATA(I,J,K)-DMIN)/(DMAX-DMIN)*
     +                                       (DHGH-DLOW)+DLOW
  104       CONTINUE
  105     CONTINUE
  106   CONTINUE
C
C Done.
C
        RETURN
C
      END


      FUNCTION FRAN()
C
C Pseudo-random-number generator.
C
        DOUBLE PRECISION X
        SAVE X
        DATA X / 2.718281828459045D0 /
        X=MOD(9821.D0*X+.211327D0,1.D0)
        FRAN=REAL(X)
        RETURN
      END


      SUBROUTINE TDRDIN (IVAL,IDEF)
C
C This routine reads a line of input to get an integer value IVAL.  If
C the line is entirely blank, the default value IDEF is returned.  If
C the line contains an illegal character, another line is obtained.
C
        CHARACTER*128 LINE
C
C Assume the default value will be returned until we find otherwise.
C
        IVAL=IDEF
C
C Get a line from standard input.
C
        READ '(A128)' , LINE
C
        IF (LINE.NE.' ') THEN
          READ (LINE,*) IVAL
        END IF
C
C Done.
C
        RETURN
C
      END


      SUBROUTINE TDRDRN (RVAL,RDEF)
C
C This routine reads a line of input to get a real number RVAL.  If
C the line is entirely blank, the default value RDEF is returned.  If
C the line contains an illegal character, another line is obtained.
C
        CHARACTER*128 LINE
C
C Assume the default value will be returned until we find otherwise.
C
        RVAL=RDEF
C
C Get a line from standard input.
C
        READ '(A128)' , LINE
C
        IF (LINE.NE.' ') THEN
          READ (LINE,*) RVAL
        END IF
C
C Done.
C
        RETURN
C
      END
