''bluatigro 26 feb 2015
WindowWidth = DisplayWidth
WindowHeight = DisplayHeight
global winx , winy , spmax
winx = WindowWidth
winy = WindowHeight
spmax = 10
dim sp( spmax , 5 )
global spx , spy , spz , spr , spkl , spmat
spx = 0
spy = 1
spz = 2
spr = 3
spkl = 4
spmat = 5
nomainwin
for i = 0 to spmax
  call sphere  i , range( -100 , 100 ) _
                 , range( -100 , 100 ) _
                 , range( 50 , 200 ) _
                 , range( 10 , 50 ) _
                 , range( 0 , 2^24-1 ) , 0
next i

open "ray" for graphics as #m
  #m "trapclose [quit]"
  for x = -100 to 100
    for y = -100 to 100
      kl = pixel( x,y,0 , x,y,1e10 )
      r = kl and 255
      g = int( kl / 256 ) and 255
      b = int( kl / 256 ^ 2 ) and 255
      #m "color " ; rgb$( r , g , b )
      #m "goto " ; winx / 2 + x ; " " ; winy / 2 - y
      #m "down"
      #m "set " ; winy / 2 + x ; " " ; winy / 2 - y
      #m "up"
    next y
  next x
  notice "ready !!"
wait
function range( low , high )
  range = rnd(0) * ( high - low ) + low
end function
sub sphere no , x,y,z , r , kl , mat
  sp( no , spx ) = x
  sp( no , spy ) = y
  sp( no , spz ) = z
  sp( no , spr ) = r
  sp( no , spkl ) = kl
  sp( no , spmat ) = mat
end sub
function rgb$( r , g , b )
  rgb$ = str$( r and 255 ) ; " " ; g and 255 ; " " ; b and 255
end function
function pixel( ox,oy,oz , dx,dy,dz )
  minst = 1e10
  found = -1
  for i = 0 to spmax
    d = hit( i , ox,oy,oz , dx,dy,dz )
    if d > 0 and d < minst then
      minst = d
      found = i
    end if
  next i
  if found < 0 then pixel = 0 : exit function
  pixel = sp( found , spkl )
end function
function hit( no , ox,oy,oz , dx,dy,dz )
  opx = sp( no , spx ) - ox
  opy = sp( no , spy ) - oy
  opz = sp( no , spz ) - oz
  b = dot( opx,opy,opz , dx,dy,dz )
  disc = b ^ 2 * dot( opx,opy,opz , dx,dy,dz ) + sp( no , spr ) ^ 2
  if disc < 0 then hit = -1 : exit function
  disc = sqr( disc )
  t = b^2 - disc
  if t > 1e-10 then hit = t : exit function
  t = b^2 + disc
  if t > 1e-10 then hit = t : exit function
  hit = -1
end function
function dot( ax,ay,az , bx,by,bz )
  dot = ax * bx + ay * by + az * bz
end function
[quit]
  close #m
end

 

nomainwin
'constants to access the point arrays
global X,Y,Z,R,Red,Green,Blue,Xp,Yp,DistToScreen,maxObjects
X=0
Y=1
Z=2
R=3
Red=4
Green=5
Blue=6

'the world is 24x24 units
'0,0 is centre of the world
'-12 far left 12 far right
'12 top -12 bottom
'objects are specified in world units
'rays are specified in world units
'pixel coordinates are converted by
'dividing the world width by the screen width
'and the world height by the screen height
'this gives a per pixel world increment value

ScrXLeft=-12
ScrXRight=12
ScrYTop=-12
ScrYBottom=12
ImgX=512 'pixel width of screen
ImgY=512 'pixel height of screen
ScrX=2*ScrXRight/ImgX 'per pixel world increment value
ScrY=2*ScrYBottom/ImgY

