vb 2017 : raytracing ?

'' bluatigro 30 aug 2017
'' frm_ray
'' needs m4x4 , dbl3d , basis3d and ray_world

public class frm_ray

    public red as new dbl3d(1,0,0)
    public green as new dbl3d(0,1,0)

    private sub set_pixel(x as int16,y as int16,kl as color)
    ''i dont have the book jet
    end sub

    private sub frm_ray_load()
        '' at begin of scene :
        spheretel = 0
        '' set drawing matix
        link(1 , 0,0,0 , 0,0,0 , xyz , 0)
        '' floting sphere
        add_sphere(0,0,0,100,red)
        '' ground
        add_sphere(0,-me.height/2-verry_big _
        ,0,verry_big,green,0.5)
    end sub

    private sub frm_paint()
        dim x as double
        dim y as double
        dim o as dbl3d
        dim d as dbl3d
        ''for every pixel on form
        for i = -me.width/2 to me.width/2
            for j = -me.height/2 to me.height/2
                 '' camara origon
                 o = dbl3d(0,0,-1000)
                 '' camara direction
                 d = dbl3d(i,j,1000)
                 set_pixel(cint(i),cint(j) _
                 ,render(o,d,7).tocolor)
            next j
        next x
    end sub

end class

'' bluatigro 29 aug 2017
'' ray world module
'' needs m4x4 , dbl3d and basis3d

public module ray_world

    public const verry_big as double = 1e9
    public const verry_small as double = 1e-9

    public spheretel as int16
    public const spheremax as int16 = 200
    public spheres(spheremax) as ray_sphere

    public light as new dbl3d(-1,1,-1)

    public class ray_mat
    '' ray material
    '' wil be extended

        public diffuse as dbl3d
        public reflection as double

    end class

    public class ray_sphere
    '' ray sphere shape
        public center as dbl3d
        public r as double
        public r2 as double
        public mat as ray_mat
        
        public sub new()
        end sub

        public sub fill(cx as double _
        ,cy as double,cz as double,r as double _
        ,m as ray_mat)
            spot(cx,cy,cz)
            center = dbl3d(cx,cy,cz)
            me.r = r
            r2 = r * r
            mat = m
        end sub

        public function hit( o as t3d , d as t3d ) as double
        '' calculate distance to sphere
        '' or return verry_big if sphere wasn't hit
            dim as double t , a , b , c , disc
            dim as t3d temp = o - center
            a = d.dot( d )
            b = 2 * temp.dot( d )
            c = temp.dot( temp ) - r2
            disc = b * b - 4 * a * c
            if disc < 0 then 
                return verry_big
            else
                dim as double e = math.sqrt( disc ) 
                dim as double demon = 2 * a
                t = ( -b - e ) / demon
                if t > verry_small then
                    return t
                end if
                t = ( -b + e ) / demon
                if t > verry_small then 
                    return t
                end if
            end if
            return verry_big
        end function

    end class 

    public sub add_sphere(x as double _
    ,y as double,z as double,r as double,kl as dbl3d _
    ,f as double = 0)
    '' add a sphere to the world 
        if spheretel >= spheremax then exit sub
        dim mat as ray_mat
        mat.diffuse = kl
        sphere(spheretel).fill(x,y,z,r,mat,f)
        spheretel += 1
    end sub

    public function render(o as dbl3d _
    ,d as dbl3d,depth as int16)as dbl3d
    '' get color of pixel
        dim i as int16
        dim isph as int16 = -1
        dim b as double
        dim sphdist as double = verry_big
        '' get clostest sphere distance
        for i = o to spheretel
            b = spheres(i).hit(o,d)
            if b < sphdist then
                sphdist = dist
                isph = i
            end if
        next i
        '' get color of clostest sphere
        '' and ca***ate shaded color
        '' wil be extended
        if i >= 0 then
            dim kl as dbl3d = spheres(isph).mat.diffuse
            dim p as dbl3d = o + d.normalize() * b
            dim n as dbl3d = p - spheres(isph).center
            dim a as double = getangle(n,light)
            return kl*(math.cos(a)/2+.5)
        end if
        '' if no sphere was hit return black
        return dbl3d()
    end function

end module

'' bluatigro 29 aug 2017
'' 3d engine module
'' needs m4x4 and Dbl3D

public module basis3D

    private dim sk(64) as Dbl3D
    private dim v(20) as m4x4
    private number as int16

    public const xyz as int16 = 0
    public const xzy as int16 = 1
    public const yxz as int16 = 2
    public const yzx as int16 = 3
    public const zxy as int16 = 4
    public const zyx as int16 = 5

    public sub link(no as int16 _
    ,x as double,y as double,z as double _
    ,xz as double,yz as double,xy as double _
    ,ax as int16,p as int16)
    '' create new matrix
    '' wil effect future drawings 
        if no < 1 or n > 20 then exit sub
        if p < 0 or p > 20 then exit sub
        if p = no then exit sub
        dim trans as m4x4
        dim rotx as m4x4
        dim roty as m4x4
        dim rotz as m4x4
        trans.translate(x,y,z)
        rotx.rotate(yz)
        roty.rotate(xz,1)
        rotz.rotate(xy,2)
        select case ax
            case xyz
                v(no)=rotx*roty*rotz*trans*v(p)
            case xzy
                v(no)=rotx*rotz*roty*trans*v(p)
            case yxz
                v(no)=roty*rotx*rotz*trans*v(p)
            case yzx
                v(no)=roty*rotz*rotx*trans*v(p)
            case zxy
                v(no)=rotz*rotx*roty*trans*v(p)
            case zyx
                v(no)=rotz*roty*rotx*trans*v(p)
            case else
                v(no)=rotx*roty*rotz*trans*v(p)
        end select
        number = no
    end sub

    public sub child(no as int16 _
    ,x as double,y as double,z as double _
    ,lim as int16,ax as int16,p as int16)
    '' create a lim matrix of a avatar
    '' wil effect furure drawings
        link(no,x,y,z,sk(lim).y,sk(lim).x,sk(lim).z),ax,p)
    end sub

    public sub skelet(lim as int16 _
    ,x as double,y as double,z as double)
    '' set the skelet lim of a avatar
        sk(lim)=Dbl3D(x,y,z)
    end sub

    public function rad( deg as double ) as double
        return deg * math.pi / 180.0
    end function

    public function pend( f as double _
    , amp as double ) as double
    '' for smoot animations of lims
        return math.sin(rad(f))*amp
    end function

    public sub spot(byref x as double _
    ,byref y as double,byref z as double)
    '' from local to world coordinates
        dim v3d as dbl3d(x,y,z)
        v3d=v(number)*v3d
        x=v3d.x
        y=v3d.y
        z=v3d.z
    end sub

    public sub local(byref x as double _
    ,byref y as double,byref z as double)
    '' from world to local coordinates
    '' only rotation matrix
        dim v3d as dbl3d(x,y,z)
        dim m44 as m4x4 
        m44 = m44.inverse(v(number))
        m44.translate(0,0,0)
        v3d=m44*v3d
        x=v3d.x
        y=v3d.y
        z=v3d.z
    end sub

