DECLARE SUB precalc ()
DECLARE SUB setpal2 ()
DECLARE FUNCTION intersect% (c1!, c0!)

DECLARE FUNCTION totaldirectlight! (V0 AS ANY, V AS ANY, N AS ANY)
DECLARE FUNCTION phong! (L AS ANY, V AS ANY, N AS ANY)
DECLARE SUB setpal ()
DECLARE SUB raytrace (V0 AS ANY, dv AS ANY, c!, depth!)
DECLARE FUNCTION magnitude! (nix!, niy!, niz!)
DECLARE SUB normalize (x!, y!, z!)
DECLARE SUB rayscan ()
'==================================================================
'RAYTRACE 5   Copyright 2000 Toshihiro Horie
'  last modified: 6/22/2000
'  released: 7/11/2000
'  Download other cool GFX demos at 
'         http://www.ocf.berkeley.edu/~horie/project.html
'  and at http://toshi.tekscode.com
'  The scene consists of:
'   * several spheres and point light sources
'==================================================================
CONST FALSE = 0
CONST TRUE = NOT FALSE
TYPE rgbtype
    red AS INTEGER
    green AS INTEGER
    blue AS INTEGER
END TYPE
TYPE vec
    x AS SINGLE
    y AS SINGLE
    z AS SINGLE
END TYPE
TYPE sph
    a AS SINGLE
    b AS SINGLE
    c AS SINGLE
    clr AS INTEGER
    r AS SINGLE
    r2 AS SINGLE
    da AS SINGLE
    db AS SINGLE
    dc AS SINGLE
    rn AS SINGLE
END TYPE
TYPE pln
    a AS SINGLE
    b AS SINGLE
    c AS SINGLE
    d AS SINGLE
END TYPE
'center of projection (aka. the eye's position)
COMMON SHARED eye AS vec
eye.x = 160: eye.y = 100: eye.z = -160
CONST numlights = 2 'actually, there is one more light than specified here
CONST numspheres = 8
DIM SHARED o(numspheres) AS sph
DIM SHARED plane(0) AS pln
DIM SHARED lights(numlights) AS vec
DIM SHARED dx2%(-160 TO 160)
precalc

FOR i% = 0 TO numlights
    lights(i%).x = eye.x - 20 + (160 * i%)
    lights(i%).y = eye.y - 120 + (190 * i%) * -(i% AND 1)
    lights(i%).z = eye.z
NEXT i%
    
'y=-40
'plane eq: Ax+By+Cz=-D
plane(0).a = 0
plane(0).b = 1
plane(0).c = 0
plane(0).d = 0

RANDOMIZE TIMER

SCREEN 13


setpal2 'for grayscale
'setpal 'for false color
DO
    'create colored spheres
    FOR i% = 0 TO numspheres
        o(i%).a = RND * 320
        o(i%).b = RND * 200
        o(i%).c = i% * 10 + 65
        o(i%).clr = i% * 10 + 2
        o(i%).r = RND * 50 + 8
        o(i%).r2 = o(i%).r * o(i%).r
        o(i%).da = eye.x - o(i%).a
        o(i%).db = eye.y - o(i%).b
        o(i%).dc = eye.z - o(i%).c
        o(i%).rn = (o(i%).da * o(i%).da) + (o(i%).db * o(i%).db) + (o(i%).dc * o(i%).dc) - o(i%).r * o(i%).r
    NEXT i%
    CLS
    CALL rayscan
    t1 = TIMER
    DO: k$ = INKEY$: LOOP WHILE TIMER - t1 < 2 AND k$ = ""
LOOP WHILE k$ = ""
END

FUNCTION intersect% (c1, c0)
'intersect% reports the number of distinct, real valued >1 intercepts
'INP: c1 and c0 are parameters for quadratic equation (A)*t^2+(C1)*t+C0 = 0
'OUT: c1 is used to return the parameter t values for the closest intercepts
discr = c1 * c1 - 4 * c0
IF discr < 0 THEN
    intersect% = 0
ELSEIF discr > 0 THEN
    sd = SQR(discr) 'always sd >=0
    c0 = .5 * (-c1 + sd)'order  !
    c1 = .5 * (-c1 - sd)'matters! because overwrites c1
   
    IF c0 <= 1 AND c1 <= 1 THEN intersect% = 0: EXIT FUNCTION
   
    'store closer intersection in c1, farther in c0
    IF c1 <= 1 THEN c1 = c0
    intersect% = 2
ELSEIF discr = 0 THEN
    intersect% = 1
    'c0 = -c1
    c1 = -c1
END IF
END FUNCTION