'Calculate how much space windows borders take
'Anatoly's tip
WindowWidth = 200
WindowHeight = 200
open "Ajusting..." for graphics_nsb as #1
#1, "home ; down ; posxy w h"
w=200-2*w : h = 200-2*h
'w and h now contain the number of pixels
'the Windows scheme/theme takes
close #1
WindowWidth  = ImgX+w
WindowHeight = ImgY+h
UpperLeftX = int((DisplayWidth-WindowWidth)/2)
UpperLeftY = int((DisplayHeight-WindowHeight)/2)
open "Raytrace" for graphics_nsb as #1
#1 "down ; fill 192 192 192 ; trapclose [quit]"



'sphere data x,y,z,radius,color
'three spheres ,red green and blue
'evenly spaced around centre of screen
'z order front to back blue, green, red

'there are five massive spheres positioned
'all round the edges and to the rear to create walls


maxObjects=7
dim sp(maxObjects,6)
sp(0,X)=0
sp(0,Y)=0
sp(0,Z)=5
sp(0,R)=4
sp(0,Red)=255
sp(0,Green)=0
sp(0,Blue)=0
sp(1,X)=-6
sp(1,Y)=0
sp(1,Z)=5
sp(1,R)=2
sp(1,Red)=0
sp(1,Green)=255
sp(1,Blue)=0
sp(2,X)=6
sp(2,Y)=0
sp(2,Z)=2
sp(2,R)=2
sp(2,Red)=0
sp(2,Green)=0
sp(2,Blue)=255
sp(3,X)=-36
sp(3,Y)=0
sp(3,Z)=20
sp(3,R)=26
sp(3,Red)=64
sp(3,Green)=128
sp(3,Blue)=255
sp(4,X)=36
sp(4,Y)=0
sp(4,Z)=20
sp(4,R)=26
sp(4,Red)=64
sp(4,Green)=128
sp(4,Blue)=255
sp(5,X)=0
sp(5,Y)=-36
sp(5,Z)=20
sp(5,R)=26
sp(5,Red)=64
sp(5,Green)=128
sp(5,Blue)=25
sp(6,X)=0
sp(6,Y)=36
sp(6,Z)=20
sp(6,R)=26
sp(6,Red)=64
sp(6,Green)=128
sp(6,Blue)=25
sp(7,X)=0
sp(7,Y)=0
sp(7,Z)=36
sp(7,R)=24
sp(7,Red)=128
sp(7,Green)=128
sp(7,Blue)=128

'point origin of camera
'centre screen, 24 world units back from the image plane
cam(X)=0
cam(Y)=0
cam(Z)=-24

'point origin of light source
'centre top, on the image plane
lit(X)=12
lit(Y)=5
lit(Z)=-10

'now cast a ray through every pixel on the screen
for Xp=0 to ImgX
    for Yp=0 to ImgY
    scan

        'point0, the camera x,y,z
        o(X)=cam(X)                     'world x camera 0
        o(Y)=cam(Y)                     'world y camera 0
        o(Z)=cam(Z)                     'world z camera -10


        'point1, the screen x,y,z
        'which is ScrXLeft (-12) + .5 * ScrX, the per pixel world increment value
        'to which we add xp * ScrX to give world x in world units
        r(X)=ScrXLeft+.5*ScrX+Xp*ScrX   'world x
        r(Y)=ScrYTop+.5*ScrY+Yp*ScrY    'world y
        r(Z)=0                          'world z image plane 0

        'flat parrallell ray no perspective
        'o(X)=r(X)
        'o(Y)=r(Y)
        'o(Z)=-10

        'subtract point0 from point1 to get vector direction d()
        d(X)=r(X)-o(X)
        d(Y)=r(Y)-o(Y)
        d(Z)=r(Z)-o(Z)

        'normalize it by dividing by its length to make a unit vector ie it sums to 1
        l=sqr(d(X)*d(X)+d(Y)*d(Y)+d(Z)*d(Z))
        d(X)=d(X)/l
        d(Y)=d(Y)/l
        d(Z)=d(Z)/l
        DistToScreen=l

        'go look for an intersect
        null=raySphereIntersect()
    next 'y
