Logo Search packages:      
Sourcecode: raster3d version File versions  Download package

render.f

      MODULE LISTS
        INTEGER, ALLOCATABLE, DIMENSION(:,:)   :: KOUNT, MOUNT
        INTEGER, ALLOCATABLE, DIMENSION(:,:)   :: TTRANS
        INTEGER ISTRANS
      END MODULE LISTS

      PROGRAM RENDER
      USE LISTS
*
*     Version 3.0  (14 Dec 2010)
*
* EAM May 1990    - add object type CYLIND (cylinder with rounded ends)
*             and CYLFLAT (cylinder with flat ends)
*           - use Phong shading on triangles if they are sequentially
*             adjoining.
* EAM Feb 1991    - port to Ultrix (minor changes) and modify output format
*             to depend on code in separate module "local".
* EAM Mar 1993    - fix embarrassingly stupid bug in cylinder shadowing
*             add object type PLANE (triangle with infinite extent)
* EAM Nov 1993    - Version 2.0 (beta test)
*             fix bug which allowed objects to shadow themselves
*             Command line options for 3 output modes
* EAM Apr 1994    - Version 2.0 release
*             TIFF output support in local.c
*             minor changes to fortran source to make IBM xlf compiler happy
* EAM Sep 1994    - EXPERIMENTAL VERSION WITH OBJECT TYPES 7, 8, 9
* EAM Jan 1995    - move DATA statement to make linux happy
* EAM Mar 1995    - fix bug in routine CYLTEST which took bites out of cylinder ends
* EAM May 1995    - fold object types 7/8/9 back into distributed code for V2.1
* EAM Jan 1996    - Version 2.2
*             Add code for transparency and faster MATERIAL bookkeeping
*             Also fix major problems with explicit surface normals
*             object type 8 expanded to describe transparency
* EAM May 1996    - antialiasing scheme 4, file indirection, 
*             minor changes to accommodate HPUX
* EAM Oct 1996  - trap and forgive shadowing error due to too small NSX,NSY
* EAM Nov 1996    - scheme 0 causes alpha blend channel in output image
*             per-tile count of transparent objects
* EAM Jan 1997    - zero length cylinders treated as spheres
*             Object types 10 + 11 (fonts and labels) accepted but ignored
*             Object type  12 reserved for other label information
*             Material properties can override object colors
*           - Material OPT(1) = 1 transparency option
*             OPT(4) = continuation lines for more material properties
* EAM Mar 1997    - Make SLOP larger, and dependent on tile size
*           - GLOW light source specified by object type 13
* EAM May 1997    - V2.3c allow multiple glow lights; make cyltest a function;
*             V2.3d remove DATA statements; more terse output; 
*             BACKFACE material option; EYEPOS = 0.0 disables perspective
* EAM Jul 1997    - add commons LISTS MATRICES NICETIES RASTER
*           - D_LINES code for quadric surfaces, ISOLATION 
* EAM Sep 1997    - fix normals of flat cylinder ends (thanks to Takaaki Fukami)
* EAM Nov 1997    - add VERTEXRGB object type to extend triangle descriptions,
*             allow # as comment delimiter in input stream
* EAM Dec 1997    - Release V2.4b
* EAM Feb 1998    - fixed bug in check against limiting radius of quadrics
* EAM Jul 1998    - Object type 16 (GLOBAL PROPERTIES); fog
* EAM Aug 1998    - render back side of transparent flat-ended cylinders
*             by duplicating object with INSIDE flag bit set.
* EAM Oct 1998    - check environmental variable R3D_LIB for input file indirection
*           - V2.4g includes some preliminary code to support Z-clipping
* EAM Nov 1998    - more work on Z-clipping
* EAM Feb 1999    - re-work output module local.c to support -jpeg and -out
* EAM May 1999    - allow explicit vertex colors for cylinders also
* EAM Jul 1999    - 2.4l preliminary work towards an after-the-fact rotation option
* EAM Sep 1999    - 2.5a command line parsing in separate routine parse()
*             label processing folded into render; routines in r3dtops.f
* EAM Jan 2000    - 2.5b general release
* EAM Feb 2000    - 2.5c (bug fixes to 2.5b)
* EAM Mar 2000    - object types VTRANSP (18) and ISOLATE2 (19)
* EAM Sep 2000    - 2.5e uncompress indirect files ending with .Z or .gz
*             discard BACKCLIP objects on input, implement BACKCLIP material
* EAM Nov 2000    - 2.5f more bug-fixes to rotation of surface normals
* EAM Feb 2001    - 2.5g command line shadows
* EAM Mar 2001    - V2.6 (alpha test)
*             bounding planes
*             revamped transparency code, remove limit of 2 stacked objects
*             break out shared maximum dimensions to paramters.incl
*             make back surface HIDING (former INMODE=4) the default for
*                 non-bounded opaque triangles 
* EAM Jul 2001    - V2.6b first release
* EAM Aug 2001    - V2.6c bug-fix for MOPT1 processing
*             PNG output format ( -DPNG_SUPPORT )
* EAM Feb 2002    - allow a little RIBBONSLOP in testing for ribbon triangles
*             fix up complicated corner cases in bounded surface algorithm
* EAM Apr 2002    - clean up auto-tiling, and allow zero NPX or NPY to trigger it
* EAM Apr 2006  - V2.6d gfortran accommodations
*               - Change AND() to iand() everywhere
*               - Change  OR() to  ior() everywhere
* EAM Mar 2008    - initialize various static storage areas
* FZ  Dec 2009    - initialize more static storage areas (valgrind runs)
* FZ  Feb 2010  - expandable memory allocation 
* JMK Dec 2010    - Bugfix for JUSTCLIPPED reported by Joe Krahn
* JMK Dec 2010    - More transparency algorithms; use OPT[2] to pass choice
* EAM Dec 2010    - Read in all header info before initializing output file
*             
*     General organization:
*
*     - read in control parameters and initial output image file
*     - read in list of objects
*       - count objects that may impinge on each tile
*       - do this for both pixel and rotated "shadow" space
*     - sort objects
*     - go through main object list in sorted order
*       - fill in short lists of objects
*     - repeat the sort etc. for the objects in shadow space
*     - that's it for the "cheap" part
*     - for each tile:
*       - for each pixel:
*         - search objects to find highest point for pixel
*       - if it's transparent find the next one down as well
*         - transform resulting (x,y,z) to shadow space
*         - find closest z' for new x',y'
*         - this tells you if the pixel is in shadow or not
*         - shade accordingly
*       - copy tile to output buffer
*
*
*     Easy-to-change constants (kept in file parameters.incl):
*
*     - maximum size of any expandable array       MAXMEM
*     - maximum number of tiles in each direction  MAXNTX*, MAXNTY*
*     - number of shadow tiles in each direction   NSX, NSY
*     - maximum number of pixels in a tile         MAXNPX*, MAXNPY*
*     - maximum number of objects                  MAXOBJ*
*     - maximum number of material specifications  MAXMAT*
*     - maximum depth of stack transparent objects MAXTRANSP
*
*     *NOTE: MAXNTX,MAXNTY,MAXNPX,MAXNPY,MAXOBJ and MAXMAT are now
*     initial array sizes for expandable arrays, rather than limits.
*
*     Input (line by line except where noted):
*
*     - TITLE    anything you like
*     - NTX,NTY  tiles in each direction
*     - NPX,NPY  pixels per tile to compute in each direction
*     - SCHEME   pixel averaging scheme (1, 2, or 3)
*       - 0 no anti-aliasing, use alpha channel
*       - 1 no anti-aliasing, no alpha channel
*       - 2 means 2x2 computing pixels for 1 output pixel
*       - 3 means 3x3 computing pixels for 2x2 output pixels
*     - 4 same as 3, but NTX,NTY expanded inside program
*     - BKGND    background colour (r,g,b in range 0 to 1)
*     - SHADOW   "shadow mode?" (T or F)
*     - IPHONG   Phong power (e.g., 20)
*     - STRAIT   straight-on (2ndary) light component (e.g., 0.1)
*     - AMBIEN   ambient light component (e.g., 0.05)
*     - SPECLR   specular reflection component (e.g., 0.30)
*     - EYEPOS   eye position along +z coordinate (e.g., 4)
*       - relative to 1=narrow dimension of screen
*       - used for perspective, EYEPOS = 0.0 disables perspective
*     - SOURCE   main light source position (x,y,z components)
*       - vector length ignored, point source is at infinity
*     - TMAT     global transformation on input objects
*       - postfix 4x4 matrix on 4 lines, as you would write it
*       - upper left 3x3 must be pure rotation
*       - lower left 1x3 is translation
*       - lower right 1x1 is global scaling (reduction)
*       - upper right 3x1 causes extra perspective (should be 0)
*       - applies to homogeneous co-ordinates (x,y,z,1)
*       - radii are only scaled down by global scaling TMAT(4,4)
*     - INMODE   object input mode (1, 2, or 3)
*       - mode 1:  all objects are triangles
*       - mode 2:  all objects are spheres
*       - mode 3:  each object will be preceded by type
*     - INFMT or INFMTS   object input format(s), 1 per line
*       - one format for modes 1 and 2, or three for mode 3
*       - each is fortran format in parentheses, or single *
*       - for 3 formats, the order of formats and details is:
*         - triangle:  x1,y1,z1,x2,y2,z2,x3,y3,z3,r,g,b
*         - sphere:    x,y,z,radius,r,g,b
*         - trcone:    x1,y1,z1,rad1,x2,y2,z2,rad2,r,g,b
*       - cylinder:  as truncated cone, but 2nd radius ignored
*     - objects
*       - modes 1,2:  each object starts on a new line
*         - read according to the single format given
*       - mode 3:  each object is preceded by a line giving type
*         - type  1: triangle (to be read with 1st format)
*         - type  2: sphere (to be read with 2nd format)
*         - type  3: cylinder with rounded ends (3rd format)
*         - type  4: [not implemented: truncated cone]
*         - type  5: cylinder with flat ends (3rd format)
*         - type  6: plane (=triangle with infinite extent) (1st format)
*         - type  7: normal vectors for previous triangle (1st format)
*         - type  8: material definition which applies to subsequent objects
*         - type  9: end previous material
*         - type 10: font selection (ignored in render)
*         - type 11: label (ignored other than to count them)
*         - type 12: (reserved for additional label processing)
*         - type 13: glow light source
*       - type 14: quadric surface (usually an ellipsoid)
*       - type 15: disable coordinate transformation of subsequent objects
*       - type 16: global properties (e.g. FOG)
*       - type 17: RGB triple for each vertex of preceding object
*       - type 18: transparency at each vertex of preceding object
*       - type 19: variant of type 15; forces unitary coordinate sytem
*         - type  0: no more objects (equivalent to eof)
*
*-----------------------------------------------------------------------------
* EAM Sep 1994
*
*     1) Object type 7 signals an extra record giving explicit vertex normals
*     for a single triangle. This extra record must directly follow the 
*     corresponding triangle and uses the same format.
*
*     2) Object type 8 signals an extra record giving extra or more explicit
*     material properties object.  Current (trial) contents of record are:
*         MPHONG  - overrides global Phong power for specular reflections
*         MSPEC   - overrides global specular scattering contribution
*         SR,SG,SB      - red/green/blue components for specular highlighting
*                   (values <0 cause highlights to match object colour)
*           CLRITY      - 0.0 (opaque) => 1.0 (transparent)
*         CLROPT  - [reserved] suboptions for transparency handling
*         OPT2    - [reserved] suboptions for bounding planes  Is it really used?
*         OPT3    - [reserved]
*         OPT4    - # of additional modifier records immediately following
*     These material properties remain in effect for subsequent objects 
*     until object type 9 appears to terminate the effect.
*
*     3) Object type 9 terminates effect of previous materials property
*
*-----------------------------------------------------------------------------
*     Object types 10 and 11 are used for specifying labels.
*     Label object types are
*       - type 10: Font_Name size alignment
*       - type 11: XYZ RGB on first line
*                label (ascii characters enclosed in quotes) on second line
*     Object type 12 is reserved to go with this, as I have a nagging
*     suspicion more information may turn out to be necessary.
*-----------------------------------------------------------------------------
*     Object type 13 specifies a "glow" light source; i.e. a non-shadowing
*     light source with finite origin GLOWSRC and illumination range GLOWRAD.
*     Specular highlights for this source specified by GLOWCOL and GPHONG.
*     0.0 < GLOW < 1.0   = contribution of glow to total lighting model.
*     13
*      GLOWSRC(3)  GLOWRAD  GLOW  GOPT GPHONG  GLOWCOL(3)
*-----------------------------------------------------------------------------
* V2.4
*     Object type 14 specifies a quadric surface
*         QQ = Ax^2 + By^2 + Cz^2 + 2Dxy + 2Eyz + 2Fxz + 2Gx + 2Hy + 2Iz + J
*     centered at XC,YC,ZC and truncated at bounding radius RC. Supporting
*     code is in file quadric.f
*     14
*       XC YC ZC RC  RED GRN BLU
*       A B C D E F G H I J
*
*     Object type 15 is a single line signaling that subsequent objects are
*     not to be transformed by the TMAT matrix in the header. This isolation
*     from TMAT is terminated by an end material record (object type 9).
*-----------------------------------------------------------------------------
* V2.6
*     Object type 4 is used internally to implement bounding planes.
*     The BOUNDING_PLANE definition is given in the input stream as
*     a modifier to a MATERIAL descriptor.  At the time it is read in
*     render creates a new object of type 4 to hold the bounding 
*     plane parameters and modifiers.
*     DETAIL(K+...) is loaded with the following parameters:
*     SUBTYPE, BPTYPE, X,Y,Z, XNORM, YNORM, ZNORM, RED, GREEN, BLUE
*-----------------------------------------------------------------------------
*
*     Object space convention:
*
*     - this is the space TMAT is supposed to map your data to
*     - centre of "virtual screen" is (0,0,0)
*     - x to right, y to top, z towards viewer
*     - the smaller of the x and y dimensions goes from -.5 to +.5
*     - z cuts off at +1 and -1 by default, but is modified by FRONTCLIP, BACKCLIP
*     - shadow box dimensions determined by NSX/NTX, NSY/NTY
*
*-----------------------------------------------------------------------------
* David Bacon's comments from Version 1.0
*
*     Bugs:
*
*     - perspective is applied to raw objects, giving wrong lighting
*     - perspective unity factor is always at z=0 in object space
*     - shadow box doesn't necessarily enclose entire view prism
*     - some ASSERT calls are commented out for extra speed
*     - SLOP parameter is empirical fix to imprecision in shadowing
*
*     Deferred priorities:
*
*     - superior pixel averaging (you should use another pass?)
*     - better assignment of triangles to tiles (do clipping?)
*     - better estimate of max. triangle elevation within tile
*
*
*     Why I don't do shadowing properly:
*
*       Although you might not notice it as a casual observer,
*     the main light source is (in effect) in different places
*     for different parts of the picture.  More precisely,
*     perspective is applied to the objects comprising the
*     picture first, and THEN the lighting from a distant
*     light source is applied.  The lighting would be correct
*     if there were no perspective, because the angle doesn't
*     change across the picture.  With a scene in perspective,
*     however, the angle of the light beam, from the eye's
*     viewpoint, should vary a little.
*       The reason I allow this error is because the use of "tiles"
*     is implicitly a parallel projection, so I have to apply the
*     perspective to the objects initially.
*       Conceivably I could "undo" this whenever taking the light
*     source point of view (this would result in using a slightly
*     different light source position for different parts of the
*     picture), but that would cause another problem:
*       Perspective distorts spheres differently depending on whether
*     you ask about each point on the surface individually or use one
*     perspective factor for all points based on the centre.  Consider
*     even a sphere in the plane where perspective is supposed to be
*     unity (z=0 for us).  It swells slightly when perspective is
*     applied point by point.
*       This is a serious problem for us because although I
*     generate non-swelled spheres by applying the perspective to
*     all of them initially, I end up asking about individual
*     points on them when wanting the light source point of view.
*     You might argue that I could simply calculate the amount
*     of swelling that has to be accounted for (by drawing tangents
*     from the eye and seeing where they intersect the constant-z
*     plane that passes through the sphere centre or something),
*     but the "swelling" is complicated in the sense that it is
*     in effect a "stretching" of the surface closest to the eye
*     and a "shrinking" in back.
*       Even if I could see how to compensate for all this, I
*     don't think it would be worth it.
*       It's probably not even worth trying to implement a better
*     approximate solution by changing the light source position
*     slightly for each object in the picture.  The ambiguity
*     as to whether the obscuring or the obscured object should have
*     the modified light source position applied would make it
*     difficult to assign objects to shadow tiles.  Trying to
*     implement full "antiperspective" for the light source would
*     just shift the problems of perspective (swelling spheres,
*     etc.) to a different locale without solving them.
*-----------------------------------------------------------------------------
*
*     Overkill:
      IMPLICIT REAL   (A-H, O-Z)
*
      INCLUDE 'VERSION.incl'
*
*     I/O units for control input, info output, label processing
      INTEGER INPUT, INPUT0, NOISE, LUNIT
      PARAMETER (INPUT0=5, NOISE=0, LUNIT=4)
*
*     Descriptor codes for the various object types
      INTEGER TRIANG, SPHERE, INTERNAL, CYLIND, CYLFLAT
      INTEGER PLANE, QUADRIC, MXTYPE, FONT, GLOWLIGHT
      INTEGER NORMS, VERTEXRGB, VERTRANSP
      PARAMETER (TRIANG   = 1)
      PARAMETER (SPHERE   = 2)
      PARAMETER (CYLIND   = 3)
      PARAMETER (INTERNAL = 4)
      PARAMETER (CYLFLAT  = 5)
      PARAMETER (PLANE    = 6)
      PARAMETER (NORMS    = 7)
      PARAMETER (MATERIAL = 8)
      PARAMETER (MATEND   = 9)
      PARAMETER (FONT     = 10, LABEL = 11)
      PARAMETER (GLOWLIGHT= 13)
      PARAMETER (QUADRIC  = 14)
      PARAMETER (ISOLATE1 = 15)
      PARAMETER (GPROP    = 16)
      PARAMETER (VERTEXRGB= 17)
      PARAMETER (VERTRANSP= 18)
      PARAMETER (ISOLATE2 = 19)
      PARAMETER (MXTYPE   = 19)
*
*     Bit definitions for FLAG array
      INTEGER    FLAT,      RIBBON,    SURFACE,   PROPS
      PARAMETER (FLAT=2,    RIBBON=4,  SURFACE=8, PROPS=16)
      INTEGER    TRANSP,    HIDDEN,    INSIDE,    MOPT1
      PARAMETER (TRANSP=32, HIDDEN=64, INSIDE=128,MOPT1=256)
      INTEGER      VCOLS,     CLIPPED,      VTRANS,      BOUNDED
      PARAMETER   (VCOLS=512, CLIPPED=1024, VTRANS=2048, BOUNDED=4096)
*
*     Bit definitions for OTMODE passed to local(1,...)
      INTEGER    ALPHACHANNEL
      PARAMETER (ALPHACHANNEL=32)
*
*     $$$$$$$$$$$$$   ARRAY SIZE LIMITS $$$$$$$$$$$$$$
*
      INCLUDE 'parameters.incl'
      INTEGER OUTSIZ
      PARAMETER (OUTSIZ = MAXNTX*MAXNPX*MAXNPY)
*
*
*     Command line options (Aug 1999) NB: nax,nay,quality MUST be integer*2
      COMMON /OPTIONS/ FONTSCALE, GAMMA, ZOOM, NSCHEME, SHADOWFLAG, XBG,
     &                 NAX, NAY, OTMODE, QUALITY, INVERT, LFLAG
      REAL             FONTSCALE, GAMMA, ZOOM
      INTEGER          NSCHEME, SHADOWFLAG, XBG
      INTEGER*4        NAX, NAY, OTMODE, QUALITY
      LOGICAL*2        INVERT, LFLAG
*
*     Title for run
      CHARACTER*132 TITLE
*
*     Number of tiles, pixels per tile
      COMMON /RASTER/ NTX,NTY,NPX,NPY
      INTEGER         NTX,NTY,NPX,NPY
*
*     Pixels per tile after anti-aliasing, output buffer line length
      INTEGER         NOX, NOY, NX
*
*     Actual image size in pixels (may include partial tiling at the edges)
*     (MUST BE INTEGER*2 for call to local()!!!)
C     INTEGER*2 NAX, NAY
*
*     One lonely tile
      REAL, ALLOCATABLE, DIMENSION(:,:,:) :: TILE
*
*     With an alpha blend channel
      REAL, ALLOCATABLE, DIMENSION(:,:)   :: ACHAN
*
*     Pixel averaging scheme
      INTEGER SCHEME
*
*     Background colour
      REAL   BKGND(3)
      INTEGER IBKGND(3)
*
*     "Shadow mode?"
      LOGICAL SHADOW
      INTEGER NSXMAX,NSYMAX
*
*     Phong power
      INTEGER IPHONG
*
*     Straight-on (secondary) light source contribution
      REAL   STRAIT
*
*     Ambient light contribution
      REAL   AMBIEN
*
*     Specular reflection component
      REAL   SPECLR
*
*     Primary light source position
      REAL   SOURCE(3)
*
*     Input transformation
      COMMON /MATRICES/ XCENT, YCENT, SCALE, EYEPOS, SXCENT, SYCENT,
     &                  TMAT, TINV, TINVT, SROT, SRTINV, SRTINVT
     &                 ,RAFTER, TAFTER
      REAL   XCENT, YCENT, SCALE, SXCENT, SYCENT
*     Transformation matrix, inverse of transpose, and transposed inverse
      REAL   TMAT(4,4), TINV(4,4), TINVT(4,4)
*     Shortest rotation from light source to +z axis
      REAL   SROT(4,4), SRTINV(4,4), SRTINVT(4,4)
*     Post-hoc transformation on top of original TMAT
      REAL   RAFTER(4,4), TAFTER(3)
      EXTERNAL DET
      REAL     DET
*
*     Distance (in +z) of viewing eye
      REAL   EYEPOS
*
*     Input mode
      INTEGER INMODE
*
*     Buffer one line of input for decoding
      CHARACTER*132 LINE
*
*     Input format(s)
      CHARACTER*80 INFMTS(MXTYPE),INFMT
*
*     Free-format input flag
      LOGICAL INFLGS(MXTYPE),INFLG
*
*     Allow very long names for file indirection
      CHARACTER*132 FULLNAME
*
*     Stuff for shading
      REAL   NL(3),NORMAL(3),LDOTN
      REAL   RGBCUR(3),RGBSHD(3),RGBFUL(3)
      REAL   SPECOL(3)
*
*     FOG parameters 
*     (fogtype  -1 = none, 0 = linear depthcuing, 1 = exponential model)
*     (fogfront  0 = front object, else fraction of front clipping plane)
*     (fogback   0 = back object,  else fraction of back clipping plane)
*
      COMMON /FOGCOM/ FOGTYPE,FOGFRONT,FOGBACK,FOGDEN,FOGLIM,FOGRGB
      INTEGER FOGTYPE
      REAL    FOGFRONT, FOGBACK, FOGDEN, FOGLIM(2), FOGRGB(3)
*
*     The s & m guys are for the shadow box in the following
*
*     Object list, consists of pointers (less 1) into detail, sdtail
      INTEGER, ALLOCATABLE, DIMENSION(:) :: LIST, MIST
*
*     Object types and flags, parallel to list
      INTEGER,   ALLOCATABLE, DIMENSION(:) :: TYPE
      INTEGER*4, ALLOCATABLE, DIMENSION(:) :: FLAG
*
*     Keep a separate list of special materials
*     and remember any special props of current material on input
CDEBUG MPARITY gets its own array because it's used in a per-pixel loop
CDEBUG (using DETAIL(LIST(MLIST(MAT))+18) cost 5% in execution time)
      INTEGER, ALLOCATABLE, DIMENSION(:) :: MLIST, MPARITY
      LOGICAL MATCOL, BACKFACE
      LOGICAL CLIPPING, MAYCLIP, JUSTCLIPPED
      REAL    RGBMAT(3)
*
*     Object details, shadow object details
      REAL, ALLOCATABLE, DIMENSION(:) :: DETAIL, SDTAIL
*
*     Input buffer for details
      REAL   BUF(100)
*
*     Number of objects in each tile's short list (m... are for shadows)
C     Moved from COMMON BLOCK LISTS to MODULE LISTS to allow dynamic
C     allocation - FZ
C     COMMON /LISTS/ KOUNT, MOUNT, TTRANS, ISTRANS
C     INTEGER KOUNT(MAXNTX,MAXNTY), MOUNT(NSX,NSY)
C     INTEGER TTRANS(MAXNTX,MAXNTY), ISTRANS
*
*     Pointer to where each tile's objects start
      INTEGER, ALLOCATABLE, DIMENSION(:,:) :: KSTART, MSTART
*
*     Pointer to where each tile's objects end
      INTEGER, ALLOCATABLE, DIMENSION(:,:) :: KSTOP, MSTOP
*
*     Short list heap
      INTEGER, ALLOCATABLE, DIMENSION(:) :: KSHORT, MSHORT
*
*     Temporary for sorting
      REAL, ALLOCATABLE, DIMENSION(:) :: ZTEMP
*
*     Where the permutation representing the sort is stored
      INTEGER, ALLOCATABLE, DIMENSION(:) :: ZINDEX
*
*     The number of "details" each object type is supposed to have
*     :       input,        object,       shadow
      INTEGER IDET(MXTYPE), KDET(MXTYPE), SDET(MXTYPE)
*
*     Support for cylinders
      EXTERNAL CYLTEST
      LOGICAL  CYLTEST, ISCYL
*
*     Support for quadrics
      REAL     QNORM(3)
      EXTERNAL QINP, QTEST
      LOGICAL  QINP, QTEST, ISQUAD
*
*     Support for transparency
      COMMON /TRANS/ NTRANSP, INDEPTH, INDTOP, TRANOVFL, ZTOP, ZHIGH,
     &               INDLIST(MAXTRANSP), ZLIST(MAXTRANSP), 
     &                 NORMLIST(3,MAXTRANSP)
      INTEGER        NTRANSP, INDEPTH, INDTOP, INDLIST, TRANOVFL
      REAL           ZTOP, ZHIGH, ZLIST, NORMLIST
      REAL    SBLEND, RGBLND(3)
*
* Support for a "glow" light source 
      REAL  GLOWSRC(3), GLOWCOL(3), GDIST(3), GLOWRAD, GLOW, GLOWMAX
      INTEGER     GOPT, GPHONG
      INTEGER,    ALLOCATABLE, DIMENSION(:) :: GLOWLIST
      INTEGER     NGLOWS
*
* Support for decompression on the fly
      EXTERNAL ungz
      INTEGER  ungz
*
* Support for BOUNDING_PLANE internal object type
      REAL     BPLANE(3), BPNORM(3), BPRGB(3)
      REAL     xn, yn, zn, xnb, ynb, znb
      INTEGER  BPTYPE, NBOUNDS, BPIND
      LOGICAL  ORTEPLIKE
      INTEGER  ORTEP
      PARAMETER (ORTEP=1)
      EXTERNAL INBOUNDS
      LOGICAL  INBOUNDS
      REAL     TEMPNORM(3)
*
*     Output buffer
      INTEGER*2, ALLOCATABLE, DIMENSION(:,:) :: OUTBUF
*
*     Copy of NOISE for ASSERT to see
      INTEGER ASSOUT
      LOGICAL VERBOSE
      COMMON /ASSCOM/ ASSOUT, VERBOSE
      SAVE /ASSCOM/
*
*     For label processing
      COMMON /LABELS/ LABOUT
      INTEGER         LABOUT
*
*     Gamma correction
      INTEGER GAMMA_MAP(256)
      PARAMETER (MAXRGB=255.0)
      LOGICAL GAMMACORRECTION
*
*     Keep track of actual coordinate limits
      COMMON /NICETIES/ TRULIM,      ZLIM,    FRONTCLIP, BACKCLIP
     &                , ISOLATION
      REAL              TRULIM(3,2), ZLIM(2), FRONTCLIP, BACKCLIP
      INTEGER           ISOLATION
*
*     Array of sizes to try allocating for expanded dynamic storage
      INTEGER NEEDMEM, TRY1(3), TRY2(3)
      INTEGER,   ALLOCATABLE, DIMENSION(:)     :: TMP1D
      INTEGER*4, ALLOCATABLE, DIMENSION(:)     :: TMP1DI4
      REAL,      ALLOCATABLE, DIMENSION(:)     :: TMP1DR
      INTEGER,   ALLOCATABLE, DIMENSION(:,:)   :: TMP2D
      INTEGER*2, ALLOCATABLE, DIMENSION(:,:)   :: TMP2DI2
      REAL,      ALLOCATABLE, DIMENSION(:,:)   :: TMP2DR
      REAL,      ALLOCATABLE, DIMENSION(:,:,:) :: TMP3DR

      LOGICAL TEST_ALLOC
      TEST_ALLOC = .TRUE.
*
*     Allocate initial space for dynamically allocatable arrays
      ALLOCATE( TILE(3,MAXNPX,MAXNPY) )
      ALLOCATE( ACHAN(MAXNPX,MAXNPY) )
      ALLOCATE( ZTEMP(MAXOBJ) )
      ALLOCATE( ZINDEX(MAXOBJ) )
      ALLOCATE( LIST(MAXOBJ), MIST(MAXOBJ) )
      ALLOCATE( TYPE(MAXOBJ) )
      ALLOCATE( FLAG(MAXOBJ) )
      ALLOCATE( MLIST(MAXMAT), MPARITY(MAXMAT) )
      ALLOCATE( DETAIL(MAXDET), SDTAIL(MAXSDT) )
      ALLOCATE( KOUNT(MAXNTX,MAXNTY), MOUNT(NSX,NSY) )
      ALLOCATE( TTRANS(MAXNTX,MAXNTY) )
      ALLOCATE( KSTART(MAXNTX,MAXNTY), MSTART(NSX,NSY) )
      ALLOCATE( KSTOP(MAXNTX,MAXNTY), MSTOP(NSX,NSY) )
      ALLOCATE( KSHORT(MAXSHR), MSHORT(MAXSSL) )
      ALLOCATE( GLOWLIST(MAXGLOWS) )
      ALLOCATE( OUTBUF(OUTSIZ,4) )

      TRULIM (1,1) = HUGE
      TRULIM (2,1) = HUGE
      TRULIM (3,1) = HUGE
      TRULIM (1,2) = -HUGE
      TRULIM (2,2) = -HUGE
      TRULIM (3,2) = -HUGE
      ZLIM(1)      = HUGE
      ZLIM(2)      = -HUGE
