1 MODULE LISTS 2 INTEGER, ALLOCATABLE, DIMENSION(:,:) :: KOUNT, MOUNT 3 INTEGER, ALLOCATABLE, DIMENSION(:,:) :: TTRANS 4 INTEGER ISTRANS 5 END MODULE LISTS 6 7 PROGRAM RENDER 8 USE LISTS 9* 10* Version 3.0 (14 Dec 2010) 11* 12* EAM May 1990 - add object type CYLIND (cylinder with rounded ends) 13* and CYLFLAT (cylinder with flat ends) 14* - use Phong shading on triangles if they are sequentially 15* adjoining. 16* EAM Feb 1991 - port to Ultrix (minor changes) and modify output format 17* to depend on code in separate module "local". 18* EAM Mar 1993 - fix embarrassingly stupid bug in cylinder shadowing 19* add object type PLANE (triangle with infinite extent) 20* EAM Nov 1993 - Version 2.0 (beta test) 21* fix bug which allowed objects to shadow themselves 22* Command line options for 3 output modes 23* EAM Apr 1994 - Version 2.0 release 24* TIFF output support in local.c 25* minor changes to fortran source to make IBM xlf compiler happy 26* EAM Sep 1994 - EXPERIMENTAL VERSION WITH OBJECT TYPES 7, 8, 9 27* EAM Jan 1995 - move DATA statement to make linux happy 28* EAM Mar 1995 - fix bug in routine CYLTEST which took bites out of cylinder ends 29* EAM May 1995 - fold object types 7/8/9 back into distributed code for V2.1 30* EAM Jan 1996 - Version 2.2 31* Add code for transparency and faster MATERIAL bookkeeping 32* Also fix major problems with explicit surface normals 33* object type 8 expanded to describe transparency 34* EAM May 1996 - antialiasing scheme 4, file indirection, 35* minor changes to accommodate HPUX 36* EAM Oct 1996 - trap and forgive shadowing error due to too small NSX,NSY 37* EAM Nov 1996 - scheme 0 causes alpha blend channel in output image 38* per-tile count of transparent objects 39* EAM Jan 1997 - zero length cylinders treated as spheres 40* Object types 10 + 11 (fonts and labels) accepted but ignored 41* Object type 12 reserved for other label information 42* Material properties can override object colors 43* - Material OPT(1) = 1 transparency option 44* OPT(4) = continuation lines for more material properties 45* EAM Mar 1997 - Make SLOP larger, and dependent on tile size 46* - GLOW light source specified by object type 13 47* EAM May 1997 - V2.3c allow multiple glow lights; make cyltest a function; 48* V2.3d remove DATA statements; more terse output; 49* BACKFACE material option; EYEPOS = 0.0 disables perspective 50* EAM Jul 1997 - add commons LISTS MATRICES NICETIES RASTER 51* - D_LINES code for quadric surfaces, ISOLATION 52* EAM Sep 1997 - fix normals of flat cylinder ends (thanks to Takaaki Fukami) 53* EAM Nov 1997 - add VERTEXRGB object type to extend triangle descriptions, 54* allow # as comment delimiter in input stream 55* EAM Dec 1997 - Release V2.4b 56* EAM Feb 1998 - fixed bug in check against limiting radius of quadrics 57* EAM Jul 1998 - Object type 16 (GLOBAL PROPERTIES); fog 58* EAM Aug 1998 - render back side of transparent flat-ended cylinders 59* by duplicating object with INSIDE flag bit set. 60* EAM Oct 1998 - check environmental variable R3D_LIB for input file indirection 61* - V2.4g includes some preliminary code to support Z-clipping 62* EAM Nov 1998 - more work on Z-clipping 63* EAM Feb 1999 - re-work output module local.c to support -jpeg and -out 64* EAM May 1999 - allow explicit vertex colors for cylinders also 65* EAM Jul 1999 - 2.4l preliminary work towards an after-the-fact rotation option 66* EAM Sep 1999 - 2.5a command line parsing in separate routine parse() 67* label processing folded into render; routines in r3dtops.f 68* EAM Jan 2000 - 2.5b general release 69* EAM Feb 2000 - 2.5c (bug fixes to 2.5b) 70* EAM Mar 2000 - object types VTRANSP (18) and ISOLATE2 (19) 71* EAM Sep 2000 - 2.5e uncompress indirect files ending with .Z or .gz 72* discard BACKCLIP objects on input, implement BACKCLIP material 73* EAM Nov 2000 - 2.5f more bug-fixes to rotation of surface normals 74* EAM Feb 2001 - 2.5g command line shadows 75* EAM Mar 2001 - V2.6 (alpha test) 76* bounding planes 77* revamped transparency code, remove limit of 2 stacked objects 78* break out shared maximum dimensions to paramters.incl 79* make back surface HIDING (former INMODE=4) the default for 80* non-bounded opaque triangles 81* EAM Jul 2001 - V2.6b first release 82* EAM Aug 2001 - V2.6c bug-fix for MOPT1 processing 83* PNG output format ( -DPNG_SUPPORT ) 84* EAM Feb 2002 - allow a little RIBBONSLOP in testing for ribbon triangles 85* fix up complicated corner cases in bounded surface algorithm 86* EAM Apr 2002 - clean up auto-tiling, and allow zero NPX or NPY to trigger it 87* EAM Apr 2006 - V2.6d gfortran accommodations 88* - Change AND() to iand() everywhere 89* - Change OR() to ior() everywhere 90* EAM Mar 2008 - initialize various static storage areas 91* FZ Dec 2009 - initialize more static storage areas (valgrind runs) 92* FZ Feb 2010 - expandable memory allocation 93* JMK Dec 2010 - Bugfix for JUSTCLIPPED reported by Joe Krahn 94* JMK Dec 2010 - More transparency algorithms; use OPT[2] to pass choice 95* EAM Dec 2010 - Read in all header info before initializing output file 96* EAM Apr 2011 - Ignore blank lines when looking for next object 97* 98* General organization: 99* 100* - read in control parameters and initial output image file 101* - read in list of objects 102* - count objects that may impinge on each tile 103* - do this for both pixel and rotated "shadow" space 104* - sort objects 105* - go through main object list in sorted order 106* - fill in short lists of objects 107* - repeat the sort etc. for the objects in shadow space 108* - that's it for the "cheap" part 109* - for each tile: 110* - for each pixel: 111* - search objects to find highest point for pixel 112* - if it's transparent find the next one down as well 113* - transform resulting (x,y,z) to shadow space 114* - find closest z' for new x',y' 115* - this tells you if the pixel is in shadow or not 116* - shade accordingly 117* - copy tile to output buffer 118* 119* 120* Easy-to-change constants (kept in file parameters.incl): 121* 122* - maximum size of any expandable array MAXMEM 123* - maximum number of tiles in each direction MAXNTX*, MAXNTY* 124* - number of shadow tiles in each direction NSX, NSY 125* - maximum number of pixels in a tile MAXNPX*, MAXNPY* 126* - maximum number of objects MAXOBJ* 127* - maximum number of material specifications MAXMAT* 128* - maximum depth of stack transparent objects MAXTRANSP 129* 130* *NOTE: MAXNTX,MAXNTY,MAXNPX,MAXNPY,MAXOBJ and MAXMAT are now 131* initial array sizes for expandable arrays, rather than limits. 132* 133* Input (line by line except where noted): 134* 135* - TITLE anything you like 136* - NTX,NTY tiles in each direction 137* - NPX,NPY pixels per tile to compute in each direction 138* - SCHEME pixel averaging scheme (1, 2, or 3) 139* - 0 no anti-aliasing, use alpha channel 140* - 1 no anti-aliasing, no alpha channel 141* - 2 means 2x2 computing pixels for 1 output pixel 142* - 3 means 3x3 computing pixels for 2x2 output pixels 143* - 4 same as 3, but NTX,NTY expanded inside program 144* - BKGND background colour (r,g,b in range 0 to 1) 145* - SHADOW "shadow mode?" (T or F) 146* - IPHONG Phong power (e.g., 20) 147* - STRAIT straight-on (2ndary) light component (e.g., 0.1) 148* - AMBIEN ambient light component (e.g., 0.05) 149* - SPECLR specular reflection component (e.g., 0.30) 150* - EYEPOS eye position along +z coordinate (e.g., 4) 151* - relative to 1=narrow dimension of screen 152* - used for perspective, EYEPOS = 0.0 disables perspective 153* - SOURCE main light source position (x,y,z components) 154* - vector length ignored, point source is at infinity 155* - TMAT global transformation on input objects 156* - postfix 4x4 matrix on 4 lines, as you would write it 157* - upper left 3x3 must be pure rotation 158* - lower left 1x3 is translation 159* - lower right 1x1 is global scaling (reduction) 160* - upper right 3x1 causes extra perspective (should be 0) 161* - applies to homogeneous co-ordinates (x,y,z,1) 162* - radii are only scaled down by global scaling TMAT(4,4) 163* - INMODE object input mode (1, 2, or 3) 164* - mode 1: all objects are triangles 165* - mode 2: all objects are spheres 166* - mode 3: each object will be preceded by type 167* - INFMT or INFMTS object input format(s), 1 per line 168* - one format for modes 1 and 2, or three for mode 3 169* - each is fortran format in parentheses, or single * 170* - for 3 formats, the order of formats and details is: 171* - triangle: x1,y1,z1,x2,y2,z2,x3,y3,z3,r,g,b 172* - sphere: x,y,z,radius,r,g,b 173* - trcone: x1,y1,z1,rad1,x2,y2,z2,rad2,r,g,b 174* - cylinder: as truncated cone, but 2nd radius ignored 175* - objects 176* - modes 1,2: each object starts on a new line 177* - read according to the single format given 178* - mode 3: each object is preceded by a line giving type 179* - type 1: triangle (to be read with 1st format) 180* - type 2: sphere (to be read with 2nd format) 181* - type 3: cylinder with rounded ends (3rd format) 182* - type 4: [not implemented: truncated cone] 183* - type 5: cylinder with flat ends (3rd format) 184* - type 6: plane (=triangle with infinite extent) (1st format) 185* - type 7: normal vectors for previous triangle (1st format) 186* - type 8: material definition which applies to subsequent objects 187* - type 9: end previous material 188* - type 10: font selection (ignored in render) 189* - type 11: label (ignored other than to count them) 190* - type 12: (reserved for additional label processing) 191* - type 13: glow light source 192* - type 14: quadric surface (usually an ellipsoid) 193* - type 15: disable coordinate transformation of subsequent objects 194* - type 16: global properties (e.g. FOG) 195* - type 17: RGB triple for each vertex of preceding object 196* - type 18: transparency at each vertex of preceding object 197* - type 19: variant of type 15; forces unitary coordinate sytem 198* - type 0: no more objects (equivalent to eof) 199* 200*----------------------------------------------------------------------------- 201* EAM Sep 1994 202* 203* 1) Object type 7 signals an extra record giving explicit vertex normals 204* for a single triangle. This extra record must directly follow the 205* corresponding triangle and uses the same format. 206* 207* 2) Object type 8 signals an extra record giving extra or more explicit 208* material properties object. Current (trial) contents of record are: 209* MPHONG - overrides global Phong power for specular reflections 210* MSPEC - overrides global specular scattering contribution 211* SR,SG,SB - red/green/blue components for specular highlighting 212* (values <0 cause highlights to match object colour) 213* CLRITY - 0.0 (opaque) => 1.0 (transparent) 214* CLROPT - [reserved] suboptions for transparency handling 215* OPT2 - [reserved] suboptions for bounding planes Is it really used? 216* OPT3 - [reserved] 217* OPT4 - # of additional modifier records immediately following 218* These material properties remain in effect for subsequent objects 219* until object type 9 appears to terminate the effect. 220* 221* 3) Object type 9 terminates effect of previous materials property 222* 223*----------------------------------------------------------------------------- 224* Object types 10 and 11 are used for specifying labels. 225* Label object types are 226* - type 10: Font_Name size alignment 227* - type 11: XYZ RGB on first line 228* label (ascii characters enclosed in quotes) on second line 229* Object type 12 is reserved to go with this, as I have a nagging 230* suspicion more information may turn out to be necessary. 231*----------------------------------------------------------------------------- 232* Object type 13 specifies a "glow" light source; i.e. a non-shadowing 233* light source with finite origin GLOWSRC and illumination range GLOWRAD. 234* Specular highlights for this source specified by GLOWCOL and GPHONG. 235* 0.0 < GLOW < 1.0 = contribution of glow to total lighting model. 236* 13 237* GLOWSRC(3) GLOWRAD GLOW GOPT GPHONG GLOWCOL(3) 238*----------------------------------------------------------------------------- 239* V2.4 240* Object type 14 specifies a quadric surface 241* QQ = Ax^2 + By^2 + Cz^2 + 2Dxy + 2Eyz + 2Fxz + 2Gx + 2Hy + 2Iz + J 242* centered at XC,YC,ZC and truncated at bounding radius RC. Supporting 243* code is in file quadric.f 244* 14 245* XC YC ZC RC RED GRN BLU 246* A B C D E F G H I J 247* 248* Object type 15 is a single line signaling that subsequent objects are 249* not to be transformed by the TMAT matrix in the header. This isolation 250* from TMAT is terminated by an end material record (object type 9). 251*----------------------------------------------------------------------------- 252* V2.6 253* Object type 4 is used internally to implement bounding planes. 254* The BOUNDING_PLANE definition is given in the input stream as 255* a modifier to a MATERIAL descriptor. At the time it is read in 256* render creates a new object of type 4 to hold the bounding 257* plane parameters and modifiers. 258* DETAIL(K+...) is loaded with the following parameters: 259* SUBTYPE, BPTYPE, X,Y,Z, XNORM, YNORM, ZNORM, RED, GREEN, BLUE 260*----------------------------------------------------------------------------- 261* 262* Object space convention: 263* 264* - this is the space TMAT is supposed to map your data to 265* - centre of "virtual screen" is (0,0,0) 266* - x to right, y to top, z towards viewer 267* - the smaller of the x and y dimensions goes from -.5 to +.5 268* - z cuts off at +1 and -1 by default, but is modified by FRONTCLIP, BACKCLIP 269* - shadow box dimensions determined by NSX/NTX, NSY/NTY 270* 271*----------------------------------------------------------------------------- 272* David Bacon's comments from Version 1.0 273* 274* Bugs: 275* 276* - perspective is applied to raw objects, giving wrong lighting 277* - perspective unity factor is always at z=0 in object space 278* - shadow box doesn't necessarily enclose entire view prism 279* - some ASSERT calls are commented out for extra speed 280* - SLOP parameter is empirical fix to imprecision in shadowing 281* 282* Deferred priorities: 283* 284* - superior pixel averaging (you should use another pass?) 285* - better assignment of triangles to tiles (do clipping?) 286* - better estimate of max. triangle elevation within tile 287* 288* 289* Why I don't do shadowing properly: 290* 291* Although you might not notice it as a casual observer, 292* the main light source is (in effect) in different places 293* for different parts of the picture. More precisely, 294* perspective is applied to the objects comprising the 295* picture first, and THEN the lighting from a distant 296* light source is applied. The lighting would be correct 297* if there were no perspective, because the angle doesn't 298* change across the picture. With a scene in perspective, 299* however, the angle of the light beam, from the eye's 300* viewpoint, should vary a little. 301* The reason I allow this error is because the use of "tiles" 302* is implicitly a parallel projection, so I have to apply the 303* perspective to the objects initially. 304* Conceivably I could "undo" this whenever taking the light 305* source point of view (this would result in using a slightly 306* different light source position for different parts of the 307* picture), but that would cause another problem: 308* Perspective distorts spheres differently depending on whether 309* you ask about each point on the surface individually or use one 310* perspective factor for all points based on the centre. Consider 311* even a sphere in the plane where perspective is supposed to be 312* unity (z=0 for us). It swells slightly when perspective is 313* applied point by point. 314* This is a serious problem for us because although I 315* generate non-swelled spheres by applying the perspective to 316* all of them initially, I end up asking about individual 317* points on them when wanting the light source point of view. 318* You might argue that I could simply calculate the amount 319* of swelling that has to be accounted for (by drawing tangents 320* from the eye and seeing where they intersect the constant-z 321* plane that passes through the sphere centre or something), 322* but the "swelling" is complicated in the sense that it is 323* in effect a "stretching" of the surface closest to the eye 324* and a "shrinking" in back. 325* Even if I could see how to compensate for all this, I 326* don't think it would be worth it. 327* It's probably not even worth trying to implement a better 328* approximate solution by changing the light source position 329* slightly for each object in the picture. The ambiguity 330* as to whether the obscuring or the obscured object should have 331* the modified light source position applied would make it 332* difficult to assign objects to shadow tiles. Trying to 333* implement full "antiperspective" for the light source would 334* just shift the problems of perspective (swelling spheres, 335* etc.) to a different locale without solving them. 336*----------------------------------------------------------------------------- 337* 338* Overkill: 339 IMPLICIT REAL (A-H, O-Z) 340* 341 INCLUDE 'VERSION.incl' 342* 343* I/O units for control input, info output, label processing 344 INTEGER INPUT, INPUT0, NOISE, LUNIT 345 PARAMETER (INPUT0=5, NOISE=0, LUNIT=4) 346* 347* Descriptor codes for the various object types 348 INTEGER TRIANG, SPHERE, INTERNAL, CYLIND, CYLFLAT 349 INTEGER PLANE, QUADRIC, MXTYPE, FONT, GLOWLIGHT 350 INTEGER NORMS, VERTEXRGB, VERTRANSP 351 PARAMETER (TRIANG = 1) 352 PARAMETER (SPHERE = 2) 353 PARAMETER (CYLIND = 3) 354 PARAMETER (INTERNAL = 4) 355 PARAMETER (CYLFLAT = 5) 356 PARAMETER (PLANE = 6) 357 PARAMETER (NORMS = 7) 358 PARAMETER (MATERIAL = 8) 359 PARAMETER (MATEND = 9) 360 PARAMETER (FONT = 10, LABEL = 11) 361 PARAMETER (GLOWLIGHT= 13) 362 PARAMETER (QUADRIC = 14) 363 PARAMETER (ISOLATE1 = 15) 364 PARAMETER (GPROP = 16) 365 PARAMETER (VERTEXRGB= 17) 366 PARAMETER (VERTRANSP= 18) 367 PARAMETER (ISOLATE2 = 19) 368 PARAMETER (MXTYPE = 19) 369* 370* Bit definitions for FLAG array 371 INTEGER FLAT, RIBBON, SURFACE, PROPS 372 PARAMETER (FLAT=2, RIBBON=4, SURFACE=8, PROPS=16) 373 INTEGER TRANSP, HIDDEN, INSIDE, MOPT1 374 PARAMETER (TRANSP=32, HIDDEN=64, INSIDE=128,MOPT1=256) 375 INTEGER VCOLS, CLIPPED, VTRANS, BOUNDED 376 PARAMETER (VCOLS=512, CLIPPED=1024, VTRANS=2048, BOUNDED=4096) 377* 378* Bit definitions for OTMODE passed to local(1,...) 379 INTEGER ALPHACHANNEL 380 PARAMETER (ALPHACHANNEL=32) 381* 382* $$$$$$$$$$$$$ ARRAY SIZE LIMITS $$$$$$$$$$$$$$ 383* 384 INCLUDE 'parameters.incl' 385 INTEGER OUTSIZ 386 PARAMETER (OUTSIZ = MAXNTX*MAXNPX*MAXNPY) 387* 388* 389* Command line options (Aug 1999) NB: nax,nay,quality MUST be integer*2 390 COMMON /OPTIONS/ FONTSCALE, GAMMA, ZOOM, NSCHEME, SHADOWFLAG, XBG, 391 & NAX, NAY, OTMODE, QUALITY, INVERT, LFLAG 392 REAL FONTSCALE, GAMMA, ZOOM 393 INTEGER NSCHEME, SHADOWFLAG, XBG 394 INTEGER*4 NAX, NAY, OTMODE, QUALITY 395 LOGICAL*2 INVERT, LFLAG 396* 397* Title for run 398 CHARACTER*132 TITLE 399* 400* Number of tiles, pixels per tile 401 COMMON /RASTER/ NTX,NTY,NPX,NPY 402 INTEGER NTX,NTY,NPX,NPY 403* 404* Pixels per tile after anti-aliasing, output buffer line length 405 INTEGER NOX, NOY, NX 406* 407* Actual image size in pixels (may include partial tiling at the edges) 408* (MUST BE INTEGER*2 for call to local()!!!) 409C INTEGER*2 NAX, NAY 410* 411* One lonely tile 412 REAL, ALLOCATABLE, DIMENSION(:,:,:) :: TILE 413* 414* With an alpha blend channel 415 REAL, ALLOCATABLE, DIMENSION(:,:) :: ACHAN 416* 417* Pixel averaging scheme 418 INTEGER SCHEME 419* 420* Background colour 421 REAL BKGND(3) 422 INTEGER IBKGND(3) 423* 424* "Shadow mode?" 425 LOGICAL SHADOW 426 INTEGER NSXMAX,NSYMAX 427* 428* Phong power 429 INTEGER IPHONG 430* 431* Straight-on (secondary) light source contribution 432 REAL STRAIT 433* 434* Ambient light contribution 435 REAL AMBIEN 436* 437* Specular reflection component 438 REAL SPECLR 439* 440* Primary light source position 441 REAL SOURCE(3) 442* 443* Input transformation 444 COMMON /MATRICES/ XCENT, YCENT, SCALE, EYEPOS, SXCENT, SYCENT, 445 & TMAT, TINV, TINVT, SROT, SRTINV, SRTINVT 446 & ,RAFTER, TAFTER 447 REAL XCENT, YCENT, SCALE, SXCENT, SYCENT 448* Transformation matrix, inverse of transpose, and transposed inverse 449 REAL TMAT(4,4), TINV(4,4), TINVT(4,4) 450* Shortest rotation from light source to +z axis 451 REAL SROT(4,4), SRTINV(4,4), SRTINVT(4,4) 452* Post-hoc transformation on top of original TMAT 453 REAL RAFTER(4,4), TAFTER(3) 454 EXTERNAL DET 455 REAL DET 456* 457* Distance (in +z) of viewing eye 458 REAL EYEPOS 459* 460* Input mode 461 INTEGER INMODE 462* 463* Buffer one line of input for decoding 464 CHARACTER*132 LINE 465* 466* Input format(s) 467 CHARACTER*80 INFMTS(MXTYPE),INFMT 468* 469* Free-format input flag 470 LOGICAL INFLGS(MXTYPE),INFLG 471* 472* Allow very long names for file indirection 473 CHARACTER*132 FULLNAME 474* 475* Stuff for shading 476 REAL NL(3),NORMAL(3),LDOTN 477 REAL RGBCUR(3),RGBSHD(3),RGBFUL(3) 478 REAL SPECOL(3) 479* 480* FOG parameters 481* (fogtype -1 = none, 0 = linear depthcuing, 1 = exponential model) 482* (fogfront 0 = front object, else fraction of front clipping plane) 483* (fogback 0 = back object, else fraction of back clipping plane) 484* 485 COMMON /FOGCOM/ FOGTYPE,FOGFRONT,FOGBACK,FOGDEN,FOGLIM,FOGRGB 486 INTEGER FOGTYPE 487 REAL FOGFRONT, FOGBACK, FOGDEN, FOGLIM(2), FOGRGB(3) 488* 489* The s & m guys are for the shadow box in the following 490* 491* Object list, consists of pointers (less 1) into detail, sdtail 492 INTEGER, ALLOCATABLE, DIMENSION(:) :: LIST, MIST 493* 494* Object types and flags, parallel to list 495 INTEGER, ALLOCATABLE, DIMENSION(:) :: TYPE 496 INTEGER*4, ALLOCATABLE, DIMENSION(:) :: FLAG 497* 498* Keep a separate list of special materials 499* and remember any special props of current material on input 500CDEBUG MPARITY gets its own array because it's used in a per-pixel loop 501CDEBUG (using DETAIL(LIST(MLIST(MAT))+18) cost 5% in execution time) 502 INTEGER, ALLOCATABLE, DIMENSION(:) :: MLIST, MPARITY 503 LOGICAL MATCOL, BACKFACE 504 LOGICAL CLIPPING, MAYCLIP, JUSTCLIPPED 505 REAL RGBMAT(3) 506* 507* Object details, shadow object details 508 REAL, ALLOCATABLE, DIMENSION(:) :: DETAIL, SDTAIL 509* 510* Input buffer for details 511 REAL BUF(100) 512* 513* Number of objects in each tile's short list (m... are for shadows) 514C Moved from COMMON BLOCK LISTS to MODULE LISTS to allow dynamic 515C allocation - FZ 516C COMMON /LISTS/ KOUNT, MOUNT, TTRANS, ISTRANS 517C INTEGER KOUNT(MAXNTX,MAXNTY), MOUNT(NSX,NSY) 518C INTEGER TTRANS(MAXNTX,MAXNTY), ISTRANS 519* 520* Pointer to where each tile's objects start 521 INTEGER, ALLOCATABLE, DIMENSION(:,:) :: KSTART, MSTART 522* 523* Pointer to where each tile's objects end 524 INTEGER, ALLOCATABLE, DIMENSION(:,:) :: KSTOP, MSTOP 525* 526* Short list heap 527 INTEGER, ALLOCATABLE, DIMENSION(:) :: KSHORT, MSHORT 528* 529* Temporary for sorting 530 REAL, ALLOCATABLE, DIMENSION(:) :: ZTEMP 531* 532* Where the permutation representing the sort is stored 533 INTEGER, ALLOCATABLE, DIMENSION(:) :: ZINDEX 534* 535* The number of "details" each object type is supposed to have 536* : input, object, shadow 537 INTEGER IDET(MXTYPE), KDET(MXTYPE), SDET(MXTYPE) 538* 539* Support for cylinders 540 EXTERNAL CYLTEST 541 LOGICAL CYLTEST, ISCYL 542* 543* Support for quadrics 544 REAL QNORM(3) 545 EXTERNAL QINP, QTEST 546 LOGICAL QINP, QTEST, ISQUAD 547* 548* Support for transparency 549 COMMON /TRANS/ NTRANSP, INDEPTH, INDTOP, TRANOVFL, ZTOP, ZHIGH, 550 & INDLIST(MAXTRANSP), ZLIST(MAXTRANSP), 551 & NORMLIST(3,MAXTRANSP) 552 INTEGER NTRANSP, INDEPTH, INDTOP, INDLIST, TRANOVFL 553 REAL ZTOP, ZHIGH, ZLIST, NORMLIST 554 REAL SBLEND, RGBLND(3) 555* 556* Support for a "glow" light source 557 REAL GLOWSRC(3), GLOWCOL(3), GDIST(3), GLOWRAD, GLOW, GLOWMAX 558 INTEGER GOPT, GPHONG 559 INTEGER, ALLOCATABLE, DIMENSION(:) :: GLOWLIST 560 INTEGER NGLOWS 561* 562* Support for decompression on the fly 563 EXTERNAL ungz 564 INTEGER ungz 565* 566* Support for BOUNDING_PLANE internal object type 567 REAL BPLANE(3), BPNORM(3), BPRGB(3) 568 REAL xn, yn, zn, xnb, ynb, znb 569 INTEGER BPTYPE, NBOUNDS, BPIND 570 LOGICAL ORTEPLIKE 571 INTEGER ORTEP 572 PARAMETER (ORTEP=1) 573 EXTERNAL INBOUNDS 574 LOGICAL INBOUNDS 575 REAL TEMPNORM(3) 576* 577* Output buffer 578 INTEGER*2, ALLOCATABLE, DIMENSION(:,:) :: OUTBUF 579* 580* Copy of NOISE for ASSERT to see 581 INTEGER ASSOUT 582 LOGICAL VERBOSE 583 COMMON /ASSCOM/ ASSOUT, VERBOSE 584 SAVE /ASSCOM/ 585* 586* For label processing 587 COMMON /LABELS/ LABOUT 588 INTEGER LABOUT 589* 590* Gamma correction 591 INTEGER GAMMA_MAP(256) 592 PARAMETER (MAXRGB=255.0) 593 LOGICAL GAMMACORRECTION 594* 595* Keep track of actual coordinate limits 596 COMMON /NICETIES/ TRULIM, ZLIM, FRONTCLIP, BACKCLIP 597 & , ISOLATION 598 REAL TRULIM(3,2), ZLIM(2), FRONTCLIP, BACKCLIP 599 INTEGER ISOLATION 600* 601* Array of sizes to try allocating for expanded dynamic storage 602 INTEGER NEEDMEM, TRY1(3), TRY2(3) 603 INTEGER, ALLOCATABLE, DIMENSION(:) :: TMP1D 604 INTEGER*4, ALLOCATABLE, DIMENSION(:) :: TMP1DI4 605 REAL, ALLOCATABLE, DIMENSION(:) :: TMP1DR 606 INTEGER, ALLOCATABLE, DIMENSION(:,:) :: TMP2D 607 INTEGER*2, ALLOCATABLE, DIMENSION(:,:) :: TMP2DI2 608 REAL, ALLOCATABLE, DIMENSION(:,:) :: TMP2DR 609 REAL, ALLOCATABLE, DIMENSION(:,:,:) :: TMP3DR 610 611 LOGICAL TEST_ALLOC 612 TEST_ALLOC = .TRUE. 613* 614* Allocate initial space for dynamically allocatable arrays 615 ALLOCATE( TILE(3,MAXNPX,MAXNPY) ) 616 ALLOCATE( ACHAN(MAXNPX,MAXNPY) ) 617 ALLOCATE( ZTEMP(MAXOBJ) ) 618 ALLOCATE( ZINDEX(MAXOBJ) ) 619 ALLOCATE( LIST(MAXOBJ), MIST(MAXOBJ) ) 620 ALLOCATE( TYPE(MAXOBJ) ) 621 ALLOCATE( FLAG(MAXOBJ) ) 622 ALLOCATE( MLIST(MAXMAT), MPARITY(MAXMAT) ) 623 ALLOCATE( DETAIL(MAXDET), SDTAIL(MAXSDT) ) 624 ALLOCATE( KOUNT(MAXNTX,MAXNTY), MOUNT(NSX,NSY) ) 625 ALLOCATE( TTRANS(MAXNTX,MAXNTY) ) 626 ALLOCATE( KSTART(MAXNTX,MAXNTY), MSTART(NSX,NSY) ) 627 ALLOCATE( KSTOP(MAXNTX,MAXNTY), MSTOP(NSX,NSY) ) 628 ALLOCATE( KSHORT(MAXSHR), MSHORT(MAXSSL) ) 629 ALLOCATE( GLOWLIST(MAXGLOWS) ) 630 ALLOCATE( OUTBUF(OUTSIZ,4) ) 631 632 TRULIM (1,1) = HUGE 633 TRULIM (2,1) = HUGE 634 TRULIM (3,1) = HUGE 635 TRULIM (1,2) = -HUGE 636 TRULIM (2,2) = -HUGE 637 TRULIM (3,2) = -HUGE 638 ZLIM(1) = HUGE 639 ZLIM(2) = -HUGE 640* 641 IDET(TRIANG) = 12 642 IDET(SPHERE) = 7 643 IDET(CYLFLAT) = 11 644 IDET(CYLIND) = 11 645 IDET(PLANE) = 12 646 IDET(NORMS ) = 9 647 IDET(MATERIAL) = 10 648 IDET(GLOWLIGHT)= 10 649 IDET(QUADRIC) = 17 650 IDET(VERTEXRGB)= 9 651 IDET(VERTRANSP)= 3 652 IDET(INTERNAL) = 20 653* 654 KDET(TRIANG) = 16 655 KDET(SPHERE) = 7 656 KDET(CYLIND) = 11 657 KDET(PLANE) = 7 658 KDET(NORMS ) = 9 659 KDET(MATERIAL) = 18 660 KDET(GLOWLIGHT)= 10 661 KDET(QUADRIC) = 17 662 KDET(VERTEXRGB)= 9 663 KDET(VERTRANSP)= 3 664 KDET(INTERNAL) = 20 665* 666 SDET(TRIANG) = 13 667 SDET(SPHERE) = 4 668 SDET(CYLIND) = 8 669 SDET(QUADRIC) = 14 670 SDET(INTERNAL) = 20 671* These object types really have no shadow details, 672* but indexing seems to require a nonzero value 673 SDET(PLANE) = 1 674 SDET(NORMS ) = 1 675 SDET(MATERIAL) = 1 676 SDET(GLOWLIGHT)= 1 677 SDET(VERTEXRGB)= 1 678 SDET(VERTRANSP)= 1 679* 680* Feb 2008 - more initializations for current gfortran 681 NSXMAX = 0 682 NSYMAX = 0 683 CLROPT = 0 684 SCHEME = 0 685 TRNSPOPT = 0 686 INFLG = .TRUE. 687 JUSTCLIPPED = .FALSE. 688 IXHI = 0 689 IXLO = 0 690 IYHI = 0 691 IYLO = 0 692 ITPASS = 0 693* 694* Copy the info (also error reporting) unit number to common 695 ASSOUT = NOISE 696 WRITE (NOISE,*) ' ' 697* 698* Initialize to level 0 of file indirection 699 INPUT = INPUT0 700* 701* Initialize unit number for label processing 702 LABOUT = LUNIT 703* 704* Initialize to no special material properties 705 MSTATE = 0 706 MATCOL = .FALSE. 707 ISOLATION = 0 708 CLIPPING = .FALSE. 709 NBOUNDS = 0 710 ORTEPLIKE = .FALSE. 711 CLRITY = 0.0 712 GLOWMAX = 0.0 713* 714* Initialize to no perspective. EYEPOS > 0 will add perspective 715 PFAC = 1.0 716 PFAC1 = 1.0 717 PFAC2 = 1.0 718 PFAC3 = 1.0 719* 720* Initialize global properties 721 FOGTYPE = -1 722 723 RGBLND(1) = 0 724 RGBLND(2) = 0 725 RGBLND(3) = 0 726* 727* EAM Aug 1999 - break out command line parsing into new routine 728 call parse 729* 730* Get title 731 100 CONTINUE 732 DO I=1,132 733 TITLE(I:I) = ' ' 734 ENDDO 735 736 READ (INPUT,'(A)',END=104,ERR=104) TITLE 737 IF (TITLE(1:1) .EQ. '#') GOTO 100 738 IF (TITLE(1:1) .EQ. '@') THEN 739 J = 1 740 K = 132 741 DO I=132,2,-1 742 IF (TITLE(I:I).EQ.' ') TITLE(I:I) = ' ' 743 IF (TITLE(I:I).NE.' ') J = I 744 IF (TITLE(I:I).EQ.'#') K = I-1 745 IF (TITLE(I:I).EQ.'!') K = I-1 746 ENDDO 747 IF (J.EQ.1) GOTO 101 748 DO WHILE (TITLE(K:K).EQ.' ') 749 K = K -1 750 ENDDO 751 OPEN (UNIT=INPUT+1,ERR=101,STATUS='OLD',FILE=TITLE(J:K)) 752 WRITE (NOISE,'(A,A)') ' + Opening input file ',TITLE(J:K) 753 INPUT = INPUT + 1 754 READ (INPUT,'(A)',ERR=101) TITLE 755 IF (TITLE(1:1) .EQ. '#') GOTO 100 756 ENDIF 757 GOTO 102 758 101 WRITE (NOISE,'(A,A)')' >> Cannot open or read file ',TITLE(2:K) 759 CALL EXIT(-1) 760 102 CONTINUE 761 K = 132 762 DO WHILE (TITLE(K:K).EQ.' ') 763 K = K -1 764 ENDDO 765 WRITE (NOISE,103) TITLE(1:K) 766 103 FORMAT('title="',A,'"') 767* 768* Get number of tiles 769 READ (INPUT,*,ERR=104,END=104) NTX,NTY 770 CALL ASSERT (NTX.GT.0, 'ntx.le.0') 771 CALL ASSERT (NTY.GT.0, 'nty.le.0') 772 GOTO 105 773104 CALL ASSERT(.FALSE., 774 & '>>> This doesnt look like a Raster3D input file! <<<') 775105 CONTINUE 776* 777* Get number of pixels per tile - 0 means autotile from values in NTX, NTY 778 READ (INPUT,*,ERR=104) NPX,NPY 779 if (npx.eq.0 .and. nax.le.0) nax = ntx 780 if (npy.eq.0 .and. nay.le.0) nay = nty 781* 782* Get pixel averaging scheme 783 READ (INPUT,*,ERR=104) SCHEME 784 if (nscheme.ge.0) scheme = nscheme 785 CALL ASSERT (SCHEME.GE.0 .AND. SCHEME.LE.4, 'bad scheme') 786* 787* Set up tiling and anti-aliasing. 788* If NAX, NAY are set, then use them to autotile 789 IF (SCHEME.LE.1) THEN 790 call autotile( nax, nay, 2 ) 791 NOX = NPX 792 NOY = NPY 793 if (nax.lt.0) nax = npx * ntx 794 if (nay.lt.0) nay = npy * nty 795 ELSEIF (SCHEME.EQ.2) THEN 796 if (nax.gt.0) nax = nax * 2 797 if (nay.gt.0) nay = nay * 2 798 if (nax.le.0) nax = NPX * NTX 799 if (nay.le.0) nay = NPY * NTY 800 call autotile( nax, nay, 2 ) 801 nax = nax / 2 802 nay = nay / 2 803 NOX = NPX/2 804 NOY = NPY/2 805 ELSEIF (SCHEME.EQ.3 .and. nscheme.ne.-4 .and. nax.le.0 ) THEN 806C Old style scheme 3 with exact tiling specified 807 nax = NPX * NTX 808 nay = NPY * NTY 809 call autotile( nax, nay, 3 ) 810 nax = (nax * 2 + 2) / 3 811 nay = (nay * 2 + 2) / 3 812 NOX = (NPX * 2) / 3 813 NOY = (NPY * 2) / 3 814 ELSE 815C Either scheme 4 or -size and anti-aliasing selected on command line 816 if (nax.gt.0) nax = (nax + (nax+1)/2) 817 if (nay.gt.0) nay = (nay + (nay+1)/2) 818 if (nax.le.0) nax = NPX*NTX + (NPX*NTX+1)/2 819 if (nay.le.0) nay = NPY*NTY + (NPY*NTY+1)/2 820 call autotile( nax, nay, 3 ) 821 nax = (nax * 2 + 2) / 3 822 nay = (nay * 2 + 2) / 3 823 NOX = (NPX * 2) / 3 824 NOY = (NPY * 2) / 3 825 SCHEME = 3 826 ENDIF 827* 828 call assert (nax.gt.0, 'nax <= 0') 829 call assert (nay.gt.0, 'nay <= 0') 830* 831 CALL ASSERT (NTX.GT.0.,'Tiling failure - ntx = 0') 832 CALL ASSERT (NTY.GT.0.,'Tiling failure - nty = 0') 833 CALL ASSERT (NPX.GT.0.,'Tiling failure - npx = 0') 834 CALL ASSERT (NPY.GT.0.,'Tiling failure - npy = 0') 835* 836* Expand arrays KOUNT, TTRANS, KSTART and KSTOP if needed, eg NTX > MAXNTX 837 if (NTX.GT.SIZE(KSTOP,1) .OR. NTY.GT.SIZE(KSTOP,2)) THEN 838* 839* Double the old allocation if that's enough, or take 150% of the new 840* If that fails, try 150% of the old or 120% of new, or just the new 841 CALL GET_TRY(SIZE(KSTOP, 1), NTX, TRY1, 2) 842 CALL GET_TRY(SIZE(KSTOP, 2), NTY, TRY2, 2) 843* 844* Try allocating memory to the arrays 845 do 900 ITRY = 1,3 846* 847* Test to see if requested allocation is valid 848 if (TRY1(ITRY) .LE. 0 .OR. TRY2(ITRY) .LE. 0 .OR. 849 & TRY1(ITRY)*TRY2(ITRY) .LE. 0 .OR. 850 & MAXMEM / TRY1(ITRY) .LT. TRY2(ITRY)) GOTO 900 851 852 ALLOCATE( TMP2D(TRY1(ITRY),TRY2(ITRY)), STAT=IERR) 853 if (ierr .NE. 0) GOTO 900 854C TMP2D = 0 855 TMP2D = KOUNT 856 CALL MOVE_ALLOC(from=TMP2D, to=KOUNT) 857 858 ALLOCATE( TMP2D(TRY1(ITRY),TRY2(ITRY)), STAT=IERR) 859 if (ierr .NE. 0) GOTO 900 860C TMP2D = 0 861 TMP2D = TTRANS 862 CALL MOVE_ALLOC(from=TMP2D, to=TTRANS) 863 864 ALLOCATE( TMP2D(TRY1(ITRY),TRY2(ITRY)),stat=ierr) 865 if (ierr .NE. 0) GOTO 900 866C TMP2D = 0 867 TMP2D = KSTART 868 CALL MOVE_ALLOC(from=TMP2D, to=KSTART) 869 870 ALLOCATE( TMP2D(TRY1(ITRY),TRY2(ITRY)),stat=ierr) 871 if (ierr .NE. 0) GOTO 900 872C TMP2D = 0 873 TMP2D = KSTOP 874 CALL MOVE_ALLOC(from=TMP2D, to=KSTOP) 875 if(TEST_ALLOC)write(NOISE,*)"Expand MAXNTX x Y to ", 876 & try1(ITRY)," x ",try2(ITRY) 877 GOTO 902 878900 CONTINUE 879 ENDIF 880* 881* These should only fail if the above dynamic allocation failed 882902 CALL ASSERT (NTX.LE.SIZE(KSTOP, 1),'Tiling failure - ntx>maxntx') 883 CALL ASSERT (NTY.LE.SIZE(KSTOP, 2),'Tiling failure - nty>maxnty') 884* 885* Expand arrays TILE, ACHAN next if needed, i.e. NPX,NPY > MAXNPX or Y 886 if (NPX.GT.SIZE(TILE,2).OR. NPY.GT.SIZE(TILE,3)) THEN 887 CALL GET_TRY(SIZE(TILE,2), NPX, TRY1, 2) 888 CALL GET_TRY(SIZE(TILE,3), NPY, TRY2, 2) 889 do 905 ITRY = 1,3 890* Test to see if requested allocation is valid 891 if (TRY1(ITRY) .LE. 0 .OR. TRY2(ITRY) .LE. 0 .OR. 892 & TRY1(ITRY)*TRY2(ITRY) .LE. 0 .OR. 893 & MAXMEM / TRY1(ITRY) .LT. TRY2(ITRY)) GOTO 905 894 895 ALLOCATE( TMP2DR(TRY1(ITRY),TRY2(ITRY)),stat=ierr) 896 if (ierr .NE. 0) GOTO 905 897 TMD2DR = 0. 898 TMP2DR = ACHAN 899 CALL MOVE_ALLOC(from=TMP2DR, to=ACHAN) 900 901 ALLOCATE( TMP3DR(3,TRY1(ITRY),TRY2(ITRY)),stat=ierr) 902 if (ierr .NE. 0) GOTO 905 903 TMD3DR = 0. 904 TMP3DR = TILE 905 CALL MOVE_ALLOC(from=TMP3DR, to=TILE) 906 if(TEST_ALLOC)write(NOISE,*)"Expand MAXNPX,Y to ", 907 & try1(ITRY)," x ",try2(ITRY) 908 GOTO 907 909905 CONTINUE 910 ENDIF 911* 912* These should only fail if the above dynamic allocation failed 913907 CALL ASSERT (NPX.LE.SIZE(TILE,2),'Tiling failure - npx>maxnpx') 914 CALL ASSERT (NPY.LE.SIZE(TILE,3),'Tiling failure - npy>maxnpy') 915* 916 IF (VERBOSE) THEN 917 WRITE (NOISE,*) 'ntx=',NTX,' nty=',NTY 918 WRITE (NOISE,*) 'npx=',NPX,' npy=',NPY 919 WRITE (NOISE,*) 'scheme=',SCHEME 920 WRITE (NOISE,*) 'nox=',NOX,' noy=',NOY 921 END IF 922 if (nax.lt.0) nax = nox*ntx 923 if (nay.lt.0) nay = noy*nty 924 NX = nox*ntx 925 LINOUT = 0 926 WRITE (NOISE,1105) 'Rendered raster size =',NPX*NTX,NPY*NTY 927 WRITE (NOISE,1105) ' Output raster size =',NAX,NAY 9281105 FORMAT(A,I7,' x',I7) 929C 930* 931* Expand array OUTBUF if needed, i.e. NTX * NOX * NOY > OUTBUF size 932* NOTE: this may avoidably fail if the above array expansions took 933* up more room than it needed, and didn't leave enough for OUTBUF. 934 NEEDMEM = NOY*NOX*NTX 935 if ( NEEDMEM .GT. SIZE(OUTBUF,1) ) THEN 936 CALL GET_TRY(SIZE(OUTBUF,1), NEEDMEM, TRY1, 1) 937 do 910 ITRY = 1,3 938* Test to see if requested allocation is valid 939 if (TRY1(ITRY).LE.0 .OR. TRY1(ITRY).GT.MAXMEM) GOTO 910 940 941 ALLOCATE( TMP2DI2(TRY1(ITRY),4), stat=ierr ) 942 if (ierr .NE. 0) GOTO 910 943 TMP2DI2 = 0 944 TMP2DI2 = OUTBUF 945 CALL MOVE_ALLOC(from=TMP2DI2, to=OUTBUF) 946 if(TEST_ALLOC)write(NOISE,*)"Expand OUTBUF to ",try1(ITRY) 947 GOTO 912 948910 CONTINUE 949 ENDIF 950912 CALL ASSERT (SIZE(OUTBUF,1).GE.NEEDMEM, 951 & 'image too large for output buffer') 952C 953C Header records and picture title 954 IF (SCHEME.EQ.0) OTMODE = ior(OTMODE,ALPHACHANNEL) 955* 956* Some derived parameters 957 XCENT = NTX*NPX/2. 958 YCENT = NTY*NPY/2. 959 SXCENT = NSX*NPX/2. 960 SYCENT = NSY*NPY/2. 961 SCALE = 2.*MIN(XCENT,YCENT) 962* This was always true; now it's explicit 963 BACKCLIP = -(SCALE+1.0) 964 FRONTCLIP = HUGE 965* Copy scheme to common, where r3dtogd can see it 966 NSCHEME = SCHEME 967* 968* Get background colour 969 READ (INPUT,*,ERR=104) BKGND 970 if (XBG.NE.0) then 971 BKGND(3) = iand(XBG,Z'00FF') 972 BKGND(2) = iand(XBG,Z'FF00')/256 973 BKGND(1) = iand(XBG,Z'FF0000')/65536 974 BKGND(3) = BKGND(3) / 255. 975 BKGND(2) = BKGND(2) / 255. 976 BKGND(1) = BKGND(1) / 255. 977 BKGND(3) = BKGND(3)**2 978 BKGND(2) = BKGND(2)**2 979 BKGND(1) = BKGND(1)**2 980 endif 981 IF (VERBOSE) THEN 982 WRITE (NOISE,1106) 'bkgnd=',BKGND 983 END IF 9841106 FORMAT(A,4F10.4,(/,4F10.4)) 985 CALL ASSERT (BKGND(1).GE.0., 'bkgnd(1) < 0') 986 CALL ASSERT (BKGND(2).GE.0., 'bkgnd(2) < 0') 987 CALL ASSERT (BKGND(3).GE.0., 'bkgnd(3) < 0') 988 CALL ASSERT (BKGND(1).LE.1., 'bkgnd(1) > 1') 989 CALL ASSERT (BKGND(2).LE.1., 'bkgnd(2) > 1') 990 CALL ASSERT (BKGND(3).LE.1., 'bkgnd(3) > 1') 991 IF (GAMMA.LT.0.99 .OR. GAMMA.GT.1.01) THEN 992 IBKGND(1) = SQRT(BKGND(1)) ** (1.0/GAMMA) * MAXRGB + 0.5 993 IBKGND(2) = SQRT(BKGND(2)) ** (1.0/GAMMA) * MAXRGB + 0.5 994 IBKGND(3) = SQRT(BKGND(3)) ** (1.0/GAMMA) * MAXRGB + 0.5 995 ELSE 996 IBKGND(1) = 255. * SQRT(BKGND(1)) + .5 997 IBKGND(2) = 255. * SQRT(BKGND(2)) + .5 998 IBKGND(3) = 255. * SQRT(BKGND(3)) + .5 999 ENDIF 1000* 1001* 1002* Get "shadows" flag 1003 READ (INPUT,*,ERR=104) SHADOW 1004 IF (SHADOWFLAG.EQ.0) SHADOW = .FALSE. 1005 IF (SHADOWFLAG.EQ.1) SHADOW = .TRUE. 1006 IF (VERBOSE) THEN 1007 WRITE (NOISE,*) 'shadow=',SHADOW 1008 END IF 1009* 1010* Get Phong power 1011 READ (INPUT,*,ERR=104) PHONG 1012 IPHONG = PHONG 1013 CALL ASSERT (IPHONG.GE.0, 'iphong < 0') 1014* A derived constant for numerical purposes in applying the 1015* Phong power in the shading algorithm. 1016* The idea is that any specular contribution less than 1017* 1E-9 (hence the 9 in 9./IPHONG) is insignificant: 1018 IF (IPHONG .NE. 0) PHOBND = 0.1**(9./IPHONG) 1019 IF (IPHONG .EQ. 0) PHOBND = 0. 1020* 1021* Get contribution of straight-on (secondary) light source 1022 READ (INPUT,*,ERR=104) STRAIT 1023 CALL ASSERT (STRAIT.GE.0., 'strait < 0') 1024 CALL ASSERT (STRAIT.LE.1., 'strait > 1') 1025* 1026* Derive contribution of primary light source 1027 PRIMAR = 1. - STRAIT 1028* 1029* Get contribution of ambient light 1030 READ (INPUT,*,ERR=104) AMBIEN 1031 CALL ASSERT (AMBIEN.GE.0., 'ambien < 0') 1032 CALL ASSERT (AMBIEN.LE.1., 'ambien > 1') 1033* 1034* Get component of specular reflection 1035 READ (INPUT,*,ERR=104) SPECLR 1036 CALL ASSERT (SPECLR.GE.0., 'speclr < 0') 1037 CALL ASSERT (SPECLR.LE.1., 'speclr > 1') 1038* 1039 IF (VERBOSE) THEN 1040 WRITE (NOISE,1104) 'iphong=',float(IPHONG),'strait=',STRAIT, 1041 & 'ambien=',AMBIEN,'speclr=',SPECLR 10421104 FORMAT(2(4X,A,F10.4)) 1043 END IF 1044* 1045* Derive component of diffuse reflection 1046 CALL ASSERT (AMBIEN+SPECLR.LE.1., 'ambien+speclr > 1') 1047 DIFFUS = 1. - (AMBIEN+SPECLR) 1048* 1049* Get distance of viewing eye 1050 READ (INPUT,*,ERR=104) EYEPOS 1051 CALL ASSERT (EYEPOS.GE.0., 'eyepos.lt.0') 1052* 1053* Get position of primary light source 1054 READ (INPUT,*,ERR=104) SOURCE 1055 SMAG = SQRT(SOURCE(1)**2 + SOURCE(2)**2 + SOURCE(3)**2) 1056 SOURCE(1) = SOURCE(1) / SMAG 1057 SOURCE(2) = SOURCE(2) / SMAG 1058 SOURCE(3) = SOURCE(3) / SMAG 1059 IF (VERBOSE) THEN 1060 WRITE (NOISE,1106) 'eyepos=',EYEPOS 1061 WRITE (NOISE,1106) 'source=',SOURCE 1062 WRITE (NOISE,1106) 'normalized source=',SOURCE 1063 END IF 1064* 1065* Get input transformation 1066 DO I=1,4 1067 READ (INPUT,*,ERR=104) (TMAT(I,J),J=1,4) 1068 END DO 1069 IF (VERBOSE) THEN 1070 WRITE (NOISE,*) 'tmat (v'' = v * tmat):' 1071 DO I=1,4 1072 WRITE (NOISE,'(4f9.4)') (TMAT(I,J),J=1,4) 1073 END DO 1074 END IF 1075* 1076* Initialize output file 1077 IERR = LOCAL(5, IBKGND(1), IBKGND(2), IBKGND(3)) 1078 IERR = LOCAL(1, NAX, NAY, OTMODE, QUALITY) 1079 IERR = LOCAL(4, TITLE) 1080* 1081* Allow command line rescaling option 1082 IF (ZOOM.LT.0.) ZOOM = -ZOOM / 100. 1083 IF (ZOOM.GT.0.) TMAT(4,4) = TMAT(4,4) / ZOOM 1084 IF (VERBOSE .AND. ZOOM.NE.0.) THEN 1085 WRITE (NOISE,1106) 'zoom factor = ',ZOOM 1086 ENDIF 1087* 1088* EAM - The original output mode was "upside down" compared 1089* to what most graphics programs expect to see. It is messy 1090* to change the evaluation order everywhere so that pixels can be 1091* streamed to stdout, so instead I invert the Y axis in TMAT and SOURCE 1092* here. 1093* The actual decision whether or not to invert is done in local.c 1094* and returned as a bit in the status word returned by local(0,...) 1095 IF (INVERT) THEN 1096 DO I = 1,4 1097 TMAT(I,2) = -TMAT(I,2) 1098 ENDDO 1099 SOURCE(2) = -SOURCE(2) 1100 ENDIF 1101* 1102* By popular demand, add a post-hoc rotation/translation option 1103* that uses matrices of the form used by O and molscript 1104* Initialized here to identity matrix; set by GPROP options. 1105 DO I = 1,4 1106 DO J = 1,4 1107 RAFTER(I,J) = 0.0 1108 ENDDO 1109 RAFTER(I,I) = 1.0 1110 ENDDO 1111 TAFTER(1) = 0.0 1112 TAFTER(2) = 0.0 1113 TAFTER(3) = 0.0 1114* 1115* Compute the rotation matrix which takes the light 1116* source to the +z axis (i.e., to the viewpoint). 1117* first make p = source cross z (and normalize p) 1118 P1 = SOURCE(2) 1119 P2 = -SOURCE(1) 1120* p3 = 0 1121 PLEN = SQRT(P1**2 + P2**2) 1122 IF (PLEN .GT. 0.0) P1 = P1 / PLEN 1123 IF (PLEN .GT. 0.0) P2 = P2 / PLEN 1124* phi is the angle between source and z (shortest route) 1125 COSPHI = SOURCE(3) 1126 SINPHI = PLEN 1127 SROT(1,1) = P1**2 + (1.-P1**2)*COSPHI 1128 SROT(1,2) = P1*P2*(1.-COSPHI) 1129 SROT(1,3) = P2*SINPHI 1130 SROT(2,1) = SROT(1,2) 1131 SROT(2,2) = P2**2 + (1.-P2**2)*COSPHI 1132 SROT(2,3) = -P1*SINPHI 1133 SROT(3,1) = -SROT(1,3) 1134 SROT(3,2) = -SROT(2,3) 1135 SROT(3,3) = COSPHI 1136 SROT(1,4) = 0.0 1137 SROT(2,4) = 0.0 1138 SROT(3,4) = 0.0 1139 SROT(4,1) = 0.0 1140 SROT(4,2) = 0.0 1141 SROT(4,3) = 0.0 1142 SROT(4,4) = 1.0 1143* 1144* Quadrics will require the inverse matrix also (and its transpose) 1145* This is also a convenient place to check legality of TMAT 1146 CALL QSETUP 1147* 1148* Get input mode 1149 READ (INPUT,*,ERR=104) INMODE 1150C WRITE (NOISE,*) 'inmode=',INMODE 1151 CALL ASSERT (INMODE.GE.1,'bad inmode') 1152* 1153* Get input format(s) 1154 IF (INMODE.EQ.1.OR.INMODE.EQ.2) THEN 1155 READ (INPUT,'(A)',ERR=104) INFMT 1156C WRITE (NOISE,*) 'infmt=',INFMT 1157 II = 0 11582 CONTINUE 1159 IF (INFMT(1:1).EQ.' ') THEN 1160 INFMT(1:79) = INFMT(2:80) 1161 INFMT(80:80) = ' ' 1162 II = II + 1 1163 IF (II.LT.80) GO TO 2 1164 ENDIF 1165 IF (INFMT(1:1).EQ.'*') THEN 1166 INFLG = .TRUE. 1167 ELSE 1168 INFLG = .FALSE. 1169 ENDIF 1170 ELSEIF (INMODE.GE.3) THEN 1171C WRITE (NOISE,*) 'infmts:' 1172 DO 4 I=1,3 1173 READ (INPUT,'(A)',ERR=104) INFMTS(I) 1174C WRITE (NOISE,*) INFMTS(I) 1175 II = 0 11763 CONTINUE 1177 IF (INFMTS(I)(1:1).EQ.' ') THEN 1178 INFMTS(I)(1:79) = INFMTS(I)(2:80) 1179 INFMTS(I)(80:80) = ' ' 1180 II = II + 1 1181 IF (II.LT.80) GO TO 3 1182 ENDIF 1183 IF (INFMTS(I)(1:1).EQ.'*') THEN 1184 INFLGS(I) = .TRUE. 1185 ELSE 1186 INFLGS(I) = .FALSE. 1187 ENDIF 11884 CONTINUE 1189 INFLGS(PLANE) = INFLGS(TRIANG) 1190 INFMTS(PLANE) = INFMTS(TRIANG) 1191 INFLGS(NORMS) = INFLGS(TRIANG) 1192 INFMTS(NORMS) = INFMTS(TRIANG) 1193 INFLGS(VERTEXRGB) = INFLGS(TRIANG) 1194 INFMTS(VERTEXRGB) = INFMTS(TRIANG) 1195 INFLGS(VERTRANSP) = INFLGS(TRIANG) 1196 INFMTS(VERTRANSP) = INFMTS(TRIANG) 1197 INFLGS(MATERIAL) = .TRUE. 1198 INFLGS(GLOWLIGHT) = .TRUE. 1199 INFLGS(QUADRIC) = .TRUE. 1200 ELSE 1201 CALL ASSERT (.FALSE.,'bad inmode') 1202 ENDIF 1203c 1204c Done with header records 1205c Do we force-close the input file, so that we can borrow headers, 1206c or should we keep going as long as the file continues? 1207c The following 4 lines implement the former, so that the initial 1208c @file command essentially means 'use his header records for me too'. 1209c 1210c IF (INPUT.NE.INPUT0) THEN 1211c CLOSE(INPUT) 1212c INPUT = INPUT0 1213c ENDIF 1214c As of V2.5e, however, we keep reading. 1215c >>> This is a change! <<< 1216c 1217c End of header processing 1218* 1219* Give them a notice to stare at while the program cranks along 1220 WRITE (NOISE,'(1X)') 1221 WRITE (NOISE,*) '%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%', 1222 & '%%%%%%%%%%%%%%%%%%%%%%%%%%' 1223 WRITE (NOISE,*) '% Raster3D ',VERSION, 1224 & ' %' 1225 WRITE (NOISE,*) '% -------------------------', 1226 & '------------- %' 1227 WRITE (NOISE,*) '% If you publish figures generated by this ', 1228 & 'program please cite %' 1229 WRITE (NOISE,*) '% Merritt & Bacon (1997) ', 1230 & 'Meth. Enzymol. 277, 505-524.',' %' 1231 WRITE (NOISE,*) '% -------------------------', 1232 & '------------- %' 1233 WRITE (NOISE,*) '% Raster3D distribution site', 1234 & ' %' 1235 WRITE (NOISE,*) '% http://www.bmsc.washington.edu/', 1236 & 'raster3d/ %' 1237 WRITE (NOISE,*) '% comments & suggestions to: ', 1238 & ' merritt@u.washington.edu %' 1239 WRITE (NOISE,*) '%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%', 1240 & '%%%%%%%%%%%%%%%%%%%%%%%%%%' 1241* 1242* If label processing is selected on command line, initialize 1243* PostScript output file. This isn't needed for Version 3, which 1244* uses libgd rather than PostScript to handle labels 1245 PSSCALE = 2.0 * MIN( NTX*NOX/2., NTY*NOY/2. ) 1246 CALL LSETUP( PSSCALE, BKGND, TITLE ) 1247 WRITE (NOISE,'(1X)') 1248* 1249* Initialize gamma correction table 1250 GAMMACORRECTION = .FALSE. 1251 IF (GAMMA.LT.0.99 .OR. GAMMA.GT.1.01) GAMMACORRECTION = .TRUE. 1252 IF (GAMMACORRECTION) THEN 1253 DO I=0,MAXRGB 1254 G = I 1255 GAMMA_MAP(I+1) = (G/MAXRGB) ** (1.0/GAMMA) * MAXRGB + 0.5 1256 ENDDO 1257 ENDIF 1258* 1259* 1260* Initialize counters 1261 DO 5 J=1,NTY 1262 DO 5 I=1,NTX 1263 KOUNT(I,J) = 0 1264 TTRANS(I,J) = 0 12655 CONTINUE 1266 DO 6 J=1,NSY 1267 DO 6 I=1,NSX 1268 MOUNT(I,J) = 0 12696 CONTINUE 1270c 1271 DO 662 I = 1,SIZE(FLAG) 1272 FLAG(I) = 0 1273 662 CONTINUE 1274 nprops = 0 1275 npropm = 0 1276 ntransp = 0 1277 nsphere = 0 1278 ncylind = 0 1279 nplanes = 0 1280 nhidden = 0 1281 ninside = 0 1282 mstate = 0 1283 nlabels = 0 1284 nglows = 0 1285 nquads = 0 1286 nclip = 0 1287 nvtrans = 0 1288 nbplanes = 0 1289 tranovfl = 0 1290c 1291c Objects in, and count up objects that may impinge on each tile 1292 NDET = 0 1293 IF (SHADOW) MDET = 0 1294 N = 0 1295c 1296c Read in next object 12977 CONTINUE 1298 IF (INMODE.EQ.1.OR.INMODE.EQ.2) THEN 1299 INTYPE = INMODE 1300 GOTO 8 1301 ENDIF 1302C 1303C READ (INPUT,*,END=50) INTYPE 1304 READ (INPUT,'(A)',END=50,ERR=47) LINE 1305c Nov 1997 - allow # comments 1306 IF (LINE(1:1) .EQ. '#') THEN 1307 GOTO 7 1308c May 1996 - allow file indirection 1309 ELSE IF (LINE(1:1) .EQ. '@') THEN 1310 J = 1 1311 K = 132 1312 L = 132 1313 DO I=132,2,-1 1314 IF (LINE(I:I).NE.' ') J = I 1315 IF (LINE(I:I).EQ.'#') K = I-1 1316 IF (LINE(I:I).EQ.'!') K = I-1 1317 IF (LINE(I:I).EQ.CHAR(0)) K = I-1 1318 IF (LINE(I:I).EQ.' ') LINE(I:I) = ' ' 1319 ENDDO 1320 IF (J.EQ.1) GOTO 7 1321 L=0 1322 DO I=J,K 1323 IF (LINE(I:I).NE.' ') L = I 1324 ENDDO 1325 K = L 1326 IF (LINE(K:K).eq.'Z' .or. LINE(K-1:K).eq.'gz') THEN 1327c ungz will uncompress into a temporary file, which ought to 1328c be deleted later. Unfortunately, that's hard to do in g77 1329c since it doesn't support dispose='DELETE'. 1330 line(k+1:k+1) = char(0) 1331 if (0 .gt. ungz( line(j:k+1), fullname )) goto 73 1332 j = 1 1333 k = 132 1334 do i=132,2,-1 1335 if (fullname(i:i).eq.' ') k = i-1 1336 if (fullname(i:i).eq.char(0)) k = i-1 1337 enddo 1338 if (verbose) 1339 & write(noise,*) 'Creating temporary file: ',fullname(j:k) 1340 open (unit=input+1,err=73,status='OLD', 1341 & file=fullname(j:k)) 1342 fullname = line(2:132) 1343 goto 72 1344 ENDIF 1345 70 CONTINUE 1346 OPEN (UNIT=INPUT+1,ERR=71,STATUS='OLD', 1347 & FILE=LINE(J:K)) 1348 FULLNAME = LINE(J:K) 1349 GOTO 72 1350 71 CONTINUE 1351 IF (LINE(K-3:K).NE.'.r3d') THEN 1352 K = K + 4 1353 LINE(K-3:K) = '.r3d' 1354 GOTO 70 1355 ENDIF 1356 CALL LIBLOOKUP( LINE(J:K), FULLNAME ) 1357 OPEN (UNIT=INPUT+1,ERR=73,STATUS='OLD',FILE=FULLNAME) 1358 72 CONTINUE 1359 DO I=132,2,-1 1360 IF (FULLNAME(I:I).EQ.' ') J = I 1361 ENDDO 1362 WRITE (NOISE,'(A,A)') ' + Opening input file ',FULLNAME(1:J) 1363 INPUT = INPUT + 1 1364 CALL ASSERT(INPUT-INPUT0.LE.MAXLEV, 1365 & 'Too many levels of indirection') 1366 GOTO 7 1367 73 WRITE (NOISE,'(A,A)') ' >> Cannot open file ',LINE(J:K) 1368 GOTO 7 1369 ELSE 1370c April 2011 - ignore blank line 1371 READ (LINE,*,END=7,ERR=74) INTYPE 1372 GOTO 76 1373 74 WRITE (NOISE,'(A,A)') ' >> Unrecognized line: ',LINE 1374 GOTO 7 1375 76 CONTINUE 1376 ENDIF 1377 IF (INTYPE.EQ.0) GO TO 50 1378 IF (INTYPE .EQ. CYLFLAT) THEN 1379 INTYPE = CYLIND 1380 FLAG(N+1) = ior( FLAG(N+1), FLAT ) 1381 ELSEIF (INTYPE .EQ. MATEND) THEN 1382 MSTATE = 0 1383 CLRITY = 0 1384 CLROPT = 0 1385 TRNSPOPT = 0 1386 MATCOL = .FALSE. 1387 ISOLATION = 0 1388 CLIPPING = .FALSE. 1389 NBOUNDS = 0 1390 ORTEPLIKE = .FALSE. 1391 GOTO 7 1392 ELSEIF (INTYPE .EQ. FONT) THEN 1393 IF (LFLAG) THEN 1394 CALL LINP( INPUT, INTYPE, .FALSE., RGBMAT ) 1395 ELSE 1396 READ (INPUT,'(A)',END=50) LINE 1397 ENDIF 1398 GOTO 7 1399 ELSEIF (INTYPE .EQ. LABEL) THEN 1400 NLABELS = NLABELS + 1 1401 IF (LFLAG) THEN 1402 IF (MSTATE.EQ.MATERIAL .AND. MATCOL) THEN 1403 CALL LINP( INPUT, INTYPE, .TRUE., RGBMAT ) 1404 ELSE 1405 CALL LINP( INPUT, INTYPE, .FALSE., RGBMAT ) 1406 ENDIF 1407 ELSE 1408 READ (INPUT,'(A)',END=50) LINE 1409 READ (INPUT,'(A)',END=50) LINE 1410 ENDIF 1411 GOTO 7 1412 ELSEIF (INTYPE .EQ. ISOLATE1) THEN 1413 ISOLATION = 1 1414 GOTO 7 1415 ELSEIF (INTYPE .EQ. ISOLATE2) THEN 1416 ISOLATION = 2 1417 GOTO 7 1418c Global Properties 1419 ELSEIF (INTYPE .EQ. GPROP) THEN 1420 READ (INPUT,'(A)',END=50) LINE 1421 L = 1 1422 DO I = 132, 1, -1 1423 IF (LINE(I:I).NE.' '.AND.LINE(I:I).NE.' ') L = I 1424 ENDDO 1425 IF (LINE(L:L+2).EQ.'FOG') THEN 1426 READ (LINE(L+4:74),*,ERR=771,END=771) 1427 & FOGTYPE, FOGFRONT, FOGBACK, FOGDEN 1428 771 CONTINUE 1429 FOGRGB(1) = BKGND(1) 1430 FOGRGB(2) = BKGND(2) 1431 FOGRGB(3) = BKGND(3) 1432 CALL CHKRGB(FOGRGB(1),FOGRGB(2),FOGRGB(3),'invalid fog color') 1433 IF (FOGTYPE.NE.1) FOGTYPE = 0 1434 IF (FOGDEN.LE.0.0) FOGDEN = 0.5 1435 ELSE IF (LINE(L:L+8).EQ.'FRONTCLIP') THEN 1436 READ (LINE(L+10:132),*,ERR=75) FRONTCLIP 1437 FRONTCLIP = FRONTCLIP * SCALE / TMAT(4,4) 1438 ELSE IF (LINE(L:L+7).EQ.'BACKCLIP') THEN 1439 READ (LINE(L+9:132),*,ERR=75) BACKCLIP 1440 BACKCLIP = BACKCLIP * SCALE / TMAT(4,4) 1441 ELSE IF (LINE(L:L+7).EQ.'ROTATION') THEN 1442 READ (LINE(L+9:132),*,ERR=773,END=773) 1443 & ((RAFTER(I,J),J=1,3),I=1,3) 1444 GOTO 774 1445 773 READ (INPUT,*,ERR=75) ((RAFTER(I,J),J=1,3),I=1,3) 1446 774 WRITE (NOISE,775) ((RAFTER(I,J),J=1,3),I=1,3) 1447 775 FORMAT('Post-rotation matrix: ',3(/,3F10.4)) 1448 D = DET(RAFTER) 1449 IF (ABS(1.0-ABS(D)).GT.0.02) WRITE (NOISE,*) 1450 & '>>> Warning: Post-rotation matrix has determinant',D 1451 IF (INVERT) THEN 1452 RAFTER(1,2) = -RAFTER(1,2) 1453 RAFTER(2,1) = -RAFTER(2,1) 1454 RAFTER(2,3) = -RAFTER(2,3) 1455 RAFTER(3,2) = -RAFTER(3,2) 1456 ENDIF 1457 CALL QSETUP 1458 ELSE IF (LINE(L:L+10).EQ.'TRANSLATION') THEN 1459 READ (LINE(L+12:132),*,ERR=776,END=776) 1460 & (TAFTER(I),I=1,3) 1461 GOTO 777 1462 776 READ (INPUT,*,ERR=75) (TAFTER(I),I=1,3) 1463 777 WRITE (NOISE,778) (TAFTER(I),I=1,3) 1464 778 FORMAT('Post-translation: ',1(/,3F10.4)) 1465 IF (INVERT) TAFTER(2) = -TAFTER(2) 1466 ELSE IF (LINE(L:L+4).EQ.'DUMMY') THEN 1467 ELSE 1468 GOTO 75 1469 ENDIF 1470 GOTO 7 1471 1472 75 CONTINUE 1473 WRITE(NOISE,'(A,A)') 1474 & '>> Unrecognized or incomplete GPROP option ',LINE 1475 GOTO 7 1476* 1477 ENDIF 1478* 1479COLD CALL ASSERT (INTYPE.GE.1.AND.INTYPE.LE.MXTYPE,'bad object') 1480 IF (INTYPE.LT.1 .OR. INTYPE.GT.MXTYPE) GOTO 47 1481 CALL ASSERT (INTYPE.NE.INTERNAL,'object type 4 not available') 1482c 1483c Read in object details, now we know what kind it is. 1484c Allow an all-zeroes record to terminate input for the 1485c benefit of those of us who might inadvertently supply 1486c a series of blank records after our real input as a 1487c side-effect of tape blocking or sloppiness or ... 14888 CONTINUE 1489 IF (INMODE.GE.3) THEN 1490 INFMT = INFMTS(INTYPE) 1491 INFLG = INFLGS(INTYPE) 1492 ENDIF 1493 IF (INFLG) THEN 1494 READ (INPUT,*,END=50) (BUF(I),I=1,IDET(INTYPE)) 1495 ELSE 1496 READ (INPUT,INFMT,END=50) (BUF(I),I=1,IDET(INTYPE)) 1497 ENDIF 1498c 15-Dec-1999 This was supposed to check for all-zero line and exit 1499c but all zeros is legal for [at least!] LABELs 1500 IF (INTYPE.EQ.LABEL) GOTO 9 1501 DO I=1,IDET(INTYPE) 1502 IF (BUF(I).NE.0.) GO TO 9 1503 ENDDO 1504 GO TO 50 15059 CONTINUE 1506* 1507* Expand array DETAIL if needed: NDET+KDET(INTYPE) > size(DETAIL) 1508 NEEDMEM = NDET+KDET(INTYPE) 1509 if ( NEEDMEM .GT. SIZE(DETAIL) ) THEN 1510 CALL GET_TRY(SIZE(DETAIL), NEEDMEM, TRY1, 1) 1511 DO 920 ITRY = 1,3 1512 if (TRY1(ITRY).LE.0 .OR. TRY1(ITRY).GT.MAXMEM) GOTO 920 1513 ALLOCATE( TMP1DR(TRY1(ITRY)), stat=ierr ) 1514 if (ierr .NE. 0) GOTO 920 1515 TMP1DR = 0. 1516 TMP1DR = DETAIL 1517 if(TEST_ALLOC)write(NOISE,*)"Expand DETAIL to ",try1(ITRY) 1518 CALL MOVE_ALLOC(from=TMP1DR, to=DETAIL) 1519 GOTO 922 1520920 CONTINUE 1521 ENDIF 1522922 CALL ASSERT (NEEDMEM.LE.SIZE(DETAIL), 1523 & 'too many object details - increase MAXDET and recompile') 1524 IF (SHADOW) THEN 1525* Expand array SDTAIL if needed 1526 NEEDMEM = MDET+SDET(INTYPE) 1527 IF ( NEEDMEM .GT. SIZE(SDTAIL) ) THEN 1528 CALL GET_TRY(SIZE(SDTAIL), NEEDMEM, TRY1, 1) 1529 DO 925 ITRY = 1,3 1530 if (TRY1(ITRY).LE.0.OR.TRY1(ITRY).GT.MAXMEM) GOTO 925 1531 ALLOCATE( TMP1DR(TRY1(ITRY)), stat=ierr ) 1532 if (ierr .NE. 0) GOTO 925 1533 TMP1DR = 0. 1534 TMP1DR = SDTAIL 1535 CALL MOVE_ALLOC(from=TMP1DR, to=SDTAIL) 1536 if(TEST_ALLOC)write(NOISE,*)"Expand SDTAIL to ", 1537 & try1(ITRY) 1538 GOTO 927 1539925 CONTINUE 1540 ENDIF 1541927 CALL ASSERT (NEEDMEM.LE.SIZE(SDTAIL), 1542 & 'too many shadow object details - increase MAXSDT and recompile') 1543 ENDIF 1544 N = N + 1 1545* Expand arrays ZTEMP, ZINDEX, LIST, MIST, TYPE, FLAG if needed 1546 IF (N.GT.SIZE(FLAG)) THEN 1547 CALL GET_TRY(SIZE(FLAG), N, TRY1, 1) 1548 DO 930 ITRY = 1,3 1549 if (TRY1(ITRY).LE.0.OR.TRY1(ITRY).GT.MAXMEM) GOTO 930 1550 ALLOCATE( TMP1DR(TRY1(ITRY)), stat=ierr ) 1551 if (ierr .NE. 0) GOTO 930 1552 TMP1DR = 0. 1553 TMP1DR = ZTEMP 1554 CALL MOVE_ALLOC(from=TMP1DR, to=ZTEMP) 1555 ALLOCATE( TMP1D(TRY1(ITRY)), stat=ierr ) 1556 if (ierr .NE. 0) GOTO 930 1557 TMP1D = 0 1558 TMP1D = ZINDEX 1559 CALL MOVE_ALLOC(from=TMP1D, to=ZINDEX) 1560 ALLOCATE( TMP1D(TRY1(ITRY)), stat=ierr ) 1561 if (ierr .NE. 0) GOTO 930 1562 TMP1D = 0 1563 TMP1D = LIST 1564 CALL MOVE_ALLOC(from=TMP1D, to=LIST) 1565 ALLOCATE( TMP1D(TRY1(ITRY)), stat=ierr ) 1566 if (ierr .NE. 0) GOTO 930 1567 TMP1D = 0 1568 TMP1D = MIST 1569 CALL MOVE_ALLOC(from=TMP1D, to=MIST) 1570 ALLOCATE( TMP1D(TRY1(ITRY)), stat=ierr ) 1571 if (ierr .NE. 0) GOTO 930 1572 TMP1D = 0 1573 TMP1D = TYPE 1574 CALL MOVE_ALLOC(from=TMP1D, to=TYPE) 1575 ALLOCATE( TMP1DI4(TRY1(ITRY)), stat=ierr ) 1576 if (ierr .NE. 0) GOTO 930 1577 TMP1DI4 = 0 1578 TMP1DI4 = FLAG 1579 CALL MOVE_ALLOC(from=TMP1DI4, to=FLAG) 1580 if(TEST_ALLOC)write(NOISE,*)"Expand MAXOBJ to ",size(FLAG) 1581 GOTO 932 1582930 CONTINUE 1583 ENDIF 1584 1585932 CALL ASSERT (N.LE.SIZE(FLAG), 1586 & 'too many objects - increase MAXOBJ and recompile') 1587C 20-Feb-1997 Save both object type and material type 1588 TYPE(N) = INTYPE 1589 IF (MSTATE.EQ.MATERIAL) FLAG(N) = FLAG(N) + 65536 * NPROPM 1590 LIST(N) = NDET 1591 IF (SHADOW) MIST(N) = MDET 1592 ISTRANS = 0 1593* From this point on, we'll use the symbolic codes for objects 1594 IF (INTYPE.EQ.TRIANG .or. INTYPE.EQ.PLANE) THEN 1595* triangle as read in 1596 X1A = BUF(1) 1597 Y1A = BUF(2) 1598 Z1A = BUF(3) 1599 X2A = BUF(4) 1600 Y2A = BUF(5) 1601 Z2A = BUF(6) 1602 X3A = BUF(7) 1603 Y3A = BUF(8) 1604 Z3A = BUF(9) 1605 RED = BUF(10) 1606 GRN = BUF(11) 1607 BLU = BUF(12) 1608 CALL CHKRGB( RED,GRN,BLU,'invalid triangle color' ) 1609 CALL ASSERT (IDET(INTYPE).EQ.12,'idet(1).ne.12') 1610 IF (MSTATE.EQ.MATERIAL) THEN 1611 FLAG(N) = ior(FLAG(N),PROPS) 1612 NPROPS = NPROPS + 1 1613 IF (CLRITY.GT.0) THEN 1614 FLAG(N) = ior(FLAG(N),TRANSP) 1615 IF (CLROPT.EQ.1) FLAG(N) = ior(FLAG(N),MOPT1) 1616 NTRANSP = NTRANSP + 1 1617 ISTRANS = 1 1618 ENDIF 1619 IF (CLIPPING) THEN 1620 FLAG(N) = ior(FLAG(N),CLIPPED) 1621 ENDIF 1622 IF (NBOUNDS.GT.0) FLAG(N) = ior(FLAG(N),BOUNDED) 1623 IF (MATCOL) THEN 1624 RED = RGBMAT(1) 1625 GRN = RGBMAT(2) 1626 BLU = RGBMAT(3) 1627 ENDIF 1628 ENDIF 1629* Isolated objects not transformed by TMAT, but still subject to inversion 1630 IF (ISOLATION.GT.0) THEN 1631 CALL ISOLATE(X1A,Y1A) 1632 CALL ISOLATE(X2A,Y2A) 1633 CALL ISOLATE(X3A,Y3A) 1634 ELSE 1635* update true coordinate limits 1636 TRULIM(1,1) = MIN( TRULIM(1,1), X1A,X2A,X3A) 1637 TRULIM(1,2) = MAX( TRULIM(1,2), X1A,X2A,X3A) 1638 TRULIM(2,1) = MIN( TRULIM(2,1), Y1A,Y2A,Y3A) 1639 TRULIM(2,2) = MAX( TRULIM(2,2), Y1A,Y2A,Y3A) 1640 TRULIM(3,1) = MIN( TRULIM(3,1), Z1A,Z2A,Z3A) 1641 TRULIM(3,2) = MAX( TRULIM(3,2), Z1A,Z2A,Z3A) 1642* modify the input, so to speak 1643 CALL TRANSF (X1A,Y1A,Z1A) 1644 CALL TRANSF (X2A,Y2A,Z2A) 1645 CALL TRANSF (X3A,Y3A,Z3A) 1646 ENDIF 1647* perspective factor for each corner 1648 IF (EYEPOS.GT.0) THEN 1649 PFAC1 = PERSP( Z1A ) 1650 PFAC2 = PERSP( Z2A ) 1651 PFAC3 = PERSP( Z3A ) 1652 END IF 1653* apply perspective 1654 X1B = X1A * PFAC1 1655 Y1B = Y1A * PFAC1 1656 Z1B = Z1A * PFAC1 1657 X2B = X2A * PFAC2 1658 Y2B = Y2A * PFAC2 1659 Z2B = Z2A * PFAC2 1660 X3B = X3A * PFAC3 1661 Y3B = Y3A * PFAC3 1662 Z3B = Z3A * PFAC3 1663* scale and translate to pixel space 1664 X1C = X1B * SCALE + XCENT 1665 Y1C = Y1B * SCALE + YCENT 1666 Z1C = Z1B * SCALE 1667 X2C = X2B * SCALE + XCENT 1668 Y2C = Y2B * SCALE + YCENT 1669 Z2C = Z2B * SCALE 1670 X3C = X3B * SCALE + XCENT 1671 Y3C = Y3B * SCALE + YCENT 1672 Z3C = Z3B * SCALE 1673* save transformed Z limits 1674 ZLIM(1) = MIN( ZLIM(1), Z1C,Z2C,Z3C ) 1675 ZLIM(2) = MAX( ZLIM(2), Z1C,Z2C,Z3C ) 1676* 1677* check for Z-clipping 1678 JUSTCLIPPED = .FALSE. 1679 IF (INTYPE.NE.PLANE) THEN 1680 IF ((MIN(Z1C,Z2C,Z3C) .GT. FRONTCLIP) 1681 & .OR.(MAX(Z1C,Z2C,Z3C) .LT. BACKCLIP)) THEN 1682 JUSTCLIPPED = .TRUE. 1683 GOTO 45 1684 ENDIF 1685 IF (CLIPPING) THEN 1686 MIND = LIST(MLIST(NPROPM)) 1687 IF ((MIN(Z1C,Z2C,Z3C) .GT. DETAIL(MIND+16)) 1688 & .OR.(MAX(Z1C,Z2C,Z3C) .LT. DETAIL(MIND+17))) THEN 1689 JUSTCLIPPED = .TRUE. 1690 GOTO 45 1691 ENDIF 1692 ENDIF 1693 ENDIF 1694* 1695* solve for coefficients of plane eqn z=Ax+By+C 1696 CALL PLANER(X1C,Y1C,Z1C, 1697 & X2C,Y2C,Z2C, 1698 & X3C,Y3C,Z3C, A,B,C,D) 1699* save results for PLANE object 1700* PLANE impinges on all tiles, but casts no shadows 1701 IF (INTYPE.EQ. PLANE) THEN 1702 NPLANES = NPLANES + 1 1703 DETAIL(NDET+1) = A 1704 DETAIL(NDET+2) = B 1705 DETAIL(NDET+3) = C 1706 DETAIL(NDET+4) = D 1707 DETAIL(NDET+5) = RED 1708 DETAIL(NDET+6) = GRN 1709 DETAIL(NDET+7) = BLU 1710 CALL ASSERT(KDET(INTYPE).EQ.7,'kdet(6).ne.7') 1711 NDET = NDET + KDET(INTYPE) 1712 DO IY = 1,NTY 1713 DO IX = 1,NTX 1714 KOUNT(IX,IY) = KOUNT(IX,IY) + 1 1715 TTRANS(IX,IY) = TTRANS(IX,IY) + ISTRANS 1716 ENDDO 1717 ENDDO 1718C write(NOISE,*)"Incr kount ",NTX," x ",NTY 1719 IF (SHADOW) THEN 1720 MDET = MDET + SDET(INTYPE) 1721 ENDIF 1722 GOTO 7 1723 ENDIF 1724* save results for normal triangles 1725 DETAIL(NDET+1) = X1C 1726 DETAIL(NDET+2) = Y1C 1727 DETAIL(NDET+3) = Z1C 1728 DETAIL(NDET+4) = X2C 1729 DETAIL(NDET+5) = Y2C 1730 DETAIL(NDET+6) = Z2C 1731 DETAIL(NDET+7) = X3C 1732 DETAIL(NDET+8) = Y3C 1733 DETAIL(NDET+9) = Z3C 1734 DETAIL(NDET+10) = A 1735 DETAIL(NDET+11) = B 1736 DETAIL(NDET+12) = C 1737 DETAIL(NDET+13) = D 1738 DETAIL(NDET+14) = RED 1739 DETAIL(NDET+15) = GRN 1740 DETAIL(NDET+16) = BLU 1741 CALL ASSERT (KDET(INTYPE).EQ.16,'kdet(1).ne.16') 1742 NDET = NDET + KDET(INTYPE) 1743* tally for tiles the object might impinge on 1744 IXLO = MIN(X1C,X2C,X3C)/NPX + 1 1745 IXHI = MAX(X1C,X2C,X3C)/NPX + 1 1746 IYLO = MIN(Y1C,Y2C,Y3C)/NPY + 1 1747 IYHI = MAX(Y1C,Y2C,Y3C)/NPY + 1 1748 IF (IXLO.LT.1 ) IXLO=1 1749 IF (IXLO.GT.NTX) GO TO 11 1750 IF (IXHI.LT.1 ) GO TO 11 1751 IF (IXHI.GT.NTX) IXHI=NTX 1752 IF (IYLO.LT.1 ) IYLO=1 1753 IF (IYLO.GT.NTY) GO TO 11 1754 IF (IYHI.LT.1 ) GO TO 11 1755 IF (IYHI.GT.NTY) IYHI=NTY 1756 DO 10 IY=IYLO,IYHI 1757 DO 10 IX=IXLO,IXHI 1758 KOUNT(IX,IY) = KOUNT(IX,IY) + 1 1759 TTRANS(IX,IY) = TTRANS(IX,IY) + ISTRANS 176010 CONTINUE 1761C write(NOISE,*)"Incr kount b ",IXLO,"-",IXHI," x ",IYLO,"-",IYHI 176211 CONTINUE 1763* repeat for shadow buffer if necessary 1764 IF (SHADOW) THEN 1765* rotate light source to z to take light source viewpoint 1766 X1R = SROT(1,1)*X1B+SROT(1,2)*Y1B+SROT(1,3)*Z1B 1767 Y1R = SROT(2,1)*X1B+SROT(2,2)*Y1B+SROT(2,3)*Z1B 1768 Z1R = SROT(3,1)*X1B+SROT(3,2)*Y1B+SROT(3,3)*Z1B 1769 X2R = SROT(1,1)*X2B+SROT(1,2)*Y2B+SROT(1,3)*Z2B 1770 Y2R = SROT(2,1)*X2B+SROT(2,2)*Y2B+SROT(2,3)*Z2B 1771 Z2R = SROT(3,1)*X2B+SROT(3,2)*Y2B+SROT(3,3)*Z2B 1772 X3R = SROT(1,1)*X3B+SROT(1,2)*Y3B+SROT(1,3)*Z3B 1773 Y3R = SROT(2,1)*X3B+SROT(2,2)*Y3B+SROT(2,3)*Z3B 1774 Z3R = SROT(3,1)*X3B+SROT(3,2)*Y3B+SROT(3,3)*Z3B 1775* scale and translate for shadow space 1776 X1S = X1R * SCALE + SXCENT 1777 Y1S = Y1R * SCALE + SYCENT 1778 Z1S = Z1R * SCALE 1779 X2S = X2R * SCALE + SXCENT 1780 Y2S = Y2R * SCALE + SYCENT 1781 Z2S = Z2R * SCALE 1782 X3S = X3R * SCALE + SXCENT 1783 Y3S = Y3R * SCALE + SYCENT 1784 Z3S = Z3R * SCALE 1785* solve plane eqn etc. 1786 CALL PLANER(X1S,Y1S,Z1S, 1787 & X2S,Y2S,Z2S, 1788 & X3S,Y3S,Z3S, A,B,C,D) 1789* save results etc. 1790 SDTAIL(MDET+1) = X1S 1791 SDTAIL(MDET+2) = Y1S 1792 SDTAIL(MDET+3) = Z1S 1793 SDTAIL(MDET+4) = X2S 1794 SDTAIL(MDET+5) = Y2S 1795 SDTAIL(MDET+6) = Z2S 1796 SDTAIL(MDET+7) = X3S 1797 SDTAIL(MDET+8) = Y3S 1798 SDTAIL(MDET+9) = Z3S 1799 SDTAIL(MDET+10) = A 1800 SDTAIL(MDET+11) = B 1801 SDTAIL(MDET+12) = C 1802 SDTAIL(MDET+13) = D 1803 CALL ASSERT (SDET(INTYPE).EQ.13,'sdet(1).ne.13') 1804 MDET = MDET + SDET(INTYPE) 1805* tally for shadow tiles the object might impinge on 1806 IXLO = MIN(X1S,X2S,X3S)/NPX + 1 1807 IXHI = MAX(X1S,X2S,X3S)/NPX + 1 1808 IYLO = MIN(Y1S,Y2S,Y3S)/NPY + 1 1809 IYHI = MAX(Y1S,Y2S,Y3S)/NPY + 1 1810 IF (IXLO.LT.1 ) IXLO=1 1811 IF (IXLO.GT.NSX) GO TO 16 1812 IF (IXHI.LT.1 ) GO TO 16 1813 IF (IXHI.GT.NSX) IXHI=NSX 1814 IF (IYLO.LT.1 ) IYLO=1 1815 IF (IYLO.GT.NSY) GO TO 16 1816 IF (IYHI.LT.1 ) GO TO 16 1817 IF (IYHI.GT.NSY) IYHI=NSY 1818 DO 15 IY=IYLO,IYHI 1819 DO 15 IX=IXLO,IXHI 1820 MOUNT(IX,IY) = MOUNT(IX,IY) + 1 182115 CONTINUE 182216 CONTINUE 1823 ENDIF 1824 ELSEIF (INTYPE.EQ.SPHERE) THEN 1825* sphere as read in 1826 XA = BUF(1) 1827 YA = BUF(2) 1828 ZA = BUF(3) 1829 RA = BUF(4) 1830 RED = BUF(5) 1831 GRN = BUF(6) 1832 BLU = BUF(7) 1833 CALL CHKRGB (RED,GRN,BLU,'invalid sphere color') 1834 CALL ASSERT (IDET(INTYPE).EQ.7,'idet(2).ne.7') 1835 IF (MSTATE.EQ.MATERIAL) THEN 1836 FLAG(N) = ior(FLAG(N),PROPS) 1837 NPROPS = NPROPS + 1 1838 IF (CLRITY.GT.0) THEN 1839 FLAG(N) = ior(FLAG(N),TRANSP) 1840 IF (CLROPT.EQ.1) FLAG(N) = ior(FLAG(N),MOPT1) 1841 NTRANSP = NTRANSP + 1 1842 ISTRANS = 1 1843 ENDIF 1844 IF (CLIPPING) THEN 1845 FLAG(N) = ior(FLAG(N),CLIPPED) 1846 ENDIF 1847 IF (NBOUNDS.GT.0) FLAG(N) = ior(FLAG(N),BOUNDED) 1848 IF (MATCOL) THEN 1849 RED = RGBMAT(1) 1850 GRN = RGBMAT(2) 1851 BLU = RGBMAT(3) 1852 ENDIF 1853 ENDIF 1854* Isolated objects not transformed by TMAT, but still subject to inversion 1855 IF (ISOLATION.GT.0) THEN 1856 CALL ISOLATE(XA,YA) 1857 ELSE 1858* update true coordinate limits 1859 TRULIM(1,1) = MIN( TRULIM(1,1), XA ) 1860 TRULIM(1,2) = MAX( TRULIM(1,2), XA ) 1861 TRULIM(2,1) = MIN( TRULIM(2,1), YA ) 1862 TRULIM(2,2) = MAX( TRULIM(2,2), YA ) 1863 TRULIM(3,1) = MIN( TRULIM(3,1), ZA ) 1864 TRULIM(3,2) = MAX( TRULIM(3,2), ZA ) 1865* modify the input, as it were 1866 CALL TRANSF (XA,YA,ZA) 1867 RA = RA / TMAT(4,4) 1868 ENDIF 1869* perspective 1870 IF (EYEPOS.GT.0) PFAC = PERSP(ZA) 1871 XB = XA * PFAC 1872 YB = YA * PFAC 1873 ZB = ZA * PFAC 1874 RB = RA * PFAC 1875* scale & translate 1876 XC = XB * SCALE + XCENT 1877 YC = YB * SCALE + YCENT 1878 ZC = ZB * SCALE 1879 RC = RB * SCALE 1880* save transformed Z limits 1881 ZLIM(1) = MIN( ZLIM(1), ZC ) 1882 ZLIM(2) = MAX( ZLIM(2), ZC ) 1883* check for Z-clipping 1884 IF (ZC .GT. FRONTCLIP .OR. ZC .LT. BACKCLIP) THEN 1885 JUSTCLIPPED = .TRUE. 1886 GOTO 45 1887 ELSE 1888 JUSTCLIPPED = .FALSE. 1889 ENDIF 1890* save results 1891 DETAIL(NDET+1) = XC 1892 DETAIL(NDET+2) = YC 1893 DETAIL(NDET+3) = ZC 1894 DETAIL(NDET+4) = RC 1895 DETAIL(NDET+5) = RED 1896 DETAIL(NDET+6) = GRN 1897 DETAIL(NDET+7) = BLU 1898 CALL ASSERT (KDET(INTYPE).EQ.7,'kdet(2).ne.7') 1899 NDET = NDET + KDET(INTYPE) 1900 nsphere = nsphere + 1 1901* tally for tiles the object might impinge on 1902 IXLO = (XC-RC)/NPX + 1 1903 IXHI = (XC+RC)/NPX + 1 1904 IYLO = (YC-RC)/NPY + 1 1905 IYHI = (YC+RC)/NPY + 1 1906 IF (IXLO.LT.1 ) IXLO=1 1907 IF (IXLO.GT.NTX) GO TO 21 1908 IF (IXHI.LT.1 ) GO TO 21 1909 IF (IXHI.GT.NTX) IXHI=NTX 1910 IF (IYLO.LT.1 ) IYLO=1 1911 IF (IYLO.GT.NTY) GO TO 21 1912 IF (IYHI.LT.1 ) GO TO 21 1913 IF (IYHI.GT.NTY) IYHI=NTY 1914 DO 20 IY=IYLO,IYHI 1915 DO 20 IX=IXLO,IXHI 1916 KOUNT(IX,IY) = KOUNT(IX,IY) + 1 1917 TTRANS(IX,IY) = TTRANS(IX,IY) + ISTRANS 191820 CONTINUE 1919C write(NOISE,*)"Incr kount c ",IXLO,"-",IXHI," x ",IYLO,"-",IYHI, 1920C & kstart(IXLO,IYLO),kstop(IXLO,IYLO),kount(IXLO,IYLO) 192121 CONTINUE 1922* repeat for shadow buffer if necessary 1923 IF (SHADOW) THEN 1924* rotate light source to z to take light source viewpoint 1925 XR = SROT(1,1)*XB+SROT(1,2)*YB+SROT(1,3)*ZB 1926 YR = SROT(2,1)*XB+SROT(2,2)*YB+SROT(2,3)*ZB 1927 ZR = SROT(3,1)*XB+SROT(3,2)*YB+SROT(3,3)*ZB 1928 RR = RB 1929* scale and translate for shadow space 1930 XS = XR * SCALE + SXCENT 1931 YS = YR * SCALE + SYCENT 1932 ZS = ZR * SCALE 1933 RS = RR * SCALE 1934* save results 1935 SDTAIL(MDET+1) = XS 1936 SDTAIL(MDET+2) = YS 1937 SDTAIL(MDET+3) = ZS 1938 SDTAIL(MDET+4) = RS 1939 CALL ASSERT (SDET(INTYPE).EQ.4,'sdet(2).ne.4') 1940 MDET = MDET + SDET(INTYPE) 1941* tally for shadow tiles the object might impinge on 1942 IXLO = (XS-RS)/NPX + 1 1943 IXHI = (XS+RS)/NPX + 1 1944 IYLO = (YS-RS)/NPY + 1 1945 IYHI = (YS+RS)/NPY + 1 1946 IF (IXLO.LT.1 ) IXLO=1 1947 IF (IXLO.GT.NSX) GO TO 26 1948 IF (IXHI.LT.1 ) GO TO 26 1949 IF (IXHI.GT.NSX) IXHI=NSX 1950 IF (IYLO.LT.1 ) IYLO=1 1951 IF (IYLO.GT.NSY) GO TO 26 1952 IF (IYHI.LT.1 ) GO TO 26 1953 IF (IYHI.GT.NSY) IYHI=NSY 1954 DO 25 IY=IYLO,IYHI 1955 DO 25 IX=IXLO,IXHI 1956 MOUNT(IX,IY) = MOUNT(IX,IY) + 1 195725 CONTINUE 195826 CONTINUE 1959 ENDIF 1960 1961c 30-Dec-99 duplicate transparent spheres, if requested, so that 1962c the inside surface can be rendered also. BUF() is still loaded 1963c with specs for the current object; we just need to set flags. 1964 IF ( iand(FLAG(N),TRANSP).NE.0 .AND. CLROPT.EQ.2 1965 & .AND. iand(FLAG(N),INSIDE).EQ.0) THEN 1966 FLAG(N+1) = INSIDE 1967 NINSIDE = NINSIDE + 1 1968 GOTO 9 1969 ENDIF 1970 1971 ELSEIF (INTYPE.EQ.CYLIND) THEN 1972* EAM May 1990 cylinder as read in 1973 X1A = BUF(1) 1974 Y1A = BUF(2) 1975 Z1A = BUF(3) 1976 R1A = BUF(4) 1977 X2A = BUF(5) 1978 Y2A = BUF(6) 1979 Z2A = BUF(7) 1980 R2A = R1A 1981 RED = BUF(9) 1982 GRN = BUF(10) 1983 BLU = BUF(11) 1984 CALL CHKRGB (RED,GRN,BLU,'invalid cylinder color') 1985 CALL ASSERT (IDET(INTYPE).EQ.11,'idet(1).ne.11') 1986* Zero length cylinder is better treated as sphere 1987* EAM 22-Nov-96 1988 IF ((iand(FLAG(N),FLAT).EQ.0) .AND. 1989 & (X1A.EQ.X2A).AND.(Y1A.EQ.Y2A).AND.(Z1A.EQ.Z2A)) THEN 1990 BUF(5) = BUF(9) 1991 BUF(6) = BUF(10) 1992 BUF(7) = BUF(11) 1993 INTYPE = SPHERE 1994 N = N-1 1995 GOTO 9 1996 ENDIF 1997 IF (MSTATE.EQ.MATERIAL) THEN 1998 FLAG(N) = ior(FLAG(N),PROPS) 1999 NPROPS = NPROPS + 1 2000 IF (CLRITY.GT.0) THEN 2001 FLAG(N) = ior(FLAG(N),TRANSP) 2002 IF (CLROPT.EQ.1) FLAG(N) = ior(FLAG(N),MOPT1) 2003 NTRANSP = NTRANSP + 1 2004 ISTRANS = 1 2005 ENDIF 2006 IF (CLIPPING) THEN 2007 FLAG(N) = ior(FLAG(N),CLIPPED) 2008 ENDIF 2009 IF (NBOUNDS.GT.0) FLAG(N) = ior(FLAG(N),BOUNDED) 2010 IF (MATCOL) THEN 2011 RED = RGBMAT(1) 2012 GRN = RGBMAT(2) 2013 BLU = RGBMAT(3) 2014 ENDIF 2015 ENDIF 2016* Isolated objects not transformed by TMAT, but still subject to inversion 2017 IF (ISOLATION.GT.0) THEN 2018 CALL ISOLATE(X1A,Y1A) 2019 CALL ISOLATE(X2A,Y2A) 2020 ELSE 2021* update true coordinate limits 2022 TRULIM(1,1) = MIN( TRULIM(1,1), X1A,X2A) 2023 TRULIM(1,2) = MAX( TRULIM(1,2), X1A,X2A) 2024 TRULIM(2,1) = MIN( TRULIM(2,1), Y1A,Y2A) 2025 TRULIM(2,2) = MAX( TRULIM(2,2), Y1A,Y2A) 2026 TRULIM(3,1) = MIN( TRULIM(3,1), Z1A,Z2A) 2027 TRULIM(3,2) = MAX( TRULIM(3,2), Z1A,Z2A) 2028* modify the input, so to speak 2029 CALL TRANSF (X1A,Y1A,Z1A) 2030 CALL TRANSF (X2A,Y2A,Z2A) 2031 R1A = R1A / TMAT(4,4) 2032 R2A = R2A / TMAT(4,4) 2033 ENDIF 2034* perspective factor for each corner 2035 IF (EYEPOS.GT.0) THEN 2036 PFAC1 = PERSP( Z1A ) 2037 PFAC2 = PERSP( Z2A ) 2038 END IF 2039* apply perspective 2040 X1B = X1A * PFAC1 2041 Y1B = Y1A * PFAC1 2042 Z1B = Z1A * PFAC1 2043 R1B = R1A * PFAC1 2044 X2B = X2A * PFAC2 2045 Y2B = Y2A * PFAC2 2046 Z2B = Z2A * PFAC2 2047 R2B = R2A * PFAC2 2048* scale and translate to pixel space 2049 X1C = X1B * SCALE + XCENT 2050 Y1C = Y1B * SCALE + YCENT 2051 Z1C = Z1B * SCALE 2052 R1C = R1B * SCALE 2053 X2C = X2B * SCALE + XCENT 2054 Y2C = Y2B * SCALE + YCENT 2055 Z2C = Z2B * SCALE 2056 R2C = R2B * SCALE 2057* save transformed Z limits 2058 ZLIM(1) = MIN( ZLIM(1), Z1C,Z2C ) 2059 ZLIM(2) = MAX( ZLIM(2), Z1C,Z2C ) 2060* check for Z-clipping 2061 IF ((MIN(Z1C,Z2C) .GT. FRONTCLIP) 2062 & .OR.(MAX(Z1C,Z2C) .LT. BACKCLIP)) THEN 2063 JUSTCLIPPED = .TRUE. 2064 GOTO 45 2065 ELSE 2066 JUSTCLIPPED = .FALSE. 2067 ENDIF 2068* save results 2069 DETAIL(NDET+1) = X1C 2070 DETAIL(NDET+2) = Y1C 2071 DETAIL(NDET+3) = Z1C 2072 DETAIL(NDET+4) = R1C 2073 DETAIL(NDET+5) = X2C 2074 DETAIL(NDET+6) = Y2C 2075 DETAIL(NDET+7) = Z2C 2076 DETAIL(NDET+8) = R2C 2077* EAM save anything else? 2078 DETAIL(NDET+9) = RED 2079 DETAIL(NDET+10) = GRN 2080 DETAIL(NDET+11) = BLU 2081 CALL ASSERT (KDET(INTYPE).EQ.11,'kdet(1).ne.11') 2082 NDET = NDET + KDET(INTYPE) 2083 ncylind = ncylind + 1 2084* tally for tiles the object might impinge on 2085 IXLO = MIN(X1C-R1C,X2C-R2C) / NPX + 1 2086 IXHI = MAX(X1C+R1C,X2C+R2C) / NPX + 1 2087 IYLO = MIN(Y1C-R1C,Y2C-R2C) / NPY + 1 2088 IYHI = MAX(Y1C+R1C,Y2C+R2C) / NPY + 1 2089 IF (IXLO.LT.1 ) IXLO=1 2090 IF (IXLO.GT.NTX) GO TO 711 2091 IF (IXHI.LT.1 ) GO TO 711 2092 IF (IXHI.GT.NTX) IXHI=NTX 2093 IF (IYLO.LT.1 ) IYLO=1 2094 IF (IYLO.GT.NTY) GO TO 711 2095 IF (IYHI.LT.1 ) GO TO 711 2096 IF (IYHI.GT.NTY) IYHI=NTY 2097 DO 710 IY=IYLO,IYHI 2098 DO 710 IX=IXLO,IXHI 2099 KOUNT(IX,IY) = KOUNT(IX,IY) + 1 2100 TTRANS(IX,IY) = TTRANS(IX,IY) + ISTRANS 2101710 CONTINUE 2102C write(NOISE,*)"Incr kount d ",IXLO,"-",IXHI," x ",IYLO,"-",IYHI 2103711 CONTINUE 2104* repeat for shadow buffer if necessary 2105 IF (SHADOW) THEN 2106* rotate light source to z to take light source viewpoint 2107 X1R = SROT(1,1)*X1B+SROT(1,2)*Y1B+SROT(1,3)*Z1B 2108 Y1R = SROT(2,1)*X1B+SROT(2,2)*Y1B+SROT(2,3)*Z1B 2109 Z1R = SROT(3,1)*X1B+SROT(3,2)*Y1B+SROT(3,3)*Z1B 2110 X2R = SROT(1,1)*X2B+SROT(1,2)*Y2B+SROT(1,3)*Z2B 2111 Y2R = SROT(2,1)*X2B+SROT(2,2)*Y2B+SROT(2,3)*Z2B 2112 Z2R = SROT(3,1)*X2B+SROT(3,2)*Y2B+SROT(3,3)*Z2B 2113* scale and translate for shadow space 2114 X1S = X1R * SCALE + SXCENT 2115 Y1S = Y1R * SCALE + SYCENT 2116 Z1S = Z1R * SCALE 2117 R1S = R1B * SCALE 2118 X2S = X2R * SCALE + SXCENT 2119 Y2S = Y2R * SCALE + SYCENT 2120 Z2S = Z2R * SCALE 2121 R2S = R2B * SCALE 2122* save results etc. 2123 SDTAIL(MDET+1) = X1S 2124 SDTAIL(MDET+2) = Y1S 2125 SDTAIL(MDET+3) = Z1S 2126 SDTAIL(MDET+4) = R1S 2127 SDTAIL(MDET+5) = X2S 2128 SDTAIL(MDET+6) = Y2S 2129 SDTAIL(MDET+7) = Z2S 2130 SDTAIL(MDET+8) = R2S 2131 CALL ASSERT (SDET(INTYPE).EQ.8,'sdet(1).ne.8') 2132 MDET = MDET + SDET(INTYPE) 2133* tally for shadow tiles the object might impinge on 2134 IXLO = MIN(X1S-R1S,X2S-R2S) / NPX + 1 2135 IXHI = MAX(X1S+R1S,X2S+R2S) / NPX + 1 2136 IYLO = MIN(Y1S-R1S,Y2S-R2S) / NPY + 1 2137 IYHI = MAX(Y1S+R1S,Y2S+R2S) / NPY + 1 2138 IF (IXLO.LT.1 ) IXLO=1 2139 IF (IXLO.GT.NSX) GO TO 716 2140 IF (IXHI.LT.1 ) GO TO 716 2141 IF (IXHI.GT.NSX) IXHI=NSX 2142 IF (IYLO.LT.1 ) IYLO=1 2143 IF (IYLO.GT.NSY) GO TO 716 2144 IF (IYHI.LT.1 ) GO TO 716 2145 IF (IYHI.GT.NSY) IYHI=NSY 2146 DO 715 IY=IYLO,IYHI 2147 DO 715 IX=IXLO,IXHI 2148 MOUNT(IX,IY) = MOUNT(IX,IY) + 1 2149715 CONTINUE 2150716 CONTINUE 2151 ENDIF 2152c 2153c 20-Aug-98 duplicate any transparent flat-ended cylinders so that 2154c the inside surface can be rendered also. BUF() is still loaded 2155c with specs for the current object; we just need to set flags. 2156 IF ( iand(FLAG(N),TRANSP).NE.0 .AND. iand(FLAG(N),FLAT).NE.0 2157 & .AND. iand(FLAG(N),INSIDE).EQ.0) THEN 2158 FLAG(N+1) = FLAT + INSIDE 2159 NINSIDE = NINSIDE + 1 2160 GOTO 9 2161 ENDIF 2162* 2163 ELSEIF (INTYPE.EQ.NORMS) THEN 2164* vertex normals as given (these belong to previous triangle) 2165 IF (JUSTCLIPPED) GOTO 46 2166 IPREV = N - 1 2167 IF ((TYPE(IPREV).EQ.VERTEXRGB).OR.(TYPE(IPREV).EQ.VERTRANSP)) 2168 & IPREV = IPREV - 1 2169 IF ((TYPE(IPREV).EQ.VERTEXRGB).OR.(TYPE(IPREV).EQ.VERTRANSP)) 2170 & IPREV = IPREV - 1 2171 CALL ASSERT (TYPE(IPREV).EQ.TRIANG,'orphan normals') 2172* Isolated objects not transformed by TMAT, but still subject to inversion 2173 IF (ISOLATION.GT.0) THEN 2174 X1C = BUF(1) 2175 Y1C = BUF(2) 2176 Z1C = BUF(3) 2177 X2C = BUF(4) 2178 Y2C = BUF(5) 2179 Z2C = BUF(6) 2180 X3C = BUF(7) 2181 Y3C = BUF(8) 2182 Z3C = BUF(9) 2183 IF (INVERT) THEN 2184 Y1C = -Y1C 2185 Y2C = -Y2C 2186 Y3C = -Y3C 2187 ENDIF 2188 ELSE 2189 X1A = BUF(1) 2190 Y1A = BUF(2) 2191 Z1A = BUF(3) 2192 X2A = BUF(4) 2193 Y2A = BUF(5) 2194 Z2A = BUF(6) 2195 X3A = BUF(7) 2196 Y3A = BUF(8) 2197 Z3A = BUF(9) 2198* Apply rotation matrix, but not translation components 2199 X1B = X1A*TMAT(1,1) + Y1A*TMAT(2,1) + Z1A*TMAT(3,1) 2200 Y1B = X1A*TMAT(1,2) + Y1A*TMAT(2,2) + Z1A*TMAT(3,2) 2201 Z1B = X1A*TMAT(1,3) + Y1A*TMAT(2,3) + Z1A*TMAT(3,3) 2202 X2B = X2A*TMAT(1,1) + Y2A*TMAT(2,1) + Z2A*TMAT(3,1) 2203 Y2B = X2A*TMAT(1,2) + Y2A*TMAT(2,2) + Z2A*TMAT(3,2) 2204 Z2B = X2A*TMAT(1,3) + Y2A*TMAT(2,3) + Z2A*TMAT(3,3) 2205 X3B = X3A*TMAT(1,1) + Y3A*TMAT(2,1) + Z3A*TMAT(3,1) 2206 Y3B = X3A*TMAT(1,2) + Y3A*TMAT(2,2) + Z3A*TMAT(3,2) 2207 Z3B = X3A*TMAT(1,3) + Y3A*TMAT(2,3) + Z3A*TMAT(3,3) 2208* Also apply post-rotation, if any 2209 X1C = RAFTER(1,1)*X1B + RAFTER(1,2)*Y1B + RAFTER(1,3)*Z1B 2210 Y1C = RAFTER(2,1)*X1B + RAFTER(2,2)*Y1B + RAFTER(2,3)*Z1B 2211 Z1C = RAFTER(3,1)*X1B + RAFTER(3,2)*Y1B + RAFTER(3,3)*Z1B 2212 X2C = RAFTER(1,1)*X2B + RAFTER(1,2)*Y2B + RAFTER(1,3)*Z2B 2213 Y2C = RAFTER(2,1)*X2B + RAFTER(2,2)*Y2B + RAFTER(2,3)*Z2B 2214 Z2C = RAFTER(3,1)*X2B + RAFTER(3,2)*Y2B + RAFTER(3,3)*Z2B 2215 X3C = RAFTER(1,1)*X3B + RAFTER(1,2)*Y3B + RAFTER(1,3)*Z3B 2216 Y3C = RAFTER(2,1)*X3B + RAFTER(2,2)*Y3B + RAFTER(2,3)*Z3B 2217 Z3C = RAFTER(3,1)*X3B + RAFTER(3,2)*Y3B + RAFTER(3,3)*Z3B 2218 ENDIF 2219C 2220C If all three Z components are negative, it's facing away from us. 2221C Old default treatment was to assume the normals were screwed up. 2222C V2.6: default appropriate cases (e.g. opaque triangles with no 2223C associated bounding planes) to hidden by assumption and 2224C thus not needing to be rendered. Mark as HIDDEN. 2225C Default treatment is overridden by CLROPT in material spec 2226C V3.0.3: Render both sides of opaque triangles 2227C 2228 IF (Z1C.GE.0 .AND. Z2C.GE.0 .AND. Z3C.GE.0) GOTO 718 2229 IF (Z1C.LT.-.01 .AND. Z2C.LT.-.01 .AND. Z3C.LT.-.01) THEN 2230 IF (CLROPT.EQ.2) THEN 2231 NINSIDE = NINSIDE + 1 2232 FLAG(IPREV) = ior( FLAG(IPREV), INSIDE ) 2233 ELSE IF ((iand(FLAG(IPREV),BOUNDED).EQ.0) 2234 & .AND.(CLRITY.NE.0 .AND. CLROPT.EQ.1)) THEN 2235 NHIDDEN = NHIDDEN + 1 2236 FLAG(IPREV) = ior( FLAG(IPREV), HIDDEN ) 2237 ELSE 2238 NINSIDE = NINSIDE + 1 2239 FLAG(IPREV) = ior( FLAG(IPREV), INSIDE ) 2240 ENDIF 2241 GOTO 718 2242 ENDIF 2243C 2244C Mixed + and - Z means the triangle "wrapped around" the edge. 2245C For solid objects the best we can do is pretend the edge is right here. 2246C For transparent objects or 2-sided surfaces we need to invert the 2247C normals also. The value of EDGESLOP is purely empirical; setting it 2248C either too low or too high makes some edges get coloured wrongly. 2249C Setting the HIDDEN flag for this record (NB: for the NORMALS, not for 2250C the triangle itself) causes the triangle to have flat shading. 2251 IF (Z1C+Z2C+Z3C .LT. 0) THEN 2252 IF (Z1C .GT. EDGESLOP) FLAG(N) = HIDDEN 2253 IF (Z2C .GT. EDGESLOP) FLAG(N) = HIDDEN 2254 IF (Z3C .GT. EDGESLOP) FLAG(N) = HIDDEN 2255 Z1C = MIN(Z1C,0.) 2256 Z2C = MIN(Z2C,0.) 2257 Z3C = MIN(Z3C,0.) 2258 ELSE 2259 IF (Z1C .LT. -EDGESLOP) FLAG(N) = HIDDEN 2260 IF (Z2C .LT. -EDGESLOP) FLAG(N) = HIDDEN 2261 IF (Z3C .LT. -EDGESLOP) FLAG(N) = HIDDEN 2262 Z1C = MAX(Z1C,0.) 2263 Z2C = MAX(Z2C,0.) 2264 Z3C = MAX(Z3C,0.) 2265 ENDIF 2266C 2267718 CONTINUE 2268 DETAIL(NDET+1) = X1C 2269 DETAIL(NDET+2) = Y1C 2270 DETAIL(NDET+3) = Z1C 2271 DETAIL(NDET+4) = X2C 2272 DETAIL(NDET+5) = Y2C 2273 DETAIL(NDET+6) = Z2C 2274 DETAIL(NDET+7) = X3C 2275 DETAIL(NDET+8) = Y3C 2276 DETAIL(NDET+9) = Z3C 2277 NDET = NDET + KDET(INTYPE) 2278 IF (SHADOW) THEN 2279 MDET = MDET + SDET(INTYPE) 2280 ENDIF 2281* 2282* Allow specification of RGB triple for each vertex of preceding 2283* triangle or cylinder. Overrides base RGB. 2284* Also overrides MATERIAL RGB, which is arguably a bug. 2285* 2286 ELSEIF (INTYPE .EQ. VERTEXRGB) THEN 2287 IF (JUSTCLIPPED) GOTO 46 2288 CALL CHKRGB(BUF(1),BUF(2),BUF(3),'invalid vertex color') 2289 CALL CHKRGB(BUF(4),BUF(5),BUF(6),'invalid vertex color') 2290 CALL CHKRGB(BUF(7),BUF(8),BUF(9),'invalid vertex color') 2291 IPREV = N - 1 2292 IF ((TYPE(IPREV).EQ.NORMS).OR.(TYPE(IPREV).EQ.VERTRANSP)) 2293 & IPREV = IPREV - 1 2294 IF ((TYPE(IPREV).EQ.NORMS).OR.(TYPE(IPREV).EQ.VERTRANSP)) 2295 & IPREV = IPREV - 1 2296c we should only see a SPHERE is if it's a collapsed cylinder 2297 IF (TYPE(IPREV).EQ.SPHERE) THEN 2298 K = LIST(IPREV) 2299 DETAIL(K+5) = BUF(1) 2300 DETAIL(K+6) = BUF(2) 2301 DETAIL(K+7) = BUF(3) 2302 GOTO 7 2303 ENDIF 2304 CALL ASSERT (TYPE(IPREV).EQ.TRIANG .OR. TYPE(IPREV).EQ.CYLIND, 2305 & 'orphan vertex colours') 2306 FLAG(IPREV) = ior( FLAG(IPREV), VCOLS ) 2307 DETAIL(NDET+1) = BUF(1) 2308 DETAIL(NDET+2) = BUF(2) 2309 DETAIL(NDET+3) = BUF(3) 2310 DETAIL(NDET+4) = BUF(4) 2311 DETAIL(NDET+5) = BUF(5) 2312 DETAIL(NDET+6) = BUF(6) 2313 DETAIL(NDET+7) = BUF(7) 2314 DETAIL(NDET+8) = BUF(8) 2315 DETAIL(NDET+9) = BUF(9) 2316 NDET = NDET + KDET(INTYPE) 2317 IF (SHADOW) THEN 2318 MDET = MDET + SDET(INTYPE) 2319 ENDIF 2320* 2321* EAM - 30-Dec-1999 2322* Allow specification of transparency at each vertex of preceding 2323* triangle or cylinder. Overrides any MATERIAL properties. 2324* 2325 ELSEIF (INTYPE .EQ. VERTRANSP) THEN 2326 IF (JUSTCLIPPED) GOTO 46 2327 IPREV = N - 1 2328 IF (TYPE(IPREV).EQ.NORMS .OR. TYPE(IPREV).EQ.VERTEXRGB) 2329 & IPREV = IPREV - 1 2330 IF (TYPE(IPREV).EQ.NORMS .OR. TYPE(IPREV).EQ.VERTEXRGB) 2331 & IPREV = IPREV - 1 2332 CALL ASSERT (TYPE(IPREV).EQ.TRIANG .OR. TYPE(IPREV).EQ.CYLIND 2333 & .OR.TYPE(IPREV).EQ.SPHERE, 2334 & 'orphan vertex transparency') 2335 NVTRANS = NVTRANS + 1 2336 IF (iand(FLAG(IPREV),TRANSP).EQ.0) NTRANSP = NTRANSP + 1 2337 FLAG(IPREV) = ior( FLAG(IPREV), TRANSP ) 2338 FLAG(IPREV) = ior( FLAG(IPREV), VTRANS ) 2339 DETAIL(NDET+1) = BUF(1) 2340 DETAIL(NDET+2) = BUF(2) 2341 DETAIL(NDET+3) = BUF(3) 2342 NDET = NDET + KDET(INTYPE) 2343 IF (SHADOW) THEN 2344 MDET = MDET + SDET(INTYPE) 2345 ENDIF 2346* 2347* Material properties are saved after enforcing legality 2348* 2349 ELSEIF (INTYPE .EQ. MATERIAL) THEN 2350* Mark this object as current material 2351 MSTATE = MATERIAL 2352 NPROPM = NPROPM + 1 2353* 2354* Expand arrays MLIST, MPARITY 2355 IF ( NPROPM .GT. SIZE(MPARITY) ) THEN 2356 CALL GET_TRY(SIZE(MPARITY), NPROPM, TRY1, 1) 2357 DO 940 ITRY = 1,3 2358 if (TRY1(ITRY).LE.0.OR.TRY1(ITRY).GT.MAXMEM) GOTO 940 2359 ALLOCATE( TMP1D(TRY1(ITRY)), stat=ierr ) 2360 if (ierr .NE. 0) GOTO 940 2361 TMP1D = 0 2362 TMP1D = MLIST 2363 CALL MOVE_ALLOC(from=TMP1D, to=MLIST) 2364 ALLOCATE( TMP1D(TRY1(ITRY)), stat=ierr ) 2365 if (ierr .NE. 0) GOTO 940 2366 TMP1D = 0 2367 TMP1D = MPARITY 2368 CALL MOVE_ALLOC(from=TMP1D, to=MPARITY) 2369 if(TEST_ALLOC)write(NOISE,*)"Expand MAXMAT to ", 2370 & try1(ITRY) 2371 GOTO 942 2372940 CONTINUE 2373 ENDIF 2374942 CALL ASSERT(NPROPM.LE.SIZE(MPARITY), 2375 & 'too many materials - increase MAXMAT and recompile') 2376 MLIST(NPROPM) = N 2377* Clear any previous material properties 2378 FLAG(N) = NPROPM*65536 2379 MATCOL = .FALSE. 2380 CLIPPING = .FALSE. 2381 NBOUNDS = 0 2382 ORTEPLIKE= .FALSE. 2383 BPRGB(1) = -1 2384* Phong power defaults to value in header 2385 IF (BUF(1).LT.0) BUF(1) = IPHONG 2386 DETAIL(NDET+1) = BUF(1) 2387* Specular reflection component defaults to value in header 2388 IF (BUF(2).LT.0 .OR. BUF(2).GT.1) BUF(2) = SPECLR 2389 DETAIL(NDET+2) = BUF(2) 2390* Negative values for specular highlighting indicate default to object 2391 CALL ASSERT (BUF(3).LE.1., 'red > 1 in material') 2392 CALL ASSERT (BUF(4).LE.1., 'grn > 1 in material') 2393 CALL ASSERT (BUF(5).LE.1., 'blu > 1 in material') 2394 DETAIL(NDET+3) = BUF(3) 2395 DETAIL(NDET+4) = BUF(4) 2396 DETAIL(NDET+5) = BUF(5) 2397 CLRITY = BUF(6) 2398 CALL ASSERT (CLRITY.GE.0., 'clarity < 0 in material') 2399 CALL ASSERT (CLRITY.LE.1., 'clarity > 1 in material') 2400 DETAIL(NDET+6) = CLRITY 2401* Transparency processing is necessarily a compromise, and several 2402* possible approximations may be useful; allow a choice among them 2403 CLROPT = BUF(7) 2404 DETAIL(NDET+7) = BUF(7) 2405* The next one is used in conjunction with bounding planes 2406* Dec 2010: no, it's used to select which transparency algorithm is used 2407 TRNSPOPT = BUF(8) 2408 DETAIL(NDET+8) = TRNSPOPT 2409* One remaining field reserved for future expansion 2410 DETAIL(NDET+9) = BUF(9) 2411* Initialize clipping planes, only used if CLIPPING is set below 2412 DETAIL(NDET+16) = FRONTCLIP 2413 DETAIL(NDET+17) = BACKCLIP 2414* Additional properties may continue on extra lines 2415 IF (INT(BUF(10)).GT.0) THEN 2416 L = 1 2417 DO I = 1,INT(BUF(10)) 2418 READ (INPUT,'(A)',END=50) LINE 2419 DO J = 132, 1, -1 2420 IF (LINE(J:J).NE.' '.AND.LINE(J:J).NE.' ') L = J 2421 ENDDO 2422 IF (LINE(L:L+4).EQ.'SOLID') THEN 2423 MATCOL = .TRUE. 2424 READ (LINE(L+6:132),*,END=720) RGBMAT(1),RGBMAT(2),RGBMAT(3) 2425 ELSE IF (LINE(L:L+7).EQ.'BACKFACE') THEN 2426 FLAG(N) = ior(FLAG(N),INSIDE) 2427 READ (LINE(L+9:132),*,END=720) RED, GRN, BLU, PHONGM, SPECM 2428 MPHONG = PHONGM 2429 IF (PHONGM.LT.0) MPHONG = IPHONG 2430 IF (SPECM.LT.0.OR.SPECM.GT.1.) SPECM = SPECLR 2431 DETAIL(NDET+11) = RED 2432 DETAIL(NDET+12) = GRN 2433 DETAIL(NDET+13) = BLU 2434 DETAIL(NDET+14) = MPHONG 2435 DETAIL(NDET+15) = SPECM 2436 ELSE IF (LINE(L:L+8).EQ.'FRONTCLIP') THEN 2437 CLIPPING = .TRUE. 2438 ZCLIP = FRONTCLIP 2439 READ (LINE(L+10:132),*,END=720) ZCLIP 2440 ZCLIP = ZCLIP * SCALE / TMAT(4,4) 2441 DETAIL(NDET+16) = ZCLIP 2442 ELSE IF (LINE(L:L+7).EQ.'BACKCLIP') THEN 2443 CLIPPING = .TRUE. 2444 ZCLIP = BACKCLIP 2445 READ (LINE(L+9:132),*,END=720) ZCLIP 2446 ZCLIP = ZCLIP * SCALE / TMAT(4,4) 2447 DETAIL(NDET+17) = ZCLIP 2448 ELSE IF (LINE(L:L+9).EQ.'ORTEP_LIKE') THEN 2449 ORTEPLIKE = .TRUE. 2450 ELSE IF (LINE(L:L+13).EQ.'BOUNDING_COLOR') THEN 2451 READ (LINE(L+15:132),*,END=720) BPRGB(1),BPRGB(2),BPRGB(3) 2452 CALL CHKRGB(BPRGB(1),BPRGB(2),BPRGB(3), 2453 & 'Invalid bounding color') 2454 ELSE IF (LINE(L:L+13).EQ.'BOUNDING_PLANE') THEN 2455 NBOUNDS = NBOUNDS + 1 2456 nbplanes = nbplanes + 1 2457* 2458* Expand arrays DETAIL and SDTAIL as needed 2459 NEEDMEM = MAX(N+NBOUNDS, NDET+NBOUNDS*KDET(INTERNAL)) 2460 if ( NEEDMEM .GT. SIZE(DETAIL) ) THEN 2461 CALL GET_TRY(SIZE(DETAIL), NEEDMEM, TRY1, 1) 2462 DO 950 ITRY = 1,3 2463 IF (TRY1(ITRY) .LE. 0 .OR. TRY1(ITRY) .GT. MAXMEM) 2464 & GOTO 950 2465 ALLOCATE( TMP1DR(TRY1(ITRY)), stat=ierr ) 2466 IF (IERR .NE. 0) GOTO 950 2467 TMP1DR = 0. 2468 TMP1DR = DETAIL 2469 if(TEST_ALLOC)write(NOISE,*)"Expand DETAIL to ", 2470 & size(TMP1DR) 2471 CALL MOVE_ALLOC(from=TMP1DR, to=DETAIL) 2472 GOTO 952 2473950 CONTINUE 2474 ENDIF 2475952 CALL ASSERT(NEEDMEM.LE.SIZE(DETAIL), 2476 & 'BP Oops') 2477 IF (SHADOW) THEN 2478 NEEDMEM = MDET+NBOUNDS*SDET(INTERNAL) 2479 if ( NEEDMEM .GT. SIZE(DETAIL) ) THEN 2480 CALL GET_TRY(SIZE(DETAIL), NEEDMEM, TRY1, 1) 2481 DO 955 ITRY = 1,3 2482 IF (TRY1(ITRY).LE.0.OR.TRY1(ITRY).GT.MAXMEM) 2483 & GOTO 955 2484 ALLOCATE( TMP1DR(TRY1(ITRY)), stat=ierr ) 2485 IF (IERR .NE. 0) GOTO 955 2486 TMP1DR = 0. 2487 TMP1DR = SDTAIL 2488 if(TEST_ALLOC)write(NOISE,*)"Expand SDTAIL to", 2489 & try1(ITRY) 2490 CALL MOVE_ALLOC(from=TMP1DR, to=SDTAIL) 2491 GOTO 957 2492955 CONTINUE 2493 ENDIF 2494957 CALL ASSERT(NEEDMEM .LE.SIZE(SDTAIL),'BP Oops') 2495 ENDIF 2496 NB = N + NBOUNDS 2497 CALL ASSERT( NB.LE.SIZE(DETAIL), 'BP Oops') 2498c OK, we've established there's room to store this bound; 2499c Flag all properties belonging to the parent material 2500 TYPE(NB) = INTERNAL 2501 FLAG(NB) = ior(FLAG(N),PROPS) 2502 IF(CLRITY.GT.0.0)THEN 2503 FLAG(NB) = ior(FLAG(NB),TRANSP) 2504 IF (CLROPT.EQ.1) FLAG(NB) = ior(FLAG(NB),MOPT1) 2505 ENDIF 2506 NBDET = KDET(MATERIAL) + (NBOUNDS-1)*KDET(INTERNAL) 2507 LIST(NB) = NDET + NBDET 2508 IF (SHADOW) THEN 2509 NBSDT = SDET(MATERIAL) + (NBOUNDS-1)*SDET(INTERNAL) 2510 MIST(NB) = MDET + NBSDT 2511 ENDIF 2512c Read in details, transform, and save 2513 READ (LINE(L+15:132),*,END=720) BPTYPE,XB,YB,ZB,xn,yn,zn 2514* Transform bounding plane along with objects 2515 CALL TRANSF(XB,YB,ZB) 2516* Rotate but don't translate normal, including post-rotation 2517 xnb = xn*TMAT(1,1) + yn*TMAT(2,1) + zn*TMAT(3,1) 2518 ynb = xn*TMAT(1,2) + yn*TMAT(2,2) + zn*TMAT(3,2) 2519 znb = xn*TMAT(1,3) + yn*TMAT(2,3) + zn*TMAT(3,3) 2520 xn = RAFTER(1,1)*xnb + RAFTER(1,2)*ynb + RAFTER(1,3)*znb 2521 yn = RAFTER(2,1)*xnb + RAFTER(2,2)*ynb + RAFTER(2,3)*znb 2522 zn = RAFTER(3,1)*xnb + RAFTER(3,2)*ynb + RAFTER(3,3)*znb 2523 temp = sqrt(xn*xn + yn*yn + zn*zn) 2524 xn = xn/temp 2525 yn = yn/temp 2526 zn = zn/temp 2527 IF (ORTEPLIKE .AND. ZN.LT.0) THEN 2528 XN = -XN 2529 YN = -YN 2530 ZN = -ZN 2531 ENDIF 2532* Save data in same arrays as for regular objects 2533* ISUBTYPE currently only used to flag ORTEP_LIKE ellipsoids 2534* BPTYPE is loaded on input but not currently used for anything 2535 ISUBTYPE = -1 2536 IF (ORTEPLIKE) ISUBTYPE = ORTEP 2537 DETAIL(NDET + NBDET + 1) = ISUBTYPE 2538 DETAIL(NDET + NBDET + 2) = BPTYPE 2539 DETAIL(NDET + NBDET + 3) = XB * SCALE + XCENT 2540 DETAIL(NDET + NBDET + 4) = YB * SCALE + YCENT 2541 DETAIL(NDET + NBDET + 5) = ZB * SCALE 2542 DETAIL(NDET + NBDET + 6) = XN 2543 DETAIL(NDET + NBDET + 7) = YN 2544 DETAIL(NDET + NBDET + 8) = ZN 2545c Most of the time BPRGB(1) is -1 to signal no special color 2546 DETAIL(NDET + NBDET + 9) = BPRGB(1) 2547 DETAIL(NDET + NBDET + 10) = BPRGB(2) 2548 DETAIL(NDET + NBDET + 11) = BPRGB(3) 2549 IF (SHADOW) THEN 2550 XR = SROT(1,1)*XB+SROT(1,2)*YB+SROT(1,3)*ZB 2551 YR = SROT(2,1)*XB+SROT(2,2)*YB+SROT(2,3)*ZB 2552 ZR = SROT(3,1)*XB+SROT(3,2)*YB+SROT(3,3)*ZB 2553 BPNORM(1) = SROT(1,1)*XN+SROT(1,2)*YN+SROT(1,3)*ZN 2554 BPNORM(2) = SROT(2,1)*XN+SROT(2,2)*YN+SROT(2,3)*ZN 2555 BPNORM(3) = SROT(3,1)*XN+SROT(3,2)*YN+SROT(3,3)*ZN 2556 SDTAIL(MDET + NBSDT + 1) = ISUBTYPE 2557 SDTAIL(MDET + NBSDT + 2) = BPTYPE 2558 SDTAIL(MDET + NBSDT + 3) = XR * SCALE + SXCENT 2559 SDTAIL(MDET + NBSDT + 4) = YR * SCALE + SYCENT 2560 SDTAIL(MDET + NBSDT + 5) = ZR * SCALE 2561 SDTAIL(MDET + NBSDT + 6) = BPNORM(1) 2562 SDTAIL(MDET + NBSDT + 7) = BPNORM(2) 2563 SDTAIL(MDET + NBSDT + 8) = BPNORM(3) 2564 ENDIF 2565 ELSE IF (LINE(L:L+6).EQ.'BUMPMAP') THEN 2566 WRITE(NOISE,*) '>> Sorry, no bumpmaps (dont you wish!)' 2567 ELSE 2568 GOTO 720 2569 ENDIF 2570 GOTO 721 2571720 WRITE(NOISE,'(A,A)') 2572 & '>> Unrecognized or incomplete MATERIAL option ',LINE 2573721 CONTINUE 2574 ENDDO 2575 ENDIF 2576* Update array pointers for material object itself 2577 DETAIL(NDET+18) = NBOUNDS 2578 NDET = NDET + KDET(INTYPE) 2579 IF (SHADOW) MDET = MDET + SDET(INTYPE) 2580* Update array pointers to allow for bounding planes and any other 2581* objects inserted by MATERIAL processing 2582 NDET = NDET + NBOUNDS*IDET(INTERNAL) 2583 IF (SHADOW) MDET = MDET + NBOUNDS*KDET(INTERNAL) 2584 N = N + NBOUNDS 2585* 2586 ELSEIF (INTYPE.EQ.GLOWLIGHT) THEN 2587 NGLOWS = NGLOWS + 1 2588* 2589* Expand array GLOWLIST 2590 IF ( NGLOWS .GT. SIZE(GLOWLIST) ) THEN 2591 CALL GET_TRY(SIZE(GLOWLIST), NGLOWS, TRY1, 1) 2592 DO 960 ITRY = 1,3 2593 if (TRY1(ITRY).LE.0.OR.TRY1(ITRY).GT.MAXMEM) GOTO 960 2594 ALLOCATE( TMP1D(TRY1(ITRY)), stat=ierr ) 2595 if (ierr .NE. 0) GOTO 960 2596 TMP1D = GLOWLIST 2597 CALL MOVE_ALLOC(from=TMP1D, to=GLOWLIST) 2598 if(TEST_ALLOC)write(NOISE,*)"Expand GLOWLIST to ", 2599 & try1(ITRY) 2600 GOTO 962 2601960 CONTINUE 2602 ENDIF 2603962 CALL ASSERT(NGLOWS.LE.SIZE(GLOWLIST),'too many glow lights') 2604 GLOWLIST(NGLOWS) = N 2605 GLOWSRC(1) = BUF(1) 2606 GLOWSRC(2) = BUF(2) 2607 GLOWSRC(3) = BUF(3) 2608 GLOWRAD = BUF(4) 2609 GLOW = BUF(5) 2610 GOPT = BUF(6) 2611 GPHONG = BUF(7) 2612 RED = BUF(8) 2613 GRN = BUF(9) 2614 BLU = BUF(10) 2615 CALL ASSERT (GLOW.GE.0,'illegal glow value') 2616 CALL ASSERT (GLOW.LE.1,'illegal glow value') 2617 IF (GLOW.GT.GLOWMAX) GLOWMAX = GLOW 2618 CALL CHKRGB (RED,GRN,BLU,'invalid glow light color') 2619* Isolated objects not transformed by TMAT, but still subject to inversion 2620 IF (ISOLATION.GT.0) THEN 2621 CALL ISOLATE(GLOWSRC(1),GLOWSRC(2)) 2622 ELSE 2623* transform coordinates and radius of glow source 2624 CALL TRANSF(GLOWSRC(1),GLOWSRC(2),GLOWSRC(3)) 2625 GLOWRAD = GLOWRAD / TMAT(4,4) 2626 ENDIF 2627 IF (EYEPOS.GT.0) PFAC = PERSP( GLOWSRC(3) ) 2628* save for rendering 2629 DETAIL(NDET+1) = GLOWSRC(1) * PFAC * SCALE + XCENT 2630 DETAIL(NDET+2) = GLOWSRC(2) * PFAC * SCALE + YCENT 2631 DETAIL(NDET+3) = GLOWSRC(3) * PFAC * SCALE 2632 DETAIL(NDET+4) = GLOWRAD * PFAC * SCALE 2633 DETAIL(NDET+5) = GLOW 2634 DETAIL(NDET+6) = GOPT 2635 DETAIL(NDET+7) = GPHONG 2636 DETAIL(NDET+8) = RED 2637 DETAIL(NDET+9) = GRN 2638 DETAIL(NDET+10) = BLU 2639 NDET = NDET + KDET(INTYPE) 2640 IF (SHADOW) THEN 2641 MDET = MDET + SDET(INTYPE) 2642 ENDIF 2643 2644 ELSEIF (INTYPE.EQ.QUADRIC) THEN 2645 NQUADS = NQUADS + 1 2646 IF (MSTATE.EQ.MATERIAL) THEN 2647 FLAG(N) = ior(FLAG(N),PROPS) 2648 NPROPS = NPROPS + 1 2649 IF (CLRITY.GT.0) THEN 2650 FLAG(N) = ior(FLAG(N),TRANSP) 2651 IF (CLROPT.EQ.1) FLAG(N) = ior(FLAG(N),MOPT1) 2652 NTRANSP = NTRANSP + 1 2653 ISTRANS = 1 2654 ENDIF 2655 IF (CLIPPING) THEN 2656 FLAG(N) = ior(FLAG(N),CLIPPED) 2657 ENDIF 2658 IF (NBOUNDS.GT.0) FLAG(N) = ior(FLAG(N),BOUNDED) 2659 IF (MATCOL) THEN 2660 BUF(5) = RGBMAT(1) 2661 BUF(6) = RGBMAT(2) 2662 BUF(7) = RGBMAT(3) 2663 ENDIF 2664 ENDIF 2665* 2666 ISQUAD = QINP( BUF(1), DETAIL(NDET+1), SHADOW, SDTAIL(MDET+1) ) 2667* 2668 IF (.NOT. ISQUAD) GOTO 45 2669 NDET = NDET + KDET(INTYPE) 2670 IF (SHADOW) THEN 2671 MDET = MDET + SDET(INTYPE) 2672 ENDIF 2673 2674C 2675C New object types go here! 2676C 2677 2678 ELSEIF (INTYPE.EQ.INTERNAL) THEN 2679 CALL ASSERT(.FALSE.,'object type 4 not available') 2680 ELSE 2681 CALL ASSERT(.FALSE.,'crash 50') 2682 ENDIF 2683 GO TO 7 2684c 2685c here to discard object due to clipping planes 2686c 268745 CONTINUE 2688 NCLIP = NCLIP + 1 268946 FLAG(N) = 0 2690 N = N - 1 2691 GO TO 7 2692 2693c 2694c 26-Aug-1999 attempt error recovery and reporting 2695c if input line is not recognized 269647 continue 2697 write (noise,'(A,A)') 'Unrecognized line: ',LINE 2698 goto 7 2699 2700* 2701* here for end of objects 2702* 270350 CONTINUE 2704 IF (INPUT.GT.INPUT0) THEN 2705 IF (VERBOSE) WRITE (NOISE,*) ' - closing indirect input file' 2706 CLOSE(INPUT) 2707 INPUT = INPUT - 1 2708 GOTO 7 2709 ENDIF 2710* 2711* help people re-center objects 2712* 2713 XA = (TRULIM(1,1) + TRULIM(1,2)) / 2.0 2714 YA = (TRULIM(2,1) + TRULIM(2,2)) / 2.0 2715 ZA = (TRULIM(3,1) + TRULIM(3,2)) / 2.0 2716 CALL TRANSF( XA, YA, ZA ) 2717 XA = TMAT(4,1) - XA * TMAT(4,4) 2718 YA = TMAT(4,2) - YA * TMAT(4,4) 2719 ZA = TMAT(4,3) - ZA * TMAT(4,4) 2720 IF (INVERT) YA = -YA 2721 WRITE (NOISE,*) 'To center objects in rendered scene, ', 2722 & 'change translation to:' 2723 WRITE (NOISE,'(3F10.4)') XA, YA, ZA 2724* 2725* Now we can set depth-cueing 2726 WRITE (LINE,577) 'Z limits (unscaled) before clipping:', 2727 & ZLIM(1)*TMAT(4,4)/SCALE,ZLIM(2)*TMAT(4,4)/SCALE 2728 WRITE (NOISE,*) LINE(2:57) 2729 WRITE (LINE,577) 'Z-clipping limits: ', 2730 & BACKCLIP*TMAT(4,4)/SCALE 2731 IF (FRONTCLIP.EQ.HUGE) THEN 2732 WRITE (LINE(48:57),'(A10)') ' none' 2733 ELSE 2734 WRITE (LINE(48:57),'(F10.2)') FRONTCLIP*TMAT(4,4)/SCALE 2735 ENDIF 2736 WRITE (NOISE,*) LINE(2:57) 2737 IF (VERBOSE) WRITE (NOISE,*) 'Scale: ', SCALE 2738 IF (VERBOSE.AND.GAMMACORRECTION) WRITE (NOISE,*) 'Gamma: ', Gamma 2739 IF (FOGTYPE .GE. 0) THEN 2740 IF (FOGBACK .EQ. 0) THEN 2741 FOGLIM(1) = ZLIM(1) 2742 ELSE 2743 FOGLIM(1) = FOGBACK * BACKCLIP 2744 ENDIF 2745 IF (FOGFRONT .EQ. 0) THEN 2746 FOGLIM(2) = ZLIM(2) 2747 ELSE IF (FRONTCLIP .LT. HUGE) THEN 2748 FOGLIM(2) = FOGFRONT * FRONTCLIP 2749 ELSE 2750 FOGLIM(2) = FOGFRONT * SCALE 2751 ENDIF 2752 IF (FOGTYPE.EQ.1) THEN 2753 WRITE(LINE,577) 'Fog (exponential) limits, density:' 2754 ELSE 2755 WRITE(LINE,577) 'Fog (linear) limits, density: ' 2756 ENDIF 2757 WRITE (LINE(38:67),578) 2758 & FOGLIM(1)*TMAT(4,4)/SCALE,FOGLIM(2)*TMAT(4,4)/SCALE,FOGDEN 2759 WRITE (NOISE,*) LINE(2:67) 2760 ENDIF 2761 577 FORMAT(1X,A,2F10.2) 2762 578 FORMAT(2F10.2,F10.2) 2763* 2764* Check list for special objects 2765* Triangle types first (vanilla/ribbon/surface) 2766 NRIB = 0 2767 NSUR = 0 2768 NTRI = 0 2769 DO 55 I = 1, N-1 2770 IF (TYPE(I).EQ.TRIANG) THEN 2771 NTRI = NTRI + 1 2772* Allow IPHONG=0 to disable special processing of triangles 2773 IF (IPHONG.EQ.0) GOTO 54 2774* Check for surface triangle (explicit normals in next record) 2775CDEBUG "I+1" should be replaced by INEXT processing 2776 IF (TYPE(I+1).EQ.NORMS.AND.iand(FLAG(I+1),HIDDEN).EQ.0) THEN 2777 FLAG(I) = ior( FLAG(I), SURFACE ) 2778 GOTO 54 2779 END IF 2780 IF (I.EQ.1) GOTO 54 2781* Check for ribbon triangles, 2782* can't possibly be one unless surrounded by other triangles 2783 IPREV = I - 1 2784 INEXT = I + 1 2785 IF (TYPE(IPREV).NE.TRIANG .OR. TYPE(INEXT).NE.TRIANG) THEN 2786 FLAG(I) = iand( FLAG(I), NOT(RIBBON) ) 2787 GOTO 54 2788 END IF 2789* trailing vertex must match one in previous triangle 2790 J = LIST(IPREV) 2791 K = LIST(I) 2792 L = LIST(INEXT) 2793 DO II = 1, 3 2794 KK = K+II 2795 if ( abs(detail(kk)-detail(j+ii+3)).gt.RIBBONSLOP 2796 & .and. abs(detail(kk)-detail(j+ii+6)).gt.RIBBONSLOP) goto 54 2797C IF (DETAIL(KK).NE.DETAIL(J+II+3) 2798C & .AND. DETAIL(KK).NE.DETAIL(J+II+6)) GOTO 54 2799 END DO 2800* leading vertex must match one in following triangle 2801 DO II = 7, 9 2802 KK = K+II 2803 IF (abs(DETAIL(KK)-DETAIL(L+II-3)).gt.RIBBONSLOP 2804 & .AND. abs(DETAIL(KK)-DETAIL(L+II-6)).gt.RIBBONSLOP) GOTO 54 2805C IF (DETAIL(KK).NE.DETAIL(L+II-3) 2806C & .AND. DETAIL(KK).NE.DETAIL(L+II-6)) GOTO 54 2807 END DO 2808 FLAG(I) = ior(FLAG(I),RIBBON) 280954 CONTINUE 2810 IF (iand(FLAG(I),RIBBON).NE.0) NRIB = NRIB + 1 2811 IF (iand(FLAG(I),SURFACE).NE.0) NSUR = NSUR + 1 2812 END IF 281355 CONTINUE 2814 IF (TYPE(N).EQ.TRIANG) NTRI = NTRI + 1 281556 CONTINUE 2816 2817* Set GLOW to maximum requested by glow light sources and bump up 2818* ambient contribution to compensate for darkening applied later 2819 AMBIEN = AMBIEN * (1. + GLOWMAX) 2820* 2821 WRITE(NOISE,*)'-------------------------------' 2822 IF (NSPHERE.NE.0) WRITE(NOISE,57) 'spheres =',NSPHERE 2823 IF (NCYLIND.NE.0) WRITE(NOISE,57) 'cylinders =',NCYLIND 2824 NTRI = NTRI - (NRIB + NSUR) 2825 IF (NPLANES.NE.0) WRITE(NOISE,57) 'planes =',NPLANES 2826 IF (NTRI.NE.0) WRITE(NOISE,57) 'plain triangles =',NTRI 2827 IF (NRIB.NE.0) WRITE(NOISE,57) 'ribbon triangles =',NRIB 2828 IF (NSUR.NE.0) WRITE(NOISE,57) 'surface triangles =',NSUR 2829 IF (NQUADS.NE.0) WRITE(NOISE,57) 'quadric surfaces =',NQUADS 2830 IF (NPROPM.NE.0) WRITE(NOISE,57) 'special materials =',NPROPM 2831 IF (NCLIP.NE.0) WRITE(NOISE,57) 'Z-clipped objects =',NCLIP 2832 IF (NTRANSP.NE.0) WRITE(NOISE,57) 'transparent objs =',NTRANSP 2833 IF (NHIDDEN.NE.0) WRITE(NOISE,57) 'hidden surfaces =',NHIDDEN 2834 IF (NINSIDE.NE.0) WRITE(NOISE,57) 'inside surfaces =',NINSIDE 2835 IF (nbplanes.NE.0) WRITE(NOISE,57)'bounding planes =',nbplanes 2836 WRITE(NOISE,57) 'total objects =',N 2837 WRITE(NOISE,*)'-------------------------------' 2838 IF (NGLOWS.GT.0) WRITE(NOISE,57) 'glow lights =',NGLOWS 2839 IF (LFLAG) THEN 2840 CALL LCLOSE( NLABELS ) 2841 WRITE(NOISE,57) 'labels =',NLABELS 2842 WRITE(NOISE,*)'-------------------------------' 2843 ELSEIF (NLABELS.NE.0) THEN 2844 WRITE(NOISE,57) 'labels (ignored) =',NLABELS 2845 WRITE(NOISE,*)'-------------------------------' 2846 ENDIF 284757 FORMAT(2X,A,I8) 2848* 2849 ZSLOP = SLOP * MAX(NPX,NPY) 2850 IF (VERBOSE) THEN 2851 WRITE (NOISE,*) 'ndet =',NDET,' SIZE(DETAIL)=',SIZE(DETAIL) 2852 IF (SHADOW) WRITE (NOISE,*) 'mdet =',MDET,' MAXSDT=', 2853 & SIZE(SDTAIL) 2854 WRITE (NOISE,*) 'EDGESLOP =',EDGESLOP 2855 WRITE (NOISE,*) ' ZSLOP =', ZSLOP 2856 ENDIF 2857* 2858* 2859* Sort objects, fill in "short lists" as indices into main list 2860* (note that it would lend itself better to "parallel 2861* processing" to form the short lists first and then 2862* sort each one - maybe even slightly more efficient in 2863* the present context!) 2864 DO 60 I = 1, N 2865 K = LIST(I) 2866 CALL ASSERT (K.GE.0,'k<0') 2867 CALL ASSERT (K.LT.NDET,'k>=ndet') 2868 IF (TYPE(I).EQ.TRIANG) THEN 2869 Z1 = DETAIL(K+3) 2870 Z2 = DETAIL(K+6) 2871 Z3 = DETAIL(K+9) 2872 ZTEMP(I) = MAX(Z1,Z2,Z3) 2873 ELSEIF (TYPE(I).EQ.SPHERE) THEN 2874 Z = DETAIL(K+3) 2875 R = DETAIL(K+4) 2876 ZTEMP(I) = Z + R 2877 ELSEIF (TYPE(I).EQ.CYLIND) THEN 2878* EAM May 1990 2879 Z1 = DETAIL(K+3) 2880 Z2 = DETAIL(K+7) 2881 R1 = DETAIL(K+4) 2882 R2 = DETAIL(K+8) 2883 ZTEMP(I) = MAX(Z1+R1,Z2+R2) 2884 ELSEIF (TYPE(I).EQ.PLANE 2885 & .OR. TYPE(I).EQ.NORMS 2886 & .OR. TYPE(I).EQ.MATERIAL 2887 & .OR. TYPE(I).EQ.VERTRANSP 2888 & .OR. TYPE(I).EQ.VERTEXRGB 2889 & .OR. TYPE(I).EQ.GLOWLIGHT 2890 & .OR. TYPE(I).EQ.INTERNAL) THEN 2891* EAM Mar 1994 (not sure this is necessary) 2892 ZTEMP(I) = SCALE + 1.0 2893 ELSEIF (TYPE(I).EQ.QUADRIC) THEN 2894 Z = DETAIL(K+3) 2895 R = DETAIL(K+4) 2896 ZTEMP(I) = Z + R 2897 ELSE 2898 CALL ASSERT(.FALSE.,'crash 60') 2899 ENDIF 290060 CONTINUE 2901 CALL HSORTD (N, ZTEMP, ZINDEX) 2902 KNTTOT = 0 2903 KNTMAX = 0 2904 DO J = 1, NTY 2905 DO I = 1, NTX 2906 KNTTOT = KNTTOT + KOUNT(I,J) 2907 IF (KOUNT(I,J).GT.KNTMAX) KNTMAX = KOUNT(I,J) 2908 ENDDO 2909 ENDDO 2910 IF (VERBOSE) THEN 2911 write (noise,*) 'max/avg length of short lists=', 2912 & kntmax,(knttot/(ntx*nty))+1 2913 WRITE (NOISE,*) 'knttot=',KNTTOT,' MAXSHR=',SIZE(KSHORT) 2914 ENDIF 2915* 2916* Expand array KSHORT 2917 IF ( KNTTOT .GT. SIZE(KSHORT) ) THEN 2918 CALL GET_TRY(SIZE(KSHORT), KNTTOT, TRY1, 1) 2919 DO 970 ITRY = 1,3 2920 if (TRY1(ITRY).LE.0.OR.TRY1(ITRY).GT.MAXMEM) GOTO 970 2921 ALLOCATE( TMP1D(TRY1(ITRY)), stat=ierr ) 2922 if (ierr .NE. 0) GOTO 970 2923 TMP1D = KSHORT 2924 CALL MOVE_ALLOC(from=TMP1D, to=KSHORT) 2925 if(TEST_ALLOC)write(NOISE,*)"Expand KSHORT to ",try1(ITRY) 2926 GOTO 972 2927970 CONTINUE 2928 ENDIF 2929972 CALL ASSERT (KNTTOT.LE.SIZE(KSHORT),'short list overflow') 2930 K = 0 2931 DO J = 1, NTY 2932 DO I = 1, NTX 2933 KSTART(I,J) = K+1 2934 KSTOP(I,J) = K 2935 K = K + KOUNT(I,J) 2936 ENDDO 2937 ENDDO 2938C write(NOISE,*)"Set kstart and kstop ",NTX," x ",NTY 2939 CALL ASSERT (K.EQ.KNTTOT,'k.ne.knttot') 2940 DO 90 I = 1, N 2941 IND = ZINDEX(N-I+1) 2942 CALL ASSERT (IND.GE.1,'ind<1') 2943 CALL ASSERT (IND.LE.N,'ind>n') 2944 K = LIST(IND) 2945 CALL ASSERT (K.GE.0,'k<0') 2946 CALL ASSERT (K.LT.NDET,'k>=ndet') 2947* impingement tests here must be same as above 2948 IF (TYPE(IND).EQ.TRIANG) THEN 2949 X1 = DETAIL(K+1) 2950 Y1 = DETAIL(K+2) 2951 Z1 = DETAIL(K+3) 2952 X2 = DETAIL(K+4) 2953 Y2 = DETAIL(K+5) 2954 Z2 = DETAIL(K+6) 2955 X3 = DETAIL(K+7) 2956 Y3 = DETAIL(K+8) 2957 Z3 = DETAIL(K+9) 2958 IXLO = MIN(X1,X2,X3)/NPX + 1 2959 IXHI = MAX(X1,X2,X3)/NPX + 1 2960 IYLO = MIN(Y1,Y2,Y3)/NPY + 1 2961 IYHI = MAX(Y1,Y2,Y3)/NPY + 1 2962 ELSEIF (TYPE(IND).EQ.SPHERE) THEN 2963 X = DETAIL(K+1) 2964 Y = DETAIL(K+2) 2965 Z = DETAIL(K+3) 2966 R = DETAIL(K+4) 2967 IXLO = (X-R)/NPX + 1 2968 IXHI = (X+R)/NPX + 1 2969 IYLO = (Y-R)/NPY + 1 2970 IYHI = (Y+R)/NPY + 1 2971 ELSEIF (TYPE(IND).EQ.CYLIND) THEN 2972 X1 = DETAIL(K+1) 2973 Y1 = DETAIL(K+2) 2974 Z1 = DETAIL(K+3) 2975 R1 = DETAIL(K+4) 2976 X2 = DETAIL(K+5) 2977 Y2 = DETAIL(K+6) 2978 Z2 = DETAIL(K+7) 2979 R2 = DETAIL(K+8) 2980 IXLO = MIN(X1-R1,X2-R2) / NPX + 1 2981 IXHI = MAX(X1+R1,X2+R2) / NPX + 1 2982 IYLO = MIN(Y1-R1,Y2-R2) / NPY + 1 2983 IYHI = MAX(Y1+R1,Y2+R2) / NPY + 1 2984 ELSEIF (TYPE(IND).EQ.PLANE) THEN 2985 IXLO = 1 2986 IXHI = NTX 2987 IYLO = 1 2988 IYHI = NTY 2989 ELSEIF (TYPE(IND).EQ.NORMS) THEN 2990 GOTO 81 2991 ELSEIF (TYPE(IND).EQ.VERTEXRGB) THEN 2992 GOTO 81 2993 ELSEIF (TYPE(IND).EQ.VERTRANSP) THEN 2994 GOTO 81 2995 ELSEIF (TYPE(IND).EQ.MATERIAL) THEN 2996 GOTO 81 2997 ELSEIF (TYPE(IND).EQ.GLOWLIGHT) THEN 2998 GOTO 81 2999 ELSEIF (TYPE(IND).EQ.QUADRIC) THEN 3000 X = DETAIL(K+1) 3001 Y = DETAIL(K+2) 3002 Z = DETAIL(K+3) 3003 R = DETAIL(K+4) 3004 IXLO = (X-R)/NPX + 1 3005 IXHI = (X+R)/NPX + 1 3006 IYLO = (Y-R)/NPY + 1 3007 IYHI = (Y+R)/NPY + 1 3008 ELSEIF (TYPE(IND).EQ.INTERNAL) THEN 3009 GOTO 81 3010 ELSE 3011 CALL ASSERT(.FALSE.,'crash 80') 3012 ENDIF 3013 IF (IXLO.LT.1 ) IXLO=1 3014 IF (IXLO.GT.NTX) GO TO 81 3015 IF (IXHI.LT.1 ) GO TO 81 3016 IF (IXHI.GT.NTX) IXHI=NTX 3017 IF (IYLO.LT.1 ) IYLO=1 3018 IF (IYLO.GT.NTY) GO TO 81 3019 IF (IYHI.LT.1 ) GO TO 81 3020 IF (IYHI.GT.NTY) IYHI=NTY 3021 DO 80 IY=IYLO,IYHI 3022 DO 80 IX=IXLO,IXHI 3023 KSTOP(IX,IY) = KSTOP(IX,IY) + 1 3024 KSHORT(KSTOP(IX,IY)) = IND 302580 CONTINUE 3026C write(NOISE,*)"Incr kstop ",IXLO,"-",IXHI," x ",IYLO,"-",IYHI, 3027C & kstart(IXLO,IYLO),kstop(IXLO,IYLO),kount(IXLO,IYLO) 302881 CONTINUE 302990 CONTINUE 3030 DO 95 J=1,NTY 3031 DO 95 I=1,NTX 3032 K1 = KSTART(I,J) 3033 K2 = KSTOP(I,J) 3034 K3 = KOUNT(I,J) 3035 if(K2-k1.ne.k3-1)then 3036 write(NOISE,*)"*** ERROR IN KOUNT,KSTART,KSTOP" 3037 write(NOISE,*)"I,J=",I,J," start,stop=",K1,K2," kount=",k3 3038 write(NOISE,*)"NTX,Y ",NTX,NTY," kount=",size(kount,1)," x", 3039 & size(kount,2)," stop=",size(kstop,1),size(kstop,2) 3040C write(NOISE,*)((kstart(III,JJJ),III=1,I),JJJ=1,J) 3041C write(NOISE,*)((kstop(III,JJJ),III=1,I),JJJ=1,J) 3042C write(NOISE,*)((kount(III,JJJ),III=1,I),JJJ=1,J) 3043 endif 3044 CALL ASSERT (K2-K1.EQ.K3-1,'k2-k1.ne.kount(i,j)-1') 3045 CALL ASSERT (K1.GE.1.AND.K1.LE.KNTTOT+1,'kstart(i,j)') 3046 CALL ASSERT (K2.GE.0.AND.K2.LE.KNTTOT,'kstop(i,j)') 304795 CONTINUE 3048* 3049* Do the short list business for shadow space too if required 3050 IF (SHADOW) THEN 3051 DO 160 I = 1, N 3052 K = MIST(I) 3053 CALL ASSERT (K.GE.0,'k.lt.0') 3054 CALL ASSERT (K.LT.MDET,'k.ge.mdet') 3055 IF (TYPE(I).EQ.TRIANG) THEN 3056 Z1 = SDTAIL(K+3) 3057 Z2 = SDTAIL(K+6) 3058 Z3 = SDTAIL(K+9) 3059 ZTEMP(I) = MAX(Z1,Z2,Z3) 3060 ELSEIF (TYPE(I).EQ.SPHERE) THEN 3061 Z = SDTAIL(K+3) 3062 R = SDTAIL(K+4) 3063 ZTEMP(I) = Z + R 3064 ELSEIF (TYPE(I).EQ.CYLIND) THEN 3065 Z1 = SDTAIL(K+3) 3066 Z2 = SDTAIL(K+7) 3067 R1 = SDTAIL(K+4) 3068 R2 = SDTAIL(K+8) 3069 ZTEMP(I) = MAX(Z1+R1,Z2+R2) 3070 ELSEIF (TYPE(I).EQ.PLANE) THEN 3071* no shadows for plane surface 3072 ELSEIF (TYPE(I).EQ.NORMS) THEN 3073* and certainly not for normals 3074 ELSEIF (TYPE(I).EQ.VERTEXRGB) THEN 3075 ELSEIF (TYPE(I).EQ.VERTRANSP) THEN 3076 ELSEIF (TYPE(I).EQ.MATERIAL) THEN 3077* or surface properties 3078 ELSEIF (TYPE(I).EQ.GLOWLIGHT) THEN 3079* you want a shadow on a light source??? 3080 ELSEIF (TYPE(I).EQ.QUADRIC) THEN 3081 Z = SDTAIL(K+3) 3082 R = SDTAIL(K+4) 3083 ZTEMP(I) = Z + R 3084 ELSEIF (TYPE(I).EQ.INTERNAL) THEN 3085 ELSE 3086 CALL ASSERT(.FALSE.,'crash 160') 3087 ENDIF 3088160 CONTINUE 3089 CALL HSORTD (N, ZTEMP, ZINDEX) 3090 MNTTOT = 0 3091 DO 170 J = 1, NSY 3092 DO 170 I = 1, NSX 3093 MNTTOT = MNTTOT + MOUNT(I,J) 3094170 CONTINUE 3095 IF (VERBOSE) WRITE (NOISE,*) 'mnttot=',MNTTOT,' MAXSSL=', 3096 & SIZE(MSHORT) 3097* 3098* Expand array MSHORT 3099 IF ( MNTTOT .GT. SIZE(MSHORT) ) THEN 3100 CALL GET_TRY(SIZE(MSHORT), MNTTOT, TRY1, 1) 3101 DO 975 ITRY = 1,3 3102 if (TRY1(ITRY).LE.0.OR.TRY1(ITRY).GT.MAXMEM) GOTO 975 3103 ALLOCATE( TMP1D(TRY1(ITRY)), stat=ierr ) 3104 if (ierr .NE. 0) GOTO 975 3105 TMP1D = MSHORT 3106 CALL MOVE_ALLOC(from=TMP1D, to=MSHORT) 3107 if(TEST_ALLOC)write(NOISE,*)"Expand MSHORT to ",try1(ITRY) 3108 GOTO 977 3109975 CONTINUE 3110 ENDIF 3111977 CALL ASSERT (MNTTOT.LE.SIZE(MSHORT), 3112 & 'shadow short list overflow') 3113 K = 0 3114 DO 175 J = 1, NSY 3115 DO 175 I = 1, NSX 3116 MSTART(I,J) = K+1 3117 MSTOP(I,J) = K 3118 K = K + MOUNT(I,J) 3119175 CONTINUE 3120 CALL ASSERT (K.EQ.MNTTOT,'k.ne.mnttot') 3121 DO 190 I = 1, N 3122 IND = ZINDEX(N-I+1) 3123 CALL ASSERT (IND.GE.1,'ind.lt.1') 3124 CALL ASSERT (IND.LE.N,'ind.gt.n') 3125 K = MIST(IND) 3126 CALL ASSERT (K.GE.0,'k.lt.0') 3127 CALL ASSERT (K.LT.MDET,'k.ge.mdet') 3128* impingement tests here must be same as above 3129 IF (TYPE(IND).EQ.TRIANG) THEN 3130 X1 = SDTAIL(K+1) 3131 Y1 = SDTAIL(K+2) 3132 Z1 = SDTAIL(K+3) 3133 X2 = SDTAIL(K+4) 3134 Y2 = SDTAIL(K+5) 3135 Z2 = SDTAIL(K+6) 3136 X3 = SDTAIL(K+7) 3137 Y3 = SDTAIL(K+8) 3138 Z3 = SDTAIL(K+9) 3139 IXLO = MIN(X1,X2,X3)/NPX + 1 3140 IXHI = MAX(X1,X2,X3)/NPX + 1 3141 IYLO = MIN(Y1,Y2,Y3)/NPY + 1 3142 IYHI = MAX(Y1,Y2,Y3)/NPY + 1 3143 ELSEIF (TYPE(IND).EQ.SPHERE) THEN 3144 X = SDTAIL(K+1) 3145 Y = SDTAIL(K+2) 3146 Z = SDTAIL(K+3) 3147 R = SDTAIL(K+4) 3148 IXLO = (X-R)/NPX + 1 3149 IXHI = (X+R)/NPX + 1 3150 IYLO = (Y-R)/NPY + 1 3151 IYHI = (Y+R)/NPY + 1 3152 ELSEIF (TYPE(IND).EQ.CYLIND) THEN 3153 X1 = SDTAIL(K+1) 3154 Y1 = SDTAIL(K+2) 3155 Z1 = SDTAIL(K+3) 3156 R1 = SDTAIL(K+4) 3157 X2 = SDTAIL(K+5) 3158 Y2 = SDTAIL(K+6) 3159 Z2 = SDTAIL(K+7) 3160 R2 = SDTAIL(K+8) 3161 IXLO = MIN(X1-R1,X2-R2) / NPX + 1 3162 IXHI = MAX(X1+R1,X2+R2) / NPX + 1 3163 IYLO = MIN(Y1-R1,Y2-R2) / NPY + 1 3164 IYHI = MAX(Y1+R1,Y2+R2) / NPY + 1 3165 ELSEIF (TYPE(IND).EQ.PLANE) THEN 3166* no shadows for plane surface 3167 GOTO 181 3168 ELSEIF (TYPE(IND).EQ.NORMS) THEN 3169 GOTO 181 3170 ELSEIF (TYPE(IND).EQ.VERTEXRGB) THEN 3171 GOTO 181 3172 ELSEIF (TYPE(IND).EQ.VERTRANSP) THEN 3173 GOTO 181 3174 ELSEIF (TYPE(IND).EQ.MATERIAL) THEN 3175 GOTO 181 3176 ELSEIF (TYPE(IND).EQ.GLOWLIGHT) THEN 3177 GOTO 181 3178 ELSEIF (TYPE(IND).EQ.QUADRIC) THEN 3179 X = SDTAIL(K+1) 3180 Y = SDTAIL(K+2) 3181 Z = SDTAIL(K+3) 3182 R = SDTAIL(K+4) 3183 IXLO = (X-R)/NPX + 1 3184 IXHI = (X+R)/NPX + 1 3185 IYLO = (Y-R)/NPY + 1 3186 IYHI = (Y+R)/NPY + 1 3187 ELSEIF (TYPE(IND).EQ.INTERNAL) THEN 3188 GOTO 181 3189 ELSE 3190 CALL ASSERT(.FALSE.,'crash 180') 3191 ENDIF 3192 IF (IXLO.LT.1 ) IXLO=1 3193 IF (IXLO.GT.NSX) GO TO 181 3194 IF (IXHI.LT.1 ) GO TO 181 3195 IF (IXHI.GT.NSX) IXHI=NSX 3196 IF (IYLO.LT.1 ) IYLO=1 3197 IF (IYLO.GT.NSY) GO TO 181 3198 IF (IYHI.LT.1 ) GO TO 181 3199 IF (IYHI.GT.NSY) IYHI=NSY 3200 DO 180 IY=IYLO,IYHI 3201 DO 180 IX=IXLO,IXHI 3202 MSTOP(IX,IY) = MSTOP(IX,IY) + 1 3203 MSHORT(MSTOP(IX,IY)) = IND 3204180 CONTINUE 3205181 CONTINUE 3206190 CONTINUE 3207 DO 195 J=1,NSY 3208 DO 195 I=1,NSX 3209 K1 = MSTART(I,J) 3210 K2 = MSTOP(I,J) 3211 K3 = MOUNT(I,J) 3212 CALL ASSERT (K2-K1.EQ.K3-1,'k2-k1.ne.mount(i,j)-1') 3213 CALL ASSERT (K1.GE.1.AND.K1.LE.MNTTOT+1,'mstart(i,j)') 3214 CALL ASSERT (K2.GE.0.AND.K2.LE.MNTTOT,'mstop(i,j)') 3215195 CONTINUE 3216 ENDIF 3217* 3218* Paint the tiles one by one 3219 DO 600 JTILE = 1, NTY 3220 DO 500 ITILE = 1, NTX 3221* bounds of this tile in pixel space 3222 XLO = (ITILE-1)*NPX 3223 XHI = ITILE*NPX 3224 YLO = (JTILE-1)*NPY 3225 YHI = JTILE*NPY 3226* initialize tile to background colour 3227 DO 200 J = 1, NPY 3228 DO 200 I = 1, NPX 3229 DO 199 IC = 1, 3 3230 TILE(IC,I,J) = BKGND(IC) 3231199 CONTINUE 3232 ACHAN(I,J) = 0.0 3233200 CONTINUE 3234* test for no objects in tile 3235 IF (KOUNT(ITILE,JTILE).EQ.0) GO TO 400 3236 NTRANSP = TTRANS(ITILE,JTILE) + nvtrans 3237 IJSTART = KSTART(ITILE,JTILE) 3238 IJSTOP = KSTOP(ITILE,JTILE) 3239* process non-empty tile 3240 DO 300 J = 1, NPY 3241 DO 300 I = 1, NPX 3242* location of the pixel in pixel space 3243 XP = XLO + 0.5 + (I-1) 3244 YP = YLO + 0.5 + (J-1) 3245* starting value of "highest z so far" 3246 ZTOP = BACKCLIP 3247* the index of the object that has it 3248 INDTOP = 0 3249* index of highest opaque object 3250 ZHIGH = ZTOP 3251* and number of transparent objects above it 3252 INDEPTH = 0 3253C Clear parity counter for all BOUNDED materials 3254 IF (NBPLANES.GT.0) THEN 3255 DO M = 1, NPROPM 3256 MPARITY(M) = 0 3257 ENDDO 3258 ENDIF 3259* find the highest pixel, using the tile's sorted list 3260C DO 240 IK = KSTART(ITILE,JTILE), KSTOP(ITILE,JTILE) 3261 DO 240 IK = IJSTART, IJSTOP 3262 IND = KSHORT(IK) 3263 K = LIST(IND) 3264 IFLAG = FLAG(IND) 3265* skip if hidden surface 3266 IF (NHIDDEN.GT.0 .AND. iand(IFLAG,HIDDEN).NE.0) goto 240 3267* further tests depend on object type 3268 IF (TYPE(IND).EQ.TRIANG) THEN 3269 X1 = DETAIL(K+1) 3270 Y1 = DETAIL(K+2) 3271 Z1 = DETAIL(K+3) 3272 X2 = DETAIL(K+4) 3273 Y2 = DETAIL(K+5) 3274 Z2 = DETAIL(K+6) 3275 X3 = DETAIL(K+7) 3276 Y3 = DETAIL(K+8) 3277 Z3 = DETAIL(K+9) 3278* cheap check for done pixel 3279 IF (Z1.LT.ZHIGH .AND. Z2.LT.ZHIGH .AND. Z3.LT.ZHIGH) 3280 & GOTO 250 3281 A = DETAIL(K+10) 3282 B = DETAIL(K+11) 3283 C = DETAIL(K+12) 3284 D = DETAIL(K+13) 3285* skip object if degenerate triangle 3286 IF (D.EQ.0) GO TO 240 3287* skip object if z not a new high 3288 ZP = A*XP + B*YP + C 3289 IF (ZP.LE.ZHIGH) GO TO 240 3290* Rigorous test to see if this point is interior to triangle 3291* NOTE: when lots of triangles are present, the following 3 lines 3292* account for the largest single chunk of rendering time (>10%)! 3293 S = (X2-X1)*(YP-Y1) - (Y2-Y1)*(XP-X1) 3294 T = (X3-X2)*(YP-Y2) - (Y3-Y2)*(XP-X2) 3295 U = (X1-X3)*(YP-Y3) - (Y1-Y3)*(XP-X3) 3296 IF ( (S.LT.0. .OR. T.LT.0. .OR. U.LT.0.) .AND. 3297 & (S.GT.0. .OR. T.GT.0. .OR. U.GT.0.) ) GO TO 240 3298* Z-clipped triangles are easy 3299 IF (iand(IFLAG,CLIPPED).NE.0) THEN 3300 MIND = LIST(MLIST(IFLAG/65536)) 3301 IF ( ZP.GT.DETAIL(MIND+16)) GOTO 240 3302 IF ( ZP.LT.DETAIL(MIND+17)) GOTO 240 3303 ENDIF 3304* Use Phong shading for surface and ribbon triangles. 3305 IF (iand(IFLAG,SURFACE+VCOLS+VTRANS).NE.0) THEN 3306 V = (Y3-Y1)*(X2-X1) - (X3-X1)*(Y2-Y1) 3307 W = (XP-X1)*(Y3-Y1) - (YP-Y1)*(X3-X1) 3308 ALPHA = W / V 3309 BETA = S / V 3310 ENDIF 3311 IF (iand(IFLAG,VCOLS+VTRANS).NE.0) THEN 3312 DETAIL(14 + LIST(IND)) = ALPHA 3313 DETAIL(15 + LIST(IND)) = BETA 3314 ENDIF 3315 IF (iand(IFLAG,SURFACE).NE.0) THEN 3316C CALL ASSERT(TYPE(IND+1).EQ.NORMS,'lost normals') 3317 A1 = DETAIL(1 + LIST(IND+1)) 3318 B1 = DETAIL(2 + LIST(IND+1)) 3319 C1 = DETAIL(3 + LIST(IND+1)) 3320 A2 = DETAIL(4 + LIST(IND+1)) 3321 B2 = DETAIL(5 + LIST(IND+1)) 3322 C2 = DETAIL(6 + LIST(IND+1)) 3323 A3 = DETAIL(7 + LIST(IND+1)) 3324 B3 = DETAIL(8 + LIST(IND+1)) 3325 C3 = DETAIL(9 + LIST(IND+1)) 3326 TEMPNORM(1) = A1 + ALPHA*(A2-A1) + BETA*(A3-A1) 3327 TEMPNORM(2) = B1 + ALPHA*(B2-B1) + BETA*(B3-B1) 3328 TEMPNORM(3) = C1 + ALPHA*(C2-C1) + BETA*(C3-C1) 3329* For ribbon triangles we take this normal for "middle" vertex, 3330* normal of previous triangle for "trailing" vertex normal, 3331* normal of next triangle for "leading" vertex normal. 3332* Then we use linear interpolation of vertex normals. 3333 ELSE IF (iand(IFLAG,RIBBON).NE.0) THEN 3334 IPREV = IND - 1 3335 INEXT = IND + 1 3336C CALL ASSERT(TYPE(IPREV).EQ.TRIANG,'lost triangle') 3337C CALL ASSERT(TYPE(INEXT).EQ.TRIANG,'lost triangle') 3338 V = (Y3-Y1)*(X2-X1) - (X3-X1)*(Y2-Y1) 3339 W = (XP-X1)*(Y3-Y1) - (YP-Y1)*(X3-X1) 3340 ALPHA = W / V 3341 BETA = S / V 3342 AT = DETAIL(10 + LIST(IPREV)) 3343 BT = DETAIL(11 + LIST(IPREV)) 3344 AL = DETAIL(10 + LIST(INEXT)) 3345 BL = DETAIL(11 + LIST(INEXT)) 3346 TEMPNORM(1) = -AT -ALPHA*(A-AT) -BETA*(AL-AT) 3347 TEMPNORM(2) = -BT -ALPHA*(B-BT) -BETA*(BL-BT) 3348 TEMPNORM(3) = 1. 3349 ELSE 3350 TEMPNORM(1) = -A 3351 TEMPNORM(2) = -B 3352 TEMPNORM(3) = 1. 3353 END IF 3354* Check bounding planes. 3355* This is different for triangles than for other shapes, as we assume 3356* that each triangle is only a facet of a larger shape that is really 3357* the 'object' being bounded. This means that rather than checking top 3358* and bottom surfaces of the current object, we have to search for 3359* them in other triangle/facets of the same bounded material. 3360 IF (iand(IFLAG,BOUNDED).NE.0) THEN 3361 MAT = IFLAG/65536 3362 M = MLIST(MAT)+1 3363 IF (.NOT.INBOUNDS( M, TYPE, LIST, DETAIL, 3364 & XP,YP,ZP,DX,DY,DZ,ZP, BPIND )) THEN 3365 GOTO 240 3366 ENDIF 3367c If this surface was above bounding plane, track parity 3368 MIND = LIST(MLIST(MAT)) 3369 IF (BPIND.NE.0) THEN 3370 IF (MPARITY(MAT).EQ.0) THEN 3371 MPARITY(MAT) = 1 3372 IF (ZP.LE.ZHIGH) GO TO 240 3373 TEMPNORM(1) = DX 3374 TEMPNORM(2) = DY 3375 TEMPNORM(3) = DZ 3376CORTEP Very ugly code to force bounding plane colors to be used 3377CORTEP but only if they are present. 3378 IF ( (BPIND.GT.0) 3379 & .AND.(DETAIL(LIST(BPIND)+9).GE.0.0)) 3380 & IND = BPIND 3381 ELSE 3382 MPARITY(MAT) = 0 3383 IF (FLAG(INDTOP)/65536.EQ.MAT) THEN 3384 IF (iand(IFLAG,TRANSP).EQ.0) THEN 3385 INDTOP = 0 3386 ZHIGH = BACKCLIP 3387 ZTOP = BACKCLIP 3388 INDEPTH= 0 3389 ELSE IF (INDEPTH.LE.1) THEN 3390 INDTOP = 0 3391 ZHIGH = BACKCLIP 3392 ZTOP = BACKCLIP 3393 INDEPTH= 0 3394 ELSE 3395 INDEPTH = INDEPTH - 1 3396 DO L = 1, INDEPTH 3397 INDLIST(L) = INDLIST(L+1) 3398 ZLIST(L) = ZLIST(L+1) 3399 NORMLIST(1,L) = NORMLIST(1,L+1) 3400 NORMLIST(2,L) = NORMLIST(2,L+1) 3401 NORMLIST(3,L) = NORMLIST(3,L+1) 3402 ENDDO 3403 ZTOP = ZLIST(1) 3404 ENDIF 3405 ENDIF 3406 GOTO 240 3407 ENDIF 3408 ENDIF 3409 ENDIF 3410* update values for object having highest z here yet 3411* 19-Feb-2002 Must wait til here to set NORMAL 3412 NORMAL(1) = TEMPNORM(1) 3413 NORMAL(2) = TEMPNORM(2) 3414 NORMAL(3) = TEMPNORM(3) 3415 IF (NTRANSP.GT.0) THEN 3416 CALL RANK( IND, ZP, NORMAL, FLAG ) 3417 ELSE 3418 ZHIGH = ZP 3419 INDTOP = IND 3420 ENDIF 3421 ELSEIF (TYPE(IND).EQ.SPHERE) THEN 3422 X = DETAIL(K+1) 3423 Y = DETAIL(K+2) 3424 Z = DETAIL(K+3) 3425 R = DETAIL(K+4) 3426* cheap check for done pixel 3427 IF (Z+R.LE.ZHIGH) GO TO 250 3428* skip object if point exterior 3429 DX = XP-X 3430 DY = YP-Y 3431 DX2 = DX**2 3432 DY2 = DY**2 3433 R2 = R**2 3434 IF (DX2+DY2 .GE. R2) GO TO 240 3435* skip object if z not a new high 3436 DZ = SQRT(R2-(DX2+DY2)) 3437C Triggered by CLROPT=2 3438 IF (iand(IFLAG,TRANSP).NE.0 .AND. 3439 & iand(IFLAG,INSIDE).NE.0) DZ = -DZ 3440 ZP = Z+DZ 3441 IF (ZP.LE.ZHIGH) GO TO 240 3442* Check bounding planes. 3443 IF (iand(IFLAG,BOUNDED).NE.0) THEN 3444 ZBACK = Z-DZ 3445 M = MLIST(IFLAG/65536)+1 3446 IF (.NOT.INBOUNDS( M, TYPE, LIST, DETAIL, 3447 & XP,YP,ZP,DX,DY,DZ,ZBACK, BPIND)) 3448 & GOTO 240 3449 ENDIF 3450* Z-clipped spheres aren't too bad 3451 IF (iand(IFLAG,CLIPPED).NE.0) THEN 3452 MIND = LIST(MLIST(IFLAG/65536)) 3453 IF (ZP.GT.DETAIL(MIND+16)) THEN 3454 ZP = Z-DZ 3455 IF (ZP.LE.ZHIGH) GOTO 240 3456 IF (ZP.GT.DETAIL(MIND+16)) GOTO 240 3457 DZ = -DZ 3458 ENDIF 3459 IF (ZP.LT.DETAIL(MIND+17)) GOTO 240 3460 ENDIF 3461* update values for object having highest z here yet 3462 NORMAL(1) = DX 3463 NORMAL(2) = DY 3464 NORMAL(3) = DZ 3465 IF (NTRANSP.GT.0) THEN 3466 CALL RANK( IND, ZP, NORMAL, FLAG ) 3467 ELSE 3468 ZHIGH = ZP 3469 INDTOP = IND 3470 ENDIF 3471 ELSEIF (TYPE(IND).EQ.CYLIND) THEN 3472* EAM May 1990 3473 X1 = DETAIL(K+1) 3474 Y1 = DETAIL(K+2) 3475 Z1 = DETAIL(K+3) 3476 R1 = DETAIL(K+4) 3477 X2 = DETAIL(K+5) 3478 Y2 = DETAIL(K+6) 3479 Z2 = DETAIL(K+7) 3480 R2 = R1 3481* EAM Mar 1993 with a better understanding of how this works 3482* add truly cheap test for cylinder entirely below current ZTOP 3483 TEMP1 = MAX(Z1+R1,Z2+R2) 3484 IF (TEMP1 .LE. ZHIGH) GOTO 250 3485* 2nd (relatively cheap) test 3486* is to check distance to cylinder axis in projection 3487 IF (X1.EQ.X2 .AND. Y1.EQ.Y2) THEN 3488 TEMP1 = 0.0 3489 ELSE 3490 TEMP1 = ((XP-X1)*(Y2-Y1) - (YP-Y1)*(X2-X1))**2 3491 & / ((X2-X1)*(X2-X1) + (Y2-Y1)*(Y2-Y1)) 3492 ENDIF 3493 IF (TEMP1 .GT. R1*R1) GOTO 240 3494* Now find Z coord in pixel space of point on surface of 3495* cylinder with these X and Y coords (ZP) 3496* Also get coords of closest point on cylinder axis (XYZA). 3497 ISCYL = CYLTEST( IFLAG, AXFRAC, 3498 & X1,Y1,Z1, X2,Y2,Z2, XP,YP,ZP, R1, XA,YA,ZA ) 3499 IF (.NOT.ISCYL) GO TO 240 3500* skip object if z not a new high 3501 IF (ZP.LE.ZHIGH) GO TO 240 3502 DX = XP - XA 3503 DY = YP - YA 3504 DZ = ZP - ZA 3505* Check bounding planes. Unfortunately we have to get the 3506* back surface first which means dummying up a call to CYLTEST 3507 IF (iand(IFLAG,BOUNDED).NE.0) THEN 3508 ZBACK = ZP 3509 ISCYL = CYLTEST( ior(IFLAG,INSIDE+TRANSP), AXFRAC, 3510 & X1,Y1,Z1, X2,Y2,Z2, XP,YP,ZBACK, R1, XA,YA,ZA) 3511 M = MLIST(IFLAG/65536)+1 3512 IF (.NOT.INBOUNDS( M, TYPE, LIST, DETAIL, 3513 & XP,YP,ZP,DX,DY,DZ,ZBACK, BPIND )) 3514 & GOTO 240 3515 ENDIF 3516* Z-clipped cylinders are messy 3517 IF (iand(IFLAG,CLIPPED).NE.0) THEN 3518 MIND = LIST(MLIST(IFLAG/65536)) 3519 IF (ZP.GT.DETAIL(MIND+16)) THEN 3520 ISCYL = CYLTEST( ior(IFLAG,INSIDE+TRANSP), AXFRAC, 3521 & X1,Y1,Z1, X2,Y2,Z2, XP,YP,ZP, R1, XA,YA,ZA ) 3522 IF (ZP.LE.ZHIGH) GOTO 240 3523 IF (ZP.GT.DETAIL(MIND+16)) GOTO 240 3524 ENDIF 3525 IF (ZP.LT.DETAIL(MIND+17)) GOTO 240 3526 DX = XP - XA 3527 DY = YP - YA 3528 DZ = ZP - ZA 3529 ENDIF 3530 NORMAL(1) = DX 3531 NORMAL(2) = DY 3532 NORMAL(3) = DZ 3533* if explicit vertex colors, need to keep fractional position 3534 IF (iand(IFLAG,ior(VCOLS,VTRANS)).NE.0) DETAIL(K+8) = AXFRAC 3535* update values for object having highest z here yet 3536 IF (NTRANSP.GT.0) THEN 3537 CALL RANK( IND, ZP, NORMAL, FLAG ) 3538 ELSE 3539 ZHIGH = ZP 3540 INDTOP = IND 3541 ENDIF 3542 ELSEIF (TYPE(IND).EQ.PLANE) THEN 3543 A = DETAIL(K+1) 3544 B = DETAIL(K+2) 3545 C = DETAIL(K+3) 3546 D = DETAIL(K+4) 3547 IF (D.EQ.0) GOTO 240 3548 ZP = A*XP + B*YP + C 3549 IF (ZP.LE.ZHIGH) GOTO 240 3550 NORMAL(1) = -A 3551 NORMAL(2) = -B 3552 NORMAL(3) = 1. 3553 IF (NTRANSP.GT.0) THEN 3554 CALL RANK( IND, ZP, NORMAL, FLAG ) 3555 ELSE 3556 ZHIGH = ZP 3557 INDTOP = IND 3558 ENDIF 3559 ELSEIF (TYPE(IND).EQ.QUADRIC) THEN 3560* First do cheap checks against projection of limiting sphere 3561 X = DETAIL(K+1) 3562 Y = DETAIL(K+2) 3563 Z = DETAIL(K+3) 3564 R = DETAIL(K+4) 3565 IF (Z+R.LE.ZHIGH) GO TO 250 3566 DX2 = (XP-X)**2 3567 DY2 = (YP-Y)**2 3568 R2 = R**2 3569 IF (DX2 + DY2 .GE. R2) GO TO 240 3570* Now find Z coord (ZP) in pixel space of point on quadric surface 3571* with these X and Y coords 3572 ISQUAD = QTEST( DETAIL(K+1), DETAIL(K+8), 3573 & XP, YP, ZP, QNORM, .FALSE., .FALSE. ) 3574 IF (.NOT.ISQUAD) GO TO 240 3575 IF (ZP.LE.ZHIGH) GO TO 240 3576* Check bounding planes. 3577 IF (iand(IFLAG,BOUNDED).NE.0) THEN 3578 M = MLIST(IFLAG/65536) 3579 IF (DETAIL(LIST(M+1)+1) .EQ. ORTEP) THEN 3580 CALL ORTEPBOUNDS( M+1, TYPE, LIST, DETAIL, XP,YP,ZP, 3581 & QNORM(1),QNORM(2),QNORM(3),ZBACK, BPIND) 3582CORTEP Very ugly code to force bounding plane colors to be used 3583CORTEP but only if they are present. 3584CORTEP An alternative would be to always set IND = BPIND, but 3585CORTEP check for presence of coloring info later, in which case 3586CORTEP IND itself needs to be temporarily saved somewhere. 3587CORTEP Or maybe just cache BPIND now and use it later if non-zero? 3588 IF ( (BPIND.GT.0) 3589 & .AND.(DETAIL(LIST(BPIND)+9).GE.0.0)) 3590 & IND = BPIND 3591 ELSE 3592 ISQUAD = QTEST( DETAIL(K+1), DETAIL(K+8), 3593 & XP, YP, ZBACK, QNORM, .FALSE., .TRUE. ) 3594 IF (.NOT.INBOUNDS( M+1, TYPE, LIST, DETAIL, XP,YP,ZP, 3595 & QNORM(1),QNORM(2),QNORM(3),ZBACK, BPIND)) 3596 & GOTO 240 3597 ENDIF 3598 ENDIF 3599C Z-clipping of quadric surfaces is encountered more frequently 3600C than for other object types, as the limiting sphere can also 3601C cause clipping. 3602C Check against limiting sphere in 3D 3603 MAYCLIP = .FALSE. 3604 DZ2 = (ZP-Z)**2 3605 IF (DX2+DY2+DZ2 .GT. R2) MAYCLIP = .TRUE. 3606 IF (iand(IFLAG,CLIPPED).NE.0) THEN 3607 MIND = LIST(MLIST(IFLAG/65536)) 3608 IF (ZP.GT.DETAIL(MIND+16)) MAYCLIP = .TRUE. 3609 ENDIF 3610 IF (MAYCLIP) THEN 3611 ISQUAD = QTEST( DETAIL(K+1), DETAIL(K+8), 3612 & XP, YP, ZP, QNORM, .FALSE., .TRUE. ) 3613 IF (.NOT.ISQUAD) GO TO 240 3614 IF (ZP.LE.ZHIGH) GO TO 240 3615 DZ2 = (ZP-Z)**2 3616 IF (DX2+DY2+DZ2 .GT. R2) GO TO 240 3617 IF (iand(IFLAG,CLIPPED).NE.0) THEN 3618 IF (ZP.GT.DETAIL(MIND+16)) GO TO 240 3619 IF (ZP.LT.DETAIL(MIND+17)) GO TO 240 3620 ENDIF 3621 ENDIF 3622 NORMAL(1) = QNORM(1) 3623 NORMAL(2) = QNORM(2) 3624 NORMAL(3) = QNORM(3) 3625* update values for object having highest z here yet 3626 IF (NTRANSP.GT.0) THEN 3627 CALL RANK( IND, ZP, NORMAL, FLAG ) 3628 ELSE 3629 ZHIGH = ZP 3630 INDTOP = IND 3631 ENDIF 3632 ELSE 3633 CALL ASSERT(.FALSE.,'crash 240') 3634 ENDIF 3635240 CONTINUE 3636250 CONTINUE 3637C Apply background fog here (added 2010) 3638 IF (FOGTYPE .GE. 0) THEN 3639 FOGDIM = FOGGY( FOGLIM(2) - ZTOP ) 3640 TILE(1,I,J) = (1.-FOGDIM)*TILE(1,I,J) + FOGDIM*FOGRGB(1) 3641 TILE(2,I,J) = (1.-FOGDIM)*TILE(2,I,J) + FOGDIM*FOGRGB(2) 3642 TILE(3,I,J) = (1.-FOGDIM)*TILE(3,I,J) + FOGDIM*FOGRGB(3) 3643 ENDIF 3644 3645* Background colour if we never found an object in this line of sight 3646 IF (INDTOP.EQ.0) GO TO 299 3647C We now know this is not a background pixel so set alpha channel to 1 3648C Modify later if it turns out the object is transparent 3649 ACHAN(I,J) = 1.0 3650C Transparency processing revamped Mar 2001 3651C If the top object is transparent we will have to come back here 3652C later and do this all again for each object in INDLIST 3653 IF (NTRANSP.NE.0) THEN 3654C CALL ASSERT(INDEPTH.GT.0,'INDEPTH = 0') 3655 ITPASS = 1 3656 ZHIGH = ZLIST(INDEPTH) 3657 INDTOP = INDLIST(INDEPTH) 3658 NORMAL(1) = NORMLIST(1,INDEPTH) 3659 NORMAL(2) = NORMLIST(2,INDEPTH) 3660 NORMAL(3) = NORMLIST(3,INDEPTH) 3661 RGBLND(1) = BKGND(1) 3662 RGBLND(2) = BKGND(2) 3663 RGBLND(3) = BKGND(3) 3664 ENDIF 3665* ZP is the "height" of the chosen pixel, 3666* and indtop tells us which object it came from: 3667 ZTOP = ZHIGH 3668 ZP = ZTOP 3669255 CONTINUE 3670C 3671C Shadowing code - look for objects that shadow the one we just found 3672C 3673 IF (SHADOW) THEN 3674* locate pixel in shadow space 3675* take out object translation & scaling 3676 XT = (XP - XCENT) / SCALE 3677 YT = (YP - YCENT) / SCALE 3678 ZT = ZP / SCALE 3679* rotate light source position to z axis 3680 XR = SROT(1,1)*XT+SROT(1,2)*YT+SROT(1,3)*ZT 3681 YR = SROT(2,1)*XT+SROT(2,2)*YT+SROT(2,3)*ZT 3682 ZR = SROT(3,1)*XT+SROT(3,2)*YT+SROT(3,3)*ZT 3683* scale and translate for shadow space 3684 XS = XR * SCALE + SXCENT 3685 YS = YR * SCALE + SYCENT 3686 ZS = ZR * SCALE 3687* determine appropriate shadow tile 3688 ISTILE = XS/NPX + 1 3689 JSTILE = YS/NPY + 1 3690* Just to get proper error message 3691 IF (JSTILE.LE.0) JSTILE = NSY + 1 - JSTILE 3692 IF (ISTILE.LE.0) ISTILE = NSX + 1 - ISTILE 3693 IF (JSTILE.GE.NSY) THEN 3694 NSYMAX = MAX(JSTILE,NSYMAX) 3695 INDSTP = 0 3696 GOTO 270 3697 END IF 3698 IF (ISTILE.GE.NSX) THEN 3699 NSXMAX = MAX(ISTILE,NSXMAX) 3700 INDSTP = 0 3701 GOTO 270 3702 END IF 3703* starting value of "highest shadow space z so far" 3704* and the index of the object that has it 3705 ZSTOP = 2.0*BACKCLIP 3706 INDSTP = 0 3707* 3708 DO 260 IK = MSTART(ISTILE,JSTILE), MSTOP(ISTILE,JSTILE) 3709 IND = MSHORT(IK) 3710 IFLAG = FLAG(IND) 3711* Ignore transparent objects except for the one being shaded 3712 IF (iand(IFLAG,TRANSP).NE.0 .AND. IND.NE.INDTOP) GOTO 260 3713 K = MIST(IND) 3714 IF (TYPE(IND).EQ.TRIANG) THEN 3715 X1 = SDTAIL(K+1) 3716 Y1 = SDTAIL(K+2) 3717 Z1 = SDTAIL(K+3) 3718 X2 = SDTAIL(K+4) 3719 Y2 = SDTAIL(K+5) 3720 Z2 = SDTAIL(K+6) 3721 X3 = SDTAIL(K+7) 3722 Y3 = SDTAIL(K+8) 3723 Z3 = SDTAIL(K+9) 3724 A = SDTAIL(K+10) 3725 B = SDTAIL(K+11) 3726 C = SDTAIL(K+12) 3727 D = SDTAIL(K+13) 3728* cheap check for done "pixel" 3729 IF (Z1.LT.ZSTOP .AND. Z2.LT.ZSTOP .AND. Z3.LT.ZSTOP) 3730 & GOTO 270 3731* skip object if degenerate triangle 3732 IF (D.EQ.0) GO TO 260 3733* skip object if z not a new high 3734 ZTEST = A*XS + B*YS + C 3735 IF (ZTEST.LE.ZSTOP) GO TO 260 3736* skip object if point exterior 3737 S = (X2-X1)*(YS-Y1)-(Y2-Y1)*(XS-X1) 3738 T = (X3-X2)*(YS-Y2)-(Y3-Y2)*(XS-X2) 3739 U = (X1-X3)*(YS-Y3)-(Y1-Y3)*(XS-X3) 3740 IF ( (S.LT.0. .OR. T.LT.0. .OR. U.LT.0.) .AND. 3741 & (S.GT.0. .OR. T.GT.0. .OR. U.GT.0.) ) GO TO 260 3742* Check bounding planes 3743 IF (iand(IFLAG,BOUNDED).NE.0) THEN 3744 MAT = IFLAG/65536 3745 M = MLIST(MAT)+1 3746 IF (.NOT.INBOUNDS( M, TYPE, MIST, SDTAIL, 3747 & XS,YS,ZTEST,DX,DY,DZ,ZTEST, BPIND )) THEN 3748 GOTO 260 3749 ENDIF 3750 MIND = MIST(MLIST(MAT)) 3751 IF (BPIND.NE.0) THEN 3752 IF (MPARITY(MAT).GE.0) THEN 3753 MPARITY(MAT) = -1 3754 ELSE 3755 MPARITY(MAT) = 0 3756 IF (FLAG(INDSTP)/65536.EQ.MAT) THEN 3757 INDSTP = 0 3758 ZSTOP = 2.*BACKCLIP 3759 ENDIF 3760 GOTO 260 3761 ENDIF 3762 ENDIF 3763 ENDIF 3764* update values for object having highest z here yet 3765 ZSTOP = ZTEST 3766 INDSTP = IND 3767 ELSEIF (TYPE(IND).EQ.SPHERE) THEN 3768 X = SDTAIL(K+1) 3769 Y = SDTAIL(K+2) 3770 Z = SDTAIL(K+3) 3771 R = SDTAIL(K+4) 3772* cheap check for done "pixel" 3773 IF (Z+R.LT.ZSTOP) GO TO 270 3774* skip object if point exterior 3775 DX = XS-X 3776 DY = YS-Y 3777 DX2 = DX**2 3778 DY2 = DY**2 3779 R2 = R**2 3780 IF (DX2+DY2 .GE. R2) GO TO 260 3781* skip object if z not a new high 3782 DZ = SQRT(R2-(DX2+DY2)) 3783 ZTEST = Z+DZ 3784 IF (ZTEST.LE.ZSTOP) GO TO 260 3785* Check bounding planes. 3786 IF (iand(IFLAG,BOUNDED).NE.0) THEN 3787 ZBACK = Z-DZ 3788 M = MLIST(IFLAG/65536)+1 3789 IF (.NOT.INBOUNDS( M, TYPE, MIST, SDTAIL, 3790 & XS,YS,ZTEST,DX,DY,DZ,ZBACK, BPIND)) 3791 & GOTO 260 3792 ENDIF 3793* update values for object having highest z here yet 3794 ZSTOP = ZTEST 3795 INDSTP = IND 3796 ELSEIF (TYPE(IND).EQ.CYLIND) THEN 3797* EAM May 1990 3798 X1 = SDTAIL(K+1) 3799 Y1 = SDTAIL(K+2) 3800 Z1 = SDTAIL(K+3) 3801 R1 = SDTAIL(K+4) 3802 X2 = SDTAIL(K+5) 3803 Y2 = SDTAIL(K+6) 3804 Z2 = SDTAIL(K+7) 3805 R2 = R1 3806* EAM Feb 93 - Test first to see if entire cylinder is below 3807* current top object in shadow space 3808 IF (MAX( Z1+R1, Z2+R2 ) .LT. ZSTOP) GOTO 270 3809* Now find Z coord (ZTEST) in pixel space of point on 3810* surface of cylinder with these X and Y coords 3811 ISCYL = CYLTEST( IFLAG, AXFRAC, 3812 & X1,Y1,Z1, X2,Y2,Z2, XS,YS,ZTEST, R1, XA,YA,ZA ) 3813 IF (.NOT.ISCYL) GO TO 260 3814* skip object if z not a new high 3815 IF (ZTEST.LE.ZSTOP) GO TO 260 3816* Check bounding planes. 3817 IF (iand(IFLAG,BOUNDED).NE.0) THEN 3818 ISCYL = CYLTEST( ior(IFLAG,INSIDE+TRANSP), AXFRAC, 3819 & X1,Y1,Z1, X2,Y2,Z2, XS,YS,ZBACK, R1, XA,YA,ZA) 3820 M = MLIST(IFLAG/65536)+1 3821 IF (.NOT.INBOUNDS( M, TYPE, MIST, SDTAIL, 3822 & XS,YS,ZTEST,DX,DY,DZ,ZBACK, BPIND )) 3823 & GOTO 260 3824 ENDIF 3825* update values for object having highest z here yet 3826 ZSTOP = ZTEST 3827 INDSTP = IND 3828 ELSEIF (TYPE(IND).EQ.QUADRIC) THEN 3829 X = SDTAIL(K+1) 3830 Y = SDTAIL(K+2) 3831 Z = SDTAIL(K+3) 3832 R = SDTAIL(K+4) 3833* cheap check against limiting sphere 3834 IF (Z+R.LT.ZSTOP) GO TO 270 3835 DX = XS-X 3836 DY = YS-Y 3837 R2 = R**2 3838 IF (DX**2 + DY**2 .GE. R2) GO TO 260 3839* Now find Z coord (ZTEST) in shadow pixel space of point on 3840* surface with these X and Y coords 3841 ISQUAD = QTEST( SDTAIL(K+1), SDTAIL(K+5), 3842 & XS, YS, ZTEST, QNORM, .TRUE., .FALSE. ) 3843CDEBUG XS, YS, ZTEST, QNORM, .TRUE., .TRUE. ) 3844CDEBUG 16-Dec-1998 I inverted the BACKSIDE = TRUE/FALSE flags from 3845CDEBUG what they "ought" to be to remove buggy shadows from a test 3846CDEBUG case parabolic hyperboloid. I don't understand why this would be 3847CDEBUG necessary, and worry a bit that it breaks something else. 3848CDEBUG 3849 IF (.NOT.ISQUAD) GO TO 260 3850* skip object if z not a new high 3851 IF (ZTEST.LE.ZSTOP) GO TO 260 3852* Check bounding planes. 3853 IF (iand(IFLAG,BOUNDED).NE.0) THEN 3854 M = MLIST(IFLAG/65536) 3855 IF (DETAIL(LIST(M)+1) .EQ. ORTEP) THEN 3856 ISQUAD = QTEST( SDTAIL(K+1), SDTAIL(K+5), 3857 & XS, YS, ZBACK, QNORM, .TRUE., .TRUE. ) 3858 IF (.NOT.INBOUNDS(M+1, 3859 & TYPE, MIST, SDTAIL, XS,YS,ZTEST, 3860 & QNORM(1),QNORM(2),QNORM(3),ZBACK, BPIND)) 3861 & GOTO 260 3862 ENDIF 3863 ENDIF 3864* Test against bounding sphere in 3D 3865* and if surface nearest to light is clipped, check back also 3866 DZ = ZTEST-Z 3867 IF (DX**2 + DY**2 + DZ**2 .GE. R2) THEN 3868 ISQUAD = QTEST( SDTAIL(K+1), SDTAIL(K+5), 3869 & XS, YS, ZTEST, QNORM, .TRUE., .TRUE. ) 3870CDEBUG XS, YS, ZTEST, QNORM, .TRUE., .FALSE. ) 3871 IF (.NOT.ISQUAD) GO TO 260 3872 IF (ZTEST.LE.ZSTOP) GO TO 260 3873 DZ = ZTEST - Z 3874 IF (DX**2 + DY**2 + DZ**2 .GE. R2) GO TO 260 3875 ENDIF 3876* update values for object having highest z here yet 3877 ZSTOP = ZTEST 3878 INDSTP = IND 3879C No more legal object types; should never happen 3880 ELSE 3881 CALL ASSERT(.FALSE.,'shadow tile error, crash 260') 3882 ENDIF 3883260 CONTINUE 3884270 CONTINUE 3885C End of search for objects that shadow this one 3886 if ((zstop+zslop).lt.zs .and. indstp.ne.0) nslow = nslow + 1 3887 ELSE 3888 ZS = 0. 3889 ZSTOP = 0. 3890 INDSTP = INDTOP 3891 ENDIF 3892* if roundoff made us miss the object, we are probably 3893* at a pixel that is very near the edge of the object 3894* from the point of view of the light source, so just 3895* treat it as if not in shadow 3896 IF (INDSTP.EQ.0) THEN 3897 ZS = 0. 3898 ZSTOP = 0. 3899 INDSTP = INDTOP 3900 ENDIF 3901* 3902* Pick up colours of object to be shaded 3903* 3904 K = LIST(INDTOP) 3905 IF (TYPE(INDTOP).EQ.TRIANG) THEN 3906 IF (iand(FLAG(INDTOP),VCOLS).NE.0) THEN 3907 ALPHA = DETAIL(K+14) 3908 BETA = DETAIL(K+15) 3909 INEXT = INDTOP + 1 3910 IF ((TYPE(INEXT).EQ.NORMS).OR.(TYPE(INEXT).EQ.VERTRANSP)) 3911 & INEXT = INEXT + 1 3912 IF ((TYPE(INEXT).EQ.NORMS).OR.(TYPE(INEXT).EQ.VERTRANSP)) 3913 & INEXT = INEXT + 1 3914 K = LIST(INEXT) 3915 CALL ASSERT(TYPE(INEXT).EQ.VERTEXRGB,'lost vertex colors') 3916 RGBCUR(1) = DETAIL(K+1) 3917 & + ALPHA*(DETAIL(K+4)-DETAIL(K+1)) 3918 & + BETA*(DETAIL(K+7)-DETAIL(K+1)) 3919 RGBCUR(2) = DETAIL(K+2) 3920 & + ALPHA*(DETAIL(K+5)-DETAIL(K+2)) 3921 & + BETA*(DETAIL(K+8)-DETAIL(K+2)) 3922 RGBCUR(3) = DETAIL(K+3) 3923 & + ALPHA*(DETAIL(K+6)-DETAIL(K+3)) 3924 & + BETA*(DETAIL(K+9)-DETAIL(K+3)) 3925 ELSE 3926 RGBCUR(1) = DETAIL(K+14) 3927 RGBCUR(2) = DETAIL(K+15) 3928 RGBCUR(3) = DETAIL(K+16) 3929 ENDIF 3930 ELSEIF (TYPE(INDTOP).EQ.SPHERE) THEN 3931 RGBCUR(1) = DETAIL(K+5) 3932 RGBCUR(2) = DETAIL(K+6) 3933 RGBCUR(3) = DETAIL(K+7) 3934 ELSEIF (TYPE(INDTOP).EQ.CYLIND) THEN 3935 IF (iand(FLAG(INDTOP),VCOLS).NE.0) THEN 3936 FRAC = DETAIL(K+8) 3937 INEXT = INDTOP + 1 3938 IF (TYPE(INEXT).EQ.VERTRANSP) INEXT = INEXT + 1 3939 K = LIST(INEXT) 3940 CALL ASSERT(TYPE(INEXT).EQ.VERTEXRGB,'lost vertex colors') 3941 RGBCUR(1) = FRAC*DETAIL(K+4) + (1.-FRAC)*DETAIL(K+1) 3942 RGBCUR(2) = FRAC*DETAIL(K+5) + (1.-FRAC)*DETAIL(K+2) 3943 RGBCUR(3) = FRAC*DETAIL(K+6) + (1.-FRAC)*DETAIL(K+3) 3944 ELSE 3945 RGBCUR(1) = DETAIL(K+9) 3946 RGBCUR(2) = DETAIL(K+10) 3947 RGBCUR(3) = DETAIL(K+11) 3948 ENDIF 3949* EAM Mar 1993 PLANE is shaded from full colour in foreground 3950* to half-saturation at horizon 3951 ELSEIF (TYPE(INDTOP).EQ.PLANE) THEN 3952 FADE = (ZP + 3.*SCALE) / (4.*SCALE) 3953 RGBCUR(1) = FADE * DETAIL(K+5) + (1.-FADE) * BKGND(1) 3954 RGBCUR(2) = FADE * DETAIL(K+6) + (1.-FADE) * BKGND(2) 3955 RGBCUR(3) = FADE * DETAIL(K+7) + (1.-FADE) * BKGND(3) 3956 ELSEIF (TYPE(INDTOP).EQ.QUADRIC) THEN 3957 RGBCUR(1) = DETAIL(K+5) 3958 RGBCUR(2) = DETAIL(K+6) 3959 RGBCUR(3) = DETAIL(K+7) 3960 ELSEIF (TYPE(INDTOP).EQ.INTERNAL) THEN 3961 RGBCUR(1) = DETAIL(K+9) 3962 RGBCUR(2) = DETAIL(K+10) 3963 RGBCUR(3) = DETAIL(K+11) 3964 ELSE 3965 WRITE(LINE,*) 'Top object claims to be type',TYPE(INDTOP) 3966 CALL ASSERT(.FALSE.,LINE) 3967 ENDIF 3968* 3969* Get shading components. 3970* 3971 3972* 3973* 11-May-1997 As of now, treat negative normal(3) as indicating the 3974* back side of a material. Default is to shrug and invert the normal. 3975* Some material have explicit BACKFACE proterties, however. 3976* 3977 BACKFACE = .FALSE. 3978 IF (NORMAL(3).LE.0) THEN 3979 NORMAL(1) = -NORMAL(1) 3980 NORMAL(2) = -NORMAL(2) 3981 NORMAL(3) = -NORMAL(3) 3982 BACKFACE = .TRUE. 3983 IF (iand(FLAG(INDTOP),PROPS).NE.0) THEN 3984 K = FLAG(INDTOP) / 65536 3985 if(K.le.0)WRITE(NOISE,*)"FLAG(",INDTOP,")=",FLAG(INDTOP) 3986 CALL ASSERT(K.GT.0,'lost material definition') 3987 IF (iand(FLAG(MLIST(K)),INSIDE).NE.0) THEN 3988 K = LIST(MLIST(K)) 3989 RGBCUR(1) = DETAIL(K+11) 3990 RGBCUR(2) = DETAIL(K+12) 3991 RGBCUR(3) = DETAIL(K+13) 3992 END IF 3993 END IF 3994 END IF 3995* 3996 ABSN = SQRT(NORMAL(1)**2 + NORMAL(2)**2 + NORMAL(3)**2) 3997c CALL ASSERT(ABS(ABSN-1.0).LT.0.02,'>> Abnormal normal') 3998 NL(1) = NORMAL(1) / ABSN 3999 NL(2) = NORMAL(2) / ABSN 4000 NL(3) = NORMAL(3) / ABSN 4001 SDIFF = NL(3) * STRAIT*DIFFUS 4002 SSP = 2.*NL(3)**2 - 1. 4003* We do the value check like this to avoid floating-point underflows 4004* in the Phonging. We also save calculation time this way, because 4005* the 0 case will occur often for reasonably high Phong powers. 4006* Note that PHOBND**IPHONG should evaluate to the cutoff value 4007* between significant and insignificant specular contribution. 4008* The contributions that are actually computed here 4009* can be smaller by a factor of STRAIT*SPECLR: 4010 IF (SSP.LT.PHOBND .OR. IPHONG.EQ.0) THEN 4011 SSPEC = 0. 4012 ELSE 4013 SSPEC = SSP**IPHONG * STRAIT*SPECLR 4014 ENDIF 4015 LDOTN = SOURCE(1)*NL(1)+SOURCE(2)*NL(2)+SOURCE(3)*NL(3) 4016 IF (LDOTN.LE.0.) THEN 4017 PDIFF = 0. 4018 PSPEC = 0. 4019 PSP = 0. 4020 ELSE 4021 PDIFF = LDOTN * PRIMAR*DIFFUS 4022 PSP = 2.*LDOTN*NL(3) - SOURCE(3) 4023* Comments as for SSPEC apply, but substitute PRIMAR for STRAIT: 4024 IF (PSP.LT.PHOBND .OR. IPHONG.EQ.0) THEN 4025 PSPEC = 0. 4026 ELSE 4027 PSPEC = PSP**IPHONG * PRIMAR*SPECLR 4028 ENDIF 4029 ENDIF 4030* 4031* experience has shown the "spots" on dark objects to be rather 4032* overpowering, especially by comparison to those on brighter 4033* objects. hence the specular reflections on dark objects are 4034* now artificially scaled down by a function which relates 4035* directly to the "brightness" of the object. 4036* this makes such objects duller, but their 4037* colour seems to come through more clearly, and they don't 4038* appear more specular than the brighter objects. 4039* the funny coefficients are ntsc: 4040 BRIGHT = 0.2 + 0.8 * SQRT(0.299*RGBCUR(1) + 4041 & 0.587*RGBCUR(2) + 4042 & 0.114*RGBCUR(3)) 4043 SSPEC = SSPEC * BRIGHT 4044 PSPEC = PSPEC * BRIGHT 4045* 4046* The usual case is white lighting, no transparency 4047* 4048 SPECOL(1) = 1.0 4049 SPECOL(2) = 1.0 4050 SPECOL(3) = 1.0 4051 SBLEND = 1.0 4052 CLRITY = 0.0 4053* 4054* Extra properties make specular highlighting calculation a 4055* bit more complex. First we have to find the MATERIAL description. 4056* 4057 IF (iand(FLAG(INDTOP),PROPS).NE.0) THEN 4058 K = FLAG(INDTOP) / 65536 4059 if(K.le.0)then 4060 WRITE(NOISE,*)"FLAG(",INDTOP-1,") =",FLAG(INDTOP-1) 4061 WRITE(NOISE,*)"FLAG(",INDTOP,") =",FLAG(INDTOP) 4062 WRITE(NOISE,*)"FLAG(",INDTOP+1,") =",FLAG(INDTOP+1) 4063 endif 4064 CALL ASSERT(K.GT.0,'lost material definition') 4065 IF (iand(FLAG(MLIST(K)),INSIDE).NE.0 .AND. 4066 & iand(FLAG(INDTOP), INSIDE).NE.0) THEN 4067 K = LIST(MLIST(K)) 4068 MPHONG = DETAIL(K+14) 4069 SPECM = DETAIL(K+15) 4070 ELSE 4071 K = LIST(MLIST(K)) 4072 MPHONG = DETAIL(K+1) 4073 SPECM = DETAIL(K+2) 4074 ENDIF 4075 SPECOL(1) = DETAIL(K+3) 4076 SPECOL(2) = DETAIL(K+4) 4077 SPECOL(3) = DETAIL(K+5) 4078 IF (SPECOL(1).LT.0) SPECOL(1) = RGBCUR(1) 4079 IF (SPECOL(2).LT.0) SPECOL(2) = RGBCUR(2) 4080 IF (SPECOL(3).LT.0) SPECOL(3) = RGBCUR(3) 4081 CLRITY = DETAIL(K+6) 4082* not currently used, as MOPT(1)=1 already marked in FLAG, 4083* but future interpretations of MOPT(1) might need it 4084 CLROPT = DETAIL(K+7) 4085 ENDIF 4086* 4087* 20-Feb-2000 Allow per-vertex transparency (obj type 18) 4088 IF (iand(FLAG(INDTOP),VTRANS).NE.0) THEN 4089 IF (TYPE(INDTOP).EQ.CYLIND) THEN 4090 K = LIST(INDTOP) 4091 FRAC = DETAIL(K+8) 4092 INEXT = INDTOP + 1 4093 IF (TYPE(INEXT).EQ.VERTEXRGB) INEXT = INEXT + 1 4094 CALL ASSERT(TYPE(INEXT).EQ.VERTRANSP,'lost vertex transp') 4095 K = LIST(INEXT) 4096 CLRITY = FRAC*DETAIL(K+1) + (1.-FRAC)*DETAIL(K+2) 4097 ELSE IF (TYPE(INDTOP).EQ.TRIANG) THEN 4098 K = LIST(INDTOP) 4099 ALPHA = DETAIL(K+14) 4100 BETA = DETAIL(K+15) 4101 INEXT = INDTOP + 1 4102 IF ((TYPE(INEXT).EQ.VERTEXRGB).OR.(TYPE(INEXT).EQ.NORMS)) 4103 & INEXT = INEXT + 1 4104 IF ((TYPE(INEXT).EQ.VERTEXRGB).OR.(TYPE(INEXT).EQ.NORMS)) 4105 & INEXT = INEXT + 1 4106 CALL ASSERT(TYPE(INEXT).EQ.VERTRANSP,'lost vertex transp') 4107 K = LIST(INEXT) 4108 CLRITY = DETAIL(K+1) 4109 & + ALPHA * (DETAIL(K+2)-DETAIL(K+1)) 4110 & + BETA * (DETAIL(K+3)-DETAIL(K+1)) 4111 CALL ASSERT(CLRITY.LE.1 .AND. CLRITY.GE.0.0, 'illegal transp') 4112 ELSE 4113 CALL ASSERT(.FALSE.,'illegal vertex transp') 4114 ENDIF 4115 ENDIF 4116 4117C 4118C This is the only computationally intensive code (as opposed to mere 4119C bookkeeping) involved in rendering transparent objects. The blend 4120C factor must be some function of the clarity/transparency, but I'm not 4121C sure exactly what the equation ought to be. The cosine function in 4122C TRNSPOPT option 0 below was chosen after purely empirical tests of the 4123C resulting image quality. If your machine bogs down incredibly due to 4124C the cosine call, then you might prefer to use TRNSPOPT 1 instead. 4125C Conversely, if you don't mind the extra computation then you could use 4126C TRNSPOPT 2, which is closer to an ideal model. 4127C Dec 2010: controlled by OPT(2) in the MATERIAL specification record 4128C 4129 IF (CLRITY.NE.0) THEN 4130 ZN = ABS(NL(3)) + MODULO(TRNSPOPT,1.0) 4131 SELECT CASE(FLOOR(TRNSPOPT)) 4132 CASE DEFAULT ! CASE 0 4133 SBLEND = .25*(1.+COS(3.1416*CLRITY*ZN))**2 4134 CASE (1) 4135 SBLEND = (1. - CLRITY*ZN)**2 4136 CASE (2) 4137 SBLEND = 1.0 - CLRITY ** ( 0.7071 / ZN ) 4138 CASE (3) 4139 SBLEND = 1.0 - CLRITY 4140 END SELECT 4141 ENDIF 4142C 4143C Final calculation of specular properties of special materials 4144C 4145 IF (iand(FLAG(INDTOP),PROPS).NE.0) THEN 4146 DIFFM = 1. - (SPECM + AMBIEN) 4147 SDIFF = SDIFF * DIFFM / DIFFUS 4148 PDIFF = PDIFF * DIFFM / DIFFUS 4149 SSPEC = 0.0 4150 PSPEC = 0.0 4151 IF (SSP .GE. PHOBND) 4152 & SSPEC = SSP**MPHONG * STRAIT*SPECM 4153 IF ((PSP .GE. PHOBND) .AND. (LDOTN.GT.0)) 4154 & PSPEC = PSP**MPHONG * PRIMAR*SPECM 4155* de-emphasize highlights from inside surface of transparent objects 4156* Could use BACKFACE flag instead of INSIDE to catch non-triangles 4157 IF (iand(FLAG(INDTOP),INSIDE).NE.0) THEN 4158 SSPEC = SSPEC * (1.-SPECM) 4159 PSPEC = PSPEC * (1.-SPECM) 4160 ENDIF 4161 END IF 4162* 4163* We now return you to your regular processing 4164* 4165 DO 280 KC = 1, 3 4166 C2ND = SBLEND*RGBCUR(KC)*(AMBIEN+SDIFF) + SSPEC*SPECOL(KC) 4167 CSUN = SBLEND*RGBCUR(KC)*PDIFF + PSPEC*SPECOL(KC) 4168 RGBSHD(KC) = C2ND 4169 RGBFUL(KC) = C2ND + CSUN 4170280 CONTINUE 4171 4172C EAM March 1997 - Support additional non-shadowing light sources 4173C which lie within figure and have a finite range of illumination. 4174 IF (NGLOWS.GT.0) THEN 4175 DO KC = 1,3 4176 RGBSHD(KC) = (1.-GLOWMAX)*RGBSHD(KC) 4177 RGBFUL(KC) = (1.-GLOWMAX)*RGBFUL(KC) 4178 ENDDO 4179* Recover glow light parameters 4180 DO NG = 1, NGLOWS 4181 IG = LIST(GLOWLIST(NG)) 4182 GLOWSRC(1) = DETAIL(IG+1) 4183 GLOWSRC(2) = DETAIL(IG+2) 4184 GLOWSRC(3) = DETAIL(IG+3) 4185 GLOWRAD = DETAIL(IG+4) 4186 GLOW = DETAIL(IG+5) 4187 GOPT = DETAIL(IG+6) 4188 GPHONG = DETAIL(IG+7) 4189 GLOWCOL(1) = DETAIL(IG+8) 4190 GLOWCOL(2) = DETAIL(IG+9) 4191 GLOWCOL(3) = DETAIL(IG+10) 4192* 4193 GDIST(1) = GLOWSRC(1) - XP 4194 GDIST(2) = GLOWSRC(2) - YP 4195 GDIST(3) = GLOWSRC(3) - ZP 4196 ABSN = SQRT(GDIST(1)**2 + GDIST(2)**2 + GDIST(3)**2) 4197 GDIST(1) = GDIST(1) / ABSN 4198 GDIST(2) = GDIST(2) / ABSN 4199 GDIST(3) = GDIST(3) / ABSN 4200 LDOTN = GDIST(1)*NL(1) + GDIST(2)*NL(2) + GDIST(3)*NL(3) 4201 IF (LDOTN.LE.0) THEN 4202 GDIFF = 0. 4203 GSPEC = 0. 4204 ELSE 4205C Might want separate diffuse param for glow; (always 1.0?) 4206C GDIFF = LDOTN * DIFFUS 4207 GDIFF = LDOTN 4208 GSP = 2.*LDOTN*NL(3) - GDIST(3) 4209 IF (GSP.LT.PHOBND .OR. GPHONG.EQ.0) THEN 4210 GSPEC = 0. 4211 ELSE 4212 GSPEC = GSP**GPHONG * SPECLR 4213 ENDIF 4214 ENDIF 4215C Limit glow effect by some function of ABSN, GLOWRAD 4216 IF (GOPT.EQ.3) THEN 4217 GFADE = MAX( 0., 1. - ABSN/GLOWRAD ) 4218 ELSE IF (GOPT.EQ.2) THEN 4219 GFADE = 1./(ABSN/GLOWRAD + 1.) 4220 ELSE IF (GOPT.EQ.1) THEN 4221 GFADE = 1./(ABSN/GLOWRAD + 1.)**2 4222 ELSE 4223 GFADE = MIN( 1., 1./(ABSN/GLOWRAD)**2 ) 4224 ENDIF 4225 DO KC = 1,3 4226C This isn't right for transparent surfaces 4227 CGLO = SBLEND*RGBCUR(KC)*GDIFF + GSPEC 4228 CGLO = GFADE * GLOWCOL(KC) * CGLO 4229 RGBSHD(KC) = RGBSHD(KC) + CGLO 4230 RGBFUL(KC) = RGBFUL(KC) + CGLO 4231 ENDDO 4232* End of this glow light 4233 ENDDO 4234 ENDIF 4235* 4236* That does it for the shading computation. 4237* ZS should still be a shadow-space co-ordinate of the pixel 4238* whose shading we were interested in, and zstop should be a 4239* shadow-space object's co-ordinate no further than that from 4240* the primary light source (modulo the empirical slop factor). 4241* 4242 IF (INDTOP.EQ.INDSTP) THEN 4243 TILE(1,I,J) = RGBFUL(1) 4244 TILE(2,I,J) = RGBFUL(2) 4245 TILE(3,I,J) = RGBFUL(3) 4246 ELSE IF (ZS+ZSLOP.GE.ZSTOP) THEN 4247 nzslop = nzslop + 1 4248 TILE(1,I,J) = RGBFUL(1) 4249 TILE(2,I,J) = RGBFUL(2) 4250 TILE(3,I,J) = RGBFUL(3) 4251 ELSE 4252 TILE(1,I,J) = RGBSHD(1) 4253 TILE(2,I,J) = RGBSHD(2) 4254 TILE(3,I,J) = RGBSHD(3) 4255 ENDIF 4256C 4257C Fog processing added July 1998; moved into transp. loop 2010. 4258C Note: Background fog is applied at beginning of this loop. 4259C Should have glow lights brighten fog? 4260 IF (FOGTYPE .GE. 0) THEN 4261 FOGDIM = FOGGY( FOGLIM(2) - ZP ) 4262 TILE(1,I,J) = (1.-FOGDIM)*TILE(1,I,J) + FOGDIM*FOGRGB(1) 4263 TILE(2,I,J) = (1.-FOGDIM)*TILE(2,I,J) + FOGDIM*FOGRGB(2) 4264 TILE(3,I,J) = (1.-FOGDIM)*TILE(3,I,J) + FOGDIM*FOGRGB(3) 4265 ENDIF 4266 4267C 4268C Transparency processing totally overhauled Feb 2001 4269C The first pass is sufficient if top object is opaque. 4270 IF (NTRANSP .EQ. 0) GOTO 299 4271 IF (INDEPTH.EQ.1 .AND. iand(FLAG(INDTOP),TRANSP).EQ.0) GOTO 299 4272 IF (ITPASS.EQ.1) THEN 4273 ACHAN(I,J) = SBLEND 4274 ELSE 4275 ACHAN(I,J) = 1. - (1.-ACHAN(I,J))*(1.-SBLEND) 4276 ENDIF 4277 RGBLND(1) = (1.-SBLEND)*RGBLND(1) + TILE(1,I,J) 4278 RGBLND(2) = (1.-SBLEND)*RGBLND(2) + TILE(2,I,J) 4279 RGBLND(3) = (1.-SBLEND)*RGBLND(3) + TILE(3,I,J) 4280C CALL ASSERT(ITPASS.LE.INDEPTH,'Ran off end of INDEPTH') 4281 IF (ITPASS.GE.INDEPTH) THEN 4282 TILE(1,I,J) = RGBLND(1) 4283 TILE(2,I,J) = RGBLND(2) 4284 TILE(3,I,J) = RGBLND(3) 4285 ZTOP = ZP 4286 ELSE 4287 ZP = ZLIST(INDEPTH-ITPASS) 4288 INDTOP = INDLIST(INDEPTH-ITPASS) 4289 NORMAL(1) = NORMLIST(1,INDEPTH-ITPASS) 4290 NORMAL(2) = NORMLIST(2,INDEPTH-ITPASS) 4291 NORMAL(3) = NORMLIST(3,INDEPTH-ITPASS) 4292 ITPASS = ITPASS + 1 4293 GOTO 255 4294 ENDIF 4295C End of transparency processing 4296C 4297299 CONTINUE 4298 4299C 4300300 CONTINUE 4301400 CONTINUE 4302* do tile averaging and save output tile in outbuf 4303C For now fold schemes 0 and 1 together; later split for efficiency? 4304 IF (SCHEME.LE.1) THEN 4305 K = (ITILE-1)*NOX 4306 DO 420 J = 1, NOY 4307 DO 415 I = 1, NOX 4308 K = K + 1 4309C CALL ASSERT (K.LE.SIZE(OUTBUF,1),'k>outsiz') 4310C 4311 DO 410 IC = 1, 3 4312 ICK = 256. * SQRT(TILE(IC,I,J)) 4313 IF (ICK.LT.0) ICK = 0 4314 IF (ICK.GT.255) ICK = 255 4315 IF (GAMMACORRECTION) ICK = GAMMA_MAP(ICK+1) 4316 OUTBUF(K,IC) = ICK 4317 410 CONTINUE 4318C 4319 IF (SCHEME.EQ.0) THEN 4320 ICK = 255. * ACHAN(I,J) 4321 IF (ICK.LT.0) ICK = 0 4322 IF (ICK.GT.255) ICK = 255 4323 OUTBUF(K,4) = ICK 4324 END IF 4325C 4326415 CONTINUE 4327 K = K + NOX*(NTX-1) 4328420 CONTINUE 4329 ELSEIF (SCHEME.EQ.2) THEN 4330 K = (ITILE-1)*NOX 4331 DO 440 J = 1, NOY 4332 DO 435 I = 1, NOX 4333 K = K + 1 4334 DO 430 IC = 1, 3 4335* I'm not quite convinced by this pixel averaging 4336* (is a corner worth too much in this setup?): 4337 TMP = (TILE(IC,2*I-1,2*J-1) + 4338 & TILE(IC,2*I ,2*J-1) + 4339 & TILE(IC,2*I-1,2*J ) + 4340 & TILE(IC,2*I ,2*J )) / 4. 4341 ICK = 256. * SQRT(TMP) 4342 IF (ICK.LT.0) ICK = 0 4343 IF (ICK.GT.255) ICK = 255 4344 IF (GAMMACORRECTION) ICK = GAMMA_MAP(ICK+1) 4345 OUTBUF(K,IC) = ICK 4346430 CONTINUE 4347435 CONTINUE 4348 K = K + NOX*(NTX-1) 4349440 CONTINUE 4350 ELSEIF (SCHEME.EQ.3) THEN 4351 NHX = NOX/2 4352 NHY = NOY/2 4353 K = (ITILE-1)*NOX 4354 DO 460 J = 1, NHY 4355 DO 455 I = 1, NHX 4356 DO 450 IC = 1, 3 4357* Bad pixel averaging?: 4358 TMP1 = (TILE(IC,3*I-2,3*J-2) + 4359 & TILE(IC,3*I-1,3*J-2)/2. + 4360 & TILE(IC,3*I-2,3*J-1)/2. + 4361 & TILE(IC,3*I-1,3*J-1)/4.) / (9./4.) 4362 TMP2 = (TILE(IC,3*I-1,3*J-2)/2. + 4363 & TILE(IC,3*I ,3*J-2) + 4364 & TILE(IC,3*I-1,3*J-1)/4. + 4365 & TILE(IC,3*I ,3*J-1)/2.) / (9./4.) 4366 TMP3 = (TILE(IC,3*I-2,3*J-1)/2. + 4367 & TILE(IC,3*I-1,3*J-1)/4. + 4368 & TILE(IC,3*I-2,3*J ) + 4369 & TILE(IC,3*I-1,3*J )/2.) / (9./4.) 4370 TMP4 = (TILE(IC,3*I-1,3*J-1)/4. + 4371 & TILE(IC,3*I ,3*J-1)/2. + 4372 & TILE(IC,3*I-1,3*J )/2. + 4373 & TILE(IC,3*I ,3*J ) ) / (9./4.) 4374 ICK1 = MIN(MAX(INT(256.*SQRT(TMP1)),0),255) 4375 ICK2 = MIN(MAX(INT(256.*SQRT(TMP2)),0),255) 4376 ICK3 = MIN(MAX(INT(256.*SQRT(TMP3)),0),255) 4377 ICK4 = MIN(MAX(INT(256.*SQRT(TMP4)),0),255) 4378 IF (GAMMACORRECTION) THEN 4379 ICK1 = GAMMA_MAP(ICK1+1) 4380 ICK2 = GAMMA_MAP(ICK2+1) 4381 ICK3 = GAMMA_MAP(ICK3+1) 4382 ICK4 = GAMMA_MAP(ICK4+1) 4383 ENDIF 4384 OUTBUF( K+1,IC) = ICK1 4385 OUTBUF( K+2,IC) = ICK2 4386 OUTBUF(NX+K+1,IC) = ICK3 4387 OUTBUF(NX+K+2,IC) = ICK4 4388450 CONTINUE 4389 K = K + 2 4390455 CONTINUE 4391 K = K + NOX*(2*NTX - 1) 4392460 CONTINUE 4393 ELSE 4394 CALL ASSERT(.FALSE.,'crash 500') 4395 ENDIF 4396500 CONTINUE 4397* Ready to write when we have completed a row of tiles 4398 K = 0 4399 DO 550 J=1,NOY 4400 LINOUT = LINOUT + 1 4401 IF (LINOUT.GT.NAY) GOTO 600 4402 IERR = LOCAL(2, OUTBUF(K+1,1), OUTBUF(K+1,2), OUTBUF(K+1,3), 4403 & OUTBUF(K+1,4) ) 4404 K = K + NX 4405C CALL ASSERT (K.LE.SIZE(OUTBUF,1),'k>outsiz') 4406550 CONTINUE 4407600 CONTINUE 4408* 4409* Report any soft failures 4410 IF ( NSXMAX.GT.0 .OR. NSYMAX.GT.0 4411 & .OR. TRANOVFL.GT.0 ) THEN 4412 WRITE(NOISE,*)' >>> WARNINGS <<<' 4413 END IF 4414 IF (NSXMAX.GT.0) WRITE(NOISE,*) 4415 & ' Possible shadow error NSXMAX=',NSXMAX 4416 IF (NSYMAX.GT.0) WRITE(NOISE,*) 4417 & ' Possible shadow error NSYMAX=',NSYMAX 4418 IF (TRANOVFL.GT.0) WRITE(NOISE,601) 4419 & ' Transparency processing truncated at MAXTRANSP=', 4420 & MAXTRANSP, ' for', TRANOVFL,' pixels' 4421601 FORMAT(A,I3,A,I10,A) 4422* 4423* Debugging information 4424 if (verbose) then 4425 if (nzslop.gt.0) write(noise,*)' NZSLOP failures=',nzslop 4426 if (nslow.gt.0) write(noise,*)' NSLOW failures=',nslow 4427 endif 4428* 4429* close up shop 4430 IERR = LOCAL(3) 4431* 4432 END 4433 4434 SUBROUTINE TRANSF (X,Y,Z) 4435* Input transformation 4436 COMMON /MATRICES/ XCENT, YCENT, SCALE, EYEPOS, SXCENT, SYCENT, 4437 & TMAT, TINV, TINVT, SROT, SRTINV, SRTINVT 4438 & ,RAFTER, TAFTER 4439 REAL XCENT, YCENT, SCALE, EYEPOS, SXCENT, SYCENT 4440* Transformation matrix, inverse of transpose, and transposed inverse 4441 REAL TMAT(4,4), TINV(4,4), TINVT(4,4) 4442* Shortest rotation from light source to +z axis 4443 REAL SROT(4,4), SRTINV(4,4), SRTINVT(4,4) 4444* Post-hoc transformation on top of original TMAT 4445 REAL RAFTER(4,4), TAFTER(3) 4446 REAL X,Y,Z 4447 REAL G(4),H(4) 4448 H(1) = X*TMAT(1,1) + Y*TMAT(2,1) + Z*TMAT(3,1) + TMAT(4,1) 4449 H(2) = X*TMAT(1,2) + Y*TMAT(2,2) + Z*TMAT(3,2) + TMAT(4,2) 4450 H(3) = X*TMAT(1,3) + Y*TMAT(2,3) + Z*TMAT(3,3) + TMAT(4,3) 4451 H(4) = X*TMAT(1,4) + Y*TMAT(2,4) + Z*TMAT(3,4) + TMAT(4,4) 4452* Apply post-hoc rotation and translation also 4453 G(1) = RAFTER(1,1)*H(1) + RAFTER(1,2)*H(2) + RAFTER(1,3)*H(3) 4454 & + TAFTER(1) 4455 G(2) = RAFTER(2,1)*H(1) + RAFTER(2,2)*H(2) + RAFTER(2,3)*H(3) 4456 & + TAFTER(2) 4457 G(3) = RAFTER(3,1)*H(1) + RAFTER(3,2)*H(2) + RAFTER(3,3)*H(3) 4458 & + TAFTER(3) 4459 CALL ASSERT (H(4).NE.0.,'infinite vector') 4460 X = G(1) / H(4) 4461 Y = G(2) / H(4) 4462 Z = G(3) / H(4) 4463 RETURN 4464 END 4465 4466 SUBROUTINE ISOLATE( X, Y ) 4467* Expand X and Y coordinates to fill image regardless of aspect ratio 4468 COMMON /OPTIONS/ FONTSCALE, GAMMA, ZOOM, NSCHEME, SHADOWFLAG, XBG, 4469 & NAX, NAY, OTMODE, QUALITY, INVERT, LFLAG 4470 REAL FONTSCALE, GAMMA, ZOOM 4471 INTEGER NSCHEME, SHADOWFLAG, XBG 4472 INTEGER*4 NAX, NAY, OTMODE, QUALITY 4473 LOGICAL*2 INVERT, LFLAG 4474 COMMON /MATRICES/ XCENT, YCENT, SCALE, EYEPOS, SXCENT, SYCENT, 4475 & TMAT, TINV, TINVT, SROT, SRTINV, SRTINVT 4476 & ,RAFTER, TAFTER 4477 REAL XCENT, YCENT, SCALE, EYEPOS, SXCENT, SYCENT 4478 COMMON /NICETIES/ TRULIM, ZLIM, FRONTCLIP, BACKCLIP 4479 & , ISOLATION 4480 REAL TRULIM(3,2), ZLIM(2), FRONTCLIP, BACKCLIP 4481 INTEGER ISOLATION 4482* 4483 IF (INVERT) Y = -Y 4484 IF (ISOLATION.EQ.2) THEN 4485 ASPECT = XCENT/YCENT 4486 IF (ASPECT.GT.1.0) X = X * ASPECT 4487 IF (ASPECT.LT.1.0) Y = Y / ASPECT 4488 ENDIF 4489 RETURN 4490 END 4491 4492 4493 4494 SUBROUTINE HSORTD (N, A, NDEX) 4495 INTEGER N 4496 REAL A(N) 4497 INTEGER NDEX(N) 4498* 4499* this formulation of heapsort is based on n. wirth, 4500* "algorithms + data structures = programs" (p. 75). 4501* 4502* the caller supplies an array, a, containing n elements, and an 4503* index array with space for n integers. 4504* a and n are considered "read-only" by the subroutine, but ndex 4505* is filled by the subroutine with the sequence of indices of a 4506* that obtain the elements of a in ascending order. this is 4507* similar to the apl unary "tree" operator. thus a(ndex(1)) is the 4508* smallest element after the sort, and a(ndex(n)) is the largest. 4509* 4510 INTEGER L, R, T 4511 DO 10 I = 1, N 451210 NDEX(I) = I 4513 L = N/2 + 1 4514 R = N 451520 IF (L .LE. 1) GO TO 30 4516 L = L - 1 4517 CALL HSIFTD (N, A, NDEX, L, R) 4518 GO TO 20 451930 IF (R .LE. 1) RETURN 4520 T = NDEX(1) 4521 NDEX(1) = NDEX(R) 4522 NDEX(R) = T 4523 R = R - 1 4524 CALL HSIFTD (N, A, NDEX, L, R) 4525 GO TO 30 4526 END 4527 SUBROUTINE HSIFTD (N, A, NDEX, L, R) 4528* used by hsortd 4529 INTEGER N 4530 REAL A(N) 4531 INTEGER NDEX(N), L, R 4532 INTEGER I, J, X 4533 I = L 4534 J = I + I 4535 X = NDEX(I) 453610 IF (J .GT. R) GO TO 20 4537 IF (J .GE. R) GO TO 15 4538 IF (A(NDEX(J)) .LT. A(NDEX(J+1))) J = J + 1 453915 IF (A(X) .GE. A(NDEX(J))) GO TO 20 4540 NDEX(I) = NDEX(J) 4541 I = J 4542 J = I + I 4543 GO TO 10 454420 NDEX(I) = X 4545 RETURN 4546 END 4547 4548 SUBROUTINE PLANER (X1,Y1,Z1,X2,Y2,Z2,X3,Y3,Z3, A,B,C,D) 4549 IMPLICIT REAL (A-Z) 4550* solve for coefficients of plane eqn z=ax+by+c 4551* and yield d=0 in case of degenerate ("edge-on") triangle 4552 D1 = Z1*(Y2-Y3) - Y1*(Z2-Z3) + Z2*Y3-Y2*Z3 4553 D2 = X1*(Z2-Z3) - Z1*(X2-X3) + X2*Z3-Z2*X3 4554 D3 = X1*(Y2*Z3-Z2*Y3) - Y1*(X2*Z3-Z2*X3) + Z1*(X2*Y3-Y2*X3) 4555 D = X1*(Y2-Y3) - Y1*(X2-X3) + X2*Y3-Y2*X3 4556 A = 0. 4557 B = 0. 4558 C = 0. 4559 IF (D.NE.0.) THEN 4560 A = D1/D 4561 B = D2/D 4562 C = D3/D 4563 ENDIF 4564 IF (ABS(A)+ABS(B)+ABS(C).GT.1E10) THEN 4565 D = 0.0 4566 ENDIF 4567 RETURN 4568 END 4569 4570 SUBROUTINE ASSERT (LOGIC, DAMMIT) 4571 LOGICAL LOGIC 4572 CHARACTER*(*) DAMMIT 4573 INTEGER ASSOUT 4574 LOGICAL VERBOSE 4575 COMMON /ASSCOM/ ASSOUT, VERBOSE 4576 SAVE /ASSCOM/ 4577 IF (LOGIC) RETURN 4578 WRITE (ASSOUT,*) 'ERROR >>>>>> ',DAMMIT 4579 CALL EXIT(-1) 4580 END 4581 4582C 4583C Find Z coord of point on surface of cylinder with known X and Y coords 4584C cylinder axis is X2 - X1, cylinder radius is R 4585C Need to find Z coord ZB. 4586C flag is 0 if cylinder had rounded ends, FLAT if it has flat ends, 4587C Also find nearest point XYZA on cylinder axis and fraction along it. 4588C 4589 FUNCTION CYLTEST ( flag, axfrac, 4590 & x1,y1,z1, x2,y2,z2, xb,yb,zb, R, xa,ya,za ) 4591 LOGICAL CYLTEST 4592c implicit NONE 4593* 4594* Bit definitions for FLAG array 4595 INTEGER FLAT, RIBBON, SURFACE, PROPS 4596 PARAMETER (FLAT=2, RIBBON=4, SURFACE=8, PROPS=16) 4597 INTEGER TRANSP, HIDDEN, INSIDE, MOPT1 4598 PARAMETER (TRANSP=32, HIDDEN=64, INSIDE=128,MOPT1=256) 4599 INTEGER VCOLS, CLIPPED, VTRANS, BOUNDED 4600 PARAMETER (VCOLS=512, CLIPPED=1024, VTRANS=2048, BOUNDED=4096) 4601 INTEGER TFI 4602 PARAMETER (TFI = TRANSP + INSIDE) 4603c 4604 INTEGER flag 4605 REAL x1,y1,z1, x2,y2,z2, xb,yb,zb 4606 REAL R 4607 REAL xa,ya,za 4608c 4609 REAL ca,cb,cg,dx,dy,dz,d2 4610 REAL A0,A1,A2,Q 4611 REAL p1,r2,dx2,dy2 4612 REAL dd1,dd2 4613c 4614c start with direction cosines * d2 4615 ca = x2 - x1 4616 cb = y2 - y1 4617 cg = z2 - z1 4618c 4619c other useful quantities 4620c (note: if d2==0 must be degenerate cylinder, really a disk) 4621 r2 = R*R 4622 d2 = ca*ca + cb*cb + cg*cg 4623 dx = xb - x1 4624 dy = yb - y1 4625 dx2 = dx**2 4626 dy2 = dy**2 4627c 4628c use these to find coefficients of quadratic equation for ZB 4629c EAM Jan 1997 test and handle dx-dy=0 4630 if (ca.eq.0. .and. cb.eq.0.) then 4631 if (z2.gt.z1) p1 = 1.0 4632 if (z2.lt.z1) p1 = -1.0 4633 goto 100 4634 end if 4635c 4636 A0 = (dx*cb - dy*ca)**2 + (dy2 + dx2)*cg*cg - r2*d2 4637 A1 = -2.0 * (dy*cg*cb + dx*ca*cg) 4638 A2 = ca*ca + cb*cb 4639 Q = A1*A1 - 4.0*A0*A2 4640 if (Q .lt. 0) then 4641C zb = -99999. 4642 cyltest = .false. 4643 return 4644 else 4645 if ( iand(flag,TFI) .eq. TFI) then 4646 dz = (-sqrt(Q) - A1) / (2.0 * A2) 4647 else 4648 dz = ( sqrt(Q) - A1) / (2.0 * A2) 4649 endif 4650 zb = z1 + dz 4651 end if 4652c 4653c now find nearest point on cylinder axis 4654c p1 is fraction along axis from x1 to x2 4655c 0 < p1 < 1 means point is on wall of cylinder 4656c 4657 dd1 = dx2 + dy2 + dz*dz 4658 dd2 = (x2-xb)**2 + (y2-yb)**2 + (z2-zb)**2 4659c 4660 p1 = (dd1 - r2) / d2 4661 if (p1 .le. 0.0) then 4662 p1 = 0.0 4663 else 4664 p1 = sqrt(p1) 4665 end if 4666c 4667 if ((dd2 .gt. (d2+r2)) .and. (dd2 .gt. dd1)) p1 = -p1 4668c 4669 if (p1 .ge. 0 .and. p1 .le. 1.0) then 4670 xa = p1*ca + x1 4671 ya = p1*cb + y1 4672 za = p1*cg + z1 4673 cyltest = .true. 4674 return 4675 end if 4676c 4677c point is either on end cap, or missed entirely 4678c 4679 100 continue 4680 if (p1 .ge. 1.0) then 4681 xa = x2 4682 ya = y2 4683 za = z2 4684 dx = xb - x2 4685 dy = yb - y2 4686 dx2 = dx**2 4687 dy2 = dy**2 4688 else if (p1 .le. 0.0) then 4689 xa = x1 4690 ya = y1 4691 za = z1 4692 end if 4693c 4694c Rounded cylinder end 4695 if (iand(flag,FLAT) .eq. 0) then 4696 if (dx2+dy2 .gt. r2) then 4697 cyltest = .false. 4698 return 4699 else 4700 if ( iand(flag,TFI) .eq. TFI) then 4701 zb = za - sqrt(r2 - (dx2+dy2)) 4702 else 4703 zb = za + sqrt(r2 - (dx2+dy2)) 4704 end if 4705 end if 4706C 4707C Flat cylinder end 4708C 4709 else 4710 if (cg .eq. 0.) then 4711C zb = -99999. 4712 cyltest = .false. 4713 return 4714 endif 4715 zb = (cg*za - ca*dx - cb*dy) / cg 4716 if (dx2 + dy2 + (zb-za)**2 .ge. r2) then 4717C zb = -99999. 4718 cyltest = .false. 4719 return 4720 endif 4721 if (p1 .ge. 1.0) then 4722 xa = xb - (x2 - x1) 4723 ya = yb - (y2 - y1) 4724 za = zb - (z2 - z1) 4725 else if (p1 .le. 0.0) then 4726 xa = xb - (x1 - x2) 4727 ya = yb - (y1 - y2) 4728 za = zb - (z1 - z2) 4729 endif 4730 end if 4731c 4732 if (p1.gt.1.0) p1 = 1.0 4733 if (p1.lt.0.0) p1 = 0.0 4734 axfrac = p1 4735c 4736 cyltest = .true. 4737 return 4738 end 4739 4740C Bookkeeping for transparency 4741C 5-Mar-2001 4742C New version that is not limited to 3-deep transparent objects. 4743C On exit: 4744C INDTOP, ZTOP contain the top object so far, and its height 4745C ZHIGH height of top opaque object 4746C 4747 subroutine rank( ind, zp, normal, flag ) 4748* 4749 implicit NONE 4750 integer ind 4751 real zp, normal(3) 4752 integer*4 flag(1) 4753* 4754 integer i,j,k 4755* 4756 INTEGER ASSOUT 4757 LOGICAL VERBOSE 4758 COMMON /ASSCOM/ ASSOUT, VERBOSE 4759* 4760 INCLUDE 'parameters.incl' 4761* 4762* Support for transparency 4763 COMMON /TRANS/ NTRANSP, INDEPTH, INDTOP, TRANOVFL, ZTOP, ZHIGH, 4764 & INDLIST(MAXTRANSP), ZLIST(MAXTRANSP), 4765 & NORMLIST(3,MAXTRANSP) 4766 INTEGER NTRANSP, INDEPTH, INDTOP, INDLIST, TRANOVFL 4767 REAL ZTOP, ZHIGH, ZLIST, NORMLIST 4768* 4769* Bit definitions for FLAG(MAXOBJ) array 4770 INTEGER FLAT, RIBBON, SURFACE, PROPS 4771 PARAMETER (FLAT=2, RIBBON=4, SURFACE=8, PROPS=16) 4772 INTEGER TRANSP, HIDDEN, INSIDE, MOPT1 4773 PARAMETER (TRANSP=32, HIDDEN=64, INSIDE=128,MOPT1=256) 4774 INTEGER VCOLS, CLIPPED, VTRANS, BOUNDED 4775 PARAMETER (VCOLS=512, CLIPPED=1024, VTRANS=2048, BOUNDED=4096) 4776* 4777* The MOPT1 flag signals an alternative mode of transparency 4778 if (iand(flag(ind),mopt1).ne.0) goto 400 4779* 4780 do i = 1, indepth 4781 if (zp.gt.zlist(i)) then 4782 if (iand(flag(ind),transp).eq.0) then 4783 indepth = i 4784 zhigh = zp 4785 goto 345 4786 else 4787 do j = indepth, i, -1 4788 indlist(j+1) = indlist(j) 4789 zlist(j+1) = zlist(j) 4790 normlist(1,j+1)= normlist(1,j) 4791 normlist(2,j+1)= normlist(2,j) 4792 normlist(3,j+1)= normlist(3,j) 4793 enddo 4794 indepth = indepth + 1 4795 goto 344 4796 endif 4797 else if (iand(flag(indlist(i)),TRANSP).eq.0) then 4798 return 4799 endif 4800 enddo 4801c If the rest of the list is transparent, add this at the end 4802 indepth = indepth + 1 4803 i = indepth 4804c 4805 344 continue 4806 if (indepth.ge.MAXTRANSP) then 4807 indepth = MAXTRANSP - 1 4808 tranovfl = tranovfl + 1 4809 return 4810 endif 4811c 4812 345 continue 4813 indlist(i) = ind 4814 zlist(i) = zp 4815 normlist(1,i) = normal(1) 4816 normlist(2,i) = normal(2) 4817 normlist(3,i) = normal(3) 4818 if (iand(flag(ind),transp).eq.0) zhigh = zp 4819 indtop = indlist(1) 4820 ztop = zlist(1) 4821 return 4822c 4823c MOPT1 version. Same routine, except this time we have the extra 4824c overhead of having to check for duplication of material. 4825 400 continue 4826 do i = 1, indepth 4827 if (zp.gt.zlist(i)) then 4828 if (flag(ind)/65536 .eq. flag(indlist(i))/65536) then 4829 goto 345 4830 endif 4831c Handle case where two MOPT1 surfaces have intervening transp obj 4832c In this case overwrite the lower MOPT1 surface 4833 do k = i, indepth-1 4834 if (flag(ind)/65536 .eq. flag(indlist(k+1))/65536) 4835 & goto 401 4836 enddo 4837 if (indepth .ge. MAXTRANSP-1) then 4838 tranovfl = tranovfl + 1 4839 else 4840 k = indepth 4841 indepth = indepth + 1 4842 endif 4843 401 continue 4844 do j = k, i, -1 4845 indlist(j+1) = indlist(j) 4846 zlist(j+1) = zlist(j) 4847 normlist(1,j+1)= normlist(1,j) 4848 normlist(2,j+1)= normlist(2,j) 4849 normlist(3,j+1)= normlist(3,j) 4850 enddo 4851 goto 345 4852 else if (iand(flag(indlist(i)),TRANSP).eq.0) then 4853 return 4854 else if (flag(ind)/65536 .eq. flag(indlist(i))/65536) then 4855 return 4856 endif 4857 enddo 4858c If the rest of the list is transparent, add this at the end 4859 indepth = indepth + 1 4860 i = indepth 4861 goto 344 4862c 4863 end 4864c 4865 4866 SUBROUTINE CHKRGB( RED, GRN, BLU, MESSAGE ) 4867 REAL RED, GRN, BLU 4868 CHARACTER*(*) MESSAGE 4869 CALL ASSERT (RED.GE.0., MESSAGE) 4870 CALL ASSERT (GRN.GE.0., MESSAGE) 4871 CALL ASSERT (BLU.GE.0., MESSAGE) 4872 CALL ASSERT (RED.LE.1., MESSAGE) 4873 CALL ASSERT (GRN.LE.1., MESSAGE) 4874 CALL ASSERT (BLU.LE.1., MESSAGE) 4875 RETURN 4876 END 4877 4878 4879 FUNCTION FOGGY( DEPTH ) 4880 REAL FOGGY, DEPTH 4881 COMMON /FOGCOM/ FOGTYPE,FOGFRONT,FOGBACK,FOGDEN,FOGLIM,FOGRGB 4882 INTEGER FOGTYPE 4883 REAL FOGFRONT, FOGBACK, FOGDEN, FOGLIM(2), FOGRGB(3) 4884 REAL FOGDIM 4885c 4886 FOGDIM = 0.5 4887 IF (FOGTYPE .EQ. 0) 4888 & FOGDIM = FOGDEN * DEPTH / (FOGLIM(2)-FOGLIM(1)) 4889 IF (FOGTYPE .GT. 0) 4890 & FOGDIM = 1. - EXP(-FOGDEN * DEPTH/(FOGLIM(2)-FOGLIM(1))) 4891 FOGDIM = MAX( 0.0, FOGDIM ) 4892 FOGDIM = MIN( 1.0, FOGDIM ) 4893 FOGGY = FOGDIM 4894 RETURN 4895 END 4896 4897 FUNCTION DET( A ) 4898 REAL DET, A(4,4) 4899 DET = A(1,1)*A(2,2)*A(3,3) + A(1,2)*A(2,3)*A(3,1) 4900 & + A(2,1)*A(3,2)*A(1,3) - A(1,1)*A(2,3)*A(3,2) 4901 & - A(3,3)*A(1,2)*A(2,1) - A(1,3)*A(2,2)*A(3,1) 4902 RETURN 4903 END 4904 4905 4906 subroutine liblookup( name, fullname ) 4907c 4908 character*(*) name 4909 character*132 fullname 4910 character*132 R3DLIB 4911c 4912 call getenv('R3D_LIB',R3DLIB) 4913c 4914 fullname = ' ' 4915 len = 0 4916 do i = 1, 132 4917 if (R3DLIB(i:i).ne.' ') len = i 4918 enddo 4919 if (len.eq.0) then 4920 fullname = name 4921 return 4922 else 4923 fullname(1:len) = R3DLIB(1:len) 4924 fullname(len+1:len+1) = '/' 4925 j = len+2 4926 endif 4927c 4928 len = 0 4929 100 continue 4930 len = len + 1 4931 if (name(len:len).ne.' ') goto 100 4932 len = len - 1 4933c 4934 fullname(j:j+len-1) = name(1:len) 4935 len = j+len-1 4936c 4937 return 4938 end 4939 4940C EAM - 21 Feb 2001 Version 2.6 4941C Test for bounding planes 4942C Returns .FALSE. if the point does not have to be rendered 4943C .TRUE. if the point is to be rendered, in which case 4944C ZP, DX, DY, DZ have been updated if on bounding plane 4945C ZP = height of point being rendered 4946C DX,DY,DZ = surface normal at point (XP,YP,ZP) 4947C BPIND is the object number of the bounding plane 4948C or 0 if the bounding planes don't trigger 4949C 4950 FUNCTION INBOUNDS( MAT, TYPE, LIST, DETAIL, 4951 & XP,YP,ZP, DX,DY,DZ, ZBACK, BPIND ) 4952 IMPLICIT NONE 4953 LOGICAL INBOUNDS 4954 INTEGER MAT, M, TYPE(1), LIST(1), BPIND 4955 REAL DETAIL(1), XP,YP,ZP, DX,DY,DZ, ZBACK 4956c 4957c Object type used for bounding planes (may change) 4958 INTEGER INTERNAL 4959 PARAMETER (INTERNAL = 4) 4960c 4961 REAL BPLANE(3), BPNORM(3) 4962 REAL XN,YN,ZN, TEMP 4963 INTEGER MIND 4964 LOGICAL TESTBACK 4965c 4966 INBOUNDS = .FALSE. 4967 BPIND = 0 4968 TESTBACK = (ZBACK.NE.ZP) 4969 M = MAT 4970 DO WHILE (TYPE(M).EQ.INTERNAL) 4971 MIND = LIST(M) 4972 BPLANE(1) = DETAIL(MIND+3) 4973 BPLANE(2) = DETAIL(MIND+4) 4974 BPLANE(3) = DETAIL(MIND+5) 4975 BPNORM(1) = DETAIL(MIND+6) 4976 BPNORM(2) = DETAIL(MIND+7) 4977 BPNORM(3) = DETAIL(MIND+8) 4978 XN = XP - BPLANE(1) 4979 YN = YP - BPLANE(2) 4980 ZN = ZP - BPLANE(3) 4981 TEMP = XN*BPNORM(1) + YN*BPNORM(2) + ZN*BPNORM(3) 4982 IF (TEMP.GT.0) THEN 4983 IF (TESTBACK) THEN 4984 ZP = ZBACK 4985 ZN = ZP - BPLANE(3) 4986 TEMP = XN*BPNORM(1) + YN*BPNORM(2) + ZN*BPNORM(3) 4987 IF (TEMP.GT.0) RETURN 4988 ENDIF 4989 DX = BPNORM(1) 4990 DY = BPNORM(2) 4991 DZ = BPNORM(3) 4992 ZP = (XN*DX + YN*DY) / (-DZ) + BPLANE(3) 4993 BPIND = M 4994 ENDIF 4995 M = M + 1 4996 ENDDO 4997 INBOUNDS = .TRUE. 4998 RETURN 4999 END 5000C 5001C Equivalent test for bounding planes in ORTEP mode (only the octant clipped 5002C by all three bounding planes is removed). Only intended for ellipsoids. 5003C This tests the AND (rather than the OR) of multiple bounding planes. 5004C 5005 SUBROUTINE ORTEPBOUNDS( MAT, TYPE, LIST, DETAIL, 5006 & XP,YP,ZP, DX,DY,DZ, ZBACK, BPIND ) 5007 IMPLICIT NONE 5008 INTEGER MAT, M, TYPE(1), LIST(1), BPIND 5009 REAL DETAIL(1), XP,YP,ZP, DX,DY,DZ, ZBACK 5010c 5011c Object type used for bounding planes (may change) 5012 INTEGER INTERNAL 5013 PARAMETER (INTERNAL = 4) 5014c 5015 REAL BPLANE(3), BPNORM(3) 5016 REAL XN,YN,ZN, TEMP, ZMAX, DXMAX, DYMAX, DZMAX 5017 INTEGER MIND 5018c 5019 M = MAT 5020 ZMAX = -1.E10 5021 DXMAX = DX 5022 DYMAX = DY 5023 DZMAX = DZ 5024 DO WHILE (TYPE(M).EQ.INTERNAL) 5025 MIND = LIST(M) 5026 BPLANE(1) = DETAIL(MIND+3) 5027 BPLANE(2) = DETAIL(MIND+4) 5028 BPLANE(3) = DETAIL(MIND+5) 5029 BPNORM(1) = DETAIL(MIND+6) 5030 BPNORM(2) = DETAIL(MIND+7) 5031 BPNORM(3) = DETAIL(MIND+8) 5032 XN = XP - BPLANE(1) 5033 YN = YP - BPLANE(2) 5034 ZN = ZP - BPLANE(3) 5035 TEMP = XN*BPNORM(1) + YN*BPNORM(2) + ZN*BPNORM(3) 5036 IF (TEMP.LT.0) THEN 5037 BPIND = 0 5038 RETURN 5039 ENDIF 5040 TEMP = (XN*BPNORM(1)+YN*BPNORM(2)) / (-BPNORM(3)) + BPLANE(3) 5041 IF (TEMP.GT.ZMAX) THEN 5042 ZMAX = TEMP 5043 DXMAX = BPNORM(1) 5044 DYMAX = BPNORM(2) 5045 DZMAX = BPNORM(3) 5046 BPIND = M 5047 ENDIF 5048 M = M + 1 5049 ENDDO 5050 ZP = ZMAX 5051 DX = DXMAX 5052 DY = DYMAX 5053 DZ = DZMAX 5054 CALL ASSERT(DZ.GT.0,'ORTEP bounds incorrectly initialized') 5055 RETURN 5056 END 5057 5058 SUBROUTINE GET_TRY(NOLD, NNEW, TRY, DIMENS) 5059* 5060* NOLD = old memory allocation 5061* NNEW = new memory allocation needed 5062* TRY = new memory allocations beyond that needed to be attempted, 5063* so that we don't have to expand the array too often. 5064* DIMENS = number of dimensions to be expanded (1 to 3) 5065* 5066* For DIMENS=1: if the new allocation is less than twice the old 5067* try to double the old allocation; otherwise add 50% to the new. 5068* If that fails, add 50% to the old or 25% of the new; 5069* If that fails, just try the new, then give up. 5070* 5071* For DIMENS = 2, expand by the square root of 2, 1.5 or 1.25 5072* for the same effect on total memory in 2-D arrays where bot 5073* dimensions are expanded. 5074* For DIMENS = 3 (not currently used) use the cube root. 5075* 5076 INTEGER NOLD, NNEW, TRY(3), N, DIMENS 5077 REAL RATIO(3, 3) 5078 DATA RATIO / 2.0, 1.5, 1.25, 5079 & 1.41, 1.22, 1.12, 5080 & 1.26, 1.14, 1.07 / 5081 N = INT( FLOAT(NOLD) * RATIO(1, DIMENS) ) 5082 if (NNEW.LT.N) THEN 5083 TRY(1) = N 5084 N = INT( FLOAT(NOLD) * RATIO(2, DIMENS) ) 5085 if (NNEW.LT.N) THEN 5086 TRY(2) = N 5087 else 5088 TRY(2) = INT( FLOAT(NNEW) * RATIO(3, DIMENS) ) 5089 ENDIF 5090 else 5091 TRY(1) = INT( FLOAT(NNEW) * RATIO(2, DIMENS) ) 5092 TRY(2) = INT( FLOAT(NNEW) * RATIO(3, DIMENS) ) 5093 endif 5094 TRY(3) = NNEW 5095 RETURN 5096 END 5097