TRT2FT. USER,NICK,DLRS2. * * JOB TO PUNCH FORTRAN SOURCE * FROM MODIFY INPUT * RESULT IS "STANDARD FORTRAN-77" WHICH * CAN BE MOVED TO NON-CDC MACHINES AND * BUILT THERE - IN THEORY. * SETJSL,7777. SETASL,7777. MODIFY,P=0,C=PUNCH,F,Z./*NOSEQ/*CREATE INPUT REWIND,PUNCH. ROUTE,PUNCH,DC=PU. ~ PARAMS COMMON C C----------------------- C-- SUNDRY PARAMETERS -- C----------------------- C REAL GTHUGE, MINEPS, EPS, MINFLT, MAXFLT, TRTPI, DEGRAD PARAMETER( GTHUGE=1E7 ) PARAMETER( MINEPS=1E-5 ) PARAMETER( EPS=1E-7 ) PARAMETER( MINFLT=-1E20 ) PARAMETER( MAXFLT=1E20 ) PARAMETER( TRTPI=3.1415926 ) PARAMETER( DEGRAD=TRTPI/180.0 ) C INTEGER EYERAY, RFLRAY, TRNRAY PARAMETER( EYERAY=1 ) PARAMETER( RFLRAY=2 ) PARAMETER( TRNRAY=3 ) C INTEGER ENTER, LEAVE PARAMETER( ENTER=1 ) PARAMETER( LEAVE=2 ) C INTEGER NOACL, SUBACL PARAMETER( NOACL=0 ) PARAMETER( SUBACL=1 ) C INTEGER BKCNST, BKVERT, BKHORZ PARAMETER( BKCNST=0 ) PARAMETER( BKVERT=1 ) PARAMETER( BKHORZ=2 ) C INTEGER LGTCON, LGTCOS, LGTBRN PARAMETER( LGTCON=0 ) PARAMETER( LGTCOS=1 ) PARAMETER( LGTBRN=2 ) C INTEGER SPHPRM, TRIPRM PARAMETER( SPHPRM=1 ) PARAMETER( TRIPRM=2 ) C INTEGER DAXX, DAXY, DAXZ PARAMETER( DAXX=1 ) PARAMETER( DAXY=2 ) PARAMETER( DAXZ=3 ) ~ RAYSTK COMMON C C---- IN-MEMORY AND GLOBALLY ACCESSIBLE DATA STRUCTURES C C---------------------------- C-- RAY STACK COMMON BLOCK -- C---------------------------- C SIZE: DIMSTK + 2 WORDS = 514 C INTEGER MAXSTK, DIMSTK PARAMETER( MAXSTK=512 ) PARAMETER( DIMSTK=MAXSTK ) C INTEGER NSRAY, CRAY INTEGER RAY(DIMSTK) COMMON /RAYSTK/ NSRAY, CRAY, RAY ~ RAYLST COMMON C C--------------------------- C-- RAY LIST COMMON BLOCK -- C--------------------------- C SIZE: 13 * DIMRAY + 2 WORDS = 6658 C INTEGER MAXRAY, DIMRAY PARAMETER( MAXRAY=512 ) PARAMETER( DIMRAY=MAXRAY ) C INTEGER NLRAY, RAYNUM, RISECT(DIMRAY), RTYPE(DIMRAY) INTEGER RDEP(DIMRAY), RRAYID(DIMRAY) COMMON /RAYLSI/ NLRAY, RAYNUM, RISECT, RTYPE, RDEP, RRAYID C REAL ROX(DIMRAY), ROY(DIMRAY), ROZ(DIMRAY) REAL RDX(DIMRAY), RDY(DIMRAY), RDZ(DIMRAY) REAL RWGTR(DIMRAY), RWGTG(DIMRAY), RWGTB(DIMRAY) COMMON /RAYLSR/ ROX, ROY, ROZ, RDX, RDY, RDZ, RWGTR, RWGTG, RWGTB ~ INSLST COMMON C C------------------------------------ C-- INTERSECTION LIST COMMON BLOCK -- C------------------------------------ C SIZE: 13 * DIMINS + 1 WORDS = 6683 C INTEGER MAXINS, DIMINS, NEWINS, SHDINS PARAMETER( MAXINS=512 ) PARAMETER( DIMINS=MAXINS+2 ) PARAMETER( NEWINS=MAXINS+1 ) PARAMETER( SHDINS=MAXINS+2 ) C INTEGER NINS, IPRIM(DIMINS), IRRAY(DIMINS), ITRAY(DIMINS) INTEGER ITTYP(DIMINS) COMMON /INSLSI/ NINS, IPRIM, IRRAY, ITRAY, ITTYP C REAL IPX(DIMINS), IPY(DIMINS), IPZ(DIMINS) REAL INX(DIMINS), INY(DIMINS), INZ(DIMINS) REAL ICOLR(DIMINS), ICOLG(DIMINS), ICOLB(DIMINS) COMMON /INSLSR/ IPX, IPY, IPZ, INX, INY, INZ, ICOLR, ICOLG, ICOLB ~ PRMLST COMMON C C--------------------------------- C-- PRIMITIVE LIST COMMON BLOCK -- C--------------------------------- C SIZE: 5 * DIMPRM + 1 WORDS = 5121 C INTEGER MAXPRM, DIMPRM PARAMETER( MAXPRM=1024 ) PARAMETER( DIMPRM=MAXPRM ) C INTEGER NPRM, MATTER(DIMPRM), PACCEL(DIMPRM), PRAYID(DIMPRM) INTEGER PRMTYP(DIMPRM), PRMIDX(DIMPRM) COMMON /PRMLST/ NPRM, MATTER, PACCEL, PRAYID, PRMTYP, PRMIDX ~ SPHLST COMMON C C----------------- C-- SPHERE LIST -- C----------------- C SIZE: 4 * DIMSPH + 1 WORDS = 4097 C INTEGER MAXSPH, DIMSPH PARAMETER( MAXSPH=1024 ) PARAMETER( DIMSPH=MAXSPH ) C INTEGER NSPH COMMON /SPHLSI/ NSPH C REAL SOX(DIMSPH), SOY(DIMSPH), SOZ(DIMSPH), SRAD(DIMSPH) COMMON /SPHLSR/ SOX, SOY, SOZ, SRAD ~ TRILST COMMON C C------------------- C-- TRIANGLE LIST -- C------------------- C 23 * DIMTRI + 1 WORDS = 5889 C INTEGER MAXTRI, DIMTRI PARAMETER( MAXTRI=256 ) PARAMETER( DIMTRI=MAXTRI ) C INTEGER NTRI, DAX(DIMTRI) COMMON /TRILSI/ NTRI, DAX C REAL TX1(DIMTRI), TY1(DIMTRI), TZ1(DIMTRI) REAL TX2(DIMTRI), TY2(DIMTRI), TZ2(DIMTRI) REAL TX3(DIMTRI), TY3(DIMTRI), TZ3(DIMTRI) REAL TNX(DIMTRI), TNY(DIMTRI), TNZ(DIMTRI), TND(DIMTRI) REAL NVX1(DIMTRI), NVY1(DIMTRI), NVZ1(DIMTRI) REAL NVX2(DIMTRI), NVY2(DIMTRI), NVZ2(DIMTRI) REAL NVX3(DIMTRI), NVY3(DIMTRI), NVZ3(DIMTRI) COMMON /TRILSR/ TX1, TY1, TZ1, TX2, TY2, TZ2, TX3, TY3, TZ3, + TNX, TNY, TNZ, NVX1, NVY1, NVZ1, NVX2, NVY2, + NVZ2, NVX3, NVY3, NVZ3, TND ~ MTRLST COMMON C C-------------------------------- C-- MATERIAL LIST COMMON BLOCK -- C-------------------------------- C SIZE: 14 * DIMMTR + 1 WORDS = 449 C INTEGER MAXMTR, DIMMTR PARAMETER( MAXMTR=32 ) PARAMETER( DIMMTR=MAXMTR ) C INTEGER NMTR COMMON /MTRLSI/ NMTR C REAL MKR(DIMMTR), MKT(DIMMTR) REAL MCAR(DIMMTR), MCAG(DIMMTR), MCAB(DIMMTR) REAL MCDR(DIMMTR), MCDG(DIMMTR), MCDB(DIMMTR) REAL MCSR(DIMMTR), MCSG(DIMMTR), MCSB(DIMMTR) REAL META(DIMMTR), MGLS(DIMMTR), MRGH(DIMMTR) COMMON /MTRLSR/ MKR, MKT, MCAR, MCAG, MCAB, MCDR, MCDG, MCDB, + MCSR, MCSG, MCSB, META, MGLS, MRGH ~ LGTLST COMMON C C----------------------------- C-- LIGHT LIST COMMON BLOCK -- C----------------------------- C SIZE: 13 * DIMLGT + 4 WORDS = 420 C INTEGER MAXLGT, DIMLGT PARAMETER( MAXLGT=32 ) PARAMETER( DIMLGT=MAXLGT ) C INTEGER NLGT, LDIR(DIMLGT) COMMON /LGTLSI/ NLGT, LDIR C REAL LCAR, LCAG, LCAB, LOX(DIMLGT), LOY(DIMLGT), LOZ(DIMLGT) REAL LCLR(DIMLGT), LCLG(DIMLGT), LCLB(DIMLGT) REAL LDX(DIMLGT), LDY(DIMLGT), LDZ(DIMLGT) REAL LTGT(DIMLGT), LTG2(DIMLGT), LRAD(DIMLGT) COMMON /LGTLSR/ LCAR, LCAG, LCAB, LOX, LOY, LOZ, LCLR, LCLG, + LCLB, LDX, LDY, LDZ, LTGT, LTG2, LRAD ~ SCREEN COMMON C C------------------------- C-- SCREEN COMMON BLOCK -- C------------------------- C SIZE: 3 * DIMPXL + 15 WORDS = 12303 C INTEGER MAXPXL, DIMPXL PARAMETER( MAXPXL=4096 ) PARAMETER( DIMPXL=MAXPXL) C INTEGER SNX, SNY, SMAXDP, OS, ACCEL, SHADOW, BKT, CURX, CURY COMMON /SCREEI/ SNX, SNY, SMAXDP, OS, ACCEL, SHADOW, BKT, + CURX, CURY C REAL FPD, APER, FOCAL, FSTOP, MINWGT REAL COLR, COLG, COLB REAL BACKR, BACKG, BACKB REAL BKRI, BKRF, BKRD REAL BKGI, BKGF, BKGD REAL BKBI, BKBF, BKBD REAL LINER(DIMPXL), LINEG(DIMPXL), LINEB(DIMPXL) COMMON /SCREER/ FPD, APER, FOCAL, FSTOP, MINWGT, + COLR, COLG, COLB, + BACKR, BACKG, BACKB, + BKRI, BKRF, BKRD, + BKGI, BKGF, BKGD, + BKBI, BKBF, BKBD, + LINER, LINEG, LINEB ~ STATS COMMON C C----------------------------- C-- STATISTICS COMMON BLOCK -- C----------------------------- C SIZE: 8 WORDS. C INTEGER MAXDPS, MAXRYS, MAXISS, TOTRAY, TOTFEL, NACC, NOACC COMMON /STATSI/ MAXDPS, MAXRYS, MAXISS, TOTRAY, TOTFEL, NACC, + NOACC C REAL NSECTS COMMON /STATSR/ NSECTS ~ ACCEL COMMON C C------------------ C-- ACCELERATION -- C------------------ C FOR TYPICAL SCENES, THE "MAX CELLS" WOULD BE BEST AROUND 20 TO 40. C (MAYBE - MY TESTS SAY 10 IS GOOD). BUT THE MEMORY COST IS QUITE HIGH. C LIKEWISE, IT WOULD BE GOOD TO ALLOW LARGISH PRIMITIVES TO BE INDEXED, C BUT THE MAXIMUM POSSIBLE LIST SIZE GOES UP AS THE CUBE OF THE LARGEST C NUMBER OF CELLS ALLOWED IN ANY DIRECTION FOR THE PRIMITIVE TO BE ACCEPTED C FOR INDEXING. SO THE MEMORY COST IS (POTENTIALLY) VERY HIGH. IT WOULD BE C POSSIBLE TO USE THE SAME LIST FOR MORE THAN ONE GRID CELL IF IT IS IDENTICAL C TO THAT FOR ANOTHER CELL. WE DON'T TRY TO DETECT THAT, THOUGH. C IF THE RATIO OF THE BIGGEST DIMENSION RANGE TO THE SMALLEST EXCEEDS C MXDMRT, ADJUST THE RANGES FOR CALCULATING THE GRID PARAMETERS. THIS SEEMS C TO BE ALWAYS A GOOD IDEA (HENCE SET TO 1.0). C C SIZE: CELDIM + 2 * LSTDIM + 13 WORDS = 1000 + 2 * 8000 + 13 = 17013 WORDS. C INTEGER MXCELX, MXCELY, MXCELZ, DMCELX, DMCELY, DMCELZ INTEGER CELDIM, MXCELS, MXPCVR, MAXLST, LSTDIM REAL MXDMRT PARAMETER( MXCELX=10 ) PARAMETER( MXCELY=10 ) PARAMETER( MXCELZ=10 ) PARAMETER( DMCELX=MXCELX ) PARAMETER( DMCELY=MXCELY ) PARAMETER( DMCELZ=MXCELZ ) PARAMETER( CELDIM=DMCELX*DMCELY*DMCELZ ) PARAMETER( MXCELS=MXCELX*MXCELY*MXCELZ ) PARAMETER( MXPCVR=2 ) PARAMETER( MAXLST=2*2*2*MXCELS ) PARAMETER( LSTDIM=MAXLST ) PARAMETER( MXDMRT=1.0 ) C INTEGER NLIST, ANX, ANY, ANZ, CELLS(CELDIM) INTEGER LPRIM(LSTDIM), LNEXT(LSTDIM) COMMON /ACCELI/ NLIST, ANX, ANY, ANZ, CELLS, LPRIM, LNEXT C REAL AXL, AXH, AYL, AYH, AZL, AZH, ADXS, ADYS, ADZS COMMON /ACCELR/ AXL, AXH, AYL, AYH, AZL, AZH, ADXS, ADYS, ADZS ~ GRID COMMON C C-------------------- C-- GRID TRAVERSAL -- C-------------------- C SIZE: 10 WORDS. C INTEGER VNINST, VASECT, VURAYI COMMON /VSCACI/ VNINST, VASECT, VURAYI C REAL VOX, VOY, VOZ, VDX, VDY, VDZ, VTMIN COMMON /VSCACR/ VOX, VOY, VOZ, VDX, VDY, VDZ, VTMIN ~ TOTSIZE COMMON C C TOTAL "COMMON BLOCK" DATA SIZE: C 514 RAY STACK C 6658 RAY LIST C 6683 INTERSECTIONS C 5121 PRIMITIVES C 4097 SPHERES C 5121 TRIANGLES C 449 MATERIALS C 420 LIGHTS C 12303 SCREEN C 8 STATS C 17013 ACCELERATION C 10 TRAVERSAL C ***** C 58397 WORDS (OF 131008 MAX AVAILABLE.) C ***** C ~ TRT2 PROGRAM TRT2(INPUT,OUTPUT,IMAGE, + TAPE5=INPUT,TAPE6=OUTPUT,TAPE7=IMAGE) C********************************************************************** C 1 2 3 4 5 6 7 V 8 C2345678901234567890123456789012345678901234567890123456789012345678901234567890 C C ******************************** C TOY RAY TRACER MARK 2 C TRT2 - FORTRAN-77 IMPLEMENTATION C ******************************** C C A SIMPLE RAY TRACER RENDERER. C THE MAIN "PURPOSE" OF THIS CODE IS TO SERVE AS A NON-TRIVIAL AND C ARGUABLY "FUN" APPLICATION TO RUN IN "RETRO-COMPUTING" ENVIRONMENTS. C OK - SO IT REALLY DOESN'T HAVE MUCH OF A PURPOSE! C C A PROTOTYPE WAS ACTUALLY WRITTEN AND DEBUGGED IN C++ (WRITTEN SO IT C COULD BE TRANSLATED PRETTY MUCH DIRECTLY TO FORTRAN-77). IT SHOULD C BE STRAIGHTFORWARD (IF TIME CONSUMING) TO MOVE IT TO OTHER LANGUAGES C - ALTHOUGH THAT WOULDN'T EXACTLY SHOW THE BEST FEATURES OF THOSE C LANUGAGES. C C IT IS ALSO AN EXPERIMENT IN HOW SIMPLE AND SMALL A "USEFUL" C RENDERER CAN BE - IN TERMS OF LINES OF CODE, ANYWAY. C C THE GOALS WERE: C 1) PORTABILITY C 2) EASE OF RECODING IN DIFFERENT LANGUAGES C 3) SMALL MEMORY FOOTPRINT C 4) STILL CAPABLE OF MAKING PRETTY PICTURES C 5) EASY EXTENSIBILITY: C A) NEW GEOMETRIC PRIMITIVE TYPES SHOULD BE STRAIGHTFORWARD C TO ADD C B) TEXTURE MAPPING, ETC., WOULD BE STRAIGHTFORWARD TO C ADD ON MACHINES WITH ENOUGH MEMORY TO MAKE IT FEASIBLE. C C CURRENT FEATURES: C 1) SPHERE PRIMITIVES (PROVING AGAIN THAT RAY TRACING IS MOSTLY C A LOAD OF BALLS). C 2) TRIANGLE PRIMITIVES WITH OPTIONAL VERTEX NORMALS C (THESE WILL BE CORRECTLY INTERPOLATED) C 3) STOCHASTIC ANTI-ALIASING (BASIC DISTRIBUTED RAY TRACING) C 4) OPTIONAL SHADOWS, OPTIONALLY APPROPRIATELY "SOFT EDGED" C 5) SPECULARLY REFLECTIVE MATERIALS (CORRECT SPECULAR REFLECTION) C 6) TRANSPARENT MATERIALS (CORRECT REFRACTION) C 7) SIMPLE SURFACE ROIUGHNESS FACTOR FOR MATERIALS. C 8) ACCELERATED RAY INTERSECTION TESTS USING A SPACE SUBDIVISION C (VOXEL) SCHEME C 9) IMAGES OUTPUT AS TEXT FILES (WELL - IT IS PORTABLE!). C 10) DIRECTIONAL LIGHTS (RAISED COSINE AND BARN DOOR TYPES) C C HAVING SPACE SUBDIVISION ACCELERATION MEANS IT ISN'T ALL C THAT SIMPLE! HOWEVER, THIS CAN MAKE A REALLY BIG DIFFERENCE C WITH LOTS OF PRIMITIVES. WELL - ASSUMING THEY CAN BE C "INDEXED" (SEE CODE). TO SAVE MEMORY, PRIMITIVES THAT C ARE TOO BIG CANNOT BE "ACCELERATED". C C BECAUSE THIS CODE MUST WORK WITH STANDARD FORTRAN-77, C RECURSION IS EXPLICITLY CODED USING A RAY STACK WITHOUT MAKING C ANY ACTUAL RECURSIVE CALLS. ALSO, ALL STORAGE MUST BE C STATICALLY ALLOCATED AT COMPILE TIME. C C AS ALWAYS WITH FLOATING POINT BASED CODE, THE PRECISION WITH WHICH C CALCULATIONS ARE DONE IS IMPORTANT. TESTS FOR APPROXIMATE EQUIVALENCE C AND "NUDGING" THINGS TO AVOID DEGENERATE CASES IS UNAVOIDABLE. PARTS OF C THE CODE WHICH MAY SUFFER FROM "PRECISION" ISSUES ARE HIGHLIGHTED WITH C THE COMMENT: PRECISION! THE PARAMETERS EPS AND MINEPS CAN TUNE SOME OF C THIS. NOTE THAT 32BIT FLOATING POINT IS *NOT* GOOD ENOUGH FOR THE GRID C TRAVERSAL CODE. IT WORKS PERFECTLY WITH 64BIT (WELL - PRETTY EXHAUSTIVE C TESTING HAS FAILED TO MAKE IT GO WRONG) BUT *NOT* WITH 32 BIT. WHAT C WILL HAPPEN WITH CDC FORMAT 60 BIT? SO FAR, SO GOOD. C C THE PRIMARY GOAL OF THE FORTRAN VERSION OF THIS CODE IS TO RUN ON C EMULATED CDC MAINFRAMES. THE TOTAL SIZE OF THE PROGRAM MUST BE KEPT C (WELL) BELOW 131000 60 BIT WORDS. C C AUTHOR: NICK GLAZZARD C C++ EXPERIMENTS AND PROTOTYPE: NOVEMBER 2003 C F-77 IMPLEMENTATION: JULY/AUGUST 2005 C TRT (MARK 1): SUMMER 1984 (LOST LONG AGO) C C---------------------------------------------------------------------- C MAIN C---------------------------------------------------------------------- C IMPLICIT CHARACTER*1 (A-Z) *CALL PARAMS *CALL SCREEN *CALL STATS *CALL PRMLST *CALL SPHLST *CALL TRILST *CALL MTRLST *CALL LGTLST *CALL RAYLST *CALL INSLST INTEGER Y, RDOK C WRITE(6,100) 100 FORMAT(1H1,'TOY RAY TRACER MARK2 V0.1 FORTRAN-77' ) WRITE(6,101) 101 FORMAT(1X,'************************************' ) C C-- INITIALIZE THE SCENE DATABASE C CALL INITDB C C-- INITIALIZE STATISTICS C CALL INITST C C-- READ THE SCENE FROM "INPUT" AND BUILD SCENE DATABASE C CALL READDB( RDOK ) IF( RDOK .EQ. 0 )THEN WRITE(6,102) 102 FORMAT(1X,'ERROR: FAILED TO READ SCENE DATABASE.') STOP ENDIF C C-- PRINT INFORMATION C WRITE(6,103) 103 FORMAT(1X,'SCREEN SPECIFICATION:') IF( ACCEL .EQ. NOACL )THEN WRITE(6,104) 104 FORMAT(6X,' ACCEL: NONE') ELSE WRITE(6,105) 105 FORMAT(6X,' ACCEL: VOXEL GRID SUBDIVISION') ENDIF IF( SHADOW .EQ. 0 )THEN WRITE(6,106) 106 FORMAT(6X,' SHADOWS: NO') ELSE WRITE(6,107)SHADOW 107 FORMAT(6X,' SHADOWS: YES,',I5,' FEELERS') ENDIF WRITE(6,108)SNX, MAXPXL 108 FORMAT(6X,' OUTPUT PIXELS:',I5,' MAX:',I5) WRITE(6,109)SNY, MAXPXL 109 FORMAT(6X,' OUTPUT LINES:',I5,' MAX:',I5) WRITE(6,110)OS 110 FORMAT(6X,' OVERSAMPLE:',I3,' (SQUARED)') WRITE(6,111)APER 111 FORMAT(6X,' APERTURE:',F9.4) WRITE(6,112)FPD 112 FORMAT(6X,' FOCAL PLANE DIST:',F9.4) WRITE(6,113)FOCAL 113 FORMAT(6X,'LENS FOCAL LENGTH:',F9.4) WRITE(6,114)FSTOP 114 FORMAT(6X,' LENS F-STOP:',F9.4) IF( BKT .EQ.BKCNST )THEN WRITE(6,115)BKRI, BKGI, BKBI 115 FORMAT(6X,'BACKGROUND: CONSTANT',F9.4,',',F9.4,',',F9.4) ELSE IF( BKT .EQ. BKVERT )THEN WRITE(6,116)BKRI, BKGI, BKBI, BKRF, BKGF, BKBF 116 FORMAT(6X,'BACKGROUND: VERTGRAD',F9.4,',',F9.4,',',F9.4, + ' TO ',F9.4,',',F9.4,',',F9.4) ELSE WRITE(6,117)BKRI, BKGI, BKBI, BKRF, BKGF, BKBF 117 FORMAT(6X,'BACKGROUND: HORZGRAD',F9.4,',',F9.4,',',F9.4, + ' TO ',F9.4,',',F9.4,',',F9.4) ENDIF C WRITE(6,118) 118 FORMAT(1X,'CONSTRAINTS:') WRITE(6,119)SMAXDP 119 FORMAT(6X,'MAX RAY TREE DEPTH: ',I5) WRITE(6,120)MINWGT 120 FORMAT(6X,'MIN SIG RAY WEIGHT: ',F9.4) C WRITE(6,121) 121 FORMAT(1X,'DATABASE STATISTICS:') WRITE(6,122)NPRM, MAXPRM 122 FORMAT(6X,'PRIMITIVES:',I5,' MAX:',I5) WRITE(6,123)NSPH, MAXSPH 123 FORMAT(6X,' SPHERES:',I5,' MAX:',I5) WRITE(6,124)NTRI, MAXTRI 124 FORMAT(6X,' TRIANGLES:',I5,' MAX:',I5) WRITE(6,125)NMTR, MAXMTR 125 FORMAT(6X,' MATERIALS:',I5,' MAX:',I5) WRITE(6,126)NLGT, MAXLGT 126 FORMAT(6X,' LIGHTS:',I5,' MAX:',I5) C C-- IF ACCELERATING, INIT THE ACCEL STRUCTURES C IF( ACCEL .EQ. SUBACL )CALL STACEL C C-- STEP OVER ALL THE LINES OF THE OUTPUT IMAGE C-- TRACE ALL RAYS NEEDED TO GENERATE THE LINE C-- SHADE THE LINE C-- OUTPUT THE LINE C WRITE(7,200)SNX, SNY 200 FORMAT(1X,2(I5)) DO 1 Y=1,SNY IF( ( ( Y / 10 ) * 10 ) .EQ. Y )THEN CALL DISPLA( '... LINE: ', Y ) ENDIF CURY = Y CALL INITLN CALL TRCLN( Y ) CALL OUTLN( Y ) 1 CONTINUE C C-- PRINT STATISTICS C WRITE(6,127) 127 FORMAT(1X,'EXECUTION STATISTICS:') WRITE(6,128)MAXRYS, MAXRAY 128 FORMAT(6X,' MAX RAYS IN TREE:',I5,' MAX:',I5) WRITE(6,129)MAXISS, MAXINS 129 FORMAT(6X,'MAX INTERSECTS PER RAY:',I5,' MAX:',I5) WRITE(6,130)NSECTS 130 FORMAT(6X,' TOTAL ISECT TESTS:',F11.0) WRITE(6,131)TOTRAY 131 FORMAT(6X,' TOTAL RAYS TRACED:',I12) WRITE(6,132)TOTFEL 132 FORMAT(6X,' TOTAL LIGHT FEELERS:',I12) C STOP END C C SUBROUTINE OUTLN( Y ) IMPLICIT CHARACTER*1 (A-Z) INTEGER Y C********************************************************************** C OUTPUT A LINE OF THE IMAGE AS FORMATTED TEXT C********************************************************************** *CALL SCREEN INTEGER TREC, TREM, IREC, S, E, I, J, K INTEGER PIXI(24) C C-- ENFORCE LEGAL COLOUR RANGE NOW. C DO 5 I=1,SNX IF( LINER(I) .LT. 0.0 )LINER(I) = 0.0 IF( LINER(I) .GT. 1.0 )LINER(I) = 1.0 IF( LINEG(I) .LT. 0.0 )LINEG(I) = 0.0 IF( LINEG(I) .GT. 1.0 )LINEG(I) = 1.0 IF( LINEB(I) .LT. 0.0 )LINEB(I) = 0.0 IF( LINEB(I) .GT. 1.0 )LINEB(I) = 1.0 5 CONTINUE C C-- OUTPUT 8 RGB VALUES PER FULL TEXT LINE (72 CHARS/LINE) C-- ORDER RGBRGB... C TREC = SNX / 8 TREM = SNX - ( TREC * 8 ) S = 1 E = 8 C C-- OUTPUT ANY FULL TEXT LINES C DO 1 IREC=1,TREC K = 1 DO 2 J=S,E PIXI(K) = INT( LINER(J) * 255.9 ) K = K + 1 PIXI(K) = INT( LINEG(J) * 255.9 ) K = K + 1 PIXI(K) = INT( LINEB(J) * 255.9 ) K = K + 1 2 CONTINUE WRITE(7,100)(PIXI(I),I=1,24) 100 FORMAT(1X,24(I3)) S = S + 8 E = E + 8 1 CONTINUE C C-- OUTPUT ANY PARTIAL TEXT LINE C IF( TREM .NE. 0 )THEN E = S + TREM - 1 K = 1 DO 3 J=S,E PIXI(K) = INT( LINER(J) * 255.9 ) K = K + 1 PIXI(K) = INT( LINEG(J) * 255.9 ) K = K + 1 PIXI(K) = INT( LINEB(J) * 255.9 ) K = K + 1 3 CONTINUE WRITE(7,100)(PIXI(I),I=1,(3*TREM)) ENDIF C RETURN END C C---------------------------------------------------------------------- C RAY TRACING ROUTINES C---------------------------------------------------------------------- C SUBROUTINE INITLN IMPLICIT CHARACTER*1 (A-Z) C********************************************************************** C INITIALIZE THE OUTPUT LINE C********************************************************************** *CALL SCREEN C INTEGER I DO 1 I=1,SNX,1 LINER(I) = 0.0 LINEG(I) = 0.0 LINEB(I) = 0.0 1 CONTINUE RETURN END C C SUBROUTINE INITST IMPLICIT CHARACTER*1 (A-Z) C********************************************************************** C INITIALIZE THE STATISTICS C********************************************************************** *CALL STATS C MAXDPS = 0 MAXRYS = 0 MAXISS = 0 TOTRAY = 0 TOTFEL = 0 NACC = 0 NOACC = 0 NSECTS = 0.0 RETURN END C C SUBROUTINE TRCLN( Y ) IMPLICIT CHARACTER*1 (A-Z) INTEGER Y C********************************************************************** C TRACE ALL RAYS NEEDED TO CREATE LINE Y C********************************************************************** *CALL PARAMS *CALL SCREEN *CALL RAYLST *CALL RAYSTK *CALL INSLST *CALL STATS C C-- FUNCTION "PROTOTYPES" C REAL RANHLF C C-- LOCAL VARIABLES C REAL LEFT, RIGHT, BOTTOM, TOP, DP, DPX, DPY, RNS, RMAG REAL FY, FX, DS, FYS, FYSJ, FXS, FXSJ, DSX, DSY, LRAD REAL ETHETA, ERAD, EX, EY INTEGER X, SUBX, SUBY, NS, OS2, UOS, IDUM C C-- CALCULATE THE PIXEL CENTER COORDINATES AT THE EDGES OF THE APERTURE C DP = APER / REAL(SNX) LEFT = ( -APER + DP ) * 0.5 RIGHT = ( APER - DP ) * 0.5 BOTTOM = LEFT * REAL(SNY) / REAL(SNX) TOP = RIGHT * REAL(SNY) / REAL(SNX) C C-- CALCULATE THE RECIPROCAL OF THE MAGNIFICATION OF THE LENS C IF( FPD .LT. ( FOCAL + 0.5 ) )FPD = FOCAL + 0.5 RMAG = ( FPD - FOCAL ) / FOCAL C C-- SCALE APERTURE EDGE COORDINATES TO FOCAL PLANE EDGE COORDINATES C LEFT = LEFT * RMAG RIGHT = RIGHT * RMAG BOTTOM = BOTTOM * RMAG TOP = TOP * RMAG C C-- CALCULATE THE EFFECTIVE LENS RADIUS C LRAD = 0.5 * FOCAL / FSTOP C C-- CALCULATE CONSTANTS FOR SUBSAMPLING C IF( OS .EQ. 0 )THEN DS = 0.0 NS = 1 RNS = 1.0 OS2 = 0 UOS = 1 ELSE DS = 1.0 / OS NS = OS * OS RNS = 1.0 / NS OS2 = OS / 2 UOS = OS ENDIF C C-- CALCULATE THE X AND Y STEPS ON THE FOCAL PLANE BETWEEN PIXEL CENTERS. C DPX = ( RIGHT - LEFT ) / REAL( SNX - 1 ) DPY = ( TOP - BOTTOM ) / REAL( SNY - 1 ) DSX = DS * DPX DSY = DS * DPY C C-- CALCULATE Y COORDINATE OF CURRENT LINE CENTER C FY = BOTTOM + ( Y - 1 ) * DPY C C-- STEP OVER THE PIXELS IN THE LINE. C DO 1 X=1,SNX CURX = X FX = LEFT + ( X - 1 ) * DPX C C-- STEP OVER Y SUBSAMPLES - AND JITTER THEM C DO 2 SUBY=1,UOS FYS = ( SUBY - OS2 - 1 ) * DSY + FY FYSJ = FYS + RANHLF(IDUM) * DSY C C-- STEP OVER X SUBSAMPLES - AND JITTER THEM. C DO 3 SUBX=1,UOS FXS = ( SUBX - OS2 - 1 ) * DSX + FX FXSJ = FXS + RANHLF(IDUM) * DSX C C-- BEGIN TRACING RAYS. THERE ARE INITIALLY NO RAYS. AND NO INTERSECTIONS. C NLRAY = 0 NSRAY = 0 NINS = 0 COLR = 0.0 COLG = 0.0 COLB = 0.0 C C-- GENERATE A SCREEN RAY FROM A RANDOMLY CHOSEN POINT ON THE LENS (EX,EY,0) C-- THROUGH THE FOCAL PLANE POSITION (FX,FY,SCREEN.FPD) C NLRAY = NLRAY + 1 RAYNUM = RAYNUM + 1 ETHETA = RANHLF(IDUM) * TRTPI ERAD = ( RANHLF(IDUM) + 0.5 ) * LRAD EX = ERAD * COS( ETHETA ) EY = ERAD * SIN( ETHETA ) RRAYID(NLRAY) = RAYNUM RTYPE(NLRAY) = EYERAY RISECT(NLRAY) = 0 RWGTR(NLRAY) = 1.0 RWGTG(NLRAY) = 1.0 RWGTB(NLRAY) = 1.0 RDEP(NLRAY) = 1 ROX(NLRAY) = EX ROY(NLRAY) = EY ROZ(NLRAY) = 0.0 RDX(NLRAY) = FXSJ - EX RDY(NLRAY) = FYSJ - EY RDZ(NLRAY) = FPD CALL NRMVEC( RDX(NLRAY), RDY(NLRAY), RDZ(NLRAY) ) TOTRAY = TOTRAY + 1 C C-- PUSH THE RAY ONTO THE RAY STACK. C NSRAY = NSRAY + 1 RAY(NSRAY) = NLRAY IF( NLRAY .GT. MAXRYS )MAXRYS = NLRAY C C-- WHILE THERE ARE RAYS ON THE RAY STACK, POP A RAY OFF C 4 CONTINUE IF( NSRAY .EQ. 0 )GOTO 5 CRAY = RAY(NSRAY) NSRAY = NSRAY - 1 C C-- FIND INTERSECTIONS OF THIS CURRENT RAY WITH PRIMITIVES. C-- WHEN THERE IS AN INTERSECTION, CALCULATE A SHADE AT THE INTERSECTION POINT C-- AND ADD THIS CONTRIBUTION TO THE PIXEL COLOUR. C-- ADD NEW REFLECTED AND TRANSMITTED RAYS TO THE RAY STACK IF APPROPRIATE. C CALL TRACE GOTO 4 5 CONTINUE C C--- ADD IN THE NEW SAMPLE COLOUR TO THE PIXEL C LINER(X) = LINER(X) + COLR LINEG(X) = LINEG(X) + COLG LINEB(X) = LINEB(X) + COLB C C-- END LOOP OVER SUBSAMPLES IN X C 3 CONTINUE C C-- END LOOP OVER SUBSAMPLES IN Y C 2 CONTINUE C C-- NORMALIZE THE RESULT COLOUR FOR THE OUTPUT PIXEL C LINER(X) = LINER(X) * RNS LINEG(X) = LINEG(X) * RNS LINEB(X) = LINEB(X) * RNS C C-- END LOOP OVER PIXELS OF OUTPUT LINE C 1 CONTINUE RETURN END C C SUBROUTINE TRACE IMPLICIT CHARACTER*1 (A-Z) C********************************************************************** C INTERSECT THE "CURRENT RAY" (RAYSTK.CRAY) WITH ALL THE PRIMITIVES IN THE SCENE. C IF THERE IS AN INTERSECTION RECORD THE INTERSECTION IN THE INTERSECTION LIST. C SET THE CURRENT RAY INTERSECTION INDEX TO THE NEW INTERSECTION. C FIND A SHADE AT THE INTERSECTION AND ADD THAT CONTRIBUTION TO THE SAMPLE COLOUR. C SPAWN NEW REFLECTED AND TRANSMITTED RAYS IF APPROPRIATE. PUSH THEM ON THE RAY STACK. C********************************************************************** *CALL PARAMS *CALL RAYLST *CALL RAYSTK *CALL INSLST *CALL PRMLST *CALL MTRLST *CALL SCREEN *CALL STATS C C-- LOCAL VARIABLES C REAL OX, OY, OZ REAL DX, DY, DZ REAL NDOTI, ETA, RETA, COS1, COS22, K2, SWGT, TVAL REAL RWGT, LRWGTR, LRWGTG, LRWGTB REAL TWGT, TWGTR, TWGTG, TWGTB REAL CWGTR, CWGTG, CWGTB INTEGER NINST, PRIM, MATL, ASECT, DEP, TTYP, IDUM C C-- GET USEFUL THINGS ABOUT THE CURRENT RAY. C OX = ROX(CRAY) OY = ROY(CRAY) OZ = ROZ(CRAY) DX = RDX(CRAY) DY = RDY(CRAY) DZ = RDZ(CRAY) CWGTR = RWGTR(CRAY) CWGTG = RWGTG(CRAY) CWGTB = RWGTB(CRAY) DEP = RDEP(CRAY) IF( DEP .GT. MAXDPS )MAXDPS = DEP C C-- TEST FOR INTERSECTIONS. IF THERE ARE ANY, WE NEED THE NEAREST. C-- IF THERE WAS AN INTERSECTION, ITS POSITION AND SURFACE NORMAL WILL BE IN C-- INSLST.XXX(INSLST.NINS+1) AFTER THIS. C NINST = NINS + 1 IF( ACCEL .EQ. NOACL )THEN CALL GENSCT( OX, OY, OZ, DX, DY, DZ, NINST, ASECT, TVAL, 1 ) ELSE CALL CELSCT( OX, OY, OZ, DX, DY, DZ, NINST, ASECT, TVAL, 1 ) ENDIF C C-- IF THERE ACTUALLY WAS AN INTERSECTION, UPDATE THE INTERSECTION COUNT. C-- CREATE NEW REFLECTED AND TRANSMITTED RAYS IF APPROPRIATE. C IF( ASECT .NE. 0 )THEN NINS = NINST IF( NINST .GT. MAXISS )MAXISS = NINST PRIM = IPRIM(NINST) MATL = MATTER(PRIM) LRWGTR = MKR(MATL) * MCDR(MATL) * CWGTR LRWGTG = MKR(MATL) * MCDG(MATL) * CWGTG LRWGTB = MKR(MATL) * MCDB(MATL) * CWGTB TWGTR = MKT(MATL) * MCDR(MATL) * CWGTR TWGTG = MKT(MATL) * MCDG(MATL) * CWGTG TWGTB = MKT(MATL) * MCDB(MATL) * CWGTB TTYP = ITTYP(NINST) ETA = META(MATL) C C-- USE MAXIMUM OF COLOUR COMPONENT WEIGHTS TO DECIDE IF TO TERMINATE C RWGT = LRWGTR IF( LRWGTG .GT. RWGT )RWGT = LRWGTG IF( LRWGTB .GT. RWGT )RWGT = LRWGTB TWGT = TWGTR IF( TWGTG .GT. TWGT )TWGT = TWGTG IF( TWGTB .GT. TWGT )TWGT = TWGTB C C-- FLIP THE NORMAL? C IF( ( DX * INX(NINST) + DY * INY(NINST) + DZ * INZ(NINST) ) + .GT. 0.0 )THEN INX(NINST) = -INX(NINST) INY(NINST) = -INY(NINST) INZ(NINST) = -INZ(NINST) ENDIF C C-- CALCULATE THE COLOUR AT THE INTERSECTION. C CALL SHADE( NINST, MATL ) SWGT = 1.0 - MKR(MATL) - MKT(MATL) IF( SWGT .LT. 0.0 )SWGT = 0.0 COLR = COLR + ( ICOLR(NINST) * SWGT * CWGTR ) COLG = COLG + ( ICOLG(NINST) * SWGT * CWGTG ) COLB = COLB + ( ICOLB(NINST) * SWGT * CWGTB ) C C-- SPAWN NEW RAYS (RECURSE, EFFECTIVELY) IF WE ARE NOT TOO DEEP. C IF( DEP .LT. SMAXDP )THEN DEP = DEP + 1 C C-- NEW REFLECTED RAY C IF( RWGT .GT. MINWGT )THEN NDOTI = 2 * ( INX(NINST) * DX + INY(NINST) * DY + + INZ(NINST) * DZ ) NLRAY = NLRAY + 1 RAYNUM = RAYNUM + 1 RRAYID(NLRAY) = RAYNUM RTYPE(NLRAY) = RFLRAY RISECT(NLRAY) = 0 RWGTR(NLRAY) = LRWGTR RWGTG(NLRAY) = LRWGTG RWGTB(NLRAY) = LRWGTB RDEP(NLRAY) = DEP ROX(NLRAY) = IPX(NINST) ROY(NLRAY) = IPY(NINST) ROZ(NLRAY) = IPZ(NINST) RDX(NLRAY) = DX - NDOTI * INX(NINST) RDY(NLRAY) = DY - NDOTI * INY(NINST) RDZ(NLRAY) = DZ - NDOTI * INZ(NINST) CALL NRMVEC( RDX(NLRAY), RDY(NLRAY), RDZ(NLRAY) ) IRRAY(NINST) = NLRAY NSRAY = NSRAY + 1 RAY(NSRAY) = NLRAY TOTRAY = TOTRAY + 1 IF( NLRAY .GT. MAXRYS )MAXRYS = NLRAY ENDIF C C-- NEW TRANSMITTED RAY C IF( TWGT .GT. MINWGT )THEN IF( TTYP .EQ. ENTER )THEN RETA = ETA ELSE RETA = 1.0 / ETA ENDIF COS1 = -( INX(NINST) * DX + INY(NINST) * DY + + INZ(NINST) * DZ ) COS22 = 1.0 - RETA * RETA * ( 1.0 - COS1 * COS1 ) IF( COS22 .GT. 0.0 )THEN K2 = RETA * COS1 - SQRT( COS22 ) NLRAY = NLRAY + 1 RAYNUM = RAYNUM + 1 RRAYID(NLRAY) = RAYNUM RTYPE(NLRAY) = TRNRAY RISECT(NLRAY) = 0 RWGTR(NLRAY) = TWGTR RWGTG(NLRAY) = TWGTG RWGTB(NLRAY) = TWGTB RDEP(NLRAY) = DEP ROX(NLRAY) = IPX(NINST) ROY(NLRAY) = IPY(NINST) ROZ(NLRAY) = IPZ(NINST) RDX(NLRAY) = RETA * DX + K2 * INX(NINST) RDY(NLRAY) = RETA * DY + K2 * INY(NINST) RDZ(NLRAY) = RETA * DZ + K2 * INZ(NINST) CALL NRMVEC( RDX(NLRAY), RDY(NLRAY), RDZ(NLRAY) ) ITRAY(NINST) = NLRAY NSRAY = NSRAY + 1 RAY(NSRAY) = NLRAY TOTRAY = TOTRAY + 1 IF( NLRAY .GT. MAXRYS )MAXRYS = NLRAY ENDIF ENDIF C C-- END OF LESS THAN MAXIMUM RAY STACK DEPTH C ENDIF C C-- IF THIS RAY DOESN'T HIT ANYTHING, ADD IN THE BACKGROUND COLOUR. C ELSE CALL CALCBK COLR = COLR + BACKR * CWGTR COLG = COLG + BACKG * CWGTG COLB = COLB + BACKB * CWGTB ENDIF C RETURN END C C SUBROUTINE CALCBK IMPLICIT CHARACTER*1 (A-Z) C********************************************************************** C CALCULATE THE BACKGROUND COLOUR FOR GRADED BACKGROUNDS C********************************************************************** *CALL PARAMS *CALL SCREEN C INTEGER CX, CY CX = CURX - 1 CY = CURY - 1 IF( BKT .EQ. BKCNST )THEN BACKR = BKRI BACKG = BKGI BACKB = BKBI ELSE IF( BKT .EQ. BKVERT )THEN BACKR = BKRI + BKRD * CY BACKG = BKGI + BKGD * CY BACKB = BKBI + BKBD * CY ELSE IF( BKT .EQ. BKHORZ )THEN BACKR = BKRI + BKRD * CX BACKG = BKGI + BKGD * CX BACKB = BKBI + BKBD * CX ENDIF IF( BACKR .LT. 0.0 )THEN BACKR = 0.0 ELSE IF( BACKR .GT. 1.0 )THEN BACKR = 1.0 ENDIF IF( BACKG .LT. 0.0 )THEN BACKG = 0.0 ELSE IF( BACKG .GT. 1.0 )THEN BACKG = 1.0 ENDIF IF( BACKB .LT. 0.0 )THEN BACKB = 0.0 ELSE IF( BACKB .GT. 1.0 )THEN BACKB = 1.0 ENDIF C RETURN END C C---------------------------------------------------------------------- C SHADING ROUTINES C---------------------------------------------------------------------- C SUBROUTINE SHADE( I, MATL ) IMPLICIT CHARACTER*1 (A-Z) INTEGER I, MATL C********************************************************************** C FIND THE COLOUR AT INTERSECTION I WHERE MATERIAL MATL IS TO BE FOUND. C********************************************************************** *CALL PARAMS *CALL INSLST *CALL SCREEN *CALL LGTLST *CALL MTRLST C REAL RANHLF C INTEGER L, F REAL DX, DY, DZ, D, S, VX, VY, VZ, HX, HY, HZ, CR, CG, CB REAL LATTR, LATTG, LATTB, P REAL NX, NY, NZ, A, LXF, LYF, LZF, LATTRF, LATTGF, LATTBF C INTEGER IDUM DATA IDUM/12345/ C C-- INITIALIZE INTERSECTION COLOUR TO AMBIENT TERM. C ICOLR(I) = MCAR(MATL) * LCAR ICOLG(I) = MCAG(MATL) * LCAR ICOLB(I) = MCAB(MATL) * LCAR C C-- STEP OVER LIGHTS. FIND CONTRIBUTION TO INTERSECTION COLOUR FOR EACH LIGHT. C DO 1 L=1,NLGT C C-- SEE IF THIS LIGHT IS PARTIALLY SHADOWED. C IF( SHADOW .EQ. 0 )THEN LATTR = 1.0 LATTG = 1.0 LATTB = 1.0 ELSE IF( ( SHADOW .EQ. 1 ) .OR. ( LRAD(L) .LT. EPS ) )THEN C-- 1 FEELER RAY = POINT LIGHT SOURCE. CALL LTFEEL( IPRIM(I), IPX(I), IPY(I), IPZ(I), + LOX(L), LOY(L), LOZ(L), LATTR, LATTG, LATTB ) ELSE C-- MULTIPLE FEELER RAYS = UNIFORM SPHERE OF LIGHT. LATTR = 0.0 LATTG = 0.0 LATTB = 0.0 DO 2 F=1,SHADOW LXF = LOX(L) + RANHLF(IDUM) * LRAD(L) LYF = LOY(L) + RANHLF(IDUM) * LRAD(L) LZF = LOZ(L) + RANHLF(IDUM) * LRAD(L) CALL LTFEEL( IPRIM(I), IPX(I), IPY(I), IPZ(I), + LXF, LYF, LZF, LATTRF, LATTGF, LATTBF ) LATTR = LATTR + LATTRF LATTG = LATTG + LATTGF LATTB = LATTB + LATTBF 2 CONTINUE LATTR = LATTR / SHADOW LATTG = LATTG / SHADOW LATTB = LATTB / SHADOW ENDIF C C-- FIND UNIT VECTOR FROM SURFACE INTERSECTION TO LIGHT. C DX = LOX(L) - IPX(I) DY = LOY(L) - IPY(I) DZ = LOZ(L) - IPZ(I) CALL NRMVEC( DX, DY, DZ ) C C-- IF THE LIGHT IS DIRECTIONAL, MODIFY LIGHT ATTENUATION BASED ON C-- ANGLE BTW LIGHT-SURFACE VECTOR AND LIGHT DIRECTION. C IF( LDIR(L) .NE. LGTCON )THEN P = DX * LDX(L) + DY * LDY(L) + DZ * LDZ(L) IF( P .LT. 0.0 )THEN P = 0.0 ELSE IF( LDIR(L) .EQ. LGTCOS )THEN C-- RAISED COSINE SORT ... P = P ** LTGT(L) ELSE IF( LDIR(L) .EQ. LGTBRN )THEN C-- BARN DOOR SORT ... A = ACOS( P ) IF( A .LT. LTGT(L) )THEN P = 1.0 ELSE IF( A .GT. LTG2(L) )THEN P = 0.0 ELSE P = ( ( A - LTGT(L) ) / ( LTG2(L) - LTGT(L) ) ) C-- NICE SMOOTH HERMITE ... P = 1.0 - ( P * P * ( 3 - 2 * P ) ) ENDIF ENDIF ENDIF LATTR = LATTR * P LATTG = LATTG * P LATTB = LATTB * P ENDIF C C-- GET THE NORMAL AND PERTURB IT IF THERE IS SOME ROUGHNESS. C NX = INX(I) NY = INY(I) NZ = INZ(I) IF( MRGH(MATL) .GT. 0.0 )THEN NX = NX + RANHLF(IDUM) * MRGH(MATL) NY = NY + RANHLF(IDUM) * MRGH(MATL) NZ = NZ + RANHLF(IDUM) * MRGH(MATL) CALL NRMVEC( NX, NY, NZ ) ENDIF C C-- DIFFUSE TERM. C D = DX * NX + DY * NY + DZ * NZ IF( D .LT. 0.0 )D = 0.0 C C-- SPECULAR TERM. C VX = - IPX(I) VY = - IPY(I) VZ = - IPZ(I) CALL NRMVEC( VX, VY, VZ ) C HX = 0.5 * ( DX + VX ) HY = 0.5 * ( DY + VY ) HZ = 0.5 * ( DZ + VZ ) C S = HX * NX + HY * NY + HZ * NZ IF( S .LT. 0.0 )THEN S = 0.0 ELSE S = S ** MGLS(MATL) ENDIF C C-- FIND COLOUR DUE TO THIS LIGHT. C CR = LCLR(L) * LATTR * ( D * MCDR(MATL) + S * MCSR(MATL) ) CG = LCLG(L) * LATTG * ( D * MCDG(MATL) + S * MCSG(MATL) ) CB = LCLB(L) * LATTB * ( D * MCDB(MATL) + S * MCSB(MATL) ) C C-- ADD TO INTERSECTION COLOUR. C ICOLR(I) = ICOLR(I) + CR ICOLG(I) = ICOLG(I) + CG ICOLB(I) = ICOLB(I) + CB 1 CONTINUE C RETURN END C C SUBROUTINE LTFEEL( SRFIDX, PX, PY, PZ, LX, LY, LZ, + LATTR, LATTG, LATTB ) IMPLICIT CHARACTER*1 (A-Z) INTEGER SRFIDX REAL PX, PY, PZ, LX, LY, LZ, LATTR, LATTG, LATTB C********************************************************************** C GENERATE A LIGHT FEELER RAY FROM SURFACE POINT (PX,PY,PZ) TO LIGHT (LX,LY,LZ). C SEE WHAT THIS INTERSECTS AND GENERATE AN APPROPRIATE (ISH) LIGHT ATTENUATION FACTOR. C THIS TRIES TO DEAL WITH TRANSPARENT OBJECTS AS WELL AS OPAQUE ONES. IT DOES SO C INCORRECTLY, AS REFRACTION CANNOT BE TAKEN INTO ACCOUNT. C********************************************************************** *CALL PARAMS *CALL STATS *CALL SCREEN *CALL INSLST *CALL PRMLST *CALL MTRLST C REAL DX, DY, DZ, DS, ACATR, ACATG, ACATB, MINCAT INTEGER PRIM, ASECT, MATL REAL TVAL, POX, POY, POZ REAL OBDX, OBDY, OBDZ,OBDS C C-- FIND DISTANCE BETWEEN SURFACE POINT & LIGHT AND A NORMALIZED VECTOR C-- POINTING FROM THE SURFACE TO THE LIGHT. IF POINT IS VERY CLOSE TO LIGHT, THERE IS C-- NO ATTENUATION. C LATTR = 1.0 LATTG = 1.0 LATTB = 1.0 TOTFEL = TOTFEL + 1 C DX = LX - PX DY = LY - PY DZ = LZ - PZ DS = DX * DX + DY * DY + DZ * DZ IF( DS .LT. MINEPS )THEN LATTR = 1.0 LATTG = 1.0 LATTB = 1.0 RETURN ENDIF C C-- ATTEMPT TO GET AROUND DX,DY,DZ APPEARING IN GENSCT/CELSCT AS C-- INDEFINITES WHEN COMPILED OPT=2 ON CDC FTN 5 C-- NOTE THAT JUST *OUTPUTTING* DX,DY,DZ IN A WRITE STATEMENT *DID* C-- FIX THE PROBLEM. IT LOOKS LIKE THE OPTIMIZER HAD DECIDED THAT C-- DS OR DX,DY,DZ ARE SUPERFLUOUS BY THIS POINT. (WRONGLY, OBVIOUSLY). C-- *EVENTUALLY* ;-) INTRODUCING OBDS FINALLY FIXED THIS ACCEPTABLY. C OBDS = SQRT( DS ) OBDX = DX / OBDS OBDY = DY / OBDS OBDZ = DZ / OBDS C C-- FIND ALL INTERSECTING PRIMITIVES. IF A PRIMITIVE IS MADE OF FULLY C-- OPAQUE MATERIAL, OR THE ACCUMULATED ATTENUATION IS VERY GREAT, C-- ABANDON SEARCH AND RETURN LATT=0. C ACATR = 1.0 ACATG = 1.0 ACATB = 1.0 C POX = PX POY = PY POZ = PZ C 1 CONTINUE IF( ACCEL .EQ. NOACL )THEN CALL GENSCT( POX, POY, POZ, OBDX, OBDY, OBDZ, SHDINS, + ASECT, TVAL, 0 ) ELSE CALL CELSCT( POX, POY, POZ, OBDX, OBDY, OBDZ, SHDINS, + ASECT, TVAL, 0 ) ENDIF C C-- SEE IF THERE IS AN INTERSECTION BETWEEN THE SURFACE POINT BEING SHADED C-- AND THE LIGHT (NOT BEYOND IT). C IF( ASECT .EQ. 0 .OR. TVAL .LT. 0.0 .OR. TVAL .GT. DS )GOTO 2 IF( ASECT .NE. 0 )THEN C C-- CHECK INTERSECTION IS NOT WITH PRIMITIVE WE ARE FINDING SHADOWS FOR. C-- (PRIMITIVES CANNOT SHADE THEMSELVES - DIFFICULT NUMERICAL ISSUES IF THEY COULD!). C IF( ASECT .NE. SRFIDX )THEN C-- PRIMITIVE INDEX WE HAVE INTERSECTED WITH. PRIM = IPRIM(SHDINS) C-- MATERIAL INDEX PRIMITIVE IS MADE FROM. MATL = MATTER(PRIM) IF( MKT(MATL) .LT. MINWGT )THEN LATTR = 0.0 LATTG = 0.0 LATTB = 0.0 RETURN ENDIF C ACATR = ACATR * ( MKT(MATL) * MCDR(MATL) ) ACATG = ACATG * ( MKT(MATL) * MCDG(MATL) ) ACATB = ACATB * ( MKT(MATL) * MCDB(MATL) ) C MINCAT = ACATR IF( ACATG .LT. MINCAT )MINCAT = ACATG IF( ACATB .LT. MINCAT )MINCAT = ACATB IF( MINCAT .LT. MINWGT )THEN LATTR = 0.0 LATTG = 0.0 LATTB = 0.0 RETURN ENDIF ENDIF C C-- FOUND AN INTERSECTION (POSSIBLY WITH SELF). C-- MOVE A TAD FURTHER ALONG THE FEELER AND SEE IF WE INTERSECT ANYTHING ELSE. C-- PRECISION! C POX = IPX(SHDINS) + EPS * OBDX POY = IPY(SHDINS) + EPS * OBDY POZ = IPZ(SHDINS) + EPS * OBDZ ENDIF GOTO 1 2 CONTINUE C LATTR = ACATR LATTG = ACATG LATTB = ACATB RETURN END C C---------------------------------------------------------------------- C INTERSECTION TESTING ROUTINES C---------------------------------------------------------------------- C SUBROUTINE GENSCT( OX, OY, OZ, DX, DY, DZ, NINST, ASECT, TVAL, + URAYID ) IMPLICIT CHARACTER*1 (A-Z) REAL OX, OY, OZ, DX, DY, DZ, TVAL INTEGER NINST, ASECT, URAYID C********************************************************************** C GENERAL INTERSECTION FINDER. GIVEN RAY WITH ORIGIN (OX,OY,OZ) AND UNIT C DIRECTION VECTOR (DX,DY,DZ), THE INDEX OF THE NEXT FREE ELEMENT IN THE C INTERSECTION LIST (NINST), FIND THE INTERSECTION CLOSEST TO (OX,OY,OZ) C FOR ALL PRIMITIVES. IF THERE IS AN INTERSECTION, RETURN C (ASECT) SET TO THE INDEX OF THE PRIMITIVE WE HIT ELSE SET IT TO 0. C IN THE CASE OF AN INTERSECTION, ITS POSITION AND THE SURFACE C NORMAL AT THAT POSITION WILL BE IN INSLST.XXX(NINST). C IF THERE IS AN INTERSECTION, RETURN THE RAY PARAMETER OF THE C INTERSECTION IN (TVAL). C IF (URAYID), USE RAY ID BASED INTERSECTION ACCELERATION. C THIS CAN'T BE DONE FOR SHADOW FEELERS. C********************************************************************** *CALL PARAMS *CALL RAYSTK *CALL INSLST *CALL RAYLST *CALL PRMLST INTEGER I, SECT, PIDX REAL T, TMIN C ASECT = 0 RISECT(CRAY) = 0 DO 1 I=1, NPRM IF( (PRAYID(I) .NE. RRAYID(CRAY)) .OR. (URAYID .EQ. 0) )THEN PIDX = PRMIDX(I) IF( PRMTYP(I) .EQ. SPHPRM )THEN CALL SPHINS( OX, OY, OZ, DX, DY, DZ, PIDX, SECT, T, I ) ELSE IF( PRMTYP(I) .EQ. TRIPRM )THEN CALL TRIINS( OX, OY, OZ, DX, DY, DZ, PIDX, SECT, T, I ) ELSE WRITE(6,100) 100 FORMAT( 1X,'INTERNAL ERROR: GENSCT.' ) WRITE(6,101) 101 FORMAT( 6X,'UNKNOWN PRIMITIVE TYPE.' ) STOP ENDIF IF( SECT .NE. 0 )THEN IF( ASECT .EQ. 0 )THEN TMIN = T CALL CPYINS( NINST, NEWINS ) RISECT(CRAY) = NINST ASECT = I PRAYID(I) = RRAYID(CRAY) ELSE IF( T .LT. TMIN )THEN TMIN = T CALL CPYINS( NINST, NEWINS ) RISECT(CRAY) = NINST ASECT = I RRAYID(I) = RRAYID(CRAY) ENDIF ENDIF ENDIF ENDIF 1 CONTINUE IF( ASECT .NE. 0 )THEN TVAL = TMIN ELSE TVAL = 0.0 ENDIF C RETURN END C C SUBROUTINE CELSCT( OX, OY, OZ, DX, DY, DZ, NINST, ASECT, TVAL, + URAYID ) IMPLICIT CHARACTER*1 (A-Z) REAL OX, OY, OZ, DX, DY, DZ, TVAL INTEGER NINST, ASECT, URAYID C********************************************************************** C FIND AN INTERSECTION USING THE ACCELERATION STRUCTURE. C THE ARGUMENTS ARE THE SAME AS FOR GENSCT(). C********************************************************************** *CALL PARAMS *CALL RAYSTK *CALL INSLST *CALL RAYLST *CALL PRMLST *CALL ACCEL *CALL GRID REAL XCELL, YCELL, ZCELL INTEGER BSECT REAL XI, YI, ZI INTEGER I, SECT, PIDX REAL T, TMIN C C-- CALCULATE THE ORIGIN COORDINATES IN THE CELL COORDINATE SYSTEM. C XCELL = ( ( OX - AXL ) / ADXS ) * ANX YCELL = ( ( OY - AYL ) / ADYS ) * ANY ZCELL = ( ( OZ - AZL ) / ADZS ) * ANZ C C-- CHECK WE ARE INSIDE THE GRID. EXCEPT FOR EYE RAYS, WE SHOULD ALWAYS BE. C-- IF NOT, THEN SEE IF THE RAY INTERSECTS THE BOUNDING BOX OF THE SCENE GEOMETRY. C IF( XCELL .LT. 0 .OR. XCELL .GE. MXCELX .OR. + YCELL .LT. 0 .OR. YCELL .GE. MXCELY .OR. + ZCELL .LT. 0 .OR. ZCELL .GE. MXCELZ )THEN CALL BOXINS( OX, OY, OZ, DX, DY, DZ, BSECT, XI, YI, ZI ) C C-- THE RAY DOESN'T HIT THE BOUNDING BOX. NO INTERSECTION. C IF( BSECT .EQ. 0 )THEN ASECT = 0 RETURN ENDIF C C-- USE THE INTERSECTION POINT TO FIND THE INITIAL GRID CELL TO TRAVERSE FROM. C XCELL = ( ( XI - AXL ) / ADXS ) * ANX YCELL = ( ( YI - AYL ) / ADYS ) * ANY ZCELL = ( ( ZI - AZL ) / ADZS ) * ANZ ENDIF C C-- SAVE RAY DATA FOR USE BY VISCEL() C VASECT = 0 VDX = DX VDY = DY VDZ = DZ VOX = OX VOY = OY VOZ = OZ VNINST = NINST VURAYI = URAYID C C-- TRAVERSE THE CELL GRID. EACH CELL VISITED IS HANDLED BY VISCEL() C-- NOTE CORRECTION TO THE DIRECTION VECTOR FOR GRID SPACE! C CALL GRDTRV( XCELL, YCELL, ZCELL, DX/ADXS, DY/ADYS, DZ/ADZS, + ANX, ANY, ANZ ) C C-- GRID TRAVERSAL ENDS EITHER WHEN WE HAVE FALLEN OFF THE SIDES OF THE GRID, C-- OR VISCEL() HAS FOUND AN INTERSECTION. RETURN (ASECT) IN EITHER CASE. C ASECT = VASECT C C-- NOW CHECK ALL THE PRIMITIVES THAT WERE TOO LARGE TO BE INDEXED IN THE C-- ACCELERATION STRUCTURE. ONE OF THESE MAY HAVE A CLOSER INTERSECTION. C IF( ASECT .NE. 0 )TMIN = VTMIN C DO 1 I=1,NPRM IF( PACCEL(I) .EQ. 0 )THEN IF( PRAYID(I) .NE. RRAYID(CRAY) .OR. URAYID .EQ. 0 )THEN PIDX = PRMIDX(I) IF( PRMTYP(I) .EQ. SPHPRM )THEN CALL SPHINS( OX, OY, OZ, DX, DY, DZ, PIDX, SECT, + T, I ) ELSE IF( PRMTYP(I) .EQ. TRIPRM )THEN CALL TRIINS( OX, OY, OZ, DX, DY, DZ, PIDX, SECT, + T, I ) ELSE WRITE(6,100) 100 FORMAT( 1X,'INTERNAL ERROR: CELSCT.' ) WRITE(6,101) 101 FORMAT( 6X,'UNKNOWN PRIMITIVE TYPE.' ) STOP ENDIF IF( SECT .NE. 0 )THEN IF( ASECT .EQ. 0 )THEN TMIN = T CALL CPYINS( NINST, NEWINS ) RISECT(CRAY) = NINST ASECT = I PRAYID(I) = RRAYID(CRAY) ELSE IF( T .LT. TMIN )THEN TMIN = T CALL CPYINS( NINST, NEWINS ) RISECT(CRAY) = NINST ASECT = I PRAYID(I) = RRAYID(CRAY) ENDIF ENDIF ENDIF ENDIF ENDIF 1 CONTINUE C IF( ASECT .NE. 0 )THEN TVAL = TMIN ELSE TVAL = 0.0 ENDIF C RETURN END C C SUBROUTINE VISCEL( X, Y, Z, QUIT ) IMPLICIT CHARACTER*1 (A-Z) INTEGER X, Y, Z, QUIT C********************************************************************** C PROCESS GRID CELL (X,Y,Z). IF THERE ARE ANY PRIMITIVES IN THIS CELL, C CHECK THEM FOR INTERSECTION WITH THE CURRENT RAY (IN VSCACC.XXX). C IF AN INTERSECTION IS FOUND, LEAVE THE PRIMITIVE INDEX IN VSCACC.ASECT C WITH THE POSITION AND NORMAL STORED IN INSLST.XXX(VSCACC.NINST) IN C THE USUAL WAY. AND TERMINATE GRID TRAVERSAL. C (IF WE WANT TO STOP GRID TRAVERSAL AT THIS CELL, RETURN (QUIT) C NON-ZERO, ELSE 0.) C********************************************************************** *CALL PARAMS *CALL RAYSTK *CALL RAYLST *CALL PRMLST *CALL ACCEL *CALL GRID *CALL INSLST INTEGER CI, LI, LP, PIDX INTEGER SECT, LASECT, CX, CY, CZ REAL XCELL, YCELL, ZCELL, T, TMIN C C-- SEE IF THERE IS A LIST OF PRIMITIVES FOR THIS CELL. C-- IF NOT CONTINUE TRAVERSAL. C CI = X + Y * ANX + Z * ( ANX * ANY ) + 1 LI = CELLS(CI) IF( LI .EQ. 0 )THEN VASECT = 0 QUIT = 0 RETURN ENDIF C C-- THERE IS A LIST OF PRIMITIVES STARTING AT LIST INDEX LI. C-- MOVE OVER THIS LIST LOOKING FOR INTERSECTIONS WITH EACH C-- PRIMITIVE ON THE LIST IN TURN. C LASECT = 0 1 CONTINUE LP = LPRIM(LI) IF( PRAYID(LP) .NE. RRAYID(CRAY) .OR. VURAYI .EQ. 0 )THEN PIDX = PRMIDX(LP) IF( PRMTYP(LP) .EQ. SPHPRM )THEN CALL SPHINS( VOX, VOY, VOZ, VDX, VDY, VDZ, PIDX, + SECT, T, LP ) ELSE IF( PRMTYP(LP) .EQ. TRIPRM )THEN CALL TRIINS( VOX, VOY, VOZ, VDX, VDY, VDZ, PIDX, + SECT, T, LP ) ELSE WRITE(6,100) 100 FORMAT( 1X,'INTERNAL ERROR: VISCEL.' ) WRITE(6,101) 101 FORMAT( 6X,'UNKNOWN PRIMITIVE TYPE.' ) STOP ENDIF IF( SECT .NE. 0 )THEN C C-- CHECK THAT THIS INTERSECTION LIES INSIDE THIS GRID CELL. C XCELL = ( ( IPX(NEWINS) - AXL ) / ADXS ) * ANX YCELL = ( ( IPY(NEWINS) - AYL ) / ADYS ) * ANY ZCELL = ( ( IPZ(NEWINS) - AZL ) / ADZS ) * ANZ CX = INT(XCELL) CY = INT(YCELL) CZ = INT(ZCELL) C IF( X .EQ. CX .AND. Y .EQ. CY .AND. Z .EQ. CZ )THEN C C-- YES - A VALID INTERSECTION. IS IT THE CLOSEST IN THIS CELL? C IF( LASECT .EQ. 0 )THEN TMIN = T CALL CPYINS( VNINST, NEWINS ) RISECT(CRAY) = VNINST LASECT = LP PRAYID(LP) = RRAYID(CRAY) ELSE IF( T .LT. TMIN )THEN TMIN = T CALL CPYINS( VNINST, NEWINS ) RISECT(CRAY) = VNINST LASECT = LP PRAYID(LP) = RRAYID(CRAY) ENDIF ENDIF ENDIF ENDIF ENDIF C LI = LNEXT(LI) IF( LI .EQ. 0 )GOTO 2 GOTO 1 2 CONTINUE C C-- NO INTERSECTION WITH PRIMITIVES IN THIS CELL. CONTINUE TRAVERSAL. C IF( LASECT .EQ. 0 )THEN VASECT = 0 QUIT = 0 ELSE C C-- INTERSECTION FOUND. TERMINATE TRAVERSAL. C VASECT = LASECT VTMIN = TMIN QUIT = 1 ENDIF C RETURN END C C SUBROUTINE SPHINS( OX, OY, OZ, DX, DY, DZ, I, SECT, T, INPRIM ) IMPLICIT CHARACTER*1 (A-Z) REAL OX, OY, OZ, DX, DY, DZ, T INTEGER I, SECT, INPRIM C********************************************************************** C SEE IF THE RAY FROM (OX,OY,OZ) DIRECTION (DX,DY,DZ) INTERSECTS C SPHERE I. THE PRIMITIVE INDEX IS INPRIM. C IF SO, SET SECT NON-ZERO AND SET T TO THE INTERSECTION PARAMETER. C OTHERWISE SET SECT TO ZERO. C IF THERE IS AN INTERSECTION SAVE THE POSITION AND NORMAL IN C INSLST.XXX(NEWINS). C********************************************************************** *CALL PARAMS *CALL SPHLST *CALL INSLST *CALL STATS REAL OCX, OCY, OCZ, RS REAL B, C, D, T0, T1, TMIN C C-- INITIALIZE. C C CALL PMDDUMP NSECTS = NSECTS + 1 SECT = 0 T = 0.0 C C-- SOLVE IMPLICIT QUADRATIC EQUATION C OCX = OX - SOX(I) OCY = OY - SOY(I) OCZ = OZ - SOZ(I) RS = SRAD(I) B = -2 * ( DX * OCX + DY * OCY + DZ * OCZ ) C = OCX * OCX + OCY * OCY + OCZ * OCZ - RS * RS D = B * B - 4 * C C C-- IF DISCRIMINANT -VE, NO INTERSECTION. C IF( D .GE. 0.0 )THEN D = SQRT( D ) C C-- OK - TWO POSSIBLE INTERSECTIONS TO BE FOUND - T0 AND T1. C-- PRECISION! C C-- T1 MUST BE .GT. T0 ... T1 WOULD BE EXITING THE SPHERE IF T0 ALSO VALID. C T1 = 0.5 * ( B + D ) IF( T1 .GT. MINEPS )THEN C C-- T1 IS AHEAD, SO WE HAVE AN INTERSECTION. C-- SO T0 MIGHT BE AHEAD AND CLOSER ... C TMIN = T1 SECT = 1 C C-- THIS WOULD BE LEAVING THE SPHERE. C ITTYP(NEWINS) = LEAVE C T0 = 0.5 * ( B - D ) IF( T0 .GT. MINEPS )THEN C C-- T0 IS AHEAD, WE HAVE A CLOSER INTERSECTION ... C TMIN = T0 C C-- THIS WOULD BE ENTERING THE SPHERE. C ITTYP(NEWINS) = ENTER ENDIF C C-- FILL IN ALL THE REQUIRED INFO FOR THE NEW INTERSECTION C IPRIM(NEWINS) = INPRIM IRRAY(NEWINS) = 0 ITRAY(NEWINS) = 0 ICOLR(NEWINS) = 0.0 ICOLG(NEWINS) = 0.0 ICOLB(NEWINS) = 0.0 IPX(NEWINS) = OX + TMIN * DX IPY(NEWINS) = OY + TMIN * DY IPZ(NEWINS) = OZ + TMIN * DZ INX(NEWINS) = ( IPX(NEWINS) - SOX(I) ) / RS INY(NEWINS) = ( IPY(NEWINS) - SOY(I) ) / RS INZ(NEWINS) = ( IPZ(NEWINS) - SOZ(I) ) / RS C C--- HAND BACK THE INTERSECTION PARAMETER C T = TMIN ENDIF ENDIF C RETURN END C C SUBROUTINE TRIINS( OX, OY, OZ, DX, DY, DZ, I, SECT, T, INPRIM ) IMPLICIT CHARACTER*1 (A-Z) REAL OX, OY, OZ, DX, DY, DZ, T INTEGER I, SECT, INPRIM C********************************************************************** C SEE IF THE RAY FROM (OX,OY,OZ) DIRECTION (DX,DY,DZ) INTERSECTS C TRIANGLE I. THE PRIMITIVE INDEX IS INPRIM. C IF SO, SET SECT NON-ZERO AND SET T TO THE INTERSECTION PARAMETER. C OTHERWISE SET SECT TO ZERO. C IF THERE IS AN INTERSECTION SAVE THE POSITION AND NORMAL C IN INSLST.XXX(NEWINS). C********************************************************************** *CALL PARAMS *CALL TRILST *CALL INSLST *CALL STATS REAL NDOTD, NDOTO, PX, PY, PZ, U0, U1, U2, V0, V1, V2 REAL ALPHA, BETA, GAMMA, NPX, NPY, NPZ C C-- INITIALIZE. BUMP INTERSECTION TEST COUNT. C NSECTS = NSECTS + 1 C C-- CHECK RAY IS NOT PARALLEL TO PLANE EMBEDDING TRIANGLE. C-- PRECISION! C NDOTD = TNX(I) * DX + TNY(I) * DY + TNZ(I) * DZ IF( ABS( NDOTD ) .GT. EPS )THEN C C-- FIND INTERSECTION PARAMETER WITH PLANE. C NDOTO = TNX(I) * OX + TNY(I) * OY + TNZ(I) * OZ T = -( NDOTO + TND(I) ) / NDOTD C C-- FIND INTERSECTION POINT WITH EMBEDDING PLANE. C PX = OX + T * DX PY = OY + T * DY PZ = OZ + T * DZ C C-- FIND CONSTANTS FOR LOCATING INTERSECTION POINT IN TERMS OF PARAMETRIC C-- COORDINATES IN A COORDINATE SYSTEM USING TWO EDGES OF THE TRIANGLE AS C-- ITS BASIS VECTORS. USE THE "BIGGEST NORMAL COMPONENT" PRECOMPUTED AS DAX C-- TO SELECT THE PLANE NORMAL TO THE LEAST DEGENERATE POSSIBLE PARAMETRIZATION. C IF( DAX(I) .EQ. DAXX )THEN U0 = PY - TY1(I) U1 = TY2(I) - TY1(I) U2 = TY3(I) - TY1(I) V0 = PZ - TZ1(I) V1 = TZ2(I) - TZ1(I) V2 = TZ3(I) - TZ1(I) ELSE IF( DAX(I) .EQ. DAXY )THEN U0 = PX - TX1(I) U1 = TX2(I) - TX1(I) U2 = TX3(I) - TX1(I) V0 = PZ - TZ1(I) V1 = TZ2(I) - TZ1(I) V2 = TZ3(I) - TZ1(I) ELSE IF( DAX(I) .EQ. DAXZ )THEN U0 = PX - TX1(I) U1 = TX2(I) - TX1(I) U2 = TX3(I) - TX1(I) V0 = PY - TY1(I) V1 = TY2(I) - TY1(I) V2 = TY3(I) - TY1(I) ELSE WRITE(6,100) 100 FORMAT( 1X,'INTERNAL ERROR: TRIINS.' ) WRITE(6,101) 101 FORMAT( 6X,'UNKNOWN PRIMITIVE TYPE.' ) STOP ENDIF C C-- FIND THE (ALPHA,BETA) COORDINATES OF THE INTERSECTION POINT. C-- PRECISION! C IF( ABS( U1 ) .LT. MINEPS )THEN BETA = U0 / U2 IF( BETA .GE. 0.0 .AND. BETA .LE. 1.0 )THEN ALPHA = ( V0 - BETA * V2 ) / V1 ELSE T = 0.0 SECT = 0 RETURN ENDIF ELSE BETA = ( V0 * U1 - U0 * V1 ) / ( V2 * U1 - U2 * V1 ) IF( BETA .GE. 0.0 .AND. BETA .LE. 1.0 )THEN ALPHA = ( U0 - BETA * U2 ) / U1 ELSE T = 0.0 SECT = 0 RETURN ENDIF ENDIF C C-- FIND THE GAMMA COORDINATE. C GAMMA = ALPHA + BETA IF( ALPHA .LT. 0.0 .OR. GAMMA .GT. 1.0 )THEN T = 0.0 SECT = 0 RETURN ENDIF C C-- WE ARE INSIDE THE TRIANGLE. C-- USE (ALPHA,BETA,GAMMA) TO INTERPOLATE THE VERTEX PROPERTIES (CURRENTLY THE NORMAL). C NPX = GAMMA * NVX1(I) + ALPHA * NVX2(I) + BETA * NVX3(I) NPY = GAMMA * NVY1(I) + ALPHA * NVY2(I) + BETA * NVY3(I) NPZ = GAMMA * NVZ1(I) + ALPHA * NVZ2(I) + BETA * NVZ3(I) CALL NRMVEC( NPX, NPY, NPZ ) C C-- SET THE INTERSECTION DATA. C IPRIM(NEWINS) = INPRIM IRRAY(NEWINS) = 0 ITRAY(NEWINS) = 0 ICOLR(NEWINS) = 0.0 ICOLG(NEWINS) = 0.0 ICOLB(NEWINS) = 0.0 IPX(NEWINS) = PX IPY(NEWINS) = PY IPZ(NEWINS) = PZ INX(NEWINS) = NPX INY(NEWINS) = NPY INZ(NEWINS) = NPZ C C-- GUESS WHETHER WE ARE ENTERING OR LEAVING. C NDOTD = NPX * DX + NPY * DY + NPZ * DZ IF( NDOTD .LT. 0.0 )THEN ITTYP(NEWINS) = ENTER ELSE ITTYP(NEWINS) = LEAVE ENDIF SECT = 1 C C-- RAY IS PARALLEL TO EMBEDDING PLANE. NO INTERSECTION POSSIBLE. C ELSE T = 0.0 SECT = 0 ENDIF C RETURN END C C SUBROUTINE CPYINS( I, J ) IMPLICIT CHARACTER*1 (A-Z) INTEGER I, J C********************************************************************** C MOVE A NEWLY FOUND INTERSECTION FROM INSLST.XXX(J) TO INSLST.XXX(I) C********************************************************************** *CALL INSLST IPRIM(I) = IPRIM(J) IRRAY(I) = IRRAY(J) ITRAY(I) = ITRAY(J) ICOLR(I) = ICOLR(J) ICOLG(I) = ICOLG(J) ICOLB(I) = ICOLB(J) IPX(I) = IPX(J) IPY(I) = IPY(J) IPZ(I) = IPZ(J) INX(I) = INX(J) INY(I) = INY(J) INZ(I) = INZ(J) ITTYP(I) = ITTYP(J) C RETURN END C C---------------------------------------------------------------------- C DATABASE ROUTINES C---------------------------------------------------------------------- C SUBROUTINE INITDB IMPLICIT CHARACTER*1 (A-Z) C********************************************************************** C INITIALIZE THE SCENE DATABASE C********************************************************************** *CALL PRMLST *CALL SPHLST *CALL MTRLST *CALL LGTLST *CALL RAYLST *CALL TRILST NPRM = 0 NSPH = 0 NTRI = 0 NMTR = 0 NLGT = 0 RAYNUM = 0 C RETURN END C C SUBROUTINE INITBK IMPLICIT CHARACTER*1 (A-Z) C********************************************************************** C INITIALISE THINGS FOR GRADED BACKGROUNDS C********************************************************************** *CALL PARAMS *CALL SCREEN REAL FN C IF( BKT .EQ. BKVERT )THEN FN = SNY - 1 ELSE IF( BKT .EQ. BKHORZ )THEN FN = SNX - 1 ELSE FN = 1.0 ENDIF BKRD = ( BKRF - BKRI ) / FN BKGD = ( BKGF - BKGI ) / FN BKBD = ( BKBF - BKBI ) / FN C RETURN END C C SUBROUTINE INITLG IMPLICIT CHARACTER*1 (A-Z) C********************************************************************** C NORMALIZE AND INVERT THE LIGHT POINTING VECTORS FOR DIRECTIONAL LIGHTS. C********************************************************************** *CALL PARAMS *CALL LGTLST INTEGER L C DO 1 L=1,NLGT IF( LDIR(L) .NE. LGTCON )THEN CALL NRMVEC( LDX(L), LDY(L), LDZ(L) ) LDX(L) = -LDX(L) LDY(L) = -LDY(L) LDZ(L) = -LDZ(L) ENDIF 1 CONTINUE C RETURN END C C SUBROUTINE INITRI( I ) IMPLICIT CHARACTER*1 (A-Z) INTEGER I C********************************************************************** C PRECOMPUTE PLANE EQUATION AND DOMINANT AXIS INFO FOR TRIANGLE (I) C********************************************************************** *CALL PARAMS *CALL TRILST REAL DXA, DYA, DZA, DXB, DYB, DZB REAL FNX, FNY, FNZ, LEN, NVX, NVY, NVZ C DXA = TX2(I) - TX1(I) DYA = TY2(I) - TY1(I) DZA = TZ2(I) - TZ1(I) DXB = TX3(I) - TX1(I) DYB = TY3(I) - TY1(I) DZB = TZ3(I) - TZ1(I) TNX(I) = DYA * DZB - DZA * DYB TNY(I) = DZA * DXB - DXA * DZB TNZ(I) = DXA * DYB - DYA * DXB TND(I) = - ( TX1(I) * TNX(I) + TY1(I) * TNY(I) + + TZ1(I) * TNZ(I) ) IF( TNX(I) .GT. TNY(I) .AND. TNX(I) .GT. TNZ(I) )THEN DAX(I) = DAXX ELSE IF( TNY(I) .GT. TNX(I) .AND. TNY(I) .GT. TNZ(I) )THEN DAX(I) = DAXY ELSE DAX(I) = DAXZ ENDIF C C-- IF THE VERTEX NORMAL LENGTH FOR A VERTEX IS TINY, SUBSTITUTE A UNIT FACET NORMAL. C-- OTHERWISE, ENSURE THE VERTEX NORMALS ARE UNIT VECTORS. C FNX = TNX(I) FNY = TNY(I) FNZ = TNZ(I) CALL NRMVEC( FNX, FNY, FNZ ) C NVX = NVX1(I) NVY = NVY1(I) NVZ = NVZ1(I) LEN = NVX * NVX + NVY * NVY + NVZ * NVZ C-- PRECISION! IF( LEN .LT. (MINEPS*MINEPS) )THEN NVX1(I) = FNX NVY1(I) = FNY NVZ1(I) = FNZ ELSE LEN = SQRT( LEN ) NVX1(I) = NVX1(I) / LEN NVY1(I) = NVY1(I) / LEN NVZ1(I) = NVZ1(I) / LEN ENDIF C NVX = NVX2(I) NVY = NVY2(I) NVZ = NVZ2(I) LEN = NVX * NVX + NVY * NVY + NVZ * NVZ C-- PRECISION! IF( LEN .LT. (MINEPS*MINEPS) )THEN NVX2(I) = FNX NVY2(I) = FNY NVZ2(I) = FNZ ELSE LEN = SQRT( LEN ) NVX2(I) = NVX2(I) / LEN NVY2(I) = NVY2(I) / LEN NVZ2(I) = NVZ2(I) / LEN ENDIF C NVX = NVX3(I) NVY = NVY3(I) NVZ = NVZ3(I) LEN = NVX * NVX + NVY * NVY + NVZ * NVZ C-- PRECISION! IF( LEN .LT. (MINEPS*MINEPS) )THEN NVX3(I) = FNX NVY3(I) = FNY NVZ3(I) = FNZ ELSE LEN = SQRT( LEN ) NVX3(I) = NVX3(I) / LEN NVY3(I) = NVY3(I) / LEN NVZ3(I) = NVZ3(I) / LEN ENDIF C RETURN END C C SUBROUTINE READDB( RDOK ) IMPLICIT CHARACTER*1 (A-Z) INTEGER RDOK C********************************************************************** C READ THE SCENE DATABASE. C THE FORMAT OF THE DATABASE IS PRETTY HORRIBLE. C RETURN 1 IN RDOK IF OK, 0 IF SOMETHING BAD HAPPENED. C C WE USE LIST DIRECTED IO. ANY INPUT THAT SUPPLIES THE REQUIRED DATA C AS SET OUT BELOW AND IS COMPATIBLE WITH THE RULES FOR LIST DIRECTED C READS IS ACCEPTABLE. THE INPUT DATA IS DESCRIBED IN TERMS OF "GROUPS" C BELOW. SOME SETS OF GROUPS MUST BE REPEATED IN THE INPUT TO SUPPLY C THE SPECIFIED NUMBER OF ITEMS. HOPEFULLY WHAT THIS MEANS WILL BE C CLEAR WHEN READING THE DESCRIPTION THAT FOLLOWS. C C GROUP 1: TRACING CONTROL PARAMETERS C ACCEL SHADOW C (INT) (INT) C ACCEL = 1 FOR ACCELERATED INTERSECTION TESTING (CELSCT) C 0 FOR THE SIMPLER "TRY EVERYTHING" (GENSCT) C SHADOW = 0 FOR NO SHADOWS C N (INTEGER) FOR NUMBER OF FEELER RAYS TO USE FOR SHADOWS. C N .GT. 1 CAN GIVE SOFT SHADOWS. C C GROUP 2: OUTPUT IMAGE PARAMETERS C NX NY OS C (INT) (INT) (INT) C NX = NUMBER OF OUTPUT PIXELS C NY = NUMBER OF OUTPUT LINES C OS = OVERSAMPLING RATE FOR JITTERED STOCHASTIC ANTI-ALIASING C C GROUP 3: CAMERA PARAMETERS C APER FPD FOCAL FSTOP C (REAL) (REAL) (REAL) (REAL) C APER = APERTURE IN FILM PLANE (E.G. IN MM) C FPD = FOCAL PLANE DISTANCE (E.G. IN MM) C DISTANCE FROM LENS CENTER ON WHICH YOU ARE FOCUSSING. C FOCAL = FOCAL LENGTH OF LENS (E.G. IN MM) C FSTOP = LENS FSTOP (USED FOR DEPTH OF FIELD - OS .GT. 1) C C GROUP 4: TRACING TERMINATION PARAMETERS C MAXDEP MINWGT C (INT) (REAL) C MAXDEP = MAXIMUM TREE DEPTH OVER WHICH TRACING A RAY WILL STOP. C MINWGT = MINIMUM WEIGHT UNDER WHICH TRACING A RAY WILL STOP. C C GROUP 5: BACKGROUND COLOUR PARAMETERS C BKTYPE IBACKR IBACKG IBACKB FBACKR FBACKG FBACKB C (INT) (REAL) (REAL) (REAL) (REAL) (REAL) REAL) C BKTYPE = 0:CONSTANT, 1:VERTICAL GRADE, 2:HORIZ GRADE C C GROUP 6: MATERIALS COUNT C NUMMTR C (INT) C NUMMTR = NUMBER OF MATERIALS BEING SPECIFIED C GROUPS 7,8,9,10 MUST FOLLOW. THERE MUST BE NUMMTR REPETITIONS C OF THIS SET OF GROUPS C C GROUP 7: AMBIENT MATERIAL COLOUR C CAR CAG CAB C (REAL) (REAL) (REAL) C CAR = AMBIENT LIGHT REFLECTED (RED COMPONENT) C CAG = AMBIENT LIGHT REFLECTED (GREEN COMPONENT) C CAB = AMBIENT LIGHT REFLECTED (BLUE COMPONENT) C ALL USUALLY IN THE RANGE 0 TO 1 C C GROUP 8: DIFFUSE MATERIAL COLOUR C CDR CDG CDB C (REAL) (REAL) (REAL) C CDR = DIFFUSE LIGHT REFLECTED (RED COMPONENT) C CDG = DIFFUSE LIGHT REFLECTED (GREEN COMPONENT) C CDB = DIFFUSE LIGHT REFLECTED (BLUE COMPONENT) C ALL USUALLY IN THE RANGE 0 TO 1 C C GROUP 9: SPECULAR MATERIAL COLOUR C CSR CSG CSB C (REAL) (REAL) (REAL) C CSR = SPECULAR LIGHT REFLECTED (RED COMPONENT) C CSG = SPECULAR LIGHT REFLECTED (GREEN COMPONENT) C CSB = SPECULAR LIGHT REFLECTED (BLUE COMPONENT) C ALL USUALLY IN THE RANGE 0 TO 1 C C GROUP 10: VARIOUS MATERIAL PARAMETERS C KR KT ETA GLOSS ROUGH C (REAL) (REAL) (REAL) (REAL) (REAL) C KR = COEFFICIENT OF REFLECTION C KT = COEFFICIENT OF TRANSMISSION C KR + KT .LE. 1.0 USUALLY C ETA = REFRACTIVE INDEX (FOR TRANSMISSION) C ETA = 1.0 FOR "AIR", "GLASS" AROUND 1.45, ETC. C GLOSS = SHARPNESS OF SPECULAR HIGHLIGHT C GLOSS = 1.0 TO 250.0 (RAISED COSINE PROFILE) C ROUGH = ROUGHNESS OF SURFACE (RANDOM NORMAL PERTURBATION) C ROUGH = 0.0 TO 1.0 C C GROUP 11: AMBIENT LIGHT C LAR LAG LAB C (REAL) (REAL) (REAL) C LAR = AMBIENT LIGHT (RED COMPONENT) C LAG = AMBIENT LIGHT (GREEN COMPONENT) C LAB = AMBIENT LIGHT (BLUE COMPONENT) C ALL USUALLY IN THE RANGE 0 TO 1 C C GROUP 12: LIGHTS COUNT C NUMLGT C (INT) C NUMLGT = NUMBER OF LIGHTS TO USE. C GROUPS 13,14 MUST FOLLOW. THERE MUST BE NUMLGT REPETITIONS C OF THIS PAIR OF GROUPS. C C GROUP 13: LIGHT PROPERTIES 1 C OX OY OZ CLR CLG CLB RADIUS C (REAL) (REAL) (REAL) (REAL) (REAL) (REAL) (REAL) C OX = X COORDINATE OF LIGHT (CAMERA SPACE) (E.G. IN MM) C OY = Y COORDINATE OF LIGHT (CAMERA SPACE) (E.G. IN MM) C OZ = Z COORDINATE OF LIGHT (CAMERA SPACE) (E.G. IN MM) C CLR = RED COMPONENT OF LIGHT COLOUR (0 TO 1 USUALLY) C CLG = GREEN COMPONENT OF LIGHT COLOUR (0 TO 1 USUALLY) C CLB = BLUE COMPONENT OF LIGHT COLOUR (0 TO 1 USUALLY) C RADIUS = SIZE OF LIGHT SOURCE (RADIUS OF SPHERE, E.G. IN MM) C C GROUP 14: LIGHT PROPERTIES 2 C DTYPE DX DY DZ ANG1 ANG2 C (INT) (REAL) (REAL) (REAL) (REAL) (REAL) C DTYPE = DIRECTIONAL TYPE OF LIGHT. C 0 = NOT DIRECTIONAL, 1 = RAISED COSINE, 2 = "BARN DOOR" C DX = X COMPONENT OF DIRECTION VECTOR (NOT NORMALISED) C DY = Y COMPONENT OF DIRECTION VECTOR (NOT NORMALISED) C DZ = Z COMPONENT OF DIRECTION VECTOR (NOT NORMALISED) C ANG1 = POWER FOR RAISED COSINE (DTYPE = 1) C = FULL POWER HALF ANGLE (DEGREES) (DTYPE = 2) C ANG2 = ZERO POWER HALF ANGLE (DEGREES) (DTYPE = 2) C C GROUP 15: SPHERE OBJECT COUNT C NUMSPH C (INT) C NUMSPH = NUMBER OF SPHERES IN SCENE. C GROUP 16 MUST FOLLOW. THERE MUST BE NUMSPH REPETITIONS C OF THIS GROUP. C C GROUP 16: SPHERE SPECIFICATION C MTRIDX OX OY OZ RAD C (INT) (REAL) (REAL) (REAL) (REAL) C MTRIDX = INDEX NUMBER OF MATERIAL FROM WHICH SPHERE IS MADE. C MATERIAL INDEX NUMBERS START AT 1 FOR THE FIRST OF C THE NUMMTR MATERIALS DEFINED (SEE ABOVE). C OX = X CENTER POSITION IN CAMERA SPACE (E.G. IN MM) C OY = Y CENTER POSITION IN CAMERA SPACE (E.G. IN MM) C OZ = Z CENTER POSITION IN CAMERA SPACE (E.G. IN MM) C RAD = RADIUS OF SPHERE (E.G. IN MM) C C GROUP 17: TRIANGLE OBJECT COUNT C NUMTRI C (INT) C NUMTRI = NUMBER OF TRIANGLES IN THE SCENE C GROUPS 18,19,20,21 MUST FOLLOW. THERE MUST BE NUMTRI C REPETITIONS OF THIS SET OF GROUPS TO DEFINE THE TRIANGLES. C C GROUP 18: TRIANGLE MATERIAL C MTRIDX C (INT) C MTRIDX = INDEX OF MATERIAL FROM WHICH THE TRIANGLE IS MADE. C MATERIAL INDEX NUMBERS START AT 1 FOR THE FIRST OF C THE NUMMTR MATERIALS DEFINED (SEE ABOVE). C C GROUP 19: FIRST TRIANGLE VERTEX SPECIFICATION C X1 Y1 Z1 NX1 NY1 NZ1 C (REAL) (REAL) (REAL) (REAL) (REAL) (REAL) C X1 = X COORDINATE IN CAMERA SPACE (E.G. IN MM) C Y1 = Y COORDINATE IN CAMERA SPACE (E.G. IN MM) C Z1 = Z COORDINATE IN CAMERA SPACE (E.G. IN MM) C NX1 = X COMPONENT OF NORMAL VECTOR AT THE VERTEX. C NY1 = Y COMPONENT OF NORMAL VECTOR AT THE VERTEX. C NZ1 = Z COMPONENT OF NORMAL VECTOR AT THE VERTEX. C IF NX1=NY1=NZ1 = 0 THEN A FACET NORMAL WILL BE CALCULATED C AND USED AT ALL VERTICES OF THE TRIANGLE. C C GROUP 20: SECOND TRIANGLE VERTEX SPECIFICATION C X2 Y2 Z2 NX2 NY2 NZ2 C (REAL) (REAL) (REAL) (REAL) (REAL) (REAL) C SEE GROUP 19. C C GROUP 21: THIRD TRIANGLE VERTEX SPECIFICATION C X3 Y3 Z3 NX3 NY3 NZ3 C (REAL) (REAL) (REAL) (REAL) (REAL) (REAL) C SEE GROUP 19. C********************************************************************** *CALL PARAMS *CALL SCREEN *CALL MTRLST *CALL LGTLST *CALL PRMLST *CALL SPHLST *CALL TRILST INTEGER I, CLINE CHARACTER*80 CDIAG C C-- SCREEN/CAMERA SPECIFICATION. C CLINE = 1 CDIAG = 'ACCELERATION SCHEME OR SHADOW FEELERS' READ(5,*,END=998,ERR=999)ACCEL, SHADOW CLINE = CLINE + 1 CDIAG = 'RESOLUTION OR OVERSAMPLING' READ(5,*,END=998,ERR=999)SNX, SNY, OS CLINE = CLINE + 1 CDIAG = 'CAMERA PARAMETERS' READ(5,*,END=998,ERR=999)APER, FPD, FOCAL, FSTOP CLINE = CLINE + 1 CDIAG = 'TRACE TERMINATION DEPTH OR WEIGHT' READ(5,*,END=998,ERR=999)SMAXDP, MINWGT C C-- CHECK HORIZONTAL RESOLUTION IS IN LIMITS C IF( SNX .GT. MAXPXL )THEN WRITE(6,100)SNX, MAXPXL 100 FORMAT(1X,'HORIZONTAL RESOLUTION ',I5,' EXCEEDS MAXIMUM ',I5) GOTO 999 ENDIF C C-- FORCE OVERSAMPLE RATE TO BE ODD. C IF( OS .GT. 0 )THEN IF( ((OS/2)*2) .EQ. OS )OS = OS + 1 ENDIF CLINE = CLINE + 1 CDIAG = 'BACKGROUND SPECIFICATION' READ(5,*,END=998,ERR=999)BKT, BKRI, BKGI, BKBI, BKRF, BKGF, BKBF C C-- CHECK GRADE TYPE C IF( BKT .LT. BKCNST .OR. BKT .GT. BKHORZ )THEN WRITE(6,101) 101 FORMAT(1X,'INVALID BACKGROUND GRADE TYPE' ) GOTO 999 ENDIF C C-- INITIALIZE BACKGROUND DATA C CALL INITBK C C-- MATERIALS SPECIFICATION. C CLINE = CLINE + 1 CDIAG = 'MATERIAL COUNT' READ(5,*,END=998,ERR=999)NMTR C C-- CHECK COUNT IS IN LIMITS C IF( NMTR .GT. MAXMTR )THEN WRITE(6,102)NMTR, MAXMTR 102 FORMAT(1X,'TOO MANY MATERIALS ',I5,' LIMIT ',I5) GOTO 999 ENDIF C C-- READ ALL THOSE MATERIALS C DO 1 I=1,NMTR CLINE = CLINE + 1 CDIAG = 'MATERIAL AMBIENT' READ(5,*,END=998,ERR=999)MCAR(I), MCAG(I), MCAB(I) CLINE = CLINE + 1 CDIAG = 'MATERIAL DIFFUSE' READ(5,*,END=998,ERR=999)MCDR(I), MCDG(I), MCDB(I) CLINE = CLINE + 1 CDIAG = 'MATERIAL SPECULAR' READ(5,*,END=998,ERR=999)MCSR(I), MCSG(I), MCSB(I) CLINE = CLINE + 1 CDIAG = 'MATERIAL OTHER' READ(5,*,END=998,ERR=999)MKR(I), MKT(I), META(I), + MGLS(I), MRGH(I) 1 CONTINUE C C-- READ AMBIENT LIGHT C CLINE = CLINE + 1 CDIAG = 'LIGHT AMBIENT' READ(5,*,END=998,ERR=999)LCAR, LCAG, LCAB C C-- LIGHT COUNT C CLINE = CLINE + 1 CDIAG = 'LIGHT COUNT' READ(5,*,END=998,ERR=999)NLGT C C-- CHECK COUNT IS IN LIMITS C IF( NLGT .GT. MAXLGT )THEN WRITE(6,109)NLGT, MAXLGT 109 FORMAT(1X,'TOO MANY LIGHTS ',I5,' LIMIT ',I5) GOTO 999 ENDIF C C-- READ ALL THOSE LIGHTS C DO 2 I=1,NLGT CLINE = CLINE + 1 CDIAG = 'LIGHT POS, COL, RAD' READ(5,*,END=998,ERR=999)LOX(I), LOY(I), LOZ(I), + LCLR(I), LCLG(I), LCLB(I), LRAD(I) CLINE = CLINE + 1 CDIAG = 'LIGHT DIR, VEC, ANG' READ(5,*,END=998,ERR=999)LDIR(I), LDX(I), LDY(I), LDZ(I), + LTGT(I), LTG2(I) IF( LTG2(I) .LE. LTGT(I) )LTG2(I) = LTGT(I) + 0.0001 LTGT(I) = LTGT(I) * DEGRAD LTG2(I) = LTG2(I) * DEGRAD 2 CONTINUE C C-- INITIALIZE THE LIGHTS DATABASE C CALL INITLG C C-- NO GEOMETRIC PRIMITIVES YET C NPRM = 0 C C-- SPHERE COUNT C CLINE = CLINE + 1 CDIAG = 'SPHERE COUNT' READ(5,*,END=998,ERR=999)NSPH C C-- CHECK COUNT IS IN LIMITS C IF( NSPH .GT. MAXSPH )THEN WRITE(6,103)NSPH, MAXSPH 103 FORMAT(1X,'TOO MANY SPHERES ',I5,' LIMIT ',I5) GOTO 999 ENDIF C C-- READ ALL THOSE SPHERES C DO 3 I=1,NSPH NPRM = NPRM + 1 IF( NPRM .GT. MAXPRM )THEN WRITE(6,104)NPRM, MAXPRM 104 FORMAT(1X,'TOO MANY PRIMITIVES ',I5,' LIMIT ',I5) GOTO 999 ENDIF CLINE = CLINE + 1 CDIAG = 'SPHERE' READ(5,*,END=998,ERR=999)MATTER(NPRM), + SOX(I), SOY(I), SOZ(I), SRAD(I) IF( MATTER(NPRM) .GT. NMTR .OR. MATTER(NPRM) .LT. 1 )THEN WRITE(6,105) MATTER(NPRM), CLINE 105 FORMAT(1X,'MTR IDX ',I5,' INVALID. LINE ',I5) GOTO 999 ENDIF C-- NOT IN ACCELERATION STRUCTURE. PACCEL(NPRM) = 0 C-- NO RAY HAS INTERSECTED IT YET. PRAYID(NPRM) = 0 C-- PRIMITIVE TYPE: SPHERE. PRMTYP(NPRM) = SPHPRM C-- INDEX OF DEFINITION IN SPHERE TABLE. PRMIDX(NPRM) = I 3 CONTINUE C C-- TRIANGLE COUNT C CLINE = CLINE + 1 CDIAG = 'TRIANGLE COUNT' READ(5,*,END=998,ERR=999)NTRI C C-- CHECK COUNT IS IN LIMITS C IF( NTRI .GT. MAXTRI )THEN WRITE(6,106)NTRI, MAXTRI 106 FORMAT(1X,'TOO MANY TRIANGLES ',I5,' LIMIT ',I5) GOTO 999 ENDIF C C-- READ ALL THOSE TRIANGLES C DO 4 I=1,NTRI NPRM = NPRM + 1 IF( NPRM .GT. MAXPRM )THEN WRITE(6,104)NPRM, MAXPRM GOTO 999 ENDIF CLINE = CLINE + 1 CDIAG = 'TRIANGLE MTR' READ(5,*,END=998,ERR=999)MATTER(NPRM) IF( MATTER(NPRM) .GT. NMTR .OR. MATTER(NPRM) .LT. 1 )THEN WRITE(6,105) MATTER(NPRM), CLINE GOTO 999 ENDIF CLINE = CLINE + 1 CDIAG = 'TRIANGLE V1' READ(5,*,END=998,ERR=999)TX1(I), TY1(I), TZ1(I), + NVX1(I), NVY1(I), NVZ1(I) CLINE = CLINE + 1 CDIAG = 'TRIANGLE V2' READ(5,*,END=998,ERR=999)TX2(I), TY2(I), TZ2(I), + NVX2(I), NVY2(I), NVZ2(I) CLINE = CLINE + 1 CDIAG = 'TRIANGLE V3' READ(5,*,END=998,ERR=999)TX3(I), TY3(I), TZ3(I), + NVX3(I), NVY3(I), NVZ3(I) C C-- INITIALIZE THE TRIANGLE DATABASE (FOR THIS TRIANGLE) C CALL INITRI( I ) C-- NOT IN ACCELERATION STRUCTURE. PACCEL(NPRM) = 0 C-- NO RAY HAS INTERSECTED IT YET. PRAYID(NPRM) = 0 C-- PRIMITIVE TYPE: SPHERE. PRMTYP(NPRM) = TRIPRM C-- INDEX OF DEFINITION IN SPHERE TABLE. PRMIDX(NPRM) = I 4 CONTINUE C C-- FINISHED SUCCESSFULLY. C RDOK = 1 RETURN C C-- ERROR ON READ C 999 CONTINUE WRITE(6,107)CDIAG, CLINE 107 FORMAT(1X,'READ ERROR DOING: ',A,' LINE ',I5) RDOK = 0 RETURN C C-- END OF FILE ON READ C 998 CONTINUE WRITE(6,108)CLINE 108 FORMAT(1X,'END OF FILE ON READ AT LINE ',I5) RDOK = 0 RETURN END C C---------------------------------------------------------------------- C GRID SPACE SUBDIVISION ACCELERATION ROUTINES C---------------------------------------------------------------------- C SUBROUTINE STACEL IMPLICIT CHARACTER*1 (A-Z) C********************************************************************** C SET UP THE ACCELERATION STRUCTURE. THIS IS A SIMPLE RECTANGULAR SPACE C SUBDIVISION SCHEME. (BUT ARE THESE THINGS EVER REALLY SIMPLE?) C********************************************************************** *CALL PARAMS *CALL STATS *CALL PRMLST *CALL ACCEL INTEGER I, J, K, P, PIDX REAL PXL, PXH, PYL, PYH, PZL, PZH REAL DXC, DYC, DZC REAL CXL, CXH, CYL, CYH, CZL, CZH INTEGER IXL, IXH, IYL, IYH, IZL, IZH REAL MAXDIR, MINDIR, XC, YC, ZC C ANX = MXCELX ANY = MXCELY ANZ = MXCELZ C C-- SET ALL "POINTERS" TO PRIMITIVE LISTS IN THE CELLS TO "NULL". C-- PRIMITIVE LISTS ARE EMPTY C DO 1 K=0,ANZ-1 DO 2 J=0,ANY-1 DO 3 I=0,ANX-1 CELLS(I+J*ANX+K*(ANX*ANY)+1) = 0 3 CONTINUE 2 CONTINUE 1 CONTINUE NLIST = 0 C C-- SURVEY THE SCENE TO GET ITS BOUNDS. C DO 4 P=1,NPRM PIDX = PRMIDX(P) IF( PRMTYP(P) .EQ. SPHPRM )THEN CALL SPHBND( PIDX, PXL, PXH, PYL, PYH, PZL, PZH ) ELSE IF( PRMTYP(P) .EQ. TRIPRM )THEN CALL TRIBND( PIDX, PXL, PXH, PYL, PYH, PZL, PZH ) ELSE WRITE(6,100) 100 FORMAT( 1X,'INTERNAL ERROR: STACEL.' ) WRITE(6,101) 101 FORMAT( 6X,'UNKNOWN PRIMITIVE TYPE.' ) STOP ENDIF IF( P .EQ. 1 )THEN AXL = PXL AXH = PXH AYL = PYL AYH = PYH AZL = PZL AZH = PZH ELSE IF( PXL .LT. AXL )AXL = PXL IF( PXH .GT. AXH )AXH = PXH IF( PYL .LT. AYL )AYL = PYL IF( PYH .GT. AYH )AYH = PYH IF( PZL .LT. AZL )AZL = PZL IF( PZH .GT. AZH )AZH = PZH ENDIF 4 CONTINUE ADXS = AXH - AXL ADYS = AYH - AYL ADZS = AZH - AZL C C-- OUTPUT USEFUL INFORMATION C WRITE(6,102) 102 FORMAT(1X,'ACCELERATION STRUCTURE INFORMATION:') WRITE(6,103)ANX, ANY, ANZ 103 FORMAT(1X,'GRID DIMENSIONS = ',I4,' X ',I4,' X ',I4) WRITE(6,104)MXPCVR, MXPCVR, MXPCVR 104 FORMAT(1X,'INDEXABLE SIZE (CELLS) = ',I4,' X ',I4,' X ',I4) WRITE(6,105)MXDMRT 105 FORMAT(1X,'ADJUST IF RATIO EXCEEDS = ',F8.3) WRITE(6,106)AXL, AXH 106 FORMAT(1X,'SCENE X RANGE = (',F8.3,':',F8.3,')') WRITE(6,107)AYL, AYH 107 FORMAT(1X,'SCENE Y RANGE = (',F8.3,':',F8.3,')') WRITE(6,108)AZL, AZH 108 FORMAT(1X,'SCENE Z RANGE = (',F8.3,':',F8.3,')') C C-- OPTIONALLY, CHOOSE THE BIGGEST AXIS TO SCALE THE WHOLE GRID C-- IF THERE IS A BIG DISPARITY IN THE DIMENSIONS OF THE AXES, THIS CAN HELP C-- ENSURE THAT CELLS ARE MADE BIG ENOUGH TO STILL INDEX PRIMITIVES. IF THE CELLS C-- ARE MADE VERY SMALL IN ONE DIRECTION, IT MAY BE THAT TOO MANY ARE COVERED C-- OTHERWISE, AND THE ACCELERATION SCHEME CAN'T BE USED. C-- IT SEEMS LIKE THIS IS ALMOST ALWAYS A GOOD IDEA ANYWAY. C MAXDIR = ADXS IF( ADYS .GT. MAXDIR )MAXDIR = ADYS IF( ADZS .GT. MAXDIR )MAXDIR = ADZS C MINDIR = ADXS IF( ADYS .LT. MINDIR )MINDIR = ADYS IF( ADZS .LT. MINDIR )MINDIR = ADZS C C--TUNE THIS FOR WHEN TO DO IT... C IF( ( MAXDIR / MINDIR ) .GT. MXDMRT )THEN XC = 0.5 * ( AXL + AXH ) YC = 0.5 * ( AYL + AYH ) ZC = 0.5 * ( AZL + AZH ) AXL = XC - 0.5 * MAXDIR AXH = XC + 0.5 * MAXDIR AYL = YC - 0.5 * MAXDIR AYH = YC + 0.5 * MAXDIR AZL = ZC - 0.5 * MAXDIR AZH = ZC + 0.5 * MAXDIR ADXS = MAXDIR ADYS = MAXDIR ADZS = MAXDIR WRITE(6,109) 109 FORMAT(1X,'SCENE RANGES ADJUSTED FOR DIMENSIONAL UNIFORMITY') WRITE(6,110)AXL, AXH 110 FORMAT(1X,'SCENE X RANGE (ADJ) = (',F8.3,':',F8.3,')') WRITE(6,111)AYL, AYH 111 FORMAT(1X,'SCENE Y RANGE (ADJ) = (',F8.3,':',F8.3,')') WRITE(6,112)AZL, AZH 112 FORMAT(1X,'SCENE Z RANGE (ADJ) = (',F8.3,':',F8.3,')') ENDIF C C-- CALCULATE THE SIZE OF EACH CELL IN EYE SPACE C DXC = ADXS / ANX DYC = ADYS / ANY DZC = ADZS / ANZ WRITE(6,113)DXC, DYC, DZC 113 FORMAT(1X,'CELL SIZES = ',F8.3,' X',F8.3,' X',F8.3) C C-- GO OVER THE PRIMITIVES AGAIN, ASSIGNING THEM TO CELLS. C-- IF THE PRIMITIVE BOUNDS COVER MORE THAN MXPCVR CELLS IN ANY DIMENSION, WE C-- DON'T USE ACCELERATION FOR THAT PRIMITIVE. THIS IS TO LIMIT THE SIZE OF C-- THE CELL LISTS. C DO 5 P=1,NPRM PIDX = PRMIDX(P) IF( PRMTYP(P) .EQ. SPHPRM )THEN CALL SPHBND( PIDX, PXL, PXH, PYL, PYH, PZL, PZH ) ELSE IF( PRMTYP(P) .EQ. TRIPRM )THEN CALL TRIBND( PIDX, PXL, PXH, PYL, PYH, PZL, PZH ) ELSE WRITE(6,100) WRITE(6,101) STOP ENDIF CXL = ( PXL - AXL ) / DXC CXH = ( PXH - AXL ) / DXC CYL = ( PYL - AYL ) / DYC CYH = ( PYH - AYL ) / DYC CZL = ( PZL - AZL ) / DZC CZH = ( PZH - AZL ) / DZC C C-- MOVE TO INTEGER GRID CELL COORDINATES. C IXL = INT(CXL) IXH = INT(CXH) IYL = INT(CYL) IYH = INT(CYH) IZL = INT(CZL) IZH = INT(CZH) C C-- ADD PRIMITIVE TO THE LIST OF PRIMITIVES FOR EACH CELL COVERED BY THE BOUNDS OF THE PRIMITIVE. C-- IF THE PRIMITIVE IS TOO BIG, DON'T USE THE ACCELERATION STRUCTURE AT ALL. C-- THIS IS TO LIMIT THE SIZE OF THE PRIMITIVE LISTS IN NASTY CASES. C IF( IXL .EQ. IXH .AND. IYL .EQ. IYH .AND. IZL .EQ. IZH )THEN PACCEL(P) = 1 NACC = NACC + 1 CALL ADDLST( IXL, IYL, IZL, P ) ELSE IF( ( IXH - IXL ) .GE. MXPCVR .OR. + ( IYH - IYL ) .GE. MXPCVR .OR. + ( IZH - IZL ) .GE. MXPCVR )THEN PACCEL(P) = 0 NOACC = NOACC + 1 ELSE PACCEL(P) = 1 NACC = NACC + 1 DO 6 K=IZL,IZH DO 7 J=IYL,IYH DO 8 I=IXL,IXH CALL ADDLST( I, J, K, P ) 8 CONTINUE 7 CONTINUE 6 CONTINUE ENDIF ENDIF 5 CONTINUE C WRITE(6,114)NACC 114 FORMAT(1X,'PRIMITIVES INDEXED = ',I5) WRITE(6,115)NOACC 115 FORMAT(1X,'PRIMITIVES NOT INDEXED = ',I5) C RETURN END C C SUBROUTINE ADDLST( I, J, K, P ) IMPLICIT CHARACTER*1 (A-Z) INTEGER I, J, K, P C********************************************************************** C ADD PRIMITIVE P TO THE LIST ASSOCIATED WITH CELL (I,J,K) IF IT ISN'T C ALREADY ON THE LIST. *** N.B.: (I,J,K) ARE ZERO BASED! *** C********************************************************************** *CALL PARAMS *CALL ACCEL INTEGER CI, LP, LC C CI = I + J * ANX + K * ( ANX * ANY ) + 1 IF( CELLS(CI) .EQ. 0 )THEN NLIST = NLIST + 1 IF( NLIST .GT. MAXLST )THEN WRITE(6,100) 100 FORMAT(1X,'ADDLST: TRYING TO ADD TOO MANY LIST ELEMENTS.') ELSE CELLS(CI) = NLIST LPRIM(NLIST) = P LNEXT(NLIST) = 0 ENDIF ELSE LP = 0 LC = CELLS(CI) 1 CONTINUE IF( LC .EQ. 0 )GOTO 2 IF( LPRIM(LC) .EQ. P )THEN RETURN ELSE LP = LC LC = LNEXT(LC) ENDIF GOTO 1 2 CONTINUE NLIST = NLIST + 1 IF( NLIST .GT. MAXLST )THEN WRITE(6,100) ELSE LNEXT(LP) = NLIST LPRIM(NLIST) = P LNEXT(NLIST) = 0 ENDIF ENDIF C RETURN END C C SUBROUTINE SPHBND( I, PXL, PXH, PYL, PYH, PZL, PZH ) IMPLICIT CHARACTER*1 (A-Z) INTEGER I REAL PXL, PXH, PYL, PYH, PZL, PZH C********************************************************************** C FIND A BOUNDING BOX FOR A SPHERE PRIMITIVE. VERY EASY! C********************************************************************** *CALL SPHLST PXL = SOX(I) - SRAD(I) PXH = SOX(I) + SRAD(I) PYL = SOY(I) - SRAD(I) PYH = SOY(I) + SRAD(I) PZL = SOZ(I) - SRAD(I) PZH = SOZ(I) + SRAD(I) C RETURN END C C SUBROUTINE TRIBND( I, PXL, PXH, PYL, PYH, PZL, PZH ) IMPLICIT CHARACTER*1 (A-Z) INTEGER I REAL PXL, PXH, PYL, PYH, PZL, PZH C********************************************************************** C FIND A BOUNDING BOX FOR A TRIANGLE PRIMITIVE. VERY EASY! C********************************************************************** *CALL TRILST PXL = TX1(I) IF( TX2(I) .LT. PXL )PXL = TX2(I) IF( TX3(I) .LT. PXL )PXL = TX3(I) C PYL = TY1(I) IF( TY2(I) .LT. PYL )PYL = TY2(I) IF( TY3(I) .LT. PYL )PYL = TY3(I) C PZL = TZ1(I) IF( TZ2(I) .LT. PZL )PZL = TZ2(I) IF( TZ3(I) .LT. PZL )PZL = TZ3(I) C PXH = TX1(I) IF( TX2(I) .GT. PXH )PXH = TX2(I) IF( TX3(I) .GT. PXH )PXH = TX3(I) C PYH = TY1(I) IF( TY2(I) .GT. PYH )PYH = TY2(I) IF( TY3(I) .GT. PYH )PYH = TY3(I) C PZH = TZ1(I) IF( TZ2(I) .GT. PZH )PZH = TZ2(I) IF( TZ3(I) .GT. PZH )PZH = TZ3(I) C RETURN END C C SUBROUTINE GRDTRV( UX, UY, UZ, VX, VY, VZ, NX, NY, NZ ) IMPLICIT CHARACTER*1 (A-Z) REAL UX, UY, UZ, VX, VY, VZ INTEGER NX, NY, NZ C********************************************************************** C VISIT ALL CELLS IN A 3D INTEGER GRID OF NX X NY X NZ CELLS THAT ARE C PIERCED BY A RAY WITH ORIGIN (UX,UY,UZ) - **POSITIVE AND WITHIN THE C GRID** - AND DIRECTION VECTOR (VX,VY,VZ) - NOT NECESSARILY NORMALIZED. C *** N.B.: CELL COORDINATES PASSED TO VISCEL() ARE ZERO BASED. *** C THIS IS AMANTIDES & WOO'S ALGORITHM. C********************************************************************** *CALL PARAMS *CALL ACCEL INTEGER X, Y, Z, SX, SY, SZ, JX, JY, JZ, QUIT REAL TX, TY, TZ, DX, DY, DZ C C-- INITIAL CELL COORDINATES C X = INT(UX) Y = INT(UY) Z = INT(UZ) C C-- VISIT THE INITIAL CELL C CALL VISCEL( X, Y, Z, QUIT ) IF( QUIT .NE. 0 )RETURN C C-- INTEGER CELL STEPS FOR THE VECTOR C IF( VX .LT. 0.0 )THEN SX = -1 ELSE SX = 1 ENDIF IF( VY .LT. 0.0 )THEN SY = -1 ELSE SY = 1 ENDIF IF( VZ .LT. 0.0 )THEN SZ = -1 ELSE SZ = 1 ENDIF C C-- FIND TX,TY,TZ - THE T VALUES THAT MOVE ALONG THE RAY TO THE BOUNDARY C-- OF THE INITIAL CELL. C-- FIND DX,DY,DZ - THE T VALUES THAT STEP THE WIDTH, HEIGHT AND DEPTH C-- OF A FULL CELL. C-- PRECISION! C IF( ABS( VX ) .LT. EPS )THEN TX = GTHUGE DX = GTHUGE ELSE IF( VX .GT. 0 )THEN TX = ( ( X + 1 ) - UX ) / VX DX = 1.0 / VX ELSE TX = ( X - UX ) / VX DX = -1.0 / VX ENDIF IF( ABS( VY ) .LT. EPS )THEN TY = GTHUGE DY = GTHUGE ELSE IF( VY .GT. 0 )THEN TY = ( ( Y + 1 ) - UY ) / VY DY = 1.0 / VY ELSE TY = ( Y - UY ) / VY DY = -1.0 / VY ENDIF IF( ABS( VZ ) .LT. EPS )THEN TZ = GTHUGE DZ = GTHUGE ELSE IF( VZ .GT. 0 )THEN TZ = ( ( Z + 1 ) - UZ ) / VZ DZ = 1.0 / VZ ELSE TZ = ( Z - UZ ) / VZ DZ = -1.0 / VZ ENDIF C C-- IF ANY OF THE T VALUES TO MOVE TO THE BOUNDARY ARE EXACTLY ZERO, C-- BUMP THEM A BIT. C-- PRECISION! C IF( TX .EQ. 0 )TX = TX + EPS IF( TY .EQ. 0 )TY = TY + EPS IF( TZ .EQ. 0 )TZ = TZ + EPS C C-- FIND JUST OUT OF BOUNDS CELL COORDS FOR THE RAY ALONG EACH AXIS C IF( SX .GT. 0 )THEN JX = NX ELSE JX = -1 ENDIF IF( SY .GT. 0 )THEN JY = NY ELSE JY = -1 ENDIF IF( SZ .GT. 0 )THEN JZ = NZ ELSE JZ = -1 ENDIF C C-- STEP THROUGH THE CELLS UNTIL WE EXIT THE ARRAY C 1 CONTINUE IF( TX .LT. TY )THEN IF( TX .LT. TZ )THEN X = X + SX IF( X .EQ. JX )RETURN TX = TX + DX ELSE Z = Z + SZ IF( Z .EQ. JZ )RETURN TZ = TZ + DZ ENDIF ELSE IF( TY .LT. TZ )THEN Y = Y + SY IF( Y .EQ. JY )RETURN TY = TY + DY ELSE Z = Z + SZ IF( Z .EQ. JZ )RETURN TZ = TZ + DZ ENDIF ENDIF C C-- VISIT THE CURRENT CELL C CALL VISCEL( X, Y, Z, QUIT ) IF( QUIT .NE. 0 )RETURN GOTO 1 C END C C SUBROUTINE BOXINS( OX, OY, OZ, DX, DY, DZ, SECT, XI, YI, ZI ) IMPLICIT CHARACTER*1 (A-Z) REAL OX, OY, OZ, DX, DY, DZ, XI, YI, ZI INTEGER SECT C********************************************************************** C SEE IF THE RAY WITH ORIGIN (OX,OY,OZ) DIRECTION (DX, DY, DZ) C (NOT NECESSARILY NORMALISED) INTERSECTS THE BOUNDING C BOX FOR THE SCENE GEOMETRIC OBJECTS FOUND BY STACEL(). IF SO, RETURN C SECT SET TO 1 AND THE INTERSECTION POINT IN C (XI, YI, ZI). OTHERWISE, RETURN SECT 0. C********************************************************************** *CALL PARAMS *CALL ACCEL REAL T1, T2, TNEAR, TFAR, T C XI = 0.0 YI = 0.0 ZI = 0.0 SECT = 0 TNEAR = MINFLT TFAR = MAXFLT C C-- PRECISION! C IF( ABS( DX ) .LT. EPS )THEN IF( OX .LT. AXL .OR. OX .GT. AXH )THEN SECT = 0 RETURN ENDIF ELSE TNEAR = ( AXL - OX ) / DX TFAR = ( AXH - OX ) / DX IF( TNEAR .GT. TFAR )THEN T = TNEAR TNEAR = TFAR TFAR = T ENDIF IF( TFAR .LT. 0.0 )THEN SECT = 0 RETURN ENDIF ENDIF C C-- PRECISION! C IF( ABS( DY ) .LT. EPS )THEN IF( OY .LT. AYL .OR. OY .GT. AYH )THEN SECT = 0 RETURN ENDIF ELSE T1 = ( AYL - OY ) / DY T2 = ( AYH - OY ) / DY IF( T1 .GT. T2 )THEN T = T1 T1 = T2 T2 = T ENDIF IF( T1 .GT. TNEAR )TNEAR = T1 IF( T2 .LT. TFAR )TFAR = T2 IF( TNEAR .GT. TFAR )THEN SECT = 0 RETURN ENDIF IF( TFAR .LT. 0.0 )THEN SECT = 0 RETURN ENDIF ENDIF C C-- PRECISION! C IF( ABS( DZ ) .LT. EPS )THEN IF( OZ .LT. AZL .OR. OZ .GT. AZH )THEN SECT = 0 RETURN ENDIF ELSE T1 = ( AZL - OZ ) / DZ T2 = ( AZH - OZ ) / DZ IF( T1 .GT. T2 )THEN T = T1 T1 = T2 T2 = T ENDIF IF( T1 .GT. TNEAR )TNEAR = T1 IF( T2 .LT. TFAR )TFAR = T2 IF( TNEAR .GT. TFAR )THEN SECT = 0 RETURN ENDIF IF( TFAR .LT. 0.0 )THEN SECT = 0 RETURN ENDIF ENDIF C C-- THERE IS AN INTERSECTION. FIND WHERE. C XI = OX + TNEAR * DX YI = OY + TNEAR * DY ZI = OZ + TNEAR * DZ SECT = 1 C RETURN END C C---------------------------------------------------------------------- C UTILITY ROUTINES C---------------------------------------------------------------------- C REAL FUNCTION RANHLF( IDUM ) IMPLICIT CHARACTER*1 (A-Z) INTEGER IDUM C********************************************************************** C RETURN A RANDOM NUMBER IN THE RANGE +/- 0.5 C********************************************************************** RANHLF = RANF() - 0.5 RETURN END C C SUBROUTINE NRMVEC( DX, DY, DZ ) IMPLICIT CHARACTER*1 (A-Z) REAL DX, DY, DZ C********************************************************************** C MAKE (DX,DY,DZ) A UNIT VECTOR C********************************************************************** *CALL PARAMS REAL LEN C LEN = SQRT( DX*DX + DY*DY + DZ*DZ ) C C-- PRECISION! C IF( LEN .GT. EPS )THEN DX = DX / LEN DY = DY / LEN DZ = DZ / LEN ENDIF C RETURN END }