VERSION 5.00
Begin VB.Form Form1 
   AutoRedraw      =   -1  'True
   Caption         =   "Form1"
   ClientHeight    =   4065
   ClientLeft      =   45
   ClientTop       =   330
   ClientWidth     =   4950
   LinkTopic       =   "Form1"
   ScaleHeight     =   4065
   ScaleWidth      =   4950
   StartUpPosition =   3  'Windows Default
   Begin VB.PictureBox Picture1 
      AutoRedraw      =   -1  'True
      Height          =   3732
      Left            =   100
      ScaleHeight     =   245
      ScaleMode       =   3  'Pixel
      ScaleWidth      =   309
      TabIndex        =   0
      Top             =   100
      Width           =   4692
   End
End
Attribute VB_Name = "Form1"
Attribute VB_GlobalNameSpace = False
Attribute VB_Creatable = False
Attribute VB_PredeclaredId = True
Attribute VB_Exposed = False
' By Fred D.Jordan fred.jordan@span.ch
' Use as you wish, but just tell me:)

Private Type Vertex
    X As Double
    Y As Double
    z As Double
End Type
    
Private Type Col
    R As Byte
    G As Byte
    B As Byte
End Type
    
Dim V As Vertex
Dim P As Vertex
Dim VM As Vertex
Dim M As Vertex
Dim N As Vertex
Dim L As Vertex
Dim R As Vertex
Dim G As Vertex
Dim xyres As Double
Dim Size As Double
Dim Light As Col
Dim Radius As Double
Dim Chess As Double
Dim ChessSize As Double
Dim h As Double


Sub compute()
    Dim Y, z As Double
    
    st = 2 * Size / xyres
    
    For z = -Size + st / 2 To Size + st / 2 Step st
        ix = 0
        iy = iy + 1
        Picture1.Refresh
        For Y = -Size + st / 2 To Size + st / 2 Step st
        P.Y = Y
        P.z = z
        traceRay
        
        ix = ix + 1
        Picture1.PSet (ix, iy), _
            RGB(Light.R, Light.G, Light.B)
        Next Y
    Next z
End Sub