next 'x

wait

[quit]
close #1
end

function raySphereIntersect()
        maxdist=1000000
        for o=0 to maxObjects
            b=2*d(X)*(o(X)-sp(o,X))+2*d(Y)*(o(Y)-sp(o,Y))+2*d(Z)*(o(Z)-sp(o,Z))
            c=(o(X)-sp(o,X))^2+(o(Y)-sp(o,Y))^2+(o(Z)-sp(o,Z))^2-sp(o,R)^2
            d = b * b - 4 * c
            if d>0 then
                t=(b*-1 - sqr(b*b-4*c))/2
                if t>DistToScreen  then
                    'store the shortest intersect of all sphere intersects
                    if t<maxdist then maxdist=t : id=o
                end if
            end if
        next
        if maxdist<1000000 then

            'establish the sphere surface intersect point
            i(X)=o(X)+d(X)*maxdist
            i(Y)=o(Y)+d(Y)*maxdist
            i(Z)=o(Z)+d(Z)*maxdist

            'get unit normal vector from sphere centre to surface intersect
            n(X)=(i(X)-sp(id,X))/sp(id,R)
            n(Y)=(i(Y)-sp(id,Y))/sp(id,R)
            n(Z)=(i(Z)-sp(id,Z))/sp(id,R)


            'get the unit normal vector from sphere surface intersect to the light
            l(X)=lit(X)-i(X)
            l(Y)=lit(Y)-i(Y)
            l(Z)=lit(Z)-i(Z)
            l=sqr(l(X)*l(X)+l(Y)*l(Y)+l(Z)*l(Z))
            l(X)=l(X)/l
            l(Y)=l(Y)/l
            l(Z)=l(Z)/l
            'the dot product of these vectors gives an indication of the light
            color=n(X)*l(X)+n(Y)*l(Y)+n(Z)*l(Z)


            'cast a ray from intersect to the light to check for shadow
            'we have done most of the ray prep
            'point0 is the intersect point1 is the light
            'so l() is our ray direction, point1-point0, normalized
            shadow=1
            for o=0 to maxObjects
                if o<>id then  'check all other spheres not the one we are currently considering
                    b=2*l(X)*(i(X)-sp(o,X))+2*l(Y)*(i(Y)-sp(o,Y))+2*l(Z)*(i(Z)-sp(o,Z))
                    c=(i(X)-sp(o,X))^2+(i(Y)-sp(o,Y))^2+(i(Z)-sp(o,Z))^2-sp(o,R)^2
                    d = b * b - 4 * c
                    if d>0 then t=(b*-1 - sqr(b*b-4*c))/2
                    if d>0 and t>0 then shadow=.5 : exit for
                end if
            next

            'add some ambient light
            if color < .3 then color = .3

            'color the pixel
            #1 "color ";sp(id,Red)*color*shadow;" ";sp(id,Green)*color*shadow;" ";sp(id,Blue)*color*shadow;" ; set ";Xp;" ";Yp

        end if

end function 

nomainwin
'constants to access the point arrays
global X,Y,Z,R,clr,Xp,Yp,DistToScreen,maxObjects
X=0
Y=1
Z=2
R=3
clr=4

global sptel
'the world is 24x24 units
'0,0 is centre of the world
'-12 far left 12 far right
'12 top -12 bottom
'objects are specified in world units
'rays are specified in world units
'pixel coordinates are converted by
'dividing the world width by the screen width
'and the world height by the screen height
'this gives a per pixel world increment value

ScrXLeft=-12
ScrXRight=12
ScrYTop=-12
ScrYBottom=12
ImgX=512 'pixel width of screen
ImgY=512 'pixel height of screen
ScrX=2*ScrXRight/ImgX 'per pixel world increment value
ScrY=2*ScrYBottom/ImgY