end module

'' bluatigro 29 aug 2017
'' matrix 4 x 4 class
'' needs Dbl3D

Public Class m4x4

    Public m(3, 3) As Double
    Private Const deg2rad as Double = Math.pi / 180

    Public Sub New()
        Dim i as Integer
        Dim j as Integer
        For i = 0 to 3
            For j = 0 to 3
                m(i, j) = 0.0
            Next j
            m(i, i) = 1.0
        Next i
    End Sub

    Public Shared Operator *(a As m4x4, b As m4x4) As m4x4
        Dim uit As m4x4
        Dim i as Int16
        Dim j as Int16
        Dim k as Int16
        For i = 0 to 3
            For j = 0 to 3
                uit.m(i, j) = 0.0
                For k = 0 to 3
                    uit.m(i, j) += a.m(i, k) * b.m(k, j)
                Next k
            Next j
        Next i
        Return uit
    End Operator

    public shared operator *(m as m4x4, v as dbl3d) as dbl3d
        return new dbl3d(v.x*m.m(0,0)+v.y*m.m(1,0) _
                        +v.z*m.m(2,0)+m.m(3,0) _
                        ,v.x*m.m(0,1)+v.y*m.m(1,1)  _
                        +v.z*m.m(2,1)+m(3,1) _
                        ,v.x*m.m(0,2)+v.y*m.m(1,2) _
                        +v.z*m.m(2,2)+m.m(3,2))
    end operator

    public function inverse(q as m4x4) as m4x4
    '' for the inverse rotation matrix
    '' it wil NOT calculate inverse translation
        dim d as double 
        dim uit as m4x4 
  
        d = q.m(0,0)*q.m(1,1)*q.m(2,2) _
          - q.m(0,0)*q.m(2,1)*q.m(1,2) _
          + q.m(1,0)*q.m(2,1)*q.m(0,2) _
          - q.m(1,0)*q.m(0,1)*q.m(2,2) _
          + q.m(2,0)*q.m(0,1)*q.m(1,2) _
          - q.m(2,0)*q.m(1,1)*q.m(0,2)
    
        uit.m(0,0)=(q.m(1,1)*q.m(2,2)-q.m(1,2)*q.m(2,1))/d
        uit.m(1,0)=(q.m(0,1)*q.m(2,2)-q.m(0,2)*q.m(2,1))/d
        uit.m(2,0)=(q.m(0,1)*q.m(1,2)-q.m(0,2)*q.m(1,1))/d
        uit.m(0,1)=(q.m(1,0)*q.m(2,2)-q.m(1,2)*q.m(2,0))/d
        uit.m(1,1)=(q.m(0,0)*q.m(2,2)-q.m(0,2)*q.m(2,0))/d
        uit.m(2,1)=(q.m(0,0)*q.m(1,2)-q.m(0,2)*q.m(1,0))/d
        uit.m(0,2)=(q.m(1,0)*q.m(2,1)-q.m(1,1)*q.m(2,0))/d
        uit.m(1,2)=(q.m(0,0)*q.m(2,1)-q.m(0,1)*q.m(2,0))/d
        uit.m(2,2)=(q.m(0,0)*q.m(1,1)-q.m(0,1)*q.m(1,2))/d
  
        return uit
  
    end function

    Public Sub translate(x As Double _
    ,y As Double ,z As Double)
        m(3, 0) = x
        m(3, 1) = y
        m(3, 2) = z
    End Sub

    Public Sub rotate(deg As Double _
    , ax As Int16 = 0)
        Dim a As Int16 = 1
        Dim b As Int16 = 2
        If ax = 1 Then
            a = 0
            b = 2
        End If
        If ax = 2 Then
            a = 0
            b = 1
        End If
        m(a, a) = Math.cos(deg * deg2rad)
        m(a, b) = -Math.sin(deg * deg2rad)
        m(b, a) = Math.sin(deg * deg2rad)
        m(b, b) = Math.cos(deg * deg2rad)
    End Sub

    Public Sub scale(x As Double _
    ,y As Double ,z As double)
        m(0, 0) = x
        m(1, 1) = y
        m(2, 2) = z
    End Sub

'' bluatigro 29 aug 2017
'' double 3d class

Public Class Dbl3D

    Public x As Double
    Public y As Double
    Public z As Double

    Public Sub New()
        x = 0
        y = 0
        z = 0
    End Sub

    Public Sub New(a As Double, b As Double, c As Double)
        x = a
        y = b
        z = c
    End Sub

    Public Function toColor() As Color
        Return Color.FromArgb(255 _
        , CInt(x * 255), CInt(y * 255), CInt(z * 255))
    End Function

    Public Shared Operator +(a As Dbl3D, b As Dbl3D) As Dbl3D
        Return New Dbl3D(a.x + b.x, a.y + b.y, a.z + b.z)
    End Operator

    Public Shared Operator -(a As Dbl3D, b As Dbl3D) As Dbl3D
        Return New Dbl3D(a.x - b.x, a.y - b.y, a.z - b.z)
    End Operator

    Public Shared Operator /(a As Dbl3D, b As Double) As Dbl3D
        Return New Dbl3D(a.x / b, a.y / b, a.z / b)
    End Operator

    Public Shared Operator *(a As Dbl3D, b As Double) As Dbl3D
        Return New Dbl3D(a.x * b, a.y * b, a.z * b)
    End Operator

    Public Function dot(b As Dbl3D) As Double
        Return x * b.x + y * b.y + z * b.z
    End Function

    Public Function cross(b As Dbl3D) As Dbl3D
        Return New Dbl3D(y * b.z - z * b.y _
                        ,z * b.x - x * b.z _
                        ,x * b.y - y * b.x)
    End Function

    Public Function length() As Double
        Return Math.sqrt(dot(Me))
    End Function

    Public Function getAngle(b As Dbl3D) As Double
        Dim la as Double length()
        Dim lb as Double b.length()
        Dim d as Double dot(b)
        Return Math.acs(d / (la * lb))
    End Function 

    Public Sub normalize() 
        Me = Me / length()
    End Sub