FUNCTION magnitude (nix, niy, niz)
magnitude = SQR(nix * nix + niy * niy + niz * niz)
END FUNCTION

SUB normalize (x, y, z)
mag = SQR(x * x + y * y + z * z)
'should make this beep in debug mode
IF mag = 0 THEN EXIT SUB
x = x / mag
y = y / mag
z = z / mag
END SUB

FUNCTION phong (L AS vec, V AS vec, N AS vec)
    DIM r AS vec
    fac = .8 / numlights 'the object absorbs 20% of the light that shines on it.
    m = 82  'specularity coefficient (higher is more specular (and shinier) )
    'R is reflection of L about N
    'R = 2(N dot L)*N - L
    b = (N.x * L.x + N.y * L.y + N.z * L.z)
    ndotL2 = b * 2
    r.x = ndotL2 * N.x - L.x
    r.y = ndotL2 * N.y - L.y
    r.z = ndotL2 * N.z - L.z

    'brightness is proportional to (R dot V)^m
     rdotv = (r.x * V.x + r.y * V.y + r.z * V.z)
     IF rdotv < 0 THEN
        'rdotv=0
        phong = (b * 125 + 2) * fac: EXIT FUNCTION
     ELSE
        'use n=8 for a somewhat shiny surface
        'Brightness = specular + diffuse + ambient
        IF rdotv > 1 THEN 'this was causing overflow
            phong = (127 + b * 125 + 2) * fac: EXIT FUNCTION
        ELSE
            phong = (rdotv ^ m * 127 + b * 125 + 2) * fac
        END IF
     END IF
END FUNCTION

SUB precalc

'nothing to precalc!

FOR i% = -160 TO 160
dx2%(i%) = i% * i%
NEXT i%

END SUB

SUB rayscan
'initialized right before raytrace is called
'finally seems to be working correctly! 6/20
DIM dv AS vec: DIM V0 AS vec
DIM N AS vec
DIM V AS vec
t1 = TIMER
'------------interlaced field 1-------------------------
'DEF SEG=&HA000
dz = z - eye.z
dz2 = dz * dz
clr = 0
    FOR y% = 20 TO 179 STEP 2
       
        dy = y% - eye.y
        dy2 = dy * dy
        dyZ2 = dy2 + dz2
        IF INKEY$ <> "" THEN EXIT SUB
        FOR x% = 0 TO 319
            dx = x% - eye.x
            isdr2 = 1! / SQR(dx2%(dx) + dyZ2)
            dxn = dx * isdr2
            dyn = dy * isdr2
            dzn = dz * isdr2
          
            closestobj% = -1
            neart = 1E+09
            FOR i% = 0 TO numspheres
  
    'a is the quadratic, b is the linear, and c is the constant term
    'o(i) is the center of the sphere we're testing against
               
                'a = 1
                b = 2 * (dxn * o(i%).da + dyn * o(i%).db + dzn * o(i%).dc)
                c = o(i%).rn
               
                IF intersect(b, c) THEN
                    IF b < neart THEN
                        neart = b
                        closestobj% = i%
                    END IF
                END IF
            NEXT i%
            IF closestobj% <> -1 THEN
                'LOCATE 1, 1: COLOR 255: PRINT neart
                 'difference between eye and collision point
                 dxi = -neart * dxn
                 dyi = -neart * dyn
                 dzi = -neart * dzn

                 'position of collision in world coordinates
                 V0.x = eye.x - dxi
                 V0.y = eye.y - dyi
                 V0.z = eye.z - dzi
                 'V=(dxi,dyi,dzi) is the vector from sphere pt. back to eye
                 'new V0=(xi,yi,zi) is point of i'sect of ray on sphere's surface
                 'new ray origin is at collision point on surface of sphere
                 'COLOR 255: PRINT USING "delta=##.## ##.## ##.## ctr=##.## ##.## ##.##"; o(closestobj%).a; o(closestobj%).b; o(closestobj%).c

                'THESE TWO SHOULD BE EQUAL!
                ' COLOR 255: PRINT USING "###.#"; SQR((V0.x - o(closestobj%).a) ^ 2 + (V0.y - o(closestobj%).b) ^ 2 + (V0.z - o(closestobj%).c) ^ 2);
                ' PRINT USING "###.#"; o(closestobj%).r2
               
                'N = normal of the sphere at the collision point
                'N is normalized by the 1/r operation.
                 N.x = (V0.x - o(closestobj%).a) / o(closestobj%).r
                 N.y = (V0.y - o(closestobj%).b) / o(closestobj%).r
                 N.z = (V0.z - o(closestobj%).c) / o(closestobj%).r
               
                 'New direction dv = reflection of V about N.
                 V.x = dxi: V.y = dyi: V.z = dzi: CALL normalize(V.x, V.y, V.z)
                 ndoti2 = 2 * (N.x * V.x + N.y * V.y + N.z * V.z)
                 dv.x = ndoti2 * N.x - V.x
                 dv.y = ndoti2 * N.y - V.y
                 dv.z = ndoti2 * N.z - V.z
                 'the new dv is normalized because symmetric
                 'reflections preserve length.

                 'add direct light contributions
                 tclr = totaldirectlight(V0, V, N)
                 clr = 0
                 IF tclr < 255 THEN
                    'recursively raytrace
                    depth = 0 'initialize ray tracer
                    CALL raytrace(V0, dv, clr, depth)
                 END IF
                 clr = clr + tclr
                 IF clr > 255 THEN clr = 255
                 'Yes, you can use POKE here, but be careful!
                    PSET (x%, y%), clr
                    PSET (x%, y% + 1), clr
            END IF
        NEXT
    NEXT