'Calculate how much space windows borders take
'Anatoly's tip
WindowWidth = 200
WindowHeight = 200
open "Ajusting..." for graphics_nsb as #1
#1, "home ; down ; posxy w h"
w=200-2*w : h = 200-2*h
'w and h now contain the number of pixels
'the Windows scheme/theme takes
close #1
WindowWidth  = ImgX+w
WindowHeight = ImgY+h
UpperLeftX = int((DisplayWidth-WindowWidth)/2)
UpperLeftY = int((DisplayHeight-WindowHeight)/2)
open "Raytrace" for graphics_nsb as #1
#1 "down ; fill 192 192 192 ; trapclose [quit]"



'sphere data x,y,z,radius,color
'three spheres ,red green and blue
'evenly spaced around centre of screen
'z order front to back blue, green, red

'there are five massive spheres positioned
'all round the edges and to the rear to create walls


maxObjects=7
dim sp(maxObjects,4)
call sphere 0,0,5 , 4 , rgb(255,0,0)
call sphere 0,5,2 , 2 , rgb(0,255,0)
call sphere 6,0,2 , 2 , rgb(0,0,255)
call sphere -1e5-12,20,26 , 1e5 , rgb(64,128,255)
call sphere 1e5+12,20,26 , 1e5 , rgb(64,128,255)
call sphere 0,-1e5-12,26 , 1e5 , rgb(64,128,255)
call sphere 0,1e5+12,26 , 1e5 , rgb(64,128,255)
call sphere 0,0,1e5+12 , 1e5 , rgb(128,128,128)

'point origin of camera
'centre screen, 24 world units back from the image plane
cam(X)=0
cam(Y)=0
cam(Z)=-24

'point origin of light source
'centre top, on the image plane
lit(X)=12
lit(Y)=5
lit(Z)=-10

'now cast a ray through every pixel on the screen
for Xp=0 to ImgX
    for Yp=0 to ImgY
    scan

        'point0, the camera x,y,z
        o(X)=cam(X)                     'world x camera 0
        o(Y)=cam(Y)                     'world y camera 0
        o(Z)=cam(Z)                     'world z camera -10


        'point1, the screen x,y,z
        'which is ScrXLeft (-12) + .5 * ScrX, the per pixel world increment value
        'to which we add xp * ScrX to give world x in world units
        r(X)=ScrXLeft+.5*ScrX+Xp*ScrX   'world x
        r(Y)=ScrYTop+.5*ScrY+Yp*ScrY    'world y
        r(Z)=0                          'world z image plane 0

        'flat parrallell ray no perspective
        'o(X)=r(X)
        'o(Y)=r(Y)
        'o(Z)=-10

        'subtract point0 from point1 to get vector direction d()
        d(X)=r(X)-o(X)
        d(Y)=r(Y)-o(Y)
        d(Z)=r(Z)-o(Z)

        'normalize it by dividing by its length to make a unit vector ie it sums to 1
        l=sqr(d(X)*d(X)+d(Y)*d(Y)+d(Z)*d(Z))
        d(X)=d(X)/l
        d(Y)=d(Y)/l
        d(Z)=d(Z)/l
        DistToScreen=l

        'go look for an intersect
        null=raySphereIntersect()
    next 'y
next 'x

wait

[quit]
close #1
end
function rgb( a , b , c )
  rgb = a + b * 256 + c * 256 ^ 2
end if
function clr.r( kl )
  clr.r = kl and 255
end function
function clr.g( kl )
  clr.g = int( kl / 256 ) and 255
end function
function clr.b( kl )
  clr.b = int( kl / 256 ^ 2 ) and 255
end function
sub sphere a , b , c , d , kl
  if sptel > maxObjects then exit sub
  sp( sptel , X ) = a
  sp( sptel , Y ) = b
  sp( sptel , Z ) = c
  sp( sptel , R ) = d
  sp( sptel , clr ) = kl  
  sptel = sptel + 1
