WindowWidth = DisplayWidth
WindowHeight = DisplayHeight
global winx , winy
gosub [math]
gosub [color]
gosub [ray]
p.y = -400
p.ny = 1
p.kl = green
stel = 0
call sphere -40,40,0 , 40 , red
call sphere 40,40,0 , 30 , yellow
ttel = 0
call tri 0,0,0 , 50,50,0 , 70,-50,0 , blue
call tri 0,0,0 , -30,30,0, -50,-30,0 , cyan
winx = WindowWidth
winy = WindowHeight
nomainwin
open "ray 2.0" for graphics as #m
  #m "trapclose [quit]"
  for x = -100 to 100
    for y = -100 to 100
      #m "goto ";winx/2+x;" ";winy/2+y
      call setcolor render(0,0,-1000,x,y,1000,0)
      #m "down"
      #m "set ";winx/2+x;" ";winy/2+y
      #m "up"
    next y
  next x
  notice "ready"
wait
[quit]
  close #m
end

''bluatigro 18 nov 2016
''math block

[math]
  global pi , golden.ratio
  pi = atn( 1 ) * 4
  golden.ratio = ( sqr( 5 ) - 1 ) / 2
return

function rad( deg )
''from degrees to radians
  rad = deg * pi / 180
end function

sub rotate byref k , byref l , deg
  s = sin( rad( deg ) )
  c = cos( rad( deg ) )
  hk = k * c - l * s
  hl = k * s + l * c
  k = hk
  l = hl
end sub

function nr$( no , m )
  nr$ = right$( "00000000" ; no , m )
end function

sub saveframe f$ , no , m
  #m "getbmp screen 0 0 " ; winx ; " " ; winy
  bmpsave "screen" , f$ + nr$( no , m ) + ".bmp"
end sub

''bluatigro 18 nov 2016
''color block
''needs math block

[color]
  global black , red , green , yellow
  global blue , magenta , cyan , white
  global pink , orange , gray , purple
  black   = rgb(   0 ,   0 ,   0 )
  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 )
  pink    = rgb( 255 , 127 , 127 )
  orange  = rgb( 255 , 127 ,   0 )
  gray    = rgb( 127 , 127 , 127 )
  purple  = rgb( 127 ,   0 , 127 )
return

sub setcolor kl
  r = int( kl and 255 )
  g = int( kl / 256 ) and 255
  b = int( kl / 256 / 256 ) and 255
  print #m , "backcolor " ; r ;" "; g ; " "; b
  print #m , "color " ; r ; " " ; g ; " " ; b
end sub

function rgb( r , g , b )
  rgb = ( int( r ) and 255 ) _
  + ( int( g ) and 255 ) * 256 _
  + ( int( b ) and 255 ) * 256 * 256
end function

function rainbow( x )
  r = sin( rad( x ) ) * 127 + 128
  g = sin( rad( x - 120 ) ) * 127 + 128
  b = sin( rad( x + 120 ) ) * 127 + 128
  rainbow = rgb( r , g , b )
end function

function mix( kl1 , f , kl2 )
  r1 = int( kl1 and 255 )
  g1 = int( kl1 / 256 ) and 255
  b1 = int( kl1 / 256 / 256 ) and 255
  r2 = int( kl2 and 255 )
  g2 = int( kl2 / 256 ) and 255
  b2 = int( kl2 / 256 / 256 ) and 255
  r = r1 + ( r2 - r1 ) * f
  g = g1 + ( g2 - g1 ) * f
  b = b1 + ( b2 - b1 ) * f
  mix = rgb( r , g , b )
end function

''bluatigro 18 nov 2016
''ray casting block
''needs vector3D , color and math block

[ray]
  global great , smal
  great = 1e7
  smal = 1e-7
  global p.x,p.y,p.z,p.nx,p.ny,p.nz,p.kl
  global smax , stel
  smax = 10
  dim s.x(smax),s.y(smax),s.z(smax),s.r(smax),s.r2(smax),s.kl(smax)
  global tmax , ttel
  tmax = 10
  dim t.x0(tmax),t.y0(tmax),t.z0(tmax) _
        ,t.x1(tmax),t.y1(tmax),t.z1(tmax) _
        ,t.x2(tmax),t.y2(tmax),t.z2(tmax) _
        ,t.xn(tmax),t.yn(tmax),t.zn(tmax),t.kl(tmax)
return

function render(ox,oy,oz,dx,dy,dz,dept)
  scan
  p = plane.hit(0,ox,oy,oz,dx,dy,dz)
  s = great
  for i = 0 to stel
    sdist = sphere.hit(i,ox,oy,oz,dx,dy,dz)
    if sdist < s then
      s = sdist
      si = i
    end if
  next i
  t = great
  for i = 0 to ttel
    tdist = tri.hit(i,ox,oy,oz,dx,dy,dz)
    if tdist < t then
      t = tdist
      ti = i
    end if
  next i
  if t < s and t < p then
    uit = t.kl(ti)
  else
    if s < p then
      call v.normalize x,y,z , dx,dy,dz
      call v.add x,y,z , ox,oy,oz , x*s,y*s,z*s
      call v.min x,y,z , x,y,z , s.x(si),s.y(si),s.z(si)
      a = angle( x,y,z , 0,1,0 )
      kl = mix( s.kl(si) , cos(a)/2+0.5 , black )
      uit = kl
    else
      if p < great then
        uit = p.kl
      else
        uit = black
      end if
    end if
  end if
  render = uit