*
      IDET(TRIANG)   = 12
      IDET(SPHERE)   = 7
      IDET(CYLFLAT)  = 11
      IDET(CYLIND)   = 11
      IDET(PLANE)    = 12
      IDET(NORMS )   = 9
      IDET(MATERIAL) = 10
      IDET(GLOWLIGHT)= 10
      IDET(QUADRIC)  = 17
      IDET(VERTEXRGB)= 9
      IDET(VERTRANSP)= 3
      IDET(INTERNAL) = 20
*
      KDET(TRIANG)   = 16
      KDET(SPHERE)   = 7
      KDET(CYLIND)   = 11
      KDET(PLANE)    = 7
      KDET(NORMS )   = 9
      KDET(MATERIAL) = 18
      KDET(GLOWLIGHT)= 10
      KDET(QUADRIC)  = 17
      KDET(VERTEXRGB)= 9
      KDET(VERTRANSP)= 3
      KDET(INTERNAL) = 20
*
      SDET(TRIANG)   = 13
      SDET(SPHERE)   = 4
      SDET(CYLIND)   = 8
      SDET(QUADRIC)  = 14
      SDET(INTERNAL) = 20
*     These object types really have no shadow details, 
*     but indexing seems to require a nonzero value
      SDET(PLANE)    = 1
      SDET(NORMS )   = 1
      SDET(MATERIAL) = 1
      SDET(GLOWLIGHT)= 1
      SDET(VERTEXRGB)= 1
      SDET(VERTRANSP)= 1
*
*     Feb 2008 - more initializations for current gfortran
      NSXMAX = 0
      NSYMAX = 0
      CLROPT = 0
      SCHEME = 0
      TRNSPOPT = 0
      INFLG = .TRUE.
      JUSTCLIPPED = .FALSE.
      IXHI = 0
      IXLO = 0
      IYHI = 0
      IYLO = 0
      ITPASS = 0
*
*     Copy the info (also error reporting) unit number to common
      ASSOUT = NOISE
      WRITE (NOISE,*) ' '
*
*     Initialize to level 0 of file indirection
      INPUT = INPUT0
*
*     Initialize unit number for label processing
      LABOUT = LUNIT
*
*     Initialize to no special material properties
      MSTATE    = 0
      MATCOL    = .FALSE.
      ISOLATION = 0
      CLIPPING  = .FALSE.
      NBOUNDS   = 0
      ORTEPLIKE = .FALSE.
      CLRITY    = 0.0
      GLOWMAX   = 0.0
*
*     Initialize to no perspective. EYEPOS > 0 will add perspective
      PFAC  = 1.0
      PFAC1 = 1.0
      PFAC2 = 1.0
      PFAC3 = 1.0
*
*     Initialize global properties
      FOGTYPE = -1

      RGBLND(1) = 0
      RGBLND(2) = 0
      RGBLND(3) = 0
*
*     EAM Aug 1999 - break out command line parsing into new routine
      call parse
*
*     Get title
  100 CONTINUE
      DO I=1,132
        TITLE(I:I) = ' '
      ENDDO

      READ (INPUT,'(A)',END=104,ERR=104) TITLE
      IF (TITLE(1:1) .EQ. '#') GOTO 100
      IF (TITLE(1:1) .EQ. '@') THEN
      J = 1
      K = 132
      DO I=132,2,-1
        IF (TITLE(I:I).EQ.'   ') TITLE(I:I) = ' '
        IF (TITLE(I:I).NE.' ') J = I
        IF (TITLE(I:I).EQ.'#') K = I-1
        IF (TITLE(I:I).EQ.'!') K = I-1
      ENDDO
      IF (J.EQ.1) GOTO 101
      DO WHILE (TITLE(K:K).EQ.' ')
        K = K -1
      ENDDO
      OPEN (UNIT=INPUT+1,ERR=101,STATUS='OLD',FILE=TITLE(J:K))
      WRITE (NOISE,'(A,A)') '  + Opening input file ',TITLE(J:K)
      INPUT = INPUT + 1
      READ (INPUT,'(A)',ERR=101) TITLE
      IF (TITLE(1:1) .EQ. '#') GOTO 100
      ENDIF
      GOTO 102
  101 WRITE (NOISE,'(A,A)')' >> Cannot open or read file ',TITLE(2:K)
      CALL EXIT(-1)
  102 CONTINUE
      K = 132
      DO WHILE (TITLE(K:K).EQ.' ')
      K = K -1
      ENDDO
      WRITE (NOISE,103) TITLE(1:K)
  103 FORMAT('title="',A,'"')
*
*     Get number of tiles
      READ (INPUT,*,ERR=104,END=104) NTX,NTY
      CALL ASSERT (NTX.GT.0, 'ntx.le.0')
      CALL ASSERT (NTY.GT.0, 'nty.le.0')
      GOTO 105
104   CALL ASSERT(.FALSE.,
     &           '>>> This doesnt look like a Raster3D input file! <<<')
105   CONTINUE
*
*     Get number of pixels per tile - 0 means autotile from values in NTX, NTY
      READ (INPUT,*,ERR=104) NPX,NPY
      if (npx.eq.0 .and. nax.le.0) nax = ntx
      if (npy.eq.0 .and. nay.le.0) nay = nty
*
*     Get pixel averaging scheme
      READ (INPUT,*,ERR=104) SCHEME
      if (nscheme.ge.0) scheme = nscheme
      CALL ASSERT (SCHEME.GE.0 .AND. SCHEME.LE.4, 'bad scheme')
*
* Set up tiling and anti-aliasing.
* If NAX, NAY are set, then use them to autotile
      IF (SCHEME.LE.1) THEN
        call autotile( nax, nay, 2 )
        NOX = NPX
        NOY = NPY
      if (nax.lt.0) nax = npx * ntx
      if (nay.lt.0) nay = npy * nty
      ELSEIF (SCHEME.EQ.2) THEN
      if (nax.gt.0) nax = nax * 2
      if (nay.gt.0) nay = nay * 2
      if (nax.le.0) nax = NPX * NTX
      if (nay.le.0) nay = NPY * NTY
        call autotile( nax, nay, 2 )
      nax = nax / 2
      nay = nay / 2
        NOX = NPX/2
        NOY = NPY/2
      ELSEIF (SCHEME.EQ.3 .and. nscheme.ne.-4 .and. nax.le.0 ) THEN
C     Old style scheme 3 with exact tiling specified
      nax = NPX * NTX
      nay = NPY * NTY
      call autotile( nax, nay, 3 )
      nax = (nax * 2 + 2) / 3
      nay = (nay * 2 + 2) / 3
      NOX = (NPX * 2) / 3
      NOY = (NPY * 2) / 3
      ELSE
C     Either scheme 4 or -size and anti-aliasing selected on command line
      if (nax.gt.0) nax = (nax + (nax+1)/2)
      if (nay.gt.0) nay = (nay + (nay+1)/2)
      if (nax.le.0) nax = NPX*NTX + (NPX*NTX+1)/2
      if (nay.le.0) nay = NPY*NTY + (NPY*NTY+1)/2
      call autotile( nax, nay, 3 )
      nax = (nax * 2 + 2) / 3
      nay = (nay * 2 + 2) / 3
      NOX = (NPX * 2) / 3
      NOY = (NPY * 2) / 3
      SCHEME = 3
      ENDIF
*
      call assert (nax.gt.0, 'nax <= 0')
      call assert (nay.gt.0, 'nay <= 0')
*
      CALL ASSERT (NTX.GT.0.,'Tiling failure - ntx = 0')
      CALL ASSERT (NTY.GT.0.,'Tiling failure - nty = 0')
      CALL ASSERT (NPX.GT.0.,'Tiling failure - npx = 0')
      CALL ASSERT (NPY.GT.0.,'Tiling failure - npy = 0')
*
*     Expand arrays KOUNT, TTRANS, KSTART and KSTOP if needed, eg NTX > MAXNTX
      if (NTX.GT.SIZE(KSTOP,1) .OR. NTY.GT.SIZE(KSTOP,2)) THEN
*
*         Double the old allocation if that's enough, or take 150% of the new
*         If that fails, try 150% of the old or 120% of new, or just the new
          CALL GET_TRY(SIZE(KSTOP, 1), NTX, TRY1, 2)
          CALL GET_TRY(SIZE(KSTOP, 2), NTY, TRY2, 2)
*
*         Try allocating memory to the arrays
          do 900 ITRY = 1,3
*        
*             Test to see if requested allocation is valid
              if (TRY1(ITRY) .LE. 0 .OR. TRY2(ITRY) .LE. 0 .OR.
     &            TRY1(ITRY)*TRY2(ITRY) .LE. 0 .OR.
     &            MAXMEM / TRY1(ITRY) .LT. TRY2(ITRY)) GOTO 900

              ALLOCATE( TMP2D(TRY1(ITRY),TRY2(ITRY)), STAT=IERR)
              if (ierr .NE. 0) GOTO 900
C             TMP2D = 0
              TMP2D = KOUNT 
              CALL MOVE_ALLOC(from=TMP2D, to=KOUNT)

              ALLOCATE( TMP2D(TRY1(ITRY),TRY2(ITRY)), STAT=IERR)
              if (ierr .NE. 0) GOTO 900
C             TMP2D = 0
              TMP2D = TTRANS
              CALL MOVE_ALLOC(from=TMP2D, to=TTRANS)

              ALLOCATE( TMP2D(TRY1(ITRY),TRY2(ITRY)),stat=ierr)
              if (ierr .NE. 0) GOTO 900
C             TMP2D = 0
              TMP2D = KSTART
              CALL MOVE_ALLOC(from=TMP2D, to=KSTART)

              ALLOCATE( TMP2D(TRY1(ITRY),TRY2(ITRY)),stat=ierr)
              if (ierr .NE. 0) GOTO 900
C             TMP2D = 0
              TMP2D = KSTOP
              CALL MOVE_ALLOC(from=TMP2D, to=KSTOP)
              if(TEST_ALLOC)write(NOISE,*)"Expand MAXNTX x Y to ",
     &              try1(ITRY)," x ",try2(ITRY)
              GOTO 902
900       CONTINUE
      ENDIF
*
*     These should only fail if the above dynamic allocation failed
902   CALL ASSERT (NTX.LE.SIZE(KSTOP, 1),'Tiling failure - ntx>maxntx')
      CALL ASSERT (NTY.LE.SIZE(KSTOP, 2),'Tiling failure - nty>maxnty')
*
*     Expand arrays TILE, ACHAN next if needed, i.e. NPX,NPY > MAXNPX or Y
      if (NPX.GT.SIZE(TILE,2).OR. NPY.GT.SIZE(TILE,3)) THEN
          CALL GET_TRY(SIZE(TILE,2), NPX, TRY1, 2)
          CALL GET_TRY(SIZE(TILE,3), NPY, TRY2, 2)
          do 905 ITRY = 1,3
*             Test to see if requested allocation is valid
              if (TRY1(ITRY) .LE. 0 .OR. TRY2(ITRY) .LE. 0 .OR.
     &            TRY1(ITRY)*TRY2(ITRY) .LE. 0 .OR.
     &            MAXMEM / TRY1(ITRY) .LT. TRY2(ITRY)) GOTO 905

              ALLOCATE( TMP2DR(TRY1(ITRY),TRY2(ITRY)),stat=ierr)
              if (ierr .NE. 0) GOTO 905
              TMD2DR = 0.
              TMP2DR = ACHAN 
              CALL MOVE_ALLOC(from=TMP2DR, to=ACHAN)

              ALLOCATE( TMP3DR(3,TRY1(ITRY),TRY2(ITRY)),stat=ierr)
              if (ierr .NE. 0) GOTO 905
              TMD3DR = 0.
              TMP3DR = TILE
              CALL MOVE_ALLOC(from=TMP3DR, to=TILE)
              if(TEST_ALLOC)write(NOISE,*)"Expand MAXNPX,Y to ",
     &            try1(ITRY)," x ",try2(ITRY)
              GOTO 907
905       CONTINUE
      ENDIF
*
*     These should only fail if the above dynamic allocation failed
907   CALL ASSERT (NPX.LE.SIZE(TILE,2),'Tiling failure - npx>maxnpx')
      CALL ASSERT (NPY.LE.SIZE(TILE,3),'Tiling failure - npy>maxnpy')
*
      IF (VERBOSE) THEN
      WRITE (NOISE,*) 'ntx=',NTX,' nty=',NTY
      WRITE (NOISE,*) 'npx=',NPX,' npy=',NPY
      WRITE (NOISE,*) 'scheme=',SCHEME
      WRITE (NOISE,*) 'nox=',NOX,' noy=',NOY
      END IF
      if (nax.lt.0) nax = nox*ntx
      if (nay.lt.0) nay = noy*nty
      NX = nox*ntx
      LINOUT = 0
      WRITE (NOISE,1105) 'Rendered raster size =',NPX*NTX,NPY*NTY
      WRITE (NOISE,1105) '  Output raster size =',NAX,NAY
1105  FORMAT(A,I7,' x',I7)
C
*
*     Expand array OUTBUF if needed, i.e. NTX * NOX * NOY > OUTBUF size 
*     NOTE: this may avoidably fail if the above array expansions took
*     up more room than it needed, and didn't leave enough for OUTBUF.
      NEEDMEM = NOY*NOX*NTX 
      if ( NEEDMEM .GT. SIZE(OUTBUF,1) ) THEN
          CALL GET_TRY(SIZE(OUTBUF,1), NEEDMEM, TRY1, 1)
          do 910 ITRY = 1,3
*             Test to see if requested allocation is valid
              if (TRY1(ITRY).LE.0 .OR. TRY1(ITRY).GT.MAXMEM) GOTO 910

              ALLOCATE( TMP2DI2(TRY1(ITRY),4), stat=ierr )
              if (ierr .NE. 0) GOTO 910
              TMP2DI2 = 0
              TMP2DI2 = OUTBUF
              CALL MOVE_ALLOC(from=TMP2DI2, to=OUTBUF)
              if(TEST_ALLOC)write(NOISE,*)"Expand OUTBUF to ",try1(ITRY)
              GOTO 912
910       CONTINUE
      ENDIF
912   CALL ASSERT (SIZE(OUTBUF,1).GE.NEEDMEM,
     &             'image too large for output buffer')
C
C     Header records and picture title
      IF (SCHEME.EQ.0) OTMODE = ior(OTMODE,ALPHACHANNEL)
*
*     Some derived parameters
      XCENT  = NTX*NPX/2.
      YCENT  = NTY*NPY/2.
      SXCENT = NSX*NPX/2.
      SYCENT = NSY*NPY/2.
      SCALE  = 2.*MIN(XCENT,YCENT)
*     This was always true; now it's explicit
      BACKCLIP  = -(SCALE+1.0)
      FRONTCLIP =  HUGE
*     Copy scheme to common, where r3dtogd can see it
      NSCHEME = SCHEME
*
*     Get background colour
      READ (INPUT,*,ERR=104) BKGND
      if (XBG.NE.0) then
      BKGND(3) = iand(XBG,X'00FF')
      BKGND(2) = iand(XBG,X'FF00')/256
      BKGND(1) = iand(XBG,X'FF0000')/65536
      BKGND(3) = BKGND(3) / 255.
      BKGND(2) = BKGND(2) / 255.
      BKGND(1) = BKGND(1) / 255.
      BKGND(3) = BKGND(3)**2
      BKGND(2) = BKGND(2)**2
      BKGND(1) = BKGND(1)**2
      endif
      IF (VERBOSE) THEN
            WRITE (NOISE,1106) 'bkgnd=',BKGND
      END IF
1106  FORMAT(A,4F10.4,(/,4F10.4))
      CALL ASSERT (BKGND(1).GE.0., 'bkgnd(1) < 0')
      CALL ASSERT (BKGND(2).GE.0., 'bkgnd(2) < 0')
      CALL ASSERT (BKGND(3).GE.0., 'bkgnd(3) < 0')
      CALL ASSERT (BKGND(1).LE.1., 'bkgnd(1) > 1')
      CALL ASSERT (BKGND(2).LE.1., 'bkgnd(2) > 1')
      CALL ASSERT (BKGND(3).LE.1., 'bkgnd(3) > 1')
      IF (GAMMA.LT.0.99 .OR. GAMMA.GT.1.01) THEN
        IBKGND(1) = SQRT(BKGND(1)) ** (1.0/GAMMA) * MAXRGB + 0.5
        IBKGND(2) = SQRT(BKGND(2)) ** (1.0/GAMMA) * MAXRGB + 0.5
        IBKGND(3) = SQRT(BKGND(3)) ** (1.0/GAMMA) * MAXRGB + 0.5
      ELSE
          IBKGND(1) = 255. * SQRT(BKGND(1)) + .5
          IBKGND(2) = 255. * SQRT(BKGND(2)) + .5
          IBKGND(3) = 255. * SQRT(BKGND(3)) + .5
      ENDIF
*
*
*     Get "shadows" flag
      READ (INPUT,*,ERR=104) SHADOW
      IF (SHADOWFLAG.EQ.0) SHADOW = .FALSE.
      IF (SHADOWFLAG.EQ.1) SHADOW = .TRUE.
      IF (VERBOSE) THEN
            WRITE (NOISE,*) 'shadow=',SHADOW
      END IF
*
*     Get Phong power
      READ (INPUT,*,ERR=104) PHONG
      IPHONG = PHONG
      CALL ASSERT (IPHONG.GE.0, 'iphong < 0')
*     A derived constant for numerical purposes in applying the
*     Phong power in the shading algorithm.
*     The idea is that any specular contribution less than
*     1E-9 (hence the 9 in 9./IPHONG) is insignificant:
      IF (IPHONG .NE. 0) PHOBND = 0.1**(9./IPHONG)
      IF (IPHONG .EQ. 0) PHOBND = 0.
*
*     Get contribution of straight-on (secondary) light source
      READ (INPUT,*,ERR=104) STRAIT
      CALL ASSERT (STRAIT.GE.0., 'strait < 0')
      CALL ASSERT (STRAIT.LE.1., 'strait > 1')
*
*     Derive contribution of primary light source
      PRIMAR = 1. - STRAIT
*
*     Get contribution of ambient light
      READ (INPUT,*,ERR=104) AMBIEN
      CALL ASSERT (AMBIEN.GE.0., 'ambien < 0')
      CALL ASSERT (AMBIEN.LE.1., 'ambien > 1')
*
*     Get component of specular reflection
      READ (INPUT,*,ERR=104) SPECLR
      CALL ASSERT (SPECLR.GE.0., 'speclr < 0')
      CALL ASSERT (SPECLR.LE.1., 'speclr > 1')
*
      IF (VERBOSE) THEN
      WRITE (NOISE,1104) 'iphong=',float(IPHONG),'strait=',STRAIT,
     &                'ambien=',AMBIEN,'speclr=',SPECLR
1104  FORMAT(2(4X,A,F10.4))
      END IF
*
*     Derive component of diffuse reflection
      CALL ASSERT (AMBIEN+SPECLR.LE.1., 'ambien+speclr > 1')
      DIFFUS = 1. - (AMBIEN+SPECLR)
*
*     Get distance of viewing eye
      READ (INPUT,*,ERR=104) EYEPOS
      CALL ASSERT (EYEPOS.GE.0., 'eyepos.lt.0')
*
*     Get position of primary light source
      READ (INPUT,*,ERR=104) SOURCE
      SMAG = SQRT(SOURCE(1)**2 + SOURCE(2)**2 + SOURCE(3)**2)
      SOURCE(1) = SOURCE(1) / SMAG
      SOURCE(2) = SOURCE(2) / SMAG
      SOURCE(3) = SOURCE(3) / SMAG
      IF (VERBOSE) THEN
      WRITE (NOISE,1106) 'eyepos=',EYEPOS
      WRITE (NOISE,1106) 'source=',SOURCE
      WRITE (NOISE,1106) 'normalized source=',SOURCE
      END IF
*
*     Get input transformation
      DO I=1,4
        READ (INPUT,*,ERR=104) (TMAT(I,J),J=1,4)
      END DO
      IF (VERBOSE) THEN
      WRITE (NOISE,*) 'tmat (v'' = v * tmat):'
      DO I=1,4
            WRITE (NOISE,'(4f9.4)') (TMAT(I,J),J=1,4)
      END DO
      END IF
*
*     Initialize output file
      IERR = LOCAL(5, IBKGND(1), IBKGND(2), IBKGND(3))
      IERR = LOCAL(1, NAX, NAY, OTMODE, QUALITY)
      IERR = LOCAL(4, TITLE)
*
*     Allow command line rescaling option
      IF (ZOOM.LT.0.) ZOOM = -ZOOM / 100.
      IF (ZOOM.GT.0.) TMAT(4,4) = TMAT(4,4) / ZOOM
      IF (VERBOSE .AND. ZOOM.NE.0.) THEN
            WRITE (NOISE,1106) 'zoom factor = ',ZOOM
      ENDIF
*
*     EAM - The original output mode was "upside down" compared
*     to what most graphics programs expect to see.  It is messy 
*     to change the evaluation order everywhere so that pixels can be 
*     streamed to stdout, so instead I invert the Y axis in TMAT and SOURCE
*     here.  
*     The actual decision whether or not to invert is done in local.c
*     and returned as a bit in the status word returned by local(0,...)
      IF (INVERT) THEN
      DO I = 1,4
        TMAT(I,2) = -TMAT(I,2)
      ENDDO
      SOURCE(2) = -SOURCE(2)
      ENDIF
*
*     By popular demand, add a post-hoc rotation/translation option
*     that uses matrices of the form used by O and molscript
*     Initialized here to identity matrix; set by GPROP options.
      DO I = 1,4
      DO J = 1,4
        RAFTER(I,J) = 0.0
      ENDDO
      RAFTER(I,I) = 1.0
      ENDDO
      TAFTER(1) = 0.0
      TAFTER(2) = 0.0
      TAFTER(3) = 0.0
*
*     Compute the rotation matrix which takes the light
*     source to the +z axis (i.e., to the viewpoint).
*     first make p = source cross z (and normalize p)
      P1 = SOURCE(2)
      P2 = -SOURCE(1)
*     p3 = 0
      PLEN = SQRT(P1**2 + P2**2)
      IF (PLEN .GT. 0.0) P1 = P1 / PLEN
      IF (PLEN .GT. 0.0) P2 = P2 / PLEN
*     phi is the angle between source and z (shortest route)
      COSPHI = SOURCE(3)
      SINPHI = PLEN
      SROT(1,1) = P1**2 + (1.-P1**2)*COSPHI
      SROT(1,2) = P1*P2*(1.-COSPHI)
      SROT(1,3) = P2*SINPHI
      SROT(2,1) = SROT(1,2)
      SROT(2,2) = P2**2 + (1.-P2**2)*COSPHI
      SROT(2,3) = -P1*SINPHI
      SROT(3,1) = -SROT(1,3)
      SROT(3,2) = -SROT(2,3)
      SROT(3,3) = COSPHI
      SROT(1,4) = 0.0
      SROT(2,4) = 0.0
      SROT(3,4) = 0.0
      SROT(4,1) = 0.0
      SROT(4,2) = 0.0
      SROT(4,3) = 0.0
      SROT(4,4) = 1.0
*
*     Quadrics will require the inverse matrix also (and its transpose)
*     This is also a convenient place to check legality of TMAT
      CALL QSETUP
*
*     Get input mode
      READ (INPUT,*,ERR=104) INMODE
C     WRITE (NOISE,*) 'inmode=',INMODE
      CALL ASSERT (INMODE.GE.1,'bad inmode')
*
*     Get input format(s)
      IF (INMODE.EQ.1.OR.INMODE.EQ.2) THEN
        READ (INPUT,'(A)',ERR=104) INFMT
C       WRITE (NOISE,*) 'infmt=',INFMT
        II = 0
2       CONTINUE
        IF (INFMT(1:1).EQ.' ') THEN
          INFMT(1:79) = INFMT(2:80)
          INFMT(80:80) = ' '
          II = II + 1
          IF (II.LT.80) GO TO 2
        ENDIF
        IF (INFMT(1:1).EQ.'*') THEN
          INFLG = .TRUE.
        ELSE
          INFLG = .FALSE.
        ENDIF
      ELSEIF (INMODE.GE.3) THEN
C       WRITE (NOISE,*) 'infmts:'
        DO 4 I=1,3
          READ (INPUT,'(A)',ERR=104) INFMTS(I)
C         WRITE (NOISE,*) INFMTS(I)
          II = 0
3         CONTINUE
          IF (INFMTS(I)(1:1).EQ.' ') THEN
            INFMTS(I)(1:79) = INFMTS(I)(2:80)
            INFMTS(I)(80:80) = ' '
            II = II + 1
            IF (II.LT.80) GO TO 3
          ENDIF
          IF (INFMTS(I)(1:1).EQ.'*') THEN
            INFLGS(I) = .TRUE.
          ELSE
            INFLGS(I) = .FALSE.
          ENDIF
4       CONTINUE
      INFLGS(PLANE)     = INFLGS(TRIANG)
      INFMTS(PLANE)     = INFMTS(TRIANG)
      INFLGS(NORMS)     = INFLGS(TRIANG)
      INFMTS(NORMS)     = INFMTS(TRIANG)
      INFLGS(VERTEXRGB) = INFLGS(TRIANG)
      INFMTS(VERTEXRGB) = INFMTS(TRIANG)
      INFLGS(VERTRANSP) = INFLGS(TRIANG)
      INFMTS(VERTRANSP) = INFMTS(TRIANG)
      INFLGS(MATERIAL)  = .TRUE.
      INFLGS(GLOWLIGHT) = .TRUE.
      INFLGS(QUADRIC)   = .TRUE.
      ELSE
        CALL ASSERT (.FALSE.,'bad inmode')
      ENDIF
c
c     Done with header records
c     Do we force-close the input file, so that we can borrow headers,
c     or should we keep going as long as the file continues?
c     The following 4 lines implement the former, so that the initial
c     @file command essentially means 'use his header records for me too'.
c
c     IF (INPUT.NE.INPUT0) THEN
c     CLOSE(INPUT)
c     INPUT = INPUT0
c     ENDIF
c     As of V2.5e, however, we keep reading.
c     >>> This is a change! <<<
c
c     End of header processing
*
*     Give them a notice to stare at while the program cranks along
      WRITE (NOISE,'(1X)')
      WRITE (NOISE,*) '%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%',
     &                '%%%%%%%%%%%%%%%%%%%%%%%%%%'
      WRITE (NOISE,*) '%                      Raster3D    ',VERSION,
     &                '                    %'
      WRITE (NOISE,*) '%            -------------------------',
     &                '-------------            %'
      WRITE (NOISE,*) '% If you publish figures generated by this ',
     &                'program please cite %'
      WRITE (NOISE,*) '%   Merritt & Bacon (1997)  ',
     &                'Meth. Enzymol. 277, 505-524.','       %'
      WRITE (NOISE,*) '%            -------------------------',
     &                '-------------            %'
      WRITE (NOISE,*) '%                  Raster3D distribution site',
     &                '                  %'
      WRITE (NOISE,*) '%          http://www.bmsc.washington.edu/',
     &                'raster3d/            %'
      WRITE (NOISE,*) '% comments & suggestions to:  ',
     &                '        merritt@u.washington.edu %'
      WRITE (NOISE,*) '%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%',
     &                '%%%%%%%%%%%%%%%%%%%%%%%%%%'
*
*     If label processing is selected on command line, initialize 
*     PostScript output file.  This isn't needed for Version 3, which
*     uses libgd rather than PostScript to handle labels
      PSSCALE = 2.0 * MIN( NTX*NOX/2., NTY*NOY/2. )
      CALL LSETUP( PSSCALE, BKGND, TITLE )
      WRITE (NOISE,'(1X)')
*
*     Initialize gamma correction table
      GAMMACORRECTION = .FALSE.
      IF (GAMMA.LT.0.99 .OR. GAMMA.GT.1.01) GAMMACORRECTION = .TRUE.
      IF (GAMMACORRECTION) THEN
            DO I=0,MAXRGB
           G = I
           GAMMA_MAP(I+1) = (G/MAXRGB) ** (1.0/GAMMA) * MAXRGB + 0.5
      ENDDO
      ENDIF
*
*
*     Initialize counters
      DO 5 J=1,NTY
      DO 5 I=1,NTX
        KOUNT(I,J) = 0
      TTRANS(I,J) = 0
5     CONTINUE
      DO 6 J=1,NSY
      DO 6 I=1,NSX
        MOUNT(I,J) = 0
6     CONTINUE
c
      DO 662 I = 1,SIZE(FLAG)
      FLAG(I) = 0
  662 CONTINUE
      nprops  = 0
      npropm  = 0
      ntransp = 0
      nsphere = 0
      ncylind = 0
      nplanes = 0
      nhidden = 0
      ninside = 0
      mstate  = 0
      nlabels = 0
      nglows  = 0
      nquads  = 0
      nclip   = 0
      nvtrans = 0
      nbplanes = 0
      tranovfl = 0
c
c     Objects in, and count up objects that may impinge on each tile
      NDET = 0
      IF (SHADOW) MDET = 0
      N = 0
c
c     Read in next object
7     CONTINUE
      IF (INMODE.EQ.1.OR.INMODE.EQ.2) THEN
        INTYPE = INMODE
      GOTO 8
      ENDIF
C
C     READ (INPUT,*,END=50) INTYPE
      READ (INPUT,'(A)',END=50,ERR=47) LINE
c     Nov 1997 - allow # comments
      IF (LINE(1:1) .EQ. '#') THEN
            GOTO 7
c     May 1996 - allow file indirection
      ELSE IF (LINE(1:1) .EQ. '@') THEN
        J = 1
      K = 132
      L = 132
      DO I=132,2,-1
        IF (LINE(I:I).NE.' ') J = I
        IF (LINE(I:I).EQ.'#') K = I-1
        IF (LINE(I:I).EQ.'!') K = I-1
        IF (LINE(I:I).EQ.CHAR(0)) K = I-1
        IF (LINE(I:I).EQ.'    ') LINE(I:I) = ' '
      ENDDO
      IF (J.EQ.1) GOTO 7
      L=0
      DO I=J,K
        IF (LINE(I:I).NE.' ') L = I
      ENDDO
      K = L
      IF (LINE(K:K).eq.'Z' .or. LINE(K-1:K).eq.'gz') THEN