end sub
function raySphereIntersect()
        maxdist=1000000
        for o=0 to maxObjects
            b=2*d(X)*(o(X)-sp(o,X))+2*d(Y)*(o(Y)-sp(o,Y))+2*d(Z)*(o(Z)-sp(o,Z))
            c=(o(X)-sp(o,X))^2+(o(Y)-sp(o,Y))^2+(o(Z)-sp(o,Z))^2-sp(o,R)^2
            d = b * b - 4 * c
            if d>0 then
                t=(b*-1 - sqr(b*b-4*c))/2
                if t>DistToScreen  then
                    'store the shortest intersect of all sphere intersects
                    if t<maxdist then maxdist=t : id=o
                end if
            end if
        next
        if maxdist < 1e7 then

            'establish the sphere surface intersect point
            i(X)=o(X)+d(X)*maxdist
            i(Y)=o(Y)+d(Y)*maxdist
            i(Z)=o(Z)+d(Z)*maxdist

            'get unit normal vector from sphere centre to surface intersect
            n(X)=(i(X)-sp(id,X))/sp(id,R)
            n(Y)=(i(Y)-sp(id,Y))/sp(id,R)
            n(Z)=(i(Z)-sp(id,Z))/sp(id,R)


            'get the unit normal vector from sphere surface intersect to the light
            l(X)=lit(X)-i(X)
            l(Y)=lit(Y)-i(Y)
            l(Z)=lit(Z)-i(Z)
            l=sqr(l(X)*l(X)+l(Y)*l(Y)+l(Z)*l(Z))
            l(X)=l(X)/l
            l(Y)=l(Y)/l
            l(Z)=l(Z)/l
            'the dot product of these vectors gives an indication of the light
            color=n(X)*l(X)+n(Y)*l(Y)+n(Z)*l(Z)


            'cast a ray from intersect to the light to check for shadow
            'we have done most of the ray prep
            'point0 is the intersect point1 is the light
            'so l() is our ray direction, point1-point0, normalized
            shadow=1
            for o=0 to maxObjects
                if o<>id then  'check all other spheres not the one we are currently considering
                    b=2*l(X)*(i(X)-sp(o,X))+2*l(Y)*(i(Y)-sp(o,Y))+2*l(Z)*(i(Z)-sp(o,Z))
                    c=(i(X)-sp(o,X))^2+(i(Y)-sp(o,Y))^2+(i(Z)-sp(o,Z))^2-sp(o,R)^2
                    d = b * b - 4 * c
                    if d>0 then t=(b*-1 - sqr(b*b-4*c))/2
                    if d>0 and t>0 then shadow=.5 : exit for
                end if
            next

            'add some ambient light
            if color < .3 then color = .3

            kl = sp( id , clr )
            'color the pixel
            #1 "color " ; clr.r( kl ) * color * shadow _
                   ;" " ; clr.g( kl ) * color * shadow _
                   ;" " ; clr.b( kl ) * color * shadow 
            #1 "set " ; Xp ; " " ; Yp

        end if

end function 
 

''bluatigro 1 mrt 2015
''rebild of internet raytracer
nomainwin
'constants to access the point arrays
global X,Y,Z,R,clr,Xp,Yp,DistToScreen,maxObjects
X = 0
Y = 1
Z = 2
R = 3
clr = 4

global sptel
'the world is 24x24 units
'0,0 is centre of the world
'-12 far left 12 far right
'12 top -12 bottom
'objects are specified in world units
'rays are specified in world units
'pixel coordinates are converted by
'dividing the world width by the screen width
'and the world height by the screen height
'this gives a per pixel world increment value

global ScrXleft , ScrXRight , ScrYTop , ScrYBottom
ScrXLeft=-12
ScrXRight=12
ScrYTop=-12
ScrYBottom=12
global ImgX , ImgY , ScrX , Scry
ImgX=512 'pixel width of screen
ImgY=512 'pixel height of screen
ScrX=2*ScrXRight/ImgX 'per pixel world increment value
ScrY=2*ScrYBottom/ImgY