Sub traceRay()
    Dim A, B, C, delta, t, xi, yi As Double
    
    ' Compute the intersection point with the sphere
    A = (P.X - V.X) ^ 2 + (P.Y - V.Y) ^ 2 + (P.z - V.z) ^ 2
    B = 2 * (P.X - V.X) * V.X + 2 * (P.Y - V.Y) * V.Y + 2 * (P.z - V.z) * V.z
    C = V.X ^ 2 + V.Y ^ 2 + V.z ^ 2 - Radius ^ 2
    
    delta = B ^ 2 - 4 * A * C
    
    ' The ray missed the sphere ...
    If (delta <= 0) Then
        t = (V.z - h) / (P.z - V.z)
        xi = t * (P.X - V.X) + V.X
        yi = t * (P.Y - V.Y) + V.Y

        ' Deduct the color of the chessed plane
        If Abs(xi) < 100000000000000# And Abs(yi) < 100000000000000# Then
            fx = Int(Abs(xi) / Chess) Mod 2
            fy = Int(Abs(yi) / Chess) Mod 2
        Else
            fx = 0
            fy = 0
        End If
    
        ' The ray hits the plane
        If (t < 0) Then
            check = fx Xor fy
            If (check = 0) Then
                Light.R = 255
                Light.G = 0
                Light.B = 0
            End If
            If (check = 1) Then
                Light.R = 255
                Light.G = 255
                Light.B = 0
            End If
            
            ' Compute the intersection point with the sphere
            A = (xi - G.X) ^ 2 + (yi - G.Y) ^ 2 + (-h - G.z) ^ 2
            B = 2 * (xi - G.X) * G.X + 2 * (yi - G.Y) * G.Y + 2 * (-h - G.z) * G.z
            C = G.X ^ 2 + G.Y ^ 2 + G.z ^ 2 - Radius ^ 2
    
            delta1 = B ^ 2 - 4 * A * C
            
            If (delta1 > 0) Then
               Light.R = 0 * Light.R / 2
               Light.G = 0 * Light.G / 2
               Light.B = 0 * Light.B / 2
            End If
            
            d = Sqr((V.X - xi) ^ 2 + (V.Y - yi) ^ 2 + (V.z + h) ^ 2)
        ' The ray hits the sky
        Else
            Light.R = 10
            Light.G = 0
            Light.B = 160
            d = Exp(P.z + Size * 0.9)
        End If
        
        addFog (d)
        
        Exit Sub
    End If
        
    t = (-B - Sqr(delta)) / (2 * A)
    
    M.X = t * (P.X - V.X) + V.X
    M.Y = t * (P.Y - V.Y) + V.Y
    M.z = t * (P.z - V.z) + V.z
    
    VM.X = M.X - V.X
    VM.Y = M.Y - V.Y
    VM.z = M.z - V.z
    
    
    
    ' Compute the normalized vectors
    N = norm(M)
    L = norm(VM)
    
    ' Compute the reflected vector direction
    R = diff(sprod(sprod(N, prod(L, N)), 2), L)
    
    ' Compute the impact location on the plane
    t = (M.z - h) / R.z
    xi = t * R.X + M.X
    yi = t * R.Y + M.Y
    
    ' Deduct the color of the chessed plane
    If Abs(xi) < 100000000000000# And Abs(yi) < 100000000000000# Then
        fx = Int(Abs(xi) / Chess) Mod 2
        fy = Int(Abs(yi) / Chess) Mod 2
    Else
        fx = 0
        fy = 0
    End If
    
    check = fx Xor fy
    If (check = 0) Then
        Light.R = 200
        Light.G = 0
        Light.B = 0
    End If
    If (check = 1) Then
        Light.R = 200
        Light.G = 200
        Light.B = 0
    End If
    
    If (t < 0) Then
        Light.R = 10
        Light.G = 0
        Light.B = 10
    End If
    
    If (Abs(xi) > ChessSize) Then
        Light.R = 10
        Light.G = 0
        Light.B = 10
    End If
    
    If (Abs(yi) > ChessSize) Then
        Light.R = 10
        Light.G = 0
        Light.B = 10
    End If
    
    addLight = Int((R.X + R.Y + R.z) * 100)
    If addLight < 0 Then
        addLight = 0
    End If
    
    
    If (Light.R + addLight > 255) Then
        Light.R = 255
    Else
        Light.R = Light.R + addLight
    End If
    
    d = Sqr((VM.X) ^ 2 + (VM.Y) ^ 2 + (VM.z) ^ 2)
    addFog (d)
End Sub
        
Sub addFog(d As Double)
    fog = Exp(-0.1 * d)
    Light.R = Light.R * fog + (1 - fog) * 255
    Light.G = Light.G * fog + (1 - fog) * 255
    Light.B = Light.B * fog + (1 - fog) * 255
End Sub
    
        
Sub initVariable()
    Radius = 4
    Chess = 0.5
    ChessSize = 10
    h = 4.5
    
    V.X = 10
    V.Y = 0
    V.z = 0
    
    G.X = -200
    G.Y = 50
    G.z = 100
    
    P.X = 2
    xyres = 300
    Size = 5
End Sub


Private Sub Form_Load()
    Show
    Refresh
    Picture1.Visible = True
    initVariable
    compute
End Sub

Private Function norm(A As Vertex) As Vertex
    module = Sqr(A.X ^ 2 + A.Y ^ 2 + A.z ^ 2)
    
    norm.X = A.X / module
    norm.Y = A.Y / module
    norm.z = A.z / module
End Function

Private Function prod(A As Vertex, B As Vertex) As Double
    prod = A.X * B.X + A.Y * B.Y + A.z * B.z
End Function

Private Function sprod(A As Vertex, s As Double) As Vertex
    sprod.X = A.X * s
    sprod.Y = A.Y * s
    sprod.z = A.z * s
End Function

Private Function diff(A As Vertex, B As Vertex) As Vertex
    diff.X = A.X - B.X
    diff.Y = A.Y - B.Y
    diff.z = A.z - B.z
End Function

Private Sub Form_MouseUp(Button As Integer, Shift As Integer, X As Single, Y As Single)
    compute
End Sub

Private Sub Form_Resize()
    Picture1.Height = ScaleWidth - 600
    Picture1.Width = ScaleHeight
    
    Width = Height
    xyres = Picture1.Width / Screen.TwipsPerPixelX
End Sub

Private Sub Picture1_Click()
    compute
End Sub