c         ungz will uncompress into a temporary file, which ought to 
c         be deleted later.  Unfortunately, that's hard to do in g77 
c         since it doesn't support dispose='DELETE'.
          line(k+1:k+1) = char(0)
          if (0 .gt. ungz( line(j:k+1), fullname )) goto 73
          j = 1
          k = 132
          do i=132,2,-1
            if (fullname(i:i).eq.' ')     k = i-1
            if (fullname(i:i).eq.char(0)) k = i-1
          enddo
          if (verbose) 
     &            write(noise,*) 'Creating temporary file: ',fullname(j:k)
          open (unit=input+1,err=73,status='OLD',
     &              file=fullname(j:k))
          fullname = line(2:132)
          goto 72
      ENDIF
   70 CONTINUE
      OPEN (UNIT=INPUT+1,ERR=71,STATUS='OLD',
     &            FILE=LINE(J:K))
      FULLNAME = LINE(J:K)
      GOTO 72
   71 CONTINUE
      IF (LINE(K-3:K).NE.'.r3d') THEN
          K = K + 4
          LINE(K-3:K) = '.r3d'
          GOTO 70
      ENDIF
      CALL LIBLOOKUP( LINE(J:K), FULLNAME )
      OPEN (UNIT=INPUT+1,ERR=73,STATUS='OLD',FILE=FULLNAME)
   72 CONTINUE
      DO I=132,2,-1
        IF (FULLNAME(I:I).EQ.' ') J = I
      ENDDO
      WRITE (NOISE,'(A,A)') '  + Opening input file ',FULLNAME(1:J)
      INPUT = INPUT + 1
      CALL ASSERT(INPUT-INPUT0.LE.MAXLEV, 
     &                  'Too many levels of indirection')
      GOTO 7
   73 WRITE (NOISE,'(A,A)') ' >> Cannot open file ',LINE(J:K)
      GOTO 7
      ELSE
      READ (LINE,*,END=50,ERR=74) INTYPE
      GOTO 76
   74 WRITE (NOISE,'(A,A)') ' >> Unrecognized line: ',LINE
      GOTO 7
   76 CONTINUE
      ENDIF
      IF (INTYPE.EQ.0) GO TO 50
      IF (INTYPE .EQ. CYLFLAT) THEN
      INTYPE = CYLIND
      FLAG(N+1) = ior( FLAG(N+1), FLAT )
      ELSEIF (INTYPE .EQ. MATEND) THEN
      MSTATE    = 0
      CLRITY    = 0
      CLROPT    = 0
      TRNSPOPT  = 0
      MATCOL    = .FALSE.
      ISOLATION = 0
      CLIPPING  = .FALSE.
      NBOUNDS   = 0
      ORTEPLIKE = .FALSE.
      GOTO 7
      ELSEIF (INTYPE .EQ. FONT) THEN
      IF (LFLAG) THEN
        CALL LINP( INPUT, INTYPE, .FALSE., RGBMAT )
      ELSE
          READ (INPUT,'(A)',END=50) LINE
      ENDIF
      GOTO 7
      ELSEIF (INTYPE .EQ. LABEL) THEN
      NLABELS = NLABELS + 1
      IF (LFLAG) THEN
        IF (MSTATE.EQ.MATERIAL .AND. MATCOL) THEN
          CALL LINP( INPUT, INTYPE, .TRUE., RGBMAT )
        ELSE
          CALL LINP( INPUT, INTYPE, .FALSE., RGBMAT )
        ENDIF
      ELSE
          READ (INPUT,'(A)',END=50) LINE
          READ (INPUT,'(A)',END=50) LINE
      ENDIF
      GOTO 7
      ELSEIF (INTYPE .EQ. ISOLATE1) THEN
      ISOLATION = 1
      GOTO 7
      ELSEIF (INTYPE .EQ. ISOLATE2) THEN
      ISOLATION = 2
      GOTO 7
c     Global Properties
      ELSEIF (INTYPE .EQ. GPROP) THEN
      READ (INPUT,'(A)',END=50) LINE
      L = 1
      DO I = 132, 1, -1
        IF (LINE(I:I).NE.' '.AND.LINE(I:I).NE.' ') L = I
      ENDDO
      IF (LINE(L:L+2).EQ.'FOG') THEN
        READ (LINE(L+4:74),*,ERR=771,END=771) 
     &          FOGTYPE, FOGFRONT, FOGBACK, FOGDEN
  771   CONTINUE
        FOGRGB(1) = BKGND(1)
        FOGRGB(2) = BKGND(2)
        FOGRGB(3) = BKGND(3)
        CALL CHKRGB(FOGRGB(1),FOGRGB(2),FOGRGB(3),'invalid fog color')
        IF (FOGTYPE.NE.1)  FOGTYPE = 0
        IF (FOGDEN.LE.0.0) FOGDEN = 0.5
      ELSE IF (LINE(L:L+8).EQ.'FRONTCLIP') THEN
        READ (LINE(L+10:132),*,ERR=75) FRONTCLIP
        FRONTCLIP = FRONTCLIP * SCALE / TMAT(4,4)
      ELSE IF (LINE(L:L+7).EQ.'BACKCLIP') THEN
        READ (LINE(L+9:132),*,ERR=75) BACKCLIP
        BACKCLIP = BACKCLIP * SCALE / TMAT(4,4)
      ELSE IF (LINE(L:L+7).EQ.'ROTATION') THEN
        READ (LINE(L+9:132),*,ERR=773,END=773) 
     &            ((RAFTER(I,J),J=1,3),I=1,3)
        GOTO 774
  773   READ (INPUT,*,ERR=75) ((RAFTER(I,J),J=1,3),I=1,3)
  774   WRITE (NOISE,775) ((RAFTER(I,J),J=1,3),I=1,3)
  775   FORMAT('Post-rotation matrix:  ',3(/,3F10.4))
        D = DET(RAFTER)
        IF (ABS(1.0-ABS(D)).GT.0.02) WRITE (NOISE,*) 
     &          '>>> Warning: Post-rotation matrix has determinant',D
        IF (INVERT) THEN
          RAFTER(1,2) = -RAFTER(1,2)
          RAFTER(2,1) = -RAFTER(2,1)
          RAFTER(2,3) = -RAFTER(2,3)
          RAFTER(3,2) = -RAFTER(3,2)
        ENDIF
        CALL QSETUP
      ELSE IF (LINE(L:L+10).EQ.'TRANSLATION') THEN
        READ (LINE(L+12:132),*,ERR=776,END=776) 
     &            (TAFTER(I),I=1,3)
        GOTO 777
  776   READ (INPUT,*,ERR=75) (TAFTER(I),I=1,3)
  777   WRITE (NOISE,778) (TAFTER(I),I=1,3)
  778   FORMAT('Post-translation:      ',1(/,3F10.4))
        IF (INVERT) TAFTER(2) = -TAFTER(2)
      ELSE IF (LINE(L:L+4).EQ.'DUMMY') THEN
      ELSE
        GOTO 75
      ENDIF
      GOTO 7

   75 CONTINUE
      WRITE(NOISE,'(A,A)') 
     &         '>> Unrecognized or incomplete GPROP option ',LINE
      GOTO 7
*
      ENDIF
*
COLD  CALL ASSERT (INTYPE.GE.1.AND.INTYPE.LE.MXTYPE,'bad object')
      IF (INTYPE.LT.1 .OR. INTYPE.GT.MXTYPE) GOTO 47
      CALL ASSERT (INTYPE.NE.INTERNAL,'object type 4 not available')
c
c     Read in object details, now we know what kind it is.
c     Allow an all-zeroes record to terminate input for the
c     benefit of those of us who might inadvertently supply
c     a series of blank records after our real input as a
c     side-effect of tape blocking or sloppiness or ...
8     CONTINUE
      IF (INMODE.GE.3) THEN
      INFMT = INFMTS(INTYPE)
      INFLG = INFLGS(INTYPE)
      ENDIF
      IF (INFLG) THEN
        READ (INPUT,*,END=50) (BUF(I),I=1,IDET(INTYPE))
      ELSE
        READ (INPUT,INFMT,END=50) (BUF(I),I=1,IDET(INTYPE))
      ENDIF
c     15-Dec-1999 This was supposed to check for all-zero line and exit
c                 but all zeros is legal for [at least!] LABELs
      IF (INTYPE.EQ.LABEL) GOTO 9
      DO I=1,IDET(INTYPE)
        IF (BUF(I).NE.0.) GO TO 9
      ENDDO
      GO TO 50
9     CONTINUE
*
*     Expand array DETAIL if needed: NDET+KDET(INTYPE) > size(DETAIL)
      NEEDMEM = NDET+KDET(INTYPE) 
      if ( NEEDMEM .GT. SIZE(DETAIL) ) THEN
          CALL GET_TRY(SIZE(DETAIL), NEEDMEM, TRY1, 1)
          DO 920 ITRY = 1,3
              if (TRY1(ITRY).LE.0 .OR. TRY1(ITRY).GT.MAXMEM) GOTO 920
              ALLOCATE( TMP1DR(TRY1(ITRY)), stat=ierr )
              if (ierr .NE. 0) GOTO 920
              TMP1DR = 0.
              TMP1DR = DETAIL
              if(TEST_ALLOC)write(NOISE,*)"Expand DETAIL to ",try1(ITRY)
              CALL MOVE_ALLOC(from=TMP1DR, to=DETAIL)
              GOTO 922
920       CONTINUE
      ENDIF
922   CALL ASSERT (NEEDMEM.LE.SIZE(DETAIL),
     & 'too many object details - increase MAXDET and recompile')
      IF (SHADOW) THEN
*         Expand array SDTAIL if needed
          NEEDMEM = MDET+SDET(INTYPE)
          IF ( NEEDMEM .GT. SIZE(SDTAIL) ) THEN
              CALL GET_TRY(SIZE(SDTAIL), NEEDMEM, TRY1, 1)
              DO 925 ITRY = 1,3
                  if (TRY1(ITRY).LE.0.OR.TRY1(ITRY).GT.MAXMEM) GOTO 925
                  ALLOCATE( TMP1DR(TRY1(ITRY)), stat=ierr )
                  if (ierr .NE. 0) GOTO 925
                  TMP1DR = 0.
                  TMP1DR = SDTAIL
                  CALL MOVE_ALLOC(from=TMP1DR, to=SDTAIL)
                  if(TEST_ALLOC)write(NOISE,*)"Expand SDTAIL to ",
     &                try1(ITRY)
                  GOTO 927
925           CONTINUE
          ENDIF
927       CALL ASSERT (NEEDMEM.LE.SIZE(SDTAIL),
     & 'too many shadow object details - increase MAXSDT and recompile')
      ENDIF
      N = N + 1
*     Expand arrays ZTEMP, ZINDEX, LIST, MIST, TYPE, FLAG if needed
      IF (N.GT.SIZE(FLAG)) THEN
          CALL GET_TRY(SIZE(FLAG), N, TRY1, 1)
          DO 930 ITRY = 1,3
              if (TRY1(ITRY).LE.0.OR.TRY1(ITRY).GT.MAXMEM) GOTO 930
              ALLOCATE( TMP1DR(TRY1(ITRY)), stat=ierr )
              if (ierr .NE. 0) GOTO 930
              TMP1DR = 0.
              TMP1DR = ZTEMP
              CALL MOVE_ALLOC(from=TMP1DR, to=ZTEMP)
              ALLOCATE( TMP1D(TRY1(ITRY)), stat=ierr )
              if (ierr .NE. 0) GOTO 930
              TMP1D = 0
              TMP1D = ZINDEX
              CALL MOVE_ALLOC(from=TMP1D, to=ZINDEX)
              ALLOCATE( TMP1D(TRY1(ITRY)), stat=ierr )
              if (ierr .NE. 0) GOTO 930
              TMP1D = 0
              TMP1D = LIST
              CALL MOVE_ALLOC(from=TMP1D, to=LIST)
              ALLOCATE( TMP1D(TRY1(ITRY)), stat=ierr )
              if (ierr .NE. 0) GOTO 930
              TMP1D = 0
              TMP1D = MIST
              CALL MOVE_ALLOC(from=TMP1D, to=MIST)
              ALLOCATE( TMP1D(TRY1(ITRY)), stat=ierr )
              if (ierr .NE. 0) GOTO 930
              TMP1D = 0
              TMP1D = TYPE
              CALL MOVE_ALLOC(from=TMP1D, to=TYPE)
              ALLOCATE( TMP1DI4(TRY1(ITRY)), stat=ierr )
              if (ierr .NE. 0) GOTO 930
              TMP1DI4 = 0
              TMP1DI4 = FLAG
              CALL MOVE_ALLOC(from=TMP1DI4, to=FLAG)
              if(TEST_ALLOC)write(NOISE,*)"Expand MAXOBJ to ",size(FLAG)
              GOTO 932
930       CONTINUE
      ENDIF

932   CALL ASSERT (N.LE.SIZE(FLAG),
     & 'too many objects - increase MAXOBJ and recompile')
C     20-Feb-1997 Save both object type and material type 
      TYPE(N) = INTYPE
      IF (MSTATE.EQ.MATERIAL) FLAG(N) = FLAG(N) + 65536 * NPROPM
      LIST(N) = NDET
      IF (SHADOW) MIST(N) = MDET
      ISTRANS  = 0
*     From this point on, we'll use the symbolic codes for objects
      IF (INTYPE.EQ.TRIANG .or. INTYPE.EQ.PLANE) THEN
*       triangle as read in
        X1A = BUF(1)
        Y1A = BUF(2)
        Z1A = BUF(3)
        X2A = BUF(4)
        Y2A = BUF(5)
        Z2A = BUF(6)
        X3A = BUF(7)
        Y3A = BUF(8)
        Z3A = BUF(9)
        RED = BUF(10)
        GRN = BUF(11)
        BLU = BUF(12)
      CALL CHKRGB( RED,GRN,BLU,'invalid triangle color' )
        CALL ASSERT (IDET(INTYPE).EQ.12,'idet(1).ne.12')
      IF (MSTATE.EQ.MATERIAL) THEN
        FLAG(N) = ior(FLAG(N),PROPS)
        NPROPS  = NPROPS + 1
        IF (CLRITY.GT.0) THEN
          FLAG(N) = ior(FLAG(N),TRANSP)
          IF (CLROPT.EQ.1) FLAG(N) = ior(FLAG(N),MOPT1)
          NTRANSP = NTRANSP + 1
          ISTRANS = 1
        ENDIF
        IF (CLIPPING) THEN
          FLAG(N) = ior(FLAG(N),CLIPPED)
        ENDIF
        IF (NBOUNDS.GT.0) FLAG(N) = ior(FLAG(N),BOUNDED)
        IF (MATCOL) THEN
          RED = RGBMAT(1)
          GRN = RGBMAT(2)
          BLU = RGBMAT(3)
        ENDIF
      ENDIF
*     Isolated objects not transformed by TMAT, but still subject to inversion
      IF (ISOLATION.GT.0) THEN
        CALL ISOLATE(X1A,Y1A)
        CALL ISOLATE(X2A,Y2A)
        CALL ISOLATE(X3A,Y3A)
      ELSE
*     update true coordinate limits
        TRULIM(1,1) = MIN( TRULIM(1,1), X1A,X2A,X3A)
        TRULIM(1,2) = MAX( TRULIM(1,2), X1A,X2A,X3A)
        TRULIM(2,1) = MIN( TRULIM(2,1), Y1A,Y2A,Y3A)
        TRULIM(2,2) = MAX( TRULIM(2,2), Y1A,Y2A,Y3A)
        TRULIM(3,1) = MIN( TRULIM(3,1), Z1A,Z2A,Z3A)
        TRULIM(3,2) = MAX( TRULIM(3,2), Z1A,Z2A,Z3A)
*       modify the input, so to speak
        CALL TRANSF (X1A,Y1A,Z1A)
        CALL TRANSF (X2A,Y2A,Z2A)
        CALL TRANSF (X3A,Y3A,Z3A)
      ENDIF
*       perspective factor for each corner
      IF (EYEPOS.GT.0) THEN
          PFAC1 = PERSP( Z1A )
          PFAC2 = PERSP( Z2A )
          PFAC3 = PERSP( Z3A )
      END IF
*       apply perspective
        X1B = X1A * PFAC1
        Y1B = Y1A * PFAC1
        Z1B = Z1A * PFAC1
        X2B = X2A * PFAC2
        Y2B = Y2A * PFAC2
        Z2B = Z2A * PFAC2
        X3B = X3A * PFAC3
        Y3B = Y3A * PFAC3
        Z3B = Z3A * PFAC3
*       scale and translate to pixel space
        X1C = X1B * SCALE + XCENT
        Y1C = Y1B * SCALE + YCENT
        Z1C = Z1B * SCALE
        X2C = X2B * SCALE + XCENT
        Y2C = Y2B * SCALE + YCENT
        Z2C = Z2B * SCALE
        X3C = X3B * SCALE + XCENT
        Y3C = Y3B * SCALE + YCENT
        Z3C = Z3B * SCALE
*     save transformed Z limits
      ZLIM(1) = MIN( ZLIM(1), Z1C,Z2C,Z3C )
      ZLIM(2) = MAX( ZLIM(2), Z1C,Z2C,Z3C )
*
*     check for Z-clipping
      JUSTCLIPPED = .FALSE.
      IF (INTYPE.NE.PLANE) THEN
        IF ((MIN(Z1C,Z2C,Z3C) .GT. FRONTCLIP)
     &    .OR.(MAX(Z1C,Z2C,Z3C) .LT. BACKCLIP)) THEN
            JUSTCLIPPED = .TRUE.
            GOTO 45
        ENDIF
        IF (CLIPPING) THEN
          MIND = LIST(MLIST(NPROPM))
          IF ((MIN(Z1C,Z2C,Z3C) .GT. DETAIL(MIND+16))
     &      .OR.(MAX(Z1C,Z2C,Z3C) .LT. DETAIL(MIND+17))) THEN
            JUSTCLIPPED = .TRUE.
            GOTO 45
          ENDIF
        ENDIF
      ENDIF
*
*       solve for coefficients of plane eqn z=Ax+By+C
        CALL PLANER(X1C,Y1C,Z1C,
     &             X2C,Y2C,Z2C,
     &             X3C,Y3C,Z3C, A,B,C,D)
*       save results for PLANE object
*     PLANE impinges on all tiles, but casts no shadows
      IF (INTYPE.EQ. PLANE) THEN
          NPLANES = NPLANES + 1
          DETAIL(NDET+1) = A
          DETAIL(NDET+2) = B
          DETAIL(NDET+3) = C
          DETAIL(NDET+4) = D
          DETAIL(NDET+5) = RED
          DETAIL(NDET+6) = GRN
          DETAIL(NDET+7) = BLU
          CALL ASSERT(KDET(INTYPE).EQ.7,'kdet(6).ne.7')
          NDET = NDET + KDET(INTYPE)
          DO IY = 1,NTY
          DO IX = 1,NTX
            KOUNT(IX,IY) = KOUNT(IX,IY) + 1
            TTRANS(IX,IY) = TTRANS(IX,IY) + ISTRANS
          ENDDO
          ENDDO
C           write(NOISE,*)"Incr kount ",NTX," x ",NTY
          IF (SHADOW) THEN
              MDET = MDET + SDET(INTYPE)
          ENDIF
          GOTO 7
      ENDIF
*       save results for normal triangles
        DETAIL(NDET+1) = X1C
        DETAIL(NDET+2) = Y1C
        DETAIL(NDET+3) = Z1C
        DETAIL(NDET+4) = X2C
        DETAIL(NDET+5) = Y2C
        DETAIL(NDET+6) = Z2C
        DETAIL(NDET+7) = X3C
        DETAIL(NDET+8) = Y3C
        DETAIL(NDET+9) = Z3C
        DETAIL(NDET+10) = A
        DETAIL(NDET+11) = B
        DETAIL(NDET+12) = C
        DETAIL(NDET+13) = D
        DETAIL(NDET+14) = RED
        DETAIL(NDET+15) = GRN
        DETAIL(NDET+16) = BLU
        CALL ASSERT (KDET(INTYPE).EQ.16,'kdet(1).ne.16')
        NDET = NDET + KDET(INTYPE)
*       tally for tiles the object might impinge on
        IXLO = MIN(X1C,X2C,X3C)/NPX + 1
        IXHI = MAX(X1C,X2C,X3C)/NPX + 1
        IYLO = MIN(Y1C,Y2C,Y3C)/NPY + 1
        IYHI = MAX(Y1C,Y2C,Y3C)/NPY + 1
        IF (IXLO.LT.1  ) IXLO=1
        IF (IXLO.GT.NTX) GO TO 11
        IF (IXHI.LT.1  ) GO TO 11
        IF (IXHI.GT.NTX) IXHI=NTX
        IF (IYLO.LT.1  ) IYLO=1
        IF (IYLO.GT.NTY) GO TO 11
        IF (IYHI.LT.1  ) GO TO 11
        IF (IYHI.GT.NTY) IYHI=NTY
        DO 10 IY=IYLO,IYHI
        DO 10 IX=IXLO,IXHI
          KOUNT(IX,IY) = KOUNT(IX,IY) + 1
        TTRANS(IX,IY) = TTRANS(IX,IY) + ISTRANS
10      CONTINUE
C       write(NOISE,*)"Incr kount b ",IXLO,"-",IXHI," x ",IYLO,"-",IYHI
11      CONTINUE
*       repeat for shadow buffer if necessary
        IF (SHADOW) THEN
*         rotate light source to z to take light source viewpoint
          X1R = SROT(1,1)*X1B+SROT(1,2)*Y1B+SROT(1,3)*Z1B
          Y1R = SROT(2,1)*X1B+SROT(2,2)*Y1B+SROT(2,3)*Z1B
          Z1R = SROT(3,1)*X1B+SROT(3,2)*Y1B+SROT(3,3)*Z1B
          X2R = SROT(1,1)*X2B+SROT(1,2)*Y2B+SROT(1,3)*Z2B
          Y2R = SROT(2,1)*X2B+SROT(2,2)*Y2B+SROT(2,3)*Z2B
          Z2R = SROT(3,1)*X2B+SROT(3,2)*Y2B+SROT(3,3)*Z2B
          X3R = SROT(1,1)*X3B+SROT(1,2)*Y3B+SROT(1,3)*Z3B
          Y3R = SROT(2,1)*X3B+SROT(2,2)*Y3B+SROT(2,3)*Z3B
          Z3R = SROT(3,1)*X3B+SROT(3,2)*Y3B+SROT(3,3)*Z3B
*         scale and translate for shadow space
          X1S = X1R * SCALE + SXCENT
          Y1S = Y1R * SCALE + SYCENT
          Z1S = Z1R * SCALE
          X2S = X2R * SCALE + SXCENT
          Y2S = Y2R * SCALE + SYCENT
          Z2S = Z2R * SCALE
          X3S = X3R * SCALE + SXCENT
          Y3S = Y3R * SCALE + SYCENT
          Z3S = Z3R * SCALE
*         solve plane eqn etc.
          CALL PLANER(X1S,Y1S,Z1S,
     &               X2S,Y2S,Z2S,
     &               X3S,Y3S,Z3S, A,B,C,D)
*         save results etc.
          SDTAIL(MDET+1) = X1S
          SDTAIL(MDET+2) = Y1S
          SDTAIL(MDET+3) = Z1S
          SDTAIL(MDET+4) = X2S
          SDTAIL(MDET+5) = Y2S
          SDTAIL(MDET+6) = Z2S
          SDTAIL(MDET+7) = X3S
          SDTAIL(MDET+8) = Y3S
          SDTAIL(MDET+9) = Z3S
          SDTAIL(MDET+10) = A
          SDTAIL(MDET+11) = B
          SDTAIL(MDET+12) = C
          SDTAIL(MDET+13) = D
          CALL ASSERT (SDET(INTYPE).EQ.13,'sdet(1).ne.13')
          MDET = MDET + SDET(INTYPE)
*         tally for shadow tiles the object might impinge on
          IXLO = MIN(X1S,X2S,X3S)/NPX + 1
          IXHI = MAX(X1S,X2S,X3S)/NPX + 1
          IYLO = MIN(Y1S,Y2S,Y3S)/NPY + 1
          IYHI = MAX(Y1S,Y2S,Y3S)/NPY + 1
          IF (IXLO.LT.1  ) IXLO=1
          IF (IXLO.GT.NSX) GO TO 16
          IF (IXHI.LT.1  ) GO TO 16
          IF (IXHI.GT.NSX) IXHI=NSX
          IF (IYLO.LT.1  ) IYLO=1
          IF (IYLO.GT.NSY) GO TO 16
          IF (IYHI.LT.1  ) GO TO 16
          IF (IYHI.GT.NSY) IYHI=NSY
          DO 15 IY=IYLO,IYHI
          DO 15 IX=IXLO,IXHI
            MOUNT(IX,IY) = MOUNT(IX,IY) + 1
15        CONTINUE
16        CONTINUE
        ENDIF
      ELSEIF (INTYPE.EQ.SPHERE) THEN
*       sphere as read in
        XA = BUF(1)
        YA = BUF(2)
        ZA = BUF(3)
        RA = BUF(4)
        RED = BUF(5)
        GRN = BUF(6)
        BLU = BUF(7)
      CALL CHKRGB (RED,GRN,BLU,'invalid sphere color')
        CALL ASSERT (IDET(INTYPE).EQ.7,'idet(2).ne.7')
      IF (MSTATE.EQ.MATERIAL) THEN
        FLAG(N) = ior(FLAG(N),PROPS)
        NPROPS  = NPROPS + 1
        IF (CLRITY.GT.0) THEN
          FLAG(N) = ior(FLAG(N),TRANSP)
          IF (CLROPT.EQ.1) FLAG(N) = ior(FLAG(N),MOPT1)
          NTRANSP = NTRANSP + 1
          ISTRANS = 1
        ENDIF
        IF (CLIPPING) THEN
          FLAG(N) = ior(FLAG(N),CLIPPED)
        ENDIF
        IF (NBOUNDS.GT.0) FLAG(N) = ior(FLAG(N),BOUNDED)
        IF (MATCOL) THEN
          RED = RGBMAT(1)
          GRN = RGBMAT(2)
          BLU = RGBMAT(3)
        ENDIF
      ENDIF
*     Isolated objects not transformed by TMAT, but still subject to inversion
      IF (ISOLATION.GT.0) THEN
        CALL ISOLATE(XA,YA)
      ELSE
*     update true coordinate limits
        TRULIM(1,1) = MIN( TRULIM(1,1), XA )
        TRULIM(1,2) = MAX( TRULIM(1,2), XA )
        TRULIM(2,1) = MIN( TRULIM(2,1), YA )
        TRULIM(2,2) = MAX( TRULIM(2,2), YA )
        TRULIM(3,1) = MIN( TRULIM(3,1), ZA )
        TRULIM(3,2) = MAX( TRULIM(3,2), ZA )
*       modify the input, as it were
          CALL TRANSF (XA,YA,ZA)
          RA = RA / TMAT(4,4)
      ENDIF
*       perspective
      IF (EYEPOS.GT.0) PFAC = PERSP(ZA)
        XB = XA * PFAC
        YB = YA * PFAC
        ZB = ZA * PFAC
        RB = RA * PFAC
*       scale & translate
        XC = XB * SCALE + XCENT
        YC = YB * SCALE + YCENT
        ZC = ZB * SCALE
        RC = RB * SCALE
*     save transformed Z limits
      ZLIM(1) = MIN( ZLIM(1), ZC )
      ZLIM(2) = MAX( ZLIM(2), ZC )
*     check for Z-clipping
      IF (ZC .GT. FRONTCLIP .OR. ZC .LT. BACKCLIP) THEN
          JUSTCLIPPED = .TRUE.
          GOTO 45
      ELSE
          JUSTCLIPPED = .FALSE.
      ENDIF
*       save results
        DETAIL(NDET+1) = XC
        DETAIL(NDET+2) = YC
        DETAIL(NDET+3) = ZC
        DETAIL(NDET+4) = RC
        DETAIL(NDET+5) = RED
        DETAIL(NDET+6) = GRN
        DETAIL(NDET+7) = BLU
        CALL ASSERT (KDET(INTYPE).EQ.7,'kdet(2).ne.7')
        NDET = NDET + KDET(INTYPE)
      nsphere = nsphere + 1
*       tally for tiles the object might impinge on
        IXLO = (XC-RC)/NPX + 1
        IXHI = (XC+RC)/NPX + 1
        IYLO = (YC-RC)/NPY + 1
        IYHI = (YC+RC)/NPY + 1
        IF (IXLO.LT.1  ) IXLO=1
        IF (IXLO.GT.NTX) GO TO 21
        IF (IXHI.LT.1  ) GO TO 21
        IF (IXHI.GT.NTX) IXHI=NTX
        IF (IYLO.LT.1  ) IYLO=1
        IF (IYLO.GT.NTY) GO TO 21
        IF (IYHI.LT.1  ) GO TO 21
        IF (IYHI.GT.NTY) IYHI=NTY
        DO 20 IY=IYLO,IYHI
        DO 20 IX=IXLO,IXHI
          KOUNT(IX,IY) = KOUNT(IX,IY) + 1
        TTRANS(IX,IY) = TTRANS(IX,IY) + ISTRANS
20      CONTINUE
C       write(NOISE,*)"Incr kount c ",IXLO,"-",IXHI," x ",IYLO,"-",IYHI,
C    &  kstart(IXLO,IYLO),kstop(IXLO,IYLO),kount(IXLO,IYLO)
21      CONTINUE
*       repeat for shadow buffer if necessary
        IF (SHADOW) THEN
*         rotate light source to z to take light source viewpoint
          XR = SROT(1,1)*XB+SROT(1,2)*YB+SROT(1,3)*ZB
          YR = SROT(2,1)*XB+SROT(2,2)*YB+SROT(2,3)*ZB
          ZR = SROT(3,1)*XB+SROT(3,2)*YB+SROT(3,3)*ZB
          RR = RB
*         scale and translate for shadow space
          XS = XR * SCALE + SXCENT
          YS = YR * SCALE + SYCENT
          ZS = ZR * SCALE
          RS = RR * SCALE
*         save results
          SDTAIL(MDET+1) = XS
          SDTAIL(MDET+2) = YS
          SDTAIL(MDET+3) = ZS
          SDTAIL(MDET+4) = RS
          CALL ASSERT (SDET(INTYPE).EQ.4,'sdet(2).ne.4')
          MDET = MDET + SDET(INTYPE)
