File: render.f

package info (click to toggle)
raster3d 3.0-3-5
links: PTS, VCS
area: main
in suites: buster
size: 5,172 kB
sloc: fortran: 9,536; ansic: 1,060; makefile: 319; sh: 258; csh: 15
file content (5096 lines) | stat: -rw-r--r-- 173,777 bytes
      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
* EAM Apr 2011	- Ignore blank lines when looking for next object
*		  
*     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 doesn''t 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
c	April 2011 - ignore blank line
	READ (LINE,*,END=7,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	V3.0.3:   Render both sides of opaque triangles
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.NE.0 .AND. 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 (don''t 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