End Class

'' bluatigro 29 aug 2017
'' ray world module
'' needs m4x4 , dbl3d and basis3d

public module ray_world

    public const verry_big as double = 1e9
    public const verry_small as double = 1e-9

    public spheretel as int16
    public const spheremax as int16 = 200
    public spheres(spheremax) as ray_sphere

    public light as new dbl3d(-1,1,-1)

    public class ray_mat
    '' ray material
    '' wil be extended

        public diffuse as dbl3d
        public reflection as double

    end class

    public class ray_sphere
    '' ray sphere shape
        public center as dbl3d
        public r as double
        public r2 as double
        public mat as ray_mat
        
        public sub new()
        end sub

        public sub fill(cx as double _
        ,cy as double,cz as double,r as double _
        ,m as ray_mat)
            spot(cx,cy,cz)
            center = dbl3d(cx,cy,cz)
            me.r = r
            r2 = r * r
            mat = m
        end sub

        public function hit( o as t3d , d as t3d ) as double
        '' calculate distance to sphere
        '' or return verry_big if sphere wasn't hit
            dim as double t , a , b , c , disc
            dim as t3d temp = o - center
            a = d.dot( d )
            b = 2 * temp.dot( d )
            c = temp.dot( temp ) - r2
            disc = b * b - 4 * a * c
            if disc < 0 then 
                return verry_big
            else
                dim as double e = math.sqrt( disc ) 
                dim as double demon = 2 * a
                t = ( -b - e ) / demon
                if t > verry_small then
                    return t
                end if
                t = ( -b + e ) / demon
                if t > verry_small then 
                    return t
                end if
            end if
            return verry_big
        end function

    end class 

    public sub add_sphere(x as double _
    ,y as double,z as double,r as double,kl as dbl3d _
    ,f as double = 0)
    '' add a sphere to the world 
        if spheretel >= spheremax then exit sub
        dim mat as ray_mat
        mat.diffuse = kl
        sphere(spheretel).fill(x,y,z,r,mat,f)
        spheretel += 1
    end sub

    public function render(o as dbl3d _
    ,d as dbl3d,depth as int16)as dbl3d
    '' get color of pixel
        dim i as int16
        dim isph as int16 = -1
        dim b as double
        dim sphdist as double = verry_big
        '' get clostest sphere distance
        for i = o to spheretel
            b = spheres(i).hit(o,d)
            if b < sphdist then
                sphdist = dist
                isph = i
            end if
        next i
        '' get color of clostest sphere
        '' and ca***ate shaded color
        '' wil be extended
        if i >= 0 then
            dim kl as dbl3d = spheres(isph).mat.diffuse
            dim p as dbl3d = o + d.normalize() * b
            dim n as dbl3d = p - spheres(isph).center
            dim a as double = n.getangle(light)
            return kl*(math.cos(a)/2+.5)
        end if
        '' if no sphere was hit return black
        return dbl3d()
    end function

end module

'' bluatigro 1 sept 2017
'' all basic stuf you need for graphics3d :
'' 3dengine , math stuf , matrix , double3d and colors

public module basis3D

Public Class Dbl3D

    Public x As Double
    Public y As Double
    Public z As Double

    Public Sub New()
        x = 0
        y = 0
        z = 0
    End Sub

    Public Sub New(a As Double, b As Double, c As Double)
        x = a
        y = b
        z = c
    End Sub

    Public Sub fill(a As Double, b As Double, c As Double)
        x = a
        y = b
        z = c
    End Sub

    Public Function toColor() As Color
        Return Color.FromArgb(255 _
        , CInt(x * 255), CInt(y * 255), CInt(z * 255))
    End Function

    Public Shared Operator +(a As Dbl3D, b As Dbl3D) As Dbl3D
        Return New Dbl3D(a.x + b.x, a.y + b.y, a.z + b.z)
    End Operator

    Public Shared Operator -(a As Dbl3D, b As Dbl3D) As Dbl3D
        Return New Dbl3D(a.x - b.x, a.y - b.y, a.z - b.z)
    End Operator

    Public Shared Operator /(a As Dbl3D, b As Double) As Dbl3D
        Return New Dbl3D(a.x / b, a.y / b, a.z / b)
    End Operator

    Public Shared Operator *(a As Dbl3D, b As Double) As Dbl3D
        Return New Dbl3D(a.x * b, a.y * b, a.z * b)
    End Operator

    Public Function dot(b As Dbl3D) As Double
        Return x * b.x + y * b.y + z * b.z
    End Function

    Public Function cross(b As Dbl3D) As Dbl3D
        Return New Dbl3D(y * b.z - z * b.y _
                        ,z * b.x - x * b.z _
                        ,x * b.y - y * b.x)
    End Function

    Public Function length() As Double
        Return Math.sqrt(dot(Me))
    End Function

    Public Function getAngle(b As Dbl3D) As Double
        Dim la as Double length()
        Dim lb as Double b.length()
        Dim d as Double dot(b)
        Return Math.acs(d / (la * lb))
    End Function 

    Public Sub normalize() 
        Me = Me / length()
    End Sub

end class