*         tally for shadow tiles the object might impinge on
          IXLO = (XS-RS)/NPX + 1
          IXHI = (XS+RS)/NPX + 1
          IYLO = (YS-RS)/NPY + 1
          IYHI = (YS+RS)/NPY + 1
          IF (IXLO.LT.1  ) IXLO=1
          IF (IXLO.GT.NSX) GO TO 26
          IF (IXHI.LT.1  ) GO TO 26
          IF (IXHI.GT.NSX) IXHI=NSX
          IF (IYLO.LT.1  ) IYLO=1
          IF (IYLO.GT.NSY) GO TO 26
          IF (IYHI.LT.1  ) GO TO 26
          IF (IYHI.GT.NSY) IYHI=NSY
          DO 25 IY=IYLO,IYHI
          DO 25 IX=IXLO,IXHI
            MOUNT(IX,IY) = MOUNT(IX,IY) + 1
25        CONTINUE
26        CONTINUE
        ENDIF

c     30-Dec-99 duplicate transparent spheres, if requested, so that
c     the inside surface can be rendered also. BUF() is still loaded
c     with specs for the current object; we just need to set flags.
      IF (  iand(FLAG(N),TRANSP).NE.0 .AND. CLROPT.EQ.2
     &  .AND. iand(FLAG(N),INSIDE).EQ.0) THEN
          FLAG(N+1) = INSIDE
          NINSIDE   = NINSIDE + 1
          GOTO 9
      ENDIF

      ELSEIF (INTYPE.EQ.CYLIND) THEN
*     EAM May 1990 cylinder as read in
        X1A = BUF(1)
        Y1A = BUF(2)
        Z1A = BUF(3)
      R1A = BUF(4)
        X2A = BUF(5)
        Y2A = BUF(6)
        Z2A = BUF(7)
      R2A = R1A
        RED = BUF(9)
        GRN = BUF(10)
        BLU = BUF(11)
      CALL CHKRGB (RED,GRN,BLU,'invalid cylinder color')
        CALL ASSERT (IDET(INTYPE).EQ.11,'idet(1).ne.11')
*     Zero length cylinder is better treated as sphere
*     EAM 22-Nov-96
      IF ((iand(FLAG(N),FLAT).EQ.0) .AND.
     &          (X1A.EQ.X2A).AND.(Y1A.EQ.Y2A).AND.(Z1A.EQ.Z2A)) THEN
            BUF(5) = BUF(9)
            BUF(6) = BUF(10)
            BUF(7) = BUF(11)
            INTYPE = SPHERE
            N = N-1
            GOTO 9
      ENDIF
      IF (MSTATE.EQ.MATERIAL) THEN
        FLAG(N) = ior(FLAG(N),PROPS)
        NPROPS  = NPROPS + 1
        IF (CLRITY.GT.0) THEN
          FLAG(N) = ior(FLAG(N),TRANSP)
          IF (CLROPT.EQ.1) FLAG(N) = ior(FLAG(N),MOPT1)
          NTRANSP = NTRANSP + 1
          ISTRANS = 1
        ENDIF
        IF (CLIPPING) THEN
          FLAG(N) = ior(FLAG(N),CLIPPED)
        ENDIF
        IF (NBOUNDS.GT.0) FLAG(N) = ior(FLAG(N),BOUNDED)
        IF (MATCOL) THEN
          RED = RGBMAT(1)
          GRN = RGBMAT(2)
          BLU = RGBMAT(3)
        ENDIF
      ENDIF
*     Isolated objects not transformed by TMAT, but still subject to inversion
      IF (ISOLATION.GT.0) THEN
        CALL ISOLATE(X1A,Y1A)
        CALL ISOLATE(X2A,Y2A)
      ELSE
*     update true coordinate limits
        TRULIM(1,1) = MIN( TRULIM(1,1), X1A,X2A)
        TRULIM(1,2) = MAX( TRULIM(1,2), X1A,X2A)
        TRULIM(2,1) = MIN( TRULIM(2,1), Y1A,Y2A)
        TRULIM(2,2) = MAX( TRULIM(2,2), Y1A,Y2A)
        TRULIM(3,1) = MIN( TRULIM(3,1), Z1A,Z2A)
        TRULIM(3,2) = MAX( TRULIM(3,2), Z1A,Z2A)
*       modify the input, so to speak
          CALL TRANSF (X1A,Y1A,Z1A)
          CALL TRANSF (X2A,Y2A,Z2A)
          R1A = R1A / TMAT(4,4)
          R2A = R2A / TMAT(4,4)
      ENDIF
*       perspective factor for each corner
      IF (EYEPOS.GT.0) THEN
          PFAC1 = PERSP( Z1A )
          PFAC2 = PERSP( Z2A )
      END IF
*       apply perspective
        X1B = X1A * PFAC1
        Y1B = Y1A * PFAC1
        Z1B = Z1A * PFAC1
        R1B = R1A * PFAC1
        X2B = X2A * PFAC2
        Y2B = Y2A * PFAC2
        Z2B = Z2A * PFAC2
        R2B = R2A * PFAC2
*       scale and translate to pixel space
        X1C = X1B * SCALE + XCENT
        Y1C = Y1B * SCALE + YCENT
        Z1C = Z1B * SCALE
        R1C = R1B * SCALE
        X2C = X2B * SCALE + XCENT
        Y2C = Y2B * SCALE + YCENT
        Z2C = Z2B * SCALE
        R2C = R2B * SCALE
*     save transformed Z limits
      ZLIM(1) = MIN( ZLIM(1), Z1C,Z2C )
      ZLIM(2) = MAX( ZLIM(2), Z1C,Z2C )
*     check for Z-clipping
      IF ((MIN(Z1C,Z2C) .GT. FRONTCLIP)
     &  .OR.(MAX(Z1C,Z2C) .LT. BACKCLIP)) THEN
          JUSTCLIPPED = .TRUE.
          GOTO 45
      ELSE
          JUSTCLIPPED = .FALSE.
      ENDIF
*       save results
        DETAIL(NDET+1) = X1C
        DETAIL(NDET+2) = Y1C
        DETAIL(NDET+3) = Z1C
        DETAIL(NDET+4) = R1C
        DETAIL(NDET+5) = X2C
        DETAIL(NDET+6) = Y2C
        DETAIL(NDET+7) = Z2C
      DETAIL(NDET+8) = R2C
* EAM   save anything else?
        DETAIL(NDET+9)  = RED
        DETAIL(NDET+10) = GRN
        DETAIL(NDET+11) = BLU
        CALL ASSERT (KDET(INTYPE).EQ.11,'kdet(1).ne.11')
        NDET = NDET + KDET(INTYPE)
      ncylind = ncylind + 1
*       tally for tiles the object might impinge on
        IXLO = MIN(X1C-R1C,X2C-R2C)  / NPX + 1
        IXHI = MAX(X1C+R1C,X2C+R2C)  / NPX + 1
        IYLO = MIN(Y1C-R1C,Y2C-R2C)  / NPY + 1
        IYHI = MAX(Y1C+R1C,Y2C+R2C)  / NPY + 1
        IF (IXLO.LT.1  ) IXLO=1
        IF (IXLO.GT.NTX) GO TO 711
        IF (IXHI.LT.1  ) GO TO 711
        IF (IXHI.GT.NTX) IXHI=NTX
        IF (IYLO.LT.1  ) IYLO=1
        IF (IYLO.GT.NTY) GO TO 711
        IF (IYHI.LT.1  ) GO TO 711
        IF (IYHI.GT.NTY) IYHI=NTY
        DO 710 IY=IYLO,IYHI
        DO 710 IX=IXLO,IXHI
          KOUNT(IX,IY) = KOUNT(IX,IY) + 1
        TTRANS(IX,IY) = TTRANS(IX,IY) + ISTRANS
710     CONTINUE
C       write(NOISE,*)"Incr kount d ",IXLO,"-",IXHI," x ",IYLO,"-",IYHI
711     CONTINUE
*       repeat for shadow buffer if necessary
        IF (SHADOW) THEN
*         rotate light source to z to take light source viewpoint
          X1R = SROT(1,1)*X1B+SROT(1,2)*Y1B+SROT(1,3)*Z1B
          Y1R = SROT(2,1)*X1B+SROT(2,2)*Y1B+SROT(2,3)*Z1B
          Z1R = SROT(3,1)*X1B+SROT(3,2)*Y1B+SROT(3,3)*Z1B
          X2R = SROT(1,1)*X2B+SROT(1,2)*Y2B+SROT(1,3)*Z2B
          Y2R = SROT(2,1)*X2B+SROT(2,2)*Y2B+SROT(2,3)*Z2B
          Z2R = SROT(3,1)*X2B+SROT(3,2)*Y2B+SROT(3,3)*Z2B
*         scale and translate for shadow space
          X1S = X1R * SCALE + SXCENT
          Y1S = Y1R * SCALE + SYCENT
          Z1S = Z1R * SCALE
          R1S = R1B * SCALE
          X2S = X2R * SCALE + SXCENT
          Y2S = Y2R * SCALE + SYCENT
          Z2S = Z2R * SCALE
          R2S = R2B * SCALE
*         save results etc.
          SDTAIL(MDET+1) = X1S
          SDTAIL(MDET+2) = Y1S
          SDTAIL(MDET+3) = Z1S
          SDTAIL(MDET+4) = R1S
          SDTAIL(MDET+5) = X2S
          SDTAIL(MDET+6) = Y2S
          SDTAIL(MDET+7) = Z2S
        SDTAIL(MDET+8) = R2S
          CALL ASSERT (SDET(INTYPE).EQ.8,'sdet(1).ne.8')
          MDET = MDET + SDET(INTYPE)
*         tally for shadow tiles the object might impinge on
          IXLO = MIN(X1S-R1S,X2S-R2S)  / NPX + 1
          IXHI = MAX(X1S+R1S,X2S+R2S)  / NPX + 1
          IYLO = MIN(Y1S-R1S,Y2S-R2S)  / NPY + 1
          IYHI = MAX(Y1S+R1S,Y2S+R2S)  / NPY + 1
          IF (IXLO.LT.1  ) IXLO=1
          IF (IXLO.GT.NSX) GO TO 716
          IF (IXHI.LT.1  ) GO TO 716
          IF (IXHI.GT.NSX) IXHI=NSX
          IF (IYLO.LT.1  ) IYLO=1
          IF (IYLO.GT.NSY) GO TO 716
          IF (IYHI.LT.1  ) GO TO 716
          IF (IYHI.GT.NSY) IYHI=NSY
          DO 715 IY=IYLO,IYHI
          DO 715 IX=IXLO,IXHI
            MOUNT(IX,IY) = MOUNT(IX,IY) + 1
715       CONTINUE
716       CONTINUE
        ENDIF
c
c     20-Aug-98 duplicate any transparent flat-ended cylinders so that
c     the inside surface can be rendered also. BUF() is still loaded
c     with specs for the current object; we just need to set flags.
      IF (  iand(FLAG(N),TRANSP).NE.0 .AND. iand(FLAG(N),FLAT).NE.0
     &  .AND. iand(FLAG(N),INSIDE).EQ.0) THEN
          FLAG(N+1) = FLAT + INSIDE
          NINSIDE   = NINSIDE + 1
          GOTO 9
      ENDIF
*
      ELSEIF (INTYPE.EQ.NORMS) THEN
*     vertex normals as given (these belong to previous triangle)
      IF (JUSTCLIPPED) GOTO 46
      IPREV = N - 1
      IF ((TYPE(IPREV).EQ.VERTEXRGB).OR.(TYPE(IPREV).EQ.VERTRANSP))
     &      IPREV = IPREV - 1
      IF ((TYPE(IPREV).EQ.VERTEXRGB).OR.(TYPE(IPREV).EQ.VERTRANSP))
     &      IPREV = IPREV - 1
      CALL ASSERT (TYPE(IPREV).EQ.TRIANG,'orphan normals')
*     Isolated objects not transformed by TMAT, but still subject to inversion
      IF (ISOLATION.GT.0) THEN
          X1C = BUF(1)
          Y1C = BUF(2)
          Z1C = BUF(3)
          X2C = BUF(4)
          Y2C = BUF(5)
          Z2C = BUF(6)
          X3C = BUF(7)
          Y3C = BUF(8)
          Z3C = BUF(9)
        IF (INVERT) THEN
          Y1C = -Y1C
          Y2C = -Y2C
          Y3C = -Y3C
        ENDIF
      ELSE
          X1A = BUF(1)
          Y1A = BUF(2)
          Z1A = BUF(3)
          X2A = BUF(4)
          Y2A = BUF(5)
          Z2A = BUF(6)
          X3A = BUF(7)
          Y3A = BUF(8)
          Z3A = BUF(9)
*     Apply rotation matrix, but not translation components
        X1B = X1A*TMAT(1,1) + Y1A*TMAT(2,1) + Z1A*TMAT(3,1)
        Y1B = X1A*TMAT(1,2) + Y1A*TMAT(2,2) + Z1A*TMAT(3,2)
        Z1B = X1A*TMAT(1,3) + Y1A*TMAT(2,3) + Z1A*TMAT(3,3)
        X2B = X2A*TMAT(1,1) + Y2A*TMAT(2,1) + Z2A*TMAT(3,1)
        Y2B = X2A*TMAT(1,2) + Y2A*TMAT(2,2) + Z2A*TMAT(3,2)
        Z2B = X2A*TMAT(1,3) + Y2A*TMAT(2,3) + Z2A*TMAT(3,3)
        X3B = X3A*TMAT(1,1) + Y3A*TMAT(2,1) + Z3A*TMAT(3,1)
        Y3B = X3A*TMAT(1,2) + Y3A*TMAT(2,2) + Z3A*TMAT(3,2)
        Z3B = X3A*TMAT(1,3) + Y3A*TMAT(2,3) + Z3A*TMAT(3,3)
*     Also apply post-rotation, if any
        X1C = RAFTER(1,1)*X1B + RAFTER(1,2)*Y1B + RAFTER(1,3)*Z1B
        Y1C = RAFTER(2,1)*X1B + RAFTER(2,2)*Y1B + RAFTER(2,3)*Z1B
        Z1C = RAFTER(3,1)*X1B + RAFTER(3,2)*Y1B + RAFTER(3,3)*Z1B
        X2C = RAFTER(1,1)*X2B + RAFTER(1,2)*Y2B + RAFTER(1,3)*Z2B
        Y2C = RAFTER(2,1)*X2B + RAFTER(2,2)*Y2B + RAFTER(2,3)*Z2B
        Z2C = RAFTER(3,1)*X2B + RAFTER(3,2)*Y2B + RAFTER(3,3)*Z2B
        X3C = RAFTER(1,1)*X3B + RAFTER(1,2)*Y3B + RAFTER(1,3)*Z3B
        Y3C = RAFTER(2,1)*X3B + RAFTER(2,2)*Y3B + RAFTER(2,3)*Z3B
        Z3C = RAFTER(3,1)*X3B + RAFTER(3,2)*Y3B + RAFTER(3,3)*Z3B
      ENDIF
C
C     If all three Z components are negative, it's facing away from us.
C     Old default treatment was to assume the normals were screwed up.
C     V2.6:     default appropriate cases (e.g. opaque triangles with no 
C             associated bounding planes) to hidden by assumption and 
C             thus not needing to be rendered. Mark as HIDDEN.
C             Default treatment is overridden by CLROPT in material spec
C
      IF (Z1C.GE.0 .AND. Z2C.GE.0 .AND. Z3C.GE.0) GOTO 718
      IF (Z1C.LT.-.01 .AND. Z2C.LT.-.01 .AND. Z3C.LT.-.01) THEN
          IF (CLROPT.EQ.2) THEN
            NINSIDE = NINSIDE + 1
            FLAG(IPREV) = ior( FLAG(IPREV), INSIDE )
          ELSE IF ((iand(FLAG(IPREV),BOUNDED).EQ.0) 
     &          .AND.(CLRITY.EQ.0 .OR. CLROPT.EQ.1)) THEN
            NHIDDEN = NHIDDEN + 1
            FLAG(IPREV) = ior( FLAG(IPREV), HIDDEN )
          ELSE
            NINSIDE = NINSIDE + 1
            FLAG(IPREV) = ior( FLAG(IPREV), INSIDE )
          ENDIF
          GOTO 718
      ENDIF
C
C     Mixed + and - Z means the triangle "wrapped around" the edge.
C     For solid objects the best we can do is pretend the edge is right here.
C     For transparent objects or 2-sided surfaces we need to invert the 
C     normals also.  The value of EDGESLOP is purely empirical; setting it 
C     either too low or too high makes some edges get coloured wrongly.  
C     Setting the HIDDEN flag for this record (NB: for the NORMALS, not for
C     the triangle itself) causes the triangle to have flat shading.
      IF (Z1C+Z2C+Z3C .LT. 0) THEN
          IF (Z1C .GT. EDGESLOP) FLAG(N) = HIDDEN
          IF (Z2C .GT. EDGESLOP) FLAG(N) = HIDDEN
          IF (Z3C .GT. EDGESLOP) FLAG(N) = HIDDEN
          Z1C = MIN(Z1C,0.)
          Z2C = MIN(Z2C,0.)
          Z3C = MIN(Z3C,0.)
      ELSE
          IF (Z1C .LT. -EDGESLOP) FLAG(N) = HIDDEN
          IF (Z2C .LT. -EDGESLOP) FLAG(N) = HIDDEN
          IF (Z3C .LT. -EDGESLOP) FLAG(N) = HIDDEN
          Z1C = MAX(Z1C,0.)
          Z2C = MAX(Z2C,0.)
          Z3C = MAX(Z3C,0.)
      ENDIF
C
718   CONTINUE
      DETAIL(NDET+1) = X1C
      DETAIL(NDET+2) = Y1C
      DETAIL(NDET+3) = Z1C
      DETAIL(NDET+4) = X2C
      DETAIL(NDET+5) = Y2C
      DETAIL(NDET+6) = Z2C
      DETAIL(NDET+7) = X3C
      DETAIL(NDET+8) = Y3C
      DETAIL(NDET+9) = Z3C
      NDET = NDET + KDET(INTYPE)
      IF (SHADOW) THEN
        MDET = MDET + SDET(INTYPE)
      ENDIF
*
*     Allow specification of RGB triple for each vertex of preceding
*     triangle or cylinder. Overrides base RGB.
*     Also overrides MATERIAL RGB, which is arguably a bug.
*
      ELSEIF (INTYPE .EQ. VERTEXRGB) THEN
      IF (JUSTCLIPPED) GOTO 46
      CALL CHKRGB(BUF(1),BUF(2),BUF(3),'invalid vertex color')
      CALL CHKRGB(BUF(4),BUF(5),BUF(6),'invalid vertex color')
      CALL CHKRGB(BUF(7),BUF(8),BUF(9),'invalid vertex color')
      IPREV = N - 1
      IF ((TYPE(IPREV).EQ.NORMS).OR.(TYPE(IPREV).EQ.VERTRANSP))
     &      IPREV = IPREV - 1
      IF ((TYPE(IPREV).EQ.NORMS).OR.(TYPE(IPREV).EQ.VERTRANSP))
     &      IPREV = IPREV - 1
c     we should only see a SPHERE is if it's a collapsed cylinder
      IF (TYPE(IPREV).EQ.SPHERE) THEN
          K = LIST(IPREV)
          DETAIL(K+5) = BUF(1)
          DETAIL(K+6) = BUF(2)
          DETAIL(K+7) = BUF(3)
          GOTO 7
      ENDIF
      CALL ASSERT (TYPE(IPREV).EQ.TRIANG .OR. TYPE(IPREV).EQ.CYLIND,
     &            'orphan vertex colours')
      FLAG(IPREV) = ior( FLAG(IPREV), VCOLS )
      DETAIL(NDET+1) = BUF(1)
      DETAIL(NDET+2) = BUF(2)
      DETAIL(NDET+3) = BUF(3)
      DETAIL(NDET+4) = BUF(4)
      DETAIL(NDET+5) = BUF(5)
      DETAIL(NDET+6) = BUF(6)
      DETAIL(NDET+7) = BUF(7)
      DETAIL(NDET+8) = BUF(8)
      DETAIL(NDET+9) = BUF(9)
      NDET = NDET + KDET(INTYPE)
      IF (SHADOW) THEN
        MDET = MDET + SDET(INTYPE)
      ENDIF
*
*     EAM - 30-Dec-1999
*     Allow specification of transparency at each vertex of preceding
*     triangle or cylinder. Overrides any MATERIAL properties.
*
      ELSEIF (INTYPE .EQ. VERTRANSP) THEN
      IF (JUSTCLIPPED) GOTO 46
      IPREV = N - 1
      IF (TYPE(IPREV).EQ.NORMS .OR. TYPE(IPREV).EQ.VERTEXRGB)
     &      IPREV = IPREV - 1
      IF (TYPE(IPREV).EQ.NORMS .OR. TYPE(IPREV).EQ.VERTEXRGB)
     &      IPREV = IPREV - 1
      CALL ASSERT (TYPE(IPREV).EQ.TRIANG .OR. TYPE(IPREV).EQ.CYLIND
     &           .OR.TYPE(IPREV).EQ.SPHERE,
     &            'orphan vertex transparency')
      NVTRANS = NVTRANS + 1
      IF (iand(FLAG(IPREV),TRANSP).EQ.0) NTRANSP = NTRANSP + 1
      FLAG(IPREV) = ior( FLAG(IPREV), TRANSP )
      FLAG(IPREV) = ior( FLAG(IPREV), VTRANS )
      DETAIL(NDET+1) = BUF(1)
      DETAIL(NDET+2) = BUF(2)
      DETAIL(NDET+3) = BUF(3)
      NDET = NDET + KDET(INTYPE)
      IF (SHADOW) THEN
        MDET = MDET + SDET(INTYPE)
      ENDIF
*
*     Material properties are saved after enforcing legality
*
      ELSEIF (INTYPE .EQ. MATERIAL) THEN
*     Mark this object as current material
      MSTATE = MATERIAL
      NPROPM = NPROPM + 1
*
*       Expand arrays MLIST, MPARITY
        IF ( NPROPM .GT. SIZE(MPARITY) ) THEN
            CALL GET_TRY(SIZE(MPARITY), NPROPM, TRY1, 1)
              DO 940 ITRY = 1,3
                  if (TRY1(ITRY).LE.0.OR.TRY1(ITRY).GT.MAXMEM) GOTO 940
                  ALLOCATE( TMP1D(TRY1(ITRY)), stat=ierr )
                  if (ierr .NE. 0) GOTO 940
                  TMP1D = 0
                  TMP1D = MLIST
                  CALL MOVE_ALLOC(from=TMP1D, to=MLIST)
                  ALLOCATE( TMP1D(TRY1(ITRY)), stat=ierr )
                  if (ierr .NE. 0) GOTO 940
                  TMP1D = 0
                  TMP1D = MPARITY
                  CALL MOVE_ALLOC(from=TMP1D, to=MPARITY)
                  if(TEST_ALLOC)write(NOISE,*)"Expand MAXMAT to ",
     &                try1(ITRY)
                  GOTO 942
940           CONTINUE
        ENDIF
942     CALL ASSERT(NPROPM.LE.SIZE(MPARITY),
     &   'too many materials - increase MAXMAT and recompile')
      MLIST(NPROPM) = N
*     Clear any previous material properties
      FLAG(N)  = NPROPM*65536
      MATCOL   = .FALSE.
      CLIPPING = .FALSE.
      NBOUNDS  = 0
      ORTEPLIKE= .FALSE.
      BPRGB(1) = -1
*     Phong power defaults to value in header
      IF (BUF(1).LT.0) BUF(1) = IPHONG
      DETAIL(NDET+1) = BUF(1)
*     Specular reflection component defaults to value in header
      IF (BUF(2).LT.0 .OR. BUF(2).GT.1) BUF(2) = SPECLR
      DETAIL(NDET+2) = BUF(2) 
*     Negative values for specular highlighting indicate default to object
        CALL ASSERT (BUF(3).LE.1., 'red > 1 in material')
        CALL ASSERT (BUF(4).LE.1., 'grn > 1 in material')
        CALL ASSERT (BUF(5).LE.1., 'blu > 1 in material')
      DETAIL(NDET+3) = BUF(3)
      DETAIL(NDET+4) = BUF(4)
      DETAIL(NDET+5) = BUF(5)
      CLRITY = BUF(6)
      CALL ASSERT (CLRITY.GE.0., 'clarity < 0 in material')
      CALL ASSERT (CLRITY.LE.1., 'clarity > 1 in material')
      DETAIL(NDET+6) = CLRITY
*     Transparency processing is necessarily a compromise, and several
*     possible approximations may be useful; allow a choice among them
      CLROPT = BUF(7)
      DETAIL(NDET+7) = BUF(7)
*     The next one is used in conjunction with bounding planes
*     Dec 2010: no, it's used to select which transparency algorithm is used
      TRNSPOPT = BUF(8)
      DETAIL(NDET+8) = TRNSPOPT
*     One remaining field reserved for future expansion
      DETAIL(NDET+9) = BUF(9)
*     Initialize clipping planes, only used if CLIPPING is set below
      DETAIL(NDET+16) = FRONTCLIP
      DETAIL(NDET+17) = BACKCLIP
*     Additional properties may continue on extra lines
      IF (INT(BUF(10)).GT.0) THEN
        L = 1
        DO I = 1,INT(BUF(10))
        READ (INPUT,'(A)',END=50) LINE
        DO J = 132, 1, -1
          IF (LINE(J:J).NE.' '.AND.LINE(J:J).NE.'     ') L = J
        ENDDO
        IF (LINE(L:L+4).EQ.'SOLID') THEN
          MATCOL = .TRUE.
          READ (LINE(L+6:132),*,END=720) RGBMAT(1),RGBMAT(2),RGBMAT(3)
        ELSE IF (LINE(L:L+7).EQ.'BACKFACE') THEN
          FLAG(N) = ior(FLAG(N),INSIDE)
          READ (LINE(L+9:132),*,END=720) RED, GRN, BLU, PHONGM, SPECM
          MPHONG = PHONGM
          IF (PHONGM.LT.0) MPHONG = IPHONG
          IF (SPECM.LT.0.OR.SPECM.GT.1.) SPECM  = SPECLR
          DETAIL(NDET+11) = RED
          DETAIL(NDET+12) = GRN
          DETAIL(NDET+13) = BLU
          DETAIL(NDET+14) = MPHONG
          DETAIL(NDET+15) = SPECM
        ELSE IF (LINE(L:L+8).EQ.'FRONTCLIP') THEN
          CLIPPING = .TRUE.
          ZCLIP = FRONTCLIP
          READ (LINE(L+10:132),*,END=720) ZCLIP
          ZCLIP = ZCLIP * SCALE / TMAT(4,4)
          DETAIL(NDET+16) = ZCLIP
        ELSE IF (LINE(L:L+7).EQ.'BACKCLIP') THEN
          CLIPPING = .TRUE.
          ZCLIP = BACKCLIP
          READ (LINE(L+9:132),*,END=720) ZCLIP
          ZCLIP = ZCLIP * SCALE / TMAT(4,4)
          DETAIL(NDET+17) = ZCLIP
        ELSE IF (LINE(L:L+9).EQ.'ORTEP_LIKE') THEN
          ORTEPLIKE = .TRUE.
        ELSE IF (LINE(L:L+13).EQ.'BOUNDING_COLOR') THEN
          READ (LINE(L+15:132),*,END=720) BPRGB(1),BPRGB(2),BPRGB(3)
          CALL CHKRGB(BPRGB(1),BPRGB(2),BPRGB(3),
     &                  'Invalid bounding color')
        ELSE IF (LINE(L:L+13).EQ.'BOUNDING_PLANE') THEN
          NBOUNDS  = NBOUNDS + 1
          nbplanes = nbplanes + 1
*
*           Expand arrays DETAIL and SDTAIL as needed
            NEEDMEM = MAX(N+NBOUNDS, NDET+NBOUNDS*KDET(INTERNAL))
            if ( NEEDMEM .GT. SIZE(DETAIL) ) THEN
                CALL GET_TRY(SIZE(DETAIL), NEEDMEM, TRY1, 1)
                DO 950 ITRY = 1,3
                    IF (TRY1(ITRY) .LE. 0 .OR. TRY1(ITRY) .GT. MAXMEM) 
     &                   GOTO 950
                    ALLOCATE( TMP1DR(TRY1(ITRY)), stat=ierr )
                    IF (IERR .NE. 0) GOTO 950
                    TMP1DR = 0.
                    TMP1DR = DETAIL
                    if(TEST_ALLOC)write(NOISE,*)"Expand DETAIL to ",
     &                 size(TMP1DR)
                    CALL MOVE_ALLOC(from=TMP1DR, to=DETAIL)
                    GOTO 952
950             CONTINUE
            ENDIF
952         CALL ASSERT(NEEDMEM.LE.SIZE(DETAIL),
     &        'BP Oops')
          IF (SHADOW) THEN
                NEEDMEM = MDET+NBOUNDS*SDET(INTERNAL)
                if ( NEEDMEM .GT. SIZE(DETAIL) ) THEN
                    CALL GET_TRY(SIZE(DETAIL), NEEDMEM, TRY1, 1)
                    DO 955 ITRY = 1,3
                        IF (TRY1(ITRY).LE.0.OR.TRY1(ITRY).GT.MAXMEM) 
     &                       GOTO 955
                        ALLOCATE( TMP1DR(TRY1(ITRY)), stat=ierr )
                        IF (IERR .NE. 0) GOTO 955
                        TMP1DR = 0.
                        TMP1DR = SDTAIL
                        if(TEST_ALLOC)write(NOISE,*)"Expand SDTAIL to",
     &                     try1(ITRY)
                        CALL MOVE_ALLOC(from=TMP1DR, to=SDTAIL)
                        GOTO 957
955                 CONTINUE
                ENDIF
957             CALL ASSERT(NEEDMEM .LE.SIZE(SDTAIL),'BP Oops')
            ENDIF
          NB = N + NBOUNDS
          CALL ASSERT( NB.LE.SIZE(DETAIL), 'BP Oops')
c         OK, we've established there's room to store this bound;
c         Flag all properties belonging to the parent material
            TYPE(NB) = INTERNAL
            FLAG(NB) = ior(FLAG(N),PROPS)
            IF(CLRITY.GT.0.0)THEN
            FLAG(NB) = ior(FLAG(NB),TRANSP)
            IF (CLROPT.EQ.1) FLAG(NB) = ior(FLAG(NB),MOPT1)
            ENDIF
            NBDET = KDET(MATERIAL) + (NBOUNDS-1)*KDET(INTERNAL)
            LIST(NB) = NDET + NBDET
            IF (SHADOW) THEN
              NBSDT = SDET(MATERIAL) + (NBOUNDS-1)*SDET(INTERNAL)
                  MIST(NB) = MDET + NBSDT
            ENDIF
c         Read in details, transform, and save
            READ (LINE(L+15:132),*,END=720) BPTYPE,XB,YB,ZB,xn,yn,zn
