global v.x,v.y,v.z
global sp.rad,sp.px,sp.py,sp.pz,sp.ex,sp.ey,sp.ez,sp.cx,sp.cy,sp.cz,sp.type
sp.rad=0
sp.px=1
sp.py=2
sp.pz=3
sp.ex=4
sp.ey=5
sp.ez=6
sp.cx=7
sp.cy=8
sp.cz=9
sp.type=10
dim bol(10,10)
global red,green,yellow,blue,mageta,cyan,white
red=rgb(255,0,0)
green=rgb(0,255,0)
yellow=rgb(255,255,0)
blue=rgb(0,0,255)
magenta=rgb(255,0,255)
cyan=rgb(0,255,255)
white=rgb(255,255,255)
lgray=rgb(200,200,200)
white2=rgb(254,254,254)
global diff , spec , refr
diff=0
spec=1
refr=2
WindowWidth = 400
WindowHeight = 400
global winx,winy
winx=WindowWidth
winy=WindowHeight
''scene
call Sphere 0,1e5, 1e5+1,40.8,81.6,   0,red,diff ''Left 
call Sphere 1,1e5, -1e5+99,40.8,81.6, 0,cyan,diff ''Rght 
call Sphere 2,1e5, 50,40.8, 1e5,      0,magenta,diff ''Back 
call Sphere 3,1e5, 50,40.8,-1e5+170,  0,0,diff ''Frnt 
call Sphere 4,1e5, 50, 1e5, 81.6,     0,green,diff ''Botm 
call Sphere 5,1e5, 50,-1e5+81.6,81.6, 0,blue,diff ''Top 
call Sphere 6,16.5, 27,16.5,47,       0,white2, spec ''Mirr 
call Sphere 7,16.5, 73,16.5,78,       0,white2, refr ''Glas 
call Sphere 8,600, 50,681.6-.27,81.6 ,white,0, diff ''Lite 


end

''sphere

sub Sphere no,r,x,y,z,e,kl,type
  bol(no,sp.px)=x
  bol(no,sp.py)=y 
  bol(no,sp.pz)=z
  bol(no,sp.rad)=r 
  call torgb e
  bol(no,sp.ex)=v.x  
  bol(no,sp.ey)=v.y  
  bol(no,sp.ez)=v.z
  call torgb kl
  bol(no,sp.cx)=v.x
  bol(no,sp.cy)=v.y
  bol(no,sp.cz)=v.z
  bol(no,sp.type)=type
end sub
function hit( no , ox,oy,oz,dx,dy,dz )
  x=bol(no,sp.px)
  y=bol(no,sp.py)
  z=bol(no,sp.pz)
  rad=bol(mo,sp.rad)
  call minus x,y,z,ox,oy,oz
  op.x=v.x:op.y=v.y:op.z=v.z
  eps = 1e-10
  b = dot(op.x,op.y,op.z,ox,oy,oz)
  det = b*b-dot(op.x,op.y,op.z,op.x,op.y,op.z)+rad*rad
  if det < 0 then
    hit = 0
    exit function
  end if
  det = sqr( det )
  t = b - det
  if t > eps then
    hit = t
    exit function
  else
    t = b + det
    if t > eps then
      hit = t
      exit function
    else
      hit = 0
    end if
  end if
end function

end function

''color math

sub torgb kl
  r = int( kl and 255 ) / 256
  g = int( kl/256 and 255 ) / 256
  b = int( kl/256^2 and 255 ) / 256
  call v3d r , g , b
end sub
function rgb( r , g , b )
  rgb = r+g*256+b*256^2
end function

''vertor math

sub v3d x , y , z
''constructor
  v.x = x
  v.y = y
  v.z = z
end sub
sub add a,b,c,d,e,f
  call v3d a+d,b+e,c+f
end sub
sub minus a,b,c,d,e,f
  call v3d a-d,b-e,c-f
end sub
sub mul a,b,c,d
  call v3d a*d,b*d,c*d
end sub
sub mult a,b,c,d,e,f
  call v3d a*d,b*e,c*f
end sub
function dot(a,b,c,d,e,f)
  dot=a*d+b*e+c*f
end function
sub cross ax,ay,az,bx,by,bz
  call v3d ay*bz-az*by,az*bx-ax*bz,ax*by-ay*bx
end sub
sub norm a,b,c
  l=sqr(a^2+b^2+c^2)
  call mul a,b,c,1/l
end sub
 

global trimax , size
trimax = 10
size = 200
dim tri(trimax,13) , pnt(255,3)
WindowWidth = DisplayWidth
WindowHeight = DisplayHeight
global winx , winy
winx = WindowWidth
winy = WindowHeight

for i = 0 to trimax
  for t = 0 to 2
    call point t,range(0-size,size) _
                ,range(0-size,size) _
                ,range(0,size)
  next t
  call tri i , 0,1,2 , int(rnd(0)*(2^24-1)) , 0
next i

open "ray triangles" for graphics as #m
  #m "trapclose [quit]"
  for x = 0 - size to size
    for y = 0 - size to size
      scan
      kl = ray( 0,0,0-size , x/size,y/size,1)
      #m "goto ";x+winx/2;" ";winy/2-y
      #m "down"
      r = kl and 255
      g = int( kl / 256 ) and 255
      b = int( kl / 256 ^ 2 ) and 255
      #m "color ";r;" ";g;" ";b
      #m "set ";x+winx/2;" ";winy/2-y
      #m "up"
     next y
  next x
  notice "ready !!"