'Calculate how much space windows borders take
'Anatoly's tip
WindowWidth = 200
WindowHeight = 200
open "Ajusting..." for graphics_nsb as #1
  #1, "home ; down ; posxy w h"
  w=200-2*w : h = 200-2*h
  'w and h now contain the number of pixels
  'the Windows scheme/theme takes
close #1
WindowWidth  = ImgX+w
WindowHeight = ImgY+h
UpperLeftX = int((DisplayWidth-WindowWidth)/2)
UpperLeftY = int((DisplayHeight-WindowHeight)/2)
open "Raytrace" for graphics_nsb as #1
  #1  "down ; fill black ; trapclose [quit]"


  global red , green , blue , gray 

  red = rgb( 255 , 0 , 0 )
  green = rgb( 0 , 255 , 0 )
  blue = rgb( 0 , 0 , 255 )
  gray = rgb( 200 , 200 , 200 )


  'sphere data x,y,z,radius,color
  'three spheres ,red green and blue
  'evenly spaced around centre of screen
  'z order front to back blue, green, red

  'there are five massive spheres positioned
  'all round the edges and to the rear to create walls


  maxObjects = 10
  dim sp( maxObjects , 4 )
  call sphere 0,0,5 , 4 , red
  call sphere 0,5,2 , 2 , green
  call sphere 6,0,2 , 2 , blue
  call sphere -1e5-12,20,26 , 1e5 , gray
  call sphere 1e5+12,20,26 , 1e5 , gray
  call sphere 0,-1e5-12,26 , 1e5 , gray
  call sphere 0,1e5+12,26 , 1e5 , gray
  call sphere 0,0,1e5+12 , 1e5 , gray

  'point origin of camera
  'centre screen, 24 world units back from the image plane
  cam(X) = 0
  cam(Y) = 0
  cam(Z) = -24

  'point origin of light source
  'centre top, on the image plane
  lit(X) = -11
  lit(Y) = -11
  lit(Z) = -11

  call render
  notice "ready !!"
wait

sub render
'now cast a ray through every pixel on the screen

  for Xp=0 to ImgX
    for Yp=0 to ImgY
      scan

      'point0, the camera x,y,z
      o(X) = cam(X)                     'world x camera 0
      o(Y) = cam(Y)                     'world y camera 0
      o(Z) = cam(Z)                     'world z camera -10


      'point1, the screen x,y,z  
      'which is ScrXLeft (-12) + .5 * ScrX, the per pixel world increment value
      'to which we add xp * ScrX to give world x in world units
      r(X) = ScrXLeft + 0.5 * ScrX + Xp * ScrX   'world x
      r(Y) = ScrYTop + 0.5 * ScrY + Yp * ScrY    'world y
      r(Z) = 0                          'world z image plane 0

      'flat parrallell ray no perspective  
      'o(X)=r(X)
      'o(Y)=r(Y)
      'o(Z)=-10

      'subtract point0 from point1 to get vector direction d()
      d(X) = r(X) - o(X)
      d(Y) = r(Y) - o(Y)
      d(Z) = r(Z) - o(Z)

      'normalize it by dividing by its length to make a unit vector ie it sums to 1
      l = sqr( d(X) * d(X) + d(Y) * d(Y) + d(Z) * d(Z) )
      d(X) = d(X) / l
      d(Y) = d(Y) / l
      d(Z) = d(Z) / l
      DistToScreen = l

      'go look for an intersect
      null = raySphereIntersect()
    next 'y
  next 'x
end sub
[quit]
  close #1
end
function nr$( no , digits )
  nr$ = right$( "0000000000" ; no , digits )
end function
sub savescreen a$
  #1  "getbmp screen 0 0 " ; ImgX ; " " ; ImgY
  bmpsave "screen" , a$ + ".bmp"