*         Transform bounding plane along with objects
            CALL TRANSF(XB,YB,ZB)
*         Rotate but don't translate normal, including post-rotation
            xnb = xn*TMAT(1,1) + yn*TMAT(2,1) + zn*TMAT(3,1)
            ynb = xn*TMAT(1,2) + yn*TMAT(2,2) + zn*TMAT(3,2)
            znb = xn*TMAT(1,3) + yn*TMAT(2,3) + zn*TMAT(3,3)
            xn = RAFTER(1,1)*xnb + RAFTER(1,2)*ynb + RAFTER(1,3)*znb
            yn = RAFTER(2,1)*xnb + RAFTER(2,2)*ynb + RAFTER(2,3)*znb
            zn = RAFTER(3,1)*xnb + RAFTER(3,2)*ynb + RAFTER(3,3)*znb
            temp = sqrt(xn*xn + yn*yn + zn*zn)
            xn = xn/temp
            yn = yn/temp
            zn = zn/temp
            IF (ORTEPLIKE .AND. ZN.LT.0) THEN
                  XN = -XN
            YN = -YN
            ZN = -ZN
            ENDIF
*         Save data in same arrays as for regular objects
*         ISUBTYPE currently only used to flag ORTEP_LIKE ellipsoids
*         BPTYPE is loaded on input but not currently used for anything 
            ISUBTYPE = -1
            IF (ORTEPLIKE) ISUBTYPE = ORTEP
            DETAIL(NDET + NBDET + 1) = ISUBTYPE
            DETAIL(NDET + NBDET + 2) = BPTYPE
            DETAIL(NDET + NBDET + 3) = XB * SCALE + XCENT
            DETAIL(NDET + NBDET + 4) = YB * SCALE + YCENT
            DETAIL(NDET + NBDET + 5) = ZB * SCALE
            DETAIL(NDET + NBDET + 6) = XN
            DETAIL(NDET + NBDET + 7) = YN
            DETAIL(NDET + NBDET + 8) = ZN
c           Most of the time BPRGB(1) is -1 to signal no special color
            DETAIL(NDET + NBDET + 9)  = BPRGB(1)
            DETAIL(NDET + NBDET + 10) = BPRGB(2)
            DETAIL(NDET + NBDET + 11) = BPRGB(3)
            IF (SHADOW) THEN
                  XR = SROT(1,1)*XB+SROT(1,2)*YB+SROT(1,3)*ZB
            YR = SROT(2,1)*XB+SROT(2,2)*YB+SROT(2,3)*ZB
            ZR = SROT(3,1)*XB+SROT(3,2)*YB+SROT(3,3)*ZB
                  BPNORM(1) = SROT(1,1)*XN+SROT(1,2)*YN+SROT(1,3)*ZN
            BPNORM(2) = SROT(2,1)*XN+SROT(2,2)*YN+SROT(2,3)*ZN
            BPNORM(3) = SROT(3,1)*XN+SROT(3,2)*YN+SROT(3,3)*ZN
            SDTAIL(MDET + NBSDT + 1) = ISUBTYPE
            SDTAIL(MDET + NBSDT + 2) = BPTYPE
            SDTAIL(MDET + NBSDT + 3) = XR * SCALE + SXCENT
            SDTAIL(MDET + NBSDT + 4) = YR * SCALE + SYCENT
            SDTAIL(MDET + NBSDT + 5) = ZR * SCALE
            SDTAIL(MDET + NBSDT + 6) = BPNORM(1)
            SDTAIL(MDET + NBSDT + 7) = BPNORM(2)
            SDTAIL(MDET + NBSDT + 8) = BPNORM(3)
            ENDIF
        ELSE IF (LINE(L:L+6).EQ.'BUMPMAP') THEN
          WRITE(NOISE,*) '>> Sorry, no bumpmaps (dont you wish!)'
        ELSE
          GOTO 720
        ENDIF
        GOTO 721
720     WRITE(NOISE,'(A,A)') 
     &          '>> Unrecognized or incomplete MATERIAL option ',LINE
721     CONTINUE
        ENDDO
      ENDIF
*       Update array pointers for material object itself
        DETAIL(NDET+18) = NBOUNDS
        NDET = NDET + KDET(INTYPE)
        IF (SHADOW) MDET = MDET + SDET(INTYPE)
*     Update array pointers to allow for bounding planes and any other
*     objects inserted by MATERIAL processing
        NDET = NDET + NBOUNDS*IDET(INTERNAL)
        IF (SHADOW) MDET = MDET + NBOUNDS*KDET(INTERNAL)
        N = N + NBOUNDS
*
      ELSEIF (INTYPE.EQ.GLOWLIGHT) THEN
      NGLOWS = NGLOWS + 1
*
*       Expand array GLOWLIST
        IF ( NGLOWS .GT. SIZE(GLOWLIST) ) THEN
            CALL GET_TRY(SIZE(GLOWLIST), NGLOWS, TRY1, 1)
              DO 960 ITRY = 1,3
                  if (TRY1(ITRY).LE.0.OR.TRY1(ITRY).GT.MAXMEM) GOTO 960
                  ALLOCATE( TMP1D(TRY1(ITRY)), stat=ierr )
                  if (ierr .NE. 0) GOTO 960
                  TMP1D = GLOWLIST
                  CALL MOVE_ALLOC(from=TMP1D, to=GLOWLIST)
                  if(TEST_ALLOC)write(NOISE,*)"Expand GLOWLIST to ",
     &                 try1(ITRY)
                  GOTO 962
960           CONTINUE
        ENDIF
962     CALL ASSERT(NGLOWS.LE.SIZE(GLOWLIST),'too many glow lights')
      GLOWLIST(NGLOWS) = N
      GLOWSRC(1) = BUF(1)
      GLOWSRC(2) = BUF(2)
      GLOWSRC(3) = BUF(3)
      GLOWRAD    = BUF(4)
      GLOW       = BUF(5)
      GOPT       = BUF(6)
      GPHONG     = BUF(7)
      RED        = BUF(8)
      GRN        = BUF(9)
      BLU        = BUF(10)
      CALL ASSERT (GLOW.GE.0,'illegal glow value')
      CALL ASSERT (GLOW.LE.1,'illegal glow value')
      IF (GLOW.GT.GLOWMAX) GLOWMAX = GLOW
      CALL CHKRGB (RED,GRN,BLU,'invalid glow light color')
*     Isolated objects not transformed by TMAT, but still subject to inversion
      IF (ISOLATION.GT.0) THEN
        CALL ISOLATE(GLOWSRC(1),GLOWSRC(2))
      ELSE
*     transform coordinates and radius of glow source
        CALL TRANSF(GLOWSRC(1),GLOWSRC(2),GLOWSRC(3))
        GLOWRAD = GLOWRAD / TMAT(4,4)
      ENDIF
      IF (EYEPOS.GT.0) PFAC = PERSP( GLOWSRC(3) )
*     save for rendering
      DETAIL(NDET+1)  = GLOWSRC(1) * PFAC * SCALE + XCENT
      DETAIL(NDET+2)  = GLOWSRC(2) * PFAC * SCALE + YCENT
      DETAIL(NDET+3)  = GLOWSRC(3) * PFAC * SCALE
      DETAIL(NDET+4)  = GLOWRAD    * PFAC * SCALE
      DETAIL(NDET+5)  = GLOW
      DETAIL(NDET+6)  = GOPT
      DETAIL(NDET+7)  = GPHONG
      DETAIL(NDET+8)  = RED
      DETAIL(NDET+9)  = GRN
      DETAIL(NDET+10) = BLU
      NDET = NDET + KDET(INTYPE)
      IF (SHADOW) THEN
        MDET = MDET + SDET(INTYPE)
      ENDIF

      ELSEIF (INTYPE.EQ.QUADRIC) THEN
      NQUADS = NQUADS + 1
      IF (MSTATE.EQ.MATERIAL) THEN
        FLAG(N) = ior(FLAG(N),PROPS)
        NPROPS  = NPROPS + 1
        IF (CLRITY.GT.0) THEN
          FLAG(N) = ior(FLAG(N),TRANSP)
          IF (CLROPT.EQ.1) FLAG(N) = ior(FLAG(N),MOPT1)
          NTRANSP = NTRANSP + 1
          ISTRANS = 1
        ENDIF
        IF (CLIPPING) THEN
          FLAG(N) = ior(FLAG(N),CLIPPED)
        ENDIF
        IF (NBOUNDS.GT.0) FLAG(N) = ior(FLAG(N),BOUNDED)
        IF (MATCOL) THEN
          BUF(5) = RGBMAT(1)
          BUF(6) = RGBMAT(2)
          BUF(7) = RGBMAT(3)
        ENDIF
      ENDIF
*
      ISQUAD = QINP( BUF(1), DETAIL(NDET+1), SHADOW, SDTAIL(MDET+1) )
*
      IF (.NOT. ISQUAD) GOTO 45
      NDET = NDET + KDET(INTYPE)
      IF (SHADOW) THEN
        MDET = MDET + SDET(INTYPE)
      ENDIF

C
C New object types go here!
C

      ELSEIF (INTYPE.EQ.INTERNAL) THEN
        CALL ASSERT(.FALSE.,'object type 4 not available')
      ELSE
        CALL ASSERT(.FALSE.,'crash 50')
      ENDIF
      GO TO 7
c
c     here to discard object due to clipping planes
c
45    CONTINUE
      NCLIP = NCLIP + 1
46    FLAG(N) = 0
      N = N - 1
      GO TO 7

c
c     26-Aug-1999 attempt error recovery and reporting 
c             if input line is not recognized
47    continue
      write (noise,'(A,A)') 'Unrecognized line: ',LINE
      goto 7

*
*     here for end of objects
*
50    CONTINUE
      IF (INPUT.GT.INPUT0) THEN
      IF (VERBOSE) WRITE (NOISE,*) ' - closing indirect input file'
      CLOSE(INPUT)
      INPUT = INPUT - 1
      GOTO 7
      ENDIF
*    
*     help people re-center objects
*
      XA = (TRULIM(1,1) + TRULIM(1,2)) / 2.0
      YA = (TRULIM(2,1) + TRULIM(2,2)) / 2.0
      ZA = (TRULIM(3,1) + TRULIM(3,2)) / 2.0
      CALL TRANSF( XA, YA, ZA )
      XA = TMAT(4,1) - XA * TMAT(4,4)
      YA = TMAT(4,2) - YA * TMAT(4,4)
      ZA = TMAT(4,3) - ZA * TMAT(4,4)
      IF (INVERT) YA = -YA
      WRITE (NOISE,*) 'To center objects in rendered scene, ',
     &                'change translation to:'
      WRITE (NOISE,'(3F10.4)') XA, YA, ZA
*
*     Now we can set depth-cueing 
      WRITE (LINE,577)    'Z limits (unscaled) before clipping:',
     &            ZLIM(1)*TMAT(4,4)/SCALE,ZLIM(2)*TMAT(4,4)/SCALE
      WRITE (NOISE,*) LINE(2:57)
      WRITE (LINE,577)    'Z-clipping limits:                  ', 
     &            BACKCLIP*TMAT(4,4)/SCALE
      IF (FRONTCLIP.EQ.HUGE) THEN
                  WRITE (LINE(48:57),'(A10)') '      none'
      ELSE
                  WRITE (LINE(48:57),'(F10.2)') FRONTCLIP*TMAT(4,4)/SCALE
      ENDIF
      WRITE (NOISE,*) LINE(2:57)
      IF (VERBOSE) WRITE (NOISE,*) 'Scale: ', SCALE
      IF (VERBOSE.AND.GAMMACORRECTION) WRITE (NOISE,*) 'Gamma: ', Gamma
      IF (FOGTYPE .GE. 0) THEN
        IF (FOGBACK .EQ. 0) THEN
          FOGLIM(1) = ZLIM(1)
      ELSE
          FOGLIM(1) = FOGBACK * BACKCLIP
      ENDIF
      IF (FOGFRONT .EQ. 0) THEN
          FOGLIM(2) = ZLIM(2)
      ELSE IF (FRONTCLIP .LT. HUGE) THEN
          FOGLIM(2) = FOGFRONT * FRONTCLIP
      ELSE
          FOGLIM(2) = FOGFRONT * SCALE
      ENDIF
      IF (FOGTYPE.EQ.1) THEN
         WRITE(LINE,577) 'Fog (exponential) limits, density:'
      ELSE
         WRITE(LINE,577) 'Fog (linear) limits, density:     '
      ENDIF
      WRITE (LINE(38:67),578) 
     &         FOGLIM(1)*TMAT(4,4)/SCALE,FOGLIM(2)*TMAT(4,4)/SCALE,FOGDEN
      WRITE (NOISE,*) LINE(2:67)
      ENDIF
  577 FORMAT(1X,A,2F10.2)
  578 FORMAT(2F10.2,F10.2)
*
*     Check list for special objects 
*     Triangle types first (vanilla/ribbon/surface)
      NRIB = 0
      NSUR = 0
      NTRI = 0
      DO 55 I = 1, N-1
      IF (TYPE(I).EQ.TRIANG) THEN
        NTRI = NTRI + 1
*       Allow IPHONG=0 to disable special processing of triangles
        IF (IPHONG.EQ.0) GOTO 54
*       Check for surface triangle (explicit normals in next record)
CDEBUG        "I+1" should be replaced by INEXT processing
        IF (TYPE(I+1).EQ.NORMS.AND.iand(FLAG(I+1),HIDDEN).EQ.0) THEN
            FLAG(I) = ior( FLAG(I), SURFACE )
          GOTO 54
        END IF
        IF (I.EQ.1) GOTO 54
*       Check for ribbon triangles,
*       can't possibly be one unless surrounded by other triangles
        IPREV = I - 1
        INEXT = I + 1
        IF (TYPE(IPREV).NE.TRIANG .OR. TYPE(INEXT).NE.TRIANG) THEN
          FLAG(I) = iand( FLAG(I), NOT(RIBBON) )
          GOTO 54
        END IF
*         trailing vertex must match one in previous triangle
        J = LIST(IPREV)
        K = LIST(I)
        L = LIST(INEXT)
        DO II = 1, 3
          KK = K+II
          if (  abs(detail(kk)-detail(j+ii+3)).gt.RIBBONSLOP
     &          .and. abs(detail(kk)-detail(j+ii+6)).gt.RIBBONSLOP) goto 54
C         IF   (DETAIL(KK).NE.DETAIL(J+II+3)
C    &      .AND. DETAIL(KK).NE.DETAIL(J+II+6)) GOTO 54
        END DO
*         leading vertex must match one in following triangle
        DO II = 7, 9
          KK = K+II
          IF   (abs(DETAIL(KK)-DETAIL(L+II-3)).gt.RIBBONSLOP
     &      .AND. abs(DETAIL(KK)-DETAIL(L+II-6)).gt.RIBBONSLOP) GOTO 54
C         IF   (DETAIL(KK).NE.DETAIL(L+II-3)
C    &      .AND. DETAIL(KK).NE.DETAIL(L+II-6)) GOTO 54
        END DO
        FLAG(I) = ior(FLAG(I),RIBBON)
54      CONTINUE
        IF (iand(FLAG(I),RIBBON).NE.0)  NRIB = NRIB + 1
        IF (iand(FLAG(I),SURFACE).NE.0) NSUR = NSUR + 1
      END IF
55    CONTINUE
      IF (TYPE(N).EQ.TRIANG) NTRI = NTRI + 1
56    CONTINUE
 
*     Set GLOW to maximum requested by glow light sources and bump up
*     ambient contribution to compensate for darkening applied later
      AMBIEN = AMBIEN * (1. + GLOWMAX)
*
      WRITE(NOISE,*)'-------------------------------'
      IF (NSPHERE.NE.0) WRITE(NOISE,57) 'spheres           =',NSPHERE
      IF (NCYLIND.NE.0) WRITE(NOISE,57) 'cylinders         =',NCYLIND
      NTRI = NTRI - (NRIB + NSUR)
      IF (NPLANES.NE.0) WRITE(NOISE,57) 'planes            =',NPLANES
      IF (NTRI.NE.0) WRITE(NOISE,57)    'plain triangles   =',NTRI
      IF (NRIB.NE.0) WRITE(NOISE,57)    'ribbon triangles  =',NRIB
      IF (NSUR.NE.0) WRITE(NOISE,57)    'surface triangles =',NSUR
      IF (NQUADS.NE.0) WRITE(NOISE,57)  'quadric surfaces  =',NQUADS
      IF (NPROPM.NE.0) WRITE(NOISE,57)  'special materials =',NPROPM
      IF (NCLIP.NE.0)   WRITE(NOISE,57) 'Z-clipped objects =',NCLIP
      IF (NTRANSP.NE.0) WRITE(NOISE,57) 'transparent objs  =',NTRANSP
      IF (NHIDDEN.NE.0) WRITE(NOISE,57) 'hidden surfaces   =',NHIDDEN
      IF (NINSIDE.NE.0) WRITE(NOISE,57) 'inside surfaces   =',NINSIDE
      IF (nbplanes.NE.0) WRITE(NOISE,57)'bounding planes   =',nbplanes
      WRITE(NOISE,57)                   'total objects     =',N
      WRITE(NOISE,*)'-------------------------------'
      IF (NGLOWS.GT.0)  WRITE(NOISE,57) 'glow lights       =',NGLOWS
      IF (LFLAG) THEN
      CALL LCLOSE( NLABELS )
            WRITE(NOISE,57)                 'labels            =',NLABELS
        WRITE(NOISE,*)'-------------------------------'
      ELSEIF (NLABELS.NE.0) THEN
        WRITE(NOISE,57) 'labels (ignored)  =',NLABELS
        WRITE(NOISE,*)'-------------------------------'
      ENDIF
57    FORMAT(2X,A,I8)
*
      ZSLOP = SLOP * MAX(NPX,NPY)
      IF (VERBOSE) THEN
        WRITE (NOISE,*) 'ndet  =',NDET,' SIZE(DETAIL)=',SIZE(DETAIL)
        IF (SHADOW) WRITE (NOISE,*) 'mdet  =',MDET,' MAXSDT=',
     &    SIZE(SDTAIL)
      WRITE (NOISE,*) 'EDGESLOP =',EDGESLOP
      WRITE (NOISE,*) '   ZSLOP =',   ZSLOP
      ENDIF
*
*
*     Sort objects, fill in "short lists" as indices into main list
*     (note that it would lend itself better to "parallel
*     processing" to form the short lists first and then
*     sort each one - maybe even slightly more efficient in
*     the present context!)
      DO 60 I = 1, N
        K = LIST(I)
        CALL ASSERT (K.GE.0,'k<0')
        CALL ASSERT (K.LT.NDET,'k>=ndet')
        IF (TYPE(I).EQ.TRIANG) THEN
          Z1 = DETAIL(K+3)
          Z2 = DETAIL(K+6)
          Z3 = DETAIL(K+9)
          ZTEMP(I) = MAX(Z1,Z2,Z3)
        ELSEIF (TYPE(I).EQ.SPHERE) THEN
          Z = DETAIL(K+3)
          R = DETAIL(K+4)
          ZTEMP(I) = Z + R
      ELSEIF (TYPE(I).EQ.CYLIND) THEN
*       EAM May 1990
        Z1 = DETAIL(K+3)
        Z2 = DETAIL(K+7)
        R1 = DETAIL(K+4)
        R2 = DETAIL(K+8)
        ZTEMP(I) = MAX(Z1+R1,Z2+R2)
      ELSEIF (TYPE(I).EQ.PLANE
     &     .OR. TYPE(I).EQ.NORMS 
     &     .OR. TYPE(I).EQ.MATERIAL
     &     .OR. TYPE(I).EQ.VERTRANSP
     &     .OR. TYPE(I).EQ.VERTEXRGB
     &     .OR. TYPE(I).EQ.GLOWLIGHT
     &     .OR. TYPE(I).EQ.INTERNAL) THEN
*       EAM Mar 1994 (not sure this is necessary)
        ZTEMP(I) = SCALE + 1.0
        ELSEIF (TYPE(I).EQ.QUADRIC) THEN
          Z = DETAIL(K+3)
          R = DETAIL(K+4)
          ZTEMP(I) = Z + R
        ELSE
          CALL ASSERT(.FALSE.,'crash 60')
        ENDIF
60    CONTINUE
      CALL HSORTD (N, ZTEMP, ZINDEX)
      KNTTOT = 0
      KNTMAX = 0
      DO J = 1, NTY
      DO I = 1, NTX
        KNTTOT = KNTTOT + KOUNT(I,J)
      IF (KOUNT(I,J).GT.KNTMAX) KNTMAX = KOUNT(I,J)
      ENDDO
      ENDDO
      IF (VERBOSE) THEN
      write (noise,*) 'max/avg length of short lists=',
     &                   kntmax,(knttot/(ntx*nty))+1
            WRITE (NOISE,*) 'knttot=',KNTTOT,' MAXSHR=',SIZE(KSHORT)
      ENDIF
*
*     Expand array KSHORT
      IF ( KNTTOT .GT. SIZE(KSHORT) ) THEN
          CALL GET_TRY(SIZE(KSHORT), KNTTOT, TRY1, 1)
          DO 970 ITRY = 1,3
              if (TRY1(ITRY).LE.0.OR.TRY1(ITRY).GT.MAXMEM) GOTO 970
              ALLOCATE( TMP1D(TRY1(ITRY)), stat=ierr )
              if (ierr .NE. 0) GOTO 970
              TMP1D = KSHORT
              CALL MOVE_ALLOC(from=TMP1D, to=KSHORT)
              if(TEST_ALLOC)write(NOISE,*)"Expand KSHORT to ",try1(ITRY)
              GOTO 972
970       CONTINUE
      ENDIF
972   CALL ASSERT (KNTTOT.LE.SIZE(KSHORT),'short list overflow')
      K = 0
      DO J = 1, NTY
      DO I = 1, NTX
        KSTART(I,J) = K+1
        KSTOP(I,J) = K
        K = K + KOUNT(I,J)
      ENDDO
      ENDDO
C     write(NOISE,*)"Set kstart and kstop ",NTX," x ",NTY
      CALL ASSERT (K.EQ.KNTTOT,'k.ne.knttot')
      DO 90 I = 1, N
        IND = ZINDEX(N-I+1)
        CALL ASSERT (IND.GE.1,'ind<1')
        CALL ASSERT (IND.LE.N,'ind>n')
        K = LIST(IND)
        CALL ASSERT (K.GE.0,'k<0')
        CALL ASSERT (K.LT.NDET,'k>=ndet')
*       impingement tests here must be same as above
        IF (TYPE(IND).EQ.TRIANG) THEN
          X1 = DETAIL(K+1)
          Y1 = DETAIL(K+2)
          Z1 = DETAIL(K+3)
          X2 = DETAIL(K+4)
          Y2 = DETAIL(K+5)
          Z2 = DETAIL(K+6)
          X3 = DETAIL(K+7)
          Y3 = DETAIL(K+8)
          Z3 = DETAIL(K+9)
          IXLO = MIN(X1,X2,X3)/NPX + 1
          IXHI = MAX(X1,X2,X3)/NPX + 1
          IYLO = MIN(Y1,Y2,Y3)/NPY + 1
          IYHI = MAX(Y1,Y2,Y3)/NPY + 1
        ELSEIF (TYPE(IND).EQ.SPHERE) THEN
          X = DETAIL(K+1)
          Y = DETAIL(K+2)
          Z = DETAIL(K+3)
          R = DETAIL(K+4)
          IXLO = (X-R)/NPX + 1
          IXHI = (X+R)/NPX + 1
          IYLO = (Y-R)/NPY + 1
          IYHI = (Y+R)/NPY + 1
      ELSEIF (TYPE(IND).EQ.CYLIND) THEN
          X1 = DETAIL(K+1)
          Y1 = DETAIL(K+2)
          Z1 = DETAIL(K+3)
        R1 = DETAIL(K+4)
          X2 = DETAIL(K+5)
          Y2 = DETAIL(K+6)
          Z2 = DETAIL(K+7)
        R2 = DETAIL(K+8)
          IXLO = MIN(X1-R1,X2-R2) / NPX + 1
          IXHI = MAX(X1+R1,X2+R2) / NPX + 1
          IYLO = MIN(Y1-R1,Y2-R2) / NPY + 1
          IYHI = MAX(Y1+R1,Y2+R2) / NPY + 1
        ELSEIF (TYPE(IND).EQ.PLANE) THEN
        IXLO = 1
        IXHI = NTX
        IYLO = 1
        IYHI = NTY
        ELSEIF (TYPE(IND).EQ.NORMS) THEN
        GOTO 81
        ELSEIF (TYPE(IND).EQ.VERTEXRGB) THEN
        GOTO 81
        ELSEIF (TYPE(IND).EQ.VERTRANSP) THEN
        GOTO 81
        ELSEIF (TYPE(IND).EQ.MATERIAL) THEN
        GOTO 81
        ELSEIF (TYPE(IND).EQ.GLOWLIGHT) THEN
        GOTO 81
      ELSEIF (TYPE(IND).EQ.QUADRIC) THEN
          X = DETAIL(K+1)
          Y = DETAIL(K+2)
          Z = DETAIL(K+3)
          R = DETAIL(K+4)
          IXLO = (X-R)/NPX + 1
          IXHI = (X+R)/NPX + 1
          IYLO = (Y-R)/NPY + 1
          IYHI = (Y+R)/NPY + 1
        ELSEIF (TYPE(IND).EQ.INTERNAL) THEN
        GOTO 81
      ELSE
          CALL ASSERT(.FALSE.,'crash 80')
        ENDIF
        IF (IXLO.LT.1  ) IXLO=1
        IF (IXLO.GT.NTX) GO TO 81
        IF (IXHI.LT.1  ) GO TO 81
        IF (IXHI.GT.NTX) IXHI=NTX
        IF (IYLO.LT.1  ) IYLO=1
        IF (IYLO.GT.NTY) GO TO 81
        IF (IYHI.LT.1  ) GO TO 81
        IF (IYHI.GT.NTY) IYHI=NTY
        DO 80 IY=IYLO,IYHI
        DO 80 IX=IXLO,IXHI
          KSTOP(IX,IY) = KSTOP(IX,IY) + 1
          KSHORT(KSTOP(IX,IY)) = IND
80      CONTINUE
C       write(NOISE,*)"Incr kstop ",IXLO,"-",IXHI," x ",IYLO,"-",IYHI,
C    &  kstart(IXLO,IYLO),kstop(IXLO,IYLO),kount(IXLO,IYLO)
81      CONTINUE
90    CONTINUE
      DO 95 J=1,NTY
      DO 95 I=1,NTX
        K1 = KSTART(I,J)
        K2 = KSTOP(I,J)
        K3 = KOUNT(I,J)
        if(K2-k1.ne.k3-1)then
           write(NOISE,*)"*** ERROR IN KOUNT,KSTART,KSTOP"
           write(NOISE,*)"I,J=",I,J," start,stop=",K1,K2," kount=",k3
           write(NOISE,*)"NTX,Y ",NTX,NTY," kount=",size(kount,1)," x",
     &       size(kount,2)," stop=",size(kstop,1),size(kstop,2)
C          write(NOISE,*)((kstart(III,JJJ),III=1,I),JJJ=1,J)
C          write(NOISE,*)((kstop(III,JJJ),III=1,I),JJJ=1,J)
C          write(NOISE,*)((kount(III,JJJ),III=1,I),JJJ=1,J)
        endif
        CALL ASSERT (K2-K1.EQ.K3-1,'k2-k1.ne.kount(i,j)-1')
        CALL ASSERT (K1.GE.1.AND.K1.LE.KNTTOT+1,'kstart(i,j)')
        CALL ASSERT (K2.GE.0.AND.K2.LE.KNTTOT,'kstop(i,j)')
95    CONTINUE
*
*     Do the short list business for shadow space too if required
      IF (SHADOW) THEN
        DO 160 I = 1, N
          K = MIST(I)
          CALL ASSERT (K.GE.0,'k.lt.0')
          CALL ASSERT (K.LT.MDET,'k.ge.mdet')
          IF (TYPE(I).EQ.TRIANG) THEN
            Z1 = SDTAIL(K+3)
            Z2 = SDTAIL(K+6)
            Z3 = SDTAIL(K+9)
            ZTEMP(I) = MAX(Z1,Z2,Z3)
          ELSEIF (TYPE(I).EQ.SPHERE) THEN
            Z = SDTAIL(K+3)
            R = SDTAIL(K+4)
            ZTEMP(I) = Z + R
        ELSEIF (TYPE(I).EQ.CYLIND) THEN
          Z1 = SDTAIL(K+3)
          Z2 = SDTAIL(K+7)
          R1 = SDTAIL(K+4)
          R2 = SDTAIL(K+8)
          ZTEMP(I) = MAX(Z1+R1,Z2+R2)
        ELSEIF (TYPE(I).EQ.PLANE) THEN
*         no shadows for plane surface
        ELSEIF (TYPE(I).EQ.NORMS) THEN
*         and certainly not for normals
        ELSEIF (TYPE(I).EQ.VERTEXRGB) THEN
        ELSEIF (TYPE(I).EQ.VERTRANSP) THEN
        ELSEIF (TYPE(I).EQ.MATERIAL) THEN
*         or surface properties
        ELSEIF (TYPE(I).EQ.GLOWLIGHT) THEN
*         you want a shadow on a light source???
          ELSEIF (TYPE(I).EQ.QUADRIC) THEN
            Z = SDTAIL(K+3)
            R = SDTAIL(K+4)
            ZTEMP(I) = Z + R
        ELSEIF (TYPE(I).EQ.INTERNAL) THEN
          ELSE
            CALL ASSERT(.FALSE.,'crash 160')
          ENDIF
160     CONTINUE
        CALL HSORTD (N, ZTEMP, ZINDEX)
        MNTTOT = 0
        DO 170 J = 1, NSY
        DO 170 I = 1, NSX
          MNTTOT = MNTTOT + MOUNT(I,J)
170     CONTINUE
        IF (VERBOSE) WRITE (NOISE,*) 'mnttot=',MNTTOT,' MAXSSL=',
     &     SIZE(MSHORT)
*
*       Expand array MSHORT
        IF ( MNTTOT .GT. SIZE(MSHORT) ) THEN
            CALL GET_TRY(SIZE(MSHORT), MNTTOT, TRY1, 1)
            DO 975 ITRY = 1,3
              if (TRY1(ITRY).LE.0.OR.TRY1(ITRY).GT.MAXMEM) GOTO 975
              ALLOCATE( TMP1D(TRY1(ITRY)), stat=ierr )
              if (ierr .NE. 0) GOTO 975
              TMP1D = MSHORT
              CALL MOVE_ALLOC(from=TMP1D, to=MSHORT)
              if(TEST_ALLOC)write(NOISE,*)"Expand MSHORT to ",try1(ITRY)
              GOTO 977