wait
[quit]
  close #m
end
function ray( ox,oy,oz , dx,dy,dz )
  id = -1
  distmin = 1e9
  for i = 0 to trimax
    dist = trihit( i , ox,oy,oz , dx,dy,dz )
    if dist < distmin and dist > 0 then
      distmin = dist
      id = i
    end if
  next i
  if id >= 0 then
    ray = tri( id , 12 )
  end if
  ray = 0
end function
function rgb( r , g , b )
  r = r and 255
  g = g and 255
  b = b and 255
  rgb = r + g * 256 + b * 256 ^ 2
end function
function range( l , h )
  range = rnd(0) * ( h - l ) + l
end function
function trihit( no , px,py,pz , dx,dy,dz )
  v0x = tri(no,0)
  v0y = tri(no,1)
  v0z = tri(no,2)
  v1x = tri(no,3)
  v1y = tri(no,4)
  v1z = tri(no,5)
  v2x = tri(no,6)
  v2y = tri(no,7)
  v2z = tri(no,8)
  call vsub e1x,e1y,e1z , v1x,v1y,v1z , v0x,v0y,v0z
  call vsub e2x,e2y,e2z , v2x,v2y,v2z , v0x,v0y,v0z
  
  call vcross hx,hy,hz , dx,dy,dz , e2x,e2y,e2z
  a = dot( e1x,e1y,e1z , hx,hy,hz )
  if abs( a ) < 1e-9 then
    trihit = -1
  end if

  f = 1 / a
  call vsub sx,sy,sz , px,py,pz , v0x,v0y,v0z
  u = f * dot( sx,sy,sz , hx,hy,hz )
  if u < 0 or u > 1 then
    trihit = -1
  end if
  
  call vcross qx,qy,qz , sx,sy,sz , e1x,e1y,e1z
  v = f * dot( dx,dy,dz , qx,qy,qz )
  if v < 0 or v > 1 then
    trihit = -1
  end if

  t = f * dot( e2x,e2y,e2z , qx,qy,qz )

  if t > 1e-9 then
    trihit = t
  else
    trihit = -1
  end if

end function
sub point no , x,y,z
  pnt(no,0)=x
  pnt(no,1)=y
  pnt(no,2)=z
end sub
sub tri no , p1 , p2 , p3 , kl , t
  tri(no,0)=pnt(p1,0)
  tri(no,1)=pnt(p1,1)
  tri(no,2)=pnt(p1,2)
  tri(no,3)=pnt(p2,0)
  tri(no,4)=pnt(p2,1)
  tri(no,5)=pnt(p2,2)
  tri(no,6)=pnt(p3,0)
  tri(no,7)=pnt(p3,1)
  tri(no,8)=pnt(p3,2)
  x1=pnt(p2,0)-pnt(p1,0)
  x2=pnt(p3,0)-pnt(p1,0)
  y1=pnt(p2,1)-pnt(p1,1)
  y2=pnt(p3,1)-pnt(p1,1)
  z1=pnt(p2,2)-pnt(p1,2)
  z2=pnt(p3,2)-pnt(p1,2)  
  call vcross x,y,z , x1,y1,z1 , x2,y2,z2
  tri(no,9)=x
  tri(no,10)=y
  tri(no,11)=z
  tri(no,12)=kl
  tri(no,13)=t
end sub
function dot( ax,ay,az , bx,by,bz )
  dot = ax * bx + ay * by + az * bz
end function
function vlen( x,y,z )
  vlen = sqr( dot( x,y,z , x,y,z ) )
end function
sub reflect byref x,byref y,byref z , ax,ay,az , bx,by,bz
  d = dot( ax,ay,az , bx,by,bz )
  call vadd cx,cy,cz , bx,by,bz , bx,by,bz
  call vscale cx,cy,cz , d
  call vsub x,y,z , cx,cy,cz , ax,ay,az
end sub
sub v3normal byref x,byref y,byref z _
, ax,ay,az , bx,by,bz , cx,cy,cz
  call vsub dx,dy,dz , bx,by,bz , ax,ay,az
  call vsub ex,ey,ez , cx,cy,cz , ax,ay,az
  call enormal x,y,z , ex,ey,ez , dx,dy,dz
end sub
sub enormal byref x,byref y,byref z , ax,ay.az , bx,by,bz
  call vcross x,y,z , ax,ay,az , bx,by,bz
  call vnormal x,y,z , x,y,z
end sub
sub vnormal byref x,byref y,byref z , a,b,c
  l = dot( a,b,c , a,b,c )
  if l > 1e-6 then
    l2 = sqr( l )
    call vmul x,y,z , a,b,c , 1/l2
  end if
end sub
sub vadd byref x,byref y,byref z , ax,ay,az , bx,by,bz
  x = ax + bx
  y = ay + by
  z = az + bz
end sub
sub vsub byref x,byref y,byref z , ax,ay,az , bx,by,bz
  x = ax - bx
  y = ay - by
  z = az - bz
end sub
sub vscale byref x,byref y,byref z , d
  x = x * d
  y = y * d
  z = z * d
end sub
sub vmul byref x , byref y , byref z , a,b,c , d
  x = a * d
  y = b * d
  z = c * d
end sub
sub vcross byref x,byref y,byref z , ax,ay,az , bx,by,bz
  x = ay * bz - az * by
  y = az * bx - ax * bz
  z = ax * by - ay * bx
end sub