end sub
function rgb( a , b , c )
  rgb = a + b * 256 + c * 256 ^ 2
end function
function clr.r( kl )
  clr.r = kl and 255
end function
function clr.g( kl )
  clr.g = int( kl / 256 ) and 255
end function
function clr.b( kl )
  clr.b = int( kl / 256 ^ 2 ) and 255
end function
sub sphere a , b , c , d , kl
  if sptel >= maxObjects then exit sub
  sp( sptel , X ) = a
  sp( sptel , Y ) = b
  sp( sptel , Z ) = c
  sp( sptel , R ) = d
  sp( sptel , clr ) = kl
  sptel = sptel + 1
end sub
function Sphere.dist( o )
  uit = 1e7
  b = 2 * d(X) * ( o(X) - sp(o,X)) _
    + 2 * d(Y) * ( o(Y) - sp(o,Y)) _
    + 2 * d(Z) * ( o(Z) - sp(o,Z))
  c = ( o(X) - sp(o,X) ) ^ 2 _
    + ( o(Y) - sp(o,Y) ) ^ 2 _
    + ( o(Z) - sp(o,Z) ) ^ 2 - sp(o,R) ^ 2
  d = b * b - 4 * c
  if d > 0 then
    uit = ( b * -1 - sqr( b * b - 4 * c ) ) / 2
  end if
  Sphere.dist = uit
end function
function raySphereIntersect()
  maxdist = 1e7
  for o = 0 to sptel
    t = Sphere.dist( o )
    if t > DistToScreen  then
    'store the shortest intersect of all sphere intersects
      if t < maxdist then maxdist = t : id = o
    end if
  next
  if maxdist < 1e7 then
  'establish the sphere surface intersect point
    i(X) = o(X) + d(X) * maxdist
    i(Y) = o(Y) + d(Y) * maxdist
    i(Z) = o(Z) + d(Z) * maxdist
  'get unit normal vector from sphere centre to surface intersect
    n(X) = ( i(X) - sp(id,X) ) / sp(id,R)
    n(Y) = ( i(Y) - sp(id,Y) ) / sp(id,R)
    n(Z) = ( i(Z) - sp(id,Z) ) / sp(id,R)
  'get the unit normal vector from sphere surface intersect to the light
    l(X) = lit(X) - i(X)
    l(Y) = lit(Y) - i(Y)
    l(Z) = lit(Z) - i(Z)
    l = sqr( l(X) * l(X) + l(Y) * l(Y) + l(Z) * l(Z) )
    l(X) = l(X) / l
    l(Y) = l(Y) / l
    l(Z) = l(Z) / l
  'the dot product of these vectors gives an indication of the light
    color = n(X) * l(X) + n(Y) * l(Y) + n(Z) * l(Z)


  'cast a ray from intersect to the light to check for shadow
  'we have done most of the ray prep
  'point0 is the intersect point1 is the light
  'so l() is our ray direction, point1-point0, normalized
    shadow = 1
    o(X) = i(X)
    o(Y) = i(Y)
    o(Z) = i(Z)
    d(X) = l(X)
    d(Y) = l(Y)
    d(Z) = l(Z)
    for o = 0 to sptel
      if o <> id then  
      'check a ll other spheres not the one we are currently considering
        t = Sphere.dist( o )
        if t > 0 and t < 1e7 then shadow = 0.5 : exit function
      end if
    next

    'add some ambient light
    if color < .3 then color = .3

    kl = sp( id , clr )
    'color the pixel
    #1 "color " ; clr.r( kl ) * color * shadow _
          ; " " ; clr.g( kl ) * color * shadow _
          ; " " ; clr.b( kl ) * color * shadow
    #1 "set " ; Xp ; " " ; Yp
  else
    #1 "color black"
    #1 "set " ; Xp ; " " ; Yp
  end if
end function