975         CONTINUE
        ENDIF
977     CALL ASSERT (MNTTOT.LE.SIZE(MSHORT),
     &  'shadow short list overflow')
        K = 0
        DO 175 J = 1, NSY
        DO 175 I = 1, NSX
          MSTART(I,J) = K+1
          MSTOP(I,J) = K
          K = K + MOUNT(I,J)
175     CONTINUE
        CALL ASSERT (K.EQ.MNTTOT,'k.ne.mnttot')
        DO 190 I = 1, N
          IND = ZINDEX(N-I+1)
          CALL ASSERT (IND.GE.1,'ind.lt.1')
          CALL ASSERT (IND.LE.N,'ind.gt.n')
          K = MIST(IND)
          CALL ASSERT (K.GE.0,'k.lt.0')
          CALL ASSERT (K.LT.MDET,'k.ge.mdet')
*         impingement tests here must be same as above
          IF (TYPE(IND).EQ.TRIANG) THEN
            X1 = SDTAIL(K+1)
            Y1 = SDTAIL(K+2)
            Z1 = SDTAIL(K+3)
            X2 = SDTAIL(K+4)
            Y2 = SDTAIL(K+5)
            Z2 = SDTAIL(K+6)
            X3 = SDTAIL(K+7)
            Y3 = SDTAIL(K+8)
            Z3 = SDTAIL(K+9)
            IXLO = MIN(X1,X2,X3)/NPX + 1
            IXHI = MAX(X1,X2,X3)/NPX + 1
            IYLO = MIN(Y1,Y2,Y3)/NPY + 1
            IYHI = MAX(Y1,Y2,Y3)/NPY + 1
          ELSEIF (TYPE(IND).EQ.SPHERE) THEN
            X = SDTAIL(K+1)
            Y = SDTAIL(K+2)
            Z = SDTAIL(K+3)
            R = SDTAIL(K+4)
            IXLO = (X-R)/NPX + 1
            IXHI = (X+R)/NPX + 1
            IYLO = (Y-R)/NPY + 1
            IYHI = (Y+R)/NPY + 1
          ELSEIF (TYPE(IND).EQ.CYLIND) THEN
            X1 = SDTAIL(K+1)
            Y1 = SDTAIL(K+2)
            Z1 = SDTAIL(K+3)
          R1 = SDTAIL(K+4)
            X2 = SDTAIL(K+5)
            Y2 = SDTAIL(K+6)
            Z2 = SDTAIL(K+7)
          R2 = SDTAIL(K+8)
            IXLO = MIN(X1-R1,X2-R2) / NPX + 1
            IXHI = MAX(X1+R1,X2+R2) / NPX + 1
            IYLO = MIN(Y1-R1,Y2-R2) / NPY + 1
            IYHI = MAX(Y1+R1,Y2+R2) / NPY + 1
        ELSEIF (TYPE(IND).EQ.PLANE) THEN
*           no shadows for plane surface
          GOTO 181
          ELSEIF (TYPE(IND).EQ.NORMS) THEN
          GOTO 181
          ELSEIF (TYPE(IND).EQ.VERTEXRGB) THEN
          GOTO 181
          ELSEIF (TYPE(IND).EQ.VERTRANSP) THEN
          GOTO 181
          ELSEIF (TYPE(IND).EQ.MATERIAL) THEN
          GOTO 181
          ELSEIF (TYPE(IND).EQ.GLOWLIGHT) THEN
          GOTO 181
          ELSEIF (TYPE(IND).EQ.QUADRIC) THEN
            X = SDTAIL(K+1)
            Y = SDTAIL(K+2)
            Z = SDTAIL(K+3)
            R = SDTAIL(K+4)
            IXLO = (X-R)/NPX + 1
            IXHI = (X+R)/NPX + 1
            IYLO = (Y-R)/NPY + 1
            IYHI = (Y+R)/NPY + 1
          ELSEIF (TYPE(IND).EQ.INTERNAL) THEN
          GOTO 181
          ELSE
            CALL ASSERT(.FALSE.,'crash 180')
          ENDIF
          IF (IXLO.LT.1  ) IXLO=1
          IF (IXLO.GT.NSX) GO TO 181
          IF (IXHI.LT.1  ) GO TO 181
          IF (IXHI.GT.NSX) IXHI=NSX
          IF (IYLO.LT.1  ) IYLO=1
          IF (IYLO.GT.NSY) GO TO 181
          IF (IYHI.LT.1  ) GO TO 181
          IF (IYHI.GT.NSY) IYHI=NSY
          DO 180 IY=IYLO,IYHI
          DO 180 IX=IXLO,IXHI
            MSTOP(IX,IY) = MSTOP(IX,IY) + 1
            MSHORT(MSTOP(IX,IY)) = IND
180       CONTINUE
181       CONTINUE
190     CONTINUE
        DO 195 J=1,NSY
        DO 195 I=1,NSX
          K1 = MSTART(I,J)
          K2 = MSTOP(I,J)
          K3 = MOUNT(I,J)
          CALL ASSERT (K2-K1.EQ.K3-1,'k2-k1.ne.mount(i,j)-1')
          CALL ASSERT (K1.GE.1.AND.K1.LE.MNTTOT+1,'mstart(i,j)')
          CALL ASSERT (K2.GE.0.AND.K2.LE.MNTTOT,'mstop(i,j)')
195     CONTINUE
      ENDIF
*
*     Paint the tiles one by one
      DO 600 JTILE = 1, NTY
      DO 500 ITILE = 1, NTX
*       bounds of this tile in pixel space
        XLO = (ITILE-1)*NPX
        XHI = ITILE*NPX
        YLO = (JTILE-1)*NPY
        YHI = JTILE*NPY
*       initialize tile to background colour
        DO 200 J = 1, NPY
        DO 200 I = 1, NPX
        DO 199 IC = 1, 3
          TILE(IC,I,J) = BKGND(IC)
199     CONTINUE
      ACHAN(I,J) = 0.0
200     CONTINUE
*       test for no objects in tile
        IF (KOUNT(ITILE,JTILE).EQ.0) GO TO 400
      NTRANSP = TTRANS(ITILE,JTILE) + nvtrans
      IJSTART = KSTART(ITILE,JTILE)
      IJSTOP  = KSTOP(ITILE,JTILE)
*       process non-empty tile
        DO 300 J = 1, NPY
        DO 300 I = 1, NPX
*         location of the pixel in pixel space
          XP = XLO + 0.5 + (I-1)
          YP = YLO + 0.5 + (J-1)
*         starting value of "highest z so far"
        ZTOP = BACKCLIP
*         the index of the object that has it
          INDTOP = 0
*       index of highest opaque object
        ZHIGH   = ZTOP
*         and number of transparent objects above it
        INDEPTH = 0
C       Clear parity counter for all BOUNDED materials
        IF (NBPLANES.GT.0) THEN
            DO M = 1, NPROPM
            MPARITY(M) = 0
            ENDDO
        ENDIF
*         find the highest pixel, using the tile's sorted list
C         DO 240 IK = KSTART(ITILE,JTILE), KSTOP(ITILE,JTILE)
          DO 240 IK = IJSTART, IJSTOP
            IND = KSHORT(IK)
            K     = LIST(IND)
          IFLAG = FLAG(IND)
*           skip if hidden surface
          IF (NHIDDEN.GT.0 .AND. iand(IFLAG,HIDDEN).NE.0) goto 240
*         further tests depend on object type
            IF (TYPE(IND).EQ.TRIANG) THEN
              X1 = DETAIL(K+1)
              Y1 = DETAIL(K+2)
              Z1 = DETAIL(K+3)
              X2 = DETAIL(K+4)
              Y2 = DETAIL(K+5)
              Z2 = DETAIL(K+6)
              X3 = DETAIL(K+7)
              Y3 = DETAIL(K+8)
              Z3 = DETAIL(K+9)
*             cheap check for done pixel
            IF (Z1.LT.ZHIGH .AND. Z2.LT.ZHIGH .AND. Z3.LT.ZHIGH)
     &              GOTO 250
              A = DETAIL(K+10)
              B = DETAIL(K+11)
              C = DETAIL(K+12)
              D = DETAIL(K+13)
*             skip object if degenerate triangle
              IF (D.EQ.0) GO TO 240
*             skip object if z not a new high
              ZP = A*XP + B*YP + C
              IF (ZP.LE.ZHIGH) GO TO 240
*             Rigorous test to see if this point is interior to triangle
*           NOTE: when lots of triangles are present, the following 3 lines
*           account for the largest single chunk of rendering time (>10%)!
              S = (X2-X1)*(YP-Y1) - (Y2-Y1)*(XP-X1)
              T = (X3-X2)*(YP-Y2) - (Y3-Y2)*(XP-X2)
              U = (X1-X3)*(YP-Y3) - (Y1-Y3)*(XP-X3)
              IF ( (S.LT.0. .OR. T.LT.0. .OR. U.LT.0.) .AND.
     &             (S.GT.0. .OR. T.GT.0. .OR. U.GT.0.) ) GO TO 240
*           Z-clipped triangles are easy
            IF (iand(IFLAG,CLIPPED).NE.0) THEN
            MIND = LIST(MLIST(IFLAG/65536))
            IF ( ZP.GT.DETAIL(MIND+16)) GOTO 240
            IF ( ZP.LT.DETAIL(MIND+17)) GOTO 240
            ENDIF
*           Use Phong shading for surface and ribbon triangles.
            IF (iand(IFLAG,SURFACE+VCOLS+VTRANS).NE.0) THEN
            V = (Y3-Y1)*(X2-X1) - (X3-X1)*(Y2-Y1)
            W = (XP-X1)*(Y3-Y1) - (YP-Y1)*(X3-X1)
            ALPHA = W / V
            BETA  = S / V
            ENDIF
            IF (iand(IFLAG,VCOLS+VTRANS).NE.0) THEN
            DETAIL(14 + LIST(IND)) = ALPHA
            DETAIL(15 + LIST(IND)) = BETA
            ENDIF
            IF (iand(IFLAG,SURFACE).NE.0) THEN
C           CALL ASSERT(TYPE(IND+1).EQ.NORMS,'lost normals')
            A1 = DETAIL(1 + LIST(IND+1))
            B1 = DETAIL(2 + LIST(IND+1))
            C1 = DETAIL(3 + LIST(IND+1))
            A2 = DETAIL(4 + LIST(IND+1))
            B2 = DETAIL(5 + LIST(IND+1))
            C2 = DETAIL(6 + LIST(IND+1))
            A3 = DETAIL(7 + LIST(IND+1))
            B3 = DETAIL(8 + LIST(IND+1))
            C3 = DETAIL(9 + LIST(IND+1))
            TEMPNORM(1) = A1 + ALPHA*(A2-A1) + BETA*(A3-A1)
            TEMPNORM(2) = B1 + ALPHA*(B2-B1) + BETA*(B3-B1)
            TEMPNORM(3) = C1 + ALPHA*(C2-C1) + BETA*(C3-C1)
*           For ribbon triangles we take this normal for "middle" vertex,
*           normal of previous triangle for "trailing" vertex normal,
*           normal of next triangle for "leading" vertex normal.
*           Then we use linear interpolation of vertex normals.
            ELSE IF (iand(IFLAG,RIBBON).NE.0) THEN
            IPREV = IND - 1
            INEXT = IND + 1
C           CALL ASSERT(TYPE(IPREV).EQ.TRIANG,'lost triangle')
C           CALL ASSERT(TYPE(INEXT).EQ.TRIANG,'lost triangle')
            V = (Y3-Y1)*(X2-X1) - (X3-X1)*(Y2-Y1)
            W = (XP-X1)*(Y3-Y1) - (YP-Y1)*(X3-X1)
            ALPHA = W / V
            BETA  = S / V
            AT = DETAIL(10 + LIST(IPREV))
            BT = DETAIL(11 + LIST(IPREV))
            AL = DETAIL(10 + LIST(INEXT))
            BL = DETAIL(11 + LIST(INEXT))
            TEMPNORM(1) = -AT -ALPHA*(A-AT) -BETA*(AL-AT)
            TEMPNORM(2) = -BT -ALPHA*(B-BT) -BETA*(BL-BT)
            TEMPNORM(3) = 1.
            ELSE
                  TEMPNORM(1) = -A
                  TEMPNORM(2) = -B
                  TEMPNORM(3) = 1.
            END IF
*           Check bounding planes.
*           This is different for triangles than for other shapes, as we assume
*           that each triangle is only a facet of a larger shape that is really
*           the 'object' being bounded. This means that rather than checking top
*           and bottom surfaces of the current object, we have to search for 
*           them in other triangle/facets of the same bounded material.
            IF (iand(IFLAG,BOUNDED).NE.0) THEN
              MAT = IFLAG/65536
            M = MLIST(MAT)+1
            IF (.NOT.INBOUNDS( M, TYPE, LIST, DETAIL,
     &                     XP,YP,ZP,DX,DY,DZ,ZP, BPIND )) THEN
               GOTO 240
            ENDIF
c           If this surface was above bounding plane, track parity
            MIND = LIST(MLIST(MAT))
            IF (BPIND.NE.0) THEN
                IF (MPARITY(MAT).EQ.0) THEN
                  MPARITY(MAT) = 1
                  IF (ZP.LE.ZHIGH) GO TO 240
                        TEMPNORM(1) = DX
                        TEMPNORM(2) = DY
                        TEMPNORM(3) = DZ
CORTEP                  Very ugly code to force bounding plane colors to be used
CORTEP                  but only if they are present.
                  IF ( (BPIND.GT.0)
     &                  .AND.(DETAIL(LIST(BPIND)+9).GE.0.0))
     &                        IND = BPIND
                ELSE
                  MPARITY(MAT) = 0
                  IF (FLAG(INDTOP)/65536.EQ.MAT) THEN
                    IF (iand(IFLAG,TRANSP).EQ.0) THEN
                      INDTOP = 0
                      ZHIGH  = BACKCLIP
                      ZTOP   = BACKCLIP
                      INDEPTH= 0
                    ELSE IF (INDEPTH.LE.1) THEN
                      INDTOP = 0
                      ZHIGH  = BACKCLIP
                      ZTOP   = BACKCLIP
                      INDEPTH= 0
                    ELSE
                      INDEPTH = INDEPTH - 1
                      DO L = 1, INDEPTH
                        INDLIST(L) = INDLIST(L+1)
                        ZLIST(L)   = ZLIST(L+1)
                        NORMLIST(1,L) = NORMLIST(1,L+1)
                        NORMLIST(2,L) = NORMLIST(2,L+1)
                        NORMLIST(3,L) = NORMLIST(3,L+1)
                      ENDDO
                      ZTOP    = ZLIST(1)
                    ENDIF
                  ENDIF
                  GOTO 240
                ENDIF
            ENDIF
            ENDIF
*             update values for object having highest z here yet
*           19-Feb-2002 Must wait til here to set NORMAL
            NORMAL(1) = TEMPNORM(1)
            NORMAL(2) = TEMPNORM(2)
            NORMAL(3) = TEMPNORM(3)
            IF (NTRANSP.GT.0) THEN
            CALL RANK( IND, ZP, NORMAL, FLAG )
            ELSE
            ZHIGH  = ZP
            INDTOP = IND
            ENDIF
            ELSEIF (TYPE(IND).EQ.SPHERE) THEN
              X = DETAIL(K+1)
              Y = DETAIL(K+2)
              Z = DETAIL(K+3)
              R = DETAIL(K+4)
*             cheap check for done pixel
              IF (Z+R.LE.ZHIGH) GO TO 250
*             skip object if point exterior
              DX = XP-X
              DY = YP-Y
              DX2 = DX**2
              DY2 = DY**2
              R2 = R**2
              IF (DX2+DY2 .GE. R2) GO TO 240
*             skip object if z not a new high
              DZ = SQRT(R2-(DX2+DY2))
C           Triggered by CLROPT=2
            IF (iand(IFLAG,TRANSP).NE.0 .AND.
     &                iand(IFLAG,INSIDE).NE.0)       DZ = -DZ
              ZP = Z+DZ
              IF (ZP.LE.ZHIGH) GO TO 240
*           Check bounding planes.
            IF (iand(IFLAG,BOUNDED).NE.0) THEN
              ZBACK = Z-DZ
            M = MLIST(IFLAG/65536)+1
            IF (.NOT.INBOUNDS( M, TYPE, LIST, DETAIL,
     &                     XP,YP,ZP,DX,DY,DZ,ZBACK, BPIND))
     &               GOTO 240
            ENDIF
*           Z-clipped spheres aren't too bad
            IF (iand(IFLAG,CLIPPED).NE.0) THEN
            MIND = LIST(MLIST(IFLAG/65536))
            IF (ZP.GT.DETAIL(MIND+16)) THEN
              ZP = Z-DZ
              IF (ZP.LE.ZHIGH) GOTO 240
              IF (ZP.GT.DETAIL(MIND+16)) GOTO 240
              DZ = -DZ
            ENDIF
            IF (ZP.LT.DETAIL(MIND+17)) GOTO 240
            ENDIF
*             update values for object having highest z here yet
              NORMAL(1) = DX
              NORMAL(2) = DY
              NORMAL(3) = DZ
            IF (NTRANSP.GT.0) THEN
            CALL RANK( IND, ZP, NORMAL, FLAG )
            ELSE
            ZHIGH  = ZP
            INDTOP = IND
            ENDIF
          ELSEIF (TYPE(IND).EQ.CYLIND) THEN
*             EAM May 1990
              X1 = DETAIL(K+1)
              Y1 = DETAIL(K+2)
              Z1 = DETAIL(K+3)
              R1 = DETAIL(K+4)
              X2 = DETAIL(K+5)
              Y2 = DETAIL(K+6)
              Z2 = DETAIL(K+7)
            R2 = R1
*           EAM Mar 1993 with a better understanding of how this works
*           add truly cheap test for cylinder entirely below current ZTOP
            TEMP1 = MAX(Z1+R1,Z2+R2)
            IF (TEMP1 .LE. ZHIGH) GOTO 250
*             2nd (relatively cheap) test
*           is to check distance to cylinder axis in projection
            IF (X1.EQ.X2 .AND. Y1.EQ.Y2) THEN
                  TEMP1 = 0.0
            ELSE
                  TEMP1 = ((XP-X1)*(Y2-Y1) - (YP-Y1)*(X2-X1))**2
     &                        / ((X2-X1)*(X2-X1) + (Y2-Y1)*(Y2-Y1))
            ENDIF
            IF (TEMP1 .GT. R1*R1) GOTO 240
*           Now find Z coord in pixel space of point on surface of
*           cylinder with these X and Y coords (ZP)
*           Also get coords of closest point on cylinder axis (XYZA).
            ISCYL = CYLTEST( IFLAG, AXFRAC,
     &                     X1,Y1,Z1, X2,Y2,Z2, XP,YP,ZP, R1, XA,YA,ZA )
            IF (.NOT.ISCYL) GO TO 240
*             skip object if z not a new high
              IF (ZP.LE.ZHIGH) GO TO 240
              DX = XP - XA
              DY = YP - YA
              DZ = ZP - ZA
*           Check bounding planes. Unfortunately we have to get the
*           back surface first which means dummying up a call to CYLTEST
            IF (iand(IFLAG,BOUNDED).NE.0) THEN
              ZBACK = ZP
            ISCYL = CYLTEST( ior(IFLAG,INSIDE+TRANSP), AXFRAC,
     &                    X1,Y1,Z1, X2,Y2,Z2, XP,YP,ZBACK, R1, XA,YA,ZA)
            M = MLIST(IFLAG/65536)+1
            IF (.NOT.INBOUNDS( M, TYPE, LIST, DETAIL,
     &                     XP,YP,ZP,DX,DY,DZ,ZBACK, BPIND ))
     &               GOTO 240
            ENDIF
*           Z-clipped cylinders are messy
            IF (iand(IFLAG,CLIPPED).NE.0) THEN
            MIND = LIST(MLIST(IFLAG/65536))
            IF (ZP.GT.DETAIL(MIND+16)) THEN
              ISCYL = CYLTEST( ior(IFLAG,INSIDE+TRANSP), AXFRAC,
     &                     X1,Y1,Z1, X2,Y2,Z2, XP,YP,ZP, R1, XA,YA,ZA )
              IF (ZP.LE.ZHIGH) GOTO 240
              IF (ZP.GT.DETAIL(MIND+16)) GOTO 240
            ENDIF
            IF (ZP.LT.DETAIL(MIND+17)) GOTO 240
                DX = XP - XA
                DY = YP - YA
                DZ = ZP - ZA
            ENDIF
              NORMAL(1) = DX
              NORMAL(2) = DY
              NORMAL(3) = DZ
*           if explicit vertex colors, need to keep fractional position
            IF (iand(IFLAG,ior(VCOLS,VTRANS)).NE.0) DETAIL(K+8) = AXFRAC
*             update values for object having highest z here yet
            IF (NTRANSP.GT.0) THEN
            CALL RANK( IND, ZP, NORMAL, FLAG )
            ELSE
            ZHIGH  = ZP
            INDTOP = IND
            ENDIF
          ELSEIF (TYPE(IND).EQ.PLANE) THEN
            A = DETAIL(K+1)
            B = DETAIL(K+2)
            C = DETAIL(K+3)
            D = DETAIL(K+4)
            IF (D.EQ.0) GOTO 240
            ZP = A*XP + B*YP + C
            IF (ZP.LE.ZHIGH) GOTO 240
            NORMAL(1) = -A
            NORMAL(2) = -B
            NORMAL(3) = 1.
            IF (NTRANSP.GT.0) THEN
            CALL RANK( IND, ZP, NORMAL, FLAG )
            ELSE
            ZHIGH  = ZP
            INDTOP = IND
            ENDIF
          ELSEIF (TYPE(IND).EQ.QUADRIC) THEN
*           First do cheap checks against projection of limiting sphere
            X = DETAIL(K+1)
            Y = DETAIL(K+2)
            Z = DETAIL(K+3)
            R = DETAIL(K+4)
            IF (Z+R.LE.ZHIGH) GO TO 250
            DX2 = (XP-X)**2
            DY2 = (YP-Y)**2
            R2 = R**2
            IF (DX2 + DY2 .GE. R2) GO TO 240
*           Now find Z coord (ZP) in pixel space of point on quadric surface
*           with these X and Y coords
            ISQUAD = QTEST( DETAIL(K+1), DETAIL(K+8), 
     &                          XP, YP, ZP, QNORM, .FALSE., .FALSE. )
            IF (.NOT.ISQUAD) GO TO 240
            IF (ZP.LE.ZHIGH) GO TO 240
*           Check bounding planes.
              IF (iand(IFLAG,BOUNDED).NE.0) THEN
              M = MLIST(IFLAG/65536)
              IF (DETAIL(LIST(M+1)+1) .EQ. ORTEP) THEN
                CALL ORTEPBOUNDS( M+1, TYPE, LIST, DETAIL, XP,YP,ZP,
     &                   QNORM(1),QNORM(2),QNORM(3),ZBACK, BPIND)
CORTEP              Very ugly code to force bounding plane colors to be used
CORTEP              but only if they are present.
CORTEP              An alternative would be to always set IND = BPIND, but
CORTEP              check for presence of coloring info later, in which case
CORTEP              IND itself needs to be temporarily saved somewhere.
CORTEP              Or maybe just cache BPIND now and use it later if non-zero?
                IF ( (BPIND.GT.0)
     &                   .AND.(DETAIL(LIST(BPIND)+9).GE.0.0))
     &                    IND = BPIND
              ELSE
                ISQUAD = QTEST( DETAIL(K+1), DETAIL(K+8), 
     &                          XP, YP, ZBACK, QNORM, .FALSE., .TRUE. )
                IF (.NOT.INBOUNDS( M+1, TYPE, LIST, DETAIL, XP,YP,ZP,
     &                   QNORM(1),QNORM(2),QNORM(3),ZBACK, BPIND))
     &                   GOTO 240
              ENDIF
              ENDIF
C           Z-clipping of quadric surfaces is encountered more frequently
C           than for other object types, as the limiting sphere can also 
C           cause clipping.
C           Check against limiting sphere in 3D
            MAYCLIP = .FALSE.
            DZ2 = (ZP-Z)**2
            IF (DX2+DY2+DZ2 .GT. R2) MAYCLIP = .TRUE.
            IF (iand(IFLAG,CLIPPED).NE.0) THEN
              MIND = LIST(MLIST(IFLAG/65536))
              IF (ZP.GT.DETAIL(MIND+16)) MAYCLIP = .TRUE.
            ENDIF
            IF (MAYCLIP) THEN
                ISQUAD = QTEST( DETAIL(K+1), DETAIL(K+8), 
     &                          XP, YP, ZP, QNORM, .FALSE., .TRUE. )
                IF (.NOT.ISQUAD) GO TO 240
                IF (ZP.LE.ZHIGH) GO TO 240
                DZ2 = (ZP-Z)**2
                IF (DX2+DY2+DZ2 .GT. R2) GO TO 240
                IF (iand(IFLAG,CLIPPED).NE.0) THEN
                  IF (ZP.GT.DETAIL(MIND+16)) GO TO 240
                  IF (ZP.LT.DETAIL(MIND+17)) GO TO 240
                ENDIF
            ENDIF
            NORMAL(1) = QNORM(1)
            NORMAL(2) = QNORM(2)
            NORMAL(3) = QNORM(3)
*             update values for object having highest z here yet
            IF (NTRANSP.GT.0) THEN
              CALL RANK( IND, ZP, NORMAL, FLAG )
            ELSE
              ZHIGH  = ZP
              INDTOP = IND
            ENDIF
            ELSE
              CALL ASSERT(.FALSE.,'crash 240')
            ENDIF
240       CONTINUE
250       CONTINUE
C         Apply background fog here (added 2010)
        IF (FOGTYPE .GE. 0) THEN
            FOGDIM = FOGGY( FOGLIM(2) - ZTOP )
            TILE(1,I,J) = (1.-FOGDIM)*TILE(1,I,J) + FOGDIM*FOGRGB(1)
            TILE(2,I,J) = (1.-FOGDIM)*TILE(2,I,J) + FOGDIM*FOGRGB(2)
            TILE(3,I,J) = (1.-FOGDIM)*TILE(3,I,J) + FOGDIM*FOGRGB(3)
        ENDIF

*         Background colour if we never found an object in this line of sight
          IF (INDTOP.EQ.0) GO TO 299
C       We now know this is not a background pixel so set alpha channel to 1
C       Modify later if it turns out the object is transparent
        ACHAN(I,J) = 1.0
C         Transparency processing revamped Mar 2001
C       If the top object is transparent we will have to come back here
C         later and do this all again for each object in INDLIST
          IF (NTRANSP.NE.0) THEN
C           CALL ASSERT(INDEPTH.GT.0,'INDEPTH = 0')
            ITPASS = 1
            ZHIGH  = ZLIST(INDEPTH)
            INDTOP = INDLIST(INDEPTH)
            NORMAL(1) = NORMLIST(1,INDEPTH)
            NORMAL(2) = NORMLIST(2,INDEPTH)
            NORMAL(3) = NORMLIST(3,INDEPTH)
            RGBLND(1) = BKGND(1)
            RGBLND(2) = BKGND(2)
            RGBLND(3) = BKGND(3)
        ENDIF
*         ZP is the "height" of the chosen pixel,
*         and indtop tells us which object it came from:
        ZTOP = ZHIGH
        ZP = ZTOP
255       CONTINUE
C
C       Shadowing code - look for objects that shadow the one we just found
C
          IF (SHADOW) THEN
*           locate pixel in shadow space
*           take out object translation & scaling
            XT = (XP - XCENT) / SCALE
            YT = (YP - YCENT) / SCALE
            ZT = ZP / SCALE
*           rotate light source position to z axis
            XR = SROT(1,1)*XT+SROT(1,2)*YT+SROT(1,3)*ZT
            YR = SROT(2,1)*XT+SROT(2,2)*YT+SROT(2,3)*ZT
            ZR = SROT(3,1)*XT+SROT(3,2)*YT+SROT(3,3)*ZT
*           scale and translate for shadow space
            XS = XR * SCALE + SXCENT
            YS = YR * SCALE + SYCENT
            ZS = ZR * SCALE
*           determine appropriate shadow tile
            ISTILE = XS/NPX + 1
            JSTILE = YS/NPY + 1
*         Just to get proper error message
            IF (JSTILE.LE.0) JSTILE = NSY + 1 - JSTILE
            IF (ISTILE.LE.0) ISTILE = NSX + 1 - ISTILE
          IF (JSTILE.GE.NSY) THEN
            NSYMAX = MAX(JSTILE,NSYMAX)
              INDSTP = 0
            GOTO 270
          END IF
          IF (ISTILE.GE.NSX) THEN
            NSXMAX = MAX(ISTILE,NSXMAX)
              INDSTP = 0
            GOTO 270
          END IF
*           starting value of "highest shadow space z so far"
*           and the index of the object that has it
            ZSTOP = 2.0*BACKCLIP
            INDSTP = 0
*
            DO 260 IK = MSTART(ISTILE,JSTILE), MSTOP(ISTILE,JSTILE)
              IND   = MSHORT(IK)
            IFLAG = FLAG(IND)
*             Ignore transparent objects except for the one being shaded
            IF (iand(IFLAG,TRANSP).NE.0 .AND. IND.NE.INDTOP) GOTO 260
              K = MIST(IND)
              IF (TYPE(IND).EQ.TRIANG) THEN
                X1 = SDTAIL(K+1)
                Y1 = SDTAIL(K+2)
                Z1 = SDTAIL(K+3)
                X2 = SDTAIL(K+4)
                Y2 = SDTAIL(K+5)
                Z2 = SDTAIL(K+6)
                X3 = SDTAIL(K+7)
                Y3 = SDTAIL(K+8)
                Z3 = SDTAIL(K+9)
                A = SDTAIL(K+10)
                B = SDTAIL(K+11)
                C = SDTAIL(K+12)
                D = SDTAIL(K+13)
*               cheap check for done "pixel"
            IF (Z1.LT.ZSTOP .AND. Z2.LT.ZSTOP .AND. Z3.LT.ZSTOP)
     &                GOTO 270
*               skip object if degenerate triangle
                IF (D.EQ.0) GO TO 260
*               skip object if z not a new high
                ZTEST = A*XS + B*YS + C
                IF (ZTEST.LE.ZSTOP) GO TO 260