'-----------interlaced field 2--------------------------
'DEF SEG=&HA000
dz = z - eye.z
dz2 = dz * dz
clr = 0
    FOR y% = 21 TO 179 STEP 2
       
        dy = y% - eye.y
        dy2 = dy * dy
        dyZ2 = dy2 + dz2
        IF INKEY$ <> "" THEN EXIT SUB
        FOR x% = 0 TO 319
            dx = x% - eye.x
            isdr2 = 1! / SQR(dx * dx + dyZ2)
            dxn = dx * isdr2
            dyn = dy * isdr2
            dzn = dz * isdr2
          
            closestobj% = -1
            neart = 1E+09
            FOR i% = 0 TO numspheres
  
    'a is the quadratic, b is the linear, and c is the constant term
    'o(i) is the center of the sphere we're testing against
               
                'a = 1
                b = 2 * (dxn * o(i%).da + dyn * o(i%).db + dzn * o(i%).dc)
                c = o(i%).rn
               
                IF intersect(b, c) THEN
                    IF b < neart THEN
                        neart = b
                        closestobj% = i%
                    END IF
                END IF
            NEXT i%
            IF closestobj% <> -1 THEN
                'LOCATE 1, 1: COLOR 255: PRINT neart
                 'difference between eye and collision point
                 dxi = -neart * dxn
                 dyi = -neart * dyn
                 dzi = -neart * dzn

                 'position of collision in world coordinates
                 V0.x = eye.x - dxi
                 V0.y = eye.y - dyi
                 V0.z = eye.z - dzi
                 'V=(dxi,dyi,dzi) is the vector from sphere pt. back to eye
                 'new V0=(xi,yi,zi) is point of i'sect of ray on sphere's surface
                 'new ray origin is at collision point on surface of sphere
                 'COLOR 255: PRINT USING "delta=##.## ##.## ##.## ctr=##.## ##.## ##.##"; o(closestobj%).a; o(closestobj%).b; o(closestobj%).c

                'THESE TWO SHOULD BE EQUAL!
                ' COLOR 255: PRINT USING "###.#"; SQR((V0.x - o(closestobj%).a) ^ 2 + (V0.y - o(closestobj%).b) ^ 2 + (V0.z - o(closestobj%).c) ^ 2);
                ' PRINT USING "###.#"; o(closestobj%).r2
               
                'N = normal of the sphere at the collision point
                'N is normalized by the 1/r operation.
                 N.x = (V0.x - o(closestobj%).a) / o(closestobj%).r
                 N.y = (V0.y - o(closestobj%).b) / o(closestobj%).r
                 N.z = (V0.z - o(closestobj%).c) / o(closestobj%).r
               
                 'New direction dv = reflection of V about N.
                 V.x = dxi: V.y = dyi: V.z = dzi: CALL normalize(V.x, V.y, V.z)
                 ndoti2 = 2 * (N.x * V.x + N.y * V.y + N.z * V.z)
                 dv.x = ndoti2 * N.x - V.x
                 dv.y = ndoti2 * N.y - V.y
                 dv.z = ndoti2 * N.z - V.z
                 'the new dv is normalized already because
                 'symmetric reflections preserve length.

                 'add direct light contributions
                 tclr = totaldirectlight(V0, V, N)
                 clr = 0
                 IF tclr < 255 THEN
                    'recursively raytrace
                    depth = 0 'initialize ray tracer
                    CALL raytrace(V0, dv, clr, depth)
                 END IF
                 clr = clr + tclr
                 IF clr > 255 THEN clr = 255
                 'Yes, you can use POKE here, but be careful!
                 PSET (x%, y%), clr
            END IF
        NEXT
    NEXT