Public Class m4x4

    Public m(3, 3) As Double

    Public Sub New()
            Dim i as Integer
            Dim j as Integer
            For i = 0 to 3
                For j = 0 to 3
                    m(i, j) = 0.0
                Next j
                m(i, i) = 1.0
            Next i
    End Sub

    Public Shared Operator *(a As m4x4, b As m4x4) As m4x4
        Dim uit As m4x4
        Dim i as Int16
        Dim j as Int16
        Dim k as Int16
        For i = 0 to 3
            For j = 0 to 3
                uit.m(i, j) = 0.0
                For k = 0 to 3
                    uit.m(i, j) += a.m(i, k) * b.m(k, j)
                Next k
            Next j
        Next i
        Return uit
    End Operator

    public shared operator *(m as m4x4, vv as dbl3d) as dbl3d
        return new dbl3d(vv.x*m.m(0,0)+vv.y*m.m(1,0) _
                        +vv.z*m.m(2,0)+m.m(3,0) _
                        ,vv.x*m.m(0,1)+vv.y*m.m(1,1)  _
                        +vv.z*m.m(2,1)+m.m(3,1) _
                        ,vv.x*m.m(0,2)+vv.y*m.m(1,2) _
                        +vv.z*m.m(2,2)+m.m(3,2))
    end operator

    public function inverse(q as m4x4) as m4x4
    '' for the inverse rotation matrix
    '' it wil NOT calculate inverse translation
        dim d as double 
        dim uit as m4x4 
  
        d = q.m(0,0)*q.m(1,1)*q.m(2,2) _
          - q.m(0,0)*q.m(2,1)*q.m(1,2) _
          + q.m(1,0)*q.m(2,1)*q.m(0,2) _
          - q.m(1,0)*q.m(0,1)*q.m(2,2) _
          + q.m(2,0)*q.m(0,1)*q.m(1,2) _
          - q.m(2,0)*q.m(1,1)*q.m(0,2)
    
        uit.m(0,0)=(q.m(1,1)*q.m(2,2)-q.m(1,2)*q.m(2,1))/d
        uit.m(1,0)=(q.m(0,1)*q.m(2,2)-q.m(0,2)*q.m(2,1))/d
        uit.m(2,0)=(q.m(0,1)*q.m(1,2)-q.m(0,2)*q.m(1,1))/d
        uit.m(0,1)=(q.m(1,0)*q.m(2,2)-q.m(1,2)*q.m(2,0))/d
        uit.m(1,1)=(q.m(0,0)*q.m(2,2)-q.m(0,2)*q.m(2,0))/d
        uit.m(2,1)=(q.m(0,0)*q.m(1,2)-q.m(0,2)*q.m(1,0))/d
        uit.m(0,2)=(q.m(1,0)*q.m(2,1)-q.m(1,1)*q.m(2,0))/d
        uit.m(1,2)=(q.m(0,0)*q.m(2,1)-q.m(0,1)*q.m(2,0))/d
        uit.m(2,2)=(q.m(0,0)*q.m(1,1)-q.m(0,1)*q.m(1,2))/d

        uit.translate(0,0,0)
  
        return uit
  
    end function

    Public Sub translate(x As Double _
    ,y As Double ,z As Double)
        m(3, 0) = x
        m(3, 1) = y
        m(3, 2) = z
    End Sub

    Public Sub rotate(deg As Double _
    , ax As Int16 = 0)
        Dim a As Int16 = 1
        Dim b As Int16 = 2
        If ax = 1 Then
            a = 0
            b = 2
        End If
        If ax = 2 Then
            a = 0
            b = 1
        End If
        m(a, a) = Math.cos(rad(deg))
        m(a, b) = -Math.sin(rad(deg))
        m(b, a) = Math.sin(rad(deg))
        m(b, b) = Math.cos(rad(deg))
    End Sub

    Public Sub scale(x As Double _
    ,y As Double ,z As double)
        m(0, 0) = x
        m(1, 1) = y
        m(2, 2) = z
    End Sub

end class

    public dim sk(64) as Dbl3D
    public dim v(20) as m4x4
    public number as int16

    public const xyz as int16 = 0
    public const xzy as int16 = 1
    public const yxz as int16 = 2
    public const yzx as int16 = 3
    public const zxy as int16 = 4
    public const zyx as int16 = 5

    public const golden_ratio as double = (math.sqrt(5)-1)/2

    ''primary colors
    public const black as new dbl3d(0,0,0)
    public const red as new dbl3d(1,0,0)
    public const green as new dbl3d(0,1,0)
    public const yellow as new dbl3d(1,1,0)
    public const blue as new dbl3d(0,0,1)
    public const magenta as new dbl3d(1,0,1)
    public const cyan as new dbl3d(0,1,1)
    public const white as new dbl3d(1,1,1)
    ''mixed colors
    public const gray as dbl3d = white / 2
    public const orange as new dbl3d(1,0.5,0)
    public const pink as new dbl3d(1,0.5,0.5)
    public const purple as dbl3d = magenta / 2
 
    public function rainbow(deg as double)as dbl3d
        return new dbl3d(sin(rad(deg))/2+.5 _
        ,sin(rad(deg-120))/2+.5 _
        ,sin(rad(deg+120))/2+.5)
    end function 

    public function rad(deg as double)as double
        return deg*math.pi/180.0
    end function

    public function degrees(r as double)as double
        return r*180.0/math.pi
    end function
  
    public sub rotate(byref k as double _
    , byref l as double,deg as double)
       dim s as double = math.sin(rad(deg))
       dim c as double = math.cos(rad(deg))
       dim hk as double = k*c-l*s
       dim hl as double = k*s+l*c
       k=hk
       l=hl
    end sub

    public function nr(no as int32,m as int16)as string
        return right("0000000000"&no,m)
    end function

    public sub link(no as int16 _
    ,x as double,y as double,z as double _
    ,xz as double,yz as double,xy as double _
    ,ax as int16,p as int16)
    '' create new matrix
    '' wil effect future drawings 
        if no < 1 or n > 20 then exit sub
        if p < 0 or p > 20 then exit sub
        if p = no then exit sub
        dim trans as m4x4
        dim rotx as m4x4
        dim roty as m4x4
        dim rotz as m4x4
        trans.translate(x,y,z)
        rotx.rotate(yz)
        roty.rotate(xz,1)
        rotz.rotate(xy,2)
        select case ax
            case xyz
                v(no)=rotx*roty*rotz*trans*v(p)
            case xzy
                v(no)=rotx*rotz*roty*trans*v(p)
            case yxz
                v(no)=roty*rotx*rotz*trans*v(p)
            case yzx
                v(no)=roty*rotz*rotx*trans*v(p)
            case zxy
                v(no)=rotz*rotx*roty*trans*v(p)
            case zyx
                v(no)=rotz*roty*rotx*trans*v(p)
            case else
                v(no)=rotx*roty*rotz*trans*v(p)
        end select
        number = no
    end sub

    public sub child(no as int16 _
    ,x as double,y as double,z as double _
    ,lim as int16,ax as int16,p as int16)
    '' create a lim matrix of a avatar
    '' wil effect furure drawings
        link(no,x,y,z,sk(lim).y,sk(lim).x,sk(lim).z),ax,p)
    end sub

    public sub skelet(lim as int16 _
    ,x as double,y as double,z as double)
    '' set the skelet lim of a avatar
        sk(lim)=Dbl3D(x,y,z)
    end sub

    public function pend( f as double _
    , amp as double ) as double
    '' for smoot animations of lims
        return math.sin(rad(f))*amp
    end function

    public sub spot(byref x as double _
    ,byref y as double,byref z as double)
    '' from local to world coordinates
        dim v3d as dbl3d(x,y,z)
        v3d=v(number)*v3d
        x=v3d.x
        y=v3d.y
        z=v3d.z
    end sub

    public sub local(byref x as double _
    ,byref y as double,byref z as double)
    '' from world to local coordinates
    '' only rotation matrix
        dim v3d as dbl3d(x,y,z)
        dim m44 as m4x4 
        m44 = m44.inverse(v(number))
        v3d=m44*v3d
        x=v3d.x
        y=v3d.y
        z=v3d.z
    end sub