*               skip object if point exterior
                S = (X2-X1)*(YS-Y1)-(Y2-Y1)*(XS-X1)
                T = (X3-X2)*(YS-Y2)-(Y3-Y2)*(XS-X2)
                U = (X1-X3)*(YS-Y3)-(Y1-Y3)*(XS-X3)
                IF ( (S.LT.0. .OR. T.LT.0. .OR. U.LT.0.) .AND.
     &               (S.GT.0. .OR. T.GT.0. .OR. U.GT.0.) ) GO TO 260
*           Check bounding planes
            IF (iand(IFLAG,BOUNDED).NE.0) THEN
              MAT = IFLAG/65536
              M = MLIST(MAT)+1
              IF (.NOT.INBOUNDS( M, TYPE, MIST, SDTAIL,
     &                       XS,YS,ZTEST,DX,DY,DZ,ZTEST, BPIND )) THEN
                  GOTO 260
              ENDIF
              MIND = MIST(MLIST(MAT))
              IF (BPIND.NE.0) THEN
                  IF (MPARITY(MAT).GE.0) THEN
                  MPARITY(MAT) = -1
                  ELSE
                  MPARITY(MAT) = 0
                  IF (FLAG(INDSTP)/65536.EQ.MAT) THEN
                    INDSTP = 0
                    ZSTOP  = 2.*BACKCLIP
                  ENDIF
                  GOTO 260
                  ENDIF
              ENDIF
            ENDIF
*               update values for object having highest z here yet
                ZSTOP = ZTEST
                INDSTP = IND
              ELSEIF (TYPE(IND).EQ.SPHERE) THEN
                X = SDTAIL(K+1)
                Y = SDTAIL(K+2)
                Z = SDTAIL(K+3)
                R = SDTAIL(K+4)
*               cheap check for done "pixel"
                IF (Z+R.LT.ZSTOP) GO TO 270
*               skip object if point exterior
                DX = XS-X
                DY = YS-Y
                DX2 = DX**2
                DY2 = DY**2
                R2 = R**2
                IF (DX2+DY2 .GE. R2) GO TO 260
*               skip object if z not a new high
                DZ = SQRT(R2-(DX2+DY2))
                ZTEST = Z+DZ
                IF (ZTEST.LE.ZSTOP) GO TO 260
*             Check bounding planes.
              IF (iand(IFLAG,BOUNDED).NE.0) THEN
                ZBACK = Z-DZ
              M = MLIST(IFLAG/65536)+1
              IF (.NOT.INBOUNDS( M, TYPE, MIST, SDTAIL,
     &                       XS,YS,ZTEST,DX,DY,DZ,ZBACK, BPIND))
     &                 GOTO 260
              ENDIF
*               update values for object having highest z here yet
                ZSTOP = ZTEST
                INDSTP = IND
            ELSEIF (TYPE(IND).EQ.CYLIND) THEN
*             EAM May 1990
                X1 = SDTAIL(K+1)
                Y1 = SDTAIL(K+2)
                Z1 = SDTAIL(K+3)
                R1 = SDTAIL(K+4)
                X2 = SDTAIL(K+5)
                Y2 = SDTAIL(K+6)
                Z2 = SDTAIL(K+7)
            R2 = R1
*           EAM Feb 93 - Test first to see if entire cylinder is below 
*           current top object in shadow space
            IF (MAX( Z1+R1, Z2+R2 ) .LT. ZSTOP) GOTO 270
*             Now find Z coord (ZTEST) in pixel space of point on
*             surface of cylinder with these X and Y coords
            ISCYL = CYLTEST( IFLAG, AXFRAC,
     &                  X1,Y1,Z1, X2,Y2,Z2, XS,YS,ZTEST, R1, XA,YA,ZA )
            IF (.NOT.ISCYL) GO TO 260
*               skip object if z not a new high
                IF (ZTEST.LE.ZSTOP) GO TO 260
*             Check bounding planes.
              IF (iand(IFLAG,BOUNDED).NE.0) THEN
              ISCYL = CYLTEST( ior(IFLAG,INSIDE+TRANSP), AXFRAC,
     &                      X1,Y1,Z1, X2,Y2,Z2, XS,YS,ZBACK, R1, XA,YA,ZA)
              M = MLIST(IFLAG/65536)+1
              IF (.NOT.INBOUNDS( M, TYPE, MIST, SDTAIL,
     &                       XS,YS,ZTEST,DX,DY,DZ,ZBACK, BPIND ))
     &               GOTO 260
              ENDIF
*               update values for object having highest z here yet
                ZSTOP = ZTEST
                INDSTP = IND
            ELSEIF (TYPE(IND).EQ.QUADRIC) THEN
                X = SDTAIL(K+1)
                Y = SDTAIL(K+2)
                Z = SDTAIL(K+3)
                R = SDTAIL(K+4)
*               cheap check against limiting sphere
                IF (Z+R.LT.ZSTOP) GO TO 270
                DX = XS-X
                DY = YS-Y
                R2 = R**2
              IF (DX**2 + DY**2 .GE. R2) GO TO 260
*             Now find Z coord (ZTEST) in shadow pixel space of point on
*             surface with these X and Y coords
            ISQUAD = QTEST( SDTAIL(K+1), SDTAIL(K+5), 
     &                          XS, YS, ZTEST, QNORM, .TRUE., .FALSE. )
CDEBUG                          XS, YS, ZTEST, QNORM, .TRUE., .TRUE. )
CDEBUG            16-Dec-1998 I inverted the BACKSIDE = TRUE/FALSE flags from
CDEBUG            what they "ought" to be to remove buggy shadows from a test
CDEBUG            case parabolic hyperboloid.  I don't understand why this would be
CDEBUG            necessary, and worry a bit that it breaks something else.
CDEBUG
            IF (.NOT.ISQUAD) GO TO 260
*               skip object if z not a new high
                IF (ZTEST.LE.ZSTOP) GO TO 260
*             Check bounding planes.
              IF (iand(IFLAG,BOUNDED).NE.0) THEN
              M = MLIST(IFLAG/65536)
              IF (DETAIL(LIST(M)+1) .EQ. ORTEP) THEN
                ISQUAD = QTEST( SDTAIL(K+1), SDTAIL(K+5),
     &                     XS, YS, ZBACK, QNORM, .TRUE., .TRUE. )
                IF (.NOT.INBOUNDS(M+1, 
     &                     TYPE, MIST, SDTAIL, XS,YS,ZTEST,
     &                       QNORM(1),QNORM(2),QNORM(3),ZBACK, BPIND))
     &                   GOTO 260
              ENDIF
              ENDIF
*           Test against bounding sphere in 3D
*           and if surface nearest to light is clipped, check back also
            DZ = ZTEST-Z
            IF (DX**2 + DY**2 + DZ**2 .GE. R2) THEN
                ISQUAD = QTEST( SDTAIL(K+1), SDTAIL(K+5),
     &                          XS, YS, ZTEST, QNORM, .TRUE., .TRUE. )
CDEBUG                          XS, YS, ZTEST, QNORM, .TRUE., .FALSE. )
                IF (.NOT.ISQUAD) GO TO 260
                IF (ZTEST.LE.ZSTOP) GO TO 260
                DZ = ZTEST - Z
                IF (DX**2 + DY**2 + DZ**2 .GE. R2) GO TO 260
            ENDIF
*               update values for object having highest z here yet
                ZSTOP = ZTEST
                INDSTP = IND
C           No more legal object types; should never happen
              ELSE
                CALL ASSERT(.FALSE.,'shadow tile error, crash 260')
              ENDIF
260         CONTINUE
270         CONTINUE
C       End of search for objects that shadow this one
          if ((zstop+zslop).lt.zs .and. indstp.ne.0) nslow = nslow + 1
          ELSE
            ZS = 0.
            ZSTOP = 0.
            INDSTP = INDTOP
          ENDIF
*         if roundoff made us miss the object, we are probably
*         at a pixel that is very near the edge of the object
*         from the point of view of the light source, so just
*         treat it as if not in shadow
          IF (INDSTP.EQ.0) THEN
            ZS = 0.
            ZSTOP = 0.
            INDSTP = INDTOP
          ENDIF
*
*         Pick up colours of object to be shaded
*
          K = LIST(INDTOP)
          IF (TYPE(INDTOP).EQ.TRIANG) THEN
          IF (iand(FLAG(INDTOP),VCOLS).NE.0) THEN
            ALPHA = DETAIL(K+14)
            BETA  = DETAIL(K+15)
            INEXT = INDTOP + 1
            IF ((TYPE(INEXT).EQ.NORMS).OR.(TYPE(INEXT).EQ.VERTRANSP))
     &            INEXT = INEXT + 1
            IF ((TYPE(INEXT).EQ.NORMS).OR.(TYPE(INEXT).EQ.VERTRANSP))
     &            INEXT = INEXT + 1
            K = LIST(INEXT)
            CALL ASSERT(TYPE(INEXT).EQ.VERTEXRGB,'lost vertex colors')
            RGBCUR(1) = DETAIL(K+1) 
     &                  + ALPHA*(DETAIL(K+4)-DETAIL(K+1))     
     &                  +  BETA*(DETAIL(K+7)-DETAIL(K+1))     
            RGBCUR(2) = DETAIL(K+2) 
     &                  + ALPHA*(DETAIL(K+5)-DETAIL(K+2))     
     &                  +  BETA*(DETAIL(K+8)-DETAIL(K+2))     
            RGBCUR(3) = DETAIL(K+3) 
     &                  + ALPHA*(DETAIL(K+6)-DETAIL(K+3))     
     &                  +  BETA*(DETAIL(K+9)-DETAIL(K+3))     
          ELSE
              RGBCUR(1) = DETAIL(K+14)
              RGBCUR(2) = DETAIL(K+15)
              RGBCUR(3) = DETAIL(K+16)
          ENDIF
          ELSEIF (TYPE(INDTOP).EQ.SPHERE) THEN
            RGBCUR(1) = DETAIL(K+5)
            RGBCUR(2) = DETAIL(K+6)
            RGBCUR(3) = DETAIL(K+7)
        ELSEIF (TYPE(INDTOP).EQ.CYLIND) THEN
          IF (iand(FLAG(INDTOP),VCOLS).NE.0) THEN
            FRAC = DETAIL(K+8)
            INEXT = INDTOP + 1
            IF (TYPE(INEXT).EQ.VERTRANSP) INEXT = INEXT + 1
            K = LIST(INEXT)
            CALL ASSERT(TYPE(INEXT).EQ.VERTEXRGB,'lost vertex colors')
              RGBCUR(1) = FRAC*DETAIL(K+4) + (1.-FRAC)*DETAIL(K+1)
              RGBCUR(2) = FRAC*DETAIL(K+5) + (1.-FRAC)*DETAIL(K+2)
              RGBCUR(3) = FRAC*DETAIL(K+6) + (1.-FRAC)*DETAIL(K+3)
          ELSE
              RGBCUR(1) = DETAIL(K+9)
              RGBCUR(2) = DETAIL(K+10)
              RGBCUR(3) = DETAIL(K+11)
          ENDIF
*       EAM Mar 1993 PLANE is shaded from full colour in foreground
*                    to half-saturation at horizon
        ELSEIF (TYPE(INDTOP).EQ.PLANE) THEN
          FADE = (ZP + 3.*SCALE) / (4.*SCALE)
          RGBCUR(1) = FADE * DETAIL(K+5) + (1.-FADE) * BKGND(1)
          RGBCUR(2) = FADE * DETAIL(K+6) + (1.-FADE) * BKGND(2)
          RGBCUR(3) = FADE * DETAIL(K+7) + (1.-FADE) * BKGND(3)
        ELSEIF (TYPE(INDTOP).EQ.QUADRIC) THEN
            RGBCUR(1) = DETAIL(K+5)
            RGBCUR(2) = DETAIL(K+6)
            RGBCUR(3) = DETAIL(K+7)
        ELSEIF (TYPE(INDTOP).EQ.INTERNAL) THEN
            RGBCUR(1) = DETAIL(K+9)
            RGBCUR(2) = DETAIL(K+10)
            RGBCUR(3) = DETAIL(K+11)
          ELSE
          WRITE(LINE,*) 'Top object claims to be type',TYPE(INDTOP)
            CALL ASSERT(.FALSE.,LINE)
          ENDIF
*
*       Get shading components.
*

*
*         11-May-1997 As of now, treat negative normal(3) as indicating the
*       back side of a material.  Default is to shrug and invert the normal.
*       Some material have explicit BACKFACE proterties, however.
*
        BACKFACE = .FALSE.
          IF (NORMAL(3).LE.0) THEN
            NORMAL(1) = -NORMAL(1)
            NORMAL(2) = -NORMAL(2)
            NORMAL(3) = -NORMAL(3)
          BACKFACE = .TRUE.
          IF (iand(FLAG(INDTOP),PROPS).NE.0) THEN
            K = FLAG(INDTOP) / 65536
                if(K.le.0)WRITE(NOISE,*)"FLAG(",INDTOP,")=",FLAG(INDTOP)
            CALL ASSERT(K.GT.0,'lost material definition')
            IF (iand(FLAG(MLIST(K)),INSIDE).NE.0) THEN
              K = LIST(MLIST(K))
              RGBCUR(1) = DETAIL(K+11)
              RGBCUR(2) = DETAIL(K+12)
              RGBCUR(3) = DETAIL(K+13)
            END IF
          END IF
        END IF
*
          ABSN = SQRT(NORMAL(1)**2 + NORMAL(2)**2 + NORMAL(3)**2)
c       CALL ASSERT(ABS(ABSN-1.0).LT.0.02,'>> Abnormal normal')
          NL(1) = NORMAL(1) / ABSN
          NL(2) = NORMAL(2) / ABSN
          NL(3) = NORMAL(3) / ABSN
          SDIFF = NL(3) * STRAIT*DIFFUS
          SSP = 2.*NL(3)**2 - 1.
*         We do the value check like this to avoid floating-point underflows 
*       in the Phonging.  We also save calculation time this way, because
*       the 0 case will occur often for reasonably high Phong powers.
*       Note that PHOBND**IPHONG should evaluate to the cutoff value
*         between significant and insignificant specular contribution. 
*         The contributions that are actually computed here
*         can be smaller by a factor of STRAIT*SPECLR:
        IF (SSP.LT.PHOBND .OR. IPHONG.EQ.0) THEN
            SSPEC = 0.
          ELSE
            SSPEC = SSP**IPHONG * STRAIT*SPECLR
          ENDIF
          LDOTN = SOURCE(1)*NL(1)+SOURCE(2)*NL(2)+SOURCE(3)*NL(3)
          IF (LDOTN.LE.0.) THEN
            PDIFF = 0.
            PSPEC = 0.
          PSP = 0.
          ELSE
            PDIFF = LDOTN * PRIMAR*DIFFUS
            PSP = 2.*LDOTN*NL(3) - SOURCE(3)
*           Comments as for SSPEC apply, but substitute PRIMAR for STRAIT:
          IF (PSP.LT.PHOBND .OR. IPHONG.EQ.0) THEN
              PSPEC = 0.
            ELSE
              PSPEC = PSP**IPHONG * PRIMAR*SPECLR
            ENDIF
          ENDIF
*
*         experience has shown the "spots" on dark objects to be rather
*         overpowering, especially by comparison to those on brighter
*         objects.  hence the specular reflections on dark objects are
*         now artificially scaled down by a function which relates
*         directly to the "brightness" of the object.
*         this makes such objects duller, but their
*         colour seems to come through more clearly, and they don't
*         appear more specular than the brighter objects.
*         the funny coefficients are ntsc:
          BRIGHT = 0.2 + 0.8 * SQRT(0.299*RGBCUR(1) +
     &                              0.587*RGBCUR(2) +
     &                              0.114*RGBCUR(3))
          SSPEC = SSPEC * BRIGHT
          PSPEC = PSPEC * BRIGHT
*
*       The usual case is white lighting, no transparency
*
        SPECOL(1) = 1.0
        SPECOL(2) = 1.0
        SPECOL(3) = 1.0
        SBLEND    = 1.0
        CLRITY    = 0.0
*
*       Extra properties make specular highlighting calculation a
*       bit more complex. First we have to find the MATERIAL description.
*
        IF (iand(FLAG(INDTOP),PROPS).NE.0) THEN
          K = FLAG(INDTOP) / 65536
            if(K.le.0)then
                  WRITE(NOISE,*)"FLAG(",INDTOP-1,") =",FLAG(INDTOP-1)
                  WRITE(NOISE,*)"FLAG(",INDTOP,") =",FLAG(INDTOP)
                  WRITE(NOISE,*)"FLAG(",INDTOP+1,") =",FLAG(INDTOP+1)
            endif
          CALL ASSERT(K.GT.0,'lost material definition')
          IF (iand(FLAG(MLIST(K)),INSIDE).NE.0  .AND.
     &            iand(FLAG(INDTOP),  INSIDE).NE.0) THEN
            K = LIST(MLIST(K))
            MPHONG  = DETAIL(K+14)
            SPECM   = DETAIL(K+15)
          ELSE
            K = LIST(MLIST(K))
            MPHONG  = DETAIL(K+1)
            SPECM   = DETAIL(K+2)
          ENDIF
          SPECOL(1) = DETAIL(K+3)
          SPECOL(2) = DETAIL(K+4)
          SPECOL(3) = DETAIL(K+5)
          IF (SPECOL(1).LT.0) SPECOL(1) = RGBCUR(1)
          IF (SPECOL(2).LT.0) SPECOL(2) = RGBCUR(2)
          IF (SPECOL(3).LT.0) SPECOL(3) = RGBCUR(3)
          CLRITY    = DETAIL(K+6)
*         not currently used, as MOPT(1)=1 already marked in FLAG,
*         but future interpretations of MOPT(1) might need it
          CLROPT    = DETAIL(K+7)
        ENDIF
*
*       20-Feb-2000 Allow per-vertex transparency (obj type 18)
        IF (iand(FLAG(INDTOP),VTRANS).NE.0) THEN
            IF (TYPE(INDTOP).EQ.CYLIND) THEN
                K = LIST(INDTOP)
              FRAC = DETAIL(K+8)
              INEXT = INDTOP + 1
              IF (TYPE(INEXT).EQ.VERTEXRGB) INEXT = INEXT + 1
              CALL ASSERT(TYPE(INEXT).EQ.VERTRANSP,'lost vertex transp')
              K = LIST(INEXT)
              CLRITY = FRAC*DETAIL(K+1) + (1.-FRAC)*DETAIL(K+2)
            ELSE IF (TYPE(INDTOP).EQ.TRIANG) THEN
              K = LIST(INDTOP)
            ALPHA = DETAIL(K+14)
            BETA  = DETAIL(K+15)
            INEXT = INDTOP + 1
            IF ((TYPE(INEXT).EQ.VERTEXRGB).OR.(TYPE(INEXT).EQ.NORMS))
     &              INEXT = INEXT + 1
            IF ((TYPE(INEXT).EQ.VERTEXRGB).OR.(TYPE(INEXT).EQ.NORMS))
     &              INEXT = INEXT + 1
            CALL ASSERT(TYPE(INEXT).EQ.VERTRANSP,'lost vertex transp')
            K = LIST(INEXT)
            CLRITY = DETAIL(K+1)
     &                   + ALPHA * (DETAIL(K+2)-DETAIL(K+1))
     &                   + BETA  * (DETAIL(K+3)-DETAIL(K+1))
            CALL ASSERT(CLRITY.LE.1 .AND. CLRITY.GE.0.0, 'illegal transp')
            ELSE
              CALL ASSERT(.FALSE.,'illegal vertex transp')
            ENDIF
        ENDIF

C
C     This is the only computationally intensive code (as opposed to mere
C     bookkeeping) involved in rendering transparent objects. The blend
C     factor must be some function of the clarity/transparency, but I'm not
C     sure exactly what the equation ought to be.  The cosine function in
C     TRNSPOPT option 0 below was chosen after purely empirical tests of the
C     resulting image quality. If your machine bogs down incredibly due to
C     the cosine call, then you might prefer to use TRNSPOPT 1 instead.
C     Conversely, if you don't mind the extra computation then you could use
C     TRNSPOPT 2, which is closer to an ideal model.
C     Dec 2010: controlled by OPT(2) in the MATERIAL specification record
C
        IF (CLRITY.NE.0) THEN
            ZN = ABS(NL(3)) + MODULO(TRNSPOPT,1.0)
            SELECT CASE(FLOOR(TRNSPOPT))
            CASE DEFAULT ! CASE 0
              SBLEND = .25*(1.+COS(3.1416*CLRITY*ZN))**2
            CASE (1)
              SBLEND = (1. - CLRITY*ZN)**2
            CASE (2)
              SBLEND = 1.0 - CLRITY ** ( 0.7071 / ZN )
            CASE (3)
              SBLEND = 1.0 - CLRITY
            END SELECT
        ENDIF
C
C       Final calculation of specular properties of special materials
C
        IF (iand(FLAG(INDTOP),PROPS).NE.0) THEN
          DIFFM     = 1. - (SPECM + AMBIEN)
          SDIFF     = SDIFF * DIFFM / DIFFUS
          PDIFF     = PDIFF * DIFFM / DIFFUS
          SSPEC     = 0.0
          PSPEC     = 0.0
          IF (SSP .GE. PHOBND) 
     &          SSPEC = SSP**MPHONG * STRAIT*SPECM
          IF ((PSP .GE. PHOBND) .AND. (LDOTN.GT.0))
     &          PSPEC = PSP**MPHONG * PRIMAR*SPECM
*         de-emphasize highlights from inside surface of transparent objects
*         Could use BACKFACE flag instead of INSIDE to catch non-triangles
          IF (iand(FLAG(INDTOP),INSIDE).NE.0) THEN
            SSPEC = SSPEC * (1.-SPECM)
            PSPEC = PSPEC * (1.-SPECM)
          ENDIF
        END IF
*
*         We now return you to your regular processing
*
          DO 280 KC = 1, 3
            C2ND = SBLEND*RGBCUR(KC)*(AMBIEN+SDIFF) + SSPEC*SPECOL(KC)
            CSUN = SBLEND*RGBCUR(KC)*PDIFF          + PSPEC*SPECOL(KC)
            RGBSHD(KC) = C2ND
            RGBFUL(KC) = C2ND + CSUN
280       CONTINUE

C EAM March 1997 - Support additional non-shadowing light sources
C which lie within figure and have a finite range of illumination.
        IF (NGLOWS.GT.0) THEN
          DO KC = 1,3
              RGBSHD(KC) = (1.-GLOWMAX)*RGBSHD(KC)
              RGBFUL(KC) = (1.-GLOWMAX)*RGBFUL(KC)
          ENDDO
*       Recover glow light parameters
        DO NG = 1, NGLOWS
              IG = LIST(GLOWLIST(NG))
            GLOWSRC(1) = DETAIL(IG+1)
            GLOWSRC(2) = DETAIL(IG+2)
            GLOWSRC(3) = DETAIL(IG+3)
            GLOWRAD    = DETAIL(IG+4)
            GLOW       = DETAIL(IG+5)
            GOPT       = DETAIL(IG+6)
            GPHONG     = DETAIL(IG+7)
            GLOWCOL(1) = DETAIL(IG+8)
            GLOWCOL(2) = DETAIL(IG+9)
            GLOWCOL(3) = DETAIL(IG+10)
*
            GDIST(1) = GLOWSRC(1) - XP
            GDIST(2) = GLOWSRC(2) - YP
            GDIST(3) = GLOWSRC(3) - ZP
            ABSN = SQRT(GDIST(1)**2 + GDIST(2)**2 + GDIST(3)**2)
            GDIST(1) = GDIST(1) / ABSN
            GDIST(2) = GDIST(2) / ABSN
            GDIST(3) = GDIST(3) / ABSN
            LDOTN = GDIST(1)*NL(1) + GDIST(2)*NL(2) + GDIST(3)*NL(3) 
            IF (LDOTN.LE.0) THEN
              GDIFF = 0.
              GSPEC = 0.
            ELSE
C           Might want separate diffuse param for glow; (always 1.0?)
C             GDIFF = LDOTN * DIFFUS
              GDIFF = LDOTN
              GSP   = 2.*LDOTN*NL(3) - GDIST(3)
              IF (GSP.LT.PHOBND .OR. GPHONG.EQ.0) THEN
                GSPEC = 0.
              ELSE
                GSPEC = GSP**GPHONG * SPECLR
              ENDIF
            ENDIF
C         Limit glow effect by some function of ABSN, GLOWRAD 
            IF (GOPT.EQ.3) THEN
            GFADE = MAX( 0., 1. - ABSN/GLOWRAD )
            ELSE IF (GOPT.EQ.2) THEN
            GFADE = 1./(ABSN/GLOWRAD + 1.)
            ELSE IF (GOPT.EQ.1) THEN
            GFADE = 1./(ABSN/GLOWRAD + 1.)**2
            ELSE
            GFADE = MIN( 1., 1./(ABSN/GLOWRAD)**2 )
            ENDIF
            DO KC = 1,3
C           This isn't right for transparent surfaces
              CGLO = SBLEND*RGBCUR(KC)*GDIFF + GSPEC
            CGLO = GFADE * GLOWCOL(KC) * CGLO
              RGBSHD(KC) = RGBSHD(KC) + CGLO
              RGBFUL(KC) = RGBFUL(KC) + CGLO
            ENDDO
*         End of this glow light
          ENDDO
          ENDIF
*
*         That does it for the shading computation.
*         ZS should still be a shadow-space co-ordinate of the pixel 
*         whose shading we were interested in, and zstop should be a
*         shadow-space object's co-ordinate no further than that from
*         the primary light source (modulo the empirical slop factor).
*
          IF (INDTOP.EQ.INDSTP) THEN
            TILE(1,I,J) = RGBFUL(1)
            TILE(2,I,J) = RGBFUL(2)
            TILE(3,I,J) = RGBFUL(3)
          ELSE IF (ZS+ZSLOP.GE.ZSTOP) THEN
          nzslop = nzslop + 1
            TILE(1,I,J) = RGBFUL(1)
            TILE(2,I,J) = RGBFUL(2)
            TILE(3,I,J) = RGBFUL(3)
          ELSE
            TILE(1,I,J) = RGBSHD(1)
            TILE(2,I,J) = RGBSHD(2)
            TILE(3,I,J) = RGBSHD(3)
          ENDIF
C
C     Fog processing added July 1998; moved into transp. loop 2010.
C       Note: Background fog is applied at beginning of this loop.
C     Should have glow lights brighten fog?
      IF (FOGTYPE .GE. 0) THEN
          FOGDIM = FOGGY( FOGLIM(2) - ZP )
          TILE(1,I,J) = (1.-FOGDIM)*TILE(1,I,J) + FOGDIM*FOGRGB(1)
          TILE(2,I,J) = (1.-FOGDIM)*TILE(2,I,J) + FOGDIM*FOGRGB(2)
          TILE(3,I,J) = (1.-FOGDIM)*TILE(3,I,J) + FOGDIM*FOGRGB(3)
      ENDIF

C
C       Transparency processing totally overhauled Feb 2001
C     The first pass is sufficient if top object is opaque.
        IF (NTRANSP .EQ. 0) GOTO 299
        IF (INDEPTH.EQ.1 .AND. iand(FLAG(INDTOP),TRANSP).EQ.0) GOTO 299
        IF (ITPASS.EQ.1) THEN
            ACHAN(I,J) = SBLEND
        ELSE
            ACHAN(I,J) = 1. - (1.-ACHAN(I,J))*(1.-SBLEND)
        ENDIF
        RGBLND(1)  = (1.-SBLEND)*RGBLND(1) + TILE(1,I,J)
        RGBLND(2)  = (1.-SBLEND)*RGBLND(2) + TILE(2,I,J)
        RGBLND(3)  = (1.-SBLEND)*RGBLND(3) + TILE(3,I,J)
C       CALL ASSERT(ITPASS.LE.INDEPTH,'Ran off end of INDEPTH')
        IF (ITPASS.GE.INDEPTH) THEN
            TILE(1,I,J) = RGBLND(1)
            TILE(2,I,J) = RGBLND(2)
            TILE(3,I,J) = RGBLND(3)
            ZTOP = ZP
        ELSE
            ZP        = ZLIST(INDEPTH-ITPASS)
            INDTOP    = INDLIST(INDEPTH-ITPASS)
            NORMAL(1) = NORMLIST(1,INDEPTH-ITPASS)
            NORMAL(2) = NORMLIST(2,INDEPTH-ITPASS)
            NORMAL(3) = NORMLIST(3,INDEPTH-ITPASS)
            ITPASS    = ITPASS + 1
            GOTO 255
        ENDIF
C     End of transparency processing
C
299   CONTINUE

C
300     CONTINUE
400     CONTINUE
*       do tile averaging and save output tile in outbuf
C     For now fold schemes 0 and 1 together; later split for efficiency?
        IF (SCHEME.LE.1) THEN
          K = (ITILE-1)*NOX
          DO 420 J = 1, NOY
          DO 415 I = 1, NOX
          K = K + 1
C         CALL ASSERT (K.LE.SIZE(OUTBUF,1),'k>outsiz')
C
          DO 410 IC = 1, 3
            ICK = 256. * SQRT(TILE(IC,I,J))
            IF (ICK.LT.0) ICK = 0
            IF (ICK.GT.255) ICK = 255
          IF (GAMMACORRECTION) ICK = GAMMA_MAP(ICK+1)
            OUTBUF(K,IC) = ICK
  410     CONTINUE
C
        IF (SCHEME.EQ.0) THEN
          ICK = 255. * ACHAN(I,J)
          IF (ICK.LT.0)   ICK = 0
          IF (ICK.GT.255) ICK = 255
          OUTBUF(K,4) = ICK
        END IF
C
415       CONTINUE
          K = K + NOX*(NTX-1)
420       CONTINUE
        ELSEIF (SCHEME.EQ.2) THEN
          K = (ITILE-1)*NOX
          DO 440 J = 1, NOY
          DO 435 I = 1, NOX
          K = K + 1
          DO 430 IC = 1, 3
*           I'm not quite convinced by this pixel averaging
*           (is a corner worth too much in this setup?):
            TMP = (TILE(IC,2*I-1,2*J-1) +
     &             TILE(IC,2*I  ,2*J-1) +
     &             TILE(IC,2*I-1,2*J  ) +
     &             TILE(IC,2*I  ,2*J  )) / 4.
            ICK = 256. * SQRT(TMP)
            IF (ICK.LT.0) ICK = 0
            IF (ICK.GT.255) ICK = 255
          IF (GAMMACORRECTION) ICK = GAMMA_MAP(ICK+1)
            OUTBUF(K,IC)   = ICK