LOCATE 1, 1: COLOR 255: PRINT USING "##.##s  "; TIMER - t1
END SUB

SUB raytrace (V0 AS vec, dv AS vec, c, depth)
IF c >= 255 THEN EXIT SUB
'this is not quite working correctly
DIM V AS vec
DIM N AS vec

'PRE: v0 is new "eye" point for ray tracing
'     dv is normalized direction vector (direction eye is looking)
'base case
IF depth > 5 THEN EXIT SUB

'normal case
closestobj% = -1
neart = 9999
'check for intersection with all spheres
FOR i% = 0 TO numspheres
    diffx = V0.x - o(i%).a
    diffy = V0.y - o(i%).b
    diffz = V0.z - o(i%).c
    c1 = 2 * (dv.x * diffx + dv.y * diffy + dv.z * diffz)
    c0 = diffx * diffx + diffy * diffy + diffz * diffz - o(i%).r2
    
    IF intersect(c1, c0) THEN
        'this c1 is actually the t value at the intersection
        'between the sphere and our ray
        IF c1 < neart THEN
            neart = c1
            closestobj% = i%
        END IF
    END IF
NEXT i%
'check for intersection with all planes (polygons)
'FOR i% = 0 TO 0
'    denom = plane(i%).a * dv.x + plane(i%).b * dv.y + plane(i%).c * dv.z
'    ' if denominator is really close to zero, intersect far away
'    IF ABS(denom) > .0001 THEN
'        numer = -plane(i%).a * v0.x + plane(i%).b * v0.y + plane(i%).c * v0.z
'        t1 = numer / denom
'        IF t1 < neart AND t1 > 0 THEN
'            z = v0.z + t * dv.z
'            x = v0.x + t * dv.x
'            'c = c + INT(z)
'            c2 = c
'            IF z > -230 AND z < 220 AND x > 0 AND x < 130 THEN
'            IF (INT(z) AND 8) < 4 THEN
'                c = 32 * -((INT(x) AND 8) < 4)
'            ELSE
'                c = 0
'            END IF
'            END IF
'            c = c2 + c
'            'IF c > 255 THEN c = 255
'            EXIT SUB
'        END IF
'    END IF
'NEXT i%




IF closestobj% <> -1 THEN
    'LOCATE 1, 1: COLOR 255: PRINT neart
     'difference between eye and collision point
     dxi = -neart * dv.x
     dyi = -neart * dv.y
     dzi = -neart * dv.z

     'position of collision in world coordinates
     V0.x = V0.x - dxi
     V0.y = V0.y - dyi
     V0.z = V0.z - dzi
     'V=(dxi,dyi,dzi) is the vector from sphere pt. back to eye
     'new V0=(xi,yi,zi) is point of i'sect of ray on sphere's surface
     'new ray origin is at collision point on surface of sphere
'     COLOR 255: PRINT USING "###.#"; SQR((V0.x - o(closestobj%).a) ^ 2 + (V0.y - o(closestobj%).b) ^ 2 + (V0.z - o(closestobj%).c) ^ 2);
'     PRINT USING "###.#"; o(closestobj%).r
   
    'N = normal of the sphere at the collision point
     N.x = (V0.x - o(closestobj%).a) / o(closestobj%).r
     N.y = (V0.y - o(closestobj%).b) / o(closestobj%).r
     N.z = (V0.z - o(closestobj%).c) / o(closestobj%).r
               
     'New direction dv = reflection of V about N.
     V.x = dxi: V.y = dyi: V.z = dzi: CALL normalize(V.x, V.y, V.z)
     ndoti2 = 2 * (N.x * V.x + N.y * V.y + N.z * V.z)
     dv.x = ndoti2 * N.x - V.x
     dv.y = ndoti2 * N.y - V.y
     dv.z = ndoti2 * N.z - V.z
     tclr = totaldirectlight(V0, V, N)
     c = c + tclr
     IF c > 255 THEN EXIT SUB
     depth = depth + 1
     CALL raytrace(V0, dv, c, depth)
END IF


END SUB

