File: render.f

package info (click to toggle)
raster3d 3.0-3-2
  • links: PTS, VCS
  • area: main
  • in suites: stretch
  • size: 4,916 kB
  • ctags: 1,557
  • sloc: fortran: 9,536; ansic: 1,060; makefile: 318; sh: 250; csh: 15
file content (5096 lines) | stat: -rw-r--r-- 173,773 bytes parent folder | download | duplicates (3)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
1001
1002
1003
1004
1005
1006
1007
1008
1009
1010
1011
1012
1013
1014
1015
1016
1017
1018
1019
1020
1021
1022
1023
1024
1025
1026
1027
1028
1029
1030
1031
1032
1033
1034
1035
1036
1037
1038
1039
1040
1041
1042
1043
1044
1045
1046
1047
1048
1049
1050
1051
1052
1053
1054
1055
1056
1057
1058
1059
1060
1061
1062
1063
1064
1065
1066
1067
1068
1069
1070
1071
1072
1073
1074
1075
1076
1077
1078
1079
1080
1081
1082
1083
1084
1085
1086
1087
1088
1089
1090
1091
1092
1093
1094
1095
1096
1097
1098
1099
1100
1101
1102
1103
1104
1105
1106
1107
1108
1109
1110
1111
1112
1113
1114
1115
1116
1117
1118
1119
1120
1121
1122
1123
1124
1125
1126
1127
1128
1129
1130
1131
1132
1133
1134
1135
1136
1137
1138
1139
1140
1141
1142
1143
1144
1145
1146
1147
1148
1149
1150
1151
1152
1153
1154
1155
1156
1157
1158
1159
1160
1161
1162
1163
1164
1165
1166
1167
1168
1169
1170
1171
1172
1173
1174
1175
1176
1177
1178
1179
1180
1181
1182
1183
1184
1185
1186
1187
1188
1189
1190
1191
1192
1193
1194
1195
1196
1197
1198
1199
1200
1201
1202
1203
1204
1205
1206
1207
1208
1209
1210
1211
1212
1213
1214
1215
1216
1217
1218
1219
1220
1221
1222
1223
1224
1225
1226
1227
1228
1229
1230
1231
1232
1233
1234
1235
1236
1237
1238
1239
1240
1241
1242
1243
1244
1245
1246
1247
1248
1249
1250
1251
1252
1253
1254
1255
1256
1257
1258
1259
1260
1261
1262
1263
1264
1265
1266
1267
1268
1269
1270
1271
1272
1273
1274
1275
1276
1277
1278
1279
1280
1281
1282
1283
1284
1285
1286
1287
1288
1289
1290
1291
1292
1293
1294
1295
1296
1297
1298
1299
1300
1301
1302
1303
1304
1305
1306
1307
1308
1309
1310
1311
1312
1313
1314
1315
1316
1317
1318
1319
1320
1321
1322
1323
1324
1325
1326
1327
1328
1329
1330
1331
1332
1333
1334
1335
1336
1337
1338
1339
1340
1341
1342
1343
1344
1345
1346
1347
1348
1349
1350
1351
1352
1353
1354
1355
1356
1357
1358
1359
1360
1361
1362
1363
1364
1365
1366
1367
1368
1369
1370
1371
1372
1373
1374
1375
1376
1377
1378
1379
1380
1381
1382
1383
1384
1385
1386
1387
1388
1389
1390
1391
1392
1393
1394
1395
1396
1397
1398
1399
1400
1401
1402
1403
1404
1405
1406
1407
1408
1409
1410
1411
1412
1413
1414
1415
1416
1417
1418
1419
1420
1421
1422
1423
1424
1425
1426
1427
1428
1429
1430
1431
1432
1433
1434
1435
1436
1437
1438
1439
1440
1441
1442
1443
1444
1445
1446
1447
1448
1449
1450
1451
1452
1453
1454
1455
1456
1457
1458
1459
1460
1461
1462
1463
1464
1465
1466
1467
1468
1469
1470
1471
1472
1473
1474
1475
1476
1477
1478
1479
1480
1481
1482
1483
1484
1485
1486
1487
1488
1489
1490
1491
1492
1493
1494
1495
1496
1497
1498
1499
1500
1501
1502
1503
1504
1505
1506
1507
1508
1509
1510
1511
1512
1513
1514
1515
1516
1517
1518
1519
1520
1521
1522
1523
1524
1525
1526
1527
1528
1529
1530
1531
1532
1533
1534
1535
1536
1537
1538
1539
1540
1541
1542
1543
1544
1545
1546
1547
1548
1549
1550
1551
1552
1553
1554
1555
1556
1557
1558
1559
1560
1561
1562
1563
1564
1565
1566
1567
1568
1569
1570
1571
1572
1573
1574
1575
1576
1577
1578
1579
1580
1581
1582
1583
1584
1585
1586
1587
1588
1589
1590
1591
1592
1593
1594
1595
1596
1597
1598
1599
1600
1601
1602
1603
1604
1605
1606
1607
1608
1609
1610
1611
1612
1613
1614
1615
1616
1617
1618
1619
1620
1621
1622
1623
1624
1625
1626
1627
1628
1629
1630
1631
1632
1633
1634
1635
1636
1637
1638
1639
1640
1641
1642
1643
1644
1645
1646
1647
1648
1649
1650
1651
1652
1653
1654
1655
1656
1657
1658
1659
1660
1661
1662
1663
1664
1665
1666
1667
1668
1669
1670
1671
1672
1673
1674
1675
1676
1677
1678
1679
1680
1681
1682
1683
1684
1685
1686
1687
1688
1689
1690
1691
1692
1693
1694
1695
1696
1697
1698
1699
1700
1701
1702
1703
1704
1705
1706
1707
1708
1709
1710
1711
1712
1713
1714
1715
1716
1717
1718
1719
1720
1721
1722
1723
1724
1725
1726
1727
1728
1729
1730
1731
1732
1733
1734
1735
1736
1737
1738
1739
1740
1741
1742
1743
1744
1745
1746
1747
1748
1749
1750
1751
1752
1753
1754
1755
1756
1757
1758
1759
1760
1761
1762
1763
1764
1765
1766
1767
1768
1769
1770
1771
1772
1773
1774
1775
1776
1777
1778
1779
1780
1781
1782
1783
1784
1785
1786
1787
1788
1789
1790
1791
1792
1793
1794
1795
1796
1797
1798
1799
1800
1801
1802
1803
1804
1805
1806
1807
1808
1809
1810
1811
1812
1813
1814
1815
1816
1817
1818
1819
1820
1821
1822
1823
1824
1825
1826
1827
1828
1829
1830
1831
1832
1833
1834
1835
1836
1837
1838
1839
1840
1841
1842
1843
1844
1845
1846
1847
1848
1849
1850
1851
1852
1853
1854
1855
1856
1857
1858
1859
1860
1861
1862
1863
1864
1865
1866
1867
1868
1869
1870
1871
1872
1873
1874
1875
1876
1877
1878
1879
1880
1881
1882
1883
1884
1885
1886
1887
1888
1889
1890
1891
1892
1893
1894
1895
1896
1897
1898
1899
1900
1901
1902
1903
1904
1905
1906
1907
1908
1909
1910
1911
1912
1913
1914
1915
1916
1917
1918
1919
1920
1921
1922
1923
1924
1925
1926
1927
1928
1929
1930
1931
1932
1933
1934
1935
1936
1937
1938
1939
1940
1941
1942
1943
1944
1945
1946
1947
1948
1949
1950
1951
1952
1953
1954
1955
1956
1957
1958
1959
1960
1961
1962
1963
1964
1965
1966
1967
1968
1969
1970
1971
1972
1973
1974
1975
1976
1977
1978
1979
1980
1981
1982
1983
1984
1985
1986
1987
1988
1989
1990
1991
1992
1993
1994
1995
1996
1997
1998
1999
2000
2001
2002
2003
2004
2005
2006
2007
2008
2009
2010
2011
2012
2013
2014
2015
2016
2017
2018
2019
2020
2021
2022
2023
2024
2025
2026
2027
2028
2029
2030
2031
2032
2033
2034
2035
2036
2037
2038
2039
2040
2041
2042
2043
2044
2045
2046
2047
2048
2049
2050
2051
2052
2053
2054
2055
2056
2057
2058
2059
2060
2061
2062
2063
2064
2065
2066
2067
2068
2069
2070
2071
2072
2073
2074
2075
2076
2077
2078
2079
2080
2081
2082
2083
2084
2085
2086
2087
2088
2089
2090
2091
2092
2093
2094
2095
2096
2097
2098
2099
2100
2101
2102
2103
2104
2105
2106
2107
2108
2109
2110
2111
2112
2113
2114
2115
2116
2117
2118
2119
2120
2121
2122
2123
2124
2125
2126
2127
2128
2129
2130
2131
2132
2133
2134
2135
2136
2137
2138
2139
2140
2141
2142
2143
2144
2145
2146
2147
2148
2149
2150
2151
2152
2153
2154
2155
2156
2157
2158
2159
2160
2161
2162
2163
2164
2165
2166
2167
2168
2169
2170
2171
2172
2173
2174
2175
2176
2177
2178
2179
2180
2181
2182
2183
2184
2185
2186
2187
2188
2189
2190
2191
2192
2193
2194
2195
2196
2197
2198
2199
2200
2201
2202
2203
2204
2205
2206
2207
2208
2209
2210
2211
2212
2213
2214
2215
2216
2217
2218
2219
2220
2221
2222
2223
2224
2225
2226
2227
2228
2229
2230
2231
2232
2233
2234
2235
2236
2237
2238
2239
2240
2241
2242
2243
2244
2245
2246
2247
2248
2249
2250
2251
2252
2253
2254
2255
2256
2257
2258
2259
2260
2261
2262
2263
2264
2265
2266
2267
2268
2269
2270
2271
2272
2273
2274
2275
2276
2277
2278
2279
2280
2281
2282
2283
2284
2285
2286
2287
2288
2289
2290
2291
2292
2293
2294
2295
2296
2297
2298
2299
2300
2301
2302
2303
2304
2305
2306
2307
2308
2309
2310
2311
2312
2313
2314
2315
2316
2317
2318
2319
2320
2321
2322
2323
2324
2325
2326
2327
2328
2329
2330
2331
2332
2333
2334
2335
2336
2337
2338
2339
2340
2341
2342
2343
2344
2345
2346
2347
2348
2349
2350
2351
2352
2353
2354
2355
2356
2357
2358
2359
2360
2361
2362
2363
2364
2365
2366
2367
2368
2369
2370
2371
2372
2373
2374
2375
2376
2377
2378
2379
2380
2381
2382
2383
2384
2385
2386
2387
2388
2389
2390
2391
2392
2393
2394
2395
2396
2397
2398
2399
2400
2401
2402
2403
2404
2405
2406
2407
2408
2409
2410
2411
2412
2413
2414
2415
2416
2417
2418
2419
2420
2421
2422
2423
2424
2425
2426
2427
2428
2429
2430
2431
2432
2433
2434
2435
2436
2437
2438
2439
2440
2441
2442
2443
2444
2445
2446
2447
2448
2449
2450
2451
2452
2453
2454
2455
2456
2457
2458
2459
2460
2461
2462
2463
2464
2465
2466
2467
2468
2469
2470
2471
2472
2473
2474
2475
2476
2477
2478
2479
2480
2481
2482
2483
2484
2485
2486
2487
2488
2489
2490
2491
2492
2493
2494
2495
2496
2497
2498
2499
2500
2501
2502
2503
2504
2505
2506
2507
2508
2509
2510
2511
2512
2513
2514
2515
2516
2517
2518
2519
2520
2521
2522
2523
2524
2525
2526
2527
2528
2529
2530
2531
2532
2533
2534
2535
2536
2537
2538
2539
2540
2541
2542
2543
2544
2545
2546
2547
2548
2549
2550
2551
2552
2553
2554
2555
2556
2557
2558
2559
2560
2561
2562
2563
2564
2565
2566
2567
2568
2569
2570
2571
2572
2573
2574
2575
2576
2577
2578
2579
2580
2581
2582
2583
2584
2585
2586
2587
2588
2589
2590
2591
2592
2593
2594
2595
2596
2597
2598
2599
2600
2601
2602
2603
2604
2605
2606
2607
2608
2609
2610
2611
2612
2613
2614
2615
2616
2617
2618
2619
2620
2621
2622
2623
2624
2625
2626
2627
2628
2629
2630
2631
2632
2633
2634
2635
2636
2637
2638
2639
2640
2641
2642
2643
2644
2645
2646
2647
2648
2649
2650
2651
2652
2653
2654
2655
2656
2657
2658
2659
2660
2661
2662
2663
2664
2665
2666
2667
2668
2669
2670
2671
2672
2673
2674
2675
2676
2677
2678
2679
2680
2681
2682
2683
2684
2685
2686
2687
2688
2689
2690
2691
2692
2693
2694
2695
2696
2697
2698
2699
2700
2701
2702
2703
2704
2705
2706
2707
2708
2709
2710
2711
2712
2713
2714
2715
2716
2717
2718
2719
2720
2721
2722
2723
2724
2725
2726
2727
2728
2729
2730
2731
2732
2733
2734
2735
2736
2737
2738
2739
2740
2741
2742
2743
2744
2745
2746
2747
2748
2749
2750
2751
2752
2753
2754
2755
2756
2757
2758
2759
2760
2761
2762
2763
2764
2765
2766
2767
2768
2769
2770
2771
2772
2773
2774
2775
2776
2777
2778
2779
2780
2781
2782
2783
2784
2785
2786
2787
2788
2789
2790
2791
2792
2793
2794
2795
2796
2797
2798
2799
2800
2801
2802
2803
2804
2805
2806
2807
2808
2809
2810
2811
2812
2813
2814
2815
2816
2817
2818
2819
2820
2821
2822
2823
2824
2825
2826
2827
2828
2829
2830
2831
2832
2833
2834
2835
2836
2837
2838
2839
2840
2841
2842
2843
2844
2845
2846
2847
2848
2849
2850
2851
2852
2853
2854
2855
2856
2857
2858
2859
2860
2861
2862
2863
2864
2865
2866
2867
2868
2869
2870
2871
2872
2873
2874
2875
2876
2877
2878
2879
2880
2881
2882
2883
2884
2885
2886
2887
2888
2889
2890
2891
2892
2893
2894
2895
2896
2897
2898
2899
2900
2901
2902
2903
2904
2905
2906
2907
2908
2909
2910
2911
2912
2913
2914
2915
2916
2917
2918
2919
2920
2921
2922
2923
2924
2925
2926
2927
2928
2929
2930
2931
2932
2933
2934
2935
2936
2937
2938
2939
2940
2941
2942
2943
2944
2945
2946
2947
2948
2949
2950
2951
2952
2953
2954
2955
2956
2957
2958
2959
2960
2961
2962
2963
2964
2965
2966
2967
2968
2969
2970
2971
2972
2973
2974
2975
2976
2977
2978
2979
2980
2981
2982
2983
2984
2985
2986
2987
2988
2989
2990
2991
2992
2993
2994
2995
2996
2997
2998
2999
3000
3001
3002
3003
3004
3005
3006
3007
3008
3009
3010
3011
3012
3013
3014
3015
3016
3017
3018
3019
3020
3021
3022
3023
3024
3025
3026
3027
3028
3029
3030
3031
3032
3033
3034
3035
3036
3037
3038
3039
3040
3041
3042
3043
3044
3045
3046
3047
3048
3049
3050
3051
3052
3053
3054
3055
3056
3057
3058
3059
3060
3061
3062
3063
3064
3065
3066
3067
3068
3069
3070
3071
3072
3073
3074
3075
3076
3077
3078
3079
3080
3081
3082
3083
3084
3085
3086
3087
3088
3089
3090
3091
3092
3093
3094
3095
3096
3097
3098
3099
3100
3101
3102
3103
3104
3105
3106
3107
3108
3109
3110
3111
3112
3113
3114
3115
3116
3117
3118
3119
3120
3121
3122
3123
3124
3125
3126
3127
3128
3129
3130
3131
3132
3133
3134
3135
3136
3137
3138
3139
3140
3141
3142
3143
3144
3145
3146
3147
3148
3149
3150
3151
3152
3153
3154
3155
3156
3157
3158
3159
3160
3161
3162
3163
3164
3165
3166
3167
3168
3169
3170
3171
3172
3173
3174
3175
3176
3177
3178
3179
3180
3181
3182
3183
3184
3185
3186
3187
3188
3189
3190
3191
3192
3193
3194
3195
3196
3197
3198
3199
3200
3201
3202
3203
3204
3205
3206
3207
3208
3209
3210
3211
3212
3213
3214
3215
3216
3217
3218
3219
3220
3221
3222
3223
3224
3225
3226
3227
3228
3229
3230
3231
3232
3233
3234
3235
3236
3237
3238
3239
3240
3241
3242
3243
3244
3245
3246
3247
3248
3249
3250
3251
3252
3253
3254
3255
3256
3257
3258
3259
3260
3261
3262
3263
3264
3265
3266
3267
3268
3269
3270
3271
3272
3273
3274
3275
3276
3277
3278
3279
3280
3281
3282
3283
3284
3285
3286
3287
3288
3289
3290
3291
3292
3293
3294
3295
3296
3297
3298
3299
3300
3301
3302
3303
3304
3305
3306
3307
3308
3309
3310
3311
3312
3313
3314
3315
3316
3317
3318
3319
3320
3321
3322
3323
3324
3325
3326
3327
3328
3329
3330
3331
3332
3333
3334
3335
3336
3337
3338
3339
3340
3341
3342
3343
3344
3345
3346
3347
3348
3349
3350
3351
3352
3353
3354
3355
3356
3357
3358
3359
3360
3361
3362
3363
3364
3365
3366
3367
3368
3369
3370
3371
3372
3373
3374
3375
3376
3377
3378
3379
3380
3381
3382
3383
3384
3385
3386
3387
3388
3389
3390
3391
3392
3393
3394
3395
3396
3397
3398
3399
3400
3401
3402
3403
3404
3405
3406
3407
3408
3409
3410
3411
3412
3413
3414
3415
3416
3417
3418
3419
3420
3421
3422
3423
3424
3425
3426
3427
3428
3429
3430
3431
3432
3433
3434
3435
3436
3437
3438
3439
3440
3441
3442
3443
3444
3445
3446
3447
3448
3449
3450
3451
3452
3453
3454
3455
3456
3457
3458
3459
3460
3461
3462
3463
3464
3465
3466
3467
3468
3469
3470
3471
3472
3473
3474
3475
3476
3477
3478
3479
3480
3481
3482
3483
3484
3485
3486
3487
3488
3489
3490
3491
3492
3493
3494
3495
3496
3497
3498
3499
3500
3501
3502
3503
3504
3505
3506
3507
3508
3509
3510
3511
3512
3513
3514
3515
3516
3517
3518
3519
3520
3521
3522
3523
3524
3525
3526
3527
3528
3529
3530
3531
3532
3533
3534
3535
3536
3537
3538
3539
3540
3541
3542
3543
3544
3545
3546
3547
3548
3549
3550
3551
3552
3553
3554
3555
3556
3557
3558
3559
3560
3561
3562
3563
3564
3565
3566
3567
3568
3569
3570
3571
3572
3573
3574
3575
3576
3577
3578
3579
3580
3581
3582
3583
3584
3585
3586
3587
3588
3589
3590
3591
3592
3593
3594
3595
3596
3597
3598
3599
3600
3601
3602
3603
3604
3605
3606
3607
3608
3609
3610
3611
3612
3613
3614
3615
3616
3617
3618
3619
3620
3621
3622
3623
3624
3625
3626
3627
3628
3629
3630
3631
3632
3633
3634
3635
3636
3637
3638
3639
3640
3641
3642
3643
3644
3645
3646
3647
3648
3649
3650
3651
3652
3653
3654
3655
3656
3657
3658
3659
3660
3661
3662
3663
3664
3665
3666
3667
3668
3669
3670
3671
3672
3673
3674
3675
3676
3677
3678
3679
3680
3681
3682
3683
3684
3685
3686
3687
3688
3689
3690
3691
3692
3693
3694
3695
3696
3697
3698
3699
3700
3701
3702
3703
3704
3705
3706
3707
3708
3709
3710
3711
3712
3713
3714
3715
3716
3717
3718
3719
3720
3721
3722
3723
3724
3725
3726
3727
3728
3729
3730
3731
3732
3733
3734
3735
3736
3737
3738
3739
3740
3741
3742
3743
3744
3745
3746
3747
3748
3749
3750
3751
3752
3753
3754
3755
3756
3757
3758
3759
3760
3761
3762
3763
3764
3765
3766
3767
3768
3769
3770
3771
3772
3773
3774
3775
3776
3777
3778
3779
3780
3781
3782
3783
3784
3785
3786
3787
3788
3789
3790
3791
3792
3793
3794
3795
3796
3797
3798
3799
3800
3801
3802
3803
3804
3805
3806
3807
3808
3809
3810
3811
3812
3813
3814
3815
3816
3817
3818
3819
3820
3821
3822
3823
3824
3825
3826
3827
3828
3829
3830
3831
3832
3833
3834
3835
3836
3837
3838
3839
3840
3841
3842
3843
3844
3845
3846
3847
3848
3849
3850
3851
3852
3853
3854
3855
3856
3857
3858
3859
3860
3861
3862
3863
3864
3865
3866
3867
3868
3869
3870
3871
3872
3873
3874
3875
3876
3877
3878
3879
3880
3881
3882
3883
3884
3885
3886
3887
3888
3889
3890
3891
3892
3893
3894
3895
3896
3897
3898
3899
3900
3901
3902
3903
3904
3905
3906
3907
3908
3909
3910
3911
3912
3913
3914
3915
3916
3917
3918
3919
3920
3921
3922
3923
3924
3925
3926
3927
3928
3929
3930
3931
3932
3933
3934
3935
3936
3937
3938
3939
3940
3941
3942
3943
3944
3945
3946
3947
3948
3949
3950
3951
3952
3953
3954
3955
3956
3957
3958
3959
3960
3961
3962
3963
3964
3965
3966
3967
3968
3969
3970
3971
3972
3973
3974
3975
3976
3977
3978
3979
3980
3981
3982
3983
3984
3985
3986
3987
3988
3989
3990
3991
3992
3993
3994
3995
3996
3997
3998
3999
4000
4001
4002
4003
4004
4005
4006
4007
4008
4009
4010
4011
4012
4013
4014
4015
4016
4017
4018
4019
4020
4021
4022
4023
4024
4025
4026
4027
4028
4029
4030
4031
4032
4033
4034
4035
4036
4037
4038
4039
4040
4041
4042
4043
4044
4045
4046
4047
4048
4049
4050
4051
4052
4053
4054
4055
4056
4057
4058
4059
4060
4061
4062
4063
4064
4065
4066
4067
4068
4069
4070
4071
4072
4073
4074
4075
4076
4077
4078
4079
4080
4081
4082
4083
4084
4085
4086
4087
4088
4089
4090
4091
4092
4093
4094
4095
4096
4097
4098
4099
4100
4101
4102
4103
4104
4105
4106
4107
4108
4109
4110
4111
4112
4113
4114
4115
4116
4117
4118
4119
4120
4121
4122
4123
4124
4125
4126
4127
4128
4129
4130
4131
4132
4133
4134
4135
4136
4137
4138
4139
4140
4141
4142
4143
4144
4145
4146
4147
4148
4149
4150
4151
4152
4153
4154
4155
4156
4157
4158
4159
4160
4161
4162
4163
4164
4165
4166
4167
4168
4169
4170
4171
4172
4173
4174
4175
4176
4177
4178
4179
4180
4181
4182
4183
4184
4185
4186
4187
4188
4189
4190
4191
4192
4193
4194
4195
4196
4197
4198
4199
4200
4201
4202
4203
4204
4205
4206
4207
4208
4209
4210
4211
4212
4213
4214
4215
4216
4217
4218
4219
4220
4221
4222
4223
4224
4225
4226
4227
4228
4229
4230
4231
4232
4233
4234
4235
4236
4237
4238
4239
4240
4241
4242
4243
4244
4245
4246
4247
4248
4249
4250
4251
4252
4253
4254
4255
4256
4257
4258
4259
4260
4261
4262
4263
4264
4265
4266
4267
4268
4269
4270
4271
4272
4273
4274
4275
4276
4277
4278
4279
4280
4281
4282
4283
4284
4285
4286
4287
4288
4289
4290
4291
4292
4293
4294
4295
4296
4297
4298
4299
4300
4301
4302
4303
4304
4305
4306
4307
4308
4309
4310
4311
4312
4313
4314
4315
4316
4317
4318
4319
4320
4321
4322
4323
4324
4325
4326
4327
4328
4329
4330
4331
4332
4333
4334
4335
4336
4337
4338
4339
4340
4341
4342
4343
4344
4345
4346
4347
4348
4349
4350
4351
4352
4353
4354
4355
4356
4357
4358
4359
4360
4361
4362
4363
4364
4365
4366
4367
4368
4369
4370
4371
4372
4373
4374
4375
4376
4377
4378
4379
4380
4381
4382
4383
4384
4385
4386
4387
4388
4389
4390
4391
4392
4393
4394
4395
4396
4397
4398
4399
4400
4401
4402
4403
4404
4405
4406
4407
4408
4409
4410
4411
4412
4413
4414
4415
4416
4417
4418
4419
4420
4421
4422
4423
4424
4425
4426
4427
4428
4429
4430
4431
4432
4433
4434
4435
4436
4437
4438
4439
4440
4441
4442
4443
4444
4445
4446
4447
4448
4449
4450
4451
4452
4453
4454
4455
4456
4457
4458
4459
4460
4461
4462
4463
4464
4465
4466
4467
4468
4469
4470
4471
4472
4473
4474
4475
4476
4477
4478
4479
4480
4481
4482
4483
4484
4485
4486
4487
4488
4489
4490
4491
4492
4493
4494
4495
4496
4497
4498
4499
4500
4501
4502
4503
4504
4505
4506
4507
4508
4509
4510
4511
4512
4513
4514
4515
4516
4517
4518
4519
4520
4521
4522
4523
4524
4525
4526
4527
4528
4529
4530
4531
4532
4533
4534
4535
4536
4537
4538
4539
4540
4541
4542
4543
4544
4545
4546
4547
4548
4549
4550
4551
4552
4553
4554
4555
4556
4557
4558
4559
4560
4561
4562
4563
4564
4565
4566
4567
4568
4569
4570
4571
4572
4573
4574
4575
4576
4577
4578
4579
4580
4581
4582
4583
4584
4585
4586
4587
4588
4589
4590
4591
4592
4593
4594
4595
4596
4597
4598
4599
4600
4601
4602
4603
4604
4605
4606
4607
4608
4609
4610
4611
4612
4613
4614
4615
4616
4617
4618
4619
4620
4621
4622
4623
4624
4625
4626
4627
4628
4629
4630
4631
4632
4633
4634
4635
4636
4637
4638
4639
4640
4641
4642
4643
4644
4645
4646
4647
4648
4649
4650
4651
4652
4653
4654
4655
4656
4657
4658
4659
4660
4661
4662
4663
4664
4665
4666
4667
4668
4669
4670
4671
4672
4673
4674
4675
4676
4677
4678
4679
4680
4681
4682
4683
4684
4685
4686
4687
4688
4689
4690
4691
4692
4693
4694
4695
4696
4697
4698
4699
4700
4701
4702
4703
4704
4705
4706
4707
4708
4709
4710
4711
4712
4713
4714
4715
4716
4717
4718
4719
4720
4721
4722
4723
4724
4725
4726
4727
4728
4729
4730
4731
4732
4733
4734
4735
4736
4737
4738
4739
4740
4741
4742
4743
4744
4745
4746
4747
4748
4749
4750
4751
4752
4753
4754
4755
4756
4757
4758
4759
4760
4761
4762
4763
4764
4765
4766
4767
4768
4769
4770
4771
4772
4773
4774
4775
4776
4777
4778
4779
4780
4781
4782
4783
4784
4785
4786
4787
4788
4789
4790
4791
4792
4793
4794
4795
4796
4797
4798
4799
4800
4801
4802
4803
4804
4805
4806
4807
4808
4809
4810
4811
4812
4813
4814
4815
4816
4817
4818
4819
4820
4821
4822
4823
4824
4825
4826
4827
4828
4829
4830
4831
4832
4833
4834
4835
4836
4837
4838
4839
4840
4841
4842
4843
4844
4845
4846
4847
4848
4849
4850
4851
4852
4853
4854
4855
4856
4857
4858
4859
4860
4861
4862
4863
4864
4865
4866
4867
4868
4869
4870
4871
4872
4873
4874
4875
4876
4877
4878
4879
4880
4881
4882
4883
4884
4885
4886
4887
4888
4889
4890
4891
4892
4893
4894
4895
4896
4897
4898
4899
4900
4901
4902
4903
4904
4905
4906
4907
4908
4909
4910
4911
4912
4913
4914
4915
4916
4917
4918
4919
4920
4921
4922
4923
4924
4925
4926
4927
4928
4929
4930
4931
4932
4933
4934
4935
4936
4937
4938
4939
4940
4941
4942
4943
4944
4945
4946
4947
4948
4949
4950
4951
4952
4953
4954
4955
4956
4957
4958
4959
4960
4961
4962
4963
4964
4965
4966
4967
4968
4969
4970
4971
4972
4973
4974
4975
4976
4977
4978
4979
4980
4981
4982
4983
4984
4985
4986
4987
4988
4989
4990
4991
4992
4993
4994
4995
4996
4997
4998
4999
5000
5001
5002
5003
5004
5005
5006
5007
5008
5009
5010
5011
5012
5013
5014
5015
5016
5017
5018
5019
5020
5021
5022
5023
5024
5025
5026
5027
5028
5029
5030
5031
5032
5033
5034
5035
5036
5037
5038
5039
5040
5041
5042
5043
5044
5045
5046
5047
5048
5049
5050
5051
5052
5053
5054
5055
5056
5057
5058
5059
5060
5061
5062
5063
5064
5065
5066
5067
5068
5069
5070
5071
5072
5073
5074
5075
5076
5077
5078
5079
5080
5081
5082
5083
5084
5085
5086
5087
5088
5089
5090
5091
5092
5093
5094
5095
5096
      MODULE LISTS
        INTEGER, ALLOCATABLE, DIMENSION(:,:)   :: KOUNT, MOUNT
        INTEGER, ALLOCATABLE, DIMENSION(:,:)   :: TTRANS
        INTEGER ISTRANS
      END MODULE LISTS

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

C
C New object types go here!
C

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

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

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

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

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

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

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

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

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

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

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



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

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

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

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

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

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


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

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


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

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

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