430       CONTINUE
435       CONTINUE
          K = K + NOX*(NTX-1)
440       CONTINUE
        ELSEIF (SCHEME.EQ.3) THEN
          NHX = NOX/2
          NHY = NOY/2
          K = (ITILE-1)*NOX
          DO 460 J = 1, NHY
          DO 455 I = 1, NHX
          DO 450 IC = 1, 3
*           Bad pixel averaging?:
            TMP1 = (TILE(IC,3*I-2,3*J-2)    +
     &              TILE(IC,3*I-1,3*J-2)/2. +
     &              TILE(IC,3*I-2,3*J-1)/2. +
     &              TILE(IC,3*I-1,3*J-1)/4.) / (9./4.)
            TMP2 = (TILE(IC,3*I-1,3*J-2)/2. +
     &              TILE(IC,3*I  ,3*J-2)    +
     &              TILE(IC,3*I-1,3*J-1)/4. +
     &              TILE(IC,3*I  ,3*J-1)/2.) / (9./4.)
            TMP3 = (TILE(IC,3*I-2,3*J-1)/2. +
     &              TILE(IC,3*I-1,3*J-1)/4. +
     &              TILE(IC,3*I-2,3*J  )    +
     &              TILE(IC,3*I-1,3*J  )/2.) / (9./4.)
            TMP4 = (TILE(IC,3*I-1,3*J-1)/4. +
     &              TILE(IC,3*I  ,3*J-1)/2. +
     &              TILE(IC,3*I-1,3*J  )/2. +
     &              TILE(IC,3*I  ,3*J  )   ) / (9./4.)
            ICK1 = MIN(MAX(INT(256.*SQRT(TMP1)),0),255)
            ICK2 = MIN(MAX(INT(256.*SQRT(TMP2)),0),255)
            ICK3 = MIN(MAX(INT(256.*SQRT(TMP3)),0),255)
            ICK4 = MIN(MAX(INT(256.*SQRT(TMP4)),0),255)
          IF (GAMMACORRECTION) THEN
            ICK1 = GAMMA_MAP(ICK1+1)
            ICK2 = GAMMA_MAP(ICK2+1)
            ICK3 = GAMMA_MAP(ICK3+1)
            ICK4 = GAMMA_MAP(ICK4+1)
          ENDIF
            OUTBUF(   K+1,IC) = ICK1
            OUTBUF(   K+2,IC) = ICK2
            OUTBUF(NX+K+1,IC) = ICK3
            OUTBUF(NX+K+2,IC) = ICK4
450       CONTINUE
          K = K + 2
455       CONTINUE
          K = K + NOX*(2*NTX - 1)
460       CONTINUE
        ELSE
          CALL ASSERT(.FALSE.,'crash 500')
        ENDIF
500   CONTINUE
*     Ready to write when we have completed a row of tiles
      K = 0
      DO 550 J=1,NOY
      LINOUT = LINOUT + 1
      IF (LINOUT.GT.NAY) GOTO 600
      IERR = LOCAL(2, OUTBUF(K+1,1), OUTBUF(K+1,2), OUTBUF(K+1,3),
     &                    OUTBUF(K+1,4) )
        K = K + NX
C       CALL ASSERT (K.LE.SIZE(OUTBUF,1),'k>outsiz')
550   CONTINUE
600   CONTINUE
*
*     Report any soft failures
      IF (  NSXMAX.GT.0 .OR. NSYMAX.GT.0 
     & .OR. TRANOVFL.GT.0  ) THEN
      WRITE(NOISE,*)'   >>> WARNINGS <<<'
      END IF
      IF (NSXMAX.GT.0) WRITE(NOISE,*)
     &                '   Possible shadow error NSXMAX=',NSXMAX      
      IF (NSYMAX.GT.0) WRITE(NOISE,*)
     &                '   Possible shadow error NSYMAX=',NSYMAX      
      IF (TRANOVFL.GT.0) WRITE(NOISE,601)
     &  '   Transparency processing truncated at MAXTRANSP=',
     &  MAXTRANSP, '  for', TRANOVFL,' pixels'
601   FORMAT(A,I3,A,I10,A)
*
*     Debugging information
      if (verbose) then
      if (nzslop.gt.0) write(noise,*)'   NZSLOP failures=',nzslop
      if (nslow.gt.0)  write(noise,*)'   NSLOW  failures=',nslow
      endif
*
*     close up shop
      IERR = LOCAL(3)
*
      END

      SUBROUTINE TRANSF (X,Y,Z)
*     Input transformation
      COMMON /MATRICES/ XCENT, YCENT, SCALE, EYEPOS, SXCENT, SYCENT,
     &                  TMAT, TINV, TINVT, SROT, SRTINV, SRTINVT
     &                 ,RAFTER, TAFTER
      REAL   XCENT, YCENT, SCALE, EYEPOS, SXCENT, SYCENT
*     Transformation matrix, inverse of transpose, and transposed inverse
      REAL   TMAT(4,4), TINV(4,4), TINVT(4,4)
*     Shortest rotation from light source to +z axis
      REAL   SROT(4,4), SRTINV(4,4), SRTINVT(4,4)
*     Post-hoc transformation on top of original TMAT
      REAL   RAFTER(4,4), TAFTER(3)
      REAL   X,Y,Z
      REAL   G(4),H(4)
      H(1) = X*TMAT(1,1) + Y*TMAT(2,1) + Z*TMAT(3,1) + TMAT(4,1)
      H(2) = X*TMAT(1,2) + Y*TMAT(2,2) + Z*TMAT(3,2) + TMAT(4,2)
      H(3) = X*TMAT(1,3) + Y*TMAT(2,3) + Z*TMAT(3,3) + TMAT(4,3)
      H(4) = X*TMAT(1,4) + Y*TMAT(2,4) + Z*TMAT(3,4) + TMAT(4,4)
*     Apply post-hoc rotation and translation also
      G(1) = RAFTER(1,1)*H(1) + RAFTER(1,2)*H(2) + RAFTER(1,3)*H(3)
     &     + TAFTER(1)
      G(2) = RAFTER(2,1)*H(1) + RAFTER(2,2)*H(2) + RAFTER(2,3)*H(3)
     &     + TAFTER(2)
      G(3) = RAFTER(3,1)*H(1) + RAFTER(3,2)*H(2) + RAFTER(3,3)*H(3)
     &     + TAFTER(3)
      CALL ASSERT (H(4).NE.0.,'infinite vector')
      X = G(1) / H(4)
      Y = G(2) / H(4)
      Z = G(3) / H(4)
      RETURN
      END

      SUBROUTINE ISOLATE( X, Y )
*     Expand X and Y coordinates to fill image regardless of aspect ratio
      COMMON /OPTIONS/ FONTSCALE, GAMMA, ZOOM, NSCHEME, SHADOWFLAG, XBG,
     &                 NAX, NAY, OTMODE, QUALITY, INVERT, LFLAG
      REAL             FONTSCALE, GAMMA, ZOOM
      INTEGER          NSCHEME, SHADOWFLAG, XBG
      INTEGER*4        NAX, NAY, OTMODE, QUALITY
      LOGICAL*2        INVERT, LFLAG
      COMMON /MATRICES/ XCENT, YCENT, SCALE, EYEPOS, SXCENT, SYCENT,
     &                  TMAT, TINV, TINVT, SROT, SRTINV, SRTINVT
     &                 ,RAFTER, TAFTER
      REAL   XCENT, YCENT, SCALE, EYEPOS, SXCENT, SYCENT
      COMMON /NICETIES/ TRULIM,      ZLIM,    FRONTCLIP, BACKCLIP
     &                , ISOLATION
      REAL              TRULIM(3,2), ZLIM(2), FRONTCLIP, BACKCLIP
      INTEGER           ISOLATION
*
      IF (INVERT) Y = -Y
      IF (ISOLATION.EQ.2) THEN
          ASPECT = XCENT/YCENT
        IF (ASPECT.GT.1.0) X = X * ASPECT
        IF (ASPECT.LT.1.0) Y = Y / ASPECT
      ENDIF
      RETURN
      END



      SUBROUTINE HSORTD (N, A, NDEX)
      INTEGER N
      REAL   A(N)
      INTEGER NDEX(N)
*
*     this formulation of heapsort is based on n. wirth,
*     "algorithms + data structures = programs" (p. 75).
*
*     the caller supplies an array, a, containing n elements, and an
*     index array with space for n integers.
*     a and n are considered "read-only" by the subroutine, but ndex
*     is filled by the subroutine with the sequence of indices of a
*     that obtain the elements of a in ascending order.  this is
*     similar to the apl unary "tree" operator.  thus a(ndex(1)) is the
*     smallest element after the sort, and a(ndex(n)) is the largest.
*
      INTEGER L, R, T
      DO 10 I = 1, N
10    NDEX(I) = I
      L = N/2 + 1
      R = N
20    IF (L .LE. 1) GO TO 30
      L = L - 1
      CALL HSIFTD (N, A, NDEX, L, R)
      GO TO 20
30    IF (R .LE. 1) RETURN
      T = NDEX(1)
      NDEX(1) = NDEX(R)
      NDEX(R) = T
      R = R - 1
      CALL HSIFTD (N, A, NDEX, L, R)
      GO TO 30
      END
      SUBROUTINE HSIFTD (N, A, NDEX, L, R)
*     used by hsortd
      INTEGER N
      REAL   A(N)
      INTEGER NDEX(N), L, R
      INTEGER I, J, X
      I = L
      J = I + I
      X = NDEX(I)
10    IF (J .GT. R) GO TO 20
      IF (J .GE. R) GO TO 15
      IF (A(NDEX(J)) .LT. A(NDEX(J+1))) J = J + 1
15    IF (A(X) .GE. A(NDEX(J))) GO TO 20
      NDEX(I) = NDEX(J)
      I = J
      J = I + I
      GO TO 10
20    NDEX(I) = X
      RETURN
      END

      SUBROUTINE PLANER (X1,Y1,Z1,X2,Y2,Z2,X3,Y3,Z3, A,B,C,D)
      IMPLICIT REAL   (A-Z)
*     solve for coefficients of plane eqn z=ax+by+c
*     and yield d=0 in case of degenerate ("edge-on") triangle
      D1 = Z1*(Y2-Y3)       - Y1*(Z2-Z3)       + Z2*Y3-Y2*Z3
      D2 = X1*(Z2-Z3)       - Z1*(X2-X3)       + X2*Z3-Z2*X3
      D3 = X1*(Y2*Z3-Z2*Y3) - Y1*(X2*Z3-Z2*X3) + Z1*(X2*Y3-Y2*X3)
      D  = X1*(Y2-Y3)       - Y1*(X2-X3)       + X2*Y3-Y2*X3
      A = 0.
      B = 0.
      C = 0.
      IF (D.NE.0.) THEN
        A = D1/D
        B = D2/D
        C = D3/D
      ENDIF
      IF (ABS(A)+ABS(B)+ABS(C).GT.1E10) THEN
            D = 0.0
      ENDIF
      RETURN
      END

      SUBROUTINE ASSERT (LOGIC, DAMMIT)
      LOGICAL LOGIC
      CHARACTER*(*) DAMMIT
      INTEGER ASSOUT
      LOGICAL VERBOSE
      COMMON /ASSCOM/ ASSOUT, VERBOSE
      SAVE /ASSCOM/
      IF (LOGIC) RETURN
      WRITE (ASSOUT,*) 'ERROR >>>>>> ',DAMMIT
      CALL EXIT(-1)
      END

C
C Find Z coord of point on surface of cylinder with known X and Y coords
C cylinder axis is X2 - X1, cylinder radius is R
C Need to find Z coord ZB.
C flag is 0 if cylinder had rounded ends, FLAT if it has flat ends,
C Also find nearest point XYZA on cylinder axis and fraction along it.
C
      FUNCTION CYLTEST  ( flag, axfrac,
     &                   x1,y1,z1,  x2,y2,z2,  xb,yb,zb,  R,  xa,ya,za )
      LOGICAL  CYLTEST
c     implicit NONE
*
*     Bit definitions for FLAG array
      INTEGER    FLAT,      RIBBON,    SURFACE,   PROPS
      PARAMETER (FLAT=2,    RIBBON=4,  SURFACE=8, PROPS=16)
      INTEGER    TRANSP,    HIDDEN,    INSIDE,    MOPT1
      PARAMETER (TRANSP=32, HIDDEN=64, INSIDE=128,MOPT1=256)
      INTEGER      VCOLS,     CLIPPED,      VTRANS,      BOUNDED
      PARAMETER   (VCOLS=512, CLIPPED=1024, VTRANS=2048, BOUNDED=4096)
      INTEGER      TFI
      PARAMETER (TFI = TRANSP + INSIDE)
c
      INTEGER flag
      REAL        x1,y1,z1, x2,y2,z2, xb,yb,zb
      REAL        R
      REAL        xa,ya,za
c
      REAL        ca,cb,cg,dx,dy,dz,d2
      REAL        A0,A1,A2,Q
      REAL        p1,r2,dx2,dy2
      REAL        dd1,dd2
c
c     start with direction cosines * d2
      ca = x2 - x1
      cb = y2 - y1
      cg = z2 - z1
c
c     other useful quantities
c     (note: if d2==0 must be degenerate cylinder, really a disk)
      r2 = R*R
      d2 = ca*ca + cb*cb + cg*cg
      dx = xb - x1
      dy = yb - y1
      dx2 = dx**2
      dy2 = dy**2
c
c     use these to find coefficients of quadratic equation for ZB
c     EAM Jan 1997 test and handle dx-dy=0
      if (ca.eq.0. .and. cb.eq.0.) then
          if (z2.gt.z1) p1 =  1.0
          if (z2.lt.z1) p1 = -1.0
          goto 100
      end if
c
      A0 = (dx*cb - dy*ca)**2 + (dy2 + dx2)*cg*cg - r2*d2
      A1 = -2.0 * (dy*cg*cb + dx*ca*cg)
      A2 = ca*ca + cb*cb
      Q  = A1*A1 - 4.0*A0*A2
      if (Q .lt. 0) then
C           zb = -99999.
            cyltest = .false.
            return
      else
            if ( iand(flag,TFI) .eq. TFI) then
                  dz = (-sqrt(Q) - A1) / (2.0 * A2)
            else
                  dz = ( sqrt(Q) - A1) / (2.0 * A2)
            endif
            zb = z1 + dz
      end if
c
c     now find nearest point on cylinder axis
c     p1 is fraction along axis from x1 to x2
c     0 < p1 < 1 means point is on wall of cylinder
c
      dd1 = dx2 + dy2 + dz*dz
      dd2 = (x2-xb)**2 + (y2-yb)**2 + (z2-zb)**2
c
      p1 = (dd1 - r2) / d2
      if (p1 .le. 0.0) then
            p1 = 0.0
      else
            p1 = sqrt(p1)
      end if
c
      if ((dd2 .gt. (d2+r2)) .and. (dd2 .gt. dd1)) p1 = -p1
c
      if (p1 .ge. 0 .and. p1 .le. 1.0) then
            xa = p1*ca + x1
            ya = p1*cb + y1
            za = p1*cg + z1
            cyltest = .true.
            return
      end if
c
c     point is either on end cap, or missed entirely
c
  100 continue
      if (p1 .ge. 1.0) then
            xa = x2
            ya = y2
            za = z2
            dx = xb - x2
            dy = yb - y2
            dx2 = dx**2
            dy2 = dy**2
      else if (p1 .le. 0.0) then
            xa = x1
            ya = y1
            za = z1
      end if
c
c Rounded cylinder end
      if (iand(flag,FLAT) .eq. 0) then
            if (dx2+dy2 .gt. r2) then
                  cyltest = .false.
                  return
            else
                if ( iand(flag,TFI) .eq. TFI) then
                  zb = za - sqrt(r2 - (dx2+dy2))
                else
                  zb = za + sqrt(r2 - (dx2+dy2))
                end if
            end if
C
C Flat cylinder end
C
      else
            if (cg .eq. 0.) then
C               zb = -99999.
                cyltest = .false.
                return
            endif
            zb = (cg*za - ca*dx - cb*dy) / cg
            if (dx2 + dy2 + (zb-za)**2 .ge. r2) then
C               zb = -99999.
                cyltest = .false.
                return
            endif
            if (p1 .ge. 1.0) then
                xa = xb - (x2 - x1)
                ya = yb - (y2 - y1)
                za = zb - (z2 - z1)
            else if (p1 .le. 0.0) then
                xa = xb - (x1 - x2)
                ya = yb - (y1 - y2)
                za = zb - (z1 - z2)
            endif
      end if
c
      if (p1.gt.1.0) p1 = 1.0
      if (p1.lt.0.0) p1 = 0.0
      axfrac = p1
c
      cyltest = .true.
      return
      end

C Bookkeeping for transparency
C 5-Mar-2001
C New version that is not limited to 3-deep transparent objects.
C On exit:
C     INDTOP, ZTOP   contain the top object so far, and its height
C     ZHIGH          height of top opaque object
C
      subroutine rank( ind, zp, normal, flag )
*
      implicit NONE
      integer  ind
      real     zp, normal(3)
      integer*4  flag(1)
*
      integer  i,j,k
*
      INTEGER ASSOUT
      LOGICAL VERBOSE
      COMMON /ASSCOM/ ASSOUT, VERBOSE
*
      INCLUDE 'parameters.incl'
*
*     Support for transparency
      COMMON /TRANS/ NTRANSP, INDEPTH, INDTOP, TRANOVFL, ZTOP, ZHIGH, 
     &               INDLIST(MAXTRANSP), ZLIST(MAXTRANSP), 
     &                 NORMLIST(3,MAXTRANSP)
      INTEGER        NTRANSP, INDEPTH, INDTOP, INDLIST, TRANOVFL
      REAL           ZTOP, ZHIGH, ZLIST, NORMLIST
*
*     Bit definitions for FLAG(MAXOBJ) array
      INTEGER    FLAT,      RIBBON,    SURFACE,   PROPS
      PARAMETER (FLAT=2,    RIBBON=4,  SURFACE=8, PROPS=16)
      INTEGER    TRANSP,    HIDDEN,    INSIDE,    MOPT1
      PARAMETER (TRANSP=32, HIDDEN=64, INSIDE=128,MOPT1=256)
      INTEGER      VCOLS,     CLIPPED,      VTRANS,      BOUNDED
      PARAMETER   (VCOLS=512, CLIPPED=1024, VTRANS=2048, BOUNDED=4096)
*
*     The MOPT1 flag signals an alternative mode of transparency
      if (iand(flag(ind),mopt1).ne.0) goto 400
*
      do i = 1, indepth
                if (zp.gt.zlist(i)) then
            if (iand(flag(ind),transp).eq.0) then
                indepth    = i
                zhigh = zp
                goto 345
            else
                do j = indepth, i, -1
                  indlist(j+1)   = indlist(j)
                  zlist(j+1)     = zlist(j)
                  normlist(1,j+1)= normlist(1,j)
                  normlist(2,j+1)= normlist(2,j)
                  normlist(3,j+1)= normlist(3,j)
                enddo
                indepth = indepth + 1
                goto 344
            endif
          else if (iand(flag(indlist(i)),TRANSP).eq.0) then
            return
          endif
      enddo
c     If the rest of the list is transparent, add this at the end
      indepth = indepth + 1
      i = indepth
c
  344 continue
      if (indepth.ge.MAXTRANSP) then
          indepth  = MAXTRANSP - 1
          tranovfl = tranovfl + 1
          return
      endif
c
  345 continue
      indlist(i) = ind
      zlist(i)   = zp
      normlist(1,i) = normal(1)
      normlist(2,i) = normal(2)
      normlist(3,i) = normal(3)
      if (iand(flag(ind),transp).eq.0) zhigh = zp
      indtop = indlist(1)
      ztop   = zlist(1)
      return
c
c     MOPT1 version. Same routine, except this time we have the extra
c     overhead of having to check for duplication of material.
  400 continue
      do i = 1, indepth
                if (zp.gt.zlist(i)) then
              if (flag(ind)/65536 .eq. flag(indlist(i))/65536) then
                goto 345
            endif
c           Handle case where two MOPT1 surfaces have intervening transp obj
c           In this case overwrite the lower MOPT1 surface
            do k = i, indepth-1
               if (flag(ind)/65536 .eq. flag(indlist(k+1))/65536) 
     &                  goto 401
            enddo
            if (indepth .ge. MAXTRANSP-1) then
                tranovfl = tranovfl + 1
            else
                k = indepth
                indepth = indepth + 1
            endif
  401       continue
            do j = k, i, -1
               indlist(j+1)   = indlist(j)
               zlist(j+1)     = zlist(j)
               normlist(1,j+1)= normlist(1,j)
               normlist(2,j+1)= normlist(2,j)
               normlist(3,j+1)= normlist(3,j)
            enddo
            goto 345
          else if (iand(flag(indlist(i)),TRANSP).eq.0) then
            return
          else if (flag(ind)/65536 .eq. flag(indlist(i))/65536) then
            return
          endif
      enddo
c     If the rest of the list is transparent, add this at the end
      indepth = indepth + 1
      i = indepth
      goto 344
c
      end
c

      SUBROUTINE CHKRGB( RED, GRN, BLU, MESSAGE )
      REAL RED, GRN, BLU
      CHARACTER*(*) MESSAGE
        CALL ASSERT (RED.GE.0., MESSAGE)
        CALL ASSERT (GRN.GE.0., MESSAGE)
        CALL ASSERT (BLU.GE.0., MESSAGE)
        CALL ASSERT (RED.LE.1., MESSAGE)
        CALL ASSERT (GRN.LE.1., MESSAGE)
        CALL ASSERT (BLU.LE.1., MESSAGE)
      RETURN
      END


      FUNCTION FOGGY( DEPTH )
      REAL     FOGGY, DEPTH
      COMMON /FOGCOM/ FOGTYPE,FOGFRONT,FOGBACK,FOGDEN,FOGLIM,FOGRGB
      INTEGER FOGTYPE
      REAL    FOGFRONT, FOGBACK, FOGDEN, FOGLIM(2), FOGRGB(3)
      REAL    FOGDIM
c
      FOGDIM = 0.5
      IF (FOGTYPE .EQ. 0) 
     &          FOGDIM = FOGDEN * DEPTH / (FOGLIM(2)-FOGLIM(1))
      IF (FOGTYPE .GT. 0)
     &          FOGDIM = 1. - EXP(-FOGDEN * DEPTH/(FOGLIM(2)-FOGLIM(1)))
      FOGDIM = MAX( 0.0, FOGDIM )
      FOGDIM = MIN( 1.0, FOGDIM )
      FOGGY  = FOGDIM
      RETURN
      END

      FUNCTION DET( A )
      REAL     DET, A(4,4)
      DET = A(1,1)*A(2,2)*A(3,3) + A(1,2)*A(2,3)*A(3,1) 
     &          + A(2,1)*A(3,2)*A(1,3) - A(1,1)*A(2,3)*A(3,2)
     &          - A(3,3)*A(1,2)*A(2,1) - A(1,3)*A(2,2)*A(3,1)
      RETURN
      END


      subroutine liblookup( name, fullname )
c
      character*(*) name
      character*132 fullname
      character*132 R3DLIB
c
      call getenv('R3D_LIB',R3DLIB)
c
      fullname = ' '
      len = 0
      do i = 1, 132
          if (R3DLIB(i:i).ne.' ') len = i
      enddo
      if (len.eq.0) then
          fullname = name
          return
      else
          fullname(1:len) = R3DLIB(1:len)
          fullname(len+1:len+1) = '/'
          j = len+2
      endif
c
      len = 0
  100 continue
      len = len + 1
      if (name(len:len).ne.' ') goto 100
      len = len - 1
c
      fullname(j:j+len-1) = name(1:len)
      len = j+len-1
c
      return
      end

C EAM - 21 Feb 2001     Version 2.6
C Test for bounding planes
C Returns .FALSE. if the point does not have to be rendered
C         .TRUE.  if the point is to be rendered, in which case
C                 ZP, DX, DY, DZ have been updated if on bounding plane
C             ZP = height of point being rendered
C             DX,DY,DZ = surface normal at point (XP,YP,ZP)
C             BPIND is the object number of the bounding plane 
C                   or 0 if the bounding planes don't trigger
C
      FUNCTION INBOUNDS( MAT, TYPE, LIST, DETAIL, 
     &                         XP,YP,ZP, DX,DY,DZ, ZBACK, BPIND )
      IMPLICIT NONE
      LOGICAL  INBOUNDS
      INTEGER             MAT, M, TYPE(1), LIST(1), BPIND
      REAL                DETAIL(1), XP,YP,ZP, DX,DY,DZ, ZBACK
c
c     Object type used for bounding planes (may change)
      INTEGER    INTERNAL
      PARAMETER (INTERNAL = 4)
c
      REAL BPLANE(3), BPNORM(3)
      REAL XN,YN,ZN, TEMP
      INTEGER MIND
      LOGICAL TESTBACK
c
      INBOUNDS = .FALSE.
      BPIND    = 0
      TESTBACK = (ZBACK.NE.ZP)
      M = MAT
      DO WHILE (TYPE(M).EQ.INTERNAL)
        MIND = LIST(M)
        BPLANE(1) = DETAIL(MIND+3)
        BPLANE(2) = DETAIL(MIND+4)
        BPLANE(3) = DETAIL(MIND+5)
        BPNORM(1) = DETAIL(MIND+6)
        BPNORM(2) = DETAIL(MIND+7)
        BPNORM(3) = DETAIL(MIND+8)
        XN = XP - BPLANE(1)
        YN = YP - BPLANE(2)
        ZN = ZP - BPLANE(3)
        TEMP = XN*BPNORM(1) + YN*BPNORM(2) + ZN*BPNORM(3)
        IF (TEMP.GT.0) THEN
          IF (TESTBACK) THEN
            ZP = ZBACK
            ZN = ZP - BPLANE(3)
            TEMP = XN*BPNORM(1) + YN*BPNORM(2) + ZN*BPNORM(3)
            IF (TEMP.GT.0) RETURN
          ENDIF
          DX = BPNORM(1)
          DY = BPNORM(2)
          DZ = BPNORM(3)
          ZP = (XN*DX + YN*DY) / (-DZ) + BPLANE(3)
          BPIND = M
        ENDIF
      M = M + 1
      ENDDO
      INBOUNDS = .TRUE.
      RETURN
      END
C
C Equivalent test for bounding planes in ORTEP mode (only the octant clipped
C by all three bounding planes is removed). Only intended for ellipsoids.
C This tests the AND (rather than the OR) of multiple bounding planes.
C
      SUBROUTINE ORTEPBOUNDS( MAT, TYPE, LIST, DETAIL, 
     &                              XP,YP,ZP, DX,DY,DZ, ZBACK, BPIND )
      IMPLICIT NONE
      INTEGER             MAT, M, TYPE(1), LIST(1), BPIND
      REAL                DETAIL(1), XP,YP,ZP, DX,DY,DZ, ZBACK
c
c     Object type used for bounding planes (may change)
      INTEGER    INTERNAL
      PARAMETER (INTERNAL = 4)
c
      REAL BPLANE(3), BPNORM(3)
      REAL XN,YN,ZN, TEMP, ZMAX, DXMAX, DYMAX, DZMAX
      INTEGER MIND
c
      M     = MAT
      ZMAX  = -1.E10
      DXMAX = DX
      DYMAX = DY
      DZMAX = DZ
      DO WHILE (TYPE(M).EQ.INTERNAL)
        MIND = LIST(M)
        BPLANE(1) = DETAIL(MIND+3)
        BPLANE(2) = DETAIL(MIND+4)
        BPLANE(3) = DETAIL(MIND+5)
        BPNORM(1) = DETAIL(MIND+6)
        BPNORM(2) = DETAIL(MIND+7)
        BPNORM(3) = DETAIL(MIND+8)
        XN = XP - BPLANE(1)
        YN = YP - BPLANE(2)
        ZN = ZP - BPLANE(3)
        TEMP = XN*BPNORM(1) + YN*BPNORM(2) + ZN*BPNORM(3)
        IF (TEMP.LT.0) THEN
            BPIND = 0
            RETURN
        ENDIF
        TEMP = (XN*BPNORM(1)+YN*BPNORM(2)) / (-BPNORM(3)) + BPLANE(3)
        IF (TEMP.GT.ZMAX) THEN
            ZMAX = TEMP
            DXMAX = BPNORM(1)
            DYMAX = BPNORM(2)
            DZMAX = BPNORM(3)
            BPIND = M
        ENDIF
        M = M + 1
      ENDDO
      ZP = ZMAX
      DX = DXMAX
      DY = DYMAX
      DZ = DZMAX
      CALL ASSERT(DZ.GT.0,'ORTEP bounds incorrectly initialized')
      RETURN
      END

        SUBROUTINE GET_TRY(NOLD, NNEW, TRY, DIMENS)
*
*       NOLD = old memory allocation
*       NNEW = new memory allocation needed
*       TRY = new memory allocations beyond that needed to be attempted,
*             so that we don't have to expand the array too often.
*       DIMENS = number of dimensions to be expanded (1 to 3)
*
*       For DIMENS=1: if the new allocation is less than twice the old 
*       try to double the old allocation; otherwise add 50% to the new.
*       If that fails, add 50% to the old or 25% of the new;
*       If that fails, just try the new, then give up.
*
*       For DIMENS = 2, expand by the square root of 2, 1.5 or 1.25
*       for the same effect on total memory in 2-D arrays where bot
*       dimensions are expanded. 
*       For DIMENS = 3 (not currently used) use the cube root.
*
        INTEGER NOLD, NNEW, TRY(3), N, DIMENS
        REAL RATIO(3, 3)
        DATA RATIO / 2.0,  1.5,  1.25,
     &               1.41, 1.22, 1.12,
     &               1.26, 1.14, 1.07 /
        N = INT( FLOAT(NOLD) * RATIO(1, DIMENS) )
        if (NNEW.LT.N) THEN
            TRY(1) = N
            N = INT( FLOAT(NOLD) * RATIO(2, DIMENS) )
            if (NNEW.LT.N) THEN
                  TRY(2) = N
              else
                  TRY(2) = INT( FLOAT(NNEW) * RATIO(3, DIMENS) )
              ENDIF
          else
              TRY(1) = INT( FLOAT(NNEW) * RATIO(2, DIMENS) )
              TRY(2) = INT( FLOAT(NNEW) * RATIO(3, DIMENS) )
          endif
          TRY(3) = NNEW
          RETURN
          END 

Generated by  Doxygen 1.6.0   Back to index