end module

'' bluatigro 1 sept 2017
'' ray world module
'' needs basis3d

public module ray_world

    public const verry_big as double = 1e9
    public const verry_small as double = 1e-9

    public spheretel as int16
    public const spheremax as int16 = 200
    public spheres(spheremax) as ray_sphere

    public light as new dbl3d(-1,1,-1)

    public class ray_mat
    '' ray material
    '' wil be extended

        public diffuse as dbl3d
        public reflection as double

    end class

public const pntmax as int16 = 256
dim pnt(pntmax) as dbl3d

public sub set_point(no as int16,x as double,y as double,z as double)
    if no < 0 or no > pntmax then exit
    spot x,y,z
    pnt(no).fill x,y,z
end sub

class ray_triangle
  dim as dbl3d punt( 2 ) , n , led 
  dim as ray_mat mat
sub fill( p1 as integer , p2 as integer _
  , p3 as integer , m as ray_mat )
  if p1 < 0 or p1 > pntmax then exit sub
  if p2 < 0 or p2 > pntmax then exit sub
  if p3 < 0 or p3 > pntmax then exit sub
  punt( 0 ) = pnt( p1 )
  punt( 1 ) = pnt( p2 )
  punt( 2 ) = pnt( p3 )
  led = ( punt( 0 ) + punt( 1 ) + punt( 2 ) ) / 3
  dim a as dbl3d = punt( 2 ) - punt( 0 )
  n = a.cross( punt( 1 ) - punt( 0 ) )
  n *= -1
  n.normalize()
  mat = m
end sub
function hit( o as t3d , d as t3d ) as double
  dim as t3d e1 = punt( 1 ) - punt( 0 )
  dim as t3d e2 = punt( 2 ) - punt( 0 )
  dim as t3d p = d \ e2
  dim as double a = dot( e1 , p )
  if abs( a ) < 1e-9 then return -1
  dim as double f = 1 / a
  dim as t3d s = o - punt( 0 )
  dim as double u = f * dot( s , p )
  if u < 0 or u > 1 then return -1
  dim as t3d q = s \ e1
  dim as double v = f * dot( d , q )
  if v < 0 or u + v > 1 then return -1
  return f * dot( e2 , q )
end function
public const trimax as int16 = 1000
public tris(trimax) as ray_triangles
public tritel as int16
public sub tri(p1 as int16,p2 as int16,p3 as int16 _
,kl as dbl3d,f as double)
    if tritel >= trimax then exit
    dim mat as ray_mat
    mat.diffuse = kl
    mat.reflection = f
    tris(tritel).fill p1,p2,p3,mat
    tritel += 1
end sub
end class
    public class ray_sphere
    '' ray sphere shape
        public center as dbl3d
        public r as double
        public r2 as double
        public mat as ray_mat
        
        public sub new()
        end sub

        public sub fill(cx as double _
        ,cy as double,cz as double,r as double _
        ,m as ray_mat)
            spot(cx,cy,cz)
            center = dbl3d(cx,cy,cz)
            me.r = r
            r2 = r * r
            mat = m
        end sub

        public function hit( o as t3d , d as t3d ) as double
        '' calculate distance to sphere
        '' or return verry_big if sphere wasn't hit
            dim as double t , a , b , c , disc
            dim as t3d temp = o - center
            a = d.dot( d )
            b = 2 * temp.dot( d )
            c = temp.dot( temp ) - r2
            disc = b * b - 4 * a * c
            if disc < 0 then 
                return verry_big
            else
                dim as double e = math.sqrt( disc ) 
                dim as double demon = 2 * a
                t = ( -b - e ) / demon
                if t > verry_small then
                    return t
                end if
                t = ( -b + e ) / demon
                if t > verry_small then 
                    return t
                end if
            end if
            return verry_big
        end function

    end class 

    public sub sphere(x as double _
    ,y as double,z as double,r as double,kl as dbl3d _
    ,f as double = 0)
    '' add a sphere to the world 
        if spheretel >= spheremax then exit sub
        dim mat as ray_mat
        mat.diffuse = kl
        mat.reflection = f
        sphere(spheretel).fill(x,y,z,r,mat)
        spheretel += 1
    end sub

    public function render(o as dbl3d _
    ,d as dbl3d,depth as int16)as dbl3d
    '' get color of pixel
        dim i as int16
        dim isph as int16 = -1
        dim itri as int16 = -1
        dim b as double
        dim sphdist as double = verry_big
        dim tridist as double = verry_big
        '' get clostest sphere distance
        for i = 0 to spheretel - 1
            b = spheres(i).hit(o,d)
            if b < sphdist then
                sphdist = b
                isph = i
            end if
        next i
        for i = 0 to tritel - 1
            b = tris(i).hit(o,d)
            if b < tridist then
                tridist = b
                itri = i
            end if
        next i
        '' get color of clostest sphere or triangle
        '' and ca***ate shaded color
        '' wil be extended
        if sphdist < tridist
          if isph >= 0 then
            dim kl as dbl3d = spheres(isph).mat.diffuse
            dim p as dbl3d = o + d.normalize() * b
            dim n as dbl3d = p - spheres(isph).center
            dim a as double = n.getangle(light)
            return kl*(math.cos(a)/2+.5)
          end if
        else
          if itri >= 0 then
            dim kl as dbl3d = tris(itri).mat.diffuse
            dim p as dbl3d = o + d.nomalize() * b
            dim n as dbl3d = tris(itri).n
            return kl*(math.cos(a)/2+.5)
          end if
        end if
        '' if no sphere was hit return black
        return black
    end function