SUB setpal
' ported from the C code in the book Computer Graphics by Foley et al.
numcolors = 256
max# = 63
DIM HSVPal(numcolors) AS rgbtype
SCREEN 13
S# = 0
S# = .9: 'saturation (0 to 1) 0=greyscale, 1=marker colors
V# = .9: 'brightness (0 to 1) 0=black, 1=brightest
p# = V# * (1# - S#)
conV# = numcolors / 360#
'if H was undefined, then we would have the greyscale.
FOR H% = 0 TO 359:       'hue (0 to 360) aka rainbow spectra
    c% = INT(H% * conV#): 'index into palette array
    Hn# = H% / 60        'Hn is in the interval [0,6)
    i% = INT(Hn#)        'i=floor(Hn)
    f# = Hn# - i%        'fractional part of Hn
    q# = V# * (1# - (S# * f#))
    t# = V# * (1# - (S# * (1 - f#)))
    'PRINT i%, c%, q#, t#
    SELECT CASE i%
        CASE 0:
            HSVPal(c%).red = INT(V# * max#)
            HSVPal(c%).green = INT(t# * max#)
            HSVPal(c%).blue = INT(p# * max#)
        CASE 1:
            HSVPal(c%).red = INT(q# * max#)
            HSVPal(c%).green = INT(V# * max#)
            HSVPal(c%).blue = INT(p# * max#)
        CASE 2:
            HSVPal(c%).red = p# * max#
            HSVPal(c%).green = V# * max#
            HSVPal(c%).blue = t# * max#
        CASE 3:
            HSVPal(c%).red = INT(p# * max#)
            HSVPal(c%).green = INT(q# * max#)
            HSVPal(c%).blue = INT(V# * max#)
        CASE 4:
            HSVPal(c%).red = INT(t# * max#)
            HSVPal(c%).green = INT(p# * max#)
            HSVPal(c%).blue = INT(V# * max#)
        CASE 5:
            HSVPal(c%).red = INT(V# * max#)
            HSVPal(c%).green = INT(p# * max#)
            HSVPal(c%).blue = INT(q# * max#)
    END SELECT
NEXT H%
OUT &H3C8, 0: OUT &H3C9, 0: OUT &H3C9, 0: OUT &H3C9, 0
FOR i% = 1 TO numcolors - 1
    OUT &H3C8, i%
    OUT &H3C9, HSVPal(i%).red
    OUT &H3C9, HSVPal(i%).green
    OUT &H3C9, HSVPal(i%).blue
NEXT i%
OUT &H3C8, 255: OUT &H3C9, 255: OUT &H3C9, 255: OUT &H3C9, 255
END SUB

SUB setpal2
OUT &H3C8, 0
FOR i% = 0 TO 255
gray% = i% \ 4
OUT &H3C9, gray%: OUT &H3C9, gray%: OUT &H3C9, gray%
NEXT
END SUB

FUNCTION totaldirectlight (V0 AS vec, V AS vec, N AS vec)
'V is normalized view vector
'v0 is new eye point (point of reflection)
DIM r AS vec
DIM L AS vec
DIM dv AS vec
'for each light
clr = 0
FOR m% = 0 TO numlights
    blocked = FALSE
    'L has the direction toward light from sphere
    L.x = lights(m%).x - V0.x
    L.y = lights(m%).y - V0.y
    L.z = lights(m%).z - V0.z
    'can get divide by zero or overflow if light is at reflection point.
    Lmag2 = L.x * L.x + L.y * L.y + L.z * L.z
    Linvmag2 = 1! / Lmag2
    lightdist = SQR(Lmag2)
    
    neart = 9.9E+09
    'check for intersection with all spheres
    FOR i% = 0 TO numspheres
        diffx = V0.x - o(i%).a
        diffy = V0.y - o(i%).b
        diffz = V0.z - o(i%).c
        c1 = 2 / lightdist * (L.x * diffx + L.y * diffy + L.z * diffz)
        c0 = diffx * diffx + diffy * diffy + diffz * diffz - o(i%).r2
        
        IF intersect(c1, c0) THEN
            'this c1 is actually the t value at the intersection
            'between the sphere and our ray
            IF c1 < neart THEN
                neart = c1
                'short circuit if light blocked by any sphere
                ' Can compare Lrealmag to neart because both neart
                '   and Lrealmag are in world coordinate units
                IF lightdist >= neart THEN blocked = TRUE: EXIT FOR
            END IF
        END IF
    NEXT i%
    IF NOT blocked THEN
        'light was closest, so add its contribution
        'normalize L vector
        L.x = L.x / lightdist
        L.y = L.y / lightdist
        L.z = L.z / lightdist
        'now L is normalized
        clr = clr + phong(L, V, N)
    ELSE
        'actually, we should recursively trace the ray here,
        'but I assume than an eclipsed ray will not make a
        'significant contribution to the brightness
    END IF

NEXT m%
IF clr > 254 THEN clr = 254
'IF clr < 0 THEN clr = 0
totaldirectlight = clr

END FUNCTION