end function

function plane.hit(no,ox,oy,oz,dx,dy,dz)
  demon = dot(p.nx,p.ny,p.nz , dx,dy,dz)
  if demon > smal then
    call v.min x,y,z , p.x,p.y,p.z , ox,oy,oz
    uit = dot(x,y,z , p.nx,p.ny,p.nz)
  else
    uit = great
  end if
  plane.hit = uit
end function

sub sphere x,y,z , d , kl
  if stel > smax then exit sub
  s.x(stel)=x
  s.y(stel)=0-y
  s.z(stel)=z
  s.r(stel)=d
  s.r2(stel)=d*d
  s.kl(stel)=kl
  stel=stel+1
end sub

function sphere.hit(no,ox,oy,oz,dx,dy,dz)
  call v.min x,y,z , ox,oy,oz , s.x(no),s.y(no),s.z(no)
  a = dot(dx,dy,dz , dx,dy,dz)
  b = 2 * dot(x,y,z , dx,dy,dz)
  c = dot(x,y,z , x,y,z)-s.r2(no)
  disc = b*b-4*a*c
  if disc < 0 then
    uit = great
  else
    e = sqr(disc)
    demon = 2*a
    t = ( 0-b - e ) / demon
    if t > smal then
      uit = t
    end if
    t = ( 0-b + e ) / demon
    if t > smal then
      uit = t
    end if
  end if
  sphere.hit = uit
end function

sub tri x0,y0,z0,x1,y1,z1,x2,y2,z2,kl
  t.x0(ttel)=x0
  t.y0(ttel)=y0
  t.z0(ttel)=z0
  t.x1(ttel)=x1
  t.y1(ttel)=y1
  t.z1(ttel)=z1
  t.x2(ttel)=x2
  t.y2(ttel)=y2
  t.z2(ttel)=z2
  t.kl(ttel)=kl
  call v.min ax,ay,az , x2,y2,z2 , x0,y0,z0
  call v.min bx,by,bz , x1,y1,z1 , x0,y0,z0
  call v.cross xn,yn,zn , ax,ay,az , bx,by,bz
  t.xn(ttel)=xn
  t.yn(ttel)=yn
  t.zn(ttel)=zn
  ttel=ttel+1
end sub

function tri.hit(no,ox,oy,oz,dx,dy,dz)
  call v.min e1x,e1y,e1z , t.x1(no),t.y1(no),t.z1(no) , t.x0(no),t.y0(no),t.z0(no)
  call v.min e2x,e2y,e2z , t.x2(no),t.y2(no),t.z2(no) , t.x0(no),t.y0(no),t.z0(no)
  call v.cross qx,qy,qz , dx,dy,dz , e2x,e2y,e2z
  a = dot(e1x,e1y,e1z , qx,qy,qz)
  if abs( a ) < smal then
    tri.hit = great
    exit function
  end if
  f = 1 / a
  call v.min x,y,z , ox,oy,oz , t.x0,t.y0,t.z0
  u = f*dot(x,y,z , qx,qy,qz)
  if u < 0 or u > 1 then
    tri.hit = great
    exit function
  end if
  call v.cross x,y,z , x,y,z , e1x,e1y,e1z
  v = f*dot(dx,dy,dz , x,y,z)
  if v < 0 or u + v > 1 then
    tri.hit = great
    exit function
  end if
  tri.hit = f*dot(e2x,e2y,e2z , x,y,z)
end function

''bluatigro 18 nov 2016
''vector3D block

sub v.min byref x,byref y,byref z,a,b,c,d,e,f
  x = a - d
  y = b - e
  z = c - f
end sub

sub v.add byref x,byref y,byref z,a,b,c,d,e,f
  x = a + d
  y = b + e
  z = c + f
end sub

sub v.cross byref x,byref y,byrefz,x1,y1,z1,x2,y2,z2
  x = y2*z1-z2*y1
  y = z2*x1-x2*z1
  z = x2*y1-y2*x1
end sub

sub v.normalize byref x,byref y,byref z,a,b,c
  l = lenght(a,b,c)
  x = a/l
  y = b/l
  z = c/l
end sub

function dot(a,b,c,d,e,f)
  dot = a * d + b * e + c * f
end function

function angle(a,b,c,d,e,f)
  dt = dot(a,b,c,d,e,f)
  la = lenght(a,b,c)
  lb = lenght(d,e,f)
  angle = acs( dt / ( la * lb ) )
end function

function lenght(a,b,c)
  lenght = sqr(a*a+b*b+c*c)
end function