end module

        spheres(spheretel).fill(x, y, z, r, mat)

Public Class Form1
    Private Sub Form1_Load(sender As Object, e As EventArgs) Handles MyBase.Load
        Me.Width = 800
        Me.Height = 600
        spheretel = 0
        tritel = 0
        sphere(0, 0, 0, 50, yellow)
        sphere(0, -Me.Height / 2 - verry_big, 0, verry_big, green)
    End Sub

    Private Sub Form1_Paint(sender As Object _
                            , e As PaintEventArgs) Handles Me.Paint
        Dim x As Double
        Dim y As Double
        Dim ix As Integer
        Dim iy As Integer
        Dim o As New Dbl3D
        Dim d As New Dbl3D
        Dim p As New Pen(Color.Red)
        For x = -Me.Width / 2 To Me.Width / 2
            For y = -Me.Height / 2 To Me.Height / 2
                o.fill(0, 0, -1000)
                d.fill(x, y, 1000)
                p.Color = render(o, d, 7).toColor
                ix = CInt(x) + Me.Width / 2
                iy = CInt(y) + Me.Height / 2
                e.Graphics.DrawLine(p, ix, iy, ix, iy + 1)
            Next
        Next
    End Sub
End Class

'' bluatigro 1 sept 2017
'' ray world module
'' needs basis3d

public module ray_world

    public const verry_big as double = 1e9
    public const verry_small as double = 1e-9

    public spheretel as int16
    public const spheremax as int16 = 200
    public spheres(spheremax) as ray_sphere

    public light as new dbl3d(-1,1,-1)

    public class ray_mat
    '' ray material
    '' wil be extended

        public diffuse as dbl3d
        public reflection as double

    end Class

    Public Const pntmax As Int16 = 256
    Public pnt(pntmax) As Dbl3D

    Public Sub set_point(no As Int16, x As Double _
            , y As Double, z As Double)
        If no < 0 Or no > pntmax Then Exit Sub
        spot(x, y, z)
        pnt(no).fill(x, y, z)
    End Sub

    Public Class ray_triangle
        Private punt(2) As Dbl3D
        Public n As New Dbl3D
        Public mat As New ray_mat
        Sub fill(p1 As Integer, p2 As Integer _
        , p3 As Integer, m As ray_mat)
            If p1 < 0 Or p1 > pntmax Then Exit Sub
            If p2 < 0 Or p2 > pntmax Then Exit Sub
            If p3 < 0 Or p3 > pntmax Then Exit Sub
            punt(0) = pnt(p1)
            punt(1) = pnt(p2)
            punt(2) = pnt(p3)
            n = cross(punt(1) - punt(0), punt(2) - punt(0))
            n *= -1
            n = normalize(n)
            mat = m
        End Sub
        Function hit(o As Dbl3D, d As Dbl3D) As Double
            Dim e1 As Dbl3D = punt(1) - punt(0)
            Dim e2 As Dbl3D = punt(2) - punt(0)
            Dim p As Dbl3D = cross(d, e2)
            Dim a As Double = dot(e1, p)
            If Math.Abs(a) < verry_small Then Return verry_big
            Dim f As Double = 1 / a
            Dim s As Dbl3D = o - punt(0)
            Dim u As Double = f * dot(s, p)
            If u < 0 Or u > 1 Then Return verry_big
            Dim q As Dbl3D = cross(s, e1)
            Dim v As Double = f * dot(d, q)
            If v < 0 Or u + v > 1 Then Return verry_big
            Return f * dot(e2, q)
        End Function
        Public Const trimax As Int16 = 1000
        Public tris(trimax) As ray_triangle
        Public tritel As Int16
        Public Sub tri(p1 As Int16, p2 As Int16, p3 As Int16 _
        , kl As Dbl3D, flect As Double)
            If tritel >= trimax Then Exit Sub
            Dim mat As New ray_mat
            mat.diffuse = kl
            mat.reflection = flect
            tris(tritel).fill(p1, p2, p3, mat)
            tritel += 1
        End Sub
    End Class
    Public Const trimax As Int16 = 1000
    Public tris(trimax) As ray_triangle
    Public tritel As Int16
    Public Class ray_sphere
        '' ray sphere shape
        Public center As Dbl3D
        Public r As Double
        Public r2 As Double
        Public mat As ray_mat

        Public Sub New()
        End Sub

        Public Sub fill(cx As Double _
        , cy As Double, cz As Double, r As Double _
        , m As ray_mat)
            spot(cx, cy, cz)
            center.fill(cx, cy, cz)
            Me.r = r
            r2 = r * r
            mat = m
        End Sub

        Public Function hit(o As Dbl3D, d As Dbl3D) As Double
            '' calculate distance to sphere
            '' or return verry_big if sphere wasn't hit
            Dim t As Double
            Dim a As Double
            Dim b As Double
            Dim c As Double
            Dim disc As Double
            Dim temp As Dbl3D = o - center
            a = dot(d, d)
            b = 2 * dot(temp, d)
            c = dot(temp, temp) - r2
            disc = b * b - 4 * a * c
            If disc < 0 Then
                Return verry_big
            Else
                Dim e As Double = Math.Sqrt(disc)
                Dim demon As Double = 2 * a
                t = (-b - e) / demon
                If t > verry_small Then
                    Return t
                End If
                t = (-b + e) / demon
                If t > verry_small Then
                    Return t
                End If
            End If
            Return verry_big
        End Function

    End Class

    Public Sub sphere(x As Double _
    , y As Double, z As Double, r As Double, kl As Dbl3D)
        '' add a sphere to the world 
        If spheretel >= spheremax Then Exit Sub
        Dim mat As New ray_mat
        mat.diffuse = kl
        mat.reflection = 0
        spheres(spheretel).fill(x, y, z, r, mat)
        spheretel += 1
    End Sub

    Public function render(o as dbl3d _
    ,d as dbl3d,depth as int16)as dbl3d
    '' get color of pixel
        dim i as int16
        dim isph as int16 = -1
        dim itri as int16 = -1
        dim b as double
        dim sphdist as double = verry_big
        dim tridist as double = verry_big
        '' get clostest sphere distance
        for i = 0 to spheretel - 1
            b = spheres(i).hit(o,d)
            if b < sphdist then
                sphdist = b
                isph = i
            end if
        next i
        for i = 0 to tritel - 1
            b = tris(i).hit(o,d)
            if b < tridist then
                tridist = b
                itri = i
            end if
        next i
        '' get color of clostest sphere or triangle
        '' and ca***ate shaded color
        '' wil be extended
        if sphdist < tridist
          if isph >= 0 then
            dim kl as dbl3d = spheres(isph).mat.diffuse
                Dim p As Dbl3D = o + normalize(d) * b
                Dim n as dbl3d = p - spheres(isph).center
                Dim a As Double = getAngle(n, light)
                Return kl*(math.cos(a)/2+.5)
          end if
        else
          if itri >= 0 then
            dim kl as dbl3d = tris(itri).mat.diffuse
                Dim p As Dbl3D = o + normalize(d) * b
                Dim n As Dbl3D = tris(itri).n
                Dim a As Double = getAngle(n, light)
                Return kl*(math.cos(a)/2+.5)
          end if
        end If
        '' if no sphere or triangle was hit return black
        Return black
    end function

end module

'' bluatigro 3 sept 2017
'' all basic stuf you need for graphics3d :
'' 3dengine , math stuf , matrix , double3d and colors

Public Module basis3D

    Public Class Dbl3D

        Public x As Double
        Public y As Double
        Public z As Double

        Public Sub New()
            x = 0
            y = 0
            z = 0
        End Sub

        Public Sub New(a As Double, b As Double, c As Double)
            x = a
            y = b
            z = c
        End Sub

        Public Sub fill(a As Double, b As Double, c As Double)
            x = a
            y = b
            z = c
        End Sub

        Public Function toColor() As Color
            Return Color.FromArgb(255 _
            , CInt(x * 255), CInt(y * 255), CInt(z * 255))
        End Function

        Public Shared Operator +(a As Dbl3D, b As Dbl3D) As Dbl3D
            Return New Dbl3D(a.x + b.x, a.y + b.y, a.z + b.z)
        End Operator

        Public Shared Operator -(a As Dbl3D, b As Dbl3D) As Dbl3D
            Return New Dbl3D(a.x - b.x, a.y - b.y, a.z - b.z)
        End Operator

        Public Shared Operator /(a As Dbl3D, b As Double) As Dbl3D
            Return New Dbl3D(a.x / b, a.y / b, a.z / b)
        End Operator

        Public Shared Operator *(a As Dbl3D, b As Double) As Dbl3D
            Return New Dbl3D(a.x * b, a.y * b, a.z * b)
        End Operator

    End Class

    Public Function normalize(a As Dbl3D) As Dbl3D
        Return a / length(a)
    End Function

    Public Function dot(a As Dbl3D, b As Dbl3D) As Double
        Return a.x * b.x + a.y * b.y + a.z * b.z
    End Function

    Public Function cross(a As Dbl3D, b As Dbl3D) As Dbl3D
        Return New Dbl3D(a.y * b.z - a.z * b.y _
                        , a.z * b.x - a.x * b.z _
                        , a.x * b.y - a.y * b.x)
    End Function

    Public Function length(a As Dbl3D) As Double
        Return Math.Sqrt(dot(a, a))
    End Function

    Public Function getAngle(a As Dbl3D, b As Dbl3D) As Double
        Dim la As Double = length(a)
        Dim lb As Double = length(b)
        Dim d As Double = dot(a, b)
        Return Math.Acos(d / (la * lb))
    End Function

    Public Class m4x4

        Public m(3, 3) As Double

        Public Sub identity()
            Dim i As Integer
            Dim j As Integer
            For i = 0 To 3
                For j = 0 To 3
                    m(i, j) = 0.0
                Next j
                m(i, i) = 1.0
            Next i
        End Sub

        Public Shared Operator *(a As m4x4, b As m4x4) As m4x4
            Dim uit As New m4x4()
            Dim i As Int16
            Dim j As Int16
            Dim k As Int16
            For i = 0 To 3
                For j = 0 To 3
                    uit.m(i, j) = 0.0
                    For k = 0 To 3
                        uit.m(i, j) += a.m(i, k) * b.m(k, j)
                    Next k
                Next j
            Next i
            Return uit
        End Operator

        Public Shared Operator *(m As m4x4, v As Dbl3D) As Dbl3D
            Return New Dbl3D(v.x * m.m(0, 0) + v.y * m.m(1, 0) _
                            + v.z * m.m(2, 0) + m.m(3, 0) _
                            , v.x * m.m(0, 1) + v.y * m.m(1, 1) _
                            + v.z * m.m(2, 1) + m.m(3, 1) _
                            , v.x * m.m(0, 2) + v.y * m.m(1, 2) _
                            + v.z * m.m(2, 2) + m.m(3, 2))
        End Operator

        Public Function inverse(q As m4x4) As m4x4
            '' for the inverse rotation matrix only
            '' it wil NOT calculate inverse translation
            Dim d As Double
            Dim uit As New m4x4()

            d = q.m(0, 0) * q.m(1, 1) * q.m(2, 2) _
              - q.m(0, 0) * q.m(2, 1) * q.m(1, 2) _
              + q.m(1, 0) * q.m(2, 1) * q.m(0, 2) _
              - q.m(1, 0) * q.m(0, 1) * q.m(2, 2) _
              + q.m(2, 0) * q.m(0, 1) * q.m(1, 2) _
              - q.m(2, 0) * q.m(1, 1) * q.m(0, 2)

            uit.m(0, 0) = (q.m(1, 1) * q.m(2, 2) - q.m(1, 2) * q.m(2, 1)) / d
            uit.m(1, 0) = (q.m(0, 1) * q.m(2, 2) - q.m(0, 2) * q.m(2, 1)) / d
            uit.m(2, 0) = (q.m(0, 1) * q.m(1, 2) - q.m(0, 2) * q.m(1, 1)) / d
            uit.m(0, 1) = (q.m(1, 0) * q.m(2, 2) - q.m(1, 2) * q.m(2, 0)) / d
            uit.m(1, 1) = (q.m(0, 0) * q.m(2, 2) - q.m(0, 2) * q.m(2, 0)) / d
            uit.m(2, 1) = (q.m(0, 0) * q.m(1, 2) - q.m(0, 2) * q.m(1, 0)) / d
            uit.m(0, 2) = (q.m(1, 0) * q.m(2, 1) - q.m(1, 1) * q.m(2, 0)) / d
            uit.m(1, 2) = (q.m(0, 0) * q.m(2, 1) - q.m(0, 1) * q.m(2, 0)) / d
            uit.m(2, 2) = (q.m(0, 0) * q.m(1, 1) - q.m(0, 1) * q.m(1, 2)) / d

            uit.translate(0, 0, 0)

            Return uit

        End Function

        Public Sub translate(x As Double _
        , y As Double, z As Double)
            m(3, 0) = x
            m(3, 1) = y
            m(3, 2) = z
        End Sub

        Public Sub rotate(deg As Double _
        , ax As Int16)
            Dim a As Int16 = 1
            Dim b As Int16 = 2
            If ax = 1 Then
                a = 0
                b = 2
            End If
            If ax = 2 Then
                a = 0
                b = 1
            End If
            m(a, a) = Math.Cos(rad(deg))
            m(a, b) = -Math.Sin(rad(deg))
            m(b, a) = Math.Sin(rad(deg))
            m(b, b) = Math.Cos(rad(deg))
        End Sub

        Public Sub scale(x As Double _
        , y As Double, z As Double)
            m(0, 0) = x
            m(1, 1) = y
            m(2, 2) = z
        End Sub

    End Class

    Private sk(64) As Dbl3D
    Private v(20) As m4x4
    Private number As Int16

    Public Const xyz As Int16 = 0
    Public Const xzy As Int16 = 1
    Public Const yxz As Int16 = 2
    Public Const yzx As Int16 = 3
    Public Const zxy As Int16 = 4
    Public Const zyx As Int16 = 5

    Public golden_ratio As Double = (Math.Sqrt(5) - 1) / 2

    ''primary colors
    Public black As New Dbl3D(0, 0, 0)
    Public red As New Dbl3D(1, 0, 0)
    Public green As New Dbl3D(0, 1, 0)
    Public yellow As New Dbl3D(1, 1, 0)
    Public blue As New Dbl3D(0, 0, 1)
    Public magenta As New Dbl3D(1, 0, 1)
    Public cyan As New Dbl3D(0, 1, 1)
    Public white As New Dbl3D(1, 1, 1)
    ''mixed colors
    Public gray As Dbl3D = white / 2
    Public orange As New Dbl3D(1, 0.5, 0)
    Public pink As New Dbl3D(1, 0.5, 0.5)
    Public purple As Dbl3D = magenta / 2

    Public Function rainbow(deg As Double) As Dbl3D
        Return New Dbl3D(Math.Sin(rad(deg)) / 2 + 0.5 _
        , Math.Sin(rad(deg - 120)) / 2 + 0.5 _
        , Math.Sin(rad(deg + 120)) / 2 + 0.5)
    End Function

    Public Function rad(deg As Double) As Double
        Return deg * Math.PI / 180.0
    End Function

    Public Function degrees(r As Double) As Double
        Return r * 180.0 / Math.PI
    End Function

    Public Sub rotate(ByRef k As Double _
    , ByRef l As Double, deg As Double)
        Dim s As Double = Math.Sin(rad(deg))
        Dim c As Double = Math.Cos(rad(deg))
        Dim hk As Double = k * c - l * s
        Dim hl As Double = k * s + l * c
        k = hk
        l = hl
    End Sub

    Public Function nr(no As Int32, m As Int16) As String
        Return Right("0000000000" & no, m)
    End Function

    Public Sub link(no As Int16 _
    , x As Double, y As Double, z As Double _
    , xz As Double, yz As Double, xy As Double _
    , ax As Int16, p As Int16)
        '' create new matrix
        '' wil effect future drawings 
        If no < 1 Or no > 20 Then Exit Sub
        If p < 0 Or p > 20 Then Exit Sub
        If p = no Then Exit Sub
        Dim trans As New m4x4
        Dim rotx As New m4x4
        Dim roty As New m4x4
        Dim rotz As New m4x4
        trans.identity()
        trans.translate(x, y, z)
        rotx.identity()
        rotx.rotate(yz, 0)
        roty.identity()
        roty.rotate(xz, 1)
        rotz.identity()
        rotz.rotate(xy, 2)
        Select Case ax
            Case xyz
                v(no) = rotx * roty * rotz * trans * v(p)
            Case xzy
                v(no) = rotx * rotz * roty * trans * v(p)
            Case yxz
                v(no) = roty * rotx * rotz * trans * v(p)
            Case yzx
                v(no) = roty * rotz * rotx * trans * v(p)
            Case zxy
                v(no) = rotz * rotx * roty * trans * v(p)
            Case zyx
                v(no) = rotz * roty * rotx * trans * v(p)
            Case Else
                v(no) = rotx * roty * rotz * trans * v(p)
        End Select
        number = no
    End Sub

    Public Sub child(no As Int16 _
    , x As Double, y As Double, z As Double _
    , lim As Int16, ax As Int16, p As Int16)
        '' create a lim matrix of a avatar
        '' wil effect furure drawings
        link(no, x, y, z, sk(lim).y, sk(lim).x, sk(lim).z, ax, p)
    End Sub

    Public Sub skelet(lim As Int16 _
    , x As Double, y As Double, z As Double)
        '' set the skelet lim of a avatar
        sk(lim) = New Dbl3D(x, y, z)
    End Sub

    Public Function pend(f As Double _
    , amp As Double) As Double
        '' for smoot animations of lims
        Return Math.Sin(rad(f)) * amp
    End Function

    Public Sub spot(ByRef x As Double _
    , ByRef y As Double, ByRef z As Double)
        '' from local to world coordinates
        Dim v3d As New Dbl3D(x, y, z)
        v3d = v(number) * v3d
        x = v3d.x
        y = v3d.y
        z = v3d.z
    End Sub

    Public Sub local(ByRef x As Double _
    , ByRef y As Double, ByRef z As Double)
        '' from world to local coordinates
        '' only rotation matrix
        Dim v3d As New Dbl3D(x, y, z)
        Dim m44 As New m4x4
        m44 = m44.inverse(v(number))
        v3d = m44 * v3d
        x = v3d.x
        y = v3d.y
        z = v3d.z
    